Changing from 3PL to 4PL gives all new random numbers
[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(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 correct[4,] <- 1
15
16 data <- rpf.sample(500, items, correct, cov=matrix(5,1,1))
17
18 if (1) {
19   spec <- mxMatrix(name="ItemSpec", nrow=3, ncol=numItems,
20                    values=sapply(items, function(m) slot(m,'spec')),
21                    free=FALSE, byrow=TRUE)
22   
23   ip.mat <- mxMatrix(name="itemParam", nrow=4, ncol=numItems,
24                      values=c(1,0,0, 1),
25                      free=c(FALSE, TRUE, FALSE, FALSE))
26   
27   m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
28   cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=TRUE)
29
30   m2 <- mxModel(model="drmmg", ip.mat, spec, m.mat, cov.mat,
31                 mxData(observed=data, type="raw"),
32                 mxExpectationBA81(mean="mean", cov="cov",
33                   ItemSpec="ItemSpec",
34                   ItemParam="itemParam"),
35                 mxFitFunctionBA81())
36   
37   m2 <- mxOption(m2, "Analytic Gradients", 'no')
38   if (1) {
39     m2 <- mxOption(m2, "Analytic Gradients", 'yes')
40     m2 <- mxOption(m2, "Verify level", '-1')
41   }
42   m2 <- mxOption(m2, "Function precision", '1.0E-7')
43   m2 <- mxOption(m2, "Calculate Hessian", "No")
44   m2 <- mxOption(m2, "Standard Errors", "No")
45   m2 <- mxRun(m2)
46   
47   omxCheckCloseEnough(m2@output$minimum, 14130.2, 1)
48   omxCheckCloseEnough(m2@matrices$cov@values[1,1], 4.357, .01)
49   
50   #print(m2@matrices$itemParam@values)
51   #print(correct.mat)
52   got <- cor(c(m2@matrices$itemParam@values),
53              c(correct))
54   omxCheckCloseEnough(got, .994, .01)
55 }
56
57 if (0) {
58   library(mirt)
59   rdata <- sapply(data, unclass)-1
60   # for flexMIRT, write CSV
61   write.table(rdata, file="irt-drm-mg.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
62   pars <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars='values')
63   pars[pars$name=="a1",'value'] <- 1
64   pars[pars$name=="a1",'est'] <- FALSE
65   pars[pars$name=="COV_11",'est'] <- TRUE
66   fit <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars=pars)
67   # LL -7064.519 * -2 = 14129.04
68   got <- coef(fit)
69   print(got$GroupPars)
70   # COV 4.551
71   got$GroupPars <- NULL
72   round(m2@matrices$itemParam@values - simplify2array(got), 2)
73   
74   # MH-RM takes forever, not run
75   pars <- confmirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars='values')
76   pars[pars$name=="a1",'value'] <- 1
77   pars[pars$name=="a1",'est'] <- FALSE
78   pars[pars$name=="COV_11",'est'] <- TRUE
79   fit <- confmirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars=pars)
80   got <- coef(fit)
81   got$GroupPars <- NULL
82   round(m2@matrices$itemParam@values - sapply(got, function(l) l[1,]), 2)
83 }