# How to measure consistency of teacher grading?

I would like to determine how consistent teacher grades are for a given item, both for that teacher and across all teachers. E.g., does teacher A consistently score item 1 as a 1,2,3, or 4; and do all teachers generally score that item the in a similar manner. I can use the items difficulty (based on additional information I have on an auto-validation part of that item) as a control variable.

What would be the best method to assess this?

My current attempt is a cross-classified model in which I judge consistency by the proportion of variance attributed to each cluster in order to determine what accounts most for the variability in scores. Here is a really naive example:

library(tidyverse)
library(lme4)
library(ordinal)

set.seed(1983)
teachers <- LETTERS[1:5]
students <- letters[1:10]
items <- 1:5

data <- expand.grid(teachers = teachers, students = students, items = items) %>%
rowwise() %>%
# naive method to make scores vary very little within teachers
mutate(scores = case_when(
teachers == "A" ~ sample(c(1:2), 1),
teachers == "B" ~ sample(c(1:2),1),
teachers == "C" ~ sample(c(2:3), 1),
teachers == "D" ~ sample(c(2:3), 1),
teachers == "E" ~ sample(c(3:4), 1)
)) %>%
mutate(scores_ordinal = as.factor(scores))

lmer_model <- lmer(scores ~ 1 + (1 | students ) + (1 | items) + (1 | teachers), data = data)
summary(lmer_model)

ordinal_model <- clmm(scores_ordinal ~ 1 + (1 | students ) + (1 | items) + (1 | teachers), data = data, link = "logit")
summary(ordinal_model)


The score can be 1-4, so the dv here is actually ordinal, and I believe the ordinal model would fit best. Here are the example results:

Cumulative Link Mixed Model fitted with the Laplace approximation

formula: scores_ordinal ~ 1 + (1 | students) + (1 | items) + (1 | teachers)
data:    data

logit flexible  250  -255.05 522.10 347(1391) 8.72e-05 3.0e+02

Random effects:
Groups   Name        Variance  Std.Dev.
students (Intercept) 1.072e-07 0.0003275
teachers (Intercept) 6.037e+00 2.4569922
items    (Intercept) 8.455e-04 0.0290768
Number of groups:  students 10,  teachers 5,  items 5

No Coefficients

Threshold coefficients:
Estimate Std. Error z value
2|1  -1.5531     1.1358  -1.367
1|3  -0.2054     1.1351  -0.181
3|4   5.0787     1.4778   3.437


If I have approached this correctly, is this output, specifically the Random effects std.dev telling me most of the variance is attributed to teachers, meaning the scores vary largely between (or within ?) teachers, a small amount for items, and very little for individual students?

If this is not the best approach, what is a viable alternative?

• It looks to me like you could use the R package rptR : cran.r-project.org/web/packages/rptR/rptR.pdf Commented Nov 9, 2023 at 21:12
• @CaroZ thanks. I will look into this. At first glance, it doesn't look like it works with ordinal data. Does the R statistic in the output represent the ratio of variance? Commented Nov 9, 2023 at 23:51

Your strategy of employing a cross-classified mixed model appears adequate for your goal of analysing teacher grading consistency for specific items, both within and between instructors. You are effectively dividing the variance in scores in this model into components related to teachers, items, and students.

• Random Effects: The standard deviations in your model output's random effects section tell you about the variability in scores attributable to each group (teachers, items, students). A higher standard deviation indicates more variation.

• Teacher variation: The substantial standard deviation for instructors indicates that there is significant variety in how teachers grade, either within each teacher or among teachers.

• Variability of Items: A lower standard deviation for items suggests less variability in scores owing to the item itself, implying that item difficulty is largely consistent between teachers.

• Variability of Students: The relatively tiny standard deviation for students shows that individual student variances have little impact on score variability.