Allow ComputeIterate to test maximum absolute change
[openmx:openmx.git] / models / passing / ifa-grm1.R
1 #options(error = utils::recover)
2 library(OpenMx)
3 library(rpf)
4
5 set.seed(9)
6
7 numItems <- 14
8 spec <- vector("list", numItems)
9 for (ix in 1:numItems) { spec[[ix]] <- rpf.grm(outcomes=sample(2:7, 1)) }
10 correct <- lapply(spec, rpf.rparam)
11
12 ability <- rnorm(500)
13 data <- rpf.sample(ability, spec, correct)
14
15 ip.mat <- mxMatrix(name="itemParam", nrow=max(sapply(correct, length)), ncol=numItems)
16 correct.mat <- ip.mat@values
17 for (ix in 1:numItems) {
18   len <- length(correct[[ix]])
19   ip.mat@free[1:len,ix] <- TRUE
20   ip.mat@values[1:len,ix] <- rpf.rparam(spec[[ix]])
21   correct.mat[1:len,ix] <- correct[[ix]]
22 }
23 ip.mat@values[!ip.mat@free] <- NA
24 correct.mat[!ip.mat@free] <- NA
25
26 m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
27 cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=FALSE)
28
29 m2 <- mxModel(model="grm1", ip.mat, m.mat, cov.mat,
30               mxData(observed=data, type="raw"),
31               mxExpectationBA81(ItemSpec=spec, ItemParam="itemParam",
32                 mean="mean", cov="cov", qpoints=31),
33               mxFitFunctionML(),
34               mxComputeOnce('expectation', context='EM'))
35 middle <- mxRun(m2)
36 omxCheckCloseEnough(sum(middle@expectation@patternLikelihood), -9742.31, .1)
37 omxCheckCloseEnough(fivenum(middle@expectation@patternLikelihood),
38                     c(-34.98313, -22.50933, -19.59691, -16.79153, -6.51683), .001)
39 omxCheckCloseEnough(sum(middle@expectation@em.expected), 7000, .01)
40 omxCheckCloseEnough(fivenum(middle@expectation@em.expected),
41                     c(0, 0, 0.00451, 1.61901, 80.99209), .01)
42
43 testDeriv <- mxModel(m2,
44               mxComputeSequence(steps=list(
45                                  mxComputeOnce('expectation', context='EM'),
46                                  mxComputeOnce('fitfunction', fit=TRUE, gradient=TRUE, hessian=TRUE, ihessian=TRUE)
47                                  )))
48 testDeriv <- mxRun(testDeriv)
49 omxCheckCloseEnough(testDeriv@fitfunction@result, 9399.954, .01)
50 omxCheckCloseEnough(fivenum(testDeriv@output$gradient),
51                     c(-14424.48407, -62.52714, -2.51876, 71.87544, 14651.19963), .01)
52 omxCheckCloseEnough(fivenum(testDeriv@output$hessian[testDeriv@output$hessian != 0]),
53                     c(-1038404.94356, -20.82415, 0.06734, 53.01333, 1038503.00369 ), .01)
54 omxCheckCloseEnough(solve(testDeriv@output$hessian), testDeriv@output$ihessian, 1e-2)
55
56 m2 <- mxModel(m2,
57               mxExpectationBA81(
58                 ItemSpec=spec, ItemParam="itemParam",
59                 mean="mean", cov="cov",
60                 qpoints=31,
61                 scores="full"),
62               mxComputeIterate(steps=list(
63                                  mxComputeOnce('expectation', context='EM'),
64                                  mxComputeNewtonRaphson(free.set='itemParam'),
65                                  mxComputeOnce('expectation'),
66                                  mxComputeOnce('fitfunction', fit=TRUE, free.set=c("mean","cov"))
67                                  )))
68 m2 <- mxRun(m2)
69
70 #print(m2@matrices$itemParam@values)
71 #print(correct.mat)
72 omxCheckCloseEnough(m2@fitfunction@result, 13969.747, .01)
73 got <- cor(c(m2@matrices$itemParam@values[ip.mat@free]),
74            c(correct.mat[ip.mat@free]))
75 omxCheckCloseEnough(got, .993, .01)
76 scores <- m2@expectation@scores.out
77 omxCheckCloseEnough(scores[1:5,1], c(0.81609, 0.74994, -0.83515, 0.79766, 0.16879), .001)
78 omxCheckCloseEnough(scores[1:5,2], c(0.43522, 0.44211, 0.4686, 0.43515, 0.3842), .001)
79 omxCheckCloseEnough(sum(abs(scores[,1] - ability) < 1*scores[,2])/500, .714, .01)
80 omxCheckCloseEnough(sum(abs(scores[,1] - ability) < 2*scores[,2])/500, .95, .01)