SIMULACION DE DATOS

n = 300 # Tamaño muestral
set.seed(123)

Evaluar en que rangos de valores del logit la probabilidad abarca todo el rango [0,1]:

logit_piloto <- seq(-20, 20, 1) ## Un rango arbitrario
p_piloto <- 1/(1+exp(-logit_piloto))
plot(logit_piloto, p_piloto, type = "o", pch = 16, col = "gray60", lwd = 2, main = "Rango de probabilidades para diferentes valores de Logit", xlab = "Logit", ylab  ="Probabilidad")
grid()

Se Observa que en el rango de logit = [-5,5] la probabilidad abarca todo el rango entre [0,1]. Por lo tanto, se usara aproximadamente dicho rango de logit

## logit muestreado aleatoriamente de distribucion normal
logit_fijo <- rnorm(n, mean = 0, sd = 2) 
b0 = -3
b1 = 1.5

De acuerdo a los coeficientes anteriores, se calcularan los valores que debería tomar “x” para que el logit abarque el rango establecido anteriormente

x = (logit_fijo -b0)/b1 ## logit = b0 + b1*x
p1 = 1/(1+exp(-logit_fijo))

A partir de los coeficientes anteriores, se genera un nuevo logit que incluya un error aleatorio con distribucion ~ N(0, 0.5)

logit_aleatorio <- b0 + b1*x + rnorm(n, mean = 0, sd = 0.5)
p2 <- 1/(1+exp(-logit_aleatorio))
plot(x, p2, main = "Logit con error aleatorio ~ N(0,0 .5)", pch = 16, col = "gray60",
     ylab = "Probabilidad") ## Se observa todo el rango de probabilidades entre [0,1]
points(x, p1, pch = 16, col = "violetred4")
grid() 

y1 = round(p2)

ESTIMACION A PARTIR DE DATOS SIMULADOS

mod <- glm(y1 ~ x, family = "binomial") 
summary(mod) 
## 
## Call:
## glm(formula = y1 ~ x, family = "binomial")
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.34339  -0.26763  -0.00928   0.17563   2.50187  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -9.1946     1.2170  -7.555 4.18e-14 ***
## x             4.5415     0.6002   7.566 3.84e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 415.23  on 299  degrees of freedom
## Residual deviance: 131.90  on 298  degrees of freedom
## AIC: 135.9
## 
## Number of Fisher Scoring iterations: 7

Los coeficientes son estadísticamente significativos pero incorectos: Los coeficientes simulados fueron b0 = -3; b1 = 3

Diagnóstico visual:

par(mfrow = c(2,2))
plot(mod)

El problema consiste en que el error aleatorio se añadio a nivel del logit. Sin embargo, las probabilidades estimadas a partir del logit en realidad son proporciones esperadas dado el valor de x, es decir E[Y|x], por lo que al redondear las probabilidades para obtener los valores de Y, en realidad no se esta incluyendo ningun error. El error deberia simularse mediante la asignacion Y=1 segun una distribucion binomial dependiente de x; es decir, cada “Yi” representa un ensayo de Bernoulli con probabilidad Yi = 1 dependiente de x.

y_binomial <- rbinom(n, 1, prob = p1 ) 

mod_binomial <- glm(y_binomial ~ x, family = "binomial")
summary(mod_binomial) ## Coeficiente de "x" significativo
## 
## Call:
## glm(formula = y_binomial ~ x, family = "binomial")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.9710  -0.7352   0.1133   0.6879   1.8581  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.2925     0.4134  -7.965 1.65e-15 ***
## x             1.7138     0.2026   8.460  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 415.55  on 299  degrees of freedom
## Residual deviance: 268.01  on 298  degrees of freedom
## AIC: 272.01
## 
## Number of Fisher Scoring iterations: 5
par(mfrow = c(2,2))
plot(mod_binomial)

su <- summary(mod_binomial)
su$coefficients[,1]
## (Intercept)           x 
##   -3.292509    1.713843
logit_estimado <- cbind(1,x) %*% su$coefficients[,1]
p_estimada <- 1/(1+exp(-logit_estimado))

logit_potencial <- b0 + b1*x^2
p_potencial <- 1/(1+exp(-logit_potencial))

par(mfrow = c(1,1))
# plot(x, p_potencial, main = "Probabilidades simuladas para el logit no lineal",
#      ylab = "Probabilidad", pch = 16, col = "darkseagreen3", cex.lab = 1.2)
# grid()
y_potencial <- rbinom(n, 1, p_potencial)

mod_potencial <- glm(y_potencial ~ x, family = "binomial")
logit_potencial_estimado <- cbind(1,x) %*% mod_potencial$coefficients
p_pot_est <- 1/(1+exp(-logit_potencial_estimado))

colores <- c("gray70", "lightblue3", "darkseagreen3", "wheat4")

par(mfrow = c(1,1))
plot(x, p1, main = "Sensibilidad a la violación del supuesto de linealidad",  col = colores[1], pch = 16,
     ylab = "Probabilidad", cex.lab = 1.2)
points(x, p_estimada, col = colores[2], pch = 16)
points(x, p_potencial, col = colores[3], pch = 16)
points(x, p_pot_est, col = colores[4], pch = 16)
grid()

leg_texto = c("Probabilidad verdadera (logit lineal)",
              "Probabilidad estimada (logit lineal)", 
              "Probabilidad verdadera (logit no lineal: x^2)",
              "Probabilidad estimada (logit no lineal: x^2)")
legend("bottomright", legend = leg_texto, col = colores, pch = 16, bty = "n", pt.cex = 1.5)

plot(x, logit_potencial, col = "darkseagreen3", main = "Logit vs X", pch = 16,
     ylab = "Logit", cex.lab = 1.2)
grid()
points(x, logit_potencial_estimado, col = "wheat4", pch = 16)
leg_texto <- c("Logit no lineal simulado", "Logit estimado")
legend("topleft", legend = leg_texto, col = c("darkseagreen3", "wheat4"),
       pch = 16, bty = "n", pt.cex = 1.5)