Move BA81 expectation output/debug stuff to @output and @debug slots
[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@debug$patternLikelihood), -9742.31, .1)
37 omxCheckCloseEnough(fivenum(middle@expectation@debug$patternLikelihood),
38                     c(-34.98313, -22.50933, -19.59691, -16.79153, -6.51683), .001)
39 omxCheckCloseEnough(sum(middle@expectation@debug$em.expected), 7000, .01)
40 omxCheckCloseEnough(fivenum(middle@expectation@debug$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, 2*9399.954, .01)
50 omxCheckCloseEnough(fivenum(testDeriv@output$gradient),
51                     2*c(-14424.48407, -62.52714, -2.51876, 71.87544, 14651.19963), .01)
52 omxCheckCloseEnough(fivenum(testDeriv@output$hessian[testDeriv@output$hessian != 0]),
53                     2*c(-1038404.94356, -20.82415, 0.06734, 53.01333, 1038503.00369 ), .01)
54 omxCheckCloseEnough(max(abs(solve(testDeriv@output$hessian) - testDeriv@output$ihessian)), 0, .001)
55
56 plan <- mxComputeSequence(steps=list(mxComputeEM('expectation',
57                               mxComputeNewtonRaphson(free.set='itemParam'),
58                               mxComputeOnce('fitfunction', fit=TRUE,
59                                             free.set=c("mean","cov"))),
60                               mxComputeOnce('fitfunction', information=TRUE, info.method="meat"),
61                               mxComputeHessianQuality(),
62                               mxComputeStandardError()))
63
64 m2 <- mxModel(m2,
65               mxExpectationBA81(
66                 ItemSpec=spec, ItemParam="itemParam",
67                 mean="mean", cov="cov",
68                 qpoints=31,
69                 scores="full"),
70               plan)
71                                   
72 m2 <- mxRun(m2)
73
74 #print(m2@matrices$itemParam@values)
75 #print(correct.mat)
76 omxCheckCloseEnough(m2@fitfunction@result, 13969.747, .01)
77 got <- cor(c(m2@matrices$itemParam@values[ip.mat@free]),
78            c(correct.mat[ip.mat@free]))
79 omxCheckCloseEnough(got, .993, .01)
80 scores <- m2@expectation@output$scores
81 omxCheckCloseEnough(scores[1:5,1], c(0.81609, 0.74994, -0.83515, 0.79766, 0.16879), .001)
82 omxCheckCloseEnough(scores[1:5,2], c(0.43522, 0.44211, 0.4686, 0.43515, 0.3842), .001)
83 omxCheckCloseEnough(sum(abs(scores[,1] - ability) < 1*scores[,2])/500, .714, .01)
84 omxCheckCloseEnough(sum(abs(scores[,1] - ability) < 2*scores[,2])/500, .95, .01)
85
86 omxCheckTrue(m2@output$infoDefinite)
87 omxCheckCloseEnough(c(m2@output$conditionNumber), 658.58, 1)
88
89 #cat(deparse(round(m2@output$standardErrors,3)))
90
91 se <- c(0.152, 0.114, 0.115, 0.261, 0.163, 0.121, 0.142,  0.109, 0.113, 0.098,
92         0.185, 0.168, 0.161, 0.132, 0.139, 0.154,  0.126, 0.113, 0.215, 0.107,
93         0.111, 0.149, 0.124, 0.12, 0.127,  0.186, 0.135, 0.134, 0.134, 0.145,
94         0.147, 0.125, 0.15, 0.12,  0.12, 0.122, 0.123, 0.152, 0.164, 0.129, 0.117,
95         0.104, 0.109,  0.139, 0.103, 0.104, 0.106, 0.109, 0.116, 0.135, 0.163, 0.115 )
96   
97 omxCheckCloseEnough(c(m2@output$standardErrors), se, .01)
98
99 i2 <- mxModel(m2,
100               mxComputeSequence(steps=list(
101                 mxComputeOnce('expectation'),
102                 mxComputeOnce('fitfunction', information=TRUE, info.method="sandwich"),
103                 mxComputeStandardError(),
104                 mxComputeHessianQuality())))
105 i2 <- mxRun(i2, silent=TRUE)
106
107 omxCheckCloseEnough(i2@output$conditionNumber, 662, 1)
108 #cat(deparse(round(i2@output$standardErrors,3)))
109 swse <- c(0.143, 0.11, 0.11, 0.238, 0.149, 0.125, 0.134, 0.106,  0.108, 0.094,
110           0.169, 0.165, 0.154, 0.128, 0.13, 0.148, 0.121,  0.119, 0.199, 0.102,
111           0.106, 0.144, 0.117, 0.115, 0.122, 0.177,  0.132, 0.13, 0.129, 0.136,
112           0.141, 0.134, 0.145, 0.116, 0.112,  0.116, 0.116, 0.151, 0.162, 0.124,
113           0.115, 0.097, 0.104, 0.125,  0.098, 0.099, 0.104, 0.107, 0.111, 0.139, 0.156, 0.11)
114 omxCheckCloseEnough(c(i2@output$standardErrors), swse, .001)