Rasch model example
Szu-Yu Chen
23 November 2020
1 Data Management
- 載入資料並改變項名稱
dta1123 <- read.table("C:/Users/ASUS/Desktop/data/starter.txt", h=T)
str(dta1123)## 'data.frame': 1350 obs. of 3 variables:
## $ sid : int 1 1 1 1 1 1 1 1 1 2 ...
## $ item: int 1 2 3 4 5 6 7 8 9 1 ...
## $ resp: int 1 1 1 1 1 1 1 1 0 1 ...names(dta1123) <- c("ID","Item","Response")
head(dta1123)## ID Item Response
## 1 1 1 1
## 2 1 2 1
## 3 1 3 1
## 4 1 4 1
## 5 1 5 1
## 6 1 6 12.資料長變寬
library(tidyverse)## -- Attaching packages ----------------------------------------------------------- tidyverse 1.3.0 --## √ ggplot2 3.3.2 √ purrr 0.3.4
## √ tibble 3.0.3 √ dplyr 1.0.2
## √ tidyr 1.1.2 √ stringr 1.4.0
## √ readr 1.3.1 √ forcats 0.5.0## -- Conflicts -------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()dtaw<- dta1123 %>%
spread(key=Item, value= Response)
head(dtaw)## ID 1 2 3 4 5 6 7 8 9
## 1 1 1 1 1 1 1 1 1 1 0
## 2 2 1 1 1 0 0 0 1 0 0
## 3 3 1 1 1 1 0 1 1 0 1
## 4 4 1 1 1 1 1 0 1 0 0
## 5 5 1 1 1 0 0 1 1 0 0
## 6 6 0 1 1 1 0 0 0 0 03.載入IRT package(ltm)
1.1 Basic Graph
每個人各試題的答對率
library(ltm)## Warning: package 'ltm' was built under R version 4.0.3## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: msm## Warning: package 'msm' was built under R version 4.0.3## Loading required package: polycor## Warning: package 'polycor' was built under R version 4.0.3plot(descript(dtaw[,2:9]),ylim=c(0,1), type='b', main="Item response curves")
grid()
2 Rasch Model
- 估計IRT參數 先將資料中ID拿掉，以item的答題反應做單參數IRT估計(Rasch)
dtawIRT<-subset(dtaw, select = -1)
m0<- rasch(dtawIRT, IRT.param = TRUE)
summary(m0)##
## Call:
## rasch(data = dtawIRT, IRT.param = TRUE)
##
## Model Summary:
## log.Lik AIC BIC
## -703.7896 1427.579 1457.685
##
## Coefficients:
## value std.err z.vals
## Dffclt.1 -2.6202 0.3778 -6.9354
## Dffclt.2 -1.5390 0.2531 -6.0810
## Dffclt.3 -1.2108 0.2273 -5.3280
## Dffclt.4 -1.1354 0.2221 -5.1118
## Dffclt.5 -0.9205 0.2091 -4.4015
## Dffclt.6 -0.8521 0.2055 -4.1453
## Dffclt.7 -1.6766 0.2655 -6.3138
## Dffclt.8 -0.1870 0.1851 -1.0098
## Dffclt.9 0.7763 0.2040 3.8051
## Dscrmn 1.1100 0.1238 8.9689
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Convergence: 0
## max(|grad|): 1.8e-05
## quasi-Newton: BFGS3 Graph
3.1 綜合ICC
All items at once
plot(m0)
3.2 個別ICC
每個item 的 ICC
1.方法1
library(eRm)## Warning: package 'eRm' was built under R version 4.0.3rm0 <- RM(dtawIRT)
plotICC(rm0, item.subset=1:9, ask=F, empICC=list("raw"), empCI=list(lty="solid"))
2.方法2
summary(m0)##
## Call:
## rasch(data = dtawIRT, IRT.param = TRUE)
##
## Model Summary:
## log.Lik AIC BIC
## -703.7896 1427.579 1457.685
##
## Coefficients:
## value std.err z.vals
## Dffclt.1 -2.6202 0.3778 -6.9354
## Dffclt.2 -1.5390 0.2531 -6.0810
## Dffclt.3 -1.2108 0.2273 -5.3280
## Dffclt.4 -1.1354 0.2221 -5.1118
## Dffclt.5 -0.9205 0.2091 -4.4015
## Dffclt.6 -0.8521 0.2055 -4.1453
## Dffclt.7 -1.6766 0.2655 -6.3138
## Dffclt.8 -0.1870 0.1851 -1.0098
## Dffclt.9 0.7763 0.2040 3.8051
## Dscrmn 1.1100 0.1238 8.9689
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Convergence: 0
## max(|grad|): 1.8e-05
## quasi-Newton: BFGScoef(m0)## Dffclt Dscrmn
## 1 -2.6201580 1.11
## 2 -1.5389723 1.11
## 3 -1.2108345 1.11
## 4 -1.1354132 1.11
## 5 -0.9205500 1.11
## 6 -0.8520519 1.11
## 7 -1.6765791 1.11
## 8 -0.1869645 1.11
## 9 0.7762629 1.11plot(m0, type="ICC", items = 1)
plot(m0, type="ICC", items = 2)
plot(m0, type="ICC", items = 3)
plot(m0, type="ICC", items = 4)
plot(m0, type="ICC", items = 5)
plot(m0, type="ICC", items = 6)
plot(m0, type="ICC", items = 7)
plot(m0, type="ICC", items = 8)
plot(m0, type="ICC", items = 9)
3.3 Person Item map
plotPImap(rm0,item.subset = "all", main = "Person-Item Map")
4 GLMM
1.根據Rasch model結果之試題難度將item依難度排序
library(dplyr)
dta1123 <- dta1123 %>% mutate(ID = factor(ID),
Item = factor(Item, levels = c("1", "7", "2", "3", "4", "5", "6", "8", "9")),
Response = factor(Response))2.以GLMM fit資料 (averaged Response was nested in person)
模式改為 Response(是否答對) ~ Items(試題難度,依單參數模型得到的難度排序)+ random effect(每個人平均答對情況, nested in person)
the random effect is very large, which indicated that the random effect cannot be overlooked.
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacksjPlot::tab_model(m01 <- glmer(Response ~ -1 + Item + (1 | ID),
data = dta1123, family = binomial),
show.obs=F, show.ngroups=F, transform=NULL,
show.se=T, show.r2=F,show.icc=F)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.289001 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?| Response | ||||
|---|---|---|---|---|
| Predictors | Log-Odds | std. Error | CI | p |
| Item [1] | 2.90 | 0.00 | 2.90 – 2.91 | <0.001 |
| Item [7] | 1.86 | 0.00 | 1.86 – 1.87 | <0.001 |
| Item [2] | 1.70 | 0.00 | 1.70 – 1.71 | <0.001 |
| Item [3] | 1.34 | 0.00 | 1.34 – 1.34 | <0.001 |
| Item [4] | 1.21 | 0.21 | 0.79 – 1.63 | <0.001 |
| Item [5] | 1.00 | 0.21 | 0.60 – 1.41 | <0.001 |
| Item [6] | 0.92 | 0.20 | 0.52 – 1.33 | <0.001 |
| Item [8] | 0.20 | 0.00 | 0.20 – 0.21 | <0.001 |
| Item [9] | -0.88 | 0.21 | -1.28 – -0.48 | <0.001 |
| Random Effects | ||||
| σ2 | 3.29 | |||
| τ00 ID | 1.16 | |||
summary(m01)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Response ~ -1 + Item + (1 | ID)
## Data: dta1123
##
## AIC BIC logLik deviance df.resid
## 1430.0 1482.1 -705.0 1410.0 1340
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.0215 -0.5795 0.3395 0.5451 3.6578
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.165 1.079
## Number of obs: 1350, groups: ID, 150
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## Item1 2.903013 0.001078 2691.885 < 2e-16 ***
## Item7 1.864884 0.001046 1783.178 < 2e-16 ***
## Item2 1.703233 0.001078 1579.343 < 2e-16 ***
## Item3 1.337304 0.001054 1268.791 < 2e-16 ***
## Item4 1.210188 0.212338 5.699 1.20e-08 ***
## Item5 1.001692 0.206531 4.850 1.23e-06 ***
## Item6 0.924094 0.204687 4.515 6.34e-06 ***
## Item8 0.203071 0.001054 192.669 < 2e-16 ***
## Item9 -0.877798 0.205296 -4.276 1.90e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## Item1 Item7 Item2 Item3 Item4 Item5 Item6 Item8
## Item7 -0.231
## Item2 -0.104 -0.231
## Item3 -0.183 0.095 -0.183
## Item4 0.000 -0.001 0.000 0.000
## Item5 0.002 0.004 0.002 0.003 0.101
## Item6 0.000 0.001 0.001 0.001 0.102 0.105
## Item8 0.183 -0.095 0.183 0.156 0.002 -0.002 0.001
## Item9 -0.001 -0.002 -0.001 -0.001 0.093 0.098 0.099 0.003
## convergence code: 0
## Model failed to converge with max|grad| = 0.289001 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?