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

View source: R/immer_hrm.R

immer_hrmR Documentation

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)

immer documentation built on May 29, 2024, 11:52 a.m.