ifa: Store latent mean & cov in regular MxMatrix
[openmx:openmx.git] / models / failing / irt-guess.R
1 library(OpenMx)
2 library(rpf)
3
4 set.seed(9)
5
6 numItems <- 8
7 i1 <- rpf.drm(numChoices=2)
8 i2 <- rpf.drm(numChoices=2, factors=2)
9 items <- vector("list", numItems)
10 for (ix in seq(1,numItems,2)) items[[ix]] <- i1
11 for (ix in seq(2,numItems,2)) items[[ix]] <- i2
12 correct <- lapply(items, rpf.rparam)
13 correct.mat <- lapply(correct, function(p) if (length(p)==4) c(p,NA) else p)
14 correct.mat <- simplify2array(correct.mat)
15
16 numPeople <- 1000
17 ability <- rbind(rnorm(numPeople), rnorm(numPeople))
18 data <- rpf.sample(ability, items, correct)
19
20 spec <- mxMatrix(name="ItemSpec", nrow=6, ncol=numItems,
21          values=sapply(items, function(m) slot(m,'spec')),
22          free=FALSE, byrow=TRUE)
23
24 ip.mat <- mxMatrix(name="itemParam", nrow=5, ncol=numItems,
25                    values=c(1,0,.5,1,NA,
26                             1,1.4,0,.5,1),
27                    free=c(rep(TRUE,3), FALSE, FALSE,
28                           rep(TRUE,4), FALSE),
29                    lbound=c(1e-6, -1e6, 1e-6, 0, NA,
30                      1e-6, 1e-6, -1e6, 1e-6, 0),
31                    ubound=c(1e3, 1e6, 1-1e-3, 1, NA,
32                      1e3, 1e3, 1e6, 1-1e-3, 1))
33
34 m2 <- mxModel(model="guess", ip.mat, spec,
35               mxData(observed=data, type="raw"),
36               mxExpectationBA81(
37                 ItemSpec="ItemSpec",
38                 ItemParam="itemParam",
39                 qpoints=19),  # use less points since accuracy is bad anyway
40               mxFitFunctionBA81())
41
42 m2 <- mxOption(m2, "Analytic Gradients", 'no')
43 if (1) {
44         m2 <- mxOption(m2, "Analytic Gradients", 'yes')
45         m2 <- mxOption(m2, "Verify level", '-1')
46 }
47 m2 <- mxOption(m2, "Function precision", '1.0E-5')
48 m2 <- mxOption(m2, "Calculate Hessian", "No")
49 m2 <- mxOption(m2, "Standard Errors", "No")
50 m2 <- mxRun(m2)
51
52 #print(m2@matrices$itemParam@values)
53 #print(correct.mat)
54 ipv <- m2@matrices$itemParam@values
55 got <- cor(c(ipv[!is.na(ipv)]),
56            c(correct.mat[!is.na(correct.mat)]))
57 omxCheckCloseEnough(got, .492, .01)   # removed prior so accuracy is bad
58 #ability <- scale(ability)
59 #omxCheckCloseEnough(m2@output$ability[,1], as.vector(ability), 4*max(m2@output$ability[,2]))
60 #omxCheckCloseEnough(cor(c(m2@output$ability[,1]), ability), .80, .01)