Simulando datos con el modelo de Regresión Binomial

Antes de pasar al ejemplo recordemos el modelo y sus parámetros.

\[log\left(\frac{p}{1-p}\right) = \beta_{0}+\beta_{1}x\]

donde,

\[p = \frac{1}{1+e^{-(\beta_{0}+\beta_{1}x)}} \] A continuación, se muestra un ejemplo de como simular observaciones de un modelo de regresión binomial:

Ejemplo

Crear una función que simule n observaciones del siguiente modelo que tiene el vector de parámetros:

\[\theta=(\begin{array} {rrr} \beta_{0}=-2,& \beta_{1}=3\\ \end{array})^\top \] A continuación se muestra la función solicitada la cual sólo tiene un argumento y que entrega como resultado un marco de datos con la información.

sim_logistic_data_1 = function(sample_size = n, beta_0 = -2, beta_1 = 3) {
  x = rnorm(n = sample_size)
  eta = beta_0 + beta_1 * x
  p = 1 / (1 + exp(-eta))
  y = rbinom(n = sample_size, size = 1, prob = p)
  data.frame(y, x)
}

Vamos a poner a prueba la función para simular 25 observaciones de la siguiente manera:

set.seed(666) 
datos_1 <- sim_logistic_data_1(sample_size = 20)
datos_1
##    y           x
## 1  1  0.75331105
## 2  1  2.01435467
## 3  0 -0.35513446
## 4  1  2.02816784
## 5  0 -2.21687445
## 6  0  0.75839618
## 7  0 -1.30618526
## 8  0 -0.80251957
## 9  0 -1.79224083
## 10 1 -0.04203245
## 11 1  2.15004262
## 12 0 -1.77023084
## 13 1  0.86465359
## 14 0 -1.72015590
## 15 0  0.13412567
## 16 0 -0.07582656
## 17 1  0.85830054
## 18 1  0.34490035
## 19 0 -0.58245269
## 20 1  0.78617038

Dado que solo tenemos una única variable predictora, podemos mostrar acontinuación gráficamente esta situación. Primero, tenga en cuenta que los datos se trazan con puntos negros. La respuesta \(y\) solo toma valores 0 y 1.

fit_glm_1 = glm(y ~ x, data = datos_1, family = binomial(link = "logit"))

plot(y ~ x, data =  datos_1, pch = 20, ylab = "Probabilidad estimada", 
     main = "Regresión logística")
grid()
curve(predict(fit_glm_1, data.frame(x), type = "response"),
      add = TRUE,col = "dodgerblue" )

Veamos nuestros coeficientes estimados.

round(coef(fit_glm_1), 1)
## (Intercept)           x 
##        -0.9         3.3

Veamos que \(\hat \beta_0 =-0.9\) y \(\hat \beta_1 =3.3\) no están cerca a los parametros tomados al inicio como, \(\beta_{0} = -2\) y \(\beta_{1} = 3\). Esto se debe a que el tamaño de la muestra es pequeño, recuerde que tomamos \(n = 20\).

Esto nos pone a pensar que muy seguramente con un número de observaciones mayor los estimadores serán mucho más proximos a los dados inicialmente.

Así, veamos que resultados obtenemos al realizar 200 observaciones.

A continuación se muestra la función solicitada y que entrega como resultado un marco de datos con la información.

sim_logistic_data_2 = function(sample_size = n, beta_0 = -2, beta_1 = 3) {
  x = rnorm(n = sample_size)
  eta = beta_0 + beta_1 * x
  p = 1 / (1 + exp(-eta))
  y = rbinom(n = sample_size, size = 1, prob = p)
  data.frame(y, x)
}

Vamos a poner a prueba la función para simular 200 observaciones de la siguiente manera:

set.seed(666) 
datos_2 <- sim_logistic_data_2(sample_size = 200)
datos_2

Nota: Por motivos de la cantidad de datos, estos no fueron incluidos en el documento, sin embargo, la idea es igual a los datos entregados anteriormente con \(n = 20\).

Acontinuación, veremos la grafica de 200 observaciones con una única variable predictora. Recuerde que los datos se trazan con puntos negros. La respuesta \(y\) solo toma valores 0 y 1.

fit_glm_2 = glm(y ~ x, data = datos_2, family = binomial(link = "logit"))

plot(y ~ x, data =  datos_2, pch = 20, ylab = "Probabilidad estimada", 
     main = "Regresión logística")
grid()
curve(predict(fit_glm_2, data.frame(x), type = "response"),
      add = TRUE,col = "red" )

Veamos nuestros coeficientes estimados.

round(coef(fit_glm_2), 1)
## (Intercept)           x 
##        -2.5         3.3

Veamos que \(\hat \beta_0 =-2.5\) y \(\hat \beta_1 =3.3\) son valores mucho más cercanos o parecidos a los parametros tomados al inicio como, \(\beta_{0} = -2\) y \(\beta_{1} = 3\).

Nota: A manera de conclusión, nos damos cuenta que a medida que \(n\) (número de observaciones) es más grande se recuperan o nos acercamos cada vez más a los verdaderos parametros.