ifa: Store latent mean & cov in regular MxMatrix
[openmx:openmx.git] / models / failing / categoricalDefinitionTest.R
1 library(OpenMx)
2 library(mvtnorm)
3
4 # should we write files, including data files and mplus input files?
5 write.files <- FALSE
6 is.windows  <- FALSE
7 if(write.files==TRUE&is.windows==TRUE)  setwd("C:/Documents and Settings/cre2w/My Documents/Research/Kevin")
8 if(write.files==TRUE&is.windows==FALSE) setwd("/Users/ryneestabrook")
9
10 # define sample size
11 n <- 250
12
13 # define loadings for simulation model when covariate=0
14 loadings1 <- c(1, 1, 1, 1, 1)
15
16 # define change in loadings dependent on covariate
17 loadings2 <- c(0, 0, 0, 0, -.5)
18
19 # define residual variance terms
20 resid <- c(.4, .4, .4, .4, .4)
21
22 # define definition variable
23 g <- runif(n, 0, 2)
24
25 # simulate data for testing
26 data <- matrix(NA,n,6)
27
28 for (i in 1:n){
29         l <- loadings1+loadings2*g[i]
30         e <- matrix(0,5,5)
31         diag(e) <- resid
32         data[i,1:5] <- rnorm(1)*l + rmvnorm(1,rep(0,5),e)
33         }
34
35 data <- data.frame(data)
36
37 names(data) <- c(paste("x",1:5, sep=""),"g")
38
39 # constrain the data to categories
40 data[data<0] <- 0
41 data[data>0] <- 1
42
43 #add the definition variable
44 data[,6] <- g
45
46 # view the correlation matrix of the data
47 cor(data)
48
49
50 # write data to file
51 if(write.files==TRUE) write.table(data, "data.dat", row.names=FALSE, col.names=FALSE)
52
53 ### Begin OpenMx ###
54
55 ## Let's try a threshold model without a definition variable##
56 A <- mxMatrix("Full", 6, 6, name="A")
57 A@values[1:5, 6] <- 1
58 A@free[1:5, 6]   <- TRUE
59 A@labels[1:5, 6] <- paste("l",1:5,sep="")
60
61 S <- mxMatrix("Symm", 6, 6, name="S")
62 diag(S@values) <- 1
63 diag(S@labels) <- c(paste("e",1:5,sep=""), "var.F")
64
65 F <- mxMatrix("Full", 5, 6, 
66     dimnames=list(names(data)[1:5], c(names(data)[1:5], "F")), 
67     name="F")
68 diag(F@values) <- 1
69
70 M <- mxMatrix("Zero", 1, 6, name="M")
71
72 # columns are variables; there are four thresholds for each variable
73 T <- mxMatrix("Full", 1, 5, TRUE, 0, 
74     dimnames = list(c(), names(data)[1:5]), 
75     name="thresholds")
76
77 catModel <- mxModel("Categorical Test",
78     mxData(data, type="raw"),
79     A, S, F, M, T,
80     mxRAMObjective("A", "S", "F", "M", thresholds="thresholds")
81 )
82
83 catResults <- mxRun(catModel)
84
85 summary(catResults)
86
87 # item difficulties
88 catResults@output$estimate[6:10]/catResults@output$estimate[1:5]
89
90 # item discriminations
91 catResults@output$estimate[1:5]
92
93
94 ## Let's try a threshold model with a definition variable##
95
96 #here's the trick. i created a dummy factor called "hack", 
97 #and labeled it's path on x5 with the name of the definition variable
98 #there's a flashier way involving mxAlgebra statements, but the required substitution isn't working
99 #look for "square bracket substitution" in the future
100 CA <- mxMatrix("Full", 7, 7, name="CA")
101 CA@values[c(1:5,7), 6] <- 1
102 CA@free[c(1:5,7), 6]   <- TRUE
103 CA@labels[c(1:5,7), 6] <- c(paste("l",1:5,sep=""), "dif")
104
105 CA@labels[5, 7] <-"data.g"
106 #trick <- mxMatrix("Full", 1, 2, TRUE, c(1,0), labels=c("alpha","beta"))
107 #def   <- mxAlgebra(alpha + beta*data.g, name = "l5")
108
109 CS <- mxMatrix("Symm", 7, 7, name="CS")
110 diag(CS@values) <- c(rep(1,6),0)
111 diag(CS@labels) <- c(paste("e",1:5,sep=""), "var.F", "var.hack")
112
113 CF <- mxMatrix("Full", 5, 7, 
114     dimnames=list(names(data)[1:5], c(names(data)[1:5], "F", "hack")), 
115     name="CF")
116 diag(CF@values) <- 1
117
118 CM <- mxMatrix("Zero", 1, 7, name="CM")
119
120 # columns are variables; there are four thresholds for each variable
121 CT <- mxMatrix("Full", 1, 5, TRUE, 0, 
122     labels = paste("t", 1:5, sep=""), 
123     dimnames = list(c(), names(data)[1:5]), 
124     name="thresholds")
125
126 catModel2 <- mxModel("Categorical Test with Definition",
127     mxData(data, type="raw"),
128     CA, CS, CF, CM, CT, 
129     mxRAMObjective("CA", "CS", "CF", "CM", thresholds="thresholds")
130 )
131
132 catModel2 <- mxOption(catModel2, "Function Precision", 1e-9)
133
134 catResults2 <- mxRun(catModel2)
135
136 summary(catResults2)
137
138
139 # item difficulties
140 catResults2@output$estimate[6:10]/catResults2@output$estimate[1:5]
141
142 # item discriminations
143 catResults2@output$estimate[1:5]
144
145 catResults2@output