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