La regresión logística exacta se utiliza para modelar variables de resultado binarias en el que los log odds de los resultados se modela como una combinación lineal de las variables de predicción. Se utiliza cuando el tamaño de la muestra es demasiado pequeña para una regresión logística regular (que utiliza el estimador basado en la máxima verosimilitud-estándar) y/o cuando algunas de las celdas formadas por el resultado y la variable predictora categórica no tienen observaciones. Las estimaciones dadas por regresión logística exacta no son asintóticas.
En esta sección se trabajará con los siguientes paquetes: “coda” y “elrm”.
###Paquetes###
library(coda)
library(elrm)
library(graphics)Supongamos que estamos interesados en los factores que influyen en si un estudiante es admitido en una escuela de ingeniería. La variable de resultado es binario (0/1): admitir o no admitir. Las variables predictoras de interés, como el género del estudiante y si el estudiante tomó el curso de Cálculo Avanzado en la escuela secundaria. Debido a que la variable de respuesta es binaria, tenemos que utilizar un modelo que maneje los valores 0/1. También, debido a la cantidad de estudiantes que participan que es pequeña, se necesita un procedimiento que pueda realizar la estimación con una muestra pequeña.
Número de estudiantes admitidos (admit), número total de solicitantes (num), desglosados por género (female=la variable de sexo femenino), y si habían tomado o no el curso de cálculo (apcalc).
###Los datos son pocos y se leen directamente###
dat <- read.table(text = "
female apcalc admit num
0 0 0 7
0 0 1 1
0 1 0 3
0 1 1 7
1 0 0 5
1 0 1 1
1 1 0 0
1 1 1 6",
header = TRUE)La variable “num”" indica el peso de frecuencia. Utilizamos esta variable para expandir el conjunto de datos y luego observar algunas tablas de frecuencia.
###repetir los datos###
###Se expande la data en función a la variable "num"###
dat <- dat[rep(1:nrow(dat), dat$num), -4]
dat## female apcalc admit
## 1 0 0 0
## 1.1 0 0 0
## 1.2 0 0 0
## 1.3 0 0 0
## 1.4 0 0 0
## 1.5 0 0 0
## 1.6 0 0 0
## 2 0 0 1
## 3 0 1 0
## 3.1 0 1 0
## 3.2 0 1 0
## 4 0 1 1
## 4.1 0 1 1
## 4.2 0 1 1
## 4.3 0 1 1
## 4.4 0 1 1
## 4.5 0 1 1
## 4.6 0 1 1
## 5 1 0 0
## 5.1 1 0 0
## 5.2 1 0 0
## 5.3 1 0 0
## 5.4 1 0 0
## 6 1 0 1
## 8 1 1 1
## 8.1 1 1 1
## 8.2 1 1 1
## 8.3 1 1 1
## 8.4 1 1 1
## 8.5 1 1 1###Algunas tablas###
xtabs(~female + apcalc, data = dat)## apcalc
## female 0 1
## 0 8 10
## 1 6 6xtabs(~female + admit, data = dat)## admit
## female 0 1
## 0 10 8
## 1 5 7xtabs(~apcalc + admit, data = dat)## admit
## apcalc 0 1
## 0 12 2
## 1 3 13xtabs(~female + apcalc + admit, data = dat)## , , admit = 0
##
## apcalc
## female 0 1
## 0 7 3
## 1 5 0
##
## , , admit = 1
##
## apcalc
## female 0 1
## 0 1 7
## 1 1 6Las tablas revelan que 30 estudiantes aplican para el programa de Ingeniería. De ellos, 15 ingresaron y 15 se les negó el ingreso. Había 18 hombres y 12 mujeres solicitantes. Dieciséis de los solicitantes habían tomado cálculo AP y 14 no. Tenga en cuenta que todas las mujeres que tomaron cálculo AP fueron admitidas, frente a sólo el 70% de los varones.
Debido al pequeño tamaño de la muestra y la presencia de celdas sin datos, la regresión logística regular no es recomendable, y ni siquiera podría ser estimable.
###Los datos son pocos y se leen directamente###
d <- read.table(text = "
female frec num
0 0.11111 18
1 0.16666 12",
header = TRUE)###Método logístico###
M1 <- glm(frec ~ female, weights=num, data = d, family = binomial)
summary(M1)##
## Call:
## glm(formula = frec ~ female, family = binomial, data = d, weights = num)
##
## Deviance Residuals:
## [1] 0 0
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.079 0.750 -2.773 0.00556 **
## female 0.470 1.078 0.436 0.66292
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1.8901e-01 on 1 degrees of freedom
## Residual deviance: -4.4411e-16 on 0 degrees of freedom
## AIC: 8.9312
##
## Number of Fisher Scoring iterations: 4###Probabilidad estimada###
F=1 ###es el valor que puede tomar female###
exp(-2.079+0.470*F)/1+exp(-2.079+0.470*F)## [1] 0.4001752Corramos un (aproximado) análisis logístico exacta con el comando elrm en el paquete elrm. Esto se basa en el muestreo MCMC (Métodos de cadenas de Markov Monte Carlo).
x <- xtabs(~admit + interaction(female, apcalc), data = dat)
x## interaction(female, apcalc)
## admit 0.0 1.0 0.1 1.1
## 0 7 5 3 0
## 1 1 1 7 6cdat <- cdat <- data.frame(female = rep(0:1, 2), apcalc = rep(0:1, each = 2), admit = x[2, ], ntrials = colSums(x))
cdat## female apcalc admit ntrials
## 0.0 0 0 1 8
## 1.0 1 0 1 6
## 0.1 0 1 7 10
## 1.1 1 1 6 6Ahora podemos calcular la regresión logística aproximada utilizando elrm y muestreo MCMC. Se ejecutarán 22,000 iteraciones con un burning de 2,000 para una cadena final del 20,000. Tenga en cuenta que para el modelo combinado de female y apcalc, utilizamos una cadena de 5 millones. Esto es porque para la inferencia, cada efecto necesita al menos 1000, pero debido a que la distribución conjunta condicional es degenerado, para el efecto female la relación de ensayos utilizables es baja, lo que significa que para lograr más de 1000, las iteraciones total debe ser extremadamente alta.
###Se usa el comando elrm para usar el muestreo MCMC.###
###modelo con efecto female únicamente.###
m.female <- elrm(formula = admit/ntrials ~ female, interest = ~female, iter = 22000, dataset = cdat, burnIn = 2000)## Generating the Markov chain ...## Progress: 0% Progress: 5% Progress: 10% Progress: 15% Progress: 20% Progress: 25% Progress: 30% Progress: 35% Progress: 40% Progress: 45% Progress: 50% Progress: 55% Progress: 60% Progress: 65% Progress: 70% Progress: 75% Progress: 80% Progress: 85% Progress: 90% Progress: 95% Progress: 100%## Generation of the Markov Chain required 5 secs
## Conducting inference ...
## Inference required 1 secs###summary del modelo que incluye estimadores y los intervalos de confianza###
summary(m.female)##
## Call:
## [[1]]
## elrm(formula = admit/ntrials ~ female, interest = ~female, iter = 22000,
## dataset = cdat, burnIn = 2000)
##
##
## Results:
## estimate p-value p-value_se mc_size
## female 0.50696 0.4904 0.00595 20000
##
## 95% Confidence Intervals for Parameters
##
## lower upper
## female -1.153438 2.13057###trace plot and histogram of sampled values from the sufficient statistic.###
plot(m.female)###modelo con efecto apcalc únicamente.###
m.apcalc <- elrm(formula = admit/ntrials ~ apcalc, interest = ~apcalc, iter = 22000, dataset = cdat, burnIn = 2000)## Generating the Markov chain ...## Progress: 0% Progress: 5% Progress: 10% Progress: 15% Progress: 20% Progress: 25% Progress: 30% Progress: 35% Progress: 40% Progress: 45% Progress: 50% Progress: 55% Progress: 60% Progress: 65% Progress: 70% Progress: 75% Progress: 80% Progress: 85% Progress: 90% Progress: 95% Progress: 100%## Generation of the Markov Chain required 5 secs
## Conducting inference ...## Warning: 'apcalc' observed value of the sufficient statistic was not
## sampled## Inference required 0 secs###summary del modelo que incluye estimadores y los intervalos de confianza###
summary(m.apcalc)##
## Call:
## [[1]]
## elrm(formula = admit/ntrials ~ apcalc, interest = ~apcalc, iter = 22000,
## dataset = cdat, burnIn = 2000)
##
##
## Results:
## estimate p-value p-value_se mc_size
## apcalc NA 0 0 20000
##
## 95% Confidence Intervals for Parameters
##
## lower upper
## apcalc NA NA###trace plot and histogram of sampled values from the sufficient statistic.###
plot(m.apcalc)Tener en cuenta que esta técnica aproximada con burnin y las iteraciones suficiente es bastante similar a las estimaciones logísticas exactas de Stata.
. La regresión logística exacta puede ser una alternativa a la regresión logistica condicional.
. La regresión logística exacta es una alternativa para pequeñas muestras.