Davin Díaz.

Ocho análisis DIF en R, paso a paso

Guía de implementación de los análisis de funcionamiento diferencial del reactivo usados en mi tesis de maestría

Davin Díaz

2024-03-01

Escritura

Ocho análisis DIF en R, paso a paso

Cuando empecé la parte analítica de mi tesis, entender qué hacía cada método de detección de funcionamiento diferencial del reactivo me costó bastante menos que conseguir que corrieran. La literatura describe los procedimientos con precisión y los paquetes de R los implementan bien, pero entre una cosa y la otra hay un trecho que casi nadie documenta, el de qué forma exacta tiene que tener la base de datos, qué argumento espera cada función, cómo se automatiza una comparación por pares sin escribirla seis veces y cómo se pasa de la salida cruda a una tabla que alguien pueda leer.

Esta guía intenta recorrer ese trecho, pues sigue los ocho análisis DIF que hice, en el orden en que los hice, con el código que usé para cada uno y una explicación de qué hace cada paso y por qué.

Alcance de este documento

Esto es una guía de implementación, no el estudio. Aquí se explica cómo ejecutar cada procedimiento en R; no se discuten los resultados sustantivos, la selección del modelo ni las decisiones de validez que se siguen de ellos.

Esas discusiones están en la tesis: Díaz García, D. E. (2024). El papel del funcionamiento diferencial del reactivo (DIF) en la validez de exámenes de alto impacto [Tesis de maestría, Universidad Nacional Autónoma de México]. Consultar en TESIUNAM.

Para la comparación entre lo que arrojó cada método y por qué divergen, está el artículo Ocho análisis, ocho respuestas distintas.

Sobre los datos

El estudio original se hizo con microdatos reales de un examen de comprensión lectora para posgrado, que no son públicos ni se distribuyen aquí: contienen información a nivel de sustentante y permiten recuperar la clave del examen.

Para que la guía se pueda seguir, el bloque de la siguiente sección genera un conjunto de datos sintéticos con la misma estructura de columnas, un modelo 2PL de fondo y DIF sembrado deliberadamente en tres reactivos. El resto del código está pensado para correr sobre esos datos; la excepción es DIFRasch(), que necesita una base con los reactivos ya calificados (ver la nota de esa sección). Los números que verás no son los de la tesis.

Paquetes utilizados

Para los análisis de la tesis se utilizaron los siguientes paquetes de R; DIFtree, WrightMap y psychotree se cargan más adelante, en la sección que los usa:

Manejo de los datos

library(data.table)
library(psych)
library(psy)
library(mirt)
library(sjPlot)
library(difR)
library(ltm)
library(lme4)
library(lordif)
library(ShinyItemAnalysis)
library(tidyverse)
library(gt)
library(kableExtra)
library(corrplot)
library(RColorBrewer)
library(eRm)
library(ggdist)

Generación de los datos sintéticos

El estudio original trabajó con tres bases derivadas de los microdatos del examen. Como esos datos no son públicos, el siguiente bloque las reconstruye por simulación, conservando la estructura exacta que espera el resto de la guía:

Las respuestas se generan bajo un modelo de dos parámetros, y el DIF se siembra deliberadamente: a los reactivos 4, 6 y 10 se les suma dificultad para los sustentantes de las áreas 3 y 4. Ese es el efecto que los métodos de las secciones siguientes deberían recuperar, y tener la respuesta correcta de antemano es justamente lo que permite comprobar que cada procedimiento está bien implementado. La siembra imita lo que muestran los datos de la tesis al comparar CFMI y CBQS contra CS y HA: ahí, el 4, el 6 y el 10 son los que señalan tres de los cuatro métodos tradicionales (Tabla 9 de la tesis), y los tres resultaron más difíciles para ciencias sociales y humanidades. Al cruzar todos los métodos, en cambio, los más consistentes fueron el 4, el 5 y el 6. La clave del examen tampoco se reproduce: en la base con las opciones elegidas, el acierto aparece como «X».

set.seed(2024)

## --- Estructura de la muestra -----------------------------------------
## Se replican los tamaños por área del estudio original.
n_area <- c(147, 417, 427, 227)
area   <- rep(1:4, times = n_area)
N      <- length(area)
J      <- 11
items  <- paste0("LR", 1:J)

## --- Parámetros de los reactivos (modelo 2PL) -------------------------
## a = discriminación, b = dificultad.
a <- c(1.02, 0.61, 0.88, 1.34, 0.95, 1.21, 0.77, 1.10, 0.54, 0.92, 1.05)
b <- c(-0.85, 0.42, -0.10, 0.61, -0.35, 0.78, 0.05, -0.52, -1.20, 0.28, 0.15)

