Chapter 7 Simulação teste linear
# library (data.table)
# banco_final <- fread('banco_final_unesp_unicamp.csv', sep = ';', dec = ',')A simulação será feita sem os itens com DIF. Mantive a calibração, mas excluí os itens problemáticos. Os itens utilizados são os do objeto itens_comuns (sem DIF).
Antes, incluir informação sobre as áreas dos itens
# as áreas dos itens
areas <- c('basico', 'clinica', 'pediatria', 'cirurgia', 'go', 'scolet')
info_area <- rep(areas, each = 20)
info_area <- data.frame(
item = names(banco_total)[-c(1,2)],
area = rep(info_area, length.out = (ncol(banco_total)-2))
)7.1 Simulação 1
Nesta análise, os itens são divididos em seis grupos de acordo com suas áreas. De cada grupo são selecionados n/6 itens, onde n é tamanho da prova. Os itens selecionados são aqueles que produzem um teste com maior área de informação. Em seguida, simula-se a resposta de todos os sujeitos aos itens da prova montada. O novo theta é comparado com o theta original. Esse processo se repete 100 vezes para cada tamanho de prova. É importante dizer que a prova é montada pensando na maior informação para a região próxima do theta com maior probabilidade na função de densidade. Não usei o intervalo de -3 a +3 porque quando acima de 2.0 existe muito pouca gente. O intervalo ficou em 3dp abaixo e 3dp acima do ponto de maior densidade.
set.seed(1234)
# obter o ponto de maior probabilidade na distribuição do theta
dens <- density(banco_final$theta)$x[which (density(banco_final$theta)$y == max (density(banco_final$theta)$y))]
# o intervalo de dificuldade terá quantos desvios?
int_dp <- 6
# obter intervalo de int_dp unidades de desvio padrão
intervalo <- c(dens - int_dp/2*sd(banco_final$theta), dens + int_dp/2*sd(banco_final$theta))
# parâmetros dos itens
pars <- data.frame(coef(calib, IRTpars = TRUE, simplify = TRUE)$unicamp$items)
# quais são os números dos itens comuns?
ordem_itens_comuns <- (names (banco_total) %in% itens_comuns)[-c(1,2)]
head (ordem_itens_comuns)## [1] TRUE FALSE TRUE TRUE TRUE TRUEpars <- pars[ordem_itens_comuns,]
# tabela dos itens com as áreas
tab_itens <- data.frame(
item = itens_comuns,
a = pars$a,
b = pars$b
)
tab_itens <- left_join(tab_itens, info_area, by = 'item')
# criar modelo do mirt com grupo único só para acelerar o processo
itens_mirt <- data.frame(
a1 = tab_itens$a,
d = -tab_itens$a * tab_itens$b
)
mod_info <- mirtCAT::generate.mirt_object(
itens_mirt,
itemtype = '2PL'
)
# pegar a área da curva de informação de cada item no intervalo que queremos
infos <- sapply(1:nrow(tab_itens), areainfo, x = mod_info, theta_lim = c(intervalo[1], intervalo[2]), simplify = 'array')
tab_itens$info <- do.call(rbind, infos[1,3,])
# ordenar os itens pela informação
tab_itens <- arrange(tab_itens, desc(info))# tamanhos das provas
n_itens <- seq (60, 102, 6)
# quantidade de replicações
replicacao <- 100analises1 <- data.frame(matrix(
ncol = 5,
nrow = length(n_itens)*replicacao
)
)
names (analises1) <- c('tamanho', 'replicacao', 'reqm', 'cor', 'theta')
for (i in 1:length(n_itens))
{
# tamanho da prova
n <- n_itens[i]
# selecionar os itens de cada área com maiores áreas de informação
area1 <- tab_itens[tab_itens$area == areas[1],][1:(n/6),]
area2 <- tab_itens[tab_itens$area == areas[2],][1:(n/6),]
area3 <- tab_itens[tab_itens$area == areas[3],][1:(n/6),]
area4 <- tab_itens[tab_itens$area == areas[4],][1:(n/6),]
area5 <- tab_itens[tab_itens$area == areas[5],][1:(n/6),]
area6 <- tab_itens[tab_itens$area == areas[6],][1:(n/6),]
# montar a prova
itens_prova <- rbind (
area1,
area2,
area3,
area4,
area5,
area6
)
# montar tabela com os parâmetros
itens_prova_mirt <- data.frame(
a1 = itens_prova$a,
d = -itens_prova$a * itens_prova$b,
g = 0,
u = 1
)
# gerar o modelo mirt com esses itens
mod_prova <- generate.mirt_object(itens_prova_mirt, '2PL')
for(rep in 1:replicacao)
{
# simular a resposta das pessoas aos itens da prova selecionada
sim_prova <- data.frame(
INEPsico::simular(
banco_final$theta,
pars = matrix(
c(
itens_prova$a,
itens_prova$b,
rep(0, nrow(itens_prova))
),
ncol = 3
)
)
)
# calcular a nota da simulação
theta_sim <- data.frame (fscores(mod_prova, response.pattern = sim_prova))
# verificar a raiz do erro quadrático médio
analises1$reqm[(replicacao*(i-1) + rep)] <- sqrt(mse(banco_final$theta, theta_sim$F1))
# verificar a correlação
analises1$cor[replicacao*(i-1) + rep] <- cor (banco_final$theta, theta_sim$F1)
# indicar o número da replicação
analises1$replicacao[replicacao*(i-1) + rep] <- rep
# indicar o tamanho amostral
analises1$tamanho[replicacao*(i-1) + rep] <- n
# indicar o theta
analises1$theta[replicacao*(i-1) + rep] <- data.frame(theta_sim$F1)
}
}O resultado é apresentado a seguir. Calculei a média da raiz do erro quadrático médio e da correlação para as replicações de cada tamanho de prova. A raiz do erro quadrático médio equivale à raiz da média do quadrado da diferença entre o valor real e o observado. Ou seja, pega o valor real, diminui do observado. Eleva ao quadrado. Calcula o somatório disso para todos os thetas e divide pelo número de observações. Depois, tira a raiz. Como eleva ao quadrado e tira a raiz, está tudo na mesma escala. A precisão foi calculada pela fórmula 1-erro^2.
resultado <- data.frame(matrix (ncol = 3, nrow = length(n_itens)))
names (resultado) <- c('n_itens', 'reqm', 'cor')
for (i in 1:length(n_itens))
{
resultado$n_itens[i] <- n_itens[i]
resultado$reqm[i] <- mean (subset (analises1, analises1$tamanho == n_itens[i])$reqm)
resultado$cor[i] <- mean (subset (analises1, analises1$tamanho == n_itens[i])$cor)
}
resultado$precisao <- 1 - resultado$reqm^2
resultado## n_itens reqm cor precisao
## 1 60 0.2743956 0.9702997 0.9247070
## 2 66 0.2629267 0.9726732 0.9308696
## 3 72 0.2538735 0.9744936 0.9355483
## 4 78 0.2441207 0.9763567 0.9404051
## 5 84 0.2358387 0.9778912 0.9443801
## 6 90 0.2287331 0.9791508 0.9476812
## 7 96 0.2209086 0.9805089 0.9511994
## 8 102 0.2157290 0.9814005 0.95346107.2 Simulação 2
Nesta simulação, selecionamos o segundo grupo de n melhores itens do banco.
analises2 <- data.frame(matrix(
ncol = 5,
nrow = length(n_itens)*replicacao
)
)
names (analises2) <- c('tamanho', 'replicacao', 'reqm', 'cor', 'theta')
for (i in 1:length(n_itens))
{
# tamanho da prova
n <- n_itens[i]
# selecionar os itens de cada área com maiores áreas de informação
# neste caso, não são os mais informativos, mas os segundos mais informativos
area1 <- tab_itens[tab_itens$area == areas[1],][((n/6)+1):(2*(n/6)),]
area2 <- tab_itens[tab_itens$area == areas[2],][((n/6)+1):(2*(n/6)),]
area3 <- tab_itens[tab_itens$area == areas[3],][((n/6)+1):(2*(n/6)),]
area4 <- tab_itens[tab_itens$area == areas[4],][((n/6)+1):(2*(n/6)),]
area5 <- tab_itens[tab_itens$area == areas[5],][((n/6)+1):(2*(n/6)),]
area6 <- tab_itens[tab_itens$area == areas[6],][((n/6)+1):(2*(n/6)),]
# montar a prova
itens_prova <- rbind (
area1,
area2,
area3,
area4,
area5,
area6
)
# montar tabela com os parâmetros
itens_prova_mirt <- data.frame(
a1 = itens_prova$a,
d = -itens_prova$a * itens_prova$b,
g = 0,
u = 1
)
# gerar o modelo mirt com esses itens
mod_prova <- generate.mirt_object(itens_prova_mirt, '2PL')
for(rep in 1:replicacao)
{
# simular a resposta das pessoas aos itens da prova selecionada
sim_prova <- data.frame(
INEPsico::simular(
banco_final$theta,
pars = matrix(
c(
itens_prova$a,
itens_prova$b,
rep(0, nrow(itens_prova))
),
ncol = 3
)
)
)
# calcular a nota da simulação
theta_sim <- data.frame (fscores(mod_prova, response.pattern = sim_prova))
# verificar a raiz do erro quadrático médio
analises2$reqm[(replicacao*(i-1) + rep)] <- sqrt(mse(banco_final$theta, theta_sim$F1))
# verificar a correlação
analises2$cor[replicacao*(i-1) + rep] <- cor (banco_final$theta, theta_sim$F1)
# indicar o número da replicação
analises2$replicacao[replicacao*(i-1) + rep] <- rep
# indicar o tamanho amostral
analises2$tamanho[replicacao*(i-1) + rep] <- n
# indicar o theta
analises2$theta[replicacao*(i-1) + rep] <- data.frame(theta_sim$F1)
}
}Resultado:
resultado <- data.frame(matrix (ncol = 3, nrow = length(n_itens)))
names (resultado) <- c('n_itens', 'reqm', 'cor')
for (i in 1:length(n_itens))
{
resultado$n_itens[i] <- n_itens[i]
resultado$reqm[i] <- mean (subset (analises2, analises2$tamanho == n_itens[i])$reqm)
resultado$cor[i] <- mean (subset (analises2, analises2$tamanho == n_itens[i])$cor)
}
resultado$precisao <- 1 - resultado$reqm^2
resultado## n_itens reqm cor precisao
## 1 60 0.3072122 0.9614824 0.9056206
## 2 66 0.3014442 0.9629022 0.9091314
## 3 72 0.2950019 0.9644590 0.9129739
## 4 78 0.2912104 0.9653632 0.9151965
## 5 84 0.2863283 0.9664617 0.9180161
## 6 90 0.2825177 0.9673144 0.9201838
## 7 96 0.2813081 0.9676125 0.9208657
## 8 102 0.2785001 0.9682164 0.92243777.3 Simulação 3
Nesta simulação, selecionamos os n ímpares melhores itens do banco.
analises3 <- data.frame(matrix(
ncol = 5,
nrow = length(n_itens)*replicacao
)
)
names (analises3) <- c('tamanho', 'replicacao', 'reqm', 'cor', 'theta')
for (i in 1:length(n_itens))
{
# tamanho da prova
n <- n_itens[i]
# selecionar os itens de cada área com maiores áreas de informação
# neste caso, vamos selecionar os primeiros itens ímpares para deixar a CAT com informativos também
area1 <- tab_itens[tab_itens$area == areas[1],][seq (1, by = 2, length.out = n/6),]
area2 <- tab_itens[tab_itens$area == areas[2],][seq (1, by = 2, length.out = n/6),]
area3 <- tab_itens[tab_itens$area == areas[3],][seq (1, by = 2, length.out = n/6),]
area4 <- tab_itens[tab_itens$area == areas[4],][seq (1, by = 2, length.out = n/6),]
area5 <- tab_itens[tab_itens$area == areas[5],][seq (1, by = 2, length.out = n/6),]
area6 <- tab_itens[tab_itens$area == areas[6],][seq (1, by = 2, length.out = n/6),]
# montar a prova
itens_prova <- rbind (
area1,
area2,
area3,
area4,
area5,
area6
)
# montar tabela com os parâmetros
itens_prova_mirt <- data.frame(
a1 = itens_prova$a,
d = -itens_prova$a * itens_prova$b,
g = 0,
u = 1
)
# gerar o modelo mirt com esses itens
mod_prova <- generate.mirt_object(itens_prova_mirt, '2PL')
for(rep in 1:replicacao)
{
# simular a resposta das pessoas aos itens da prova selecionada
sim_prova <- data.frame(
INEPsico::simular(
banco_final$theta,
pars = matrix(
c(
itens_prova$a,
itens_prova$b,
rep(0, nrow(itens_prova))
),
ncol = 3
)
)
)
# calcular a nota da simulação
theta_sim <- data.frame (fscores(mod_prova, response.pattern = sim_prova))
# verificar a raiz do erro quadrático médio
analises3$reqm[(replicacao*(i-1) + rep)] <- sqrt(mse(banco_final$theta, theta_sim$F1))
# verificar a correlação
analises3$cor[replicacao*(i-1) + rep] <- cor (banco_final$theta, theta_sim$F1)
# indicar o número da replicação
analises3$replicacao[replicacao*(i-1) + rep] <- rep
# indicar o tamanho amostral
analises3$tamanho[replicacao*(i-1) + rep] <- n
# indicar o theta
analises3$theta[replicacao*(i-1) + rep] <- data.frame(theta_sim$F1)
}
}Resultado:
resultado <- data.frame(matrix (ncol = 3, nrow = length(n_itens)))
names (resultado) <- c('n_itens', 'reqm', 'cor')
for (i in 1:length(n_itens))
{
resultado$n_itens[i] <- n_itens[i]
resultado$reqm[i] <- mean (subset (analises3, analises3$tamanho == n_itens[i])$reqm)
resultado$cor[i] <- mean (subset (analises3, analises3$tamanho == n_itens[i])$cor)
}
resultado$precisao <- 1 - resultado$reqm^2
resultado## n_itens reqm cor precisao
## 1 60 0.2798497 0.9686696 0.9216841
## 2 66 0.2708924 0.9706426 0.9266173
## 3 72 0.2609067 0.9726886 0.9319277
## 4 78 0.2548861 0.9738986 0.9350331
## 5 84 0.2479853 0.9752227 0.9385033
## 6 90 0.2423142 0.9763064 0.9412839
## 7 96 0.2375142 0.9772452 0.9435870
## 8 102 0.2329574 0.9780815 0.94573097.4 Simulação 4
Nesta simulação, selecionamos os n pares melhores itens do banco.
analises4 <- data.frame(matrix(
ncol = 5,
nrow = length(n_itens)*replicacao
)
)
names (analises4) <- c('tamanho', 'replicacao', 'reqm', 'cor', 'theta')
for (i in 1:length(n_itens))
{
# tamanho da prova
n <- n_itens[i]
# selecionar os itens de cada área com maiores áreas de informação
# neste caso, vamos selecionar os primeiros itens ímpares para deixar a CAT com informativos também
area1 <- tab_itens[tab_itens$area == areas[1],][seq (2, by = 2, length.out = n/6),]
area2 <- tab_itens[tab_itens$area == areas[2],][seq (2, by = 2, length.out = n/6),]
area3 <- tab_itens[tab_itens$area == areas[3],][seq (2, by = 2, length.out = n/6),]
area4 <- tab_itens[tab_itens$area == areas[4],][seq (2, by = 2, length.out = n/6),]
area5 <- tab_itens[tab_itens$area == areas[5],][seq (2, by = 2, length.out = n/6),]
area6 <- tab_itens[tab_itens$area == areas[6],][seq (2, by = 2, length.out = n/6),]
# montar a prova
itens_prova <- rbind (
area1,
area2,
area3,
area4,
area5,
area6
)
# montar tabela com os parâmetros
itens_prova_mirt <- data.frame(
a1 = itens_prova$a,
d = -itens_prova$a * itens_prova$b,
g = 0,
u = 1
)
# gerar o modelo mirt com esses itens
mod_prova <- generate.mirt_object(itens_prova_mirt, '2PL')
for(rep in 1:replicacao)
{
# simular a resposta das pessoas aos itens da prova selecionada
sim_prova <- data.frame(
INEPsico::simular(
banco_final$theta,
pars = matrix(
c(
itens_prova$a,
itens_prova$b,
rep(0, nrow(itens_prova))
),
ncol = 3
)
)
)
# calcular a nota da simulação
theta_sim <- data.frame (fscores(mod_prova, response.pattern = sim_prova))
# verificar a raiz do erro quadrático médio
analises4$reqm[(replicacao*(i-1) + rep)] <- sqrt(mse(banco_final$theta, theta_sim$F1))
# verificar a correlação
analises4$cor[replicacao*(i-1) + rep] <- cor (banco_final$theta, theta_sim$F1)
# indicar o número da replicação
analises4$replicacao[replicacao*(i-1) + rep] <- rep
# indicar o tamanho amostral
analises4$tamanho[replicacao*(i-1) + rep] <- n
# indicar o theta
analises4$theta[replicacao*(i-1) + rep] <- data.frame(theta_sim$F1)
}
}Resultado:
resultado <- data.frame(matrix (ncol = 3, nrow = length(n_itens)))
names (resultado) <- c('n_itens', 'reqm', 'cor')
for (i in 1:length(n_itens))
{
resultado$n_itens[i] <- n_itens[i]
resultado$reqm[i] <- mean (subset (analises4, analises4$tamanho == n_itens[i])$reqm)
resultado$cor[i] <- mean (subset (analises4, analises4$tamanho == n_itens[i])$cor)
}
resultado$precisao <- 1 - resultado$reqm^2
resultado## n_itens reqm cor precisao
## 1 60 0.2856978 0.9671389 0.9183768
## 2 66 0.2773147 0.9690566 0.9230965
## 3 72 0.2694112 0.9707747 0.9274176
## 4 78 0.2628051 0.9721931 0.9309335
## 5 84 0.2551876 0.9737286 0.9348793
## 6 90 0.2483627 0.9750609 0.9383159
## 7 96 0.2438717 0.9759508 0.9405266
## 8 102 0.2395258 0.9767718 0.9426274