Package 'psychotools'

Anchor Methods for the Detection of Uniform DIF the Rasch Model

Description

The anchor function provides a variety of anchor methods for the detection of uniform differential item functioning (DIF) in the Rasch model between two pre-specified groups. These methods can be divided in an anchor class that determines characteristics of the anchor method and an anchor selection that determines the ranking order of candidate anchor items. The aim of the anchor function is to provide anchor items for DIF testing, e.g. with anchortest.

Usage

anchor(object, ...)
## Default S3 method:
anchor(object, object2,
  class = c("constant", "forward"), select = NULL,
  length = NULL, range = c(0.1, 0.8), ...)
## S3 method for class 'formula'
anchor(formula, data = NULL, subset = NULL,
  na.action = NULL, weights = NULL, model = raschmodel, ...)

Arguments

object, object2	Fitted model objects of class “raschmodel” estimated via conditional maximum likelihood using raschmodel.
...	further arguments passed over to an internal call of anchor.default in the formula method. In the default method, these additional arguments are currently not being used.
class	character. Available anchor classes are the constant anchor class implying a constant anchor length defined by length and the iterative forward anchor class, for an overview see Kopf et al. (2015a).
select	character. Several anchor selection strategies are available: "MTT", "MPT", "MT", "MP", "NST", "AO", "AOP" based on Kopf et al. (2015b) as well as "Gini", "CLF", "GiniT", "CLFT" based on Strobl et al. (2021). The latter four can only be combined with class = "constant" and length = 1. The default is select = "Gini" unless either length > 1 where "MPT" is used or class = "forward" where "MTT" is used. For more details see below.
length	integer. It pre-defines a maximum anchor length. Per default, the forward anchor grows up to the proportion of currently presumed DIF-free items specified in range and the constant anchor class selects one anchor item, unless an explicit limiting number is defined in length by the user.
range	numeric vector of length 2. The first element is the percentage of first anchor candidates to be excluded for consideration when the forward anchor class is used and the second element determines a percentage of currently presumed DIF-free items up to which the anchor from the forward anchor class is allowed to grow.
formula	formula of type y ~ x where y specifies a matrix of dichotomous item responses and x the grouping variable, e.g., gender, for which DIF should be tested for.
data	a data frame containing the variables of the specified formula.
subset	logical expression indicating elements or rows to keep: missing values are taken as false.
na.action	a function which indicates what should happen when the data contain missing values (NAs).
weights	an optional vector of weights (interpreted as case weights).
model	an IRT model fitting function with a suitable itempar method, by default raschmodel.

Details

The anchor methods provided consist of an anchor class that determines characteristics of the anchor and an anchor selection that determines the ranking order of candidate anchor items.

In the constant anchor class, the anchor length is pre-defined by the user within the argument length, defaulting to a length of one. In contrast, the iterative forward class starts with a single anchor item and includes items in the anchor as long as the anchor length is shorter than a certain percentage of the number of items that do not display statistically significant DIF (default: 0.8). Furthermore, a percentage of first anchor candidates is excluded from consideration (default: 0.1) and the user is allowed to set a maximum number of anchor items using the argument length. A detailed description of the anchor classes can be found in Kopf et al. (2015a).

In more recent work Strobl et al. (2021) suggest a simpler yet powerful anchor method based on inequality criteria like the Gini coefficient. A similar approach based on the component loss function (CLF) was suggested by Muthén & Asparouhov (2014). These criteria can be shown to attain their optimium for a single-anchor, thus correponding to a constant class of length 1. Due to the simple structure in combination with good empirical performance the Gini-based selection was made the default in version 0.7-0 of the package.

Both anchor classes require an explicit anchor selection strategy (as opposed to the all-other anchor class which is therefore not included in the function anchor). The anchor selection strategy determines the ranking order of candidate anchor items. In case of two groups, each item $j, j = 1, \ldots, k$ (where $k$ denotes the number of items in the test) obtains a criterion value $c_j$ that is defined by the anchor selection strategy. The ranking order is determined by the rank of the criterion value rank$(c_j)$.

The criterion values $c_j$ for item $j$ from the different anchor selection strategies are provided in the following equations: $d_j$ denotes the difference of the item parameters, $t_j$ the corresponding test statistic, and $p_j$ the resulting p-values. In all cases, the anchor items are given in parentheses. Furthermore, $\mathrm{Gini}(\cdot)$ denotes the Gini inequality index, $\mathrm{CLF}(\cdot)$ the component loss function (sum of square root values), $1(\cdot)$ the indicator function, $\lceil 0.5\cdot k \rceil$ the empirical 50% quantile, and $A_\mathrm{purified}$ the anchor after purification steps. More detailed descriptions are available in Strobl et al. (2021) and Kopf et al. (2015b).

Gini selection (of item parameter differences) by Strobl et al. (2021):

$c_j^\mathrm{Gini} = - \mathrm{Gini} (\{ |d_1(j)|, \ldots, |d_k(j)| \})$

GiniT selection (of test statistics) similar to Strobl et al. (2021):

$c_j^\mathrm{GiniT} = - \mathrm{Gini} (\{ |t_1(j)|, \ldots, |t_k(j)| \})$

CLF selection (of item parameter differences) by Muthén & Asparouhov (2014):

$c_j^\mathrm{CLF} = \mathrm{CLF} (\{ |d_1(j)|, \ldots, |d_k(j)| \})$

CLFT selection (of test statistics) similar to Muthén & Asparouhov (2014):

$c_j^\mathrm{CLFT} = \mathrm{CLF} (\{ |t_1(j)|, \ldots, |t_k(j)| \})$

All-other selection by Woods (2009), here abbreviated AO:

$c_j^\mathrm{AO} = | t_j (\{1,\ldots,k\}\setminus j) |$

All-other purified selection by Wang et al. (2012), here abbreviated AOP:

$c_j^\mathrm{AOP} = | t_j ( A_\mathrm{purified} ) |$

Number of significant threshold selection based on Wang et al. (2004), here abbreviated NST:

$c_j^\mathrm{NST} = \sum_{l \in \{1,\ldots,k\} \setminus j} 1 \left\{ p_j ( \{l\} ) \leq \alpha \right\}) |$

Mean test statistic selection by Shih et al. (2009), here abbreviated MT:

$c_j^\mathrm{MT} = \frac{1}{k-1} \sum_{l \in \{1,\ldots,k\} \setminus j} \left| t_j ( \{l\}) \right|$

Mean p-value selection by Kopf et al. (2015b), here abbreviated MP:

$c_j^\mathrm{MP} = - \frac{1}{k-1} \sum_{l \in \{1,\ldots,k\} \setminus j} p_j ( \{l\} )$

Mean test statistic threshold selection by Kopf et al. (2015b), here abbreviated MTT:

$c_j^\mathrm{MTT} = \sum_{l \in \{1,\ldots,k\} \setminus j} 1 \left\{ \left| t_j ( \{l\} ) \right| > \left( \left| \frac{1}{k-1} \sum_{l \in \{ 1, \ldots, k \} \setminus j} t_j ( \{l\} ) \right| \right)_{\left( \lceil 0.5\cdot k\rceil \right)} \right\}$

Mean p-value threshold selection by Kopf et al. (2015b), here abbreviated MPT:

$c_j^\mathrm{MPT} = - \sum_{l \in \{1,\ldots,k\} \setminus j} 1 \left\{ p_j ( \{l\} ) > \left( \frac{1}{k-1} \sum_{l \in \{ 1, \ldots, k \} \setminus j} p_j ( \{l\} ) \right)_{ \left( \lceil 0.5\cdot k\rceil \right)} \right\}$

Kopf et al. (2015b) recommend to combine the class = "constant" with select = "MPT" and the class = "forward" with select = "MTT", respectively.

The all-other anchor class (that assumes that DIF is balanced i.e. no group has an advantage in the test) is here not considered as explicit anchor selection and, thus, not included in the anchor function (but in the anchortest function). Note that the all-other anchor class requires strong prior knowledge that DIF is balanced.

Value

An object of class anchor, i.e. a list including

anchor_items	the integer index for the selected anchor items
ranking_order	a ranking order (integer index) of the candidate anchor items by their criterion values
criteria	the criterion values obtained in the anchor selection for each item (unsorted)

References

Kopf J, Zeileis A, Strobl C (2015a). A Framework for Anchor Methods and an Iterative Forward Approach for DIF Detection. Applied Psychological Measurement, 39(2), 83–103. doi:10.1177/0146621614544195

Kopf J, Zeileis A, Strobl C (2015b). Anchor Selection Strategies for DIF Analysis: Review, Assessment, and New Approaches. Educational and Psychological Measurement, 75(1), 22–56. doi:10.1177/0013164414529792

Muthén B, Asparouhov T (2014). IRT Studies of Many Groups: The Alignment Method. Frontiers in Psychology, 5, 978. doi:10.3389/fpsyg.2014.00978

Shih CL, Wang WC (2009). Differential Item Functioning Detection Using the Multiple Indicators, Multiple Causes Method with a Pure Short Anchor. Applied Psychological Measurement, 33(3), 184–199.

Strobl C, Kopf J, Kohler L, von Oertzen T, Zeileis A (2021). Anchor Point Selection: Scale Alignment Based on an Inequality Criterion. Applied Psychological Measurement, 45(3), 214–230. doi:10.1177/0146621621990743

Wang WC (2004). Effects of Anchor Item Methods on the Detection of Differential Item Functioning within the Family of Rasch Models. Journal of Experimental Education, 72(3), 221–261.

Wang WC, Shih CL, Sun GW (2012). The DIF-Free-then-DIF Strategy for the Assessment of Differential Item Functioning. Educational and Psychological Measurement, 72(4), 687–708.

Woods C (2009). Empirical Selection of Anchors for Tests of Differential Item Functioning. Applied Psychological Measurement, 33(1), 42–57.

See Also

anchortest

Examples

## Verbal aggression data
data("VerbalAggression", package = "psychotools")

## Gini anchor (Strobl et al. 2021) for gender DIF in the self-to-blame situations
anchor(resp2[, 1:12] ~ gender , data = VerbalAggression)

## alternatively: based on fitted raschmodel objects
raschmodels <- with(VerbalAggression, lapply(levels(gender), function(i) 
  raschmodel(resp2[gender == i, 1:12])))
anchor(raschmodels[[1]], raschmodels[[2]])

if(requireNamespace("multcomp")) {

## four anchor items from constant anchor class using MPT-selection (Kopf et al. 2015b)
anchor(object = raschmodels[[1]], object2 = raschmodels[[2]], 
  class = "constant", select = "MPT", length = 4)

## iterative forward anchor class using MTT-selection (Kopf et al. 2015b)
set.seed(1)
fanchor <- anchor(object = raschmodels[[1]], object2 = raschmodels[[2]],
  class = "forward", select = "MTT", range = c(0.05, 1))
fanchor

## the same using the formula interface
set.seed(1)
fanchor2 <- anchor(resp2[, 1:12] ~ gender , data = VerbalAggression,
  class = "forward", select = "MTT", range = c(0.05, 1))

## criteria really the same?
all.equal(fanchor$criteria, fanchor2$criteria, check.attributes = FALSE)
}

---

Anchor methods for the detection of uniform DIF in the Rasch model

Description

The anchortest function provides a Wald test (see, e.g., Glas, Verhelst, 1995) for the detection of uniform differential item functioning (DIF) in the Rasch model between two pre-specified groups. A variety of anchor methods is available to build a common scale necessary for the comparison of the item parameters in the Rasch model.

Usage

anchortest(object, ...)
## Default S3 method:
anchortest(object, object2,
  class = c("constant", "forward", "all-other", "fixed"), select = NULL,
  test = TRUE, adjust = "none", length = NULL, range = c(0.1, 0.8), ...)
## S3 method for class 'formula'
anchortest(formula, data = NULL, subset = NULL,
  na.action = NULL, weights = NULL, model = raschmodel, ...)

Arguments

