Extended Rasch Model with R

Part 1: Do we need the Rasch Model?

Rasch Model

to measure latent traits (e.g., ability or attitude)
origin in psychology
now applied in various disciplines:
health, education, social sciences, economics, finance, . . .

distinction: latent variables - manifest observations

latent traits are usually assessed through the responses of a
sample of subjects to a set of items (questions)

Example

Aim: obtain a measurement for Stress

usual procedure:

Sum up number of Yes-answers
obtain scores for persons (high value - highly stressed)
do some statistics (eg. means for female/ males)

Is this justifiable?

Questions

How about the sum scores?
same stress for persons with same score?
YES to items 1 and 2 the same as YES to items 3 and 4?
do all items measure the same latent trait? (unidimensional)
do we mix up several latent traits?
What about the scale properties?
is the sumscore interval scaled? (sum of binary variables)
the central limit theorem justifiable? (the i.i.d. assumption)
How about linearity?
consider training: progress 1 –> 3 the same as 2 –> 4
Distinction manifest observations from latent trait?
manifest sum score is equated with location on latent trait

need for some thoughts on measurement!

Measurement

Measurement is the mapping from a set of predefined objects and their empirically observable relations onto a set of numbers and their relations

Specific Objectivity (Rasch, 1960)

Part 2 : Praticum

The R package eRm (extended Rasch modelling)

library(eRm)

Call the data and head the data

load("/home/yg1hw/eRm file/stress.RData")
head(stress)

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    0    1    1    0    1    0
## [2,]    0    0    1    1    1    0
## [3,]    0    0    0    0    0    1
## [4,]    0    0    1    0    0    0
## [5,]    1    0    1    1    0    1
## [6,]    1    0    1    0    1    1

Fitting the RM

rm.res <- RM(stress)
rm.res

## 
## Results of RM estimation: 
## 
## Call:  RM(X = stress) 
## 
## Conditional log-likelihood: -202.1232 
## Number of iterations: 13 
## Number of parameters: 5 
## 
## Item (Category) Difficulty Parameters (eta):
##                  I2          I3       I4        I5         I6
## Estimate -0.1101525 -0.06109055 0.241353 0.6812941 -0.3995985
## Std.Err   0.2027044  0.20317040 0.207776 0.2203618  0.2014671

default is: RM(datamatrix, sum0 = TRUE, other options)
sum0 defines constraints (for estimability): TRUE . . . sum zero, FALSE . . . first item set to 0
the output gives dificulty parameters

Constraints and the Design Matrix

model.matrix(rm.res)

##         eta 1 eta 2 eta 3 eta 4 eta 5
## beta I1    -1    -1    -1    -1    -1
## beta I2     1     0     0     0     0
## beta I3     0     1     0     0     0
## beta I4     0     0     1     0     0
## beta I5     0     0     0     1     0
## beta I6     0     0     0     0     1

model.matrix(RM(stress, sum0 = FALSE))

##         eta 1 eta 2 eta 3 eta 4 eta 5
## beta I1     0     0     0     0     0
## beta I2     1     0     0     0     0
## beta I3     0     1     0     0     0
## beta I4     0     0     1     0     0
## beta I5     0     0     0     1     0
## beta I6     0     0     0     0     1

summary(rm.res)

## 
## Results of RM estimation: 
## 
## Call:  RM(X = stress) 
## 
## Conditional log-likelihood: -202.1232 
## Number of iterations: 13 
## Number of parameters: 5 
## 
## Item (Category) Difficulty Parameters (eta): with 0.95 CI:
##    Estimate Std. Error lower CI upper CI
## I2   -0.110      0.203   -0.507    0.287
## I3   -0.061      0.203   -0.459    0.337
## I4    0.241      0.208   -0.166    0.649
## I5    0.681      0.220    0.249    1.113
## I6   -0.400      0.201   -0.794   -0.005
## 
## Item Easiness Parameters (beta) with 0.95 CI:
##         Estimate Std. Error lower CI upper CI
## beta I1    0.352      0.201   -0.043    0.747
## beta I2    0.110      0.203   -0.287    0.507
## beta I3    0.061      0.203   -0.337    0.459
## beta I4   -0.241      0.208   -0.649    0.166
## beta I5   -0.681      0.220   -1.113   -0.249
## beta I6    0.400      0.201    0.005    0.794

Extracting Information the item parameter estimates

coef(rm.res)

##     beta I1     beta I2     beta I3     beta I4     beta I5     beta I6 
##  0.35180555  0.11015250  0.06109055 -0.24135303 -0.68129407  0.39959849

the variance-covariance matrix of item parameter estimates

vcov(rm.res)

##              [,1]         [,2]         [,3]         [,4]         [,5]
## [1,]  0.041089085 -0.007883439 -0.008375267 -0.009639680 -0.007580398
## [2,] -0.007883439  0.041278211 -0.008399888 -0.009630899 -0.007669699
## [3,] -0.008375267 -0.008399888  0.043170866 -0.009727776 -0.008336458
## [4,] -0.009639680 -0.009630899 -0.009727776  0.048559311 -0.009799404
## [5,] -0.007580398 -0.007669699 -0.008336458 -0.009799404  0.040589011

