Estimating polytomous IRT models in R
Tarid Wongvorachan
September 26th, 2020
In this practice, we will be using the science data set, which is a survey of the Consumer Protection and Perceptions of Science and Technology (1992). The data has seven items in total,
- Comfort: Science and technology are making our lives healthier, easier and more comfortable.
- Environment: Scientific and technological research cannot play an important role in protecting the environment and repairing it.
- Work: The application of science and new technology will make work more interesting.
- Future: Thanks to science and technology, there will be more opportunities for the future generations.
- Technology: New technology does not depend on basic scientific research.
- Industry: Scientific and technological research do not play an important role in industrial development.
- Benefit: The benefits of science are greater than any harmful effect it may have.
Importing Science data set
setwd("D:/Class Materials & Work/Summer 2020 practice/IRT_Poly")
science <- read.csv("science.csv", header=TRUE)
head(science)## Comfort Environment Work Future Technology Industry Benefit
## 1 4 4 4 3 4 3 2
## 2 3 4 3 3 3 3 3
## 3 3 2 2 2 4 4 3
## 4 3 3 2 2 4 4 3
## 5 3 1 4 4 2 3 1
## 6 4 3 4 3 3 4 3library(mirt)Partial & Generalized Partial Credit Models
We will be relying on mirt function to estimate item parameters for the Partial Credit Model. In the estimation process, we first begin with defining the latent trait “liking.science” based on the seven items in the data set.
Although we set itemtype as Rasch, the mirt function recognizes that the items have more than two response categories, and thus uses the Partial Credit Model to estimate item parameters.
model.pcm <- 'liking.science = 1-7'
results.pcm <- mirt(data=science, model=model.pcm, itemtype="Rasch", SE=TRUE, verbose=FALSE)
coef.pcm <- coef(results.pcm, IRTpars=TRUE, simplify=TRUE)We can save the estimated item parameters as a new data frame called items.pcm and print the results using the following commands:
items.pcm <- as.data.frame(coef.pcm$items)
print(items.pcm)## a b1 b2 b3
## Comfort 1 -2.269529 -2.2818104 1.1951986
## Environment 1 -1.471985 -0.5899431 0.2537860
## Work 1 -1.386737 -0.8029415 1.5160908
## Future 1 -2.002460 -1.1970343 0.9105307
## Technology 1 -1.977071 -0.6711969 0.3402364
## Industry 1 -1.969342 -1.4911110 0.1271229
## Benefit 1 -1.886126 -0.7456167 1.0713703In the output, we see a data frame of the estimated item parameters. The names of the items are printed on the left, and the first column shows the discrimination parameter, which is constrained to one for all items due to the PCM assumption (similar to Rasch’s).
The column b1 to b3 indicate the estimated thresholds or step parameters. There are three columns because the data set has four response options; Thus, the PCM estimates three threshold parameters (or steps), 1 to 2, 2 to 3, and 3 to 4.
Using the same commands above, we can also estimate item parameters for the Generalized Partial Credit Model (G-PCM), which allows a-parameter (item discrimination) of the model to be freely estimated for the items. Instead of itemtype="Rasch", we need to use itemtype="gpcm" to tell R that we want to estimate the “Generalized” Partial Credit Model this time.
model.gpcm <- 'liking.science = 1-7'
results.gpcm <- mirt(data=science, model=model.gpcm, itemtype="gpcm", SE=TRUE, verbose=FALSE)## EM cycles terminated after 500 iterations.coef.gpcm <- coef(results.gpcm, IRTpars=TRUE, simplify=TRUE)
items.gpcm <- as.data.frame(coef.gpcm$items)
print(items.gpcm)## a b1 b2 b3
## Comfort 0.87538955 -3.250347 -2.8567249 1.5202484
## Environment -0.03501584 32.393912 13.6362126 -3.5801431
## Work 0.83702216 -2.037048 -1.0364378 2.0640667
## Future 2.18522588 -2.086705 -0.9827835 0.8367083
## Technology -0.03801980 42.679630 14.3641172 -5.8043006
## Industry 0.12953982 -12.168178 -10.1545248 0.5402662
## Benefit 0.73353673 -2.875516 -1.0956900 1.6135753In the printed output, we can see that the column “a” is not fixed to 1 anymore. However, some of our threshold parameters have turned into very extreme numbers. This might be an indication that the G-PCM is probably not a good option for the science data set.
Rating Scale Model
This time, we will estimate item parameters using the Rating Scale Model (RSM). Similar to the PCM, this model has three threshold parameters for all items (because of four possible response options).
However, we will use the argument itemtype="rsm" instead for item parameter estimation. We then save the estimated item parameters in the items.rsm data frame and print the output using the print function.
model.rsm <- 'liking.science = 1-7'
results.rsm <- mirt(data=science, model=model.rsm, itemtype="rsm", verbose=FALSE)
coef.rsm <- coef(results.rsm, simplify=TRUE)
items.rsm <- as.data.frame(coef.rsm$items)
print(items.rsm)## a1 b1 b2 b3 c
## Comfort 1 -2.104223 -1.342868 0.5128852 0.0000000
## Environment 1 -2.104223 -1.342868 0.5128852 -0.3239308
## Work 1 -2.104223 -1.342868 0.5128852 -0.7130688
## Future 1 -2.104223 -1.342868 0.5128852 -0.2494849
## Technology 1 -2.104223 -1.342868 0.5128852 -0.2353408
## Industry 1 -2.104223 -1.342868 0.5128852 0.2559425
## Benefit 1 -2.104223 -1.342868 0.5128852 -0.5210050From the output above, the a-parameter is fixed to one as RSM is the polytomous form of the Rasch’s Model. The following b1 to b3 columns are threshold parameters (or steps), and the c column is the location parameter, which is the estimated mean of all threshold parameter for each item; Thus, Work is the easiest item, and Industry is the hardest item.
Graded Response Model
For Graded Response Model (GRM), we will use the argument itemtype = "graded" instead.
model.grm <- 'liking.science = 1-7'
results.grm <- mirt(data=science, model=model.grm, itemtype="graded", SE=TRUE, verbose=FALSE)
coef.grm <- coef(results.grm, IRTpars=TRUE, simplify=TRUE)
items.grm <- as.data.frame(coef.grm$items)
print(items.grm)## a b1 b2 b3
## Comfort 1.11185062 -4.438832 -2.4181197 1.3462905
## Environment 0.06229055 -40.634260 -13.4096480 11.5906085
## Work 1.17055741 -2.463243 -0.7617422 1.9050371
## Future 2.14668567 -2.342952 -0.9922735 0.8778317
## Technology 0.05882352 -51.633868 -16.2840234 12.6841326
## Industry 0.35143873 -10.532348 -5.1825245 1.0036080
## Benefit 1.13606611 -2.981573 -0.8887580 1.5019365The output indicates that this model has a poor fit as seen from extreme values in many item parameters.
Nominal Response Model
For Nominal Resmonse Model, we will use the argument itemtype = "nominal". However, do note that this model is not appropriate to the nature of our data set, but we will still implement it as our example.
model.nrm <- 'liking.science = 1-7'
results.nrm <- mirt(data=science, model=model.nrm, itemtype="nominal", SE=TRUE, verbose=FALSE)
coef.nrm <- coef(results.nrm, IRTpars=TRUE, simplify=TRUE)
items.nrm <- as.data.frame(coef.nrm$items)
print(items.nrm)## a1 a2 a3 a4 c1 c2
## Comfort 0.2479192 -0.4393110 -0.3188852 0.5102770 -2.1664458 -0.28557185
## Environment 0.6834829 -0.8405704 -0.7274799 0.8845673 -1.0695872 0.12548444
## Work 0.2551751 -0.6666461 -0.3092425 0.7207136 -0.8110968 0.27818439
## Future 0.2005983 -0.4946839 -0.4579307 0.7520163 -1.5027516 0.13920603
## Technology 0.6254204 -0.7270765 -0.6622835 0.7639396 -1.4773630 0.26284009
## Industry 0.9480691 -0.9591746 -1.1987801 1.2098856 -1.7992223 -0.09618228
## Benefit 0.2973991 -0.5917818 -0.4148605 0.7092432 -1.2440319 0.35709035
## c3 c4
## Comfort 1.8583455 0.5936722
## Environment 0.6538681 0.2902346
## Work 1.1068668 -0.5739544
## Future 1.2190123 0.1445333
## Technology 0.8337212 0.3808017
## Industry 1.0822432 0.8131614
## Benefit 1.0589349 -0.1719934Visualizing Polytomous IRT Models
We will use the plot-method function in the mirt package to visualize characteristics of the model (and items). Unlike dichotomous IRT models, polytomous IRT models have Option Characteristic Curves (OCC), which works as an extention from Item Characteristic Curves (ICC).
Because the items have more than two response categories, OCCs have multiple curves in the same plot. Each curve represents the probability of selecting a particular response option as a function of the latent trait.
In the plot function, we first enter the mirt-class object that stores the item parameters (i.e., result.pcm), and specify the type of plot using type = 'trace'.
Assuming we want to see the OCCs for items 1 and 3, we include which.items = c(1,3) in the script. Without specifying the item number, the default command is to plot OCCs for all items in a single plot.
main = "" eliminates printing a title above the plot, and we can choose the settings for the plot inside with par.settings. To make the curves thicker, we use lwd = 2 (or higher, if you want). In addition, we use lty=1:4 to request different types of lines for each OCC.
The auto.key command adjusts the legend settings. We want the legend to show which lines represent the OCCs. The following figure shows the OCCs for all seven items in the science data.
plot(results.pcm, type = 'trace',
main = "", par.settings = simpleTheme(lty=1:4,lwd=2),
auto.key=list(points=FALSE,lines=TRUE, columns=4)) 
In the figure, each item has four OCCs representing the four response options. The response categories are labeled as “P1” to “P4”. The OCCs are similar across all items, with P1 on the far left and P4 on the far right of the plot. The other OCCs in the middle also are ordered properly.
Next, we will draw a plot of item information function (IIF) to show the amount of information that each item can explain across different levels of the latent trait (Theta). We will use the same plot-method function with type = 'infotrace'. In the par.settings option, we only specify the thickness of the IIFs.
plot(results.pcm, type = 'infotrace',
main = "", par.settings = simpleTheme(lwd=2))
In addition to IIFs for particular items, we can show the Test Information Function (TIF) to display the total amount of information, as well as the distribution of conditional standard error of measurement (cSEM). We will draw a plot of TIF and cSEM using the range of -4 to 4 for the latent trait, with type = 'info' for the TIF and type = 'SE' for the cSEM plot.
plot(results.pcm, type = 'info', theta_lim = c(-4,4), lwd=2) 
plot(results.pcm, type = 'SE', theta_lim = c(-4,4), lwd=2) 
sessionInfo()## R version 4.0.0 (2020-04-24)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] mirt_1.32.1 lattice_0.20-41
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.4.6 cluster_2.1.0 knitr_1.28
## [4] magrittr_1.5 splines_4.0.0 MASS_7.3-51.6
## [7] Deriv_4.0.1 rlang_0.4.6 dcurver_0.9.1
## [10] stringr_1.4.0 highr_0.8 tools_4.0.0
## [13] parallel_4.0.0 grid_4.0.0 nlme_3.1-148
## [16] mgcv_1.8-31 xfun_0.14 vegan_2.5-6
## [19] htmltools_0.4.0 yaml_2.2.1 digest_0.6.25
## [22] permute_0.9-5 GPArotation_2014.11-1 Matrix_1.2-18
## [25] evaluate_0.14 rmarkdown_2.1 stringi_1.4.6
## [28] compiler_4.0.0Reference: https://rpubs.com/okanbulut/lecture2