## --- DIF sembrado ------------------------------------------------------
items_dif <- c(4, 6, 10)
magnitud  <- 0.6          # en logits

## --- Simulación de respuestas -----------------------------------------
theta <- rnorm(N, mean = 0, sd = 1)
resp  <- matrix(NA_integer_, nrow = N, ncol = J, dimnames = list(NULL, items))

for (j in seq_len(J)) {
  b_j <- rep(b[j], N)
  if (j %in% items_dif) {
    b_j <- b_j + ifelse(area %in% c(3, 4), magnitud, 0)
  }
  p <- 1 / (1 + exp(-a[j] * (theta - b_j)))
  resp[, j] <- rbinom(N, size = 1, prob = p)
}

## --- Covariables administrativas --------------------------------------
posgrado <- sprintf("%02d", sample(1:12, N, replace = TRUE))
genero   <- sample(1:2, N, replace = TRUE)
nivel    <- sample(1:2, N, replace = TRUE, prob = c(0.25, 0.75))
edad     <- round(rnorm(N, mean = 27, sd = 4))

## --- Las tres bases que usa el resto de la guía -----------------------

## 1. Solo los reactivos, calificados
reactivos_lectura <- as.data.frame(resp)

## 2. Reactivos más covariables de agrupación
lectura_dicotomicos <- data.frame(area, posgrado, edad, genero, nivel, resp)

## 3. Base completa con la opción elegida. Se reconstruye a partir del
##    acierto: si el sustentante acertó se le asigna la clave, si no, un
##    distractor al azar. La clave es ficticia a propósito: todos los
##    aciertos se marcan con "X" y los errores con una opción de la A a la D,
##    para que nadie la confunda con la del examen original.
clave    <- rep("X", J)
opciones <- c("A", "B", "C", "D")

crudas <- vapply(seq_len(J), function(j) {
  distractores <- setdiff(opciones, clave[j])
  ifelse(resp[, j] == 1L, clave[j], sample(distractores, N, replace = TRUE))
}, character(N))
colnames(crudas) <- items

lectura <- data.frame(
  area,
  posgrado,
  periodo      = "02",
  ID           = sprintf("%05d", seq_len(N)),
  edad,
  genero,
  nacionalidad = NA_character_,
  nivel,
  lic_proced   = NA_character_,
  lic_UNAM     = NA_character_,
  plantel_UNAM = NA_character_,
  crudas,
  stringsAsFactors = FALSE
)

dim(lectura); dim(lectura_dicotomicos); dim(reactivos_lectura)

Una función para tablas con formato APA

Buena parte del trabajo de reportar resultados es repetitivo. Conviene definir una sola vez la función que da formato a las tablas y reutilizarla en todo el documento:

## Tabla estilo APA
apa <- function(x, title = " ",...) {
  gt(x,...) %>%
    tab_options(
      table.border.top.color = "white",
      heading.title.font.size = px(16),
      column_labels.border.top.width = 3,
      column_labels.border.top.color = "black",
      column_labels.border.bottom.width = 3,
      column_labels.border.bottom.color = "black",
      table_body.border.bottom.color = "black",
      table.border.bottom.color = "white",
      table.width = pct(100),
      table.background.color = "white"
    ) %>%
    cols_align(align="center") %>%
    tab_style(
      style = list(
        cell_borders(
          sides = c("top", "bottom"),
          color = "white",
          weight = px(1)
        ),
        cell_text(
          align="center"
        ),
        cell_fill(color = "white", alpha = NULL)
      ),
      locations = cells_body(
        columns = everything(),
        rows = everything()
      )
    ) %>%
    #title setup
    tab_header(
      title = gt::html("<i>", title, "</i>")
    ) %>%
    opt_align_table_header(align = "left")
}

Análisis descriptivos para examen de comprensión lectora

Distribución de puntajes totales

Para explorar de forma descriptiva la distribución de puntajes de los sustentantes dentro de cada una de las áreas de posgrado, se comenzó por una gráfica de lluvia (raincloud plot). La ventaja de este tipo de gráfica es que superpone tres representaciones de la misma variable: la densidad de la distribución en forma de campana, un diagrama de caja con la mediana y los cuartiles, y la nube de puntos individuales. Así se ve a la vez la forma general de la distribución y la dispersión caso por caso, sin tener que elegir entre una y otra.

Para elaborar este tipo de gráfico a partir del puntaje total de la prueba se siguieron los siguientes pasos:

Los pasos recién descritos corresponden uno a uno con las líneas del siguiente bloque de código utilizado para generar el gráfico.

