Pass design as a simple matrix, not MxMatrix
[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 eip.mat <- mxAlgebra(ItemParam, name="EItemParam")
43
44 m.mat <- mxMatrix(name="mean", nrow=1, ncol=3, values=0, free=FALSE)
45 cov.mat <- mxMatrix(name="cov", nrow=3, ncol=3, values=diag(3), free=FALSE)
46
47 m1 <- mxModel(model="bifactor",
48           ip.mat, m.mat, cov.mat, eip.mat,
49           mxData(observed=data, type="raw"),
50           mxExpectationBA81(mean="mean", cov="cov",
51              ItemSpec=items,
52              design=design,
53              EItemParam="EItemParam",
54             qpoints=29,
55             scores="full"),
56           mxFitFunctionBA81(ItemParam="ItemParam"),
57               mxComputeIterate(steps=list(
58                 mxComputeOnce("EItemParam"),
59                 mxComputeOnce('expectation', context='EM'),
60                            mxComputeNewtonRaphson(free.set='ItemParam'),
61 #                mxComputeGradientDescent(free.set='ItemParam'),
62                 mxComputeOnce('expectation'),
63                         mxComputeGradientDescent(start=TRUE, useGradient=TRUE,
64                                                  free.set=c('mean','cov'))
65 #                mxComputeOnce(start=TRUE, 'fitfunction')
66                                  )))
67
68         m1 <- mxOption(m1, "Analytic Gradients", 'Yes')
69         m1 <- mxOption(m1, "Verify level", '-1')
70 m1 <- mxOption(m1, "Function precision", '1.0E-5')
71
72 m1 <- mxRun(m1, silent=TRUE)
73 #print(correct.mat)
74 #print(m1@matrices$ItemParam@values)
75 got <- cor(c(m1@matrices$ItemParam@values), c(correct.mat))
76 omxCheckCloseEnough(got, .977, .01)
77 scores <- m1@expectation@scores.out
78 omxCheckCloseEnough(cor(c(scores[,1]), c(theta[1,])), .758, .01)
79 omxCheckCloseEnough(cor(c(scores[,2]), c(theta[2,])), .781, .01)
80 omxCheckCloseEnough(cor(c(scores[,3]), c(theta[3,])), .679, .01)
81
82 omxCheckCloseEnough(sum(abs(scores[,2] - theta[2,]) < 2*scores[,5]), 933, 5)
83 omxCheckCloseEnough(sum(abs(scores[,2] - theta[2,]) < 3*scores[,5]), 1000, 2)