Add test for fitting a latent distribution with fixed item parameters
[openmx:openmx.git] / models / passing / ifa-drm-mg.R
1 library(OpenMx)
2 library(rpf)
3
4 set.seed(9)
5
6 numItems <- 30
7 i1 <- rpf.drm(multidimensional=TRUE)
8 items <- list()
9 items[1:numItems] <- i1
10 correct <- matrix(NA, 4, numItems)
11 for (x in 1:numItems) correct[,x] <- rpf.rparam(i1)
12 correct[1,] <- 1
13 correct[3,] <- logit(0)
14 correct[4,] <- logit(1)
15
16 ip.mat <- mxMatrix(name="itemParam", nrow=4, ncol=numItems,
17                    values=c(1,0, logit(0), logit(1)),
18                    free=c(FALSE, TRUE, FALSE, FALSE))
19
20 latent <- mxModel('latent',
21                   mxMatrix(nrow=1, ncol=1, free=TRUE, values=0, name="mean"),
22                   mxMatrix(type="Symm", nrow=1, ncol=1, free=TRUE, values=1, name="cov"),
23                   mxDataDynamic("cov", expectation="latentTest.expectation"),
24                   mxExpectationNormal(covariance="cov", means="mean"),
25                   mxFitFunctionML())
26
27 # We generate data in the "wrong" order to preserve compatibility
28 # with the previous version of the test.
29 data <- rpf.sample(500, items, correct, cov=matrix(5,1,1))
30
31 ldata <- rpf.sample(300, items, correct, mean=.5, cov=matrix(5,1,1))
32
33 ip.fix <- ip.mat
34 ip.fix$free[,] <- FALSE
35 ip.fix$values[,] <- correct
36
37 if (0) {
38   plan <- mxComputeEM('expectation', 'scores',
39                       mxComputeGradientDescent(paste('latent',c('mean','cov'), sep="."),
40                                                fitfunction="latent.fitfunction"),
41                       verbose=0L)
42 }
43
44 plan <- mxComputeIterate(list(
45   mxComputeGradientDescent(paste('latent',c('mean','cov'), sep="."),
46                            fitfunction="latent.fitfunction"),
47   mxComputeOnce('fitfunction', 'fit')))
48
49 m1 <- mxModel(model="latentTest", ip.fix, latent,
50               mxData(observed=ldata, type="raw"),
51               mxExpectationBA81(mean="latent.mean", cov="latent.cov",
52                                 ItemSpec=items, ItemParam="itemParam"),
53               mxFitFunctionML(),
54               plan)
55 m1Fit <- mxRun(m1)
56 omxCheckCloseEnough(m1Fit$output$estimate, c(.46, 4.44), .1)
57 omxCheckCloseEnough(m1Fit$output$fit, 8251.09, .01)
58
59 latent.plan <- NULL
60 if (1) {
61         latent.plan <- mxComputeGradientDescent('latent.cov', fitfunction="latent.fitfunction")
62 } else {
63         latent.plan <- 
64             mxComputeSequence(list(mxComputeOnce('expectation'),
65                                    mxComputeOnce('expectation', "latentDistribution", "copy"),
66                                    mxComputeOnce('fitfunction', "set-starting-values")),
67                               freeSet='latent.cov')
68 }
69
70 latent$mean$free[,] <- FALSE
71 latent$data$expectation <- "drmmg.expectation"
72
73 m2 <- mxModel(model="drmmg", ip.mat, latent,
74               mxData(observed=data, type="raw"),
75               mxExpectationBA81(mean="latent.mean", cov="latent.cov",
76                                 ItemSpec=items, ItemParam="itemParam"),
77               mxFitFunctionML(),
78               mxComputeEM('expectation', 'scores',
79                           mxComputeSequence(list(
80                               mxComputeNewtonRaphson(freeSet='itemParam'),
81                               latent.plan)),
82                           verbose=0L))
83
84 m2 <- mxRun(m2)
85
86 omxCheckCloseEnough(m2$output$fit, 14129.94, .01)
87 omxCheckCloseEnough(m2$submodels$latent$matrices$cov$values[1,1], 4.377, .01)
88                 
89 emstat <- m2$compute$output
90 omxCheckCloseEnough(emstat$EMcycles, 36, 1)
91 #omxCheckCloseEnough(emstat$totalMstep, 763, 20) # includes the latent distribution
92
93                                         #print(m2$matrices$itemParam$values)
94                                         #print(correct.mat)
95 mask <- is.finite(correct)
96 got <- cor(c(m2$matrices$itemParam$values[mask]),
97            c(correct[mask]))
98 omxCheckCloseEnough(got, .994, .01)
99
100 if (1) {
101   ip.mat <- mxMatrix(name="itemParam", nrow=4, ncol=numItems,
102                      values=c(1,0, logit(0), logit(1)),
103                      free=c(TRUE, TRUE, FALSE, FALSE))
104   ip.mat$labels[1,] <- 'a1'
105   
106   m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
107   cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=FALSE)
108
109   m2 <- mxModel(model="drmmg", ip.mat, m.mat, cov.mat,
110                 mxData(observed=data, type="raw"),
111                 mxExpectationBA81(mean="mean", cov="cov",
112                                   ItemSpec=items, ItemParam="itemParam"),
113                 mxFitFunctionML(),
114                 mxComputeSequence(steps=list(
115                   mxComputeOnce('expectation', 'scores'),
116                   mxComputeOnce('fitfunction', c('gradient', 'hessian', 'ihessian')),
117                   mxComputeReportDeriv()
118                 )))
119   deriv <- mxRun(m2, silent=TRUE)
120   omxCheckCloseEnough(deriv$output$ihessian, solve(deriv$output$hessian), 1e-4)
121   
122   if (0) {
123     m3 <- mxModel(m2, mxComputeSequence(list(
124       mxComputeOnce('fitfunction', 'fit'),
125       mxComputeNumericDeriv(parallel=FALSE, iterations=2L),
126       mxComputeReportDeriv())))
127     deriv <- mxRun(m3)
128     stop("ok")
129   }
130   
131   m2 <- mxModel(model="drmmg", ip.mat, m.mat, cov.mat,
132                 mxData(observed=data, type="raw"),
133                 mxExpectationBA81(mean="mean", cov="cov",
134                                   ItemSpec=items, ItemParam="itemParam"),
135                 mxFitFunctionML(),
136                 mxComputeEM('expectation', 'scores',
137                             mxComputeNewtonRaphson(freeSet='itemParam')))
138
139   m2 <- mxRun(m2)
140   emstat <- m2$compute$output
141   omxCheckCloseEnough(emstat$EMcycles, 38, 1)
142   omxCheckCloseEnough(emstat$totalMstep, 94, 5)
143   omxCheckCloseEnough(m2$fitfunction$result, 14129.04, .01)
144   omxCheckCloseEnough(m2$matrices$itemParam$values[1,], rep(2.133, numItems), .002)
145   # correct values are from flexMIRT
146   est <- c(-0.838622, -1.02653, -0.0868472, -0.251784, 0.953364,  0.735258, 0.606918,
147            1.04239, 0.466055, -2.05196, -0.0456446,  -0.320668, -0.362073, 2.02502,
148            0.635298, -0.0731132, -2.05196,  -0.0456446, -1.17429, 0.880002, -0.838622,
149            -0.838622, 1.02747,  0.424094, -0.584298, 0.663755, 0.663755, 0.064287, 1.38009,
150            1.01259 )
151   omxCheckCloseEnough(m2$matrices$itemParam$values[2,], est, .002)
152 }
153
154 if (0) {
155   library(mirt)
156   rdata <- sapply(data, unclass)-1
157   # for flexMIRT, write CSV
158   #write.table(rdata, file="ifa-drm-mg.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
159   pars <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars='values')
160   pars[pars$name=="a1",'value'] <- 1
161   pars[pars$name=="a1",'est'] <- FALSE
162   pars[pars$name=="COV_11",'est'] <- TRUE
163   fit <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars=pars)
164   # LL -7064.519 * -2 = 14129.04
165   got <- coef(fit)
166   print(got$GroupPars)
167   # COV 4.551
168   got$GroupPars <- NULL
169   round(m2$matrices$itemParam$values - simplify2array(got), 2)
170   
171   # MH-RM takes forever, not run
172   pars <- confmirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars='values')
173   pars[pars$name=="a1",'value'] <- 1
174   pars[pars$name=="a1",'est'] <- FALSE
175   pars[pars$name=="COV_11",'est'] <- TRUE
176   fit <- confmirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars=pars)
177   got <- coef(fit)
178   got$GroupPars <- NULL
179   round(m2$matrices$itemParam$values - sapply(got, function(l) l[1,]), 2)
180 }