confindence intervals for the item parameter estimates

confint(rm.res, "beta")

##                2.5 %     97.5 %
## beta I1 -0.043119498  0.7467306
## beta I2 -0.287140876  0.5074459
## beta I3 -0.337116114  0.4592972
## beta I4 -0.648586506  0.1658804
## beta I5 -1.113195206 -0.2493929
## beta I6  0.004730141  0.7944668

the conditional log likelihood

logLik(rm.res)

## Conditional log Lik.: -202.1232 (df=5)

Plot ICCs

plotjointICC(rm.res, legend=TRUE, cex=0.75, xlim = c(-5, 5))

Plot single ICC (Example-Item 3)

plotICC(rm.res, item.subset = 3, ask = FALSE)

Plot ICCs

plotICC(rm.res, item.subset = 1:6, ask = F, empICC = list("raw"),empCI = list(lty = "solid"))

Plot Person Item Map

plotPImap(rm.res,sorted=TRUE)

Person Parameter Estimation

pp <- person.parameter(rm.res)
print(pp)

## 
## Person Parameters:
## 
##  Raw Score     Estimate Std.Error
##          0 -2.644767051        NA
##          1 -1.653469544 1.1039389
##          2 -0.717574135 0.8777577
##          3 -0.002996393 0.8301053
##          4  0.713595524 0.8800763
##          5  1.655491598 1.1076270
##          6  2.653641029        NA

If NAs in the data, different person parameters are estimated for every NA-pattern group

Methods for Person Parameter Estimation Results

logLik(pp)

## Unconditional (joint) log Lik.: -17.37912 (df=5)

confint(pp)

## $NAgroup1
##          2.5 %    97.5 %
## P1  -1.6299728 1.6239800
## P2  -1.6299728 1.6239800
## P3  -3.8171500 0.5102109
## P4  -3.8171500 0.5102109
## P5  -1.0113223 2.4385134
## P6  -1.0113223 2.4385134
## P7  -2.4379477 1.0027994
## P8  -3.8171500 0.5102109
## P9  -1.0113223 2.4385134
## P10 -3.8171500 0.5102109
## P11 -2.4379477 1.0027994
## P13 -2.4379477 1.0027994
## P14 -2.4379477 1.0027994
## P15 -3.8171500 0.5102109
## P16 -1.6299728 1.6239800
## P18 -2.4379477 1.0027994
## P19 -0.5154175 3.8264007
## P20 -3.8171500 0.5102109
## P21 -0.5154175 3.8264007
## P22 -3.8171500 0.5102109
## P23 -1.6299728 1.6239800
## P24 -3.8171500 0.5102109
## P25 -0.5154175 3.8264007
## P26 -1.6299728 1.6239800
## P28 -1.6299728 1.6239800
## P30 -1.0113223 2.4385134
## P31 -2.4379477 1.0027994
## P32 -3.8171500 0.5102109
## P33 -1.6299728 1.6239800
## P34 -1.6299728 1.6239800
## P35 -1.6299728 1.6239800
## P36 -2.4379477 1.0027994
## P37 -3.8171500 0.5102109
## P38 -1.0113223 2.4385134
## P40 -2.4379477 1.0027994
## P42 -3.8171500 0.5102109
## P43 -3.8171500 0.5102109
## P44 -1.6299728 1.6239800
## P45 -2.4379477 1.0027994
## P46 -1.0113223 2.4385134
## P47 -1.0113223 2.4385134
## P48 -2.4379477 1.0027994
## P49 -2.4379477 1.0027994
## P50 -1.0113223 2.4385134
## P51 -2.4379477 1.0027994
## P52 -1.6299728 1.6239800
## P53 -0.5154175 3.8264007
## P54 -0.5154175 3.8264007
## P55 -1.0113223 2.4385134
## P56 -1.0113223 2.4385134
## P58 -1.6299728 1.6239800
## P59 -3.8171500 0.5102109
## P60 -1.6299728 1.6239800
## P61 -3.8171500 0.5102109
## P62 -3.8171500 0.5102109
## P63 -3.8171500 0.5102109
## P64 -1.0113223 2.4385134
## P65 -1.6299728 1.6239800
## P66 -2.4379477 1.0027994
## P67 -3.8171500 0.5102109
## P68 -2.4379477 1.0027994
## P69 -1.6299728 1.6239800
## P70 -2.4379477 1.0027994
## P71 -3.8171500 0.5102109
## P72 -1.0113223 2.4385134
## P73 -1.6299728 1.6239800
## P74 -1.0113223 2.4385134
## P75 -0.5154175 3.8264007
## P76 -3.8171500 0.5102109
## P77 -1.0113223 2.4385134
## P78 -0.5154175 3.8264007
## P79 -3.8171500 0.5102109
## P80 -1.6299728 1.6239800
## P82 -2.4379477 1.0027994
## P83 -3.8171500 0.5102109
## P84 -2.4379477 1.0027994
## P85 -3.8171500 0.5102109
## P88 -3.8171500 0.5102109
## P89 -2.4379477 1.0027994
## P91 -1.0113223 2.4385134
## P92 -3.8171500 0.5102109
## P93 -1.6299728 1.6239800
## P95 -1.6299728 1.6239800
## P96 -1.0113223 2.4385134
## P97 -2.4379477 1.0027994
## P98 -1.0113223 2.4385134

