immer_HRM: Hierarchical Rater Model (Patz et al., 2002)

In alexanderrobitzsch/immer: Item Response Models for Multiple Ratings

Hierarchical Rater Model (Patz et al., 2002)

Description

Estimates the hierarchical rater model (HRM; Patz et al., 2002; see Details) with Markov Chain Monte Carlo using a Metropolis-Hastings algorithm.

Usage

immer_hrm(dat, pid, rater, iter, burnin, N.save=3000, prior=NULL,  est.a=FALSE,
           est.sigma=TRUE,  est.mu=FALSE, est.phi="a", est.psi="a",
           MHprop=NULL, theta_like=seq(-10,10,len=30), sigma_init=1, print_iter=20 )

## S3 method for class 'immer_hrm'
summary(object, digits=3, file=NULL, ...)

## S3 method for class 'immer_hrm'
plot(x,...)

## S3 method for class 'immer_hrm'
logLik(object,...)

## S3 method for class 'immer_hrm'
anova(object,...)

## S3 method for class 'immer_hrm'
IRT.likelihood(object,...)

## S3 method for class 'immer_hrm'
IRT.posterior(object,...)

Arguments

dat	Data frame with item responses
pid	Person identifiers
rater	Rater identifiers
iter	Number of iterations
burnin	Number of burnin iterations
N.save	Maximum number of samples to be saved.
prior	Parameters for prior distributions
est.a	Logical indicating whether a parameter should be estimated.
est.sigma	Logical indicating whether \sigma parameter should be estimated.
est.mu	Optional logical indicating whether the mean \mu of the trait \theta should be estimated.
est.phi	Type of \phi _{ir} parameters to be estimated. If est.phi="a", then \phi_{ir} is estimated for all items and all raters. If est.phi="r", then \phi_{ir}=\phi_r is rater specific, while for est.phi="i" it is item specific (\phi_{ir}=\phi_i). In case of est.phi="n", no \phi parameters are estimated and all \phi parameters are fixed at 0.
est.psi	Type of \psi_{ir} parameters to be estimated. The arguments follow the same conventions as est.phi, but also allows est.psi="e" (exchangeable) which means \psi_{ir}=\psi, i.e assuming the same \psi parameter for all items and raters.
MHprop	Parameters for Metropolis Hastings sampling. The standard deviation of the proposal distribution is adaptively computed (Browne & Draper, 2000).
theta_like	Grid of \theta values to be used for likelihood approximation
sigma_init	Initial value for sigma
print_iter	Integer indicating that after each print_iterth iteration output on the console should be displayed.
object	Object of class immer_hrm
digits	Number of digits after decimal to print
file	Name of a file in which the output should be sunk
x	Object of class immer_hrm
...	Further arguments to be passed. See sirt::plot.mcmc.sirt for options in plot.

Details

The hierarchical rater model is defined at the level of persons

P( \xi _{pi}=\xi | \theta_p ) \propto \exp ( \xi \cdot a_i \cdot \theta_p - b_{ik} )

where \theta_p is normally distributed with mean \mu and standard deviation \sigma.

At the level of ratings, the model is defined as

P( X_{pir}=x | \theta_p, \xi_{pi} ) \propto \exp( - ( x - \xi_{pi} - \phi_{ir} )^2 / ( 2 \cdot \psi_{ir} ) )

Value

A list with following entries

person	Data frame with estimated person parameters
item	Data frame with estimated item parameters
rater_pars	Data frame with estimated rater parameters
est_pars	Estimated item and trait distribution parameters arranged in vectors and matrices.
sigma	Estimated standard deviation \sigma of trait \theta
mu	Estimated mean \mu of trait \theta
mcmcobj	Object of class mcmc.list for coda package.
summary.mcmcobj	Summary of all parameters
EAP.rel	EAP reliability
ic	Parameters for information criteria
f.yi.qk	Individual likelihood evaluated at theta_like
f.qk.yi	Individual posterior evaluated at theta_like
theta_like	Grid of \theta values for likelihood approximation
pi.k	Discretized \theta distribution
like	Log-likelihood value
MHprop	Updated parameters in Metropolis-Hastings sampling

References

Browne, W. J., & Draper, D. (2000). Implementation and performance issues in the Bayesian and likelihood fitting of multilevel models. Computational Statistics, 15, 391-420.

Patz, R. J., Junker, B. W., Johnson, M. S., & Mariano, L. T. (2002). The hierarchical rater model for rated test items and its application to large-scale educational assessment data. Journal of Educational and Behavioral Statistics, 27(4), 341-384.

See Also

Simulate the HRM with immer_hrm_simulate.

Examples

## Not run: 
library(sirt)
library(TAM)

