Attachment | Size |
---|---|
Screen Shot 2014-09-30 at 12.48.51 PM.png | 44.35 KB |
Hi all,
I've run a model and am a bit puzzled that it won't give me the information criteria (AIC and BIC). Everything else is computed, including chi, k, and df, which are the sufficient statistics. Am I missing something? I've included a screenshot. It seems the AIC should be
1557.625 + ncol(d)(ncol(d)-1)-2453 = 1707.625,
where ncol(d) = 33
That is a bit mysterious. Can you provide an example in R code?
Sure. Hopefully this is intelligible. I'll try to comment as I go:
paths is a vector of strings that name variables, as is cors (e.g., paths = c("var1", "var2", "var2", "var3"))
paths = matrix(paths, ncol=2, byrow=T)
cors = matrix(cors, ncol=2, byrow=T)
f.use (and factor.mat) is a matrix where the first column lists the "to" variables, and the second column lists the "from"
variables. f.use just subsets depending on which factors I choose to model.
f.use = factor.mat[which(factors%in%which.factors),]
manifestVars=f.use[,1]
latentVars=which.factors
factor loadings
mxFLoadings = mxPath(from=f.use[,2], to=f.use[,1], arrows=1, free=T)
paths
mxPaths = mxPath(from=paths[,1], paths[,2], arrows=1, free=T)
correlations
mxCors = mxPath(from=cors[,1], cors[,2], arrows=2, free=T)
residual variances
mxResids_man = mxPath(from=manifestVars, to=manifestVars, arrows=2, free=T, values=.7)
mxResids_lat = mxPath(from=latentVars, to=latentVars, arrows=2, free=F, values=1)
dataset
dat = mxData(data.matrix(d[manifestVars, manifestVars]),type="cor", numObs=80)
put together
mod = mxModel("Test Model", type="RAM",
manifestVars=manifestVars, latentVars=latentVars,
mxFLoadings, mxPaths, mxCors, mxResids_man, mxResids_lat,
dat)
below, I added the "calculate hessian" and "standard errors" arguments because I thought it might help, but it didn't
mod = mxOption(mod, "Calculate Hessian", "Yes")
mod = mxOption(mod, "Standard Errors", "Yes")
fit = mxRun(mod, intervals=T)
I don't care about comments. I just need code that works. You wrote "paths = matrix(paths, ncol=2, byrow=T)" but this fails because paths is not defined yet.
Yeah...the dataset had PHI, but I've included fake data (and code) below. I'm not getting AIC/BIC from the following:
generate a random correlation matrix (requires fifer package)
install_github("fifer", "dustinfife")
require(fifer)
set seed to make it reproducible.
set.seed(1)
d = data.frame(random.correlation(12))
names(d) = c(intersperse(c("IL_4_", "IL_5_", "IL_13_"), 1:3), "DNA1", "DNA2", "DNA3") ### intersperse requires fifer package
row.names(d) = names(d)
create factor matrix
factor.mat = cbind(names(d)[1:9], c("Th1", "Th1", "Th1", "Th2", "Th2", "Th2", "Th3", "Th3", "Th3"))
name manifest/latent variables
manifestVars = names(d)
latentVars = unique(factor.mat[,2])
factor loadings
mxFLoadings = mxPath(from=factor.mat[,2], to=factor.mat[,1], arrows=1, free=T)
paths
mxThPaths = mxPath(from=c("Th1", "Th2"), to=c("Th2", "Th3"), arrows=1, free=T)
mxAbPaths = mxPath(from=c("DNA1", "DNA2"), to=c("DNA2", "DNA3"), arrows=1, free=T)
mxCause = mxPath(from=c("DNA1", "DNA2"), to=c("Th2", "Th3"), arrows=1, free=T)
correlations
mxCors = mxPath(from=c("Th1", "Th2", "Th3"), to=c("DNA1", "DNA2", "DNA3"), arrows=2, free=T)
residual variances
mxResids_man = mxPath(from=manifestVars, to=manifestVars, arrows=2, free=T, values=.7)
mxResids_lat = mxPath(from=latentVars, to=latentVars, arrows=2, free=c(F,T,T), values=1)
dataset
dat = mxData(data.matrix(make.symmetric(d)),type="cor", numObs=80)
run the model
mod = mxModel("Test Model", type="RAM",
manifestVars=manifestVars, latentVars=latentVars,
mxFLoadings, mxThPaths, mxAbPaths, mxCause, mxCors, mxResids_man, mxResids_lat, mxResidCors,
dat)
mod = mxOption(mod, "Calculate Hessian", "Yes")
mod = mxOption(mod, "Standard Errors", "Yes")
fit = mxRun(mod, intervals=T)
summary(fit)
Are we suppose to treat 'cov' and 'cor' the same way?
I'm not sure what you mean, but changing it to "cov" does fix it (after you remove the "mxResidCors" bit in the mxModel function...sorry about that). So why would cov estimate BIC/AIC but not cor?
This is an apparent bug in
summary
for type='cor' data. Not many people use cor data, so we never caught it. The bug is now fixed in trunk, and the fix will be included in the next Beta release. In the meantime, I guess use 'cov' data. Of course note that the degrees of freedom penalty for AIC and BIC will differ across type='cor' and type='cov' because their degrees of freedom are not the same. The parameters penalty versions of these are unaffected by cov vs cor.Thanks Mike! This is monumental--I believe this is the first time I've reported a "bug" that was actually a bug :)