attention: confint(pp) gives values for all subjects if there are NAs in the data, confidence intervals are printed for each NA group

Plot of Person Parameter Estimates

plot(pp)

Part 3: Testing the Rasch Model

Using eRm : LRtest()

lr <- LRtest(rm.res)
print(lr)

## 
## Andersen LR-test: 
## LR-value: 6.448 
## Chi-square df: 5 
## p-value:  0.265

factors as split criteria

sex <- rep(c("male", "female"), each = 50)
LRtest(rm.res, splitcr = sex)

## 
## Andersen LR-test: 
## LR-value: 5.635 
## Chi-square df: 5 
## p-value:  0.343

more than two groups (e.g., based on rawscores)

gr3 <- cut(rowSums(stress), 3)
LRtest(rm.res, splitcr = gr3)

## Warning in LRtest.Rm(rm.res, splitcr = gr3): 
## The following items were excluded due to inappropriate response patterns
## within subgroups:
## I1 I3
## 
## Full and subgroup models are estimated without these items!

## 
## Andersen LR-test: 
## LR-value: 9.055 
## Chi-square df: 6 
## p-value:  0.171

using eRm: Waldtest()

wt <- Waldtest(rm.res)
wt

## 
## Wald test on item level (z-values):
## 
##         z-statistic p-value
## beta I1      -0.045   0.964
## beta I2       0.189   0.850
## beta I3      -0.564   0.573
## beta I4      -0.892   0.372
## beta I5      -0.660   0.509
## beta I6       2.463   0.014

Graphical Procedure

plotGOF(lr)

lr <- LRtest(rm.res, se = T)
plotGOF(lr, conf = list(), xlim = c(-1.5, 1.5), ylim = c(-1.5, 1.5))

plotGOF(lr, ctrline = list(), xlim = c(-1.5, 1.5), ylim = c(-1.5, 1.5))

both require objects obtained from person.parameter()

pp <- person.parameter(rm.res)
itemfit(pp)

## 
## Itemfit Statistics: 
##      Chisq df p-value Outfit MSQ Infit MSQ Outfit t Infit t
## I1  78.028 85   0.691      0.907     0.948    -0.74   -0.51
## I2  86.465 85   0.435      1.005     0.992     0.08   -0.04
## I3  77.243 85   0.713      0.898     0.918    -0.78   -0.82
## I4  87.174 85   0.414      1.014     1.005     0.14    0.08
## I5  73.947 85   0.798      0.860     0.863    -0.72   -1.16
## I6 106.633 85   0.056      1.240     1.190     1.85    1.87

Graphical Procedure:plotPWmap()

plotPWmap(rm.res)

diferent colours, plotting symbols etc. possible

plotPWmap(rm.res, pmap = TRUE, imap=TRUE,
          item.subset = "all", person.subset = 1:5,
          mainitem = "Item Map", mainperson = "Person Map",
          mainboth="Item/Person Map",
          latdim = "Latent Dimension",
          tlab = "Infit t statistic",
          pp = NULL, cex.gen = 0.6, cex.pch=1,
          person.pch = 16, item.pch = 16,
          personCI = NULL, itemCI = list(), horiz=FALSE)

pers <- person.parameter(rm.res)
perstheta<-data.frame(pers$thetapar)
library(data.table)
setnames(perstheta, "NAgroup1", "theta")#Rename theta
itembeta<-data.frame(pers$betapar)
setnames(itembeta, "pers.betapar", "beta")#Rename beta
WLEestimates.mod1 <- perstheta
thresholds.mod1 <- itembeta
library(WrightMap)
library(RColorBrewer)
wrightMap(WLEestimates.mod1, thresholds.mod1, main.title = "This is  Wright Map"
          , axis.persons = "This is the person distribution"
          , axis.items = "This are my survey questions")

##                beta
## beta I1  0.35180555
## beta I2  0.11015250
## beta I3  0.06109055
## beta I4 -0.24135303
## beta I5 -0.68129407
## beta I6  0.39959849

Martin-Löf Test

MLoef(rm.res, splitcr = "median")

## 
## Martin-Loef-Test (split criterion: median)
## LR-value: 8.009 
## Chi-square df: 8 
## p-value: 0.433

MLoef(rm.res, splitcr = c(1, 1, 1, 1, 0, 0))

## 
## Martin-Loef-Test (split criterion: user-defined)
## LR-value: 7.934 
## Chi-square df: 7 
## p-value: 0.338