## Figura 1 - Distribución de puntajes por área
lectura_dicotomicos%>%
  mutate(Total = rowSums(reactivos_lectura))%>%
  drop_na(genero)%>%
  mutate(area = factor(area))%>%
  ggplot(aes(x = area, y = Total, color = area)) +
  ggdist::stat_halfeye(adjust =.8, width = .5, justification = -.5, aes(fill=area)) +
  geom_boxplot(width = .25, outlier.shape = NA) +
  geom_point(size = 1.3, alpha = .4, position = position_jitter(seed =1, width = .4))  +
  coord_cartesian(xlim = c(1.2, NA)) +
  theme_light() +
  labs(title = "Figura X. Distribución de puntajes en examen de lectura y gramática por área de posgrado",
       x = "Área", y = "Puntaje total")+
  scale_fill_manual(values = c("#44BEFE", "#098BD1", "#0D547A", "#092A3C"))+
  scale_color_manual(values = c("#44BEFE", "#098BD1", "#0D547A", "#092A3C"))

Porcentaje de respuestas correctas por reactivo

Una vez observados los resultados a nivel de puntajes totales por área, para visualizar las diferencias a nivel de porcentaje de respuestas correctas por reactivo por área, se llevaron a cabo los siguientes pasos:



# Figura 2 -  Porcentaje de respuestas correctas por área



reactivos_lectura%>%
  cbind(lectura$area)%>%
  mutate(area = factor(`lectura$area`))%>%
  select(-`lectura$area`)%>%
  drop_na()%>%
  group_by(area)%>%
  summarise(across(everything(), mean))%>%
  pivot_longer(!area,
               names_to = "reactivo", values_to = "porcentaje")%>%
  mutate(porcentaje = round(porcentaje*100,1))%>%
  ggplot(aes(x = factor(reactivo, levels = names(reactivos_lectura)), y =porcentaje, fill = area)) + geom_col(position = "dodge") +
  theme_light() +
  labs(title = "Figura X. Porcentaje de respuestas correctas por reactivo en examen de lectura y gramática por área de posgrado",
       x = "Reactivo", y = "Porcentaje de respuestas correctas") +
  scale_fill_manual(values = c("#44BEFE", "#098BD1", "#0D547A", "#092A3C"))

Análisis DIF con distintos métodos

Método de Mantel-Haenszel

Uno de los métodos más utilizados para detectar DIF tanto por su simplicidad como por su practicidad es el método de Mantel-Haenszel. En R, es posible realizar este análisis mediante la función difMH() perteneciente al paquete difR de Magis et al. (2010). Esta función permite la comparación directa entre dos grupos para la detección de DIF. Los principales parámetros que requiere son:

Dados estos requerimientos, es necesario transformar las bases de datos con las que se ha trabajado hasta ahora. Para ejemplificar el uso básico de esta función, se presenta a continuación el análisis entre las áreas 1 y 2. Dentro del bloque de código se realizan las siguientes operaciones:

 lecturaDIF <- lectura%>%
  dplyr::select(area)%>%
  bind_cols(reactivos_lectura)


  areas <- lecturaDIF%>%
    filter(area %in% c(1,2))

  reactivos <- areas%>%
    select(!area)

  MH<-difMH(reactivos, areas$area, focal.name= 2)

  coef <- round(-2.35 * log(MH$alphaMH), 4)

  DIF <- as.vector(symnum(coef, c(-Inf,-1.5,-1, 1, 1.5, Inf), symbols = c("--","-","", "+", "++")))

  tabladif<-paste(coef, DIF, sep = " ")%>%as.data.frame()%>%
    set_names("Area1 Vs Area2")%>%
        bind_cols(reactivos = names(reactivos),. )

Para leer el signo hay que tener presente cuál es el grupo focal. alphaMH es la razón entre los momios de acertar del grupo de referencia y los del focal a igual puntaje total, así que un MH D-DIF negativo indica que el reactivo es más difícil para el grupo focal, aquí el área 2, y uno positivo, que lo es para el de referencia. Es la convención de ETS, y conviene revisar el focal.name antes de interpretar cualquier tabla de este método. En la tesis (p. 66) esta regla aparece invertida: con el primer grupo como referencia, un valor positivo indica que el reactivo es más sencillo para el segundo grupo, no para el primero.

Si bien este procedimiento es relativamente sencillo, hacerlo para cada uno de los pares se vuelve un proceso repetitivo que puede bien ser automatizado. Para ello, en el siguiente bloque de código replicamos el código de la sección anterior y lo encerramos dentro de un bucle del tipo for(). Más concretamente, en el siguiente bloque se realizan los siguientes pasos:




pairs_list <- list(c(1,2), c(1,3), c(1,4), c(2,3), c(2,4), c(3,4))


