File: categoricalDefinitionTest.R

package info (click to toggle)
r-cran-openmx 2.21.1%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 14,412 kB
  • sloc: cpp: 36,577; ansic: 13,811; fortran: 2,001; sh: 1,440; python: 350; perl: 21; makefile: 5
file content (131 lines) | stat: -rw-r--r-- 3,889 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
# ===========
# = HISTORY =
# ===========
# 2017-04-14 04:57PM TBATES
# 	update for   mxExpectationRAM("A", "S", "F", "M", thresholds="thresholds"),mxFitFunctionML()
# 	fix "Function precision" bug
# has a bug (no dimnames) at 104
library(OpenMx)
library(mvtnorm)

# define sample size
n <- 250
# define loadings for simulation model when covariate=0
loadings1 <- c(1, 1, 1, 1, 1)
# define change in loadings dependent on covariate
loadings2 <- c(0, 0, 0, 0, -.5)
# define residual variance terms
resid <- c(.4, .4, .4, .4, .4)
# define definition variable
g <- runif(n, 0, 2)
# simulate data for testing
data <- matrix(NA, n, 6)
for (i in 1:n){
	l <- loadings1 + loadings2 * g[i]
	e <- matrix(0, 5, 5)
	diag(e) <- resid
	data[i,1:5] <- rnorm(1)*l + rmvnorm(1,rep(0,5),e)
}
data <- data.frame(data)
names(data) <- c(paste("x",1:5, sep=""),"g")
# constrain the data to categories
data[data < 0] <- 0
data[data > 0] <- 1
# add the definition variable
data[, 6] <- g

# view the correlation matrix of the data
cor(data)

# ================
# = Begin OpenMx =
# ================

## Let's try a threshold model without a definition variable##
A <- mxMatrix("Full", 6, 6, name="A")
A$values[1:5, 6] <- 1
A$free[1:5, 6]   <- TRUE
A$labels[1:5, 6] <- paste("l",1:5,sep="")

S <- mxMatrix("Symm", 6, 6, name="S")
diag(S$values) <- 1
diag(S$labels) <- c(paste("e",1:5,sep=""), "var.F")

F <- mxMatrix("Full", 5, 6, dimnames=list(names(data)[1:5], c(names(data)[1:5], "F")), name="F")
diag(F$values) <- 1

M <- mxMatrix("Zero", 1, 6, name="M")

# columns are variables; there are four thresholds for each variable
T <- mxMatrix("Full", 1, 5, TRUE, 0, dimnames = list(c(), names(data)[1:5]), name="thresholds")

catModel <- mxModel("Categorical Test", A, S, F, M, T,
	mxData(data, type="raw"),
  mxExpectationRAM("A", "S", "F", "M", thresholds="thresholds"),
	mxFitFunctionML()
)

catResults <- mxRun(catModel)

summary(catResults)

# item difficulties
catResults$output$estimate[6:10]/catResults$output$estimate[1:5]

# item discriminations
catResults$output$estimate[1:5]


## Let's try a threshold model with a definition variable##

# here's the trick. i created a dummy factor called "hack", 
# and labeled it's path on x5 with the name of the definition variable
# there's a flashier way involving mxAlgebra statements, but the required substitution isn't working
# look for "square bracket substitution" in the future
CA <- mxMatrix("Full", 7, 7, name="CA")
CA$values[c(1:5,7), 6] <- 1
CA$free[c(1:5,7), 6]   <- TRUE
CA$labels[c(1:5,7), 6] <- c(paste("l",1:5,sep=""), "dif")

CA$labels[5, 7] <-"data.g"
#trick <- mxMatrix("Full", 1, 2, TRUE, c(1,0), labels=c("alpha","beta"))
#def   <- mxAlgebra(alpha + beta*data.g, name = "l5")

CS <- mxMatrix("Symm", 7, 7, name="CS")
diag(CS$values) <- c(rep(1,6),0)
diag(CS$labels) <- c(paste("e",1:5,sep=""), "var.F", "var.hack")

CF <- mxMatrix("Full", 5, 7, dimnames=list(names(data)[1:5], c(names(data)[1:5], "F", "hack")), name="CF")
diag(CF$values) <- 1

# =================================================
# = FAILS to RUN DUE TO LACKING DIMNAMES ON MEANS =
# =================================================
CM <- mxMatrix("Zero", 1, 7, name = "CM",)

# columns are variables; there are four thresholds for each variable
CT <- mxMatrix("Full", 1, 5, TRUE, 0, 
    labels = paste("t", 1:5, sep=""), 
    dimnames = list(c(), names(data)[1:5]), 
    name="thresholds")

catModel2 <- mxModel("Categorical Test with Def", CA, CS, CF, CM, CT, 
    mxData(data, type="raw"),
	  mxExpectationRAM("CA", "CS", "CF", "CM", thresholds="thresholds"),
		mxFitFunctionML()		
)

catModel2 <- mxOption(catModel2, "Function precision", 1e-9)

catResults2 <- mxRun(catModel2)

summary(catResults2)


# item difficulties
catResults2$output$estimate[6:10]/catResults2$output$estimate[1:5]

# item discriminations
catResults2$output$estimate[1:5]

catResults2$output