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

Description Usage Arguments Details Value References See Also Examples

Description

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

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
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 σ parameter should be estimated.

est.mu

Optional logical indicating whether the mean μ of the trait θ should be estimated.

est.phi

Type of φ _{ir} parameters to be estimated. If est.phi="a", then φ_{ir} is estimated for all items and all raters. If est.phi="r", then φ_{ir}=φ_r is rater specific, while for est.phi="i" it is item specific (φ_{ir}=φ_i). In case of est.phi="n", no φ parameters are estimated and all φ parameters are fixed at 0.

est.psi

Type of ψ_{ir} parameters to be estimated. The arguments follow the same conventions as est.phi, but also allows est.psi="e" (exchangeable) which means ψ_{ir}=ψ, i.e assuming the same ψ 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 θ 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( ξ _{pi}=ξ | θ_p ) \propto \exp ( ξ \cdot a_i \cdot θ_p - b_{ik} )

where θ_p is normally distributed with mean μ and standard deviation σ.

At the level of ratings, the model is defined as

P( X_{pir}=x | θ_p, ξ_{pi} ) \propto \exp( - ( x - ξ_{pi} - φ_{ir} )^2 / ( 2 \cdot ψ_{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 σ of trait θ

mu

Estimated mean μ of trait θ

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 θ values for likelihood approximation

pi.k

Discretized θ 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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
## 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 Aug. 16, 2020, 10:11 p.m.