Update regression tests
[openmx:openmx.git] / models / failing / irt-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(1,numItems),
30                    rep(2,numItems/2), rep(3, 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 spec <- mxMatrix(name="ItemSpec", nrow=3, ncol=numItems,
36          values=sapply(items, function(m) slot(m,'spec')),
37          free=FALSE, byrow=TRUE)
38
39 design <- mxMatrix(name="Design", nrow=maxDim, ncol=numItems,
40                    values=design)
41
42 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
43                    values=c(1.414, 1, 0, 0, 1),
44                    lbound=c(1e-6, 1e-6, -1e6, 0, 0),
45                    free=c(rep(TRUE, 3), FALSE, FALSE))
46 ip.mat@free.group <- 'param'
47
48 #ip.mat@values[2,1] <- correct.mat[2,1]
49 #ip.mat@free[2,1] <- FALSE
50
51 eip.mat <- mxMatrix(name="EItemParam", nrow=maxParam, ncol=numItems,
52                    values=c(1.414, 1, 0, 0, 1), free=TRUE)
53
54 m.mat <- mxMatrix(name="mean", nrow=1, ncol=3, values=0, free=FALSE)
55 cov.mat <- mxMatrix(name="cov", nrow=3, ncol=3, values=diag(3), free=FALSE)
56
57 m1 <- mxModel(model="bifactor",
58           spec, design,
59           ip.mat, m.mat, cov.mat, eip.mat,
60           mxData(observed=data, type="raw"),
61           mxExpectationBA81(mean="mean", cov="cov",
62              ItemSpec="ItemSpec",
63              Design="Design",
64              EItemParam="EItemParam",
65             qpoints=29,
66             scores="full"),
67           mxFitFunctionBA81(ItemParam="ItemParam"),
68               mxComputeIterate(steps=list(
69                 mxComputeAssign(from="ItemParam", to="EItemParam"),
70                 mxComputeOnce(expectation='expectation', context='E'),
71                 mxComputeGradientDescent(free.group='param'),
72                 mxComputeOnce(expectation='expectation', context='M'),
73                 mxComputeOnce(fitfunction='fitfunction'))))
74
75         m1 <- mxOption(m1, "Analytic Gradients", 'Yes')
76         m1 <- mxOption(m1, "Verify level", '-1')
77 m1 <- mxOption(m1, "Function precision", '1.0E-5')
78
79 m1 <- mxRun(m1, silent=TRUE)
80 #print(correct.mat)
81 #print(m1@matrices$ItemParam@values)
82 got <- cor(c(m1@matrices$ItemParam@values), c(correct.mat))
83 omxCheckCloseEnough(got, .966, .01)
84 omxCheckCloseEnough(cor(c(m1@output$ability[1,]), c(theta[1,])), .747, .01)
85 omxCheckCloseEnough(cor(c(m1@output$ability[3,]), c(theta[2,])), .781, .01)
86 omxCheckCloseEnough(cor(c(m1@output$ability[5,]), c(theta[3,])), .679, .01)
87
88 omxCheckCloseEnough(sum(abs(m1@output$ability[3,] - theta[2,]) < 2*m1@output$ability[4,]), 929)
89 omxCheckCloseEnough(sum(abs(m1@output$ability[3,] - theta[2,]) < 3*m1@output$ability[4,]), 998)