Estimates the hierarchical rater model (HRM; Patz et al., 2002; see Details) with Markov Chain Monte Carlo using a Metropolis-Hastings algorithm.
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 immer_HRM
summary(object,digits=3,...)
# S3 method for immer_HRM
plot(x,...)
# S3 method for immer_HRM
logLik(object,...)
# S3 method for immer_HRM
anova(object,...)
# S3 method for immer_HRM
IRT.likelihood(object,...)
# S3 method for immer_HRM
IRT.posterior(object,...)
Data frame with item responses
Person identifiers
Rater identifiers
Number of iterations
Number of burnin iterations
Maximum number of samples to be saved.
Parameters for prior distributions
Logical indicating whether
Logical indicating whether
Optional logical indicating whether the mean
Type of est.phi="a"
, then est.phi="r"
, then est.phi="i"
it is item specific
(est.phi="n"
, no
Type of est.phi
, but also
allows est.psi="e"
(exchangeable) which means
Parameters for Metropolis Hastings sampling. The standard deviation of the proposal distribution is adaptively computed (Browne & Draper, 2000).
Grid of
Object of class immer_HRM
Number of digits after decimal to print
Object of class immer_HRM
Further arguments to be passed. See
sirt::plot.mcmc.sirt
for options in plot
.
A list with following entries
Data frame with estimated person parameters
Data frame with estimated item parameters
Data frame with estimated rater parameters
Estimated item and trait distribution parameters arranged in vectors and matrices.
Estimated standard deviation
Estimated mean
Object of class mcmc.list
for coda package.
Summary of all parameters
EAP reliability
Parameters for information criteria
Individual likelihood evaluated at theta_like
Individual posterior evaluated at theta_like
Grid of
Discretized
Log-likelihood value
Updated parameters in Metropolis-Hastings sampling
The hierarchical rater model is defined at the level of persons
At the level of ratings, the model is defined as
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.
Simulate the HRM with simulate_HRM
.
# NOT RUN {
library(sirt)
library(TAM)
#############################################################################
# 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 <- 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::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::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 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::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)
# }
Run the code above in your browser using DataLab