We introduce a shiny web application to facilitate the construction of Item Factor Analysis (a.k.a. Item Response Theory) models using the OpenMx package. The web application assists with importing data, outcome recoding, and model specification. However, the app does not conduct any analysis but, rather, generates an analysis script. Generated Rmarkdown output serves dual purposes: to analyze a data set and demonstrate good programming practices. The app can be used as a teaching tool or as a starting point for custom analysis scripts.
OpenMx, a modular package originally designed for structural equation modeling (Neale et al. in press), recently gained the ability to handle Item Factor Analysis (a.k.a. Item Response Theory, Modern Test Theory) models (Pritikin et al. 2015). Although a goal of OpenMx is to cater to the statistical power user and facilitate analyses that are difficult to conduct in other software, the development team is always on the lookout for ways to ease the learning curve for novice users as well. Here we introduce a new shiny (RStudio and Inc. 2014) web application to generate OpenMx code in Rmarkdown format (Allaire et al. 2014). We believe this code generator substantially lowers the barrier to entry for novice users of Item Factor Analysis (IFA) and encourages a culture of literate programming (Knuth 1984) and reproducible science (Peng 2011; Nosek et al. 2015). The generated code can be customized at many levels. This flexibility enables the production of custom analyses and reports as users grow more sophisticated in their modeling expectations.
Item analysis is concerned with items that are scored correct/incorrect
or on an ordinal scale. Many psychological surveys use an ordinal scale.
For example, participants may be asked to respond to an item like, “I
feel I am in charge of the situation in which I live.” on a 5-point
Likert scale from agree to disagree. Whether dichotomous or ordinal,
the conditional likelihood of response
Several models are readily available to plug in as the response
probability function
The dichotomous response probability function can model items when there
are exactly two possible outcomes. It is defined as,
where
The graded response model is a response probability function for 2 or
more outcomes (Samejima 1969; Cai 2010b). For outcomes k in 0 to K,
slope vector
The nominal model is a response probability function for items with 3 or
more outcomes (e.g., Thissen et al. 2010). It can be defined as,
Modern test theory employs item models,
To further inform your intuition about item models, it can be helpful to
place yourself in the position of the optimization algorithm. Enter the
following commands to launch the model explorer tool and browse to the
output web server address. It is possible to do this without RStudio,
but RStudio makes everything easier so we recommend using RStudio.
Note that the port number (3726
printed below) may be different on
your computer.
> library(ifaTools)
> itemModelExplorer()
Listening on http://127.0.0.1:3726
Your browser should show a screen similar to Figure 2. Try experimenting with all the controls. Early in the development of item models, model parameters closely corresponded to the psychological concepts of difficulty and discrimination (Birnbaum 1968). For example, difficult items are only answered correctly by the brightest examinees while most examinees may correctly answer easy items. Discrimination quantifies how much we learn from a given response. Well-designed items discriminate examinee skill. The causes of poor item discrimination are many. An item may be hurt by awkward wording, by asking examinees something that is somewhat off-topic, or by asking the same question in slightly different ways.
Some item model parameters still retain a close connection to difficulty
and discrimination. For example, the dichotomous model’s
The way to interpret these plots is to start by picking a value for the
latent trait. Suppose we know that examinee Alice has a latent ability
of 2 logit units. If we trace across the plot where the
Much can be gleaned about item models by inspection of these plots.
However, it is worth conveying a few additional details specific to
particular item models. The dichotomous model’s plogis
function. The
By default, the nominal model uses trend for the T.a
and T.c
matrices (Equation (5)). This parameterization is also
known as the Fourier basis. The effect is that the alf
and gam
parameters control the lowest frequency variation to the highest
frequency variation. To develop an intuition for how this works, set all
parameters to zero then set a
, alf1
and gam2
to 1. Experiment with
the gam
parameters before you experiment with the alf
parameters.
Refer to Thissen et al. (2010) for discussion of the possibilities of this item
model. Custom T.a
and T.c
matrices are not available in the model
explorer app, but can be specified in R code.
The “Show/Edit parameters
” checkbox has a special didactic purpose.
Open two copies of the item model explorer. On one copy, un-check the
“Show/Edit parameters
” checkbox to hide the parameters and click the
“Draw new parameters
” button. On the second copy of the item model
explorer, adjust the item model parameters to try to match the plot
produced on the first item model explorer. You can check your answers by
checking the “Show/Edit parameters
” checkbox. When you play this game,
you are doing part of what the optimization algorithm does when it fits
models to data. Note that there is no need to practice this skill. The
computer will do it for you.
Enter the following commands to launch the model builder tool and browse
to the output web server address. As before, it is possible to do this
without RStudio, but RStudio makes everything easier so we recommend
using RStudio. Note that the port number (3726
printed below) may be
different on your computer.
> library(ifaTools)
> modelBuilder()
Listening on http://127.0.0.1:3726
g341-19.csv
data.Your browser should show a screen similar to
Figure 3. Take care not to hit the Reload
button
because that will reset the app. Learn how to save your work (detailed
below) before you experiment with the Reload
button. Across the top
are tabs that organize the major functions of the model builder app. On
the left side is a control panel for the currently selected tab Data
.
Example data sets are available at the bottom of the control panel. You
are welcome to experiment with these, but we will focus on the process
of loading an external data set. If you prefer to follow along with a
video then browse to http://youtu.be/xHeb5_CWnCk for dichotomous data
and http://youtu.be/iwtpleltteQ for polytomous data.
Click on the “Choose File
” buttonChoose CSV File
’.” The exact appearance may differ
on your system.g341-19.csv
, a
dichotomous data set that is available in the
ifaTools package
(Pritikin 2015a). The first 6 lines will appear in the “Unparsed content
”
section (see Figure 4).Parsed content
” section reports an
error. By reading the error carefully, you will find mention that
“duplicate ‘row.names’ are not allowed.” Since “Row names?
” is enabled
in the control panel, the model builder app expects the first column of
data to be dedicated to row names. A row name is typically the
examinee’s name or numeric identifier. A glance at the unparsed content
reveals that no row names have been given in this data set.
g341-19.csv
data set when parsed
incorrectly as a single column.Click the “Row names?
” checkbox in the control panel to disable row
names. Immediately (without reloading the data), the error message in
the “Parsed content
” section will be replaced by some of the data
organized into a single column. The column name will read
X010111111100
. A column name like that should raise suspicion. Since
the “Header?
” checkbox is enabled in the control panel, the model
builder app expects the first line of the data to contain column names.
Therefore, the first line of data is misinterpreted.
Click the “Header?
” checkbox in the control panel to disable column
names. The column in the “Parsed content
” section will now be labeled
V1
. Click on the “Item summary
” control as an alternate way to
verify that the data is loaded and parsed accurately. The main content
area includes two elements, a selection for the “Row frequency column
”
and a table of items by Outcomes
and Missing
(see
Figure 5). The “Row frequency column
” selection
is used when you have already reduced your data to unique rows and row
counts. The example data set LSAT6
is in this format. For our current
data set, leave “Row frequency column
” set to
The information conveyed in the item table should rouse suspicion. There
is only 1 row (or 1 item) with 721 outcomes. What happened is that the
parsing is still not working and all the items are treated as a single
column. For example, the first row “0 1 0 1 1 1 1 1 1 1 0 0” is not
treated as 12 separate numbers but as a single uninterpreted unit. To
fix the parsing, we need to select Whitespace
as the Separator
in
the control panel. With this last error corrected, the table is updated
with 12 items labeled V1, V2, …, V12 and all with 2 outcomes. With the
data correctly loaded, click on the “Outcomes
”tab on the top bar.
The control panel on the “Outcomes
” tab packs a lot of functionality
(Figure 6). The first two selectors are mutually
exclusive and permit working with all items having the same outcomes or
with specific items, respectively. The outcome set “V1” is currently
selected as seen in the control panel on the left side. In these data,
all items have the same possible outcomes (0 or 1). Therefore, there is
only one outcome set. The name “V1” does not just refer to the item “V1”
but to all items, because all items have the same outcomes.
For clarity, it is often helpful to rename outcomes. The “Recode from
”
selector should have “0” selected. Change the to
selector to
<Rename>
, enter “incorrect” for the “New name
” value, and click the
“Add mapping
” button. This will create a recoding rule that will show
up in the “Recode Table
” output (Figure 7).
Similarly, rename the “1” outcome to “correct” and again click the
“Add mapping
” button. At this point, you should have 2 rules in the
“Recode Table
” output.
Add mapping
” button. (right) The outcome reorder tool
with 1 reordering rule.
Click on the “Reorder
” tab. Here you will find the outcomes sorted in
lexical order. Drag one of the outcomes to reverse the order (as in
Figure 7). Similar to the operation of renaming
outcomes, this reordering step is not strictly necessary but is often
helpful to keep things straight in our minds. With dichotomous outcomes,
there are not that many ways that things can go wrong. However, it is
good practice to use self-explanatory labels and standardized ordering.
This is especially true when there are more than 2 outcomes to worry
about.
We will not use the “Reverse
” tab and other control panel elements in
the present example. In survey design, it is common for outcomes to have
a canonical order with some items reverse scored. The “Reverse
” tab is
used to reverse the outcome order of some subset of items without
dragging the outcomes around with the “Reorder
” tool. This can save a
lot of dragging when there are more than 2 outcomes. The “Add outcome
”
tool permits the addition of outcomes that are not represented in the
data. This might happen when a measure is in development and we are
fitting a preliminary model just to get a vague idea of how the item is
working. To fit an item that lacks data on some outcomes, it is usually
necessary to use the nominal response model with a less than full rank
parameterization (similar to Thissen et al. 1989). In addition to renaming, the
recode mappings tool can merge or collapse outcomes.
For example, we might have an outcome set consisting of “Agree,” “Agree
somewhat,” “Not sure,” “Disagree somewhat,” and “Disagree.” It is
straightforward to merge “Agree somewhat” to “Agree” and “Disagree
somewhat” to “Disagree,” resulting in only 3 outcomes. Of course, it is
not always obvious which outcomes to merge. The recode tool can also
recode an outcome to <NA>
, which causes that outcome to be interpreted
as missing. Finally, the “Discard
” button at the bottom of the control
panel allows us to discard any rule that we created. So feel free to
experiment.
Click on the “Model
” tab on the top bar. The first decision we need to
make is how many latent factors to include in our model
(Figure 8). If we are not sure, a good starting
point is 1. By default, the first latent factor is called teacup
. In
case there was any doubt, “teacup” is not a very good name for a latent
trait. We deliberately picked ridiculous factors names to encourage
users to pick names that make sense in the context of the data under
analysis. For example, a good factor name for a math test might be
math
. If you cannot think of a good factor name, you could use “latent
trait,” but this name only works well for a single factor model. You
really should make an effort to think of descriptive names before you
start using trait1, trait2, etc. If you are not sure how many factors to
use then use 1 for now. We will revisit this question later.
The “Reorder
” tab allows you to change the order of your items. This
can be helpful because of the way that item model and parameter editing
works. To get familiar with how item editing works, click on the
“Parameters
” tab. The main content area of the “Parameters
” tab
contains a lot of information. The first thing to notice is that all of
the tables have the same column labels. Each item is assigned to a
column. Using the control panel, we will focus on only a subset of
items. Change the first selector “Edit items from:
” from V1
to V7
.
This will hide the first 6 items, making the tables in the main content
area look more manageable (Figure 9). The first
two selectors facilitate item range selection. To reveal all items
again, use the “Focus all items
” button. Item selection is important
to understand because the remainder of the controls in the control panel
operate on only the selected (visible) items.
With only items V7
to V12
visible, just to demonstrate how it is
done, let us place an equality constraint on the slope or latent factor
loading. Type “slope” into the Label
textbox and click the
“Set Label
” button. The label slope
should appear in all columns of
the first row of the Labels
table in the main content area. Now let us
switch to the first 6 items. This can be accomplished in a variety of
ways. One way is to change the first selector from V7
to V6
and the
second selector from V12
to V1
.
With only items V1
to V6
visible, select drm
from the “Model
”
selector. The value drm
is an abbreviation for the 4PL dichotomous
response model (Loken and Rulison 2010), which has four parameters when there is one
factor. The g
and u
rows should appear in all of the tables in the
main content area. Parameter g
is the pseudo-guessing lower bound and
u
is the upper bound. The upper bound has been used to better model
high ability examinees in a Computerized Adaptive Testing context
(Magis 2013). Even high ability examinees may occasionally miss an easy
item. Here we will fix the upper bound to 1 (meaning that an examinee
with sufficiently high ability will never answer incorrectly). Since the
bound parameters are expressed on a logit scale, we will use
u
from the “Parameter
” selector and Inf
from the “Free
” selector (since Starting values
” for u
should change to Inf
and the corresponding
“Is free
” row should change from TRUE
to FALSE
. With this
constraint, the 4PL dichotomous response model is equivalent to the 3PL
model (Birnbaum 1968).
The pseudo-guessing lower bound g
represents the chance that an
examinee will get the item correct by guessing. Typically, the expected
guessed-correct probability for a 3-alternative item is
As a compromise between fixing and freeing, a Bayesian prior can be used
with the mode of the prior set to the expected probability. While some
researchers are uneasy about the use of priors (Gelman 2008), the
practice is not new (e.g, Baker and Kim 2004 7). The prior ensures that
a parameter will converge even when there is not enough data to estimate
it, but at the same time, the model retains some flexibility to adapt to
unexpected data. To set a prior, drag the “Prior mode
” slider and
click the “Set
” button. Let us imagine that these items had 4
alternatives. Select g
from the Parameter
selector, move the
“Prior mode
” slider to 4, and click the nearby Set
button. Two
tables will change. Each cell of the g
row of the Labels
table will
be assigned a unique label. This is necessary because Bayesian priors
are associated with labels, not with particular parameters. In addition,
the “Bayesian prior mode
” table will show g
row. The logit usage reflects that the parameter is estimated on the
logit scale, but it is much easier for humans to understand a
probability expressed as a fraction rather than raw logit units.
We will not use the “Nominal Tc
” selector for these data.
“Nominal Tc
” only applies to items with more than 2 possible outcomes
when using the nominal response model (Thissen et al. 2010). Before proceeding,
go ahead and click the “Focus all items
” button. Your screen should
look like Figure 10, except for different starting
values. Click on the “Exclude
” tab. This is an easy way to exclude
some of the items from analysis. Finally, click on the “Summary
” tab.
Here you will find a summary of your model settings. Note that the
number of outcomes may differ from the number of outcomes reported in
the summary table found on under the “Data
” top bar page due to
recoding.
We are done setting up our model. Before proceeding, it would be wise to
save our model configuration so we can come back at a later time and
make small adjustments without going through the whole exhaustive
process again. Click Settings
on the top bar. In the main content
area, you will find a preview of what the settings file will look like.
Click the “Download
” button and move the file to a suitable location
on your computer.
To verify that you can reliably restore the saved settings, open a new
browser tab to the same address by pasting the URL address from the
current tab (without closing the current one). You should get a screen
like Figure 3. Again go through the procedure of
loading the data (Figures 4 and
5). Once your data is loaded, click Settings
on
the top bar and load the file that you recently saved. If all goes well,
you should see a screen similar to Figure 11. Go
ahead and look back through the editors under the Outcomes
and Model
sections of the top bar. You should find all your hard work faithfully
preserved. Now you can close either of the 2 browser tabs that you have
open. They both have the same status.
With our model set up and saved, click Analysis
on the top bar. This
screen looks and functions similarly to the Settings
screen. However,
the control panel offers a few options specific to code generation. The
“Functional form for dichotomous bound prior density
” selector chooses
the distributional form of the Bayesian prior. Logit-normal
is the
more recent scheme (Cai et al. 2011). The “Information matrix method
” control
is set to Oakes
by default. In a simulation study included with
OpenMx, the Oakes method (Oakes 1999) exhibited accuracy comparable to
central difference with Richardson extrapolation and time linear in the
number of parameters. Click the “Download
” button and save the
Rmarkdown
code. The Rmarkdown
file and your data need to be in the
same directory. Either move the Rmarkdown
file to your data directory,
or alternately, you can specify a full path in the read.csv statement
(lines 16–17). Open the file in RStudio and click the “Knit HTML
”
button.
---
title: "g341-19"
date: "14-Nov-2014"
output: html_document
---
```{r}
options(width = 120, scipen = 2, digits = 2) ## \label{e1:options}
suppressPackageStartupMessages(library(OpenMx)) ## \label{e1:pkg1}
suppressPackageStartupMessages(library(rpf))
suppressPackageStartupMessages(library(ifaTools))
library(xtable)
options(xtable.type = 'html') ## \label{e1:pkg2}
# Adjust the path in the next statement to load your data
data <- read.csv(file = 'g341-19.csv', header = FALSE, sep = ' ', ## \label{e1:load1}
stringsAsFactors = FALSE, check.names = FALSE) ## \label{e1:load2}
colnames(data) <- mxMakeNames(colnames(data), unique = TRUE) ## \label{e1:colnames}
factors <- "ability" ## \label{e1:factor1}
numFactors <- length(factors)
spec <- list()
spec[1:6] <- rpf.drm(factors = numFactors)
spec[7:12] <- rpf.grm(factors = numFactors, outcomes = 2) ## \label{e1:factor2}
names(spec) <- c("V1", "V2", "V3", "V4", "V5", "V6", "V7", "V8", "V9", "V10",
"V11", "V12")
missingColumns <- which(is.na(match(names(spec), colnames(data)))) ## \label{e1:colmatch}
if (length(missingColumns)) {
stop(paste('Columns missing in the data:',
omxQuotes(names(spec)[missingColumns])))
}
data[names(spec)] <- mxFactor(data[names(spec)], levels = 0:1, ## \label{e1:mxfactor}
labels = c("incorrect", "correct"))
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)) ## \label{e1:item}
imat$free['p4',1:6] <- FALSE
imat$values['p4',1:6] <- Inf
imat$labels['ability',7:12] <- 'slope'
imat$labels['p3',1:1] <- 'V1_g'
imat$labels['p3',2:2] <- 'V2_g'
imat$labels['p3',3:3] <- 'V3_g'
imat$labels['p3',4:4] <- 'V4_g'
imat$labels['p3',5:5] <- 'V5_g'
imat$labels['p3',6:6] <- 'V6_g'
hasLabel <- !is.na(imat$labels)
for (lab1 in unique(imat$labels[hasLabel])) {
imat$values[hasLabel & imat$labels == lab1] <-
sample(imat$values[hasLabel & imat$labels == lab1], 1)
}
data <- compressDataFrame(data)
itemModel <- mxModel(model = 'itemModel', imat,
mxData(observed = data, type = 'raw',
numObs = sum(data[['freq']]), sort = FALSE),
mxExpectationBA81(ItemSpec = spec, weightColumn = 'freq'),
mxFitFunctionML())
priorLabels <- c("V1_g", "V2_g", "V3_g", "V4_g", "V5_g", "V6_g")
priorParam <- mxMatrix(name = 'priorParam', nrow = 1,
ncol = length(priorLabels), free = TRUE, labels = priorLabels)
priorParam$values <- imat$values[ match(priorParam$labels, imat$labels) ]
priorMode <- c(priorParam$values)
priorMode[1:6] <- logit(1/4)
priorModel <- univariatePrior('logit-norm', priorLabels, priorMode)
container <- mxModel(model = 'container', itemModel, priorModel, ## \label{e1:container}
mxFitFunctionMultigroup(groups = c('itemModel.fitfunction',
'univariatePrior.fitfunction')))
emStep <- mxComputeEM('itemModel.expectation', 'scores',
mxComputeNewtonRaphson(), verbose = 2L, ## \label{e1:verbose}
information = 'oakes1999', infoArgs = list(fitfunction = 'fitfunction'))
computePlan <- mxComputeSequence(list(EM = emStep,
HQ = mxComputeHessianQuality(),
SE = mxComputeStandardError()))
m1Fit <- mxRun(container) ## \label{e1:run}
m1Grp <- as.IFAgroup(m1Fit$itemModel, minItemsPerScore = 1L) ## \label{e1:ifagroup}
```
The details of the generated report are likely to evolve. There may be
some differences between the generated code in this article and the most
recent version, but there should be a correspondence between the basic
elements. The first chunk of code builds the model and optimizes it. We
adjust the output of long lines and numbers (line 8) and load packages
(lines 9–13). The data set is loaded (Figure 5)
in lines 16–17. Latent factors are configured
(Figure 8) in lines 20–24. We strongly encourage
the use of informative column (item) names, but line 18 will take care
of assigning syntactically valid column names if informative names are
not available. The script is crafted such that it can work on other data
sets as long as the required columns are found (line 28). mxFactor
does the recoding and reordering
(Figures 6–7). mxFactor
offers a number of important safety and convenience features beyond what
is available from factor
or ordered
(line 34). The item mxMatrix
(line 43) contains most of the information in the item tables
(Figure 10). Everything goes into the container
model (line 72). The model is optimized (line 83). Since we did not
disable "Show model fitting progress
" (reflected by verbose = 2L
at line 77), we obtain some diagnostics upon knitting the Rmarkdown
to
HTML.
[0] ComputeEM: Welcome, tolerance=1e-09 accel=varadhan2008 info=2
[0] ComputeEM: msteps 2 initial fit 37185.0001 ## \label{e2:traj1}
[0] ComputeEM[2]: msteps 11 fit 34167.9816 rel change 0.0882995805
[0] ComputeEM[3]: msteps 5 fit 33699.978 rel change 0.0138873556
[0] ComputeEM[4]: msteps 14 fit 33549.9723 rel change 0.00447111437
[0] ComputeEM[5]: msteps 5 fit 33455.9478 rel change 0.00281039684
[0] ComputeEM[6]: msteps 3 fit 33454.4705 rel change 4.41596231e-05
[0] ComputeEM[7]: msteps 3 fit 33453.8021 rel change 1.99793343e-05
[0] ComputeEM[8]: msteps 3 fit 33453.2067 rel change 1.77968988e-05
[0] ComputeEM[9]: msteps 2 fit 33453.2062 rel change 1.57420931e-08
[0] ComputeEM[10]: msteps 2 fit 33453.206 rel change 5.03007605e-09
[0] ComputeEM[11]: msteps 2 fit 33453.2059 rel change 2.89615064e-09
[0] ComputeEM[12]: msteps 2 fit 33453.2059 rel change 6.61823582e-10 ## \label{e2:traj2}
[0] ComputeEM: cycles 12/500 total mstep 54 fit 33453.205893 ## \label{e2:fit}
[0] ComputeEM: Oakes1999 method=simple perturbation=0.001000
[0] ComputeEM: deriv of gradient for 0
[0] ComputeEM: deriv of gradient for 1
[0] ...
[0] ComputeEM: deriv of gradient for 24
Given that the starting values are random, the initial fit and
trajectory (lines 88–99) should differ when you try optimizing the same
model but the optimum (line 100) should be the same to within
The function as.IFAgroup
(line 85) is a bridge between OpenMx and
rpf (Pritikin 2015b). The
rpf name is an acronym for response probability function and contains
many IFA-specific diagnostic functions. In addition, rpf can be used
to analyze IFA models optimized by packages other than OpenMx. This
modularity facilitates the comparison of parameter estimates between
packages. While most of the code discussed so far is related to
OpenMx, the remainder of this report will mostly involve rpf.
An item factor model was fit with ` length(factors)`
factors (` factors`), -2LL = $` m1Fit$output$fit`$.
The condition number of the information
matrix was ` m1Fit$output$conditionNumber`.
This is a boilerplate report of model fit. It renders as, “An item
factor model was fit with 1 factors (ability),
Although IFA models consider an examinee’s response pattern as the unit of analysis, the sum-score is often of chief practical importance. For example, students taking the Standardized Aptitude Test for college admission only receive their sum-score and do not even know which items they answered correctly or incorrectly (unless they earned the maximum sum-score). The observation that the sum-score is important has lead to the development of a family of diagnostic tests that examine how well an IFA model predicts sum-scores.
```{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()
```
The first plot provides an overview of how all the items work together
to predict sum-scores (Figure 12). sumScoreEAPTest
also conducts a statistical test to produce a
```{r,results='asis'}
ct <- ChenThissen1997(m1Grp)
print(xtable(ct$pval,
paste('Log p-value of local dependence between item pairs.')))
```
V1 | V2 | V3 | V4 | V5 | V6 | V7 | V8 | V9 | V10 | V11 | |
---|---|---|---|---|---|---|---|---|---|---|---|
V2 | 0.1 | ||||||||||
V3 | |||||||||||
V4 | 1.7 | 1.2 | |||||||||
V5 | 5.1 | ||||||||||
V6 | 0.2 | 0.1 | 0.3 | ||||||||
V7 | 0.1 | 6.4 | 0.6 | 2.4 | 0.3 | 5.0 | |||||
V8 | 0.1 | 3.9 | |||||||||
V9 | 3.6 | 10.7 | 1.1 | 2.6 | 1.8 | 5.8 | 37.1 | 0.3 | |||
V10 | 9.1 | 0.3 | 0.1 | 10.2 | 16.2 | ||||||
V11 | |||||||||||
V12 | 3.9 | 0.8 | 0.6 |
A test of local dependence is important to examine, as it is an
assumption of IFA (Yen 1993). Table 1 exhibits the log
```{r,results='asis'}
sfit <- SitemFit(m1Grp)
tbl <- t(sapply(sfit, function(r)
c(n = r$n, df = r$df, stat = r$statistic, pval = r$pval)))
print(xtable(tbl, paste0('Sum-score item-wise fit.'))
```
n | df | statistic | ||
---|---|---|---|---|
V1 | 2844 | 8 | 6.58 | |
V2 | 2844 | 8 | 7.30 | |
V3 | 2844 | 8 | 7.17 | |
V4 | 2844 | 8 | 10.12 | |
V5 | 2844 | 8 | 19.00 | |
V6 | 2844 | 8 | 8.50 | |
V7 | 2844 | 9 | 33.45 | |
V8 | 2844 | 9 | 5.48 | |
V9 | 2844 | 9 | 34.42 | |
V10 | 2844 | 10 | 12.61 | |
V11 | 2844 | 8 | 43.06 | |
V12 | 2844 | 8 | 20.20 |
Sum-score item fit is another tool for model assessment
(Orlando and Thissen 2000; Kang and Chen 2008). Each item is tested against the null
hypothesis that the other items accurately predict the outcome of the
item of interest (Table 2). Again SitemFit
help for details on the
options for missing data.
```{r,fig.height=2}
map1 <- itemResponseMap(m1Grp, factor = 1)
ggplot(map1, aes(x = score, y = item, label = outcome)) +
geom_text(size = 4, position = position_jitter(h = .25))
```
An item response map is a crude tool for diagnosing certain model misspecifications (Figure 13). Each outcome is assigned the average latent score of all the examinees who picked that outcome. Usually we advocate the use of the actual outcome labels (e.g., “incorrect” and “correct”) instead of numbers. For this plot, however, we make an exception because the numbers make it easy to inspect whether the outcomes are in ascending order. If descending order is observed then it is worth checking whether the item needs to be reverse scored or to consider whether the item was misinterpreted by some examinees. If the response data were manually collected then the data entry process should also be checked for errors.
```{r,fig.height=3}
pl <- lapply(names(sfit), function(item) { SitemPlot(sfit, item) })
for (px in 1:length(pl)) {
print(pl[[px]])
}
Two approaches are available to plot response probability functions
against a latent trait. The same ingredients that go into the production
of Table 2 can also be plotted
(Figure 14). A similar plot can be obtained by plotting
the outcomes probabilities against the latent trait. This is known as an
item characteristic curve plot (Figure 14). The main
advantage of SitemPlot
over iccPlot
is that SitemPlot
is
one-dimensional regardless of the number of latent factors. With
iccPlot
, you must pick a basis vector in the latent space. The tiny
numbers across the probability = 1 line of Figures 14
and 14 are the sample size at that point on the
basis <- rep(0, length(factors))
basis[1] <- 1
plotInformation(m1Grp, width = 5, basis = basis)
```
Figure 15 exhibits item information by latent score.
Similar to iccPlot
, this plot requires the selection of a basis vector
when there is more than 1 latent factor. Notice that items V7-V12 peak
at the same height (near 0.31). This is due to our equality constraint
on the slope or factor loading on these items. By placing this
constraint, we assume a priori that each of these items contributes
exactly the same amount of information.
```{r}
summary(m1Fit)
```
Summary of container
free parameters:
name matrix row col Estimate Std.Error
1 itemModel.item[1,1] itemModel.item ability V1 1.82 0.278
2 itemModel.item[2,1] itemModel.item p2 V1 -0.51 0.230
3 V1_g itemModel.item p3 V1 -1.14 0.208
4 itemModel.item[1,2] itemModel.item ability V2 1.24 0.119
5 itemModel.item[2,2] itemModel.item p2 V2 2.58 0.140
6 V2_g itemModel.item p3 V2 -1.27 0.337
7 itemModel.item[1,3] itemModel.item ability V3 1.56 0.261
8 itemModel.item[2,3] itemModel.item p2 V3 -1.03 0.272
9 V3_g itemModel.item p3 V3 -1.16 0.192
10 itemModel.item[1,4] itemModel.item ability V4 1.36 0.161
11 itemModel.item[2,4] itemModel.item p2 V4 0.41 0.158
12 V4_g itemModel.item p3 V4 -1.10 0.277
13 itemModel.item[1,5] itemModel.item ability V5 1.41 0.196
14 itemModel.item[2,5] itemModel.item p2 V5 -0.47 0.203
15 V5_g itemModel.item p3 V5 -1.03 0.203
16 itemModel.item[1,6] itemModel.item ability V6 1.50 0.130
17 itemModel.item[2,6] itemModel.item p2 V6 1.84 0.119
18 V6_g itemModel.item p3 V6 -1.43 0.317
19 slope itemModel.item ability V7 1.12 0.037
20 itemModel.item[2,7] itemModel.item p2 V7 3.50 0.097
21 itemModel.item[2,8] itemModel.item p2 V8 1.57 0.056
22 itemModel.item[2,9] itemModel.item p2 V9 2.70 0.075
23 itemModel.item[2,10] itemModel.item p2 V10 2.27 0.066
24 itemModel.item[2,11] itemModel.item p2 V11 -0.28 0.047
25 itemModel.item[2,12] itemModel.item p2 V12 0.15 0.047
observed statistics: 720
estimated parameters: 25
degrees of freedom: 695
fit value ( -2lnL units ): 33453
number of observations: 2844
Information Criteria:
| df Penalty | Parameters Penalty | Sample-Size Adjusted
AIC: 32063 33503 NA
BIC: 27926 33652 33573
To get additional fit indices, see help(mxRefModels)
timestamp: 2016-02-24 10:14:59
Wall clock time (HH:MM:SS.hh): 00:00:02.72
OpenMx version number: 2.3.1.254
Need help? See help(mxSummary)
Exhibited above is the OpenMx provided summary of model fit. IFA models are exponential family models so we obtain AIC and BIC. More fit statistics are available if we provide the saturated and independence reference models. Reference models will be requested in our next example.
Since many things are common between dichotomous and polytomous items,
we will move quickly through the process of model set up and result
interpretation. Click on the “Choose File
” button and select
preschool.csv
, a data set from Thissen and Steinberg (1988) available in the
ifaTools package. Click the “Row names?
” checkbox in the control
panel to disable row names. The format of these data are closer to what
is expected by default than our first example so less fiddling is
required. Click on the “Item summary
” tab. Here it appears that there
are 3 items, but the freq
column is not an item. freq
indicates how
many times a row appeared in the original data. These data are
compressed; only unique rows are provided with frequency counts. To
instruct the model builder to interpret the freq
column as frequency
counts, select freq
from the “Row frequency column
” selector
(Figure 16).
This data set is from a preschool test of numerical knowledge. Each item
is actually a combination of 2 dichotomous items. Similar questions were
asked regarding the number 3 and the number 4 and the pattern of
responses mapped to an outcome code. The outcomes should be renamed with
the recoding tool under the “Outcomes
” tab on the top bar (recall
Figure 7). Outcomes 0, 1, 2, and 3 should be
renamed to “neither,” “3 only,” “4 only,” and “both correct,”
respectively, using the “Recode
” tab under the Outcomes
top bar
page. After renaming, reorder the items into the correct order
(Figure 17).
Click Model
on the top bar. On the “Factors
” tab, we will name the
single latent factor “math.” Switch to the Parameters
tab. Here we
select nrm
from the “Model
” selector. The acronym “nrm” stands for
the nominal response model (Thissen et al. 2010). This parameterization of the
nominal model can accommodate basis matrices
With math
and alf1
because they have the same effect on the model and would cause the model
to be unidentified. Fix alf1
to 1. Select alf1
from the
“Parameter
” selector and select 1 from the “Free
” selector. Since we
have worked with this data set already, we know a few things that can
give us a more parsimonious model. The alf2
parameters can be set
equal since both items exhibit poor discrimination between neither
,
3 only
, and 4 only
but good discrimination between these outcomes
and both correct
. Select alf2
with the “Parameter
” selector and
set the label to eq1
. Since both items are equally difficult, we can
equate gam1
. Select gam1
with the “Parameter
” selector and set the
label to eq2
. To avoid overfitting with the highest frequency basis
vector, fix gam3
to 0. Select gam3
with the Parameter
selector and
select 0
with the Free
selector. Figure 17
exhibits the final parameter settings.
Click Analysis
on the top bar. Ensure that “Fit reference models
” is
selected, and download the analysis script. The Rmarkdown
file and
your data need to be in the same directory. Either move the Rmarkdown
file to your data directory, or alternately, you can specify a full path
in the read.csv statement (line 162). Open the file in RStudio and
click the “Knit HTML
” button. Although this is a simple model, it can
take almost 100 E-M cycles to converge. Therefore, we omit reproduction
of the diagnostic output issued during model fit.
---
title: "preschool"
date: "18-Nov-2014"
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(file = 'preschool.csv', stringsAsFactors = FALSE,
check.names = FALSE) ## \label{e3:load}
colnames(data) <- mxMakeNames(colnames(data), unique = TRUE)
data[['freq']] <- as.numeric(data[['freq']]) ## \label{e3:freq}
factors <- "math"
numFactors <- length(factors)
spec <- list()
spec[1:2] <- rpf.nrm(factors = numFactors, outcomes = 4,
T.a = 'trend', T.c = 'trend') ## \label{e3:nominal}
names(spec) <- c("Match", "Identify")
missingColumns <- which(is.na(match(names(spec), colnames(data))))
if (length(missingColumns)) {
stop(paste('Columns missing in the data:',
omxQuotes(names(spec)[missingColumns])))
}
data[names(spec)] <- mxFactor(data[names(spec)], levels = 0:3,
labels = c("neither", "3 only", "4 only", "both correct"))
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))
imat$free['p2',] <- FALSE
imat$values['p2',1:2] <- 1
imat$free['p7',] <- FALSE
imat$values['p7',1:2] <- 0
imat$labels['p3',] <- 'eq1'
imat$labels['p5',] <- 'eq2'
hasLabel <- !is.na(imat$labels)
for (lab1 in unique(imat$labels[hasLabel])) {
imat$values[hasLabel & imat$labels == lab1] <- ## \label{e3:equality}
sample(imat$values[hasLabel & imat$labels == lab1], 1)
}
itemModel <- mxModel(model = 'itemModel', imat,
mxData(observed = data, type = 'raw', numObs = sum(data[['freq']]),
sort = FALSE), ## \label{e3:sort}
mxExpectationBA81(ItemSpec = spec, weightColumn = 'freq'),
mxFitFunctionML())
emStep <- mxComputeEM('itemModel.expectation', 'scores',
mxComputeNewtonRaphson(), verbose = 2L,
information = 'oakes1999',
infoArgs = list(fitfunction = 'fitfunction'))
computePlan <- mxComputeSequence(list(emStep,
mxComputeHessianQuality(),
mxComputeStandardError()))
m1Fit <- mxRun(mxModel(itemModel, computePlan))
m1Grp <- as.IFAgroup(m1Fit, minItemsPerScore = 1L)
```
Although response pattern frequencies are typically natural numbers, fractional frequencies are not prohibited (line 164). A Fourier basis is used for both nominal model transformation matrices (line 170). Customization is limited in the model builder app, but you can use any matrices by editing the generated code. Starting values must respect equality constraints (line 197). By default OpenMx, sorts data prior to optimization. Since our data are already compressed, no benefit would be obtained by sorting (line 202).
An item factor model was fit with ` length(factors)`
factors (` factors`), -2LL = $` m1Fit$output$fit`$.
The condition number of the information matrix was
` round(m1Fit$output$conditionNumber)`.
The boilerplate renders as, “An item factor model was fit with 1 factors
(math),
n | df | statistic | ||
---|---|---|---|---|
Match | 592 | 7 | 9.46 | |
Identify | 592 | 7 | 9.33 |
```{r}
summary(m1Fit, refModels = mxRefModels(m1Fit, run = TRUE)) ## \label{e3:refModels}
```
Summary of itemModel
free parameters:
name matrix row col Estimate Std.Error
1 itemModel.item[1,1] item math Match 0.82 0.37
2 eq1 item p3 Match -1.18 0.27
3 itemModel.item[4,1] item p4 Match 0.50 0.19
4 eq2 item p5 Match 0.18 0.05
5 itemModel.item[6,1] item p6 Match -0.78 0.20
6 itemModel.item[1,2] item math Identify 0.79 0.35
7 itemModel.item[4,2] item p4 Identify -0.25 0.36
8 itemModel.item[6,2] item p6 Identify -1.40 0.21
observed statistics: 15
estimated parameters: 8
degrees of freedom: 7
fit value ( -2lnL units ): 2767
saturated fit value ( -2lnL units ): 2758
number of observations: 592
chi-square: X2 ( df=7 ) = 9.2, p = 0.24
Information Criteria:
| df Penalty | Parameters Penalty | Sample-Size Adjusted
AIC: 2753 2783 NA
BIC: 2723 2819 2793
CFI: 0.98
TLI: 0.97 (also known as NNFI)
RMSEA: 0.023 [95% CI (0, 0.064)]
Prob(RMSEA <= 0.05): 0.88
timestamp: 2016-02-24 10:15:07
Wall clock time (HH:MM:SS.hh): 00:00:04.55
OpenMx version number: 2.3.1.254
Need help? See help(mxSummary)
Although the outcomes are not strictly ordered for Identify
in the
item outcome map (Figure 18), other measures
of model fit look reasonable. The sum-score item fit tests are not
statistically significant at the 0.01 level
(Table 3). This indicates good item-level fit. Since
we requested a saturated and independence model (mxRefModels
; line
223), CFI (Comparative Fit Index), TLI (Tucker Lewis Index), and RMSEA
(Root Mean Square Error of Approximation) are available in the OpenMx
summary and suggest adequate model fit. These relative indices of fit
are appropriate for these data because there are observations for all
possible response patterns. However, be forewarned that as the
multinomial table becomes more sparse, these indices become inaccurate.
For sparse data, a more accurate assessment of model fit is available
using other methods (Bartholomew and Tzamourani 1999; Cai and Hansen 2013).
> container2 <- container
> container2$itemModel$item$labels['ability', ] <- NA
> m3 <- addExploratoryFactors(container2, 0)
> m3 <- mxRun(m3, silent = TRUE)
> mxCompare(m3, m1)
base comparison ep minus2LL df AIC diffLL diffdf p
1 container1 <NA> 30 33369 690 31989 NA NA NA
2 container1 container1 25 33454 695 32064 85 5 7.7e-17
> m4 <- addExploratoryFactors(container2, 1)
> m4 <- mxRun(m4, silent = TRUE)
> mxCompare(m4, m2)
base comparison ep minus2LL df AIC diffLL diffdf p
1 container2 <NA> 41 33325 679 31967 NA NA NA
2 container2 container2 36 33339 684 31971 14 5 0.013
> grid.arrange(plotTwoFactors(m2$itemModel$item$values[1:2, ]) +
+ labs(title = "a."), plotTwoFactors(m4$itemModel$item$values[1:2, ]) +
+ labs(title = "b."), ncol = 2)
plotTwoFactors
is
given in the Appendix.
A Rasch model is obtained when all slope parameters are constrained to
be equal and the variance is fixed to 1.0, or equivalently, all slopes
are fixed to 1.0 with free variance (Rasch 1960/1993). If your interest is
Rasch models with a single latent factor then you can take advantage of
Rasch residual-based fit statistics. Infit and outfit are available from
pf.1dim.fit
.
A common problem is that we do not know how many latent factors to
employ to most accurately model our data. Fortunately, there is a method
item factor analysis (Bock et al. 1988) analogous to factor analysis of
continuous indicators (Lovie and Lovie 1996). We will employ the likelihood ratio
test for inference. The likelihood ratio test is asymptotically
consistent for sparse multinomial distributions (Haberman 1977).
However, in finite samples, we should not expect that the null
distribution is well calibrated. In brief, the
> m1 <- addExploratoryFactors(container, 0)
> m1 <- mxRun(m1, silent = TRUE)
> m2 <- addExploratoryFactors(container, 1)
> m2 <- mxRun(m2, silent = TRUE)
> mxCompare(m2, m1)
base comparison ep minus2LL df AIC diffLL diffdf p
1 container2 <NA> 36 33339 684 31971 NA NA NA
2 container2 container1 25 33454 695 32064 115 11 1.5e-19
Here we find that there is reasonably good support in favor of a two
factor solution. However, the slope of items 7-12 are constrained equal.
Maybe this constraint was a mistake. It is possible that these items are
well modeled by a single factor when all the slopes are freed. We cannot
directly compare m2
against a single factor model without the slope
constraint because these models are not nested. However, we can make a
number of similar comparisons.
We find that there is a dramatic improvement in fit whether we relax the constraint on items 7-12 or we add another factor. Without knowing more about how the data were collected, parsimony favors a single factor model without constraints on the slopes. We can further check this idea by comparison of two factor models with and without the slope constraint (Figure 19).
A
Recall that the optimization algorithm uses equal interval quadrature to
evaluate the integral in Equation (2). It is important to
understand how the quadrature grid influences model optimization
accuracy and time. Let
Figure 20 exhibits a simulation study of the influence of
quadrature on model accuracy. All comparisons are against a 41 point
quadrature of width 5.0. Before computing the Euclidean distance
(
Item factor analysis with more than two factors requires patience and expertise. Model optimization time becomes an uncomfortable hindrance to experimentation. An optimization algorithm better suited to many latent factors, such as the Metropolis-Hastings Robbins-Monro algorithm (Cai 2010b), is not yet available in OpenMx. The model builder offers as many as five factors because additional factors do not always increase estimation time. Suppose all items load on a general factor. In the special case that each item loads on at most one additional factor, many additional factors will not increase estimation time. One important use for this kind of factor structure is to account for local dependence (DeMars 2006). For example, a reading comprehension test might have 3-4 items that relate to a single passage. The items within each passage will likely exhibit local dependence. One way to account for this kind of test structure is to add passage specific latent factors. Since the passages are disjoint, all of the passage specific factors will count as a single factor with respect to estimation time (Cai 2010a).
We gave detailed instructions on how to set up IFA models for analysis of both dichotomous and polytomous data using the model builder app. We hope this will ease the learning curve for the construction of IFA models in OpenMx. The model builder app offers limited flexibility by design to reduce the number of options for novice users. For example, there is no facility for construction of multiple group models. This may be construed as a disadvantage, but we argue that keeping the app as simple as possible is important for newcomers to IFA. Learning OpenMx can be a daunting prospect. OpenMx, rpf, and ifaTools are free software. The source code is available for everybody to view, modify, and use. If you find this software useful, we hope you will cite us in your publications.
Factors are plotted in a coordinate system determined by a varimax rotation (line 2). Promax axes are superimposed (line 9).
plotTwoFactors <- function(slope) {
lvm <- varimax(toFactorLoading(slope))$loadings ## \label{e4:varimax}
if (any(abs(lvm[lvm < 0]) > .001)) stop("Got negative loadings")
lvm[lvm<0] <- 0
df <- as.data.frame(lvm[, 1:2])
df$name <- rownames(df)
pl <- ggplot(df, aes_string(x = rownames(slope)[1],
y = rownames(slope)[2], label = "name")) + geom_text(size = 3)
pm <- promax(lvm[, 1:2])$rotmat ## \label{e4:promax}
for (dx in 1:ncol(pm)) {
d1 <- .5 * pm[, dx] / sqrt(sum(pm[, dx]^2))
pl <- pl + geom_segment(x = .5, y = .5, xend = d1[1] + .5,
yend = d1[2] + .5, arrow = arrow(length = unit(.5, "cm")))
}
pl + xlim(0, 1) + ylim(0, 1)
}
OpenMx, shiny, Rmarkdown, ifaTools, rpf
MissingData, Psychometrics, WebTechnologies
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
N. Pritikin & M. Schmidt, "Model Builder for Item Factor Analysis with OpenMx", The R Journal, 2016
BibTeX citation
@article{RJ-2016-013, author = {N. Pritikin, Joshua and M. Schmidt, Karen}, title = {Model Builder for Item Factor Analysis with OpenMx}, journal = {The R Journal}, year = {2016}, note = {https://rjournal.github.io/}, volume = {8}, issue = {1}, issn = {2073-4859}, pages = {182-203} }