object, object2	Fitted model objects of class “raschmodel” estimated via conditional maximum likelihood using raschmodel.
...	further arguments passed over to an internal call of anchor.default in the formula method. In the default method, these additional arguments are currently not being used.
class	character. Available anchor classes are the constant anchor class implying a constant anchor length defined by length, the iterative forward anchor class that iteratively includes items in the anchor and the all-other anchor class, for an overview see Kopf et al. (2015a). Additionally, the class can be fixed, then select needs to be the numeric index of the fixed selected anchor items.
select	character or numeric. Several anchor selection strategies are available, for details see anchor. Alternatively, for class = "fixed", select needs to be the numeric index of the fixed selected anchor items. Defaults are set such that class = "constant" is combined with select = "Gini" while class = "forward" is combined with select = "MTT". And if select is numeric, then class = "fixed" is used.
test	logical. Should the Wald test be returned for the intended anchor method as final DIF test?
adjust	character. Should the final DIF test be adjusted for multiple testing? For the type of adjustment, see summary.glht and p.adjust.
length	integer. It pre-defines a maximum anchor length. Per default, the forward anchor grows up to the proportion of currently presumed DIF-free items specified in range and the constant anchor class selects four anchor items, unless an explicit limiting number is defined in length by the user.
range	numeric vector of length 2. The first element is the percentage of first anchor candidates to be excluded for consideration when the forward anchor class is used and the second element determines a percentage of currently presumed DIF-free items up to which the anchor from the forward anchor class is allowed to grow.
formula	formula of type y ~ x where y specifies a matrix of dichotomous item responses and x the grouping variable, e.g., gender, for which DIF should be tested for.
data	a data frame containing the variables of the specified formula.
subset	logical expression indicating elements or rows to keep: missing values are taken as false.
na.action	a function which indicates what should happen when the data contain missing values (NAs).
weights	an optional vector of weights (interpreted as case weights).
model	an IRT model fitting function with a suitable itempar method, by default raschmodel.

Details

To conduct the Wald test (see, e.g., Glas, Verhelst, 1995) for uniform DIF in the Rasch model, the user needs to specify an anchor method. The anchor methods can be divided in an anchor class that determines characteristics of the anchor method and an anchor selection that determines the ranking order of candidate anchor items.

Explicit anchor selection strategies are used in the constant anchor class and in the iterative forward anchor class, for a detailed description see anchor. Since $k-1$ parameters are free in the estimation, only $k-1$ estimated standard errors result. Thus, the first anchor item obtains no DIF test result and we report $k-1$ test results. This decision is applied only to those methods that rely on an explicit anchor selection strategy.

In the constant anchor class, the anchor length is pre-defined by the user within the argument length. The default is a single anchor item. The iterative forward class starts with a single anchor item and includes items in the anchor as long as the anchor length is shorter than a certain percentage of the number of items that do not display statistically significant DIF. The default proportion is set to 0.8 in the argument range. Alternatively, the user is allowed to set a maximum number of anchor items using the argument length. Both anchor classes require an explicit anchor selection strategy as opposed to the all-other anchor class.

The all-other anchor class is here not considered as explicit anchor selection and, thus, only included in the anchortest function. For the all-other anchor class, the strategy is set to "none", since all items except for the item currently studied for DIF are used as anchor. Thus, no explicit anchor selection strategy is required and we report $k$ test results. Note that the all-other anchor class requires strong prior knowledge that DIF is balanced.

See Strobl et al. (2021) and Kopf et al. (2015ab) for a detailed introduction. For convenience a trivial "fixed" anchor class is provided where the selected anchor is given directly (e.g., as chosen by a practitioner or by some other anchor selection method).

Value

An object of class anchor, i.e. a list including

anchor_items	the anchor items for DIF analysis.
ranking_order	a ranking order of candidate anchor items.
criteria	the criterion values obtained by the respective anchor selection.
anchored_item_parameters	the anchored item parameters using the anchor items.
anchored_covariances	the anchored covariance matrices using the anchor items.
final_tests	the final Wald test for uniform DIF detection if intended.

References

Glas CAW, Verhelst ND (1995). “Testing the Rasch Model.” In Fischer GH, Molenaar IW (eds.), Rasch Models: Foundations, Recent Developments, and Applications, chapter 5. Springer-Verlag, New York.

Kopf J, Zeileis A, Strobl C (2015a). A Framework for Anchor Methods and an Iterative Forward Approach for DIF Detection. Applied Psychological Measurement, 39(2), 83–103. doi:10.1177/0146621614544195

Kopf J, Zeileis A, Strobl C (2015b). Anchor Selection Strategies for DIF Analysis: Review, Assessment, and New Approaches. Educational and Psychological Measurement, 75(1), 22–56. doi:10.1177/0013164414529792

Strobl C, Kopf J, Kohler L, von Oertzen T, Zeileis A (2021). Anchor Point Selection: Scale Alignment Based on an Inequality Criterion. Applied Psychological Measurement, 45(3), 214–230. doi:10.1177/0146621621990743

Wang WC (2004). Effects of Anchor Item Methods on the Detection of Differential Item Functioning within the Family of Rasch Models. Journal of Experimental Education, 72(3), 221–261.

Woods C (2009). Empirical Selection of Anchors for Tests of Differential Item Functioning. Applied Psychological Measurement, 33(1), 42–57.

See Also

anchor

Examples

if(requireNamespace("multcomp")) {

o <- options(digits = 4)

## Verbal aggression data
data("VerbalAggression", package = "psychotools")

## Rasch model for the self-to-blame situations; gender DIF test
raschmodels <- with(VerbalAggression, lapply(levels(gender), function(i) 
  raschmodel(resp2[gender == i, 1:12])))

## single anchor from Gini selection (default)
gini1 <- anchortest(object = raschmodels[[1]], object2 = raschmodels[[2]])
gini1
summary(gini1)

## four anchor items from constant anchor class using MPT selection
const1 <- anchortest(object = raschmodels[[1]], object2 = raschmodels[[2]],
  class = "constant", select = "MPT", length = 4)
const1
summary(const1)

## iterative forward anchor class using MTT selection
set.seed(1)
forw1 <- anchortest(object = raschmodels[[1]], object2 = raschmodels[[2]], 
  class = "forward", select = "MTT", test = TRUE,
  adjust = "none", range = c(0.05,1))
forw1

## DIF test with fixed given anchor (arbitrarily selected to be items 1 and 2)
anchortest(object = raschmodels[[1]], object2 = raschmodels[[2]], select = 1:2)

options(digits = o$digits)
}

---

Coercing Item Response Data

Description

Coercing "itemresp" data objects to other classes.

Usage

## S3 method for class 'itemresp'
as.list(x, items = NULL, mscale = TRUE, df = FALSE, ...)

Arguments

x	an object of class "itemresp".
items	character, integer, or logical for subsetting the items.
mscale	logical. Should the measurement scale labels be used for creating factor levels? If FALSE, the values 0, 1, ... are used.
df	logical. Should a data frame of factors be returned? If FALSE, a plain list of factors is returned.
...	currently not used.

Details

The as.list method coerces item response data to a list (or data frame) of factors with factor levels either taken from the mscale(x) or as the values 0, 1, ....

The as.data.frame method returns a data frame with a single column of class "itemresp".

Furthermore, as.matrix, as.integer, as.double all return a matrix with the item responses coded as values 0, 1, ...

The as.character method simply calls format.itemresp.

is.itemresp can be used to check wether a given object is of class "itemresp".

See Also

itemresp

Examples

## item responses from binary matrix
x <- cbind(c(1, 0, 1, 0), c(1, 0, 0, 0), c(0, 1, 1, 1))
xi <- itemresp(x)
## change mscale
mscale(xi) <- c("-", "+")
xi

## coercion to list of factors with levels taken from mscale
as.list(xi)
## same but levels taken as integers 0, 1
as.list(xi, mscale = FALSE)
## only for first two items
as.list(xi, items = 1:2)
## result as data.frame
as.list(xi, df = TRUE)

## data frame with single itemresp column
as.data.frame(xi)

## integer matrix
as.matrix(xi)

## character vector
as.character(xi)

## check class of xi
is.itemresp(xi)

---

Bradley-Terry Model Fitting Function

Description

btmodel is a basic fitting function for simple Bradley-Terry models.

Usage

btmodel(y, weights = NULL, type = c("loglin", "logit"), ref = NULL,
  undecided = NULL, position = NULL, start = NULL, vcov = TRUE, estfun =
  FALSE, ...)

Arguments

y	paircomp object with the response.
weights	an optional vector of weights (interpreted as case weights).
type	character. Should an auxiliary log-linear Poisson model or logistic binomial be employed for estimation? The latter is not available if undecided effects are estimated.
ref	character or numeric. Which object parameter should be the reference category, i.e., constrained to zero?
undecided	logical. Should an undecided parameter be estimated?
position	logical. Should a position effect be estimated?
start	numeric. Starting values when calling glm.fit.
vcov	logical. Should the estimated variance-covariance be included in the fitted model object?
estfun	logical. Should the empirical estimating functions (score/gradient contributions) be included in the fitted model object?
...	further arguments passed to functions.

Details

btmodel provides a basic fitting function for Bradley-Terry models, intended as a building block for fitting Bradley-Terry trees and Bradley-Terry mixtures in the psychotree package, respectively. While btmodel is intended for individual paired-comparison data, the eba package provides functions for aggregate data.

btmodel returns an object of class "btmodel" for which several basic methods are available, including print, plot, summary, coef, vcov, logLik, estfun and worth.

Value

btmodel returns an S3 object of class "btmodel", i.e., a list with components as follows.

y	paircomp object with the response
coefficients	estimated parameters on log-scale (without the first parameter which is always constrained to be 0),
vcov	covariance matrix of the parameters in the model,
loglik	log-likelihood of the fitted model,
df	number of estimated parameters,
weights	the weights used (if any),
n	number of observations (with non-zero weights),
type	character for model type (see above),
ref	character for reference category (see above),
undecided	logical for estimation of undecided parameter (see above),
position	logical for estimation of position effect (see above),
labels	character labels of the objects compared,
estfun	empirical estimating function (also known as scores or gradient contributions).

See Also

pcmodel, gpcmodel, rsmodel, raschmodel, nplmodel, the eba package

Examples

o <- options(digits = 4)

## data
data("GermanParties2009", package = "psychotools")

## Bradley-Terry model
bt <- btmodel(GermanParties2009$preference)
summary(bt)
plot(bt)

options(digits = o$digits)

---

Generic Conspiracist Beliefs Scale (2016 Data)

Description

Responses of 2449 persons to 15 five-point likert-rated items (0 = disagree to 4 = agree) measuring belief in conspiracy theories as well as responses on 2 covariates.

Usage

data("ConspiracistBeliefs2016", package = "psychotools")

Format

A data frame containing 2449 observations on 3 variables.

resp

Item response matrix with 15 items (see details below).

area

Factor coding the area one lived in as a child ("rural", "suburban", "urban").

gender

Factor coding gender ("male", "female", "other").

Details

The Open Source Psychometrics Project published this dataset collected online in 2016. Persons responded to the Generic Conspiracist Beliefs (GCB) Scale (Brotherton, French & Pickering, 2013) as well as other additional questions primarily for personal amusement. At the end of the test but before the results were displayed, users were asked if they would allow their responses to be saved for research. Only users who agreed are part of this dataset. Individuals with age lower than 13 years were not recorded. Moreover, two persons stating their age to be 5555 years or higher as well as 44 persons with missing data in area or gender were excluded from this dataset. The 15 items of the GCB Scale are:

Q1:	The government is involved in the murder of innocent citizens and/or well-known public figures, and keeps this a secret.
Q2:	The power held by heads of state is second to that of small unknown groups who really control world politics.
Q3:	Secret organizations communicate with extraterrestrials, but keep this fact from the public.
Q4:	The spread of certain viruses and/or diseases is the result of the deliberate, concealed efforts of some organization.
Q5:	Groups of scientists manipulate, fabricate, or suppress evidence in order to deceive the public.
Q6:	The government permits or perpetrates acts of terrorism on its own soil, disguising its involvement.
Q7:	A small, secret group of people is responsible for making all major world decisions, such as going to war.
Q8:	Evidence of alien contact is being concealed from the public.
Q9:	Technology with mind-control capacities is used on people without their knowledge.
Q10:	New and advanced technology which would harm current industry is being suppressed.
Q11:	The government uses people as patsies to hide its involvement in criminal activity.
Q12:	Certain significant events have been the result of the activity of a small group who secretly manipulate world events.
Q13:	Some UFO sightings and rumors are planned or staged in order to distract the public from real alien contact.
Q14:	Experiments involving new drugs or technologies are routinely carried out on the public without their knowledge or consent.
Q15:	A lot of important information is deliberately concealed from the public out of self-interest.

Additional information can be found online (see below) via inspecting the codebook contained in ‘GCBS.zip’.

Source

https://openpsychometrics.org/_rawdata/.

References

