ifa: Add options for EAP scores
[openmx:openmx.git] / models / failing / irt-2d.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(7)
11
12 numItems <- 5
13 numPeople <- 500
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=min(ix,maxDim),
20                                multidimensional=TRUE)
21         correct[[ix]] <- rpf.rparam(items[[ix]])
22         if (ix>1) correct[[ix]][[4]] <- 0   # no guessing, for now
23 }
24 correct[[1]][[3]] <- 0   # no guessing, for now
25 correct[[1]][[5]] <- 1
26 correct.mat <- simplify2array(correct)
27
28 maxParam <- max(vapply(items, function(i) rpf.numParam(i), 0))
29 maxOutcomes <- max(vapply(items, function(i) i@outcomes, 0))
30
31 design <- matrix(c(1, 1,1,1,2,
32                   NA,2,2,2,1), byrow=TRUE, nrow=2)
33
34 ability <- array(rnorm(numPeople * 2), dim=c(2, numPeople))
35 #cov <- matrix(c(1, .68, .68, 1), nrow=2)
36 #ability <- rmvnorm(numPeople, sigma=cov)
37 data <- rpf.sample(ability, items, correct, design)
38
39 spec <- mxMatrix(name="ItemSpec", nrow=6, ncol=numItems,
40          values=sapply(items, function(m) slot(m,'spec')),
41          free=FALSE, byrow=TRUE)
42
43 design <- mxMatrix(name="Design", nrow=maxDim, ncol=numItems,
44                   values=design)
45
46 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
47                    values=c(1, 1.4, 0, 0, 1),
48                   lbound=c(1e-6, -1e6, 0,  0, 0,
49                     rep(c(1e-6, 1e-6, -1e6, 0, 0), 4)),
50          free=c(TRUE, TRUE, FALSE, FALSE, FALSE,
51           rep(c(TRUE, TRUE, TRUE, FALSE, FALSE), 4)))
52 ip.mat@values[4,1] <- 1
53
54 m1 <- mxModel(model="2dim",
55           spec, design,
56           ip.mat,
57           mxData(observed=data, type="raw"),
58           mxExpectationBA81(
59              ItemSpec="ItemSpec",
60              Design="Design",
61              ItemParam="ItemParam",
62             qpoints=29,
63             scores="full"),
64           mxFitFunctionBA81())
65
66 m1 <- mxOption(m1, "Analytic Gradients", 'no')
67 if (1) {
68         m1 <- mxOption(m1, "Analytic Gradients", 'yes')
69         m1 <- mxOption(m1, "Verify level", '-1')
70 }
71 m1 <- mxOption(m1, "Function precision", '1.0E-5')
72 m1 <- mxOption(m1, "Calculate Hessian", "No")
73 m1 <- mxOption(m1, "Standard Errors", "No")
74
75 m1 <- mxRun(m1, silent=TRUE)
76
77 print(m1@output$latent.mean)
78 print(m1@output$latent.cov)
79                                         #print(m1@matrices$ItemParam@values)
80                                         #print(correct.mat)
81 omxCheckCloseEnough(cor(c(m1@matrices$ItemParam@values),
82                         c(correct.mat)), .916, .01)
83 max.se <- max(m1@output$ability[c(2,4),])
84 omxCheckCloseEnough(m1@output$ability[c(1,3),], ability, max.se*2.5)
85 omxCheckCloseEnough(.57, cor(c(m1@output$ability[c(1,3),]), c(ability)), .01)