Compute condition number of information matrix
[openmx:openmx.git] / models / nightly / fm-example2-4.R
1 # flexMIRT 1.88 Example 2-4 "Graded Model Combined Calibration and Scoring Syntax with Recoding"
2
3 library(OpenMx)
4 library(rpf)
5
6 set.seed(1)
7 m2.data <- suppressWarnings(try(read.table("models/nightly/data/NCSsim.dat"), silent=TRUE))
8 if (is(m2.data, "try-error")) m2.data <- read.table("data/NCSsim.dat")
9
10 m2.spec <- list()
11 m2.spec[1:18] <- rpf.grm(outcomes=5)
12 m2.numItems <- length(m2.spec)
13
14 for (c in 1:m2.numItems) {
15   m2.data[[c]] <- mxFactor(m2.data[[c]], levels=1:m2.spec[[c]]@outcomes)
16 }
17
18 m2.maxParam <-max(sapply(m2.spec, rpf.numParam))
19
20 ip.mat <- mxMatrix(name="ItemParam", nrow=m2.maxParam, ncol=m2.numItems,
21                    values=c(1, seq(1,-1,length.out=4)), free=TRUE)
22
23 #  m2.fmfit <- read.flexmirt("~/2012/sy/fm/ms-rasch-prm.txt")
24 # cat(deparse(round(m2.fmfit$G1$param,6)))
25 #  ip.mat@values <- m2.fmfit$G1$param
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=FALSE)
29
30 m2 <- mxModel(model="m2", m.mat, cov.mat, ip.mat,
31               mxData(observed=m2.data, type="raw"),
32               mxExpectationBA81(mean="mean", cov="cov",
33                                 ItemSpec=m2.spec,
34                                 ItemParam="ItemParam"),
35               mxFitFunctionML(),
36               mxComputeEM('expectation',
37                           mxComputeNewtonRaphson(free.set='ItemParam'),
38                           mxComputeOnce('fitfunction', fit=TRUE, free.set=c("mean", "cov"))))
39 #  m2 <- mxOption(m2, "Number of Threads", 1)
40 m2 <- mxRun(m2, silent=TRUE)
41 omxCheckCloseEnough(m2@fitfunction@result, 140199.13, .01)
42
43 #print(m2@matrices$ItemParam@values - fmfit)
44 print(m2@output$backendTime)
45
46 grp <- list(spec=m2.spec,
47             param=m2@matrices$ItemParam@values,
48             free=apply(ip.mat@free, 2, sum),
49             mean=m2@matrices$mean@values, cov=m2@matrices$cov@values)
50 colnames(grp$param) <- paste("i", 1:dim(grp$param)[2], sep="")
51 #got <- chen.thissen.1997(grp, m2.data)
52 #got$std