Generic ComputeEM with Ramsay (1975) acceleration
[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 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
36                    values=c(1.414, 1, 0, 0, 1),
37                    free=c(rep(TRUE, 3), FALSE, FALSE))
38
39 #ip.mat@values[2,1] <- correct.mat[2,1]
40 #ip.mat@free[2,1] <- FALSE
41
42 m.mat <- mxMatrix(name="mean", nrow=1, ncol=3, values=0, free=FALSE)
43 cov.mat <- mxMatrix(name="cov", nrow=3, ncol=3, values=diag(3), free=FALSE)
44
45 m1 <- mxModel(model="bifactor",
46           ip.mat, m.mat, cov.mat,
47           mxData(observed=data, type="raw"),
48           mxExpectationBA81(mean="mean", cov="cov",
49              ItemSpec=items,
50              design=design,
51              ItemParam="ItemParam",
52             qpoints=29),
53               mxFitFunctionML(),
54               mxComputeOnce('expectation', context='EM'))
55 m1 <- mxRun(m1)
56
57 omxCheckCloseEnough(sum(m1@expectation@patternLikelihood), -12629.4, .1)
58 omxCheckCloseEnough(fivenum(m1@expectation@patternLikelihood),
59                     c(-15.7575854, -14.9684791, -14.0992631, -12.3467773, -3.5902924 ), 1e-4)
60 omxCheckCloseEnough(sum(m1@expectation@em.expected), 20000, 1)
61 omxCheckCloseEnough(fivenum(m1@expectation@em.expected),
62                     c(0, 0, 1.8e-06, 0.0034365, 43.2895967), 1e-4)
63
64 m1 <- mxModel(m1,
65               mxComputeIterate(steps=list(
66                 mxComputeOnce('expectation', context='EM'),
67                 mxComputeOnce('fitfunction', fit=TRUE, gradient=TRUE, hessian=TRUE)
68               )))
69 m1 <- mxRun(m1)
70 omxCheckCloseEnough(m1@fitfunction@result, 11850.68, .01)
71 omxCheckCloseEnough(fivenum(m1@output$gradient), c(-369.32879, -14.47296, 13.1165, 50.07066, 323.04627 ), .01)
72 omxCheckCloseEnough(fivenum(m1@output$hessian[m1@output$hessian != 0]),
73                     c(-53.666201, -7.6857353, -6.0121325, 89.8735155, 192.6600613 ), 1e-4)
74
75 m1 <- mxModel(m1,
76               mxData(observed=data, type="raw"),
77               mxExpectationBA81(mean="mean", cov="cov",
78                                 ItemSpec=items,
79                                 design=design,
80                                 ItemParam="ItemParam",
81                                 qpoints=29, scores="full"),
82               mxComputeEM('expectation',
83                           mxComputeNewtonRaphson(free.set='ItemParam'),
84                           mxComputeOnce('fitfunction', free.set=c("mean","cov"), fit=TRUE)
85                                  ))
86
87 m1 <- mxRun(m1, silent=TRUE)
88 #print(correct.mat)
89 #print(m1@matrices$ItemParam@values)
90 got <- cor(c(m1@matrices$ItemParam@values), c(correct.mat))
91 omxCheckCloseEnough(got, .977, .01)
92 scores <- m1@expectation@scores.out
93 omxCheckCloseEnough(cor(c(scores[,1]), c(theta[1,])), .758, .01)
94 omxCheckCloseEnough(cor(c(scores[,2]), c(theta[2,])), .781, .01)
95 omxCheckCloseEnough(cor(c(scores[,3]), c(theta[3,])), .679, .01)
96
97 omxCheckCloseEnough(sum(abs(scores[,2] - theta[2,]) < 2*scores[,5]), 933, 5)
98 omxCheckCloseEnough(sum(abs(scores[,2] - theta[2,]) < 3*scores[,5]), 1000, 2)