Fix Sandwich; add tests
[openmx:openmx.git] / models / nightly / ifa-bifactor.R
1 # echo 0 > /proc/self/coredump_filter  # normally 023
2 # R --vanilla --no-save -f models/failing/bock-aitkin-1981.R
3 # R -d gdb --vanilla --no-save -f models/failing/bock-aitkin-1981.R
4
5 #options(error = browser)
6 require(OpenMx)
7 require(rpf)
8 library(mvtnorm)
9
10 set.seed(5)
11
12 numItems <- 20
13 numPersons <- 1000
14 maxDim <- 2
15
16 items <- vector("list", numItems)
17 correct <- vector("list", numItems)
18 for (ix in 1:numItems) {
19         items[[ix]] <- rpf.drm(factors=maxDim)
20         correct[[ix]] <- rpf.rparam(items[[ix]])
21         correct[[ix]][[4]] <- 0   # no guessing, for now
22         correct[[ix]][[5]] <- 1   # upper cound
23 }
24 correct.mat <- simplify2array(correct)
25
26 maxParam <- max(vapply(items, rpf.numParam, 0))
27 maxOutcomes <- max(vapply(items, function(i) i@outcomes, 0))
28
29 design <- matrix(c(rep(1L,numItems),
30                    rep(2L,numItems/2), rep(3L, numItems/2)), byrow=TRUE, nrow=2)
31
32 theta <- t(rmvnorm(numPersons, mean=rnorm(3, sd=.25)))
33 data <- rpf.sample(theta, items, correct, design)
34
35 if (0) {
36   rdata <- sapply(data, unclass)-1
37   # for flexMIRT, write CSV
38   write.table(rdata, file="ifa-bifactor.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
39 }
40
41 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
42                    values=c(1.414, 1, 0, 0, 1),
43                    free=c(rep(TRUE, 3), FALSE, FALSE))
44
45 #ip.mat@values[2,1] <- correct.mat[2,1]
46 #ip.mat@free[2,1] <- FALSE
47
48 m.mat <- mxMatrix(name="mean", nrow=1, ncol=3, values=0, free=FALSE)
49 cov.mat <- mxMatrix(name="cov", nrow=3, ncol=3, values=diag(3), free=FALSE)
50
51 m1 <- mxModel(model="bifactor",
52           ip.mat, m.mat, cov.mat,
53           mxData(observed=data, type="raw"),
54           mxExpectationBA81(mean="mean", cov="cov",
55              ItemSpec=items,
56              design=design,
57              ItemParam="ItemParam",
58             qpoints=29),
59               mxFitFunctionML(),
60               mxComputeOnce('expectation', context='EM'))
61 m1 <- mxRun(m1)
62
63 omxCheckCloseEnough(sum(m1@expectation@patternLikelihood), -12629.4, .1)
64 omxCheckCloseEnough(fivenum(m1@expectation@patternLikelihood),
65                     c(-15.7575854, -14.9684791, -14.0992631, -12.3467773, -3.5902924 ), 1e-4)
66 omxCheckCloseEnough(sum(m1@expectation@em.expected), 20000, 1)
67 omxCheckCloseEnough(fivenum(m1@expectation@em.expected),
68                     c(0, 0, 1.8e-06, 0.0034365, 43.2895967), 1e-4)
69
70 m1 <- mxModel(m1,
71               mxComputeSequence(steps=list(
72                 mxComputeOnce('expectation', context='EM'),
73                 mxComputeOnce('fitfunction', fit=TRUE, gradient=TRUE, hessian=TRUE)
74               )))
75 m1 <- mxRun(m1)
76 omxCheckCloseEnough(m1@fitfunction@result, 2*11850.68, .01)
77 omxCheckCloseEnough(fivenum(m1@output$gradient), 2*c(-369.32879, -14.47296, 13.1165, 50.07066, 323.04627 ), .01)
78 omxCheckCloseEnough(fivenum(m1@output$hessian[m1@output$hessian != 0]),
79                     2*c(-53.666201, -7.6857353, -6.0121325, 89.8735155, 192.6600613 ), 1e-4)
80
81 m1 <- mxModel(m1,
82               mxData(observed=data, type="raw"),
83               mxExpectationBA81(mean="mean", cov="cov",
84                                 ItemSpec=items,
85                                 design=design,
86                                 ItemParam="ItemParam",
87                                 qpoints=29, scores="full"),
88               mxComputeEM('expectation',
89                           mxComputeNewtonRaphson(free.set='ItemParam'),
90                           mxComputeOnce('fitfunction', free.set=c("mean","cov"), fit=TRUE)
91                                  ))
92
93 m1 <- mxRun(m1, silent=TRUE)
94 #print(correct.mat)
95 #print(m1@matrices$ItemParam@values)
96 omxCheckCloseEnough(m1@output$minimum, 20859.87, .01)
97 got <- cor(c(m1@matrices$ItemParam@values), c(correct.mat))
98 omxCheckCloseEnough(got, .977, .01)
99 scores <- m1@expectation@scores.out
100 omxCheckCloseEnough(cor(c(scores[,1]), c(theta[1,])), .758, .01)
101 omxCheckCloseEnough(cor(c(scores[,2]), c(theta[2,])), .781, .01)
102 omxCheckCloseEnough(cor(c(scores[,3]), c(theta[3,])), .679, .01)
103
104 omxCheckCloseEnough(sum(abs(scores[,2] - theta[2,]) < 2*scores[,5]), 933, 5)
105 omxCheckCloseEnough(sum(abs(scores[,2] - theta[2,]) < 3*scores[,5]), 1000, 2)
106
107 i1 <- mxModel(m1,
108               mxComputeSequence(steps=list(
109                 mxComputeOnce('expectation', context="EM"),
110                 mxComputeOnce('fitfunction', information=TRUE, info.method="meat"),
111                 mxComputeStandardError(),
112                 mxComputeHessianQuality())))
113 i1 <- mxRun(i1, silent=TRUE)
114
115 #cat(deparse(round(i1@output$standardErrors,3)))
116 se <- c(0.195, 0.275, 0.129, 0.12, 0.123, 0.083, 0.336, 0.196,  0.137, 0.157,
117         0.148, 0.126, 0.262, 0.223, 0.194, 0.18, 0.215,  0.148, 0.223, 0.314,
118         0.286, 0.135, 0.135, 0.103, 0.246, 0.361,  0.139, 0.121, 0.118, 0.079,
119         0.123, 0.138, 0.095, 0.141, 0.159,  0.111, 0.161, 0.136, 0.094, 0.144, 0.155,
120         0.096, 0.141, 0.147,  0.101, 0.154, 0.219, 0.169, 0.212, 0.177, 0.172, 0.152, 0.197,
121         0.125, 0.179, 0.178, 0.107, 0.151, 0.135, 0.101)
122 omxCheckCloseEnough(i1@output$conditionNumber, 59, 1)
123 omxCheckCloseEnough(c(i1@output$standardErrors), se, .001)
124 em.meat <- i1@output$hessian
125
126 i1 <- mxModel(m1,
127               mxComputeSequence(steps=list(
128                 mxComputeOnce('expectation'),
129                 mxComputeOnce('fitfunction', information=TRUE, info.method="meat"))))
130 i1 <- mxRun(i1, silent=TRUE)
131 omxCheckCloseEnough(max(abs(i1@output$hessian - em.meat)), 0, .001)
132
133 i2 <- mxModel(m1,
134               mxComputeSequence(steps=list(
135                 mxComputeOnce('expectation'),
136                 mxComputeOnce('fitfunction', information=TRUE, info.method="sandwich"),
137                 mxComputeStandardError(),
138                 mxComputeHessianQuality())))
139 i2 <- mxRun(i2, silent=TRUE)
140 omxCheckCloseEnough(i2@output$conditionNumber, 142, 1)
141
142 swse <- c(0.247, 0.41, 0.142, 0.134, 0.125, 0.092, 0.378, 0.211,  0.174, 0.185,
143           0.192, 0.133, 0.273, 0.256, 0.203, 0.229, 0.247,  0.172, 0.29, 0.292,
144           0.267, 0.137, 0.145, 0.105, 0.337, 0.476,  0.17, 0.13, 0.126, 0.085,
145           0.141, 0.155, 0.105, 0.159, 0.207,  0.119, 0.169, 0.15, 0.109, 0.176,
146           0.18, 0.119, 0.183, 0.182,  0.125, 0.192, 0.225, 0.171, 0.209, 0.228,
147           0.175, 0.2, 0.26, 0.135,  0.214, 0.218, 0.136, 0.151, 0.148, 0.117)
148
149 #cat(deparse(round(i2@output$standardErrors,3)))
150 omxCheckCloseEnough(c(i2@output$standardErrors), swse, .001)