Rework specification of free param groups per Tim's suggestion
[openmx:openmx.git] / models / passing / ifa-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         eip.mat <- mxAlgebra(itemParam, name="EItemParam", fixed=TRUE)
28
29         m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
30         cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=TRUE)
31
32         m2 <- mxModel(model="drmmg", ip.mat, spec, m.mat, cov.mat, eip.mat,
33                       mxData(observed=data, type="raw"),
34                       mxExpectationBA81(mean="mean", cov="cov",
35                                         ItemSpec="ItemSpec",
36                                         EItemParam="EItemParam"),
37                       mxFitFunctionBA81(ItemParam="itemParam"),
38                                         # integrate FitFuncBA81 into FitFuncML
39                       mxComputeIterate(steps=list(
40                                          mxComputeOnce("EItemParam"),
41                                          mxComputeOnce('expectation', context='E'),
42                                         #                                  mxComputeGradientDescent(free.set='itemParam'),
43                                          mxComputeNewtonRaphson(free.set='itemParam'),
44                                          mxComputeOnce('expectation', context='M'),
45                                          mxComputeOnce('fitfunction')
46                                          )))
47         
48         if (0) {
49                 fm <- read.flexmirt("/home/joshua/irt/ifa-drm-mg/ifa-drm-mg-prm.txt")  # fix TODO
50                 cModel <- m2
51                 cModel@matrices$itemParam@values[2,] <- fm$G1$param[2,]
52                 cModel@matrices$cov@values <- fm$G1$cov
53                 cModel <- mxModel(cModel,
54                                   mxExpectationBA81(mean="mean", cov="cov",
55                                                     ItemSpec="ItemSpec",
56                                                     EItemParam="EItemParam", scores="full"),
57                                   mxComputeSequence(steps=list(
58                                                       mxComputeOnce('expectation', context='M'),
59                                                       mxComputeOnce('fitfunction'))))
60                 cModel <- mxRun(cModel)
61                 cModel@matrices$cov@values - fm$G1$cov
62                 cModel@output$minimum
63         }
64
65         if(1) {
66                 m2 <- mxOption(m2, "Analytic Gradients", 'Yes')
67                 m2 <- mxOption(m2, "Verify level", '-1')
68                 m2 <- mxOption(m2, "Function precision", '1.0E-7')
69                 m2 <- mxRun(m2)
70                 
71                 omxCheckCloseEnough(m2@fitfunction@result, 14129.94, .01)
72                 omxCheckCloseEnough(m2@matrices$cov@values[1,1], 4.377, .01)
73                 
74                                         #print(m2@matrices$itemParam@values)
75                                         #print(correct.mat)
76                 got <- cor(c(m2@matrices$itemParam@values),
77                            c(correct))
78                 omxCheckCloseEnough(got, .994, .01)
79         }
80 }
81
82 if (0) {
83   library(mirt)
84   rdata <- sapply(data, unclass)-1
85   # for flexMIRT, write CSV
86   #write.table(rdata, file="ifa-drm-mg.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
87   pars <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars='values')
88   pars[pars$name=="a1",'value'] <- 1
89   pars[pars$name=="a1",'est'] <- FALSE
90   pars[pars$name=="COV_11",'est'] <- TRUE
91   fit <- mirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars=pars)
92   # LL -7064.519 * -2 = 14129.04
93   got <- coef(fit)
94   print(got$GroupPars)
95   # COV 4.551
96   got$GroupPars <- NULL
97   round(m2@matrices$itemParam@values - simplify2array(got), 2)
98   
99   # MH-RM takes forever, not run
100   pars <- confmirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars='values')
101   pars[pars$name=="a1",'value'] <- 1
102   pars[pars$name=="a1",'est'] <- FALSE
103   pars[pars$name=="COV_11",'est'] <- TRUE
104   fit <- confmirt(rdata, 1, itemtype="2PL", D=1, quadpts=49, pars=pars)
105   got <- coef(fit)
106   got$GroupPars <- NULL
107   round(m2@matrices$itemParam@values - sapply(got, function(l) l[1,]), 2)
108 }