---
title: "data"
date: "08-Mar-2019"
output: html_document
---
```{r}
options(width=120, scipen=2, digits=2)
suppressPackageStartupMessages(library(OpenMx))
suppressPackageStartupMessages(library(rpf))
suppressPackageStartupMessages(library(ifaTools))
library(xtable)
options(xtable.type='html')
# Adjust the path in the next statement to load your data
data <- read.csv(na.strings=c(), file='fake.csv',stringsAsFactors=FALSE,check.names=FALSE)
colnames(data) <- mxMakeNames(colnames(data), unique=TRUE)
data[['case']] <- NULL # excluded
data[['twin']] <- NULL # excluded
data[['vetsaid']] <- NULL # excluded
data[['age']] <- NULL # excluded
data[['afqtrt']] <- NULL # excluded
data[['afqtms']] <- NULL # excluded
data[['afqtpct']] <- NULL # excluded
data[['afqtpcttran']] <- NULL # excluded
data[['zyg10']] <- NULL # excluded
data[['afqtwg']] <- NULL # excluded
factors <- "afqt"
numFactors <- length(factors)
spec <- list()
spec[1:100] <- list(rpf.grm(factors=numFactors,outcomes=2))
names(spec) <- c("afqt1n", "afqt2n", "afqt3n", "afqt4n", "afqt5n", "afqt6n",
"afqt7n", "afqt8n", "afqt9n", "afqt10n", "afqt11n", "afqt12n",
"afqt13n", "afqt14n", "afqt15n", "afqt16n", "afqt17n", "afqt18n",
"afqt19n", "afqt20n", "afqt21n", "afqt22n", "afqt23n", "afqt24n",
"afqt25n", "afqt26n", "afqt27n", "afqt28n", "afqt29n", "afqt30n",
"afqt31n", "afqt32n", "afqt33n", "afqt34n", "afqt35n", "afqt36n",
"afqt37n", "afqt38n", "afqt39n", "afqt40n", "afqt41n", "afqt42n",
"afqt43n", "afqt44n", "afqt45n", "afqt46n", "afqt47n", "afqt48n",
"afqt49n", "afqt50n", "afqt51n", "afqt52n", "afqt53n", "afqt54n",
"afqt55n", "afqt56n", "afqt57n", "afqt58n", "afqt59n", "afqt60n",
"afqt61n", "afqt62n", "afqt63n", "afqt64n", "afqt65n", "afqt66n",
"afqt67n", "afqt68n", "afqt69n", "afqt70n", "afqt71n", "afqt72n",
"afqt73n", "afqt74n", "afqt75n", "afqt76n", "afqt77n", "afqt78n",
"afqt79n", "afqt80n", "afqt81n", "afqt82n", "afqt83n", "afqt84n",
"afqt85n", "afqt86n", "afqt87n", "afqt88n", "afqt89n", "afqt90n",
"afqt91n", "afqt92n", "afqt93n", "afqt94n", "afqt95n", "afqt96n",
"afqt97n", "afqt98n", "afqt99n", "afqt100n")
missingColumns <- which(is.na(match(names(spec), colnames(data))))
if (length(missingColumns)) {
stop(paste('Columns missing in the data:', omxQuotes(names(spec)[missingColumns])))
}
for (col in c("afqt1n", "afqt2n", "afqt6n", "afqt8n", "afqt25n", "afqt49n"
)) {
data[[col]] <- mxFactor(data[[col]], levels=0:1)
}
for (col in c("afqt3n", "afqt4n", "afqt5n", "afqt7n", "afqt9n", "afqt10n",
"afqt11n", "afqt12n", "afqt13n", "afqt14n", "afqt15n", "afqt16n",
"afqt17n", "afqt18n", "afqt19n", "afqt20n", "afqt21n", "afqt22n",
"afqt23n", "afqt24n", "afqt26n", "afqt27n", "afqt28n", "afqt29n",
"afqt30n", "afqt31n", "afqt32n", "afqt33n", "afqt34n", "afqt35n",
"afqt36n", "afqt37n", "afqt38n", "afqt39n", "afqt40n", "afqt41n",
"afqt42n", "afqt43n", "afqt44n", "afqt45n", "afqt46n", "afqt47n",
"afqt48n", "afqt50n", "afqt51n", "afqt52n", "afqt53n", "afqt54n",
"afqt55n", "afqt56n", "afqt57n", "afqt58n", "afqt59n", "afqt60n",
"afqt61n", "afqt62n", "afqt63n", "afqt64n", "afqt65n", "afqt66n",
"afqt67n", "afqt68n", "afqt69n", "afqt70n", "afqt71n", "afqt72n",
"afqt73n", "afqt74n", "afqt75n", "afqt76n", "afqt77n", "afqt78n",
"afqt79n", "afqt80n", "afqt81n", "afqt82n", "afqt83n", "afqt84n",
"afqt85n", "afqt86n", "afqt87n", "afqt88n", "afqt89n", "afqt90n",
"afqt91n", "afqt92n", "afqt93n", "afqt94n", "afqt95n", "afqt96n",
"afqt97n", "afqt98n", "afqt99n", "afqt100n")) {
data[[col]] <- mxFactor(data[[col]], levels=c("0", "1"), exclude="NA")
}
#set.seed(1) # uncomment to get the same starting values every time
startingValues <- mxSimplify2Array(lapply(spec, rpf.rparam))
rownames(startingValues) <- paste0('p', 1:nrow(startingValues))
rownames(startingValues)[1:numFactors] <- factors
imat <- mxMatrix(name='item', values=startingValues, free=!is.na(startingValues))
# Remove non-factor columns (if any)
data <- data[,sapply(data, is.factor)]
# Remove all-missing rows
data <- data[rowSums(!is.na(data[,names(spec)])) != 0,]
itemModel <- mxModel(model='itemModel', imat,
mxData(observed=data, type='raw'),
mxExpectationBA81(ItemSpec=spec),
mxFitFunctionML())
emStep <- mxComputeEM('itemModel.expectation', 'scores',
mxComputeNewtonRaphson(), verbose=2L,
information='oakes1999', infoArgs=list(fitfunction='fitfunction'))
computePlan <- mxComputeSequence(list(EM=emStep,
HQ=mxComputeHessianQuality(),
SE=mxComputeStandardError()))
itemModel <- mxModel(itemModel, computePlan)
m1Fit <- mxRun(itemModel)
m1Grp <- as.IFAgroup(m1Fit, minItemsPerScore=1L)
if(0) {
fake <- rpf.sample(1300, grp=m1Grp)
write.csv(file="fake.csv", fake, row.names = FALSE)
}
```
An item factor model was fit with `r length(factors)`
factors (`r ifelse(length(factors), factors, '-')`), -2LL=$`r m1Fit$output$fit`$.
The condition number of the information matrix was `r round(m1Fit$output$conditionNumber)`.
```{r,fig.height=2}
got <- sumScoreEAPTest(m1Grp)
df <- data.frame(score=as.numeric(names(got[['observed']])),
expected=got[['expected']], observed=got[['observed']])
df <- melt(df, id='score', variable.name='source', value.name='n')
ggplot(df, aes(x=score, y=n, color=source)) + geom_line()
```
```{r,fig.height=3}
basis <- rep(0, length(factors))
basis[1] <- 1
plotInformation(m1Grp, width=5, basis=basis)
```
```{r}
summary(m1Fit, refModels=mxRefModels(itemModel, run = TRUE))
```
```{r,results='asis'}
citation('OpenMx')
```