Overview

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)

Procedure

Outline

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.

Edit src/libifa-rpf.cpp

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.

numSpec

Each 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; }

numParam

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

paramInfo

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

prob

The 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 logprob

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

deriv1

This 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

}

deriv2

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

dTheta

Compute 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

}

rescale

Currently 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");
}

Edit other, mostly R files

Add R functions to create item specification

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

Edit DESCRIPTION

Any new files created in the previous step ought to be added to the “Collate” statement in the DESCRIPTION file.

Edit R/Classes.R

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.

Putting things together and testing

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.

Installing and testing from R

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.

Add formal testing files

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/

Testing estimation

Forthcoming.

Generating documentation

Forthcoming.