ifa: Store latent mean & cov in regular MxMatrix
[openmx:openmx.git] / models / failing / irt-drm-mg.R
1 library(OpenMx)
2 library(rpf)
3
4 set.seed(9)
5
6 numItems <- 30
7 i1 <- rpf.drm(numChoices=4, multidimensional=TRUE)
8 items <- list()
9 items[1:numItems] <- i1
10 correct <- matrix(NA, 4, numItems)
11 for (x in 1:numItems) correct[,x] <- rpf.rparam(i1)
12 correct[1,] <- 1
13 correct[3,] <- 0
14
15 data <- rpf.sample(500, items, correct, cov=matrix(5,1,1))
16
17 if (1) {
18   spec <- mxMatrix(name="ItemSpec", nrow=6, ncol=numItems,
19                    values=sapply(items, function(m) slot(m,'spec')),
20                    free=FALSE, byrow=TRUE)
21   
22   ip.mat <- mxMatrix(name="itemParam", nrow=4, ncol=numItems,
23                      values=c(1,0,0, 1),
24                      free=c(FALSE, TRUE, FALSE, FALSE))
25   
26   m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
27   cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=TRUE)
28
29   m2 <- mxModel(model="drmmg", ip.mat, spec, m.mat, cov.mat,
30                 mxData(observed=data, type="raw"),
31                 mxExpectationBA81(mean="mean", cov="cov",
32                   ItemSpec="ItemSpec",
33                   ItemParam="itemParam"),
34                 mxFitFunctionBA81())
35   
36   m2 <- mxOption(m2, "Analytic Gradients", 'no')
37   if (1) {
38     m2 <- mxOption(m2, "Analytic Gradients", 'yes')
39     m2 <- mxOption(m2, "Verify level", '-1')
40   }
41   m2 <- mxOption(m2, "Function precision", '1.0E-7')
42   m2 <- mxOption(m2, "Calculate Hessian", "No")
43   m2 <- mxOption(m2, "Standard Errors", "No")
44   m2 <- mxRun(m2)
45   
46   omxCheckCloseEnough(m2@output$minimum, 14097, 1)
47   omxCheckCloseEnough(m2@matrices$cov@values[1,1], 4.518, .01)
48   
49   #print(m2@matrices$itemParam@values)
50   #print(correct.mat)
51   got <- cor(c(m2@matrices$itemParam@values),
52              c(correct))
53   omxCheckCloseEnough(got, .996, .01)
54 }
55
56 if (0) {
57   library(mirt)
58   rdata <- sapply(data, unclass)-1
59   # for flexMIRT, write CSV
60   write.table(rdata, file="irt-drm-mg.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
61   pars <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars='values')
62   pars[pars$name=="a1",'value'] <- 1
63   pars[pars$name=="a1",'est'] <- FALSE
64   pars[pars$name=="COV_11",'est'] <- TRUE
65   fit <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars=pars)
66   # LL -6905.250 * -2 = 13810.5
67   got <- coef(fit)
68   print(got$GroupPars)
69   # COV 4.913
70   got$GroupPars <- NULL
71   round(m2@matrices$itemParam@values - simplify2array(got), 2)
72   
73   # MH-RM takes forever, not run
74   pars <- confmirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars='values')
75   pars[pars$name=="a1",'value'] <- 1
76   pars[pars$name=="a1",'est'] <- FALSE
77   pars[pars$name=="COV_11",'est'] <- TRUE
78   fit <- confmirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars=pars)
79   got <- coef(fit)
80   got$GroupPars <- NULL
81   round(m2@matrices$itemParam@values - sapply(got, function(l) l[1,]), 2)
82 }