Chapter 3 TRI 2PL
3.1 Análise de DIF
Faremos uma análise de DIF com o método de regressão logística para verificar a invariância da medida entre as faculdades.
Primeiro, mostrarei a ideia por trás da análise tomando um item como ilustração. Em seguida, apresentarei o resultado.
A primeira coisa a se fazer é preparar o banco. É preciso filtrar o banco pelo ano de aplicação, e selecionar as variáveis ID, grupo e os itens. A análise se deu por ano, por considerar que existem pessoas que não compareceram nos dois anos e principalmente porque eu quis evitar que problemas em um ano contaminassem a análise do outro ano. Como são 120 itens por ano, considero que esse número já é suficiente para dar uma boa estimativa do theta da pessoa. Por isso, podemos analisar ano a ano. E de quebra, é mais simples analisar assim do que fazer os cruzamentos.
# selecionar somente o ano
banco.dif <- select (banco_total, ID2, grupo, starts_with('Q2019')) %>%
na.omit()Agora precisamos calibrar os itens para calcular o theta das pessoas. Sim, precisa calibrar de novo porque será utilizado o método de múltiplos grupos, em que cada grupo será uma faculdade.
# calibrar
calib1 <- multipleGroup(
banco.dif[,-c(1:2)],
1,
group = banco.dif$grupo,
itemtype = '2PL',
TOL = .01,
invariance = c('free_means', 'free_var', names(banco.dif)[-c(1:2)])
)##
Iteration: 1, Log-Lik: -86433.803, Max-Change: 1.03349
Iteration: 2, Log-Lik: -84832.437, Max-Change: 0.21569
Iteration: 3, Log-Lik: -84813.951, Max-Change: 0.03188
Iteration: 4, Log-Lik: -84813.433, Max-Change: 0.01167
Iteration: 5, Log-Lik: -84813.275, Max-Change: 0.00581# calcular o theta
theta1 <- fscores(calib1)A análise consiste em comparar três modelos:
M1: item ~ intercepto + theta
M2: item ~ intercepto + theta + grupo
M3: item ~ intercepto + theta + grupo + theta*grupo
Se houver diferença signifativa na inserção de algum parâmetro, ou se o modelo explica bem mais (R2) os dados, então consideramos que existe DIF. O modelo de regressão logística é criado por meio da função lrm.
# M1 do item 1
m1 <- lrm(
Q2019_001 ~ theta1,
data = banco.dif
)
m1## Logistic Regression Model
##
## lrm(formula = Q2019_001 ~ theta1, data = banco.dif)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 1198 LR chi2 0.14 R2 0.000 C 0.518
## 0 648 d.f. 1 g 0.025 Dxy 0.036
## 1 550 Pr(> chi2) 0.7048 gr 1.025 gamma 0.037
## max |deriv| 5e-14 gp 0.006 tau-a 0.018
## Brier 0.248
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -0.1632 0.0580 -2.81 0.0049
## F1 -0.0225 0.0594 -0.38 0.7049
## # M2 do item 1
m2 <- lrm(
Q2019_001 ~ theta1 + grupo,
data = banco.dif
)
m2## Logistic Regression Model
##
## lrm(formula = Q2019_001 ~ theta1 + grupo, data = banco.dif)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 1198 LR chi2 0.19 R2 0.000 C 0.514
## 0 648 d.f. 2 g 0.029 Dxy 0.028
## 1 550 Pr(> chi2) 0.9087 gr 1.029 gamma 0.029
## max |deriv| 2e-14 gp 0.007 tau-a 0.014
## Brier 0.248
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -0.1777 0.0879 -2.02 0.0433
## F1 -0.0234 0.0595 -0.39 0.6945
## grupo=unicamp 0.0257 0.1172 0.22 0.8266
## # M3 do item 1
m3 <- lrm(
Q2019_001 ~ theta1*grupo,
data = banco.dif
)
m3## Logistic Regression Model
##
## lrm(formula = Q2019_001 ~ theta1 * grupo, data = banco.dif)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 1198 LR chi2 0.29 R2 0.000 C 0.530
## 0 648 d.f. 3 g 0.030 Dxy 0.060
## 1 550 Pr(> chi2) 0.9625 gr 1.030 gamma 0.061
## max |deriv| 7e-14 gp 0.007 tau-a 0.030
## Brier 0.248
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -0.1787 0.0880 -2.03 0.0423
## F1 -0.0461 0.0947 -0.49 0.6264
## grupo=unicamp 0.0253 0.1172 0.22 0.8290
## F1 * grupo=unicamp 0.0376 0.1218 0.31 0.7575
## Note que não necessariamente o theta vai ter coeficiente significativo. Isso pode acontecer nos casos em que a discriminação do item é baixa.
A comparação dos modelos é feita de duas maneiras:
Pelo teste da razão de verossimilhança, que verifica se o parâmetro acrescido a um modelo tem valor diferente de zero. Isso significa que a comparação é realizada entre dois modelos. O valor-p corresponde à área da cauda à direita dessa medida em uma distribuição qui-quadrado. O número de graus de liberdade será igual à diferença entre o número de parâmetros dos modelos comparados. Como são quatro cadernos e um é a referência, são três parâmetros adicionados em cada modelo
Pelo tamanho do efeito, uma medida de ajuste que se assemelha ao R2 calculado em regressões lineares.
Em estudo de simulação, Jodoin e Gierl (2001) estabeleceram que um item apresenta DIF caso o teste da razão de verossimilhança apresente significância ao nível de 0,05 e a medida do tamanho do efeito (diferença entre o pseudo R2 dos modelos) seja maior ou igual a 0,035. Adotaremosa esse critério.
Teste de razão de verossimilhança:
# desvio
deviance1 <- m1$deviance[2]
deviance2 <- m2$deviance[2]
deviance3 <- m3$deviance[2]
# diferença de graus de liberdade
df12 <- df23 <- 3
df13 <- 6
# qui-quadrado da razão de verossimilhança
chi12 <- 1 - pchisq(deviance1 - deviance2, df12) %>%
round(4)
chi13 <- 1 - pchisq(deviance1 - deviance3, df13) %>%
round(4)
chi23 <- 1 - pchisq(deviance2 - deviance3, df23) %>%
round(4)
chi12## [1] 0.9972chi13## [1] 0.9999chi23## [1] 0.9924Tamanho do efeito:
m2$stats[10] - m1$stats[10]## R2
## 5.353027e-05m3$stats[10] - m1$stats[10]## R2
## 0.0001599024m3$stats[10] - m2$stats[10]## R2
## 0.0001063722Para assinalarmos o item com DIF, devemos comparar os modelos dois a dois na combinação dos critérios.
(chi12 < .05) & (m2$stats[10] - m1$stats[10] > .035)## R2
## FALSE(chi13 < .05) & (m3$stats[10] - m1$stats[10] > .035)## R2
## FALSE(chi23 < .05) & (m3$stats[10] - m2$stats[10] > .035)## R2
## FALSEA função rundif do pacote lordif gera uma tabela com todas as estatísticas das comparações entre os modelos, o que facilita nossa vida.
# verificar as colunas que são os itens
itens <- names (banco.dif)[substr(names(banco.dif), 1, 5) == 'Q2019']
dif <- rundif(
item = itens,
resp = banco.dif[,itens],
theta = as.numeric(theta1),
gr = banco.dif$grupo,
criterion = 'CHISQR',
alpha = .05,
wt = NULL
)
head (dif$stats)## item ncat chi12 chi13 chi23 beta12
## 1 Q2019_001 1 0.8266 0.9308 0.7574 0.0393
## 2 Q2019_002 1 0.0000 0.0000 0.0091 0.0229
## 3 Q2019_003 1 0.0167 0.0209 0.1569 0.0125
## 4 Q2019_004 1 0.0001 0.0005 0.7970 0.0719
## 5 Q2019_005 1 0.5933 0.0724 0.0259 0.0018
## 6 Q2019_006 1 0.0112 0.0003 0.0016 0.1011
## pseudo12.McFadden pseudo13.McFadden pseudo23.McFadden
## 1 0.0000 0.0001 0.0001
## 2 0.0294 0.0337 0.0043
## 3 0.0035 0.0047 0.0012
## 4 0.0102 0.0102 0.0000
## 5 0.0002 0.0032 0.0030
## 6 0.0041 0.0105 0.0063
## pseudo12.Nagelkerke pseudo13.Nagelkerke
## 1 0.0001 0.0002
## 2 0.0498 0.0569
## 3 0.0059 0.0080
## 4 0.0174 0.0175
## 5 0.0003 0.0049
## 6 0.0073 0.0186
## pseudo23.Nagelkerke pseudo12.CoxSnell pseudo13.CoxSnell
## 1 0.0001 0.0000 0.0001
## 2 0.0071 0.0365 0.0417
## 3 0.0021 0.0044 0.0059
## 4 0.0001 0.0124 0.0124
## 5 0.0046 0.0002 0.0036
## 6 0.0112 0.0053 0.0135
## pseudo23.CoxSnell df12 df13 df23
## 1 0.0001 1 2 1
## 2 0.0052 1 2 1
## 3 0.0015 1 2 1
## 4 0.0001 1 2 1
## 5 0.0034 1 2 1
## 6 0.0082 1 2 1Agora sim, verificar os critérios
itens_dif1 <- unique(
c(
which ((dif$stats$chi12 <= .05) & (dif$stats$pseudo12.Nagelkerke >= .035)),
which ((dif$stats$chi13 <= .05) & (dif$stats$pseudo13.Nagelkerke >= .035)),
which ((dif$stats$chi23 <= .05) & (dif$stats$pseudo23.Nagelkerke >= .035))
)
)
itens_dif1## [1] 2 70 14 25 51 59 112Para facilitar, criei uma função para fazer essa análise.
analise.dif <- function (banco, grupo)
{
# calibrar
calib1 <- multipleGroup(
banco,
1,
group = grupo,
itemtype = '2PL',
TOL = .01,
invariance = c('free_means', 'free_var', names(banco))
)
# calcular o theta
theta1 <- fscores(calib1)
# verificar as colunas que são os itens
itens <- names (banco)
dif1 <- rundif(
item = itens,
resp = banco[,itens],
theta = as.numeric(theta1),
gr = grupo,
criterion = 'CHISQR',
alpha = .05,
wt = NULL
)
itens_dif1 <- unique(
c(
which ((dif1$stats$chi12 <= .05) & (dif1$stats$pseudo12.Nagelkerke >= .035)),
which ((dif1$stats$chi13 <= .05) & (dif1$stats$pseudo13.Nagelkerke >= .035)),
which ((dif1$stats$chi23 <= .05) & (dif1$stats$pseudo23.Nagelkerke >= .035))
)
)
itens_dif1
if (length(itens_dif1) == 0)
stop(return(cat('Não há itens marcados com DIF')))
# purificar a medida
# calibrar
calib2 <- multipleGroup(
banco,
1,
group = grupo,
itemtype = '2PL',
TOL = .01,
invariance = c('free_means', 'free_var', names(banco)[-c(itens_dif1)])
)
# calcular o theta
theta2 <- fscores(calib2)
dif2 <- rundif(
item = itens,
resp = banco[,itens],
theta = as.numeric(theta2),
gr = grupo,
criterion = 'CHISQR',
alpha = .05,
wt = NULL
)
itens_dif2 <- unique(
c(
which ((dif2$stats$chi12 <= .05) & (dif2$stats$pseudo12.Nagelkerke >= .035)),
which ((dif2$stats$chi13 <= .05) & (dif2$stats$pseudo13.Nagelkerke >= .035)),
which ((dif2$stats$chi23 <= .05) & (dif2$stats$pseudo23.Nagelkerke >= .035))
)
)
# verificar se há algum novo item com DIF. Se der TRUE, os resultados são iguais (não tem outro item com DIF)
while (all.equal(itens_dif1, itens_dif2) != TRUE) {
itens_dif1 <- itens_dif2
# calibrar
calib2 <- multipleGroup(
banco,
1,
group = grupo,
itemtype = '2PL',
TOL = .01,
invariance = c('free_means', 'free_var', names(banco)[-c(itens_dif1)])
)
# calcular o theta
theta2 <- fscores(calib2)
dif2 <- rundif(
item = itens,
resp = banco[,itens],
theta = as.numeric(theta2),
gr = grupo,
criterion = 'CHISQR',
alpha = .05,
wt = NULL
)
itens_dif2 <- unique(
c(
which ((dif2$stats$chi12 <= .05) & (dif2$stats$pseudo12.Nagelkerke >= .035)),
which ((dif2$stats$chi13 <= .05) & (dif2$stats$pseudo13.Nagelkerke >= .035)),
which ((dif2$stats$chi23 <= .05) & (dif2$stats$pseudo23.Nagelkerke >= .035))
)
)
}
message(
cat(
paste0(
'Os itens marcados com DIF são: ',
paste0 (itens_dif2, collapse = ', '
)
)
)
)
return(itens_dif2)
}Os resultados das análises são:
# DIF 2019
banco.dif <- select (banco_total, ID2, grupo, starts_with('Q2019')) %>%
na.omit()
dif19 <- analise.dif(
banco.dif[,-c(1,2)],
banco.dif$grupo
)
dif19 <- names(banco.dif)[(dif19+2)]
# DIF 2020
banco.dif <- select (banco_total, ID2, grupo, starts_with('Q2020')) %>%
na.omit()
dif20 <- analise.dif(
banco.dif[,-c(1,2)],
banco.dif$grupo
)
dif20 <- names(banco.dif)[(dif20+2)]
# DIF 2021
banco.dif <- select (banco_total, ID2, grupo, starts_with('Q2021')) %>%
na.omit()
dif21 <- analise.dif(
banco.dif[,-c(1,2)],
banco.dif$grupo
)
dif21 <- names(banco.dif)[(dif21+2)]dif19## [1] "Q2019_002" "Q2019_014" "Q2019_070" "Q2019_059"
## [5] "Q2019_112"dif20## [1] "Q2020_004" "Q2020_042" "Q2020_068" "Q2020_084"
## [5] "Q2020_114" "Q2020_079"dif21## [1] "Q2021_096"Com esses itens marcados com DIF, precisamos calibrar considerando que eles possuem parâmetros diferentes entre os grupos.
3.2 Calibração
Aqui o modelo escolhido é o de 2PL. Para a calibração foi utilizado o algoritmo EM (expectation-maximization) com 0,01 de tolerância como critério de convergência.
# selecionar os itens comuns, ou seja, os que apresentaram invariância (e tirar tamvém as duas primeiras variáveis: ID e grupo)
itens_comuns <- names (banco_total)[!(names (banco_total) %in% c(dif19, dif20, dif21))][-c(1,2)]calib <- multipleGroup(
banco_total[,-c(1,2)],
1,
group = banco_total$grupo,
itemtype = '2PL',
TOL = .01,
invariance = c('free_means', 'free_var', itens_comuns)
)##
## Call:
## multipleGroup(data = banco_total[, -c(1, 2)], model = 1, group = banco_total$grupo,
## itemtype = "2PL", invariance = c("free_means", "free_var",
## itens_comuns), TOL = 0.01)
##
## Full-information item factor analysis with 1 factor(s).
## Converged within 0.01 tolerance after 101 EM iterations.
## mirt version: 1.34
## M-step optimizer: nlminb
## EM acceleration: Ramsay
## Number of rectangular quadrature: 61
## Latent density type: Gaussian
##
## Log-likelihood = -247355.9
## Estimated parameters: 1434
## AIC = 496195.9; AICc = 496648.7
## BIC = 500695.4; SABIC = 498337.7pars_unesp <- data.frame(coef(calib, IRTpars = TRUE, simplify = TRUE)$unesp$items)
pars_unicamp <- data.frame(coef(calib, IRTpars = TRUE, simplify = TRUE)$unicamp$items)
DT::datatable(pars_unesp)DT::datatable(pars_unicamp)3.3 Propriedades das provas
Os valores de dificuldade dos itens variaram para a Unesp de -22.72 a 18.58. Para a Unicamp, variaram de -22.72 a 18.58. Abaixo, a distribuição dos itens de acordo com sua dificuldade, a distribuição da população e a curva de informação do teste considerando todos os itens.
Esses gráficos foram feitos com as respostas da calibração. Ou seja, o sujeito que fez 2019 e 2020 é considerado uma vez. Mais à frente no relatório, ele terá dois escores.
Para a distribuição dos itens, consideramos intervalos de um desvio padrão. Por exemplo, a barra centrada em 0,5 contém os itens que estão entre 0,0 e 1,0; a barra centrada em -1,5 contém os itens entre -2,0 e -1,0.