## One factor IRT Model (2pl).  100 items - data faked to be similar to AFQT data from VETSA study

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)
}

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()

basis <- rep(0, length(factors))
basis[1] <- 1
plotInformation(m1Grp, width=5, basis=basis)
summary(m1Fit, refModels=mxRefModels(itemModel, run = TRUE))