Graded Response Model (GRM)

A tutorial in conducting Graded Response Model using mirt in R.

Ramadhan Dwi Marvianto (Unit Pengembangan Alat Psikodiagnostika UGM)
2023-06-15

Data Preparation

Library

library(mirt)       # Load data set and conduct GRM
library(dplyr)      # Tidy the codes
library(psych)      # Find skewness of the item distribution
library(rmarkdown)  # Display tables in good ones
library(knitr)      # Display tables in good ones
library(lavaan)     # Conduct CFA analysis

Load Data

data("Science")

Assumption: Unidimensionality

Descriptive Statistic

#### Creating a Table
desc.stat = matrix(nrow = 4, ncol = 7) %>% as.data.frame()
colnames(desc.stat) <- c("Mean", "SD", "Skewness", "Option 1",
                        "Option 2", "Option 3", "Option 4")
rownames(desc.stat) <- colnames(Science)

#### Assigning value to the Descriptive Statistic Table
##### Mean
desc.stat$Mean = apply(Science, 2, mean) %>% round(3)

##### SD
desc.stat$SD = apply(Science, 2, sd) %>% round(3)

##### Skewness
desc.stat$Skewness = apply(Science, 2, skew) %>% round(3)

##### Options Proportion
desc.stat[1,4:7] = round(table(Science[1])/sum(table(Science[1])),3)
desc.stat[2,4:7] = round(table(Science[2])/sum(table(Science[2])),3)
desc.stat[3,4:7] = round(table(Science[3])/sum(table(Science[3])),3)
desc.stat[4,4:7] = round(table(Science[4])/sum(table(Science[4])),3)

##### Displaying the table
paged_table(desc.stat)

Pearson Correlation Matrix

cor(Science) %>% round(3) %>% as.data.frame() %>% paged_table()

Confirmatory Factor Analysis

##### Model Specification
mod.cfa = 'F =~ Comfort + Work + Future + Benefit'

##### Running Analysis
my.cfa <- cfa(mod.cfa, data = Science, std.lv = T)

##### Model fit
fitMeasures(my.cfa, c("chisq", "df", "pvalue", 
                      "rmsea", "srmr", "gfi", "cfi", "tli")) %>% 
  as.data.frame() %>% t() %>% as.data.frame() %>% round(3) %>% paged_table()

Model Comparison of Polytomous IRT

Partial Credit Model (PCM)

#### Running analysis
my.pcm = mirt(Science, 1, itemtype = "Rasch", verbose = F)

#### Global fit
( res.pcm = M2(my.pcm, type = "C2") %>% round(3) %>% as.data.frame() %>% 
    paged_table() )

Generalzied Partial Credit Model (GPCM)

#### Running analysis
my.gpcm = mirt(Science, 1, itemtype = "gpcm", verbose = F)

#### Global fit
( res.gpcm = M2(my.gpcm, type = "C2") %>% round(3) %>% as.data.frame() %>% 
    paged_table() )

Graded Response Model (GRM)

#### Running analysis
my.grm = mirt(Science, 1, itemtype = "graded", verbose = F)

#### Global fit
( res.grm = M2(my.grm, type = "C2") %>% round(3) %>% as.data.frame() %>% 
    paged_table() )

Model Comparison

#### PCM vs GPCM
( res_01 = anova(my.pcm, my.gpcm) %>% as.data.frame() %>% round(3) %>% 
    paged_table() )
#### GPCM vs GRM
( res_02 = anova(my.gpcm, my.grm) %>% as.data.frame() %>% round(3) %>% 
    paged_table() )
#### Creating a Table
IRT_comparison = matrix(nrow = 3, ncol = 10) %>% as.data.frame()
colnames(IRT_comparison) = c(colnames(res.pcm), colnames(res_01)[4])
rownames(IRT_comparison) = c("PCM", "GPCM", "GRM")

#### Assigning value to the table
IRT_comparison[1,] = c(res.pcm, res_01[1,4])
IRT_comparison[2,] = c(res.gpcm, res_01[2,4])
IRT_comparison[3,] = c(res.grm, res_02[2,4])

#### Displaying the table
paged_table(IRT_comparison)

Assumption: Local Depedence using Q3

Q3 summary statistics:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -0.320  -0.249  -0.225  -0.190  -0.205   0.085 

        Comfort   Work Future Benefit
Comfort   1.000 -0.203 -0.252   0.085
Work     -0.203  1.000 -0.208  -0.242
Future   -0.252 -0.208  1.000  -0.320
Benefit   0.085 -0.242 -0.320   1.000

Graded Response Model’s parameter

Item location and discrmination

