Add sanity check for minItemsPerScore
[openmx:openmx.git] / models / passing / ifa-allna.R
1 #options(error = browser)
2 require(OpenMx)
3 require(rpf)
4
5 numItems <- 4
6 maxDim <- 2
7
8 items <- list()
9 items[1:numItems] <- rpf.grm(factors=maxDim)
10 correct.mat <- sapply(items, rpf.rparam)
11
12 maxParam <- max(vapply(items, rpf.numParam, 0))
13 maxOutcomes <- max(vapply(items, function(i) i$outcomes, 0))
14
15 design <- matrix(c(rep(1L,numItems),
16                    rep(2L,numItems/2), rep(3L, numItems/2)), byrow=TRUE, nrow=2)
17
18 data <- rpf.sample(4, items, correct.mat, design)
19
20 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
21                    values=c(1.414, 1, 0), free=TRUE)
22 colnames(ip.mat) <- colnames(data)
23 m.mat <- mxMatrix(name="mean", nrow=1, ncol=3, values=0, free=FALSE)
24 colnames(m.mat) <- paste("f", 1:3, sep="")
25 cov.mat <- mxMatrix(name="cov", nrow=3, ncol=3, values=diag(3), free=FALSE)
26 dimnames(cov.mat) <- list(paste("f", 1:3, sep=""), paste("f", 1:3, sep=""))
27
28 mkmodel <- function(data) {
29   mxModel(model="bifactor",
30           ip.mat, m.mat, cov.mat,
31           mxData(observed=data, type="raw"),
32           mxExpectationBA81(mean="mean", cov="cov", debugInternal=TRUE,
33                             ItemSpec=items, design=design, ItemParam="ItemParam", qpoints=29),
34           mxFitFunctionML(),
35           mxComputeOnce('expectation', 'scores'))
36 }
37
38 tdata <- data
39 tdata[1,] <- NA
40 omxCheckError(mxRun(mkmodel(tdata)), "1:Data row 1 has no information about ability 1
41 2:Data row 1 has no information about ability 2
42 3:Data row 1 has no information about ability 3")
43
44 tdata <- data
45 tdata[2,1:2] <- NA
46 omxCheckError(mxRun(mkmodel(tdata)), "Data row 2 has no information about ability 2")
47
48 tdata <- data
49 tdata[3,3:4] <- NA
50 omxCheckError(mxRun(mkmodel(tdata)), "Data row 3 has no information about ability 3")
51
52 #-------------------------------------------------
53
54 mcar <- function(ret, mcar) {
55   size <- prod(dim(ret))
56   mask <- rep(FALSE, size)
57   mask[sample.int(size, size * mcar)] <- TRUE
58   shaped.mask <- array(mask, dim=dim(ret))
59   ret[shaped.mask] <- NA
60   ret
61 }
62
63 set.seed(1)
64
65 numItems <- 12
66 items <- list()
67 items[1:numItems] <- rpf.grm()
68 correct.mat <- sapply(items, rpf.rparam)
69
70 maxParam <- max(vapply(items, rpf.numParam, 0))
71
72 slicen <- 50
73 data <- rpf.sample(slicen, items, correct.mat)
74 for (m in seq(.1, .9, .05)) {
75   data <- rbind(data, mcar(rpf.sample(slicen, items, correct.mat), m))
76 }
77
78 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems, values=correct.mat)
79 colnames(ip.mat) <- colnames(data)
80
81 m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0)
82 colnames(m.mat) <- 'f1'
83 cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, dimnames=list('f1','f1'))
84
85 result <- expand.grid(ips=0:numItems, v=NA)
86 for (r in 1:nrow(result)) {
87   m1 <- mxModel(model="perScore", ip.mat, m.mat, cov.mat,
88                 mxData(observed=data, type="raw"),
89                 mxExpectationBA81(mean="mean", cov="cov", ItemSpec=items, ItemParam="ItemParam",
90                                   naAction="pass", scores="full", minItemsPerScore=result$ips[r]),
91                 mxComputeOnce('expectation'))
92   fit1 <- mxRun(m1, silent=TRUE)
93   v <- var(fit1$expectation$output$scores[,1], na.rm=TRUE)
94   result$v[r] <- v
95 }
96
97 c1 <- coef(lm(v ~ ips, result))
98 #cat(deparse(round(c1,4)))
99 omxCheckCloseEnough(c1, c(0.5763, 0.0144), .001)
100
101 # ------------------------------
102
103 m1 <- mxModel(model="perScore", ip.mat, m.mat, cov.mat,
104               mxData(observed=data, type="raw"),
105               mxExpectationBA81(mean="mean", cov="cov", ItemSpec=items, ItemParam="ItemParam",
106                                 naAction="pass", scores="full", minItemsPerScore=as.integer(numItems+1)),
107               mxComputeOnce('expectation'))
108 omxCheckError(mxRun(m1, silent=TRUE),
109               "perScore.expectation: minItemsPerScore (=13) cannot be larger than the number of items (=12)")