All Code & Exploration were inspired by :

Paek, I., & Cole, K. (2019). Using R for Item Response Theory Model Applications (1st ed.). Routledge. https://doi.org/10.4324/9781351008167

setwd("/Users/isaiahmireles/Desktop/Misconceptions")
cohorts <- read.csv("cohorts.csv")

1 Student Counts (ct) per term & (term, exam) :

library(tidyverse)
cohorts |> group_by(term) |> summarize(ct = n())
cohorts |> group_by(term, exam) |> summarize(ct = n())
lst <- cohorts |>
  group_by(term, exam) |>
  group_split()

2 S24 Midterm

S24M <- lst[[4]]
# only the Qs
S24MQ <- S24M |> select(matches("^Q\\d+$"))
dim(S24MQ)
## [1] 149  34
S24MQ |> head()
library(eRm)
Rasch_CML <- RM(S24MQ)
help(RM)

  • eRm package uses : conditional maximum likelihood (CML) estimation. Therefore it does not require the specification of the latent trait distribution ( \(\theta\) ).

    • Above I fit a Rasch Model

      • Here we assume \(\sum_{i=1}^{34}b_i=0\) ; meaning, the sum of item difficulty is 0.

        • Typical item Difficulty

          • \(b_i\approx0\) – about avr. difficulty

        • Easy Item Difficulty

          • Making \(b_i < 0\) – below avr. difficulty

        • Hard Item Difficulty

          • Making \(b_i > 0\) – above avr. difficulty

2.1 Inspect Model Object :

Rasch_CML |> attributes()
## $names
##  [1] "X"           "X01"         "model"       "loglik"      "npar"       
##  [6] "iter"        "convergence" "etapar"      "se.eta"      "hessian"    
## [11] "betapar"     "se.beta"     "W"           "call"       
## 
## $class
## [1] "dRm" "Rm"  "eRm"

3 Did our model converge normally?

Rasch_CML$conv
## [1] 1
  • Yes! Finished within the default stopping criteria without being flagged by particular problems.

4 Item Difficulty & Standard Error

item_difficulty <- cbind(-Rasch_CML$betapar, Rasch_CML$se.beta)
colnames(item_difficulty) <- c("est", "SE")
item_difficulty <- item_difficulty |> as.data.frame(); item_difficulty

5 Difficulty Confidence Interval

difficulty_CI <- cbind(-confint(Rasch_CML)[,2], -confint(Rasch_CML)[,1])
colnames(difficulty_CI) <- c("Low", "High")
difficulty_CI <- difficulty_CI |> as.data.frame(); difficulty_CI
difficulty <- cbind(difficulty_CI, est=item_difficulty[,1], item = names(Rasch_CML$betapar))

# summary(Rasch_CML)
# better method^
difficulty |>
  ggplot(aes(
    x = est,
    y = reorder(item, est),
    xmin = Low,
    xmax = High
  )) +
  geom_pointrange() +
  geom_vline(
    xintercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Item Difficulty CI-95%",
    x = "Item Difficulty",
    y = "Item"
  ) +
  theme_minimal()

# set margin sizes
par(
  mfrow = c(2, 1),
  mar = c(3, 4, 2, 1)
)

hist(
  item_difficulty$est,
  main = "Item Difficulty Dist"
)

boxplot(
  item_difficulty$est,
  horizontal = TRUE
)

# reset 
par(mfrow = c(1, 1))

  • As we can see difficulty is normally distributed

    • Most problems are typical difficulty with few being extremely easy, few being extremely hard.

6 What percent of people got each item correct?

percent_correct <- S24MQ |>
  summarize(
    across(
      everything(),
      ~ mean(.x, na.rm = TRUE) * 100
    )
  )
percent_correct

7 Do these measures correlate?

percent_correct_long <- percent_correct |>
  pivot_longer(
    cols = everything(),
    names_to = "item",
    values_to = "percent_correct"
  )

difficulty_percent <- difficulty |>
  mutate(
    item = str_remove(item, "^beta ")
  ) |>
  left_join(
    percent_correct_long,
    by = "item"
  )
