Item Response Theory practice
Tarid Wongvorachan
September 23th, 2020
We will be using the package ltm with two test data sets, one dichotomous and one polytomous.
An example with Dichotomous data
setwd("D:/Class Materials & Work/Summer 2020 practice/IRT")
#Loading required packages
library(ltm)
library(mirt)We will be using the data from Law School Admission Test (LSAT), N = 1000, 5 items. The data can be called with data(LSAT) or data <- data.frame(LSAT).
As an initial step, we can use ltm::descript for descriptive statistics.
data(LSAT)
descript(LSAT)##
## Descriptive statistics for the 'LSAT' data-set
##
## Sample:
## 5 items and 1000 sample units; 0 missing values
##
## Proportions for each level of response:
## 0 1 logit
## Item 1 0.076 0.924 2.4980
## Item 2 0.291 0.709 0.8905
## Item 3 0.447 0.553 0.2128
## Item 4 0.237 0.763 1.1692
## Item 5 0.130 0.870 1.9010
##
##
## Frequencies of total scores:
## 0 1 2 3 4 5
## Freq 3 20 85 237 357 298
##
##
## Point Biserial correlation with Total Score:
## Included Excluded
## Item 1 0.3620 0.1128
## Item 2 0.5668 0.1532
## Item 3 0.6184 0.1728
## Item 4 0.5344 0.1444
## Item 5 0.4354 0.1216
##
##
## Cronbach's alpha:
## value
## All Items 0.2950
## Excluding Item 1 0.2754
## Excluding Item 2 0.2376
## Excluding Item 3 0.2168
## Excluding Item 4 0.2459
## Excluding Item 5 0.2663
##
##
## Pairwise Associations:
## Item i Item j p.value
## 1 1 5 0.565
## 2 1 4 0.208
## 3 3 5 0.113
## 4 2 4 0.059
## 5 1 2 0.028
## 6 2 5 0.009
## 7 1 3 0.003
## 8 4 5 0.002
## 9 3 4 7e-04
## 10 2 3 4e-04From the output above, inspection of non significant results can be used to reveal ‘problematic’ items in pairwise association. Latent variable models assume that t the high associations between items can be explained by a set of latent variables, so any pair of items that is not related to each other violates this assumption.
Additionally, Item 1 seems to be the easist item as seen from its highest proportion of correct response.
We will fit the original Rasch model, which fixes the item discrimination (aka a parameter) of all items to 1, to the data. The ltm::rasch() assumes equal a-parameter across items with an estimated value. In order to impose the constraint = 1, the constraint argument is used.
This argument accepts a two-column matrix where the first column denotes the parameter and the second column indicates the value at which the corresponding parameter should be fixed.
fit1 <- rasch(LSAT, constraint = cbind(length(LSAT) + 1, 1))
summary(fit1)##
## Call:
## rasch(data = LSAT, constraint = cbind(length(LSAT) + 1, 1))
##
## Model Summary:
## log.Lik AIC BIC
## -2473.054 4956.108 4980.646
##
## Coefficients:
## value std.err z.vals
## Dffclt.Item 1 -2.8720 0.1287 -22.3066
## Dffclt.Item 2 -1.0630 0.0821 -12.9458
## Dffclt.Item 3 -0.2576 0.0766 -3.3635
## Dffclt.Item 4 -1.3881 0.0865 -16.0478
## Dffclt.Item 5 -2.2188 0.1048 -21.1660
## Dscrmn 1.0000 NA NA
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Convergence: 0
## max(|grad|): 6.3e-05
## quasi-Newton: BFGSThe results of the descriptive analysis are also validated by the model fit, where items 3 and 1 are the most difficult and the easiest respectively (the lower the b-parameter value, the easier).
The parameter estimates can be transformed to probability estimates using the coef() method:
coef(fit1, prob = TRUE, order = TRUE)## Dffclt Dscrmn P(x=1|z=0)
## Item 1 -2.8719712 1 0.9464434
## Item 5 -2.2187785 1 0.9019232
## Item 4 -1.3880588 1 0.8002822
## Item 2 -1.0630294 1 0.7432690
## Item 3 -0.2576109 1 0.5640489The last column denotes the probability of a positive response to the ith item for the average individual. The argument order = T indicates the output to sort the items according to the difficulty estimates.
In order to check the fit of the model to the data, the argument GoF.rasch() and margins() can be used. The former argument performs a parametric Bootstrap goodness-of-fit test using Pearson’s Chi-square statistics, while the latter examines the two- and three-way chi-square residual analysis.
GoF.rasch(fit1, B = 199)##
## Bootstrap Goodness-of-Fit using Pearson chi-squared
##
## Call:
## rasch(data = LSAT, constraint = cbind(length(LSAT) + 1, 1))
##
## Tobs: 30.6
## # data-sets: 200
## p-value: 0.295Based on 200 data-sets, the non significant p-value suggests an acceptable fit of the model. The null hypothesis states that the observed data and the model fit with each other (no differences).
Now, for two-way margin analysis.
margins(fit1)##
## Call:
## rasch(data = LSAT, constraint = cbind(length(LSAT) + 1, 1))
##
## Fit on the Two-Way Margins
##
## Response: (0,0)
## Item i Item j Obs Exp (O-E)^2/E
## 1 2 4 81 98.69 3.17
## 2 1 5 12 18.45 2.25
## 3 3 5 67 80.04 2.12
##
## Response: (1,0)
## Item i Item j Obs Exp (O-E)^2/E
## 1 3 5 63 51.62 2.51
## 2 2 4 156 139.78 1.88
## 3 3 4 108 99.42 0.74
##
## Response: (0,1)
## Item i Item j Obs Exp (O-E)^2/E
## 1 2 4 210 193.47 1.41
## 2 2 3 135 125.07 0.79
## 3 1 4 53 47.24 0.70
##
## Response: (1,1)
## Item i Item j Obs Exp (O-E)^2/E
## 1 2 4 553 568.06 0.40
## 2 3 5 490 501.43 0.26
## 3 2 3 418 427.98 0.23From the above output, using the 3.5 rule of thumb, the value of all two-way combinations are below the cut-off (same way of how statistical hypothesis works) and therefore indicate a good fit to the two-way margins.
Next, we will examine the fit to the three-way margins:
margins(fit1, type = "three-way", nprint = 2) #nprint returns 2 highest residual values for each combinations##
## Call:
## rasch(data = LSAT, constraint = cbind(length(LSAT) + 1, 1))
##
## Fit on the Three-Way Margins
##
## Response: (0,0,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 2 3 4 48 66.07 4.94 ***
## 2 1 3 5 6 13.58 4.23 ***
##
## Response: (1,0,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 4 70 82.01 1.76
## 2 2 4 5 28 22.75 1.21
##
## Response: (0,1,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 5 3 7.73 2.90
## 2 3 4 5 37 45.58 1.61
##
## Response: (1,1,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 3 4 5 48 36.91 3.33
## 2 1 2 4 144 126.35 2.47
##
## Response: (0,0,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 3 5 41 34.58 1.19
## 2 2 4 5 64 72.26 0.94
##
## Response: (1,0,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 4 190 174.87 1.31
## 2 1 2 3 126 114.66 1.12
##
## Response: (0,1,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 5 42 34.35 1.70
## 2 1 4 5 46 38.23 1.58
##
## Response: (1,1,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 3 4 5 397 416.73 0.93
## 2 2 3 4 343 361.18 0.91
##
## '***' denotes a chi-squared residual greater than 3.5The three-way margins suggest a problematic fit for two triplets of items, both containing item 3.
We can try fitting the unconstrained version of Rasch model to see the difference. This time, no need for the constraint argument:
fit2 <- rasch(LSAT)
summary(fit2)##
## Call:
## rasch(data = LSAT)
##
## Model Summary:
## log.Lik AIC BIC
## -2466.938 4945.875 4975.322
##
## Coefficients:
## value std.err z.vals
## Dffclt.Item 1 -3.6153 0.3266 -11.0680
## Dffclt.Item 2 -1.3224 0.1422 -9.3009
## Dffclt.Item 3 -0.3176 0.0977 -3.2518
## Dffclt.Item 4 -1.7301 0.1691 -10.2290
## Dffclt.Item 5 -2.7802 0.2510 -11.0743
## Dscrmn 0.7551 0.0694 10.8757
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Convergence: 0
## max(|grad|): 2.9e-05
## quasi-Newton: BFGSThe output suggests that the discrimination parameter is different from 1. The difference can be tested with a likelihood ratio test using anova()
anova(fit1, fit2)##
## Likelihood Ratio Table
## AIC BIC log.Lik LRT df p.value
## fit1 4956.11 4980.65 -2473.05
## fit2 4945.88 4975.32 -2466.94 12.23 1 <0.001By comparing model summary of the constrained and unconstrained version, the latter is more suitable for the LSAT data due to its smaller AIC and BIC. We can double check thie result by testing the unconstrained model fitto the three-way margins, which yields a problematic fit with the constrained model.
margins(fit2, type = "three-way", nprint = 2)##
## Call:
## rasch(data = LSAT)
##
## Fit on the Three-Way Margins
##
## Response: (0,0,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 3 5 6 9.40 1.23
## 2 3 4 5 30 25.85 0.67
##
## Response: (1,0,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 2 4 5 28 22.75 1.21
## 2 2 3 4 81 74.44 0.58
##
## Response: (0,1,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 5 3 7.58 2.76
## 2 1 3 4 5 9.21 1.92
##
## Response: (1,1,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 2 4 5 51 57.49 0.73
## 2 3 4 5 48 42.75 0.64
##
## Response: (0,0,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 3 5 41 33.07 1.90
## 2 2 3 4 108 101.28 0.45
##
## Response: (1,0,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 2 3 4 210 218.91 0.36
## 2 1 2 4 190 185.56 0.11
##
## Response: (0,1,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 3 5 23 28.38 1.02
## 2 1 4 5 46 42.51 0.29
##
## Response: (1,1,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 4 520 526.36 0.08
## 2 1 2 3 398 393.30 0.06The new three-way margins suggests a good fit with the unconstrained Rasch model.
Finally, we investigate two more possible extensions of the unconstrained Rasch model, the two-parameter logistic model that assumes a different discrimination parameter per item, and Rasch model that incorporates a guessing parameter.
Two-parameter model can also be fitted with ltm() The formula of ltm() is two-sided, where its left is either a data frame or a matrix, and its right allows only z1 and/or z2 latent variables with the latter one serves in the case of interaction.
#2PL (unconstrained a and b)
fit3 <- ltm(LSAT ~ z1)
summary(fit3)##
## Call:
## ltm(formula = LSAT ~ z1)
##
## Model Summary:
## log.Lik AIC BIC
## -2466.653 4953.307 5002.384
##
## Coefficients:
## value std.err z.vals
## Dffclt.Item 1 -3.3597 0.8669 -3.8754
## Dffclt.Item 2 -1.3696 0.3073 -4.4565
## Dffclt.Item 3 -0.2799 0.0997 -2.8083
## Dffclt.Item 4 -1.8659 0.4341 -4.2982
## Dffclt.Item 5 -3.1236 0.8700 -3.5904
## Dscrmn.Item 1 0.8254 0.2581 3.1983
## Dscrmn.Item 2 0.7229 0.1867 3.8721
## Dscrmn.Item 3 0.8905 0.2326 3.8281
## Dscrmn.Item 4 0.6886 0.1852 3.7186
## Dscrmn.Item 5 0.6575 0.2100 3.1306
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Convergence: 0
## max(|grad|): 0.024
## quasi-Newton: BFGS#compare
anova(fit2, fit3)##
## Likelihood Ratio Table
## AIC BIC log.Lik LRT df p.value
## fit2 4945.88 4975.32 -2466.94
## fit3 4953.31 5002.38 -2466.65 0.57 4 0.967Next, we can fit the Rasch model with c parameter by tpm(). The syntax is very similar to rasch(), and it also allows one fit three-parameter logistic model as well.
#Rasch w/ c param
fit4 <- tpm(LSAT, type = "rasch", max.guessing = 1)
summary(fit4)##
## Call:
## tpm(data = LSAT, type = "rasch", max.guessing = 1)
##
## Model Summary:
## log.Lik AIC BIC
## -2466.731 4955.461 5009.447
##
## Coefficients:
## value std.err z.vals
## Gussng.Item 1 0.0830 0.8652 0.0959
## Gussng.Item 2 0.1962 0.3545 0.5535
## Gussng.Item 3 0.0081 0.0817 0.0994
## Gussng.Item 4 0.2565 0.4776 0.5372
## Gussng.Item 5 0.4957 0.4839 1.0242
## Dffclt.Item 1 -3.1765 1.5090 -2.1050
## Dffclt.Item 2 -0.7723 1.0281 -0.7513
## Dffclt.Item 3 -0.2707 0.2375 -1.1395
## Dffclt.Item 4 -1.0332 1.4031 -0.7364
## Dffclt.Item 5 -1.4332 1.9055 -0.7521
## Dscrmn 0.8459 0.1851 4.5711
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Optimizer: optim (BFGS)
## Convergence: 0
## max(|grad|): 0.073anova(fit2, fit4)##
## Likelihood Ratio Table
## AIC BIC log.Lik LRT df p.value
## fit2 4945.88 4975.32 -2466.94
## fit4 4955.46 5009.45 -2466.73 0.41 5 0.995The max.guessing argument specifies the upper bound for the guessing ( c ) parameters. FOr both 2PL and Rasch w/ c param, the data suggests no statistical differences with the unconstrained Rasch ( p = 0.967 and p = 0.995 respectively), as well as larger AIC and BIC; Therefore, they are not required.
We can then produce the Item Characteristic, the Item Information, and the Test Information Curves of the unconstrained Rasch with plot().
par(mfrow = c(2, 2)) #mar to adjust the margin
#Item characterisrics curve by default
plot(fit2, legend = TRUE, cx = "bottomright", lwd = 2,
cex.main = 1.5, cex.lab = 1.3, cex = 0.6)
#Item Information Curve
plot(fit2, type = "IIC", annot = F, items = NULL, lwd = 2,
cex.main = 1.5, cex.lab = 1.3) #items = NULL to plot all items
plot(fit2, type = "IIC", items = 0, lwd = 2,
cex.main = 1.5, cex.lab = 1.3) #items = NULL to plot Item Information Curve
plot(0:1, 0:1, type = "n", ann = FALSE, axes = FALSE)
info1 <- information(fit2, c(-4, 0)) #Item information on the negative area of theta (-4 to 0)
info2 <- information(fit2, c(0, 4)) #Item information on the positive area of theta (0 to 4)
text(0.5, 0.5, labels = paste("Total Information:", round(info1$InfoTotal, 3),
"\n\nInformation in (-4, 0):", round(info1$InfoRange, 3),
paste("(", round(100 * info1$PropRange, 2), "%)", sep = ""),
"\n\nInformation in (0, 4):", round(info2$InfoRange, 3),
paste("(", round(100 * info2$PropRange, 2), "%)", sep = "")))
The argument information() produces the area under the Test or Item Information Curves in a specified interval. From the test information curve,the items asked in LSAT mainly provide information for respondents with low ability.
In particular, the amount of Test Information for ability levels in the interval (-4 - 0) is almost 60%, and the item that seems to distinguish between respondents with higher ability levels is the third one.
Finally, the ability estimates can be obtained using the factor.scores().
factor.scores(fit2)##
## Call:
## rasch(data = LSAT)
##
## Scoring Method: Empirical Bayes
##
## Factor-Scores for observed response patterns:
## Item 1 Item 2 Item 3 Item 4 Item 5 Obs Exp z1 se.z1
## 1 0 0 0 0 0 3 2.364 -1.910 0.790
## 2 0 0 0 0 1 6 5.468 -1.439 0.793
## 3 0 0 0 1 0 2 2.474 -1.439 0.793
## 4 0 0 0 1 1 11 8.249 -0.959 0.801
## 5 0 0 1 0 0 1 0.852 -1.439 0.793
## 6 0 0 1 0 1 1 2.839 -0.959 0.801
## 7 0 0 1 1 0 3 1.285 -0.959 0.801
## 8 0 0 1 1 1 4 6.222 -0.466 0.816
## 9 0 1 0 0 0 1 1.819 -1.439 0.793
## 10 0 1 0 0 1 8 6.063 -0.959 0.801
## 11 0 1 0 1 1 16 13.288 -0.466 0.816
## 12 0 1 1 0 1 3 4.574 -0.466 0.816
## 13 0 1 1 1 0 2 2.070 -0.466 0.816
## 14 0 1 1 1 1 15 14.749 0.049 0.836
## 15 1 0 0 0 0 10 10.273 -1.439 0.793
## 16 1 0 0 0 1 29 34.249 -0.959 0.801
## 17 1 0 0 1 0 14 15.498 -0.959 0.801
## 18 1 0 0 1 1 81 75.060 -0.466 0.816
## 19 1 0 1 0 0 3 5.334 -0.959 0.801
## 20 1 0 1 0 1 28 25.834 -0.466 0.816
## 21 1 0 1 1 0 15 11.690 -0.466 0.816
## 22 1 0 1 1 1 80 83.310 0.049 0.836
## 23 1 1 0 0 0 16 11.391 -0.959 0.801
## 24 1 1 0 0 1 56 55.171 -0.466 0.816
## 25 1 1 0 1 0 21 24.965 -0.466 0.816
## 26 1 1 0 1 1 173 177.918 0.049 0.836
## 27 1 1 1 0 0 11 8.592 -0.466 0.816
## 28 1 1 1 0 1 61 61.235 0.049 0.836
## 29 1 1 1 1 0 28 27.709 0.049 0.836
## 30 1 1 1 1 1 298 295.767 0.593 0.862By default, factor.scores() produces ability estimates for the observed response patterns (every combination available); if ability estimates are required for non observed or specific response patterns, these could be specified using the resp.patterns argument, such as:
factor.scores(fit2, resp.patterns = rbind(c(0,1,1,0,0), c(0,1,0,1,0)))##
## Call:
## rasch(data = LSAT)
##
## Scoring Method: Empirical Bayes
##
## Factor-Scores for specified response patterns:
## Item 1 Item 2 Item 3 Item 4 Item 5 Obs Exp z1 se.z1
## 1 0 1 1 0 0 0 0.944 -0.959 0.801
## 2 0 1 0 1 0 0 2.744 -0.959 0.801An example with Polytomous (or ordinal) data
The data we consider here come from the Environment section of the 1990 British Social Attitudes Survey, N = 291, 6 items, 3 ordinal response options. The dataframe can be loaded with data(Environment).
data(Environment)
descript(Environment)##
## Descriptive statistics for the 'Environment' data-set
##
## Sample:
## 6 items and 291 sample units; 0 missing values
##
## Proportions for each level of response:
## very concerned slightly concerned not very concerned
## LeadPetrol 0.6151 0.3265 0.0584
## RiverSea 0.8007 0.1753 0.0241
## RadioWaste 0.7457 0.1924 0.0619
## AirPollution 0.6495 0.3196 0.0309
## Chemicals 0.7491 0.1924 0.0584
## Nuclear 0.5155 0.3265 0.1581
##
##
## Frequencies of total scores:
## 6 7 8 9 10 11 12 13 14 15 16 17 18
## Freq 96 51 37 27 26 18 13 7 6 6 1 1 2
##
##
## Cronbach's alpha:
## value
## All Items 0.8215
## Excluding LeadPetrol 0.8218
## Excluding RiverSea 0.7990
## Excluding RadioWaste 0.7767
## Excluding AirPollution 0.7751
## Excluding Chemicals 0.7790
## Excluding Nuclear 0.8058
##
##
## Pairwise Associations:
## Item i Item j p.value
## 1 1 2 0.001
## 2 1 3 0.001
## 3 1 4 0.001
## 4 1 5 0.001
## 5 1 6 0.001
## 6 2 3 0.001
## 7 2 4 0.001
## 8 2 5 0.001
## 9 2 6 0.001
## 10 3 4 0.001From the descriptive output, the first response, “very concerned”, has the highest frequency. The p-values for the pairwise associations indicate significant associations between all items.
An alternative method to explore the degree of association between pairs of items can be done with rcor.test() for non-parametric correlation coefficients.
rcor.test(Environment, method = "kendall")##
## LeadPetrol RiverSea RadioWaste AirPollution Chemicals Nuclear
## LeadPetrol ***** 0.385 0.260 0.457 0.305 0.279
## RiverSea <0.001 ***** 0.399 0.548 0.403 0.320
## RadioWaste <0.001 <0.001 ***** 0.506 0.623 0.484
## AirPollution <0.001 <0.001 <0.001 ***** 0.504 0.382
## Chemicals <0.001 <0.001 <0.001 <0.001 ***** 0.463
## Nuclear <0.001 <0.001 <0.001 <0.001 <0.001 *****
##
## upper diagonal part contains correlation coefficient estimates
## lower diagonal part contains corresponding p-valuesThe rcor.test() provides two options for nonparametric calculation, Kendall’s Tau and Spearman’s rho in method argument.
Initially, we will fit the constrained version of Graded Response Model (GRM) that assumes equal a parameter across items (similar to Rasch model). The model is fitted by grm() as follows:
fit5 <- grm(Environment, constrained = TRUE)
summary(fit5)##
## Call:
## grm(data = Environment, constrained = TRUE)
##
## Model Summary:
## log.Lik AIC BIC
## -1106.193 2238.386 2286.139
##
## Coefficients:
## $LeadPetrol
## value
## Extrmt1 0.395
## Extrmt2 1.988
## Dscrmn 2.218
##
## $RiverSea
## value
## Extrmt1 1.060
## Extrmt2 2.560
## Dscrmn 2.218
##
## $RadioWaste
## value
## Extrmt1 0.832
## Extrmt2 1.997
## Dscrmn 2.218
##
## $AirPollution
## value
## Extrmt1 0.483
## Extrmt2 2.448
## Dscrmn 2.218
##
## $Chemicals
## value
## Extrmt1 0.855
## Extrmt2 2.048
## Dscrmn 2.218
##
## $Nuclear
## value
## Extrmt1 0.062
## Extrmt2 1.266
## Dscrmn 2.218
##
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Convergence: 0
## max(|grad|): 0.0049
## quasi-Newton: BFGSIf standard errors for the parameter estimates are required, you can add the argument Hessian = T to the function grm().
Similarly to our dichotomous case, the fit of the model can be checked using margins() for two-way margins.
margins(fit5)##
## Call:
## grm(data = Environment, constrained = TRUE)
##
## Fit on the Two-Way Margins
##
## LeadPetrol RiverSea RadioWaste AirPollution Chemicals Nuclear
## LeadPetrol - 10.03 9.98 5.19 7.85 16.93
## RiverSea - 5.06 17.12 2.56 7.14
## RadioWaste - 6.78 20.60 12.09
## AirPollution - 4.49 4.57
## Chemicals - 3.85
## Nuclear -The upper diagonal part of the outputcontains the residuals, and the lower diagonal part indicates the pairs for which the residuals exceed the threshold value. The output above looks good. Next, we will try with the three-way margins.
margins(fit5, type = "three")##
## Call:
## grm(data = Environment, constrained = TRUE)
##
## Fit on the Three-Way Margins
##
## Item i Item j Item k (O-E)^2/E
## 1 1 2 3 28.52
## 2 1 2 4 34.26
## 3 1 2 5 29.91
## 4 1 2 6 42.74
## 5 1 3 4 33.03
## 6 1 3 5 66.72
## 7 1 3 6 65.31
## 8 1 4 5 25.48
## 9 1 4 6 34.46
## 10 1 5 6 39.49
## 11 2 3 4 29.63
## 12 2 3 5 37.74
## 13 2 3 6 32.50
## 14 2 4 5 27.08
## 15 2 4 6 36.77
## 16 2 5 6 19.49
## 17 3 4 5 38.99
## 18 3 4 6 26.91
## 19 3 5 6 39.62
## 20 4 5 6 22.25Both the two- and three-way residuals show a good fit of the constrained model to the data, but checking the fit of the model in the margins does not correspond to an overall goodness-of-fit test. As a result, we will fit the unconstrained version of the GRM as well.
fit6 <- grm(Environment) #unconstrained GRM
summary(fit6)##
## Call:
## grm(data = Environment)
##
## Model Summary:
## log.Lik AIC BIC
## -1090.404 2216.807 2282.927
##
## Coefficients:
## $LeadPetrol
## value
## Extrmt1 0.487
## Extrmt2 2.584
## Dscrmn 1.378
##
## $RiverSea
## value
## Extrmt1 1.058
## Extrmt2 2.499
## Dscrmn 2.341
##
## $RadioWaste
## value
## Extrmt1 0.779
## Extrmt2 1.793
## Dscrmn 3.123
##
## $AirPollution
## value
## Extrmt1 0.457
## Extrmt2 2.157
## Dscrmn 3.283
##
## $Chemicals
## value
## Extrmt1 0.809
## Extrmt2 1.868
## Dscrmn 2.947
##
## $Nuclear
## value
## Extrmt1 0.073
## Extrmt2 1.427
## Dscrmn 1.761
##
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Convergence: 0
## max(|grad|): 0.003
## quasi-Newton: BFGSWe can use a likelihood ratio test to check if the unconstrained version GRM is better than its constrained one.
anova(fit5, fit6)##
## Likelihood Ratio Table
## AIC BIC log.Lik LRT df p.value
## fit5 2238.39 2286.14 -1106.19
## fit6 2216.81 2282.93 -1090.40 31.58 5 <0.001The LRT indicates that the unconstrained GRM is preferable for the Environment data. We can plot the Item Characteristic Curve (ICC) of all 6 items, Item Information Curve (IIC), and Test Information Curve (TIC) below.
par(mar=c(2,6,2,2), mfrow = c(2, 2))
plot(fit6, lwd = 2, cex = 0.6, legend = TRUE, cx = "left",
xlab = "Latent Trait", cex.main = 0.7, cex.lab = 0.7, cex.axis = 1,)
plot(fit6, type = "IIC", lwd = 2, cex = 0.6, legend = TRUE, cx = "topleft",
xlab = "Latent Trait", cex.main = 0.7, cex.lab = 0.7, cex.axis = 1)
plot(fit6, type = "IIC", items = 0, lwd = 2, xlab = "Latent Trait",
cex.main = 0.7, cex.lab = 0.7, cex.axis = 1)
info3 <- information(fit6, c(-4, 0))
info4 <- information(fit6, c(0, 4))
text(-1.9, 8, labels = paste("Information in (-4, 0):",
paste(round(100 * info3$PropRange, 1), "%", sep = ""),
"\n\nInformation in (0, 4):",
paste(round(100 * info4$PropRange, 1), "%", sep = "")), cex = 0.5)
From the Item Characteristic Curve, we observe that there is low probability of endorsing the g the first option, “very concerned”, for relatively high latent trait levels, which means that the questions asked are not considered as major environmental issues by the respondent.
The Test Information Curve also tells us that the test provides 89% of the total information for high latent trait levels.
Finally, the Item Information Curve indicates that items in LeadPetrol and Nuclear provide little information in the whole latent trait continuum. We can check this in detail using information().
information(fit6, c(-4, 4))##
## Call:
## grm(data = Environment)
##
## Total Information = 26.97
## Information in (-4, 4) = 26.7 (98.97%)
## Based on all the itemsinformation(fit6, c(-4, 4), items = c(1, 6)) #for item 1 and 6##
## Call:
## grm(data = Environment)
##
## Total Information = 5.36
## Information in (-4, 4) = 5.17 (96.38%)
## Based on items 1, 6We observe that item 1 and 6 provide only the 20.36% of the total information (from the total of 26.91); Thus, thus they could probably be excluded from a similar future study.
Finally, a useful comparison is to plot the ICC of each response option separately.
par(mar=c(2,6,2,2), mfrow = c(2, 2))
plot(fit6, category = 1, lwd = 2, cex = 0.7, legend = TRUE, cx = -4.5,
cy = 0.85, xlab = "Latent Trait", cex.main = 0.7, cex.lab = 0.7,
cex.axis = 1)
for (ctg in 2:3)
{
plot(fit6, category = ctg, lwd = 2, cex = 0.5, annot = FALSE,
xlab = "Latent Trait", cex.main = 0.7, cex.lab = 0.7,
cex.axis = 1)
}
From the plot, the response option for RadioWaste and Chemicals have nearly identical characteristic curves for all categories, indicating that these two items are probably regarded to have the same effect on the environment.
Extra features of ltm().
Aside from fitting common models to dichotomous and polytomous data, the package can either fits latent trait models with one or two latent variables, which allows for situations in which correlation between the latent variables is plausible.
We will use the data from 1990 Workplace Industrial Relation Survey (WIRS), which has possible relationships between latent variables.
The data is dichotomous, with N = 1005, and 6 items, and it investigates formal and informal latent variables.
However, the fit of two-factor model to the three-way margins is not successful, with some triplets of items having high residual values. Here, we can do a simple two-factor model with an interaction to take relationships between the variables into account.
data(WIRS)
fit7 <- ltm(WIRS ~ z1 + z2) #Normal two-factor model
fit8 <- ltm(WIRS ~ z1 * z2) #Two-factor model w/Interaction
anova(fit7, fit8)##
## Likelihood Ratio Table
## AIC BIC log.Lik LRT df p.value
## fit7 6719.08 6807.51 -3341.54
## fit8 6634.42 6752.33 -3293.21 96.66 6 <0.001The p value suggests that the model with interaction provides a better fit to the data. We can also double check the fit with three-way margins, which yields acceptable values for all items.
fit8##
## Call:
## ltm(formula = WIRS ~ z1 * z2)
##
## Coefficients:
## (Intercept) z1 z2 z1:z2
## Item 1 -1.097 0.271 -3.317 -1.160
## Item 2 0.944 1.091 2.173 2.955
## Item 3 -1.458 1.811 -0.314 0.543
## Item 4 -1.465 1.116 0.556 0.269
## Item 5 -1.083 2.133 -0.466 1.147
## Item 6 -2.753 1.698 -0.935 1.072
##
## Log.Lik: -3293.212margins(fit8, type = "three-way", nprint = 2)##
## Call:
## ltm(formula = WIRS ~ z1 * z2)
##
## Fit on the Three-Way Margins
##
## Response: (0,0,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 4 171 163.95 0.30
## 2 4 5 6 494 500.20 0.08
##
## Response: (1,0,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 4 5 6 102 94.96 0.52
## 2 1 2 4 187 191.93 0.13
##
## Response: (0,1,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 3 4 5 85 78.37 0.56
## 2 2 3 4 84 77.69 0.51
##
## Response: (1,1,0)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 4 5 34 29.20 0.79
## 2 4 5 6 79 86.03 0.57
##
## Response: (0,0,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 4 15 21.38 1.90
## 2 1 5 6 16 19.30 0.56
##
## Response: (1,0,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 4 46 39.35 1.12
## 2 4 5 6 8 10.77 0.71
##
## Response: (0,1,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 3 4 5 43 51.00 1.25
## 2 2 4 5 23 26.52 0.47
##
## Response: (1,1,1)
## Item i Item j Item k Obs Exp (O-E)^2/E
## 1 1 2 4 23 26.48 0.46
## 2 1 4 6 18 16.32 0.17Aside from improving the fit, the model also offers an intuitive interpretation. Estimates of the result are under additive (multiple) parameterization.
Also, if we multiply the value of the second factor (z2) and the interaction term (z1:z2) to rotate the factor, we can see that z1 has high factor loading from item 3-6 (formal factor), whereas z2 has high factor loading on item 1 (informal factor), and item 2 has high loading on both z1 and z2 (general factor).
As for interaction, majority of the items have negative sign, which indicates negative relationship between factors. Focusing on single factor may reduce the chance of using anotehr factor.
We can also compute factor scores for most combinations of response pattern with ltm::factor.scores
factor.scores(fit8, resp.patterns = NULL, method = "MI")##
## Call:
## ltm(formula = WIRS ~ z1 * z2)
##
## Scoring Method: Multiple Imputation
## # Imputations: 5
##
## Factor-Scores for observed response patterns:
## Item 1 Item 2 Item 3 Item 4 Item 5 Item 6 Obs Exp z1 se.z1 z2
## 1 0 0 0 0 0 0 132 135.280 -0.967 0.654 0.423
## 2 0 0 0 0 0 1 5 2.638 -0.054 0.688 -0.299
## 3 0 0 0 0 1 0 17 11.274 0.006 0.639 -0.210
## 4 0 0 0 1 0 0 13 13.770 -0.556 0.668 0.455
## 5 0 0 1 0 0 0 12 8.129 0.066 0.650 -0.161
## 6 0 0 1 0 0 1 1 0.718 0.606 0.599 -0.269
## 7 0 0 1 0 1 0 2 3.830 0.640 0.575 -0.223
## 8 0 0 1 0 1 1 2 0.889 1.110 0.581 -0.227
## 9 0 0 1 1 0 0 1 1.716 0.430 0.630 -0.129
## 10 0 0 1 1 1 0 1 1.584 0.963 0.587 -0.175
## 11 0 1 0 0 0 0 172 167.744 -0.267 0.555 0.502
## 12 0 1 0 0 0 1 4 6.983 0.320 0.533 0.217
## 13 0 1 0 0 1 0 38 42.974 0.352 0.492 0.361
## 14 0 1 0 0 1 1 11 6.877 0.836 0.491 0.304
## 15 0 1 0 1 0 0 45 42.406 0.052 0.501 0.778
## 16 0 1 0 1 0 1 1 2.658 0.597 0.494 0.434
## 17 0 1 0 1 1 0 22 25.657 0.604 0.447 0.651
## 18 0 1 0 1 1 1 8 7.131 1.090 0.478 0.578
## 19 0 1 1 0 0 0 21 26.141 0.317 0.503 0.385
## 20 0 1 1 0 0 1 2 3.042 0.815 0.510 0.250
## 21 0 1 1 0 1 0 28 27.509 0.828 0.477 0.371
## 22 0 1 1 0 1 1 11 11.960 1.326 0.518 0.400
## 23 0 1 1 1 0 0 13 12.734 0.551 0.469 0.635
## 24 0 1 1 1 0 1 3 2.545 1.052 0.511 0.447
## 25 0 1 1 1 1 0 35 31.049 1.043 0.466 0.655
## 26 0 1 1 1 1 1 30 30.698 1.599 0.539 0.737
## 27 1 0 0 0 0 0 65 70.989 -0.160 0.731 -0.797
## 28 1 0 0 0 0 1 11 11.166 0.098 0.748 -1.137
## 29 1 0 0 0 1 0 34 36.303 0.246 0.704 -0.924
## 30 1 0 0 0 1 1 10 9.346 0.339 1.198 -1.165
## 31 1 0 0 1 0 0 20 9.346 0.244 0.700 -0.705
## 32 1 0 0 1 0 1 1 1.743 0.621 0.716 -0.811
## 33 1 0 0 1 1 0 5 7.187 0.695 0.678 -0.712
## 34 1 0 0 1 1 1 2 2.081 1.095 0.691 -0.668
## 35 1 0 1 0 0 0 24 24.990 0.506 0.692 -0.894
## 36 1 0 1 0 0 1 8 6.235 0.803 0.768 -0.941
## 37 1 0 1 0 1 0 31 24.560 0.916 0.675 -0.796
## 38 1 0 1 0 1 1 4 8.343 1.330 0.685 -0.692
## 39 1 0 1 1 0 0 2 5.343 0.896 0.672 -0.751
## 40 1 0 1 1 0 1 1 1.570 1.276 0.681 -0.717
## 41 1 0 1 1 1 0 8 8.246 1.333 0.656 -0.653
## 42 1 0 1 1 1 1 7 3.831 1.781 0.679 -0.567
## 43 1 1 0 0 0 0 55 56.423 -0.358 0.770 -0.466
## 44 1 1 0 0 0 1 10 6.477 0.004 0.827 -0.681
## 45 1 1 0 0 1 0 26 20.048 0.300 0.692 -0.383
## 46 1 1 0 0 1 1 3 5.179 0.820 0.612 -0.311
## 47 1 1 0 1 0 0 5 7.083 0.160 0.630 -0.254
## 48 1 1 0 1 1 0 4 4.584 0.717 0.571 -0.236
## 49 1 1 0 1 1 1 2 1.378 1.192 0.576 -0.210
## 50 1 1 1 0 0 0 13 10.508 0.413 0.623 -0.353
## 51 1 1 1 0 0 1 1 1.941 0.895 0.597 -0.326
## 52 1 1 1 0 1 0 6 9.258 0.947 0.579 -0.279
## 53 1 1 1 0 1 1 5 3.964 1.437 0.594 -0.226
## 54 1 1 1 1 0 0 3 2.565 0.778 0.574 -0.254
## 55 1 1 1 1 0 1 2 0.662 1.250 0.588 -0.245
## 56 1 1 1 1 1 0 4 5.156 1.288 0.575 -0.210
## 57 1 1 1 1 1 1 3 4.170 1.817 0.622 -0.163
## se.z2
## 1 0.844
## 2 0.504
## 3 0.483
## 4 0.955
## 5 0.481
## 6 0.371
## 7 0.369
## 8 0.320
## 9 0.423
## 10 0.341
## 11 0.724
## 12 0.535
## 13 0.603
## 14 0.549
## 15 0.789
## 16 0.636
## 17 0.746
## 18 0.712
## 19 0.622
## 20 0.506
## 21 0.582
## 22 0.581
## 23 0.762
## 24 0.637
## 25 0.755
## 26 0.756
## 27 0.561
## 28 0.696
## 29 0.580
## 30 1.180
## 31 0.477
## 32 0.504
## 33 0.446
## 34 0.406
## 35 0.541
## 36 0.605
## 37 0.473
## 38 0.399
## 39 0.459
## 40 0.414
## 41 0.374
## 42 0.306
## 43 0.614
## 44 0.753
## 45 0.445
## 46 0.336
## 47 0.473
## 48 0.344
## 49 0.293
## 50 0.398
## 51 0.323
## 52 0.311
## 53 0.274
## 54 0.347
## 55 0.289
## 56 0.289
## 57 0.265Reference:
Rizopoulos, D. (2006). Ltm: An R package for latent variable modeling and item response theory analyses. Journal of Statistical Software, 17(5). https://doi.org/10.18637/jss.v017.i05
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 ltm_1.1-1 polycor_0.7-10
## [5] msm_1.6.8 MASS_7.3-51.6
##
## 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 Deriv_4.0.1
## [7] rlang_0.4.6 dcurver_0.9.1 stringr_1.4.0
## [10] highr_0.8 tools_4.0.0 parallel_4.0.0
## [13] grid_4.0.0 nlme_3.1-148 vegan_2.5-6
## [16] mgcv_1.8-31 xfun_0.14 htmltools_0.4.0
## [19] yaml_2.2.1 survival_3.1-12 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.0 expm_0.999-4 mvtnorm_1.1-0