coef(my.grm, IRT = T, simplify = T)$items %>% round(3) %>% 
  as.data.frame() %>% paged_table()

Item fit

#### No Correction
itemfit(x = my.grm, fit_stats = "S_X2", which.items = 1:4)
     item  S_X2 df.S_X2 RMSEA.S_X2 p.S_X2
1 Comfort 4.774       6          0  0.573
2    Work 8.756       9          0  0.460
3  Future 7.543       8          0  0.479
4 Benefit 9.928      11          0  0.537
#### Correction using False Discovery Rate (Benjamini and Hocherg)
itemfit(x = my.grm, fit_stats = "S_X2", which.items = 1:4, p.adjust = "BH")
     item  S_X2 df.S_X2 RMSEA.S_X2 p.S_X2
1 Comfort 4.774       6          0  0.573
2    Work 8.756       9          0  0.573
3  Future 7.543       8          0  0.573
4 Benefit 9.928      11          0  0.573
#### Correction of Bonferonni
itemfit(x = my.grm, fit_stats = "S_X2", which.items = 1:4, p.adjust = "bonferroni")
     item  S_X2 df.S_X2 RMSEA.S_X2 p.S_X2
1 Comfort 4.774       6          0      1
2    Work 8.756       9          0      1
3  Future 7.543       8          0      1
4 Benefit 9.928      11          0      1

Option response function (ORF)

update(itemplot(my.grm, 1, type = "trace"), main = colnames(Science)[1])
update(itemplot(my.grm, 2, type = "trace"), main = colnames(Science)[2])
update(itemplot(my.grm, 3, type = "trace"), main = colnames(Science)[3])
update(itemplot(my.grm, 4, type = "trace"), main = colnames(Science)[4])

Category characteristic curves or cumulative probability curve

update(itemplot(my.grm, 1, type = "threshold"), main = colnames(Science)[1])
update(itemplot(my.grm, 2, type = "threshold"), main = colnames(Science)[2])
update(itemplot(my.grm, 3, type = "threshold"), main = colnames(Science)[3])
update(itemplot(my.grm, 4, type = "threshold"), main = colnames(Science)[4])

Item information function (IIF)

update(itemplot(my.grm, 1, type = "infoSE"), main = colnames(Science)[1])
update(itemplot(my.grm, 2, type = "infoSE"), main = colnames(Science)[2])
update(itemplot(my.grm, 3, type = "infoSE"), main = colnames(Science)[3])
update(itemplot(my.grm, 4, type = "infoSE"), main = colnames(Science)[4])

ORF and IIF in a single plot

update(itemplot(my.grm, 1, type = "infotrace"), main = colnames(Science)[1])
update(itemplot(my.grm, 2, type = "infotrace"), main = colnames(Science)[2])
update(itemplot(my.grm, 3, type = "infotrace"), main = colnames(Science)[3])
update(itemplot(my.grm, 4, type = "infotrace"), main = colnames(Science)[4])

Test information function (TIF) and conditional standard error

plot(my.grm, type = 'infoSE', theta_lim = c(-6, 6), 
     main="Test Information Function and Conditional Standard Error")

Conditional reliability

plot(my.grm, type = 'rxx', theta_lim = c(-6, 6), 
     main="Conditional Reliabilty Curve" )

Test expected score

plot(my.grm, type = 'score', theta_lim = c(-6, 6), 
     main = "Test Expected Score")

Reliability estimation

Marginal reliability

marginal_rxx(my.grm) %>% round(3)
[1] 0.669

Empirical reliability

( EAP.est <- fscores(object = my.grm, method = "EAP", full.scores.SE = T) %>% 
  round(3) %>% empirical_rxx() %>% round(3) )
   F1 
0.667

Summed Score Expected A Priori

fscores(my.grm, method = 'EAPsum', full.scores = FALSE) %>% round(3)
      df       X2       p.X2    rxx_F1
stats 10 16.31202 0.09104204 0.6241528
   Sum.Scores     F1 SE_F1 observed expected std.res
4           4 -2.749 0.629        2    0.124   5.324
5           5 -2.431 0.617        1    0.790   0.236
6           6 -2.081 0.610        2    2.760   0.457
7           7 -1.718 0.602        1    7.252   2.322
8           8 -1.364 0.598       11   15.893   1.227
9           9 -1.012 0.604       32   29.674   0.427
10         10 -0.649 0.610       58   48.363   1.386
11         11 -0.287 0.605       70   68.485   0.183
12         12  0.082 0.600       91   80.553   1.164
13         13  0.487 0.613       56   65.668   1.193
14         14  0.934 0.617       36   42.637   1.016
15         15  1.384 0.622       20   22.331   0.493
16         16  1.854 0.654       12    7.470   1.657