difficulty_percent |>
  ggplot(aes(
    x = est,
    y = percent_correct
  )) +
  geom_point() +
  geom_text(
    aes(label = item),
    nudge_y = 2,
    size = 3
  ) +
  labs(
    title = "Item Difficulty vs. Percent Correct",
    x = "Item Difficulty",
    y = "Percent Correct"
  ) +
  theme_minimal()

  • clear, non-linear relationship
cor(difficulty_percent$est, difficulty_percent$percent_correct)^2
## [1] 0.9154302

  • about 92% of variation can be described via regressing est and percent_correct

    • Therefore despite non-linear trend, we may adequately summarize relationship linearly

8 Item Characteristic Curve (ICC)

# include data pts (emperical ICC)
plotICC(Rasch_CML, empICC = list("raw"))

  • usually for this we go per-item and analyze each. Above is just a demonstration.

9 Person Specific latent trait

Person.Rasch_CML <- person.parameter(Rasch_CML)
Person.Rasch_CML 
## 
## Person Parameters:
## 
##  Raw Score     Estimate Std.Error
##         13 -0.673932618 0.4142515
##         15 -0.336827848 0.4077035
##         16 -0.171289853 0.4062696
##         17 -0.006417529 0.4059753
##         18  0.158657094 0.4067802
##         19  0.324814862 0.4087111
##         20  0.493033349 0.4117921
##         21  0.664332670 0.4161253
##         22  0.839784956 0.4218264
##         23  1.020667015 0.4290985
##         24  1.208586752 0.4382412
##         25  1.405523315 0.4496752
##         26  1.613988463 0.4640264
##         27  1.837478272 0.4822548
##         28  2.081086749 0.5059079
##         29  2.352518700 0.5375958
##         30  2.664566213 0.5821649
##         31  3.041136296 0.6498887
##         32  3.535665240 0.7677078
##         33  4.318303521 1.0440940
##         34  5.170735537        NA
?person.parameter

cohorts |> 
  filter(exam=="Midterm", term=="S24") |> 
  select(student_id,Total.Score) |> 
  arrange(Total.Score) |> 
  head(3)
cohorts |> 
  filter(exam=="Midterm", term=="S24") |> 
  select(student_id,Total.Score) |> 
  arrange(desc(Total.Score)) |> 
  head(3)

  • Here we see for a given raw score the est. ability ( \(\theta\) ) score and its estimated Standard error.

  • Notice our estimates are only from the lowest to highest possible score obtained by students in data (S24M)

summary(S24M$Total.Score)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.00   23.00   27.00   26.33   30.00   34.00
hist(S24M$Total.Score)

  • Bottom 25% of class (Quartile 1) have an ability of 1.02 or below

  • Mean Equates to an ability of 1.61

  • Median Equates to 1.84

  • Score of 20 equates to 0.49 ability

    • As we can see our ability does match our emperical total score performance. based on the assumption of our model, it suggests that we have a very capable sample (ie this was an easy exam for them)*
  • Notice ability rises increasingly fast near the maximum possible score

10 Person Specific Ability

Person_ability <- coef(Person.Rasch_CML) |> as.data.frame()
colnames(Person_ability) <- c("ability")
Person_ability |> arrange(ability)
coef(Person.Rasch_CML) |> hist(xlab = "Ability", main="")

summary(Person_ability$ability)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.6739  1.0207  1.8375  1.8939  2.6646  5.1707
summary(S24M$Total.Score)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.00   23.00   27.00   26.33   30.00   34.00

  • Estimated Person specific ability

    • Based on Raw Score – same score means same ability

    • Here we see the typical ability to actually be centered around 2 (1.89) – that would equate to about a raw score of 27, which from earlier makes more sense on the true typical behavior of students

11 Raw Score vs. Ability

mean_raw <- mean(S24M$Total.Score, na.rm = TRUE)

mean_ability <- mean(
  Person_ability$ability,
  na.rm = TRUE
)

plot(
  Person.Rasch_CML,
  main = "Raw Score vs. Person Ability"
)

abline(
  v = mean_raw,
  lty = 2
)