Brotherton R, French CC, Pickering AD (2013). Measuring Belief in Conspiracy Theories: The Generic Conspiracist Beliefs Scale. Frontiers in Psychology, 4, 279.

Open Source Psychometrics Project (2016). Data From: The Generic Conspiracist Beliefs Scale [Dataset]. Retrieved from https://openpsychometrics.org/_rawdata/.

See Also

gpcmodel

Examples

## overview
data("ConspiracistBeliefs2016", package = "psychotools")
str(ConspiracistBeliefs2016)

## response
plot(itemresp(ConspiracistBeliefs2016$resp))
## covariates
summary(ConspiracistBeliefs2016[, -1])

---

Extract/Set Covariates

Description

A generic function for extracting/setting covariates for an object.

Usage

covariates(object, ...)
  covariates(object) <- value

Arguments

object	an object.
...	arguments passed to methods.
value	an object.

Examples

## method for "paircomp" data
pc <- paircomp(rbind(
  c(1,  1,  1), # a > b, a > c, b > c
  c(1,  1, -1), # a > b, a > c, b < c
  c(1, -1, -1), # a > b, a < c, b < c
  c(1,  1,  1)))
covariates(pc)
covariates(pc) <- data.frame(foo = factor(c(1, 2, 2), labels = c("foo", "bar")))
covariates(pc)

---

Response Curve Plots for IRT Models

Description

Base graphics plotting function for response curve plot visualization of IRT models.

Usage

curveplot(object, ref = NULL, items = NULL, names = NULL,
    layout = NULL, xlim = NULL, ylim = c(0, 1), col = NULL,
    lty = NULL, main = NULL, xlab = "Latent trait",
    ylab = "Probability", add = FALSE, ...)

Arguments

object	a fitted model object of class "raschmodel", "rsmodel", "pcmodel", "nplmodel" or "gpcmodel".
ref	argument passed over to internal calls of predict.
items	character or numeric, specifying the items for which response curves should be visualized.
names	character, specifying labels for the items.
layout	matrix, specifying how the response curve plots of different items should be arranged.
xlim, ylim	numeric, specifying the x and y axis limits.
col	character, specifying the colors of the response curve lines. The length of col should be the maximum number of available categories.
lty	numeric, specifying the line type of the response curve lines. The length of lty should either be one or the maximum number of available categories. In the first case, a single line type is used for all category response curves. In the latter case, separate line types for each category response curve are used.
main	character, specifying the overall title of the plot.
xlab, ylab	character, specifying the x and y axis labels.
add	logical. If TRUE, new response curves are added to an existing plot. Only possible when a single item is visualized.
...	further arguments passed to internal calls of matplot.

Details

The response curve plot visualization illustrates the predicted probabilities as a function of the ability parameter $\theta$ under a certain IRT model. This type of visualization is sometimes also called item/category operating curves or item/category characteristic curves.

See Also

regionplot, profileplot, infoplot, piplot

Examples

## load verbal aggression data
data("VerbalAggression", package = "psychotools")

## fit Rasch, rating scale and partial credit model to verbal aggression data
rmmod <- raschmodel(VerbalAggression$resp2)
rsmod <- rsmodel(VerbalAggression$resp)
pcmod <- pcmodel(VerbalAggression$resp)

## curve plots of the dichotomous RM
plot(rmmod, type = "curves")

## curve plots under the RSM for the first six items of the data set
plot(rsmod, type = "curves", items = 1:6)

## curve plots under the PCM for the first six items of the data set with
## custom labels
plot(pcmod, type = "curves", items = 1:6, names = paste("Item", 1:6))

## compare the predicted probabilities under the RSM and the PCM for a single
## item
plot(rsmod, type = "curves", item = 1)
plot(pcmod, type = "curves", item = 1, lty = 2, add = TRUE)
legend(x = "topleft", y = 1.0, legend = c("RSM", "PCM"), lty = 1:2, bty = "n")


if(requireNamespace("mirt")) {
## fit 2PL and generaliced partial credit model to verbal aggression data
twoplmod <- nplmodel(VerbalAggression$resp2)
gpcmod <- gpcmodel(VerbalAggression$resp)

## curve plots of the dichotomous 2PL
plot(twoplmod, type = "curves", xlim = c(-6, 6))

## curve plots under the GPCM for the first six items of the data set
plot(gpcmod, type = "curves", items = 1:6, xlim = c(-6, 6))
}

---

Extract Discrimination Parameters of Item Response Models

Description

A class and generic function for representing and extracting the discrimination parameters of a given item response model.

Usage

discrpar(object, ...)
  ## S3 method for class 'raschmodel'
discrpar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'rsmodel'
discrpar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'pcmodel'
discrpar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'nplmodel'
discrpar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'gpcmodel'
discrpar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)

Arguments

object	a fitted model object whose discrimination parameters should be extracted.
ref	a restriction to be used. Not used for models estimated via CML as the discrimination parameters are fixed to 1 in raschmodels, rsmodels and pcmodels. For models estimated via MML (nplmodels and gpcmodels), the parameters are by default identified via the distributional parameters of the person parameters (mean and variance of a normal distribution). Nevertheless, a restriction on the ratio scale can be applied.
alias	logical. If TRUE (the default), the aliased parameters are included in the return vector (and in the variance-covariance matrix if vcov = TRUE). If FALSE, these parameters are removed. For raschmodels, rsmodels and pcmodels where all discrimination parameters are fixed to 1, this means that an empty numeric vector and an empty variance-covariance matrix is returned if alias is FALSE.
vcov	logical. If TRUE (the default), the variance-covariance matrix of the discrimination parameters is attached as attribute vcov.
...	further arguments which are currently not used.

Details

discrpar is both, a class to represent discrimination parameters of item response models as well as a generic function. The generic function can be used to extract the discrimination parameters of a given item response model.

For objects of class discrpar, several methods to standard generic functions exist: print, coef, vcov. coef and vcov can be used to extract the discrimination parameters and their variance-covariance matrix without additional attributes.

Value

A named vector with discrimination parameters of class discrpar and additional attributes model (the model name), ref (the items or parameters used as restriction/for normalization), alias (either TRUE or a named numeric vector with the aliased parameters not included in the return value), and vcov (the estimated and adjusted variance-covariance matrix).

See Also

personpar, itempar, threshpar, guesspar, upperpar

Examples

o <- options(digits = 4)

## load verbal aggression data
data("VerbalAggression", package = "psychotools")

## fit Rasch model to verbal aggression data
rmod <- raschmodel(VerbalAggression$resp2)

## extract the discrimination parameters
dp1 <- discrpar(rmod)

## extract the standard errors
sqrt(diag(vcov(dp1)))

if(requireNamespace("mirt")) {
## fit 2PL to verbal aggression data
twoplmod <- nplmodel(VerbalAggression$resp2)

## extract the discrimination parameters
dp2 <- discrpar(twoplmod)

## this time with the first discrimination parameter being the reference
discrpar(twoplmod, ref = 1)

## extract the standard errors
sqrt(diag(vcov(dp2)))
}

options(digits = o$digits)

---

Calculation of the Elementary Symmetric Functions and Their Derivatives

Description

Calculation of elementary_symmetric_functions (ESFs), their first and, in the case of dichotomous items, second derivatives with sum or difference algorithm for the Rasch, rating scale and partial credit model.

Usage

elementary_symmetric_functions(par, order = 0L, log = TRUE,
  diff = FALSE, engine = NULL)

Arguments

par	numeric vector or a list. Either a vector of item difficulty parameters of dichotomous items (Rasch model) or a list of item-category parameters of polytomous items (rating scale and partial credit model).
order	integer between 0 and 2, specifying up to which derivative the ESFs should be calculated. Please note, second order derivatives are currently only possible for dichtomous items in an R implementation engine == "R".
log	logical. Are the parameters given in par on log scale? Primarily used for internal recursive calls of elementary_symmetric_functions.
diff	logical. Should the first and second derivatives (if requested) of the ESFs calculated with sum (FALSE) or difference algorithm (TRUE).
engine	character, either "C" or "R". If the former, a C implementation is used to calculcate the ESFs and their derivatives, otherwise ("R") pure R code is used.

Details

Depending on the type of par, the elementary symmetric functions for dichotomous (par is a numeric vector) or polytomous items (par is a list) are calculated.

For dichotomous items, the summation and difference algorithm published in Liou (1994) is used. For calculating the second order derivatives, the equations proposed by Jansens (1984) are employed.

For polytomous items, the summation and difference algorithm published by Fischer and Pococny (1994) is used (see also Fischer and Pococny, 1995).

Value

elementary_symmetric_function returns a list of length 1 + order.

If order = 0, then the first (and only) element is a numeric vector with the ESFs of order 0 to the maximum score possible with the given parameters.

If order = 1, the second element of the list contains a matrix, with the rows corresponding to the possible scores and the columns corresponding to the derivatives with respect to the i-th parameter of par.

For dichotomous items and order = 2, the third element of the list contains an array with the second derivatives with respect to every possible combination of two parameters given in par. The rows of the individual matrices still correspond to the possibles scores (orders) starting from zero.

References

Liou M (1994). More on the Computation of Higher-Order Derivatives of the Elementary Symmetric Functions in the Rasch Model. Applied Psychological Measurement, 18, 53–62.

Jansen PGW (1984). Computing the Second-Order Derivatives of the Symmetric Functions in the Rasch Model. Kwantitatieve Methoden, 13, 131–147.

Fischer GH, and Ponocny I (1994). An Extension of the Partial Credit Model with an Application to the Measurement of Change. Psychometrika, 59(2), 177–192.

Fischer GH, and Ponocny I (1995). “Extended Rating Scale and Partial Credit Models for Assessing Change.” In Fischer GH, and Molenaar IW (eds.). Rasch Models: Foundations, Recent Developments, and Applications.

Examples

## zero and first order derivatives of 100 dichotomous items
 di <- rnorm(100)
 system.time(esfC <- elementary_symmetric_functions(di, order = 1))
 
 ## again with R implementation
 system.time(esfR <- elementary_symmetric_functions(di, order = 1,
 engine = "R"))

 ## are the results equal?
 all.equal(esfC, esfR)


 ## calculate zero and first order elementary symmetric functions
 ## for 10 polytomous items with three categories each.
 pi <- split(rnorm(20), rep(1:10, each = 2))
 x <- elementary_symmetric_functions(pi)

 ## use difference algorithm instead and compare results
 y <- elementary_symmetric_functions(pi, diff = TRUE)
 all.equal(x, y)

---

Popularity of First Names

Description

Preferences of 192 respondents choosing among six boys names with respect to their popularity.

Usage

data("FirstNames")

Format

A data frame containing 192 observations on 11 variables.

preference

Paired comparison of class paircomp. All 15 pairwise choices among six boys names: Tim, Lucas, Michael, Robin, Benedikt, and Julius.

ordered.pref

Ordered paired comparison of class paircomp. Same as preference, but within-pair order is recognized.

gender

Factor coding gender.

age

Integer. Age of the respondents in years.

education

Ordered factor. Level of education: 1 Hauptschule with degree (Secondary General School), 2 and 3 Realschule without and with degree (Intermediate Secondary School), 4 and 5 Gymnasium without and with degree (High School), 6 and 7 Studium without and with degree (University).

children

Integer. Number of children.

state

Factor. State of Germany where participant grew up.

state.reg

Factor. The region (south, north-west, east) each state belongs to.

fname

Factor. Participant's fist name(s). (Umlaute in Jörg and Jürgen have been transliterated to Joerg and Juergen for portability of the data.)

interviewer

Factor. Interviewer id.

gender.int

Factor coding interviewer's gender.

Details

A survey was conducted at the Department of Psychology, Universität Tübingen, in June 2009. The sample was stratified by gender and age (younger versus older than 30 years) with 48 participants in each group. The interviewers were Psychology Master's students who collected the data for course credits.

Participants were presented with 15 pairs of boys names in random order. On each trial, their task was to choose the name they would rather give to their own child. The pairs of boys names were read to the participants one at a time. A given participant compared each pair in one order only, hence the NA's in ordered.pref.

