ifa: Add options for EAP scores
[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 }
23 correct.mat <- simplify2array(correct)
24
25 maxParam <- max(vapply(items, rpf.numParam, 0))
26 maxOutcomes <- max(vapply(items, function(i) i@outcomes, 0))
27
28 design <- matrix(c(rep(1,numItems),
29                    rep(2,numItems/2), rep(3, numItems/2)), byrow=TRUE, nrow=2)
30
31 theta <- t(rmvnorm(numPersons, mean=rnorm(3, sd=.25)))
32 data <- rpf.sample(theta, items, correct, design)
33
34 spec <- mxMatrix(name="ItemSpec", nrow=6, ncol=numItems,
35          values=sapply(items, function(m) slot(m,'spec')),
36          free=FALSE, byrow=TRUE)
37
38 design <- mxMatrix(name="Design", nrow=maxDim, ncol=numItems,
39                    values=design)
40
41 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
42                    values=c(1.414, 1, 0, 0, 1),
43                    lbound=c(1e-6, 1e-6, -1e6, 0, 0),
44                    free=c(rep(TRUE, 3), FALSE, FALSE))
45
46 #ip.mat@values[2,1] <- correct.mat[2,1]
47 #ip.mat@free[2,1] <- FALSE
48
49 m1 <- mxModel(model="bifactor",
50           spec, design,
51           ip.mat,
52           mxData(observed=data, type="raw"),
53           mxExpectationBA81(
54              ItemSpec="ItemSpec",
55              Design="Design",
56              ItemParam="ItemParam",
57             qpoints=29,
58             scores="full"),
59           mxFitFunctionBA81()
60 )
61
62 m1 <- mxOption(m1, "Analytic Gradients", 'no')
63 if (1) {
64         m1 <- mxOption(m1, "Analytic Gradients", 'yes')
65         m1 <- mxOption(m1, "Verify level", '-1')
66 }
67 m1 <- mxOption(m1, "Function precision", '1.0E-5')
68 m1 <- mxOption(m1, "Calculate Hessian", "No")
69 m1 <- mxOption(m1, "Standard Errors", "No")
70
71 m1 <- mxRun(m1, silent=TRUE)
72 #print(correct.mat)
73 #print(m1@matrices$ItemParam@values)
74 got <- cor(c(m1@matrices$ItemParam@values), c(correct.mat))
75 omxCheckCloseEnough(got, .966, .01)
76 omxCheckCloseEnough(cor(c(m1@output$ability[1,]), c(theta[1,])), .774, .01)
77 omxCheckCloseEnough(cor(c(m1@output$ability[3,]), c(theta[2,])), .657, .01)
78 omxCheckCloseEnough(cor(c(m1@output$ability[5,]), c(theta[3,])), .683, .01)
79
80 omxCheckCloseEnough(sum(abs(m1@output$ability[3,] - theta[2,]) < 2*m1@output$ability[4,]), 919)
81 omxCheckCloseEnough(sum(abs(m1@output$ability[3,] - theta[2,]) < 3*m1@output$ability[4,]), 987)