immer (version 0.4-0)

immer_HRM: 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) )

## S3 method for class 'immer_HRM':
summary(object,digits=3,...)

## 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 i

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

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

object

Object of class immer_HRM

digits

Number of digits after decimal to print

x

Object of class immer_HRM

...

Further arguments to be passed. See sirt::plot.mcmc.sirt for options in plot.

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

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} ) )$$

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 simulate_HRM.

Examples

library(sirt)
library(TAM)
	
#############################################################################
# SIMULATED EXAMPLE 1: Simulated data using the "simulate_HRM" 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 <- 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 + 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 <- simulate_HRM( 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_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_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_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 in sirt package

#--- 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
res <- immer_HRM( dat[ ,  paste0("k",1:5) ] , pid=dat$idstud , rater=dat$rater , 
             iter=iter ,  burnin=burnin  )
summary(res)
plot(res , layout=1 , ask=TRUE)
plot(res , layout=2 , ask=TRUE)