ifa: Snapshot
[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.gpcm(3)
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 <- t(simplify2array(correct))
26
27 good.data <- rpf.sample(250, items, correct)
28 data <- mcar(good.data, 1/3)
29 #head(data)
30
31 spec <- mxMatrix(name="ItemSpec", nrow=3, ncol=numItems,
32          values=sapply(items, function(m) slot(m,'spec')),
33          free=FALSE, byrow=TRUE)
34
35 ip.mat <- mxMatrix(name="itemParam", nrow=3, ncol=numItems,
36                    values=c(1,0,0),
37                    lbound=c(1e-6, -1e6, -1e6),
38                    free=TRUE)
39
40 m2 <- mxModel(model="test3", ip.mat, spec,
41               mxData(observed=data, type="raw"),
42               mxExpectationBA81(
43                 ItemSpec="ItemSpec",
44                 ItemParam="itemParam"),
45               mxFitFunctionBA81())
46 m2 <- mxOption(m2, "Analytic Gradients", 0)
47 if (1) {
48         m2 <- mxOption(m2, "Analytic Gradients", 'yes')
49         m2 <- mxOption(m2, "Verify level", '-1')
50 }
51 m2 <- mxOption(m2, "Calculate Hessian", "No")
52 m2 <- mxOption(m2, "Standard Errors", "No")
53 m2 <- mxOption(m2, "Function precision", '1.0E-5')
54 m2 <- mxRun(m2)
55
56 got <- cor(c(m2@matrices$itemParam@values),
57            c(t(correct.mat)))
58 omxCheckCloseEnough(got, .894, .01)