for(i in 1:length(pairs_list)){

  areas <- lecturaDIF%>%
    filter(area %in% c(pairs_list[[i]][1],pairs_list[[i]][2]))

  reactivos <- areas%>%
    select(!area)

  MH<-difMH(reactivos, areas$area, focal.name= pairs_list[[i]][2])

  coef <- round(-2.35 * log(MH$alphaMH), 4)

  DIF <- as.vector(symnum(coef, c(-Inf,-1.5,-1, 1, 1.5, Inf), symbols = c("--","-","", "+", "++")))

  tabladif<-paste(coef, DIF, sep = " ")%>%as.data.frame()%>%
    set_names(noquote(paste0("Area",pairs_list[[i]][1]," VS Area", pairs_list[[i]][2])))

if(exists("TablaDIF_CL_MH") == T){TablaDIF_CL_MH <- bind_cols(TablaDIF_CL_MH, tabladif)}
if(exists("TablaDIF_CL_MH") == F){TablaDIF_CL_MH <- bind_cols(reactivos = names(reactivos_lectura), tabladif)}

print(names(TablaDIF_CL_MH))

}



TablaDIF_CL_MH%>%
kbl(caption = "Presencia de DIF Área vs Área en Examen de Lectura con Método Mantel-Haenszel")%>%
  kable_classic(full_width = F, html_font = "Cambria")%>%
  add_header_above(c(" " = 1, "FMI y CBQS" = 1, "FMI y CS" = 1, "FMI y HA" = 1,
                     "CBQS y CS" = 1, "CBQS y HA" = 1, "CS y HA" = 1))

Método de Regresión Logística

Para el método de Regresión logística, se utiliza un procedimiento similar al anterior, incluyendo un bucle que llama la función difLogistic del mismo paquete difR. La principal diferencia con el código anterior es que, en lugar de calcular manualmente a partir del valor de Alfa de Mantel-Haenszel, en este caso tomamos directamente el valor p para determinar el nivel de DIF, complementando esta prueba de significancia con el valor de delta que la misma prueba arroja. Otro aspecto distinto al utilizar este método de regresión logística es que tenemos la posibilidad de buscar DIF de tipo uniforme o no uniforme al especificar el parámetro type = "udif" o type = "nudif".





for(i in 1:length(pairs_list)){

  areas <- lecturaDIF%>%
    filter(area %in% c(pairs_list[[i]][1],pairs_list[[i]][2]))

  reactivos <- areas%>%
    select(!area)%>%as.data.frame()

  LR<-difLogistic(reactivos, areas$area, type = "udif", focal.name= pairs_list[[i]][2])

  coef <- round(LR$p.value, 4)
  DIF <- as.vector(symnum(coef, abs(c(0,0.001,0.01,0.05,1)),
                          symbols = c("***","**", "*", "")))
  delta <- round(LR$deltaR2,3)
  tabladif<- paste(delta, DIF, sep = " ")%>%as.data.frame()%>%
    set_names(noquote(paste0("Area",pairs_list[[i]][1]," VS Area", pairs_list[[i]][2])))

if(exists("TablaDIF_CL_LRU") == T){TablaDIF_CL_LRU <- bind_cols(TablaDIF_CL_LRU, tabladif)}
if(exists("TablaDIF_CL_LRU") == F){TablaDIF_CL_LRU <- bind_cols(reactivos = names(reactivos_lectura), tabladif)}

print(names(TablaDIF_CL_LRU))

}


TablaDIF_CL_LRU%>%
kbl(caption = "Análisis DIF con método de Regresión Logística para DIF uniforme área contra área en examen de comprensión lectora")%>%
  kable_classic(full_width = F, html_font = "Cambria")%>%
  add_header_above(c(" " = 1, "FMI y CBQS" = 1, "FMI y CS" = 1, "FMI y HA" = 1,
                     "CBQS y CS" = 1, "CBQS y HA" = 1, "CS y HA" = 1))

Regresión logística recursiva exploratoria

Además del método mostrado en la sección anterior, también es posible utilizar un método exploratorio recursivo de regresión logística gracias a la función DIFtree() perteneciente al paquete del mismo nombre. Para utilizar esta función, necesitamos realizar algunos pasos previos. Dentro del siguiente bloque de código se ejecutan todas las siguientes acciones:

## Regresión logística recursiva exploratoria
## DIFtree
library(DIFtree)

### Preparacion de los datos
covar_lectura <- lectura%>%
  select(area,posgrado, genero, nivel)%>%
  bind_cols(reactivos_lectura)%>%
  drop_na()%>%
  select(area,posgrado, genero, nivel)%>%
  mutate_all(as.numeric)%>%
  as.data.frame()


diftree_data <- lectura%>%
  select(area,posgrado, genero, nivel)%>%
  bind_cols(reactivos_lectura)%>%
  drop_na()%>%
  select(LR1:LR11)%>%
  as.matrix()