The names were selected to fall within the upper (Tim, Lucas), mid (Michael, Robin) and lower (Benedikt, Julius) range of the top 100 of the most popular boys names in Germany in the years from 1990 to 1999 (https://www.beliebte-vornamen.de/3778-1990er-jahre.htm). The names have either front (e, i) or back (o, u) vowels in the stressed syllables. Phonology of the name and attractiveness of a person have been shown to be related (Perfors, 2004; Hartung et al., 2009).

References

Hartung F, Klenovsak D, Santiago dos Santos L, Strobl C, Zaefferer D (2009). Are Tims Hot and Toms Not? Probing the Effect of Sound Symbolism on Perception of Facial Attractiveness. Presented at the 31th Annual Meeting of the Cognitive Science Society, July 27–August 1, Amsterdam, The Netherlands.

Perfors A (2004). What's in a Name? The Effect of Sound Symbolism on Perception of Facial Attractiveness. Presented at the 26th Annual Meeting of the Cognitive Science Society, August 5–7, Chicago, USA.

See Also

paircomp

Examples

data("FirstNames", package = "psychotools")
summary(FirstNames$preference)
covariates(FirstNames$preference)

---

Choice among German Political Parties

Description

Preferences of 192 respondents choosing among five German political parties and abstention from voting.

Usage

data("GermanParties2009")

Format

A data frame containing 192 observations on 6 variables.

preference

Paired comparison of class paircomp. All 15 pairwise choices among five German parties and abstention from voting.

ordered.pref

Ordered paired comparison of class paircomp. Same as preference, but within-pair order is recognized.

gender

Factor coding gender.

age

Integer. Age of the respondents in years.

education

Ordered factor. Level of education: 1 no degree, 2 Hauptschule (Secondary General School), 3 Realschule (Intermediate Secondary School), 4 Gymnasium (High School), 5 Studium (University)

crisis

Factor. Do you feel affected by the economic crisis?

interviewer

Factor. Interviewer id.

Details

A survey was conducted at the Department of Psychology, Universität Tübingen, in June 2009, three months before the German election. The sample was stratified by gender and age (younger versus older than 30 years) with 48 participants in each group.

The parties to be compared were Die Linke (socialists), Die Grünen (ecologists), SPD (social democrats), CDU/CSU (conservatives), and FDP (liberals). In addition, there was the option of abstaining from voting (coded as none).

Participants were presented with 15 pairs of options in random order. On each trial, their task was to choose the party they would rather vote for at an election for the German parliament. A given participant compared each pair in one order only, hence the NA's in ordered.pref.

In order to minimize response biases, the pairs of options were read to the participants one at a time. Participants made their choices by crossing either “First Option” or “Second Option” on an anonymous response sheet.

The interviewers were Psychology Master's students who collected the data for course credits. Since they mainly interviewed people they knew, the results are not representative of the political opinions in Germany. As far as the winner of the survey (Die Grünen) is concerned, however, the results agree with the outcome of the election for the Tübingen voters.

The results of the election on September 27, 2009 (number of so-called Zweitstimmen in percent) were:

	Germany	Tübingen
Die Linke	11.9	8.5
Die Grünen	10.7	27.9
SPD	23.0	21.1
CDU/CSU	33.8	23.0
FDP	14.6	13.9
Others	6.0	5.7

The voter turnout was 70.8 percent in Germany and 80.5 percent in Tübingen.

See Also

paircomp

Examples

data("GermanParties2009", package = "psychotools")
summary(GermanParties2009$preference)

---

Generalized Partial Credit Model Fitting Function

Description

gpcmodel is a basic fitting function for generalized partial credit models providing a wrapper around mirt and multipleGroup relying on marginal maximum likelihood (MML) estimation via the standard EM algorithm.

Usage

gpcmodel(y, weights = NULL, impact = NULL, type = c("GPCM", "PCM"),
  grouppars = FALSE, vcov = TRUE, nullcats = "downcode", 
  start = NULL, method = "BFGS", maxit = 500, reltol = 1e-5, ...)

Arguments

y	item response object that can be coerced (via as.matrix) to a numeric matrix with scores 0, 1, ... Typically, either already a matrix, data frame, or dedicated object of class itemresp.
weights	an optional vector of weights (interpreted as case weights)

.

impact	an optional factor allowing for grouping the subjects (rows). If specified, a multiple-group model is fitted to account for impact (see details below). By default, no impact is modelled, i.e., a single-group model is used.
type	character string, specifying the type of model to be estimated. In addition to the default GPCM (generalized partial credit model) it is also possible to estimate a standard PCM (partial credit model) by marginal maximum likelihood (MML).
grouppars	logical. Should the estimated distributional group parameters of a multiple group model be included in the model parameters?
vcov	logical or character specifying the type of variance-covariance matrix (if any) computed for the final model. The default vcov = TRUE corresponds to vcov = "Oakes", see mirt for further options. If set to vcov = FALSE (or vcov = "none"), vcov() will return a matrix of NAs only.
nullcats	character string, specifying how items with null categories (i.e., categories not observed) should be treated. Currently only "downcode" is available, i.e., all categories above a null category are shifted down to close the observed gap(s).
start	an optional vector or list of starting values (see examples below).
method, maxit, reltol	control parameters for the optimizer employed by mirt for the EM algorithm.
...	further arguments passed to mirt or multipleGroup, respectively.

Details

gpcmodel provides a basic fitting function for generalized partial credit models (GPCMs) providing a wrapper around mirt (and multipleGroup, respectively) relying on MML estimation via the standard EM algorithm (Bock & Aitkin, 1981). Models are estimated under the slope/intercept parametrization, see e.g. Chalmers (2012). The probability of person $i$ falling into category $x_{ij}$ of item $j$ out of all categories $p_{j}$ is modelled as:

$P(X_{ij} = x_{ij}|\theta_{i},a_{j},\boldsymbol{d_{j}}) = \frac{\exp{(x_{ij}a_{j}\theta_{i} + d_{jx_{ij}})}}{\displaystyle\sum_{k = 0} ^ {p_{j}}\exp{(ka_{j}\theta_{i} + d_{jk})}}$

Note that all $d_{j0}$ are fixed at 0. A reparametrization of the intercepts to the classical IRT parametrization, see e.g. Muraki (1992), is provided via threshpar.

If an optional impact variable is supplied, a multiple-group model of the following form is being fitted: Item parameters are fixed to be equal across the whole sample. For the first group of the impact variable the person parameters are fixed to follow the standard normal distribution. In the remaining impact groups, the distributional parameters (mean and variance of a normal distribution) of the person parameters are estimated freely. See e.g. Baker & Kim (2004, Chapter 11), Debelak & Strobl (2019), or Schneider et al. (2022) for further details. To improve convergence of the model fitting algorithm, the first level of the impact variable should always correspond to the largest group. If this is not the case, levels are re-ordered internally.

If grouppars is set to TRUE the freely estimated distributional group parameters (if any) are returned as part of the model parameters.

Instead of the default GPCM, a standard partial credit model (PCM) can also be estimated via MML by setting type = "PCM". In this case all slopes are restricted to be equal across all items.

gpcmodel returns an object of class "gpcmodel" for which several basic methods are available, including print, plot, summary, coef, vcov, logLik, estfun, discrpar, itempar, threshpar, and personpar.

Value

gpcmodel returns an S3 object of class "gpcmodel", i.e., a list of the following components:

coefficients	estimated model parameters in slope/intercept parametrization,
vcov	covariance matrix of the model parameters,
data	modified data, used for model-fitting, i.e., centralized so that the first category is zero for all items, treated null categories as specified via argument "nullcats" and without observations with zero weight,
items	logical vector of length ncol(y), indicating which items were used during estimation,
categories	list of length ncol(y), containing integer vectors starting from one to the number of categories minus one per item,
n	number of observations (with non-zero weights),
n_org	original number of observations in y,
weights	the weights used (if any),
na	logical indicating whether the data contain NAs,
nullcats	currently always NULL as eventual items with null categories are handled via "downcode",
impact	either NULL or the supplied impact variable with the levels reordered in decreasing order (if this has not been the case prior to fitting the model),
loglik	log-likelihood of the fitted model,
df	number of estimated (more precisely, returned) model parameters,
code	convergence code from mirt,
iterations	number of iterations used by mirt,
reltol	convergence threshold passed to mirt,
grouppars	the logical grouppars value,
type	the type of model restriction specified,
mirt	the mirt object fitted internally,
call	original function call.

References

Baker FB, Kim SH (2004). Item Response Theory: Parameter Estimation Techniques. Chapman & Hall/CRC, Boca Raton.

Bock RD, Aitkin M (1981). Marginal Maximum Likelihood Estimation of Item Parameters: Application of an EM Algorithm. Psychometrika, 46(4), 443–459.

Chalmers RP (2012). mirt: A Multidimensional Item Response Theory Package for the R Environment. Journal of Statistical Software, 48(6), 1–29. doi:10.18637/jss.v048.i06

Debelak R, Strobl C (2019). Investigating Measurement Invariance by Means of Parameter Instability Tests for 2PL and 3PL Models. Educational and Psychological Measurement, 79(2), 385–398. doi:10.1177/0013164418777784

Muraki E (1992). A Generalized Partial Credit Model: Application of an EM Algorithm. Applied Psychological Measurement, 16(2), 159–176.

Schneider L, Strobl C, Zeileis A, Debelak R (2022). An R Toolbox for Score-Based Measurement Invariance Tests in IRT Models. Behavior Research Methods, forthcoming. doi:10.3758/s13428-021-01689-0

See Also

pcmodel, rsmodel, nplmodel, raschmodel, btmodel

Examples

if(requireNamespace("mirt")) {

o <- options(digits = 4)

## mathematics 101 exam results
data("MathExam14W", package = "psychotools")

## generalized partial credit model
gpcm <- gpcmodel(y = MathExam14W$credit)
summary(gpcm)

## how to specify starting values as a vector of model parameters
st <- coef(gpcm)
gpcm <- gpcmodel(y = MathExam14W$credit, start = st)
## or a list containing a vector of slopes and a list of intercept vectors
## itemwise
set.seed(0)
st <- list(a = rlnorm(13, 0, 0.0625), d = replicate(13, rnorm(2, 0, 1), FALSE))
gpcm <- gpcmodel(y = MathExam14W$credit, start = st)

## visualizations
plot(gpcm, type = "profile")
plot(gpcm, type = "regions")
plot(gpcm, type = "piplot")
plot(gpcm, type = "curves", xlim = c(-6, 6))
plot(gpcm, type = "information", xlim = c(-6, 6))
## visualizing the IRT parametrization
plot(gpcm, type = "curves", xlim = c(-6, 6), items = 1)
abline(v = threshpar(gpcm)[[1]])
abline(v = itempar(gpcm)[1], lty = 2)

options(digits = o$digits)
}

---

Extract Guessing Parameters of Item Response Models

Description

A class and generic function for representing and extracting the so-called guessing parameters of a given item response model.

Usage

guesspar(object, ...)
  ## S3 method for class 'raschmodel'
guesspar(object, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'rsmodel'
guesspar(object, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'pcmodel'
guesspar(object, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'nplmodel'
guesspar(object, alias = TRUE, logit = FALSE, vcov = TRUE, ...)
  ## S3 method for class 'gpcmodel'
guesspar(object, alias = TRUE, vcov = TRUE, ...)

Arguments

object	a fitted model object whose guessing parameters should be extracted.
alias	logical. If TRUE (the default), the aliased parameters are included in the return vector (and in the variance-covariance matrix if vcov = TRUE). If FALSE, these parameters are removed. For raschmodels, rsmodels, pcmodels and gpcmodels, where all guessing parameters are fixed to 0, this means that an empty numeric vector and an empty variance-covariace matrix is returned if alias is FALSE.
logit	logical. If a nplmodel of type "3PL" or "4PL" model has been fit, the guessing parameters were estimated on the logit scale. If logit = FALSE, these estimates and the variance-covariance (if requested) are retransformed using the logistic function and the delta method.
vcov	logical. If TRUE (the default), the variance-covariance matrix of the guessing parameters is attached as attribute vcov.
...	further arguments which are currently not used.

Details

guesspar is both, a class to represent guessing parameters of item response models as well as a generic function. The generic function can be used to extract the guessing parameters of a given item response model.

For objects of class guesspar, several methods to standard generic functions exist: print, coef, vcov. coef and vcov can be used to extract the guessing parameters and their variance-covariance matrix without additional attributes.

Value

A named vector with guessing parameters of class guesspar and additional attributes model (the model name), alias (either TRUE or a named numeric vector with the aliased parameters not included in the return value), logit (indicating whether the estimates are on the logit scale or not), and vcov (the estimated and adjusted variance-covariance matrix).

See Also

personpar, itempar, threshpar, discrpar, upperpar

Examples

if(requireNamespace("mirt")) {

o <- options(digits = 3)

## load simulated data
data("Sim3PL", package = "psychotools")

## fit 2PL to data simulated under the 3PL
twoplmod <- nplmodel(Sim3PL$resp)

## extract the guessing parameters (all fixed at 0)
gp1 <- guesspar(twoplmod)

## fit 3PL to data simulated under the 3PL
threeplmod <- nplmodel(Sim3PL$resp, type = "3PL")

## extract the guessing parameters
gp2 <- guesspar(threeplmod)

## extract the standard errors
sqrt(diag(vcov(gp2)))

## extract the guessing parameters on the logit scale
gp2_logit <- guesspar(threeplmod, logit = TRUE)

## along with the delta transformed standard errors
sqrt(diag(vcov(gp2_logit)))

options(digits = o$digits)
}

---

Information Plots for IRT Models

Description

Base graphics plotting function for information plot visualization of IRT models.

Usage

infoplot(object, what = c("categories", "items", "test"),
    ref = NULL, items = NULL, names = NULL, layout = NULL, xlim = NULL,
    ylim = NULL, col = NULL, lty = NULL, lwd = NULL, main = NULL, legend = TRUE,
    xlab = "Latent trait", ylab = "Information", add = FALSE, ...)

Arguments

object	a fitted model object of class "raschmodel", "rsmodel", "pcmodel", "nplmodel" or "gpcmodel".
what	character, specifying the type of information to visualize.
ref	argument passed over to internal calls of predict.
items	character or numeric, specifying the items for which information curves should be visualized.
names	character, specifying labels for the items.
layout	matrix, specifying how the item or category information curves of different items should be arranged. If null and what is set to "items", the item information curves are overlayed within a single plot.
xlim, ylim	numeric, specifying the x and y axis limits.
col	character, specifying the colors of the test, item or category information curves.
lty	numeric, specifying the line type of the information curves.
lwd	numeric, specifying the line width of the information curves.
main	character, specifying the overall title of the plot.
legend	logical, specifying if a legend is drawn when multiple item information curves are overlayed. The labels in the legend correspond to the item names (which can be specified in the argument names).
xlab, ylab	character, specifying the x and y axis labels.
add	logical. If TRUE, new information curves are added to an existing plot. Only possible for a test or a single item information curve.
...	further arguments passed to internal calls of matplot.

Details

The information plot visualization illustrates the test, item or category information as a function of the ability parameter $\theta$ under a certain IRT model. Further details on the computation of the displayed information can be found on the help page of the function predict.pcmodel.

See Also

curveplot, regionplot, profileplot, piplot

Examples

## load verbal aggression data
data("VerbalAggression", package = "psychotools")

## fit Rasch and partial credit model to verbal aggression data
rmmod <- raschmodel(VerbalAggression$resp2)
pcmod <- pcmodel(VerbalAggression$resp)

## category information plots for all items under the dichotomous RM
plot(rmmod, type = "information", what = "categories")

## category information plots for all items under the PCM
plot(pcmod, type = "information", what = "categories")

## overlayed item information plots for the first six items of the
## data set under the PCM
plot(pcmod, type = "information", what = "items", items = 1:6)

## a comparison of the item information for the first six items under the
## dichotomous RM and the PCM
plot(pcmod, type = "information", what = "items", items = 1:6,
  xlim = c(-5, 5))
plot(rmmod, type = "information", what = "items", items = 1:6,
  lty = 2, add = TRUE)
legend(x = "topright", legend = c("PCM", "RM"), lty = 1:2, bty = "n")

## a comparison of the test information based on all items of the
## data set under the dichotomous RM and the PCM
plot(pcmod, type = "information", what = "test", items = 1:6, xlim = c(-5, 5))
plot(rmmod, type = "information", what = "test", items = 1:6, lty = 2,
  add = TRUE)
legend(x = "topright", legend = c("PCM", "RM"), lty = 1:2, bty = "n")

if(requireNamespace("mirt")) {
## fit 2PL to verbal aggression data
twoplmod <- nplmodel(VerbalAggression$resp2)

## category information plots for all items under the dichotomous 2PL
plot(twoplmod, type = "information", what = "categories")
}

---

Extract Item Parameters of Item Response Models

Description

A class and generic function for representing and extracting the item parameters of a given item response model.

Usage

itempar(object, ...)
  ## S3 method for class 'raschmodel'
itempar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'rsmodel'
itempar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'pcmodel'
itempar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'nplmodel'
itempar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'gpcmodel'
itempar(object, ref = NULL, alias = TRUE, vcov = TRUE, ...)
  ## S3 method for class 'btmodel'
itempar(object, ref = NULL, alias = TRUE, vcov = TRUE, log = FALSE, ...)

Arguments

object	a fitted model or tree object whose item parameters should be extracted.
ref	a vector of labels or position indices of item parameters or a contrast matrix which should be used as restriction/for normalization. If NULL (the default) for all models except models estimated via MML, all items are used (sum zero restriction). For models estimated via MML (nplmodels and gpcmodels), the parameters are by default identified via the distributional parameters of the person parameters (mean and variance of a normal distribution). Nevertheless, a restriction on the interval scale can be applied.
alias	logical. If TRUE (the default), the aliased parameter is included in the return vector (and in the variance-covariance matrix if vcov = TRUE). If FALSE, it is removed. If the restriction given in ref depends on several parameters, the first parameter of the restriction specified is (arbitrarily) chosen to be removed if alias is FALSE.
vcov	logical. If TRUE (the default), the (transformed) variance-covariance matrix of the item parameters is attached as attribute vcov. If FALSE, an NA-matrix is attached.
log	logical. Whether to return the estimated model parameters on the logit (TRUE) or preference scale (FALSE).
...	further arguments which are currently not used.

Details

itempar is both, a class to represent item parameters of item response models as well as a generic function. The generic function can be used to extract the item parameters of a given item response model.

For Rasch models and n-parameter logistic models, itempar returns the estimated item difficulty parameters $\hat{\beta}_{j}$ under the restriction specified in argument ref. For rating scale models, itempar returns computed item location parameters $\hat{\beta}_{j}$ under the restriction specified in argument ref. These are computed from the estimated item-specific parameters $\hat{\xi}_{j}$ (who mark the location of the first category of an item on the latent theta axis). For partial credit models and generalized partial credit models, itempar returns ‘mean’ absolute item threshold parameters, $\hat{\beta}_{j} = \frac{1}{p_{j}} \sum_{k = 1}^{p_{j}}\hat{\delta}_{jk}$, i.e., a single parameter per item is returned which results as the mean of the absolute item threshold parameters $\hat{\delta}_{jk}$ of this item. Based upon these ‘mean’ absolute item threshold parameters $\hat{\beta}_{j}$, the restriction specified in argument ref is applied. For all models, the variance-covariance matrix of the returned item parameters is adjusted according to the multivariate delta rule.

For objects of class itempar, several methods to standard generic functions exist: print, coef, vcov. coef and vcov can be used to extract the estimated calculated item parameters and their variance-covariance matrix without additional attributes. Based on this Wald tests or confidence intervals can be easily computed, e.g., via confint.

Two-sample item-wise Wald tests for DIF in the item parameters can be carried out using the function anchortest.

Value

A named vector with item parameters of class itempar and additional attributes model (the model name), ref (the items or parameters used as restriction/for normalization), alias (either FALSE or a named character vector with the removed aliased parameter, and vcov (the adjusted covariance matrix of the estimates if vcov = TRUE or an NA-matrix otherwise).

See Also

personpar, threshpar, discrpar, guesspar, upperpar

Examples

o <- options(digits = 4)

## load verbal aggression data
data("VerbalAggression", package = "psychotools")

## fit a Rasch model to dichotomized verbal aggression data
raschmod <- raschmodel(VerbalAggression$resp2)

## extract item parameters with sum zero or use last two items as anchor
ip1 <- itempar(raschmod)
ip2a <- itempar(raschmod, ref = 23:24) # with position indices
ip2b <- itempar(raschmod, ref = c("S4WantShout", "S4DoShout")) # with item label

ip1
ip2a

all.equal(ip2a, ip2b)

## extract vcov
vc1 <- vcov(ip1)
vc2 <- vcov(ip2a)

## adjusted standard errors,
## smaller with more items used as anchors
sqrt(diag(vc1))
sqrt(diag(vc2))

## Wald confidence intervals
confint(ip1)
confint(ip2a)

options(digits = o$digits)

---

Data Structure for Item Response Data

Description

A class for representing data from questionnaires along with methods for many generic functions.

Usage

itemresp(data, mscale = NULL, labels = NULL, names = NULL)

Arguments

data	matrix or data frame. A matrix or data frame with integer values or factors where the rows correspond to subjects and the columns to items. See below for details.
mscale	integer or character. A list of vectors (either integer or character) giving the measurement scale. See below for details. By default guessed from data.
labels	character. A vector of character labels for the items. By default, the column names of data are used or, if these are not available, the string "item" along with numbers 1, 2, ... is used.
names	character. A vector of names (or IDs) for the subjects. By default, no subject names are used.

Details

itemresp is designed for item response data of $n$ subjects for $k$ items.

The item responses should be coded in a matrix data with $n$ rows (subjects) and $k$ columns (items). Alternatively, data can be a data frame with $n$ rows (subjects) and $k$ variables (items), which can be either factors or integer valued vectors.

mscale provides the underlying measurement scale either as integer or character vector(s). If all items are measured on the same scale, mscale can be a vector. Alternatively, it can be provided as a named list of vectors for each item. If the list contains one unnamed element, this element will be used as the measurement scale for items that have not been named. Integers or characters not present in mscale but in data will be replaced by NA. All items must be measured with at least 2 categories. By default, mscale is set to the full range of observed values for all integer items (see example below) and the corresponding levels for all factor items in data.

Methods to standard generic functions include: str, length (number of subjects), dim (number of subjects and items), is.na (only TRUE if all item responses are NA for a subject), print (see print.itemresp for details), summary and plot (see summary.itemresp for details), subsetting via [ and subset (see subset.itemresp for details), is.itemresp and various coercion functions to other classes (see as.list.itemresp for details).

Extracting/replacing properties is available through: labels for the item labels, mscale for the measurement scale, names for subject names/IDs.

Value

itemresp returns an object of class "itemresp" which is a matrix (data transformed to integers 0, 1, ...) plus an attribute "mscale" as a named list for each item (after being checked and potentially suitably coerced or transformed to all integer or all character).

See Also

print.itemresp, summary.itemresp, as.list.itemresp, subset.itemresp

Examples

## binary responses to three items, coded as matrix
x <- cbind(c(1, 0, 1, 0), c(1, 0, 0, 0), c(0, 1, 1, 1))
## transformed to itemresp object
xi <- itemresp(x)

## printing (see also ?print.itemresp)
print(xi)
print(xi, labels = TRUE)

## subsetting/indexing (see also ?subset.itemresp)
xi[2]
xi[c(TRUE, TRUE, FALSE, FALSE)]
subset(xi, items = 1:2)
dim(xi)
length(xi)

## summary/visualization (see also ?summary.itemresp)
summary(xi)
plot(xi)

## query/set measurement scale labels
## extract mscale (tries to collapse to vector)
mscale(xi)
## extract as list
mscale(xi, simplify = FALSE)
## replacement by list
mscale(xi) <- list(item1 = c("no", "yes"),
  item2 = c("nay", "yae"), item3 = c("-", "+"))
xi
mscale(xi)
## replacement with partially named list plus default
mscale(xi) <- list(item1 = c("n", "y"), 0:1)
mscale(xi)
## replacement by vector (if number of categories constant)
mscale(xi) <- c("-", "+")
mscale(xi, simplify = FALSE)

## query/set item labels and subject names
labels(xi)
labels(xi) <- c("i1", "i2", "i3")
names(xi)
names(xi) <- c("John", "Joan", "Jen", "Jim")
print(xi, labels = TRUE)

## coercion (see also ?as.list.itemresp)
## to integer matrix
as.matrix(xi)
## to data frame with single itemresp column
as.data.frame(xi)
## to list of factors
as.list(xi)
## to data frame with factors
as.list(xi, df = TRUE)


## polytomous responses with missing values and unequal number of
## categories in a data frame
d <- data.frame(
  q1 = c(-2, 1, -1, 0, NA, 1, NA),
  q2 = c(3, 5, 2, 5, NA, 2, 3),
  q3 = factor(c(1, 2, 1, 2, NA, 3, 2), levels = 1:3,
    labels = c("disagree", "neutral", "agree")))
di <- itemresp(d)
di

## auto-completion of mscale: full range (-2, ..., 2) for q1, starting
## from smallest observed (negative) value (-2) to the same (positive)
## value (2), full (positive) range for q2, starting from smallest
## observed value (2) to largest observed value (5), missing category of
## 4 is detected, for q3 given factor levels are used
mscale(di)

## set mscale for q2 and add category 1, q1 and q3 are auto-completed:
di <- itemresp(d, mscale = list(q2 = 1:5))

## is.na.itemresp - only true for observation 5 (all missing)
is.na(di)

## illustration for larger data set
data("VerbalAggression", package = "psychotools")
r <- itemresp(VerbalAggression$resp[, 1:12])
str(r)
head(r)
plot(r)
summary(r)
prop.table(summary(r), 1)

## dichotomize response
r2 <- r
mscale(r2) <- c(0, 1, 1)
plot(r2)

## transform to "likert" package
if(require("likert")) {
lik <- likert(as.data.frame(as.list(r)))
lik
plot(lik)
}

---

Set Labels

Description

A generic function for setting labels for an object.

Usage

labels(object) <- value

Arguments

object	an object.
value	an object.

Examples

## method for "paircomp" data
pc <- paircomp(rbind(
  c(1,  1,  1), # a > b, a > c, b > c
  c(1,  1, -1), # a > b, a > c, b < c
  c(1, -1, -1), # a > b, a < c, b < c
  c(1,  1,  1)))
labels(pc)
labels(pc) <- c("ah", "be", "ce")
pc

---

Mathematics 101 Exam Results

Description

Responses of 729 students to 13 items in a written exam of introductory mathematics along with several covariates.

Usage

data("MathExam14W")

Format

A data frame containing 729 observations on 9 variables.

solved

Item response matrix (of class itemresp) with values 1/0 coding solved correctly/other.

credits

Item response matrix (of class itemresp) with values 2/1/0 coding solved correctly/incorrectly/not attempted.

nsolved

Integer. The number of items solved correctly.

tests

Integer. The number of online test exercises solved correctly prior to the written exam.

gender

Factor indicating gender.

study

Factor indicating two different types of business/economics degrees. Either the 3-year bachelor program (571) or the 4-year diploma program (155).

semester

Integer. The number of semesters enrolled in the given university program.

attempt

Factor. The number of times the course/exam has been attempted (including the current attempt).

group

Factor indicating whether the students were in the first or second batch (with somewhat different items) in the exam.

Details

The data provides individual end-term exam results from a Mathematics 101 course for first-year business and economics students at Universität Innsbruck. The format of the course comprised biweekly online tests (26 numeric exercises, conducted in OpenOLAT) and a written exam at the end of the semester (13 single-choice exercises with five answer alternatives). The course covers basics of analysis, linear algebra, financial mathematics, and probability calculus (where the latter is not assessed in this exam).

In this exam, 729 students participated (out of 941 registered in the course). To avoid cheating, all students received items with essentially the same questions but different numbers (using the exams infrastructure of Zeileis et al. 2014). Also, due to the large number of students two groups of students had to be formed which received partially different items. The items which differed (namely 1, 5, 6, 7, 8, 9, 11, 12) varied in the setup/story, but not in the mathematical skills needed to solve the exercises. Prior to the exam, the students could select themselves either into the first group (early in the morning) or the second group (starting immediately after the end of the first group).

Correctly solved items yield 100 percent of the associated points. Items without correct solution can either be unanswered (0 percent) or receive an incorrect answer (minus 25 percent) to discourage random guessing. In the examples below, the items are mostly only considered as binary. Typically, students with 8 out of 13 correct answers passed the course.

Source

Department of Statistics, Universität Innsbruck

References

Zeileis A, Umlauf N, Leisch F (2014). Flexible Generation of E-Learning Exams in R: Moodle Quizzes, OLAT Assessments, and Beyond. Journal of Statistical Software, 58(1), 1–36. doi:10.18637/jss.v058.i01

See Also

itemresp, raschmodel, pcmodel, anchortest

Examples

## load data and exclude extreme scorers
data("MathExam14W", package = "psychotools")
MathExam14W <- transform(MathExam14W,
  points = 2 * nsolved - 0.5 * rowSums(credits == 1)
)
me <- subset(MathExam14W, nsolved > 0 & nsolved < 13)


## item response data:
## solved (correct/other) or credits (correct/incorrect/not attempted)
par(mfrow = c(1, 2))
plot(me$solved)
plot(me$credits)

## PCA
pr <- prcomp(me$solved, scale = TRUE)
names(pr$sdev) <- 1:10
plot(pr, main = "", xlab = "Number of components")
biplot(pr, col = c("transparent", "black"), main = "",
  xlim = c(-0.065, 0.005), ylim = c(-0.04, 0.065))


## points achieved (and 50% threshold)
par(mfrow = c(1, 1))
hist(MathExam14W$points, breaks = -4:13 * 2 + 0.5,
  col = "lightgray", main = "", xlab = "Points")
abline(v = 12.5, lwd = 2, col = 2)


## Rasch and partial credit model
ram <- raschmodel(me$solved)
pcm <- pcmodel(me$credits)

## various types of graphics displays
plot(ram, type = "profile")
plot(pcm, type = "profile", add = TRUE, col = "blue")
plot(ram, type = "piplot")
plot(pcm, type = "piplot")
plot(ram, type = "region")
plot(pcm, type = "region")
plot(ram, type = "curves")
plot(pcm, type = "curves")



## test for differential item function with automatic anchoring
## passing vs. not passing students
at1 <- anchortest(solved ~ factor(nsolved <= 7), data = me,
  adjust = "single-step")
at1
plot(at1$final_tests)
## -> "good" students discriminate somewhat more
## (quad/payflow/lagrange are slightly more difficult)

## group 1 vs. group 2
at2 <- anchortest(solved ~ group, data = me, adjust = "single-step")
at2
plot(at2$final_tests)
## -> quad/payflow/planning easier for group 1
## -> hesse slightly easier for group 2

## bring out differences between groups 1 and 2
## by (anchored) item difficulty profiles
ram1 <- raschmodel(subset(me, group == "1")$solved)
ram2 <- raschmodel(subset(me, group == "2")$solved)
plot(ram1, parg = list(ref = at2$anchor_items), ylim = c(-2, 3))
plot(ram2, parg = list(ref = at2$anchor_items), add = TRUE, col = "blue")
legend("topleft", c("Group 1", "Group 2"), pch = 21,
  pt.bg = c("lightgray", "blue"), bty = "n")

---

Memory Deficits in Psychiatric Patients

Description

Response frequencies of 96 patients who took part in a pair-clustering experiment to assess their memory deficits.

Usage

data("MemoryDeficits")

Format

A data frame containing 576 observations on 7 variables.

ID

Participant ID.

group

Factor with four levels specifying patient or control group of participant.

trial

Trial number from 1 to 6.

E1

Number of pairs recalled adjacently.

E2

Number of pairs recalled non-adjacently.

E3

Number of single pair members recalled.

E4

Number of non-recalled pairs.

Details

Riefer, Knapp, Batchelder, Bamber and Manifold (2002) report a study on memory deficits in schizophrenic (n = 29) and organic alcoholic (n = 21) patients who were compared to two matched control groups (n = 25, n = 21). Participants were presented with 20 pairs of semantically related words. In a later memory test, they freely recalled the presented words. This procedure was repeated for a total of six study and test trials. Responses were classified into four categories: both words in a pair are recalled adjacently (E1) or non-adjacently (E2), one word in a pair is recalled (E3), neither word in a pair is recalled (E4).

Source

The data were made available by William H. Batchelder.

References

Riefer DM, Knapp BR, Batchelder WH, Bamber D, Manifold V (2002). Cognitive Psychometrics: Assessing Storage and Retrieval Deficits in Special Populations with Multinomial Processing Tree Models. Psychological Assessment, 14, 184–201.

Examples

data("MemoryDeficits", package = "psychotools")
aggregate(cbind(E1, E2, E3, E4) ~ trial + group, MemoryDeficits, sum)

---

Multinomial Processing Tree (MPT) Model Fitting Function

Description

mptmodel is a basic fitting function for multinomial processing tree (MPT) models.

Usage

mptmodel(y, weights = NULL, spec, treeid = NULL,
  optimargs = list(control = list(reltol = .Machine$double.eps^(1/1.2),
                                  maxit = 1000),
                   init = NULL),
  start = NULL, vcov = TRUE, estfun = FALSE, ...)

Arguments

y	matrix of response frequencies.
weights	an optional vector of weights (interpreted as case weights).
spec	an object of class mptspec: typically result of a call to mptspec. A symbolic description of the model to be fitted.
treeid	a factor that identifies each tree in a joint multinomial model.
optimargs	a list of arguments passed to the optimization function (optim).
start	a vector of starting values for the parameter estimates between zero and one.
vcov	logical. Should the estimated variance-covariance be included in the fitted model object?
estfun	logical. Should the empirical estimating functions (score/gradient contributions) be included in the fitted model object?
...	further arguments passed to functions.

Details

mptmodel provides a basic fitting function for multinomial processing tree (MPT) models, intended as a building block for fitting MPT trees in the psychotree package. While mptmodel is intended for individual response frequencies, the mpt package provides functions for aggregate data.

MPT models are specified using the mptspec function. See the documentation in the mpt package for details.

mptmodel returns an object of class "mptmodel" for which several basic methods are available, including print, plot, summary, coef, vcov, logLik, estfun and predict.

Value

mptmodel returns an S3 object of class "mptmodel", i.e., a list with components as follows:

y	a matrix with the response frequencies,
coefficients	estimated parameters (for extraction, the coef function is preferred),
loglik	log-likelihood of the fitted model,
npar	number of estimated parameters,
weights	the weights used (if any),
nobs	number of observations (with non-zero weights),
ysum	the aggregate response frequencies,
fitted, goodness.of.fit, ...	see mpt in the mpt package.

See Also

btmodel, pcmodel, gpcmodel, rsmodel, raschmodel, nplmodel, mptspec, the mpt package

Examples

o <- options(digits = 4)

## data
data("SourceMonitoring", package = "psychotools")

## source-monitoring MPT model
mpt1 <- mptmodel(SourceMonitoring$y, spec = mptspec("SourceMon"))
summary(mpt1)
plot(mpt1)

options(digits = o$digits)

---

Extract/Replace Measurement Scale

Description

Generic functions for extracting and replacing the measurement scale from an object.

Usage

mscale(object, ...)
  mscale(object) <- value

Arguments

object	an object.
...	arguments passed to methods.
value	an object describing the measurement scale.

Examples

## methods for "paircomp" data
pc <- paircomp(rbind(
  c(2,  1,  0),
  c(1,  1, -1),
  c(1, -2, -1),
  c(0,  0,  0)))
pc

## extract
mscale(pc)

## replace (collapse to >/=/< scale)
mscale(pc) <- sign(mscale(pc))
pc


## similar for "itemresp" data
ir <- itemresp(cbind(
  c(-1, 0, 1, 1, 0),
  c(0, 1, 2, 1, 2),
  c(1, 2, 1, 1, 3)))
ir

## extract
mscale(ir)

## replace (single scale for all items)
mscale(ir) <- 1:3
ir

---

Parametric Logistic Model (n-PL) Fitting Function

Description

nplmodel is a basic fitting function for n-PL type parametric logistic IRT models (2PL, 3PL, 3PLu, 4PL, Rasch/1PL), providing a wrapper around mirt and multipleGroup relying on marginal maximum likelihood (MML) estimation via the standard EM algorithm.

Usage

nplmodel(y, weights = NULL, impact = NULL,
  type = c("2PL", "3PL", "3PLu", "4PL", "1PL", "RM"),
  grouppars = FALSE, vcov = TRUE, 
  start = NULL, method = "BFGS", maxit = 500, reltol = 1e-5, ...)

Arguments

y	item response object that can be coerced (via as.matrix) to a numeric matrix with scores 0, 1. Typically, either already a matrix, data frame, or dedicated object of class itemresp.
weights	an optional vector of weights (interpreted as case weights).
impact	an optional factor allowing for grouping the subjects (rows). If specified, a multiple-group model is fitted to account for impact (see details below). By default, no impact is modelled, i.e., a single-group model is used.
type	character string, specifying the type of parametric logistic IRT model to be estimated (see details below).
grouppars	logical. Should the estimated distributional group parameters of a multiple group model be included in the model parameters?
vcov	logical or character specifying the type of variance-covariance matrix (if any) computed for the final model. The default vcov = TRUE corresponds to vcov = "Oakes", see mirt for further options. If set to vcov = FALSE (or vcov = "none"), vcov() will return a matrix of NAs only.
start	an optional vector or list of starting values (see examples below).
method, maxit, reltol	control parameters for the optimizer employed by mirt for the EM algorithm.
...	further arguments passed to mirt or multipleGroup, respectively.

Details

nplmodel (plmodel for backward compatibility with earlier psychotools versions) provides a basic fitting function for n-PL type parametric logistic IRT models (2PL, 3PL, 3PLu, 4PL, Rasch/1PL) providing a wrapper around mirt and multipleGroup relying on MML estimation via the standard EM algorithm (Bock & Aitkin, 1981). Models are estimated under the slope/intercept parametrization, see e.g. Chalmers (2012). The probability of person $i$ ‘solving’ item $j$ is modelled as:

$P(X_{ij} = 1|\theta_{i},a_{j},d_{j},g_{j},u_{j}) = g_{j} + \frac{(u_{j} - g_{j})}{1 + \exp{(-(a_{j}\theta_{i} + d_{j}))}}$

A reparametrization of the intercepts to the classical IRT parametrization, $b_{j} = -\frac{d_{j}}{a_{j}}$, is provided via the corresponding itempar method.

If an optional impact variable is supplied, a multiple-group model of the following form is being fitted: Item parameters are fixed to be equal across the whole sample. For the first group of the impact variable the person parameters are fixed to follow the standard normal distribution. In the remaining impact groups, the distributional parameters (mean and variance of a normal distribution) of the person parameters are estimated freely. See e.g. Baker & Kim (2004, Chapter 11), Debelak & Strobl (2019), or Schneider et al. (2022) for further details. To improve convergence of the model fitting algorithm, the first level of the impact variable should always correspond to the largest group. If this is not the case, levels are re-ordered internally.

If grouppars is set to TRUE the freely estimated distributional group parameters (if any) are returned as part of the model parameters.

By default, type is set to "2PL". Therefore, all so-called guessing parameters are fixed at 0 and all upper asymptotes are fixed at 1. "3PL" results in all upper asymptotes being fixed at 1 and "3PLu" results in all all guessing parameters being fixed at 0. "4PL" results in a full estimated model as specified above. Finally, if type is set to "1PL" (or equivalently "RM"), an MML-estimated Rasch model is being fitted. This means that all slopes are restricted to be equal across all items, all guessing parameters are fixed at 0 and all upper asymptotes are fixed at 1.

Note that internally, the so-called guessing parameters and upper asymptotes are estimated on the logit scale (see also mirt). Therefore, most of the basic methods below include a logit argument, which can be set to TRUE or FALSE allowing for a retransformation of the estimates and their variance-covariance matrix (if requested) using the logistic function and the delta method if logit = FALSE.

nplmodel returns an object of class "nplmodel" for which several basic methods are available, including print, plot, summary, coef, vcov, logLik, estfun, discrpar, itempar, threshpar, guesspar, upperpar, and personpar.

Finally, if type is set to "1PL", a Rasch model is estimated. Here, a common slope parameter is estimated for all items, whereas the person parameters are assumed to follow a standard normal distribution. Please note that this variant of the Rasch model differs from the one used by mirt, which sets all slope parameters to 1, and estimates the variance of the person parameters instead. Both variants are mathematically equivalent and hence should lead to the same intercept parameter estimates. For numerical reasons, nplmodel and mirt can lead to slightly different item parameter estimates, though, under their respective default settings, in particular when some items are very easy or very difficult and the common slope parameter is large. A distinct advantage of the variant used by nplmodel is that it allows a direct comparison of the slope and intercept parameters with that estimated in more complex IRT models, such as the 2PL model.

Value

nplmodel returns an S3 object of class "nplmodel", i.e., a list of the following components:

coefficients	estimated model parameters in slope/intercept parametrization,
vcov	covariance matrix of the model parameters,
data	modified data, used for model-fitting, i.e., without observations with zero weight,
items	logical vector of length ncol(y), indicating which items were used during estimation,
n	number of observations (with non-zero weights),
n_org	original number of observations in y,
weights	the weights used (if any),
na	logical indicating whether the data contain NAs,
impact	either NULL or the supplied impact variable with the levels reordered in decreasing order (if this has not been the case prior to fitting the model),
loglik	log-likelihood of the fitted model,
df	number of estimated (more precisely, returned) model parameters,
code	convergence code from mirt,
iterations	number of iterations used by mirt,
reltol	convergence threshold passed to mirt,
grouppars	the logical grouppars value,
type	the type of model restriction specified,
mirt	the mirt object fitted internally,
call	original function call.

References

Baker FB, Kim SH (2004). Item Response Theory: Parameter Estimation Techniques. Chapman & Hall/CRC, Boca Raton.

Bock RD, Aitkin M (1981). Marginal Maximum Likelihood Estimation of Item Parameters: Application of an EM Algorithm. Psychometrika, 46(4), 443–459.

Chalmers RP (2012). mirt: A Multidimensional Item Response Theory Package for the R Environment. Journal of Statistical Software, 48(6), 1–29. doi:10.18637/jss.v048.i06

Debelak R, Strobl C (2019). Investigating Measurement Invariance by Means of Parameter Instability Tests for 2PL and 3PL Models. Educational and Psychological Measurement, 79(2), 385–398. doi:10.1177/0013164418777784

Schneider L, Strobl C, Zeileis A, Debelak R (2022). An R Toolbox for Score-Based Measurement Invariance Tests in IRT Models. Behavior Research Methods, forthcoming. doi:10.3758/s13428-021-01689-0

See Also

raschmodel, gpcmodel, rsmodel, pcmodel, btmodel

Examples

if(requireNamespace("mirt")) {

o <- options(digits = 4)

## mathematics 101 exam results
data("MathExam14W", package = "psychotools")

## 2PL
twopl <- nplmodel(y = MathExam14W$solved)
summary(twopl)

## how to specify starting values as a vector of model parameters
st <- coef(twopl)
twopl <- nplmodel(y = MathExam14W$solved, start = st)
## or a list containing a vector of slopes and a vector of intercepts
set.seed(0)
st <- list(a = rlnorm(13, 0, 0.0625), d = rnorm(13, 0, 1))
twopl <- nplmodel(y = MathExam14W$solved, start = st)

## visualizations
plot(twopl, type = "profile")
plot(twopl, type = "regions")
plot(twopl, type = "piplot")
plot(twopl, type = "curves", xlim = c(-6, 6))
plot(twopl, type = "information", xlim = c(-6, 6))
## visualizing the IRT parametrization
plot(twopl, type = "curves", xlim = c(-6, 6), items = 1)
abline(v = itempar(twopl)[1])
abline(h = 0.5, lty = 2)

## 2PL accounting for gender impact
table(MathExam14W$gender)
mtwopl <- nplmodel(y = MathExam14W$solved, impact = MathExam14W$gender,
  grouppars = TRUE)
summary(mtwopl)
plot(mtwopl, type = "piplot")
## specifying starting values as a vector of model parameters, note that in
## this example impact is being modelled and therefore grouppars must be TRUE
## to get all model parameters
st <- coef(mtwopl)
mtwopl <- nplmodel(y = MathExam14W$solved, impact = MathExam14W$gender,
  start = st)
## or a list containing a vector of slopes, a vector of intercepts and a vector
## of means and a vector of variances as the distributional group parameters
set.seed(1)
st <- list(a = rlnorm(13, 0, 0.0625), d = rnorm(13, 0, 1), m = 0, v = 1)
mtwopl <- nplmodel(y = MathExam14W$solved, impact = MathExam14W$gender,
  start = st)

## MML estimated Rasch model (1PL)
rm <- nplmodel(y = MathExam14W$solved, type = "1PL")
summary(rm)

options(digits = o$digits)
}

---

Pair Clustering Data in Klauer (2006)

Description

Response frequencies of 63 participants who took part in a pair-clustering experiment.

Usage

data("PairClustering")

Format

A data frame containing 126 observations on 8 variables.

ID

Participant ID.

trial

Trial number, 1 or 2.

E1

Number of pairs recalled adjacently.

E2

Number of pairs recalled non-adjacently.

E3

Number of single pair members recalled.

E4

Number of non-recalled pairs.

F1

Number of recalled singleton words.

F2

Number of non-recalled singleton words.

Details

Klauer (2006) reports a pair-clustering experiment with 63 participants, who were presented with ten pairs of related words and five unrelated singleton words. In a later memory test, they freely recalled the presented words. This procedure was repeated for two study and test trials. For pairs, responses were classified into four categories: both words in a pair are recalled adjacently (E1) or non-adjacently (E2), one word in a pair is recalled (E3), neither word in a pair is recalled (E4); for singletons, into two categories: word recalled (F1), word not recalled (F2).

Source

Stahl C, Klauer KC (2007). HMMTree: A Computer Program for Latent-Class Hierarchical Multinomial Processing Tree Models. Behavior Research Methods, 39, 267–273.

References

Klauer KC (2006). Hierarchical Multinomial Processing Tree Models: A Latent-Class Approach. Psychometrika, 71, 1–31.

Examples

data("PairClustering", package = "psychotools")
aggregate(cbind(E1, E2, E3, E4, F1, F2) ~ trial, PairClustering, sum)

---

Data Structure for Paired Comparisons

Description

A class for representing data from paired comparison experiments along with methods for many generic functions.

Usage

paircomp(data,
    labels = NULL, mscale = NULL, ordered = FALSE, covariates = NULL)

Arguments

data	matrix. A matrix with integer values where the rows correspond to subjects and the columns to paired comparisons between objects. See below for details.
labels	character. A vector of character labels for the objects. By default a suitable number of letters is used.
mscale	integer. A vector of integers giving the measurement scale. See below for details. By default guessed from data.
ordered	logical. Does data contain both orderings of each comparison?
covariates	data.frame. An optional data.frame with object covariates, i.e., it must have the same number of rows as the length of labels. May be NULL (default).

Details

paircomp is designed for holding paired comparisons of $k$ objects measured for $n$ subjects.

The comparisons should be coded in an integer matrix data with $n$ rows (subjects) and $k \choose 2$ columns (unless ordered = TRUE, see below). The columns must be ordered so that objects are sequentially compared with all previous objects, i.e.: 1:2, 1:3, 2:3, 1:4, 2:4, 3:4, etc. Each column represents the results of a comparison for two particular objects. Positive values signal that the first object was preferred, negative values that the second was preferred, zero signals no preference. Larger absolute values signal stronger preference.

mscale provides the underlying measurement scale. It must be a symmetric sequence of integers of type (-i):i where i must be at least 1. However, it may exclude 0 (i.e., forced choice).

If ordered = TRUE, the order of comparison matters and thus data is assumed to have twice as many columns. The second half of columns then corresponds to the comparisons 2:1, 3:1, 3:2, 4:1, 4:2, 4:3, etc.

Value

paircomp returns an object of class "paircomp" which is a matrix (essentially data) with all remaining arguments of paircomp as attributes (after being checked and potentially suitably coerced or transformed).

See Also

subset.paircomp, print.paircomp

Examples

## a simple paired comparison
pc <- paircomp(rbind(
  c(1,  1,  1), # a > b, a > c, b > c
  c(1,  1, -1), # a > b, a > c, b < c
  c(1, -1, -1), # a > b, a < c, b < c
  c(1,  1,  1)))

## basic methods
pc
str(pc)
summary(pc)
pc[2:3]
c(pc[2], pc[c(1, 4)])

## methods to extract/set attributes
labels(pc)
labels(pc) <- c("ah", "be", "ce")
pc
mscale(pc)
covariates(pc)
covariates(pc) <- data.frame(foo = factor(c(1, 2, 2), labels = c("foo", "bar")))
covariates(pc)
names(pc)
names(pc) <- LETTERS[1:4]
pc

## reorder() and subset() both select a subset of
## objects and/or reorders the objects
reorder(pc, c("ce", "ah"))


## include paircomp object in a data.frame
## (i.e., with subject covariates)
dat <- data.frame(
  x = rnorm(4),
  y = factor(c(1, 2, 1, 1), labels = c("hansi", "beppi")))
dat$pc <- pc
dat


## formatting with long(er) labels and extended scale
pc2 <- paircomp(rbind(
  c(4,  1,  0),
  c(1,  2, -1),
  c(1, -2, -1),
  c(0,  0,  -3)),
  labels = c("Nordrhein-Westfalen", "Schleswig-Holstein", "Baden-Wuerttemberg"))
## default: abbreviate
print(pc2)
print(pc2, abbreviate = FALSE)
print(pc2, abbreviate = FALSE, width = FALSE)


## paired comparisons with object covariates
pc3 <- paircomp(rbind(
  c(2,  1,  0),
  c(1,  1, -1),
  c(1, -2, -1),
  c(0,  0,  0)),
  labels = c("New York", "Rio", "Tokyo"),
  covariates = data.frame(hemisphere = factor(c(1, 2, 1), labels = c("North", "South"))))
covariates(pc3)

---

Partial Credit Model Fitting Function

Description

pcmodel is a basic fitting function for partial credit models.

Usage

pcmodel(y, weights = NULL, nullcats = c("keep", "downcode", "ignore"),
  start = NULL, reltol = 1e-10, deriv = c("sum", "diff"),
  hessian = TRUE, maxit = 100L, full = TRUE, ...)

Arguments

y	item response object that can be coerced (via as.matrix) to a numeric matrix with scores 0, 1, ... Typically, either already a matrix, data frame, or dedicated object of class itemresp.
weights	an optional vector of weights (interpreted as case weights).
deriv	character. If "sum" (the default), the first derivatives of the elementary symmetric functions are calculated with the sum algorithm. Otherwise ("diff") the difference algorithm (faster but numerically unstable) is used.
nullcats	character string, specifying how items with null categories (i.e., categories not observed) should be treated (see details below).
start	an optional vector of starting values.
hessian	logical. Should the Hessian of the final model be computed? If set to FALSE, the vcov method can only return NAs and consequently no standard errors or tests are available in the summary.
reltol, maxit, ...	further arguments passed to optim.
full	logical. Should a full model object be returned? If set to FALSE, no variance-covariance matrix and no matrix of estimating functions are computed.

Details

pcmodel provides a basic fitting function for partial credit models, intended as a building block for fitting partial credit trees. It estimates the partial credit model suggested by Masters (1982) under the cumulative threshold parameterization, i.e., the item-category parameters $\eta_{jk} = \sum_{\ell = 1}^{k}\delta_{jk}$ are estimated by the the function pcmodel.

Null categories, i.e., categories which have not been used, can be problematic when estimating a partial credit model. Several strategies have been suggested to cope with null categories. pcmodel allows to select from three possible strategies via the argument nullcats. If nullcats is set to "keep" (the default), the strategy suggested by Wilson & Masters (1993) is used to handle null categories. That basically means that the integrity of the response framework is maintained, i.e., no category scores are changed. This is not the case, when nullcats is set to "downcode". Then all categories above a null category are shifted down to close the existing gap. In both cases ("keep" and "downcode") the number of estimated parameters is reduced by the number of null categories. When nullcats is set to "ignore", these are literally ignored and a threshold parameter is estimated during the optimization nevertheless. This strategy is used by the related package eRm when fitting partial credit models via eRm::PCM.

pcmodel returns an object of class "pcmodel" for which several basic methods are available, including print, plot, summary, coef, vcov, logLik, discrpar, itempar, estfun, threshpar, and personpar.

Value

pcmodel returns an S3 object of class "pcmodel", i.e., a list the following components:

coefficients	a named vector of estimated item-category parameters (without the first item-category parameter which is constrained to 0),
vcov	covariance matrix of the parameters in the model,
data	modified data, used for model-fitting, i.e., cleaned for items without variance, centralized so that the first category is zero for all items, treated null categories as specified via argument "nullcats" and without observations with zero weight. Be careful, this is different than for objects of class "raschmodel" or "btmodel", where data contains the original data,
items	logical vector of length ncol(dat), indicating which items have variance (TRUE), i.e., are identified and have been used for the estimation or not (FALSE),
categories	list of length ncol(y), containing integer vectors starting from one to the number of categories minus one per item,
n	number of observations (with non-zero weights),
n_org	original number of observations in y,
weights	the weights used (if any),
na	logical indicating whether the data contain NAs,
nullcats	either NULL or, if there have been null categories, a list of length ncol(y) with logical vectors specifying which categories are null categories (TRUE) or not (FALSE),
esf	list of elementary symmetric functions and their derivatives for estimated parameters,
loglik	log-likelihood of the fitted model,
df	number of estimated parameters,
code	convergence code from optim,
iterations	number of iterations used by optim,
reltol	tolerance passed to optim,
call	original function call.

References

Masters GN (1992). A Rasch Model for Partial Credit Scoring. Psychometrika, 47(2), 149–174.

Wilson M, Masters GN (1993). The Partial Credit Model and Null Categories. Psychometrika, 58(1), 87–99.

See Also

gpcmodel, rsmodel, raschmodel, nplmodel, btmodel

Examples

o <- options(digits = 4)

## Verbal aggression data
data("VerbalAggression", package = "psychotools")

## Partial credit model for the other-to-blame situations
pcm <- pcmodel(VerbalAggression$resp[, 1:12])
summary(pcm)

## visualizations
plot(pcm, type = "profile")
plot(pcm, type = "regions")
plot(pcm, type = "piplot")
plot(pcm, type = "curves")
plot(pcm, type = "information")

## Get data of situation 1 ('A bus fails to
## stop for me') and induce a null category in item 2.
pcd <- VerbalAggression$resp[, 1:6, drop = FALSE]
pcd[pcd[, 2] == 1, 2] <- NA

## fit pcm to these data, comparing downcoding and keeping strategy
pcm_va_keep  <- pcmodel(pcd, nullcats = "keep")
pcm_va_down  <- pcmodel(pcd, nullcats = "downcode")

plot(x = coef(pcm_va_keep), y = coef(pcm_va_down),
     xlab = "Threshold Parameters (Keeping)",
     ylab = "Threshold Parameters (Downcoding)",
     main = "Comparison of two null category strategies (I2 with null category)", 
     pch = rep(as.character(1:6), each = 2)[-3])
abline(b = 1, a = 0)

options(digits = o$digits)

---

Extract Person Parameters of Item Response Models

Description

A class and generic function for representing and estimating the person parameters of a given item response model.

Usage

personpar(object, ...)
## S3 method for class 'raschmodel'
personpar(object, personwise = FALSE, ref = NULL,
  vcov = TRUE, interval = NULL, tol = 1e-8, ...)
## S3 method for class 'rsmodel'
personpar(object, personwise = FALSE, ref = NULL,
  vcov = TRUE, interval = NULL, tol = 1e-8, ...)
## S3 method for class 'pcmodel'
personpar(object, personwise = FALSE, ref = NULL,
  vcov = TRUE, interval = NULL, tol = 1e-8, ...)
## S3 method for class 'nplmodel'
personpar(object, personwise = FALSE, vcov = TRUE,
  interval = NULL, tol = 1e-6, method = "EAP", ...)
## S3 method for class 'gpcmodel'
personpar(object, personwise = FALSE, vcov = TRUE,
  interval = NULL, tol = 1e-6, method = "EAP", ...)

Arguments

object	a fitted model object for which person parameters should be returned/estimated.
personwise	logical. Should the distributional parameters of the latent person ability distribution be computed (default) or the person-wise (individual) person parameters? See below for details.
ref	a vector of labels or position indices of item parameters or a contrast matrix which should be used as restriction/for normalization. This argument will be passed over to internal calls of itempar.
vcov	logical. Should a covariance matrix be returned/estimated for the person parameter estimates? See also details below.
interval	numeric vector of length two, specifying an interval for uniroot or fscores to calculate the person parameter estimates.
tol	numeric tolerance passed to uniroot or fscores.
method	type of estimation method being passed to fscores.
...	further arguments which are passed to optim in case of vcov = TRUE or fscores.

Details

personpar is both a class to represent person parameters of item response models as well as a generic function. The generic function can be used to return/estimate the person parameters of a given item response model.

By default, the function personpar() reports the distribution parameters of the assumed person ability distribution. For models estimated by marginal maximum likelihood estimation (MML) this is the mean/variance of the underlying normal distribution, whereas for models estimated by conditional maximum likelihood estimation (CML) this is a discrete distribution with one estimation for each observed raw score in the data.

Alternatively, when setting personwise = TRUE, the person parameter for each person/subject in the underlying data set can be extracted. In the CML case, this simply computes the raw score for each person and then extracts the corresponding person parameter. In the MML case, this necessitates (numerically) integrating out the individual person parameters (also known as factor scores or latent trait estimates) based on the underlying normal distribution.

More specifically, the following algorithms are employed for obtaining the distributional person parameters:

• In the MML case – i.e., for nplmodels and gpcmodels – the distributional parameters are already part of the model specification. In a single-group specification and in the reference group of a multi-group specification the mean/variance parameters are fixed to 0/1. In the multi-group case the remaining mean/variance parameters were already estimated along with all other model parameters and simply need to be extracted. Analogously, the corresponding variance-covariance matrix just needs to be extracted and has zero covariances in the cells corresponding to fixed parameters.

• In the CML case – i.e., raschmodels, rsmodels, and pcmodels – the distributional parameters are estimated via uniroot() with the estimation equations given by Hoijtink & Boomsma (1995) as well as Andersen (1995). This approach is fast and estimates for all possible raw scores are available. If the covariance matrix of the estimated person parameters is requested (vcov = TRUE), an additional call of optim is employed to obtain the Hessian numerically. With this approach, person parameters are available only for observed raw scores.

As already explained above, obtaining the person-wise (individual) person paremeters (or ability estimates or factor scores) is straightforward in the CML case. In the MML case, fscores is used, see Chalmers (2012) for further details. If personwise = TRUE, the associated variance-covariance matrix is not provided and simply a matrix with NAs is returned. (The same is done for vcov = FALSE.)

For objects of class personpar, several methods to standard generic functions exist: print, coef, vcov. coef and vcov can be used to extract the person parameters and covariance matrix without additional attributes. Based on this Wald tests or confidence intervals can be easily computed, e.g., via confint.

Value

A named vector with person parmeters of class personpar and additional attributes "model" (the model name), "vcov" (the covariance matrix of the estimates if vcov = TRUE or an NA-matrix otherwise) and "type" (the type of the parameters, depending on personwise).

References

Andersen EB (1995). Polytomous Rasch Models and Their Estimation. In Fischer GH, Molenaar IW (eds.). Rasch Models: Foundations, Recent Developments, and Applications. Springer, New York.

Chalmers RP (2012). mirt: A Multidimensional Item Response Theory Package for the R Environment. Journal of Statistical Software, 48(6), 1–29.

Hoijtink H, Boomsma A (1995). On Person Parameter Estimation in the Dichotomous Rasch Model. In Fischer GH, Molenaar IW (eds.). Rasch Models: Foundations, Recent Developments, and Applications. Springer, New York.

See Also

itempar, threshpar, discrpar, guesspar, upperpar

Examples

o <- options(digits = 3)

## load verbal aggression data
data("VerbalAggression", package = "psychotools")

## fit a Rasch model to dichotomized verbal aggression data and
ram <- raschmodel(VerbalAggression$resp2)

## extract person parameters
## (= parameters of the underlying ability distribution)
rap <- personpar(ram)
rap

## extract variance-covariance matrix and standard errors
vc <- vcov(rap)
sqrt(diag(vc))

## Wald confidence intervals
confint(rap)

## now match each person to person parameter with the corresponding raw score
personpar(ram, personwise = TRUE)[1:6]

## person parameters for RSM/PCM fitted to original polytomous data
rsm <- rsmodel(VerbalAggression$resp)
pcm <- pcmodel(VerbalAggression$resp)
cbind(personpar(rsm, vcov = FALSE), personpar(pcm, vcov = FALSE))

if(requireNamespace("mirt")) {
## fit a 2PL accounting for gender impact and
twoplm <- nplmodel(VerbalAggression$resp2, impact = VerbalAggression$gender)

## extract the person parameters
## (= mean/variance parameters from the normal ability distribution)
twoplp <- personpar(twoplm)
twoplp

## extract the standard errors
sqrt(diag(vcov(twoplp)))

## Wald confidence intervals
confint(twoplp)

## now look at the individual person parameters
## (integrated out over the normal ability distribution)
personpar(twoplm, personwise = TRUE)[1:6]
}

options(digits = o$digits)

---