abline(
  h = mean_ability,
  lty = 2
)

  • average Score, Ability Cross is provided alongside the equated Score, Ability Graph

  • Graph extends from worse to best score obtained in data

12 Person item Map (Wright Map)

plotPImap(Rasch_CML, sorted=TRUE)

  • notice person ability are highly varied – maybe there is some multiple modality

  • Dots are item difficulty

    • Difficulty ranges from -3 to 3 while person ability appears to range from -1 to 4

13 Model Data Fit : GOF measures

  • GOF : Goodness of Fit

    • It is suggested we use various techniques not one to evaluate model(s)

13.1 Andersen’s LR test

# Andersen’s LR test
alrt <- LRtest(Rasch_CML)
## Warning in LRtest.Rm(Rasch_CML): 
## The following items were excluded due to inappropriate response patterns within
## subgroups:
## Q1 Q2 Q7 Q9 Q20 Q25
## 
## Full and subgroup models are estimated without these items!
alrt
## 
## Andersen LR-test: 
## LR-value: 49.913 
## Chi-square df: 27 
## p-value:  0.005
?LRtest.Rm

  • This test basically asks : After splitting by median raw score values, are the item-difficulty parameters essentially the same for the lower-scoring and higher-scoring groups?

  • Therefore, you reject the null hypothesis of parameter invariance. Item difficulties are not sufficiently similar across the two ability groups, so there is evidence that the simple Rasch model does NOT fit the data well globally

    • Notice however many questions were taken away

13.2 Martin-Löf (MLoef) LR test

MLoef(Rasch_CML, splitcr = "median")
## 
## Martin-Loef-Test (split criterion: median)
## LR-value: 72.539 
## Chi-square df: 288 
## p-value: 1
?MLoef

  • splits data based on subgroups of items and conducts a test of a homogeneous item set

  • according to this test, the items are consistent with measuring a single latent trait

In other words :

  • The exam may still be measuring one main ability.

  • But some items may function differently across student subgroups.

13.3 Wald test

alrt
## 
## Andersen LR-test: 
## LR-value: 49.913 
## Chi-square df: 27 
## p-value:  0.005
  • So because of our results from earlier^, its a good idea to take a look at each particular question
Waldtest(Rasch_CML)
## Warning in Waldtest.Rm(Rasch_CML): 
## The following items were excluded due to inappropriate response patterns within
## subgroups:
## Q1 Q2 Q7 Q9 Q20 Q25
## 
## Subgroup models are estimated without these items!
## 
## Wald test on item level (z-values):
## 
##          z-statistic p-value
## beta Q3       -0.273   0.785
## beta Q4        1.658   0.097
## beta Q5        0.799   0.424
## beta Q6       -0.766   0.444
## beta Q8       -0.126   0.900
## beta Q10      -1.141   0.254
## beta Q11       0.805   0.421
## beta Q12      -0.567   0.571
## beta Q13       0.247   0.805
## beta Q14       1.141   0.254
## beta Q15       0.778   0.437
## beta Q16       0.389   0.697
## beta Q17      -0.143   0.886
## beta Q18      -0.839   0.402
## beta Q19      -1.191   0.234
## beta Q21       0.576   0.564
## beta Q22       3.589   0.000
## beta Q23      -0.352   0.725
## beta Q24       4.259   0.000
## beta Q26      -0.135   0.892
## beta Q27      -1.114   0.265
## beta Q28      -1.114   0.265
## beta Q29      -0.396   0.692
## beta Q30       0.112   0.911
## beta Q31      -1.192   0.233
## beta Q32       0.713   0.476
## beta Q33      -1.737   0.082
## beta Q34       0.087   0.931
?Waldtest

  • Q22, 24 have unusual z scores

    • suggesting these as the possible location of group invariance

  • It is suggested to use Bonferroni correction because of the possible family-wise Type I error rate

wald <- Waldtest(Rasch_CML)
## Warning in Waldtest.Rm(Rasch_CML): 
## The following items were excluded due to inappropriate response patterns within
## subgroups:
## Q1 Q2 Q7 Q9 Q20 Q25
## 
## Subgroup models are estimated without these items!
wald$coef.table <- cbind(
  wald$coef.table,
  p.adj.bonf = p.adjust(wald$coef.table[, "p-value"],
                        method = "bonferroni")
)