### Ejecución de análisis
nudiftree_lectura<- DIFtree(Y = diftree_data, X = covar_lectura,
                            model = "Logistic", type = "udif")

tabla_logtree <- summary(nudiftree_lectura)

tabla_logtree[[2]]%>%
  dplyr::rename(Tipo = type, Particiones = nosplits)%>%
  apa("Tabla 9.  Análisis DIF exploratorio del examen de lectura y gramática con
      método recursivo de regresión logística")

Los números de tabla que aparecen en los títulos del código no siguen la numeración de la tesis, donde este resultado es la Tabla 6 y la Tabla 9 corresponde a la comparación de dos grupos con los cuatro métodos tradicionales. Con type = "udif", como en el bloque, se reproduce la Tabla 6; con type = "dif", que admite también DIF no uniforme, el 6 y el 10 salen como no uniformes y el 5 deja de tener DIF.

Métodos DIF paramétricos

Exploración de los datos con Mapa de Wright

Dentro del marco de la TRI, una visualización útil es el llamado Mapa de Wright, que permite visualizar bajo una misma escala tanto a los reactivos como a los sustentantes. Para elaborarlo, el paquete WrightMap es una herramienta muy útil, aunque requiere de algunos pasos adicionales que involucran otro tipo de paquetes como el paquete mirt para generar el modelo TRI o algunas funciones adicionales básicas de R para la generación del gráfico. Para elaborarlo, necesitamos realizar los siguientes pasos:



## Mapa de Wright con modelo Rasch
library(ShinyItemAnalysis)
library(WrightMap)
library(RColorBrewer)


modeloRasch_CL <- mirt(data=reactivos_lectura,model = 1, itemtype="Rasch",
                       SE=TRUE, verbose=FALSE)

thetaRaschScores <-fscores(modeloRasch_CL)%>%
  cbind(area = lectura$area)%>%
  as_tibble()


thetaRaschCL1 <- thetaRaschScores%>%
  filter(area == 1)%>%
  select(F1)

thetaRaschCL2 <- thetaRaschScores%>%
  filter(area == 2)%>%
  select(F1)

thetaRaschCL3 <- thetaRaschScores%>%
  filter(area == 3)%>%
  select(F1)

thetaRaschCL4 <- thetaRaschScores%>%
  filter(area == 4)%>%
  select(F1)


coef.Rasch_CL <- coef(modeloRasch_CL, IRTpars=TRUE, simplify=TRUE)

b_RaschCL <- coef.Rasch_CL[[1]][,2]

item_color_palette <- matrix(rep(c("#44BEFE","#098BD1","#0D547A","#092A3C")),
                             10, byrow = TRUE, ncol = 4)

split.screen(figs = matrix(c(    0, .10,0,1,
                                 .10, .20,0,1,
                                 .20, .30,0,1,
                                 .30, .40,0,1,
                                 .40,  1, 0,1), ncol = 4, byrow = TRUE))

personHist(thetaRaschCL1,yRange = c(-1.2,1.5), dim.lab.cex = 1,
           dim.names = "CFMI", dim.color = "#44BEFE",
           show.axis.logits = FALSE, axis.persons = "Sustentantes", breaks = 10)

screen(2)
personHist(thetaRaschCL2,yRange = c(-1.2,1.5), dim.lab.cex = 1,
           dim.names = "CBQS", dim.color = "#098BD1",
           show.axis.logits = FALSE, axis.persons = NULL, breaks = 10)

screen(3)
personHist(thetaRaschCL3,yRange = c(-1.2,1.5), dim.lab.cex = 1,
           dim.names = "CS", dim.color = "#0D547A",
           show.axis.logits = FALSE, axis.persons = NULL, breaks = 10)

screen(4)
personHist(thetaRaschCL4, yRange = c(-1.2,1.5), dim.lab.cex = 1,
           dim.names = "HA", dim.color = "#092A3C",
           show.axis.logits = FALSE,axis.persons = NULL, breaks = 10)

screen(5)
itemModern(b_RaschCL, yRange = c(-1.2,1.5), thr.sym.cex = 2, thr.sym.col.bg = "#092A3C")
mtext("Figura 3. Mapa de Wright del examen de lectura y gramática con distribución por área", side = 3, font = 2, line = 1)

close.screen(all.screens = TRUE)

Análisis Raschtree para detección exploratoria de DIF

Los métodos anteriores requieren decidir de antemano qué variable define los grupos. El raschtree invierte el planteamiento: recibe un conjunto de covariables y busca por sí mismo qué particiones de la muestra producen diferencias en los parámetros de los reactivos. En lugar de responder “¿hay DIF por área?”, responde “¿por qué variable, y en qué punto de corte, se parten los datos de forma que los reactivos dejen de comportarse igual?”.

La preparación tiene una particularidad que conviene notar: la función espera que la matriz de respuestas sea una sola columna dentro del data frame, no once columnas sueltas. Por eso el primer paso guarda los reactivos como una matriz en la columna resp, y el segundo elimina las columnas individuales para que no aparezcan también como covariables:


## Figura 5 - Análisis Raschtree

library(psychotree)
lectura_dicotomicos$resp <- as.matrix(lectura_dicotomicos[, 6:16])

## Eliminar variables no requeridas
lecturaDIF1 <- lectura_dicotomicos[, -(6:16)]

## Creación de modelo psychotree lectura
raschtreeCL1 <- psychotree::raschtree(resp ~ area  + genero + nivel, data = lecturaDIF1)

## Gráfica del modelo
plot(raschtreeCL1, title = "Figura 5")

Es la fórmula que produjo la Figura 5 de la tesis. La Figura 6 añade la edad (resp ~ area + edad + genero + nivel), y agregar el posgrado no cambia el árbol.

Análisis DIF con método de diferencia de Logits con un modelo Rasch

Este método ajusta un modelo de Rasch y luego estima los parámetros de dificultad por separado para cada grupo y los compara con la prueba de Wald. El DIF se expresa como la diferencia entre ambos, en logits, lo que tiene la ventaja de ser directamente interpretable: un valor de 0.5 significa que el reactivo es medio logit más difícil para un grupo que para el otro, controlando el nivel de habilidad. En la función de abajo, la diferencia es la dificultad del área de código menor menos la del área de código mayor, de modo que un valor negativo indica que el reactivo es más difícil para la segunda.

La contrapartida es que todo el procedimiento descansa sobre el supuesto de Rasch —discriminación igual en todos los reactivos—, y si ese modelo no ajusta a los datos, la comparación se vuelve poco informativa. Conviene por tanto revisar el ajuste del modelo antes de leer estos resultados; el porqué está desarrollado en el artículo de comparación.

Nota sobre esta función

DIFRasch() espera que base contenga las columnas administrativas de area a plantel_UNAM seguidas de los reactivos ya calificados como acierto o error, porque la selección select(!area:plantel_UNAM) descarta las primeras por posición. La base lectura de esta guía guarda la opción elegida, no el acierto, así que para usar la función hay que construir antes una base equivalente con los reactivos calificados.


library(eRm)
DIFRasch <- function(base, referencia, focal){

  areas <- base%>%
    filter(area %in% c(referencia,focal))

  reactivos <- areas%>%
    select(!area:plantel_UNAM)

  modeloR<-RM(reactivos)
  ## Waldtest() ordena los grupos por su código: coef es la dificultad del área
  ## de código menor menos la del mayor, sin importar cuál se pase como focal
  subgroup_diffs <- Waldtest(modeloR, splitcr = areas$area)
  subgroup_1_diffs <- subgroup_diffs$betapar1
  subgroup_2_diffs <- subgroup_diffs$betapar2
  coef <- round((subgroup_1_diffs*-1)-(subgroup_2_diffs*-1),3)

  DIF <- symnum(coef, c(-Inf,-1,-.5, .5, 1, Inf), symbols = c("--","-","", "+", "++"))

  tabla_dif<-cbind(coef, DIF)
}

Análisis DIF con método de Ji cuadrada de Lord

El método de Ji cuadrada de Lord compara directamente los parámetros estimados de cada reactivo entre el grupo de referencia y el grupo focal, bajo un modelo de TRI que se especifica explícitamente. Esa especificación es lo que lo distingue de los métodos anteriores: al pasar model = "2PL" se pide que la comparación se haga sobre un modelo de dos parámetros, en lugar de heredar los supuestos del modelo de Rasch.

Como el procedimiento se repite para cada par de áreas, conviene encapsularlo en una función en lugar de un bucle. La función recibe la base, el grupo de referencia y el grupo focal, y devuelve el estadístico de cada reactivo acompañado de su nivel de significancia:



## DIF X2 de Lord
## Función personalizada para análisis DIF con X2 de Lord
DIF_LORD <- function(base, referencia, focal){

  areas <- base%>%
    filter(area %in% c(referencia,focal))

  reactivos <- areas%>%
    select(!area)

  Lord<-difLord(reactivos, areas$area, focal.name= focal, model = "2PL")
  coef <- round(Lord$p.value, 4)
  DIF <- as.vector(symnum(coef, abs(c(0,0.001,0.01,0.05,1)),
                          symbols = c("***","**", "*", "")))
  chi <- round(Lord$LordChi,3)
  tabladif<- paste(chi, DIF, sep = " ")
}

Análisis DIF con áreas combinadas

Dividir la muestra en cuatro áreas deja grupos pequeños, y varios de los métodos anteriores pierden potencia —o directamente no convergen— con pocos casos por grupo. Una alternativa es colapsar las cuatro áreas en dos agrupaciones defendibles conceptualmente: por un lado las áreas de ciencias exactas y de la salud, por el otro las de ciencias sociales y humanidades.

Esto no sustituye al análisis por pares, sino que lo complementa: da una partición con suficientes casos en cada lado como para que los métodos basados en TRI trabajen con holgura, a costa de perder la resolución entre áreas específicas.

En la tesis, las funciones de las secciones anteriores se llamaron con (lectura_2grps, 2, 1), es decir, con CS y HA como referencia y CFMI y CBQS como grupo focal. La Tabla 9 de la tesis se calculó además con los 1,218 casos, sin el drop_na() del bloque; con él quedan 1,201 y algunos valores cambian un poco (el MH D-DIF del reactivo 10, por ejemplo, pasa de 1.000 a 0.990). Para usar DIFRasch() con esta base, que solo tiene area y los reactivos, hay que cambiar su selección por select(LR1:LR11).

## DIF con 2 grupos. Comparación de métodos
lectura_2grps <- lecturaDIF%>%
  mutate(area = case_when(area %in% 1:2 ~ 1, area %in% 3:4 ~ 2))%>%
  mutate(area = as.factor(area))%>%
  select(area,LR1:LR11)%>%
  drop_na()

Análisis DIF con metodologías robustas mediante TRI

Los métodos de esta última sección son los más exigentes, y también los que hacen explícito algo que los anteriores dejan implícito: antes de comparar grupos hay que decidir qué modelo representa los datos. El bloque siguiente hace las dos cosas en orden.

Primero, la selección del modelo. Se estiman los modelos de Rasch, 2PL y 3PL, se comparan sus índices de bondad de ajuste mediante el estadístico M2, y se contrasta además el modelo unidimensional elegido contra uno exploratorio de dos dimensiones mediante una prueba de razón de verosimilitud. Sobre el modelo resultante, que en la tesis fue el 2PL unidimensional, se revisa el supuesto de independencia local con los residuales Q3; en la tesis, ningún par de reactivos pasó de .2 en valor absoluto. El Anexo de la tesis calcula el Q3 sobre el modelo de dos dimensiones, pero la Figura 7 corresponde al 2PL unidimensional, que es el que usa este bloque. Además, el modelo de dos dimensiones no convergió, así que la comparación de la Tabla 11 conviene leerla con reserva; antes de leer esa prueba con tus propios datos, comprueba que CL_mirt2 haya convergido (extract.mirt(CL_mirt2, "converged")).

Después, el análisis DIF propiamente dicho. Se estima un modelo de dos parámetros por grupo con multipleGroup() y se contrasta la invarianza de los parámetros a1 y d con la función DIF(). El procedimiento se ejecuta dos veces: la primera de forma exploratoria, para identificar los reactivos que se comportan de manera más estable y que pueden servir como anclas; la segunda fijando esos reactivos ancla mediante el argumento invariance, de modo que la comparación entre grupos se haga sobre una escala común. Sin ese segundo paso, los parámetros de cada grupo están en escalas que no son directamente comparables y el resultado es difícil de interpretar. En la tesis, el primer paso con las cuatro áreas dejó como único reactivo libre de DIF al 2, que es el que fija el código (LR2); en la comparación de CFMI y CBQS contra CS y HA, las anclas fueron el 1, el 2 y el 8.


## Análisis DIF TRI robustos

## Selección de modelo TRI

## Generación de cada modelo
CL_irt1<-mirt(reactivos_lectura, model = 1, itemtype = "Rasch")
CL_irt2<-mirt(reactivos_lectura, model = 1, itemtype = "2PL")
CL_irt3<-mirt(reactivos_lectura, model = 1, itemtype = "3PL")
CL_mirt2<-mirt(reactivos_lectura, model = 2, itemtype = "2PL") # Modelo bidimensional




### Generación de índices de bondad de ajuste

fit_CL_irt1 <- M2(CL_irt1, na.rm = TRUE)
fit_CL_irt2 <- M2(CL_irt2, na.rm = TRUE)
fit_CL_irt3 <- M2(CL_irt3, na.rm = TRUE)


nombres_modelos1 <- c("Modelo Rasch", "Modelo 2PL", "Modelo 3PL")

## Tabla final
rbind(fit_CL_irt1, fit_CL_irt2, fit_CL_irt3)%>%
  select(M2:RMSEA, CFI)%>%
  cbind(nombres_modelos1, .)%>%
  gt(rowname_col = "nombres_modelos1")%>%
  tab_header(title = "Tabla 3. Comparación de modelos TRI unidimensionales en examen de lectura y gramática")%>%
  tab_stubhead(label = "Modelos")%>%
  fmt_number(columns = c(2,4:6),
             decimals = 3)

## Análisis de Razón de Verosimilitud entre modelo unidimensional y bidimensional
anova_irt1_vs_mirt2 <-anova(CL_irt2, CL_mirt2)

## Generación de la tabla
Modelos <- c("Modelo 1", "Modelo 2")
anova_irt1_vs_mirt2%>%
  select(AIC, BIC, logLik, X2, df, p)%>%
  cbind(Modelos, .)%>%
  gt()%>%
  tab_header(title = "Tabla 4. Comparativa de Modelo unidimensional y multidimensional")%>%  fmt_number(columns = 2:5,
             decimals = 2)%>%
  tab_footnote(footnote = "Modelo unidimensional",
               locations = cells_body(columns = Modelos, rows = 1))%>%
  tab_footnote(footnote = "Modelo exploratorio de dos dimensiones",
               locations = cells_body(columns = Modelos, rows = 2))


## Evaluación de supuesto de independencia local con modelo elegido
residuals(CL_irt2, type = "LDG2")

residuals(CL_irt2, type = "Q3")%>%
  corrplot(method = "number", type="lower", col=brewer.pal(n=11, name="RdBu"),tl.srt=90)


## Generación de modelo de dos parámetros específicos para cada grupo usando multipleGroup del paquete mirt
modelo2PL_CL_A <- multipleGroup(data=reactivos_lectura,model = 1, group = factor(lectura$area),
                                itemtype="2PL", SE=TRUE)

#Análisis DIF a partir del modelo de dos parámetros anterior
dif_CL1 <-DIF(modelo2PL_CL_A, which.par = c("a1", "d"))
dif_CL1%>%arrange(X2)

## A partir del análisis anterior, se identifican reactivos ancla y se vuelve a ejecutar el análisis fijando los parámetros en dichos reactivos

modelo2PL_CL_ANCH <- multipleGroup(data=reactivos_lectura,model = 1, group = factor(lectura$area),
                                   itemtype="2PL", SE=TRUE, invariance = c(c("LR2"),
                                                                           'free_means', 'free_var'))

#Análisis DIF ya con los reactivos ancla

dif_CL2<-DIF(modelo2PL_CL_ANCH, c("a1", "d"), items2test = c(1,3:11), maxiter = 1000)


nombres_reactivos2 <- rownames(dif_CL2)
## generación de la tabla final

dif_CL2%>%
  cbind(nombres_reactivos2,.)%>%
  rename(Reactivos =nombres_reactivos2)%>%
  gt()%>%
  tab_header(title = "Análisis DIF con método de Razón de Verosimilitud en examen de lectura y gramática")%>%
  fmt_number(columns = c(3:7,9),
             decimals = 3)

Comprobación con los datos sintéticos

Como el DIF de esta guía se sembró a propósito en los reactivos 4, 6 y 10, se puede verificar qué tanto lo recupera cada método. El bloque siguiente lo hace con Mantel-Haenszel y con la regresión logística.

## Recuperación del DIF sembrado, comparando areas 1 y 2 (sin DIF)
## contra areas 3 y 4 (con DIF)
grupos_2 <- ifelse(area %in% c(1, 2), 1, 2)

MH_check <- difMH(reactivos_lectura, grupos_2, focal.name = 2)
LR_check <- difLogistic(as.data.frame(reactivos_lectura), grupos_2,
                        type = "udif", focal.name = 2)

data.frame(
  Reactivo  = items,
  Sembrado  = ifelse(seq_len(J) %in% items_dif, "si", "no"),
  MH_p      = round(MH_check$p.value, 4),
  RegLog_p  = round(LR_check$p.value, 4)
) %>%
  kbl(caption = "Recuperación del DIF sembrado con datos sintéticos") %>%
  kable_classic(full_width = FALSE, html_font = "Cambria")

Para leer completo

Esta guía cubre el cómo. La discusión de qué significan los resultados, por qué los métodos divergen y qué se puede concluir de ello está repartida entre dos lugares:

Díaz García, D. E. (2024). El papel del funcionamiento diferencial del reactivo (DIF) en la validez de exámenes de alto impacto [Tesis de maestría, Universidad Nacional Autónoma de México]. Programa de Maestría y Doctorado en Psicología, Residencia en Evaluación Educativa, Facultad de Psicología. Consultar en TESIUNAM

Nota sobre la elaboración

Los análisis y el código de esta guía provienen de mi tesis de maestría y son de mi autoría. La redacción de esta versión y el bloque de generación de datos sintéticos se elaboraron con asistencia de IA, y fueron revisados y verificados por mí. La responsabilidad por el contenido es enteramente mía.

Referencias