Pass item models in a list instead of an MxMatrix
[openmx:openmx.git] / models / passing / ifa-drm-mg2.R
1 library(OpenMx)
2 library(rpf)
3
4 set.seed(9)
5
6 numItems <- 20
7 i1 <- rpf.grm(outcomes=2)
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
13 data.g1 <- rpf.sample(500, items, correct)
14 data.g2 <- rpf.sample(500, items, correct, mean=-1, cov=matrix(5,1,1))
15 data.g3 <- rpf.sample(500, items, correct, mean=1, cov=matrix(.5,1,1))
16
17 if (0) {
18   # for flexMIRT, write CSV
19   write.table(sapply(data.g1, unclass)-1, file="drm-mg2-g1.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
20   write.table(sapply(data.g2, unclass)-1, file="drm-mg2-g2.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
21   write.table(sapply(data.g3, unclass)-1, file="drm-mg2-g3.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
22 }
23
24 mkgroup <- function(model.name, data, latent.free) {
25   ip.mat <- mxMatrix(name="ItemParam", nrow=2, ncol=numItems,
26                      values=c(1,0), free=TRUE)
27   
28   for (ix in 1:numItems) {
29     for (px in 1:2) {
30       name <- paste(c('p',ix,',',px), collapse='')
31       ip.mat@labels[px,ix] <- name
32     }
33   }
34   
35   eip.mat <- mxAlgebra(ItemParam, name="EItemParam", fixed=TRUE)
36   
37   dims <- 1
38   m.mat <- mxMatrix(name="mean", nrow=1, ncol=dims, values=0, free=latent.free)
39   cov.mat <- mxMatrix(name="cov", nrow=dims, ncol=dims, values=diag(dims),
40                       free=latent.free)
41   
42   m1 <- mxModel(model=model.name, ip.mat, eip.mat, m.mat, cov.mat,
43                 mxData(observed=data, type="raw"),
44                 mxExpectationBA81(
45                   ItemSpec=items,
46                   EItemParam="EItemParam",
47                   mean="mean", cov="cov"),
48                 mxFitFunctionBA81(ItemParam="ItemParam", rescale=FALSE))
49   m1
50 }
51
52 g1 <- mkgroup("g1", data.g1, FALSE)
53 g2 <- mkgroup("g2", data.g2, TRUE)
54 g3 <- mkgroup("g3", data.g3, TRUE)
55
56 groups <- paste("g", 1:3, sep="")
57
58 # load flexmirt fit and compare TODO
59
60 if (1) {
61   grpModel <- mxModel(model="groupModel", g1, g2, g3,
62                       mxFitFunctionMultigroup(paste(groups, "fitfunction", sep=".")),
63                       mxComputeIterate(steps=list(
64                         mxComputeOnce(paste(groups, "EItemParam", sep=".")),
65                         mxComputeOnce(paste(groups, 'expectation', sep='.'), context='EM'),
66                         mxComputeNewtonRaphson(free.set=paste(groups,'ItemParam',sep=".")),
67                         mxComputeOnce(paste(groups, 'expectation', sep=".")),
68                         mxComputeOnce('fitfunction')
69                       )))
70   
71   #grpModel <- mxOption(grpModel, "Number of Threads", 1)
72   
73   # NPSOL options:
74   grpModel <- mxOption(grpModel, "Analytic Gradients", 'Yes')
75   grpModel <- mxOption(grpModel, "Verify level", '-1')
76   grpModel <- mxOption(grpModel, "Function precision", '1.0E-7')
77   
78   grpModel <- mxRun(grpModel)
79   #grpModel@output$minimum TODO
80   omxCheckCloseEnough(grpModel@submodels$g2@matrices$mean@values, -.834, .01)
81   omxCheckCloseEnough(grpModel@submodels$g2@matrices$cov@values, 3.93, .01)
82   omxCheckCloseEnough(grpModel@submodels$g3@matrices$mean@values, .933, .01)
83   omxCheckCloseEnough(grpModel@submodels$g3@matrices$cov@values, .444, .01)
84 }
85
86 if (0) {
87   library(mirt)
88   rdata <- sapply(data, unclass)-1
89   # for flexMIRT, write CSV
90   #write.table(rdata, file="ifa-drm-mg.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
91   pars <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars='values')
92   pars[pars$name=="a1",'value'] <- 1
93   pars[pars$name=="a1",'est'] <- FALSE
94   pars[pars$name=="COV_11",'est'] <- TRUE
95   fit <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars=pars)
96   # LL -7064.519 * -2 = 14129.04
97   got <- coef(fit)
98   print(got$GroupPars)
99   # COV 4.551
100   got$GroupPars <- NULL
101   round(m2@matrices$itemParam@values - simplify2array(got), 2)
102 }