# look at only sign.
wald$coef.table[
  wald$coef.table[, "p.adj.bonf"] < 0.05,
  ,
  drop = FALSE
]
##          z-statistic      p-value   p.adj.bonf
## beta Q22    3.588739 3.322818e-04 0.0093038901
## beta Q24    4.259341 2.050308e-05 0.0005740863
  • notice still Q22, Q24 are significant
plotGOF(alrt, conf=list())

  • Item difficulty estimates from the two groups are plotted with their confidence ellipses

    • Notice Q22, 24 (problematic) questions are the ones not intersecting with the line

13.4 item fit indices*

itemfit(Person.Rasch_CML)
## 
## Itemfit Statistics: 
##       Chisq  df p-value Outfit MSQ Infit MSQ Outfit t Infit t Discrim
## Q1  103.937 147   0.997      0.702     0.988   -0.353   0.069   0.185
## Q2   74.320 147   1.000      0.502     0.943   -0.719  -0.048   0.261
## Q3   93.158 147   1.000      0.629     0.894   -0.711  -0.347   0.354
## Q4  176.393 147   0.049      1.192     1.158    1.807   2.054   0.258
## Q5  144.118 147   0.552      0.974     0.995   -0.210  -0.034   0.356
## Q6  185.745 147   0.017      1.255     1.017    1.936   0.247   0.374
## Q7  113.789 147   0.981      0.769     0.957    0.179   0.150   0.064
## Q8  141.447 147   0.614      0.956     1.012   -0.199   0.161   0.383
## Q9   72.886 147   1.000      0.492     0.908   -0.498  -0.066   0.268
## Q10 121.417 147   0.939      0.820     0.882   -1.581  -1.623   0.528
## Q11 168.023 147   0.113      1.135     0.925    0.431  -0.025   0.192
## Q12 106.208 147   0.995      0.718     0.848   -0.892  -0.961   0.469
## Q13 223.114 147   0.000      1.508     0.902    1.191  -0.399   0.296
## Q14 153.111 147   0.348      1.035     1.013    0.214   0.137   0.263
## Q15 188.145 147   0.012      1.271     1.097    0.898   0.648   0.112
## Q16 122.490 147   0.930      0.828     0.926   -0.271  -0.258   0.312
## Q17 115.728 147   0.973      0.782     1.002   -0.279   0.099   0.180
## Q18 104.386 147   0.997      0.705     0.879   -1.305  -1.017   0.486
## Q19  90.399 147   1.000      0.611     0.890   -1.339  -0.670   0.444
## Q20  43.721 147   1.000      0.295     0.842   -1.139  -0.299   0.394
## Q21 129.776 147   0.843      0.877     0.973   -0.566  -0.233   0.345
## Q22 207.358 147   0.001      1.401     1.225    2.509   2.315   0.096
## Q23 122.152 147   0.933      0.825     0.919   -1.058  -0.909   0.428
## Q24 241.299 147   0.000      1.630     1.396    4.379   4.708   0.013
## Q25  97.628 147   0.999      0.660     0.921   -1.034  -0.421   0.362
## Q26 133.376 147   0.783      0.901     1.037   -0.157   0.258   0.195
## Q27  89.453 147   1.000      0.604     0.831   -1.699  -1.341   0.527
## Q28 116.760 147   0.969      0.789     0.991   -0.792  -0.026   0.314
## Q29 139.757 147   0.652      0.944     0.990   -0.387  -0.107   0.388
## Q30  97.341 147   0.999      0.658     0.945   -0.847  -0.209   0.343
## Q31 103.137 147   0.998      0.697     0.866   -1.746  -1.411   0.503
## Q32 142.973 147   0.578      0.966     1.057   -0.162   0.677   0.296
## Q33 125.845 147   0.896      0.850     0.880   -1.523  -1.679   0.470
## Q34 151.668 147   0.379      1.025     1.009    0.274   0.150   0.366
?itemfit