Compute condition number of information matrix
[openmx:openmx.git] / models / nightly / ifa-fisher-info.R
1 # Here is some code to compute the expected Fisher information matrix,
2 # Equation 9 in Bock & Aitken (1981).
3 #
4 # The Appendix of Tian, Cai, Thissen, & Xin (2012) has a good discussion.
5
6 library(rpf)
7 library(OpenMx)
8 library(numDeriv)
9
10 sample.size <- 1000
11 numItems <- 5
12 spec <- list()
13 spec[1:numItems] <- rpf.grm()
14 param <- structure(c(1.76438, -1.02925, 0.676778, -0.639092, 0.832068,  1.44757, 2.01329, -1.52509, 1.00949, 0.10616),
15                    .Dim = c(2L, 5L ), .Dimnames = list(NULL, c("i1", "i2", "i3", "i4", "i5")))
16
17 if (0) {
18   data.fm <- sapply(rpf.sample(sample.size, spec, param), unclass)-1
19   write.table(data.fm, file="info.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
20 }
21
22 ip.mat <- mxMatrix(name="ItemParam", nrow=dim(param)[1], ncol=numItems,
23                    values=param, free=TRUE)
24 m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE)
25 cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=FALSE)
26
27 compute.gradient <- function(v) {
28   pat <- data.frame(t(v))
29   for (c in colnames(pat)) {
30     pat[[c]] <- ordered(pat[[c]], levels=0:1)
31   }
32   
33   m1 <- mxModel(model="latent",
34                 ip.mat, m.mat, cov.mat,
35                 mxData(observed=pat, type="raw"),
36                 mxExpectationBA81(mean="mean", cov="cov",
37                                   ItemSpec=spec,
38                                   ItemParam="ItemParam"),
39                 mxFitFunctionML(),
40                 mxComputeSequence(steps=list(
41                   mxComputeOnce('expectation', context="EM"),
42                   mxComputeOnce('fitfunction', gradient=TRUE))))
43   m1.fit <- mxRun(m1, silent=TRUE)
44   
45   ref <- m1.fit@expectation@patternLikelihood
46   
47   # The gradient is computed in log likelihood units
48   # but we need probability units. In log likelihood units,
49   # the gradient matches numDeriv to around 1e-8. After
50   # converting units, the precision goes way down to
51   # around 0.01.
52   grad <- exp(ref) - exp(ref + m1.fit@output$gradient)
53   
54   list(l=exp(ref), g=grad)
55 }
56
57 fisher.info <- function() {
58   numPat <- prod(vapply(spec, function(s) s@outcomes, 1))
59   patLik <- c()
60   grad <- matrix(NA, nrow=numPat, sum(ip.mat@free))
61   for (h in 1:numPat) {
62     index <- h - 1
63     v <- rep(NA, numItems)
64     for (i in 1:numItems) {
65       v[i] <- index %% spec[[i]]@outcomes
66       index <- index %/% spec[[i]]@outcomes
67     }
68     stuff <- compute.gradient(v)
69     patLik[h] <- stuff$l
70     grad[h,] <- stuff$g
71   }
72   list(l=patLik, g=grad)
73 }
74
75 fi <- fisher.info()
76 omxCheckCloseEnough(sum(fi$l), 1, .01)
77
78 se <- sqrt(diag(solve(sample.size * t(fi$g) %*% diag(1/fi$l) %*% fi$g)))
79 fm.se <- c(.25, .13, .11, .07, .14, .10, .32, .18, .13, .08)
80 omxCheckCloseEnough(se, fm.se, .04)