Update regression tests
[openmx:openmx.git] / models / failing / irt-missingdata.R
1 #options(error = utils::recover)
2 library(OpenMx)
3 library(rpf)
4
5 mcar <- function(data, pct) {
6         size <- prod(dim(data))
7         erase <- rep(TRUE, size * pct)
8         mask <- c(erase, rep(FALSE, size - length(erase)))[order(runif(size))]
9         shaped.mask <- array(mask, dim=dim(data))
10         data[shaped.mask] <- NA
11         data <- data[apply(is.na(data), 1, sum) != dim(data)[2],]  # remove when all items are missing
12         data
13 }
14
15 set.seed(8)
16
17 numItems <- 5
18 i1 <- rpf.nrm(3, T.c=diag(2))
19 items <- vector("list", numItems)
20 correct <- vector("list", numItems)
21 for (ix in 1:numItems) {
22   items[[ix]] <- i1
23   correct[[ix]] <- rpf.rparam(i1)
24 }
25 correct.mat <- simplify2array(correct)
26 correct.mat[2,] <- 1
27 correct.mat[3,] <- 0
28
29 good.data <- rpf.sample(250, items, correct.mat)
30 data <- mcar(good.data, 1/3)
31 #head(data)
32
33 spec <- mxMatrix(name="ItemSpec", nrow=19, ncol=numItems,
34          values=sapply(items, function(m) slot(m,'spec')),
35          free=FALSE, byrow=TRUE)
36
37 ip.mat <- mxMatrix(name="itemParam", nrow=5, ncol=numItems,
38                    values=c(1,1,0,0,0),
39                    lbound=c(1e-6, rep(NA,4)),
40                    free=c(TRUE,FALSE,FALSE,TRUE,TRUE))
41 ip.mat@free.group <- 'param'
42
43 eip.mat <- mxMatrix(name="EItemParam", nrow=5, ncol=numItems,
44                    values=c(1,1,0,0,0), free=TRUE)
45
46 m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
47 cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=FALSE)
48
49 m2 <- mxModel(model="test3", ip.mat, spec, m.mat, cov.mat, eip.mat,
50               mxData(observed=data, type="raw"),
51               mxExpectationBA81(mean="mean", cov="cov",
52                                 ItemSpec="ItemSpec",
53                                 EItemParam="EItemParam"),
54               mxFitFunctionBA81(ItemParam="itemParam"),
55               mxComputeIterate(steps=list(
56                 mxComputeAssign(from="itemParam", to="EItemParam"),
57                 mxComputeOnce(expectation='expectation', context='E'),
58                 mxComputeGradientDescent(free.group='param'),
59                 mxComputeOnce(expectation='expectation', context='M'),
60                 mxComputeOnce(fitfunction='fitfunction'))))
61
62         m2 <- mxOption(m2, "Analytic Gradients", 'Yes')
63         m2 <- mxOption(m2, "Verify level", '-1')
64 m2 <- mxOption(m2, "Function precision", '1.0E-5')
65 m2 <- mxRun(m2)
66
67 got <- cor(c(m2@matrices$itemParam@values[c(1,4,5),]),
68            c(correct.mat[c(1,4,5),]))
69 omxCheckCloseEnough(got, .936, .01)