#############################################################################
# EXAMPLE 1: Simulated data using the immer::immer_hrm_simulate() function
#############################################################################

# define data generating parameters
set.seed(1997)
N <- 500  # number of persons
I <- 4    # number of items
R <- 3    # number of raters
K <- 3    # maximum score
sigma <- 2  # standard deviation
theta <- stats::rnorm( N, sd=sigma )  # abilities
# item intercepts
b <- outer( seq( -1.5, 1.5, len=I), seq( -2, 2, len=K), "+" )
# item loadings
a <- rep(1,I)
# rater severity parameters
phi <- matrix( c(-.3, -.2, .5), nrow=I, ncol=R, byrow=TRUE )
phi <- phi + stats::rnorm( phi, sd=.3 )
phi <- phi - rowMeans(phi)
# rater variability parameters
psi <- matrix( c(.1, .4, .8), nrow=I, ncol=R, byrow=TRUE )
# simulate HRM data
data <- immer::immer_hrm_simulate( theta, a, b, phi=phi, psi=psi )
pid <- data$pid
rater <- data$rater
dat <- data[, - c(1:2) ]

#----------------------------------------------------------------
#*** Model 1: estimate HRM with equal item slopes
iter <- 3000
burnin <- 500
mod1 <- immer::immer_hrm( dat, pid, rater, iter=iter,  burnin=burnin )
summary(mod1)
plot(mod1,layout=2,ask=TRUE)

# relations among convergence diagnostic statistics
par(mfrow=c(1,2))
plot( mod1$summary.mcmcobj$PercVarRatio, log(mod1$summary.mcmcobj$effSize), pch=16)
plot( mod1$summary.mcmcobj$PercVarRatio, mod1$summary.mcmcobj$Rhat, pch=16)
par(mfrow=c(1,1))

# extract individual likelihood
lmod1 <- IRT.likelihood(mod1)
str(lmod1)
# extract log-likelihood value
logLik(mod1)

# write coda files into working directory
sirt::mcmclist2coda(mod1$mcmcobj, name="hrm_mod1")

#----------------------------------------------------------------
#*** Model 2: estimate HRM with estimated item slopes
mod2 <- immer::immer_hrm( dat, pid, rater, iter=iter,  burnin=burnin,
            est.a=TRUE, est.sigma=FALSE)
summary(mod2)
plot(mod2,layout=2,ask=TRUE)

# model comparison
anova( mod1, mod2 )

#----------------------------------------------------------------
#*** Model 3: Some prior specifications
prior <- list()
# prior on mu
prior$mu$M <- .7
prior$mu$SD <- 5
# fixed item parameters for first item
prior$b$M <-  matrix( 0, nrow=4, ncol=3 )
prior$b$M[1,] <- c(-2,0,2)
prior$b$SD <-  matrix( 10, nrow=4, ncol=3 )
prior$b$SD[1,] <- 1E-4
# estimate model
mod3 <- immer::immer_hrm( dat, pid, rater, iter=iter,  burnin=burnin, prior=prior)
summary(mod3)
plot(mod3)

#----------------------------------------------------------------
#*** Model 4: Multi-faceted Rasch models in TAM package

# create facets object
facets <- data.frame( "rater"=rater )

#-- Model 4a: only main rater effects
form <- ~ item*step + rater
mod4a <- TAM::tam.mml.mfr( dat, pid=pid, facets=facets, formulaA=form)
summary(mod4a)

#-- Model 4b: item specific rater effects
form <- ~ item*step + item*rater
mod4b <- TAM::tam.mml.mfr( dat, pid=pid, facets=facets, formulaA=form)
summary(mod4b)

#----------------------------------------------------------------
#*** Model 5: Faceted Rasch models with sirt::rm.facets

#--- Model 5a: Faceted Rasch model with only main rater effects
mod5a <- sirt::rm.facets( dat, pid=pid, rater=rater )
summary(mod5a)

#--- Model 5b: allow rater slopes for different rater discriminations
mod5b <- sirt::rm.facets( dat, pid=pid, rater=rater, est.a.rater=TRUE )
summary(mod5b)

#############################################################################
# EXAMPLE 2: data.ratings1 (sirt package)
#############################################################################

data(data.ratings1, package="sirt")
dat <- data.ratings1

# set number of iterations and burnin iterations
set.seed(87)
iter <- 1000
burnin <- 500
# estimate model
mod <- immer::immer_hrm( dat[,  paste0("k",1:5) ], pid=dat$idstud, rater=dat$rater,
             iter=iter, burnin=burnin )
summary(mod)
plot(mod, layout=1, ask=TRUE)
plot(mod, layout=2, ask=TRUE)

## End(Not run)

alexanderrobitzsch/immer documentation built on March 27, 2024, 5:48 a.m.