ifa: Use R_GetCCallable to access libirt-rpf
[openmx:openmx.git] / models / failing / irt-bifactor.R
1 # echo 0 > /proc/self/coredump_filter  # normally 023
2 # R --vanilla --no-save -f models/failing/bock-aitkin-1981.R
3 # R -d gdb --vanilla --no-save -f models/failing/bock-aitkin-1981.R
4
5 #options(error = browser)
6 require(OpenMx)
7 require(rpf)
8
9 set.seed(5)
10
11 numItems <- 6
12 numPersons <- 1000
13 maxDim <- 2
14
15 items <- vector("list", numItems)
16 correct <- vector("list", numItems)
17 for (ix in 1:numItems) {
18         items[[ix]] <- rpf.drm(dimensions=maxDim)
19         correct[[ix]] <- rpf.rparam(items[[ix]])
20         correct[[ix]][[4]] <- 0   # no guessing, for now
21 }
22 correct.mat <- simplify2array(correct)
23
24 maxParam <- max(vapply(items, function(i) i@numParam, 0))
25 maxOutcomes <- max(vapply(items, function(i) i@numOutcomes, 0))
26
27 design <- matrix(c(rep(1,numItems),
28                    2,2,2,3,3,3), byrow=TRUE, nrow=2)
29
30 data <- rpf.sample(numPersons, items, correct, design)
31
32 spec <- mxMatrix(name="ItemSpec", nrow=6, ncol=numItems,
33          values=sapply(items, function(m) slot(m,'spec')),
34          free=FALSE, byrow=TRUE)
35
36 design <- mxMatrix(name="Design", nrow=maxDim, ncol=numItems,
37                    values=design)
38
39 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
40                    values=c(
41                      rep(1.414,numItems),
42                      rep(1,numItems),
43                      rep(0,numItems),
44                      rep(0,numItems)),
45          free=c(rep(TRUE, numItems*3),
46                  rep(FALSE, numItems*1)),
47          byrow=TRUE)
48
49 #ip.mat@values[2,1] <- correct.mat[2,1]
50 #ip.mat@free[2,1] <- FALSE
51
52 m1 <- mxModel(model="bifactor",
53           spec, design,
54           ip.mat,
55           mxData(observed=data, type="raw"),
56           mxExpectationBA81(
57              ItemSpec="ItemSpec",
58              Design="Design",
59              ItemParam="ItemParam"),
60           mxFitFunctionBA81()
61 )
62
63 m1 <- mxOption(m1, "Analytic Gradients", 'no')
64 if (1) {
65         m1 <- mxOption(m1, "Analytic Gradients", 'yes')
66         m1 <- mxOption(m1, "Verify level", '-1')
67 }
68 m1 <- mxOption(m1, "Function precision", '1.0E-5')
69 m1 <- mxOption(m1, "Calculate Hessian", "No")
70 m1 <- mxOption(m1, "Standard Errors", "No")
71
72 m1 <- mxRun(m1, silent=TRUE)
73 print(correct.mat)
74 print(m1@matrices$ItemParam@values)
75 got <- cor(c(m1@matrices$ItemParam@values), c(correct.mat))
76 omxCheckCloseEnough(got, .957, .01) # caution, not converged yet