This document describes the process of adding custom item models to rpf as was done by a contributor to the rpf package (Carl F. Falk) with assistance from the package author (Joshua N. Pritikin) for the logistic function of a monotonic polynomial (LMP) item model (Falk & Cai, in press).
This document assumes experience in editing R and C++ code and installing R packages from source. In addition, use of versioning control software (e.g., git) is also helpful. Not mentioned are platform-specific tools or compilers that need to be in place in order to compile/install packages from source.
The rpf package is modular in the sense that one only needs to write the underlying functions specific to the item model, such as the traceline (ICC, IRF, etc.), derivatives w.r.t. item parameters and the latent trait(s), and a few utility functions. As a byproduct, estimation of a measurement model that utilizes the item model may be possible via use of an external package (e.g., OpenMx)
Obtain rpf source code from CRAN or another source such as GitHub. I used the latter as I wanted to use git to track my own changes and then later merge with the rpf package.
Most work needs to be done in C++ using the first file mentioned in the subsequent section. Testing appears to be easiest from R, however, with some functions in the rpf package yielding access to the underlying C++ functions. It is therefore possible to also glean some insight into what the underyling C++ functions ought to do, and their input by examining rpf documentation and experimenting with associated functions.
The bulk of work that needs to be done to implement an item model takes place in the src/libifa-rpf.cpp file.
For example, librpf_model[] towards the end of the
file contains available item models with associated C++
functions required to implement each item model. As a
concrete example, the following code chunk
effectively lists all functions used by the
unidimensional LMP item model:
{"lmp",
irt_rpf_1dim_lmp_numSpec,
irt_rpf_1dim_lmp_numParam,
irt_rpf_1dim_lmp_paramInfo,
irt_rpf_1dim_lmp_prob,
irt_rpf_logprob_adapter,
irt_rpf_1dim_lmp_deriv1,
irt_rpf_1dim_lmp_deriv2,
irt_rpf_1dim_lmp_dTheta,
irt_rpf_1dim_lmp_rescale,
}
To add a novel item model, one would need to add an
analogous chunk of code to librpf_model[] with
the given functions pertaining to those used by the
novel item model. Each type of function takes input
in the same way regardless of the item model. This
is a feature that helps make rpf modular.
It is easiest to understand the purpose of each function by focusing on the suffixes for the functions above. These suffixes are referenced below. The prefixes often tell us which item model the function pertains to. Below also effectively contains stubs that could be used to start constructing a new item model.
numSpecEach item model (e.g., when created by R) contains
a vector that contains a specification for the item
model. This will contain information such as the
number of categories, factors, and other item-specific
information. I will discuss how this specification is
constructed from R input at a later section of this
document. I assume the numSpec function
merely returns the number of elements of this vector.
For example:
static int
irt_rpf_1dim_lmp_numSpec(const double *spec)
{ return RPF_ISpecCount; }
numParamBased only on information in the item specification, return an integer that indicates how many parameters the item has. For example, the LMP model has a user-specified integer k that determines the number of parameters and appears as the last element in the specification vector.
static int
irt_rpf_1dim_lmp_numParam(const double *spec)
{
int k = spec[RPF_ISpecCount];
return(2+2*k);
}
paramInfoReturns information (in pointers) regarding a single
item parameter based on the item specification *spec,
the index of the item parameter param. This generally assumes
that the item parameters always appear in a particular
order for the item model, which can be determined based
only on the item specification. For instance, a
two-parameter logistic item model will always have a
slope and intercept term, which appear in that order,
and the number of slopes depends on the number of factors.
Returned information may include the type of
parameter **type which is a text description of
the item parameter (e.g., “Slope” or “Intercept”), and
the upper and lower bounds for the item parameter, *upper
and *lower, which may be set to nan("unset") if
bounds are not applicable. Below is only a stub, not the full
function for the LMP model.
static void
irt_rpf_1dim_lmp_paramInfo(const double *spec, const int param,
const char **type, double *upper, double *lower)
{
\\ Code implementing paramInfo goes here
}
probThe traceline (ICC, IRF, etc.) function for the item model.
This function should compute the probability of response
to each category (stored in *out in that order),
based on a fixed vector of item parameters,*param, and
at a single point of the latent space, *th. That is,
if the item model is dichotomous, *out will contain
two values, and if the latent trait has three dimensions,
*th will have three values. In unidimensional models,
*th will have only a single value.
static void
irt_rpf_1dim_lmp_prob(const double *spec,
const double *param, const double *th,
double *out)
{
\\ Code implementing prob goes here
}
logprob_adapter or logprobI assume this function computes the log of the traceline,
which in some cases can be done simply by obtaining
information from the prob function corresponding
to the item model and is what logprob_adapter does.
irt_rpf_logprob_adapter(const double *spec,
const double *param, const double *th,
double *out)
{
\\ Code implementing logprob goes here
}
deriv1This function computes first and second-order complete data derivatives the negative log-likelihood for a single item with respect to item parameters. e.g., for use in optimization at the M-Step of the EM algorithm.
This function assumes complete data on the latent trait as is often the case after the E-Step of the EM algorithm has effectively computed expected counts for each possible category at each quadrature node.
Note that derivatives are at a single point of the latent
trait, *where which is analogous to *th from
the traceline function. This may represent a single
quadrature node, or if complete data were available it
might represent a single respondent's score on the
latent trait.
*weight is a vector that determines how much
weight to give to each response category. If input to
the function were to assume a single respondent, this
vector is effectively an indicator function that will
have a “1” in the proper location indicating the
person's response. In the case of EM, the weight vector
may contain expected counts at this particular quadrature
node for each of the item's possible response categories.
*out will contain the derivatives for use
elsewhere. Note that this vector may not be empty
when entering the function, e.g., in the case that
deriv1 is repeatedly called to compute
derivatives summing across multiple
quadrature points or respondents. From the rpf
documentation, the order of elements in *out is…
“For p parameters, the first p values are the first
derivative and the next p(p+1)/2 columns are the lower
triangle of the second derivative.”
static void
irt_rpf_1dim_lmp_deriv1(const double *spec,
const double *param,
const double *where,
const double *weight, double *out)
{
// Code implementing deriv1 goes here
}
deriv2This function implements derivatives or computations
that may occur once any time derivatives are taken,
regardless of how many quadrature nodes or respondents.
Whereas deriv1 for example, may be called repeatedly
(once for for each quadrature node, deriv2 would
only be called once when computing derivatives across
such quadrature nodes. That is, currently the
rpf.dLL function in R from rpf may do deriv1
multiple times followed by deriv2
before returning derivatives to the user.
static void irt_rpf_1dim_lmp_deriv2(const double *spec,
const double *param,
double *out)
{
\\ Code implementing deriv2 goes here
}
dThetaCompute derivatives of the traceline w.r.t. the latent trait, e.g., for use in computing item information.
These derivative computations happen at a single point
of the latent trait, *where. Derivatives
for each category response function should appear,
in order from least to gretest, in *grad (first-order),
and *hess (second-order).
The vector *dir is a basis vector used in the
case of multidimensional models.
static void irt_rpf_1dim_lmp_dTheta(const double *spec, const double *param,
const double *where, const double *dir,
double *grad, double *hess)
{
\\ Code implementing dTheta goes here
}
rescaleCurrently this function is not yet used and could go unimplemented. However, the intuition behind it was based on Schilling & Bock (2005), equations 22 and 23 for implementation of adaptive quadrature.
static void
irt_rpf_1dim_lmp_rescale(const double *spec, double *param, const int *paramMask,
const double *mean, const double *cov)
{
error("Rescale for LMP model not implemented");
}
The the R/ subfolder, each item model or class of item models has their own associated R file. For example, “grm.R” or lmp.R" for the Graded Response model and the LMP item model.
A function common to each of these files is one that will construct the item specification based on user input. An example for the Graded Response Models is:
rpf.grm <- function(outcomes=2, factors=1, multidimensional=TRUE) {
if (!multidimensional && factors > 1) {
stop("More than 1 dimension must use a multidimensional model")
}
m <- NULL
id <- -1
if (!multidimensional) {
stop("The old parameterization is no longer available")
} else {
id <- rpf.id_of("grm")
m <- new("rpf.mdim.grm",
outcomes=outcomes,
factors=factors)
}
m@spec <- c(id, m@outcomes, m@factors)
m
}
Note that input to the function may be item specific,
however, what the function returns ought to be in a
standard format. Specifically, a list containing the id
of the item model (as determined by rpf.id_of, which
requires a String that matches that added to the end
of librpf_model[]), the number of categories (or
outcomes), and number of factors is required. Additional
optional elements may be appended to the specification
as required by the item model. For example, LMP requires
that k be added to the end of the specification to
determine the order of a polynomial and number of item
parameters, and does not currently have a multidimensional
version.
rpf.lmp <- function(q=1, multidimensional=FALSE) {
if(!(q%%1==0)){
stop("q must be an integer >= 0")
}
if(multidimensional){
stop("Multidimensional LMP model is not yet supported")
}
m <- NULL
id <- -1
id <- rpf.id_of("lmp")
m <- new("rpf.1dim.lmp",
outcomes=2,
factors=1)
m@spec <- c(id, 2, m@factors, q)
m
}
Other functions in such an .R file appear to be
optional, but may include those that specify how random
parameters be generated for the item model, e.g.,
rpf.rparam, or how an item model may be modified
to create a similar item model rpf.modify
Note also that these functions ought to have documentation for the item model that is in a format usable by roxygen2
Any new files created in the previous step ought to be added to the “Collate” statement in the DESCRIPTION file.
A class ought to be added to this file for the item model. This will allow generic functions or wrappers to access item specific C++ code. It is difficult to provide specific recommendations on how to add the class other than to look at other examples. For instance, the following code snippet defines the LMP dichotomous response model, “rpf.lmp.drm”, as a subclass of the unidimensional response model superclass, “rpf.1dim”.
##' Unidimensional logistic function of a monotonic polynomial
##'
##' @export
##' @name Class rpf.1dim.lmp
##' @rdname rpf.1dim.lmp-class
##' @aliases rpf.1dim.lmp-class
##'
setClass("rpf.1dim.lmp", contains='rpf.1dim')
Thus, the above is all that was required to add a class for the LMP item model. Alternative item models may use the multidimensional superclass, “rpf.mdim” or a specialized class for graded response models.
To avoid headaches in debugging, it may be useful
to write code incrementally and test/debug before
continuing to write more. For instance, after
writing stubs for C++ code, see if the code will
compile, perhaps by using R CMD check rpf
or some such from the command line of the OS.
Not all pieces need to be in place before testing
in R. For instance, it may make practical sense
to implement the traceline function in C++, with
suffix, prob, and test via R (see below)
before continuing to write deriv1
and dTheta functions and so on.
Examination of the accompanying R functions,
their documentation in the rpf package,
and their output for well-known item models
may also provide insight into how the underlyling
C++ functions operate.
To test the underlying C++ from R, all
modifications to R files as mentioned in the
previous section ought to be done. However,
only C++ functions that you want to test
need to be implemented. Compile and install
the package for testing, e.g.,
R CMD INSTALL rpf.
Fire up R and begin testing.
library(rpf)
lmp.item<-rpf.lmp(q=2) # create item w/ 5th order polynomial
par<-c(.69,.71,-.5,-8.48,.52,-3.32) # item parameters
theta<-seq(-3,3,.1) # grid for latent trait
## Test the traceline or "prob" C++ function
P<-rpf.prob(lmp.item, par, theta)
## Prettier plots than this are of course possible
plot(theta, P[2,], type="l", ylim=c(0,1), xlab="Theta", ylab="P(Theta)")
The plot above should look like Example 2 from Falk and Cai (in press).
Access to derivatives w.r.t. item parameters and
latent traits is also possible. Note that rpf.dLL
appears to call both deriv1 and deriv2
C++ functions.
## Derivatives of negative log-likelihood at arbitrary point with arbitrary weights
## Rounding only for easy reading for tutorial
round(rpf.dLL(lmp.item, par, theta[1], weight=c(5,7)),2)
## [1] 19.73 -6.11 -58.72 0.11 25.32 0.40 28.35 -2.67 0.83
## [10] -84.40 7.95 1171.18 0.16 -0.02 -0.15 0.11 36.40 -3.43
## [19] -190.32 0.29 36.38 0.58 -0.05 -4.09 0.01 0.22 0.40
For latent traits, with easy ways of testing accuracy against numerical derivatives - at least for a unidimensional model.
## Analytical derivatives "deriv1" followed by "deriv2"
rpf.dTheta(lmp.item, par, where=-.5, dir=1)
## $gradient
## [1] -0.4404183 0.4404183
##
## $hessian
## [1] -0.3160515 0.3160515
## Numerical derivatives
library(numDeriv)
dTheta.wrap<-function(theta, spec, par, cat=1){
rpf.prob(spec, par, theta)[cat]
}
## should match first element of gradient from rpf.dTheta
grad(dTheta.wrap, -.5, spec=lmp.item, par=par)
## [1] -0.4404183
## should match first element of hessian from rpf.dTheta
hessian(dTheta.wrap, -.5, spec=lmp.item, par=par)
## [,1]
## [1,] -0.3160515
One could also test numerical derivatives of the log-likelihood w.r.t. item parameters in a similar manner. However, this may take a little extra work, and note that numerical derivatives may not always work well, especially when testing during estimation and when parameter estimates are far from the MLE.
If any formal testing is done that the user wishes to be reproducible, which is highly encouraged and some may argue is necessary, these formal tests can currently be added to inst/tests/
Forthcoming.
Forthcoming.