# Here is some code to compute the expected Fisher information matrix, # Equation 9 in Bock & Aitken (1981). # # The Appendix of Tian, Cai, Thissen, & Xin (2012) has a good discussion. library(rpf) library(OpenMx) library(numDeriv) mxOption(NULL, 'loglikelihoodScale', -1) sample.size <- 1000 numItems <- 5 spec <- list() spec[1:numItems] <- rpf.grm() param <- structure(c(1.76438, -1.02925, 0.676778, -0.639092, 0.832068, 1.44757, 2.01329, -1.52509, 1.00949, 0.10616), .Dim = c(2L, 5L ), .Dimnames = list(NULL, c("i1", "i2", "i3", "i4", "i5"))) if (0) { data.fm <- sapply(rpf.sample(sample.size, spec, param), unclass)-1 write.table(data.fm, file="info.csv", quote=FALSE, row.names=FALSE, col.names=FALSE) } ip.mat <- mxMatrix(name="ItemParam", nrow=dim(param)[1], ncol=numItems, values=param, free=TRUE) m.mat <- mxMatrix(name="mean", nrow=1, ncol=1, values=0, free=FALSE) cov.mat <- mxMatrix(name="cov", nrow=1, ncol=1, values=1, free=FALSE) compute.gradient <- function(v) { pat <- data.frame(t(v)) for (c in colnames(pat)) { pat[[c]] <- ordered(pat[[c]], levels=0:1) } m1 <- mxModel(model="latent", ip.mat, m.mat, cov.mat, mxData(observed=pat, type="raw"), mxExpectationBA81(mean="mean", cov="cov", debugInternal=TRUE, ItemSpec=spec, ItemParam="ItemParam"), mxFitFunctionML(), mxComputeSequence(steps=list( mxComputeOnce('expectation', 'scores'), mxComputeOnce('fitfunction', 'gradient'), mxComputeReportDeriv()))) m1.fit <- mxRun(m1, silent=TRUE) ref <- m1.fit\$expectation\$debug\$patternLikelihood # The gradient is computed in log likelihood units # but we need probability units. In log likelihood units, # the gradient matches numDeriv to around 1e-8. After # converting units, the precision goes way down to # around 0.01. grad <- exp(ref) - exp(ref + m1.fit\$output\$gradient) list(l=exp(ref), g=grad) } fisher.info <- function() { numPat <- prod(vapply(spec, function(s) s\$outcomes, 1)) patLik <- c() grad <- matrix(NA, nrow=numPat, sum(ip.mat\$free)) for (h in 1:numPat) { index <- h - 1 v <- rep(NA, numItems) for (i in 1:numItems) { v[i] <- index %% spec[[i]]\$outcomes index <- index %/% spec[[i]]\$outcomes } stuff <- compute.gradient(v) patLik[h] <- stuff\$l grad[h,] <- stuff\$g } list(l=patLik, g=grad) } fi <- fisher.info() omxCheckCloseEnough(sum(fi\$l), 1, .01) se <- sqrt(diag(solve(sample.size * t(fi\$g) %*% diag(1/fi\$l) %*% fi\$g))) fm.se <- c(.25, .13, .11, .07, .14, .10, .32, .18, .13, .08) omxCheckCloseEnough(se, fm.se, .04)