| log.dosis | expuestos | muertos |
|---|---|---|
| 1.6907 | 59 | 6 |
| 1.7242 | 60 | 13 |
| 1.7552 | 62 | 18 |
| 1.7842 | 56 | 28 |
| 1.8113 | 63 | 52 |
| 1.8369 | 59 | 53 |
| 1.8610 | 62 | 61 |
| 1.8839 | 60 | 60 |
Modelos Lineales Generalizados (GLM)
Ajuste de un modelo lineal generalizado, estimación y contraste de hipótesis de las funciones base, y parámetros de un modelo lineal generalizado
Modelos lineares: Traje a la medida (distribución)
Función Link para algunas distribuciones
Distribución binomial
Datos binomiales no siguen una distribución normal: observamos un pico en 0 y un pico en 1, pero nada enmedio (imagen de la derecha).
Al modelar esta variable con lm() y graficar los valores esperados con el modelo contra los observados no siguen una línea.
Modelos lineares: Traje a la medida (link) para Datos binomiales
Modelo logístico para la mortalidad de escarabajos
En un estudio de toxicología, se está interesado en la tasa de mortalidad de escarabajos expuestos a disulfuro de carbono gaseoso. La Tabla muestra el número de escarabajos muertos después de cinco horas de exposición de este tóxico a diferentes concentraciones del químico. Podemos ver que a medida que aumenta la dosis, la mortalidad de escarabajos va aumentando.
Una alternativa es utilizar una transformación que describa la relación de forma adecuada y que garantice que la respuesta se encuentre dentro del espacio del parámetro. Es decir, que la estimación de la proporción de escarabajos muertos esté entre 0 y 1.
Sea: \(y_{ij}=1\) si el escarabajo \(j\) expuesto a las dosis \(i\), muere; \(0\) si lo contrario.
Por lo que para la i-ésima dosis, el número de escarabajos que muere es igual a \(y^*_i=\sum^{n_1}_{j=1}y_{ij}\). Además asumiendo independencia, podemos decir que:
\[y^*\sim binomial(n_i,\pi_i), i=1,...,n_i\]
donde \(0<\pi_1<1\) depende de la dosis. Definiendo \(y_i=y^*_i/n_i\) (la proporción de escarabajos muertos para la dosis \(i\)), tenemos que:
\(E(y_i|dosis_i)=\pi_i\) y \(V(y_i|dosis_i)=\pi_i(1-\pi_i)/n_i\)
Dado que \(0\leq \pi_i \leq 1\), se podría proponer que:
\[\pi=g^{-1}(dosis,\beta)=\frac{exp(\beta_0+log(dosis_i)\beta_1)}{1+exp(\beta_0+log(dosis_i)\beta_1)}=\frac{exp(\eta_i)}{1+exp(\eta_i)}\]
Por lo que \(g(\pi_i)=log(\frac{\pi_i}{1-\pi_i})=\eta_i\). Esta función recibe el nombre de función “logit” y su forma la podemos ver gráficamente en la Figura. Como vemos, está función está acotada entre 0 y 1, garantizando que las probabilidades estimadas siempre estén en este rango, y tiene forma de ‘S’ o sigmoidal.
Ajuste e interpretación de parámetros del modelo logístico
Para los datos de mortalidad de escarabajos, tenemos el siguiente modelo:
\(y^*_i\sim binomial(\eta_i,\pi_i)\), donde \(\pi_i=\frac{exp(\beta_0+\beta_1*log(dosis_i))}{1+exp(\beta_0+\beta_1*log(dosis_i))}\)
La estimación de los parámetros del GLM por MLE se puede hacer utilizando la función glm(). En esta función debemos determinar la distribución de probabilidad asociada a los datos y la función de enlace con el argumento family. Para el modelo logístico, usamos family=binomial (por defecto la función de enlace es logit):
Call:
glm(formula = cbind(dead, n - dead) ~ logdose, family = binomial)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -60.717 5.181 -11.72 <2e-16 ***
logdose 34.270 2.912 11.77 <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: 284.202 on 7 degrees of freedom
Residual deviance: 11.232 on 6 degrees of freedom
AIC: 41.43
Number of Fisher Scoring iterations: 4Aquí vemos que el efecto del logarítmo de la dosis es positivo y significativo. Es decir que, un incremento en la dosis significa un aumento en la probabidad de que el escarabajo muera. La relación la podemos ver de forma gráfica en la figura.
Datos de mortalidad de escarabajos. Estimación de la probabilidad de que el escarabajo muera en función del logarítmo de la dosis de disulfuro de carbono gaseoso.
La interpretación de los coeficientes del modelo logístico se hacen a través de los odds ratio. Primero, un odd está definido como: \[odd=\frac{P(Y=1)}{P(Y=0)}\]
Es decir, es una razón entre la probabilidad de éxito sobre la probabilidad de fracaso. Por ejemplo, si el odd es igual a dos, estaríamos diciendo que la probabilidad de éxito es el doble de la probabilidad de fracaso. Para el caso del modelo logístico simple, tenemos que:
\[odd(x)=\frac{P(Y=1|x)}{P(Y=0|x)}=\frac{\frac{exp(\beta_0+x*\beta_1)}{1+exp(\beta_0+x*\beta_1)}}{1-\frac{exp(\beta_0+x*\beta_1)}{1+exp(\beta_0+x*\beta_1)}}=exp(\beta_0+x*\beta_1)\]
Ahora, los odds ratio están definido como \(OR=\frac{odd(x=a+1)}{odd(x=a)}\). Para el caso del modelo logístico simple tenemos que el OR es:
\[OR(x)=\frac{odd(x=a+1)}{odd(x=a)}=\frac{exp[\beta_0+(a+1)*\beta_1]}{exp(\beta_0+a*\beta_1)}=exp(\beta_1)\]
Por lo tanto, por cada cambio unitario en \(x\), los odds de morir incrementan por un factor de \(exp(\beta_1)\). Por un cambio en \(x\) de \(a\) a \(a+\delta\), tenemos que \(OR=exp(\delta*\beta_1)\) Para los datos de los escarabajos un cambio unitario en el log(dosis) es muy grande, por lo que podríamos interpretar \(\beta_1\) usando un aumento en la log dosis mas pequeño, por ejemplo de \(\delta=0.02\)
\[exp(0.01\times34.270)=1.984\]
Es decir, si el log de la concentración de aumenta en
0.02, entonces la probabilidad de que el escarabajo muera aumenta casi el doble.
Ajuste e interpretación de parámetros del modelo Poisson
Para el ensayo clínico en epilepsia, el modelo propuesto es:
\(y_i \sim Poisson(\lambda_i)\), donde \(\lambda_i=exp(\beta_0+\beta_1treat_i+\beta_2base_i+\beta_2base_itreat_i)\)
Para el ajuste del modelo usamos la función glm() con el argumento family=poisson.
Call:
glm(formula = seizure.rate ~ treatment + base + base:treatment,
family = poisson, data = epilepsy4)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.2958397 0.0569022 22.773 < 2e-16 ***
treatmentProgabide -0.2492267 0.0767807 -3.246 0.00117 **
base 0.0214484 0.0009171 23.387 < 2e-16 ***
treatmentProgabide:base 0.0004227 0.0010777 0.392 0.69489
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 2521.75 on 235 degrees of freedom
Residual deviance: 989.84 on 232 degrees of freedom
AIC: 1763.8
Number of Fisher Scoring iterations: 5Los resultados del ajuste muestran que el número de ataques pre-tratamiento tiene un efecto positivo sobre la media del número de ataques epilépticos. Es decir, esta aumenta entre más ataques previos tenga el paciente. Mientras que, el tratamiento parece tener un efecto negativo. Es decir, el prograbide reduce la media del número de ataques epilépticos con respecto al placebo. Finalmente, el efecto interacción es muy pequeño sugiriendo que el efecto del tratamiento no se ve afectado por la frecuencia de los ataques epilépticos que tiene el paciente. La estimación de la media se puede observa de forma gráfica en la Figura.
Datos de epilepsia. Estimación de la media de ataques epilépticos (por dos semanas) luego de cuatro semanas de tratamiento. La línea negra representa la media de ataques para el grupo placebo, mientras que, la linea roja la media de ataques para el tratamiento con Progabide
En el modelo Poisson simple, tenemos que \(E(Y|x)=exp(\beta_0+x\beta_1)\). Si aumentamos a \(x\) en \(\delta\) unidades, tenemos que \(E(Y|x+\delta)=exp(\beta_0+(x+\delta)\beta_1)\). Ahora si calculamos el logaritmo de la razón de los valores esperados obtenemos que:
\[log[\frac{E(Y|x=a+\delta)}{E(Y|x=a)}]=\delta\beta_1\]
Por lo tanto, \(exp(\delta\beta_1)=\frac{E(Y|x=a+\delta)}{E(Y|x=a)}\). Es decir que, \(exp(\delta\beta_1)\) es una tasa de crecimiento del valor esperado de Y por un aumento de \(x\) en \(\delta\) unidades.
Por ejemplo, si el número de ataques epilépticos pre-tratamiento aumenta 20 casos, entonces el valor esperado de ataques epilépticos post-tramiento aumenta un 52% \((exp(20\times0.021)=1.522)\). Además, el tramiento con Progabida reduce en un poco más del 40% el número de episodios de epilepsia \((exp(-0.3635857)=0.695)\). Finalmente, el efecto de los casos pre-tramiento es el mismo para los pacientes en el grupo tratamiento y el control \((exp(0.0009)\approx1)\)
Mortalidad de escarabajos
La estimación de la probabilidad de que el escarabajo muera se puede observar en la Figura
Intervalo del 95% de confianza para la probabilidad de que el escarabajo muera.
Ataques epilépticos
Los intervalos del 95% de confianza para los parámetros del modelo son:
Aquí vemos que el efecto del tratamiento está entre -0.4 y -0.09. Es decir, que la reducción del número de ataques epilépticos está entre el 10% y el 49%. Mientras que, como ya se mencionó antes, el efecto interacción no tiene un aporte significativo.
La media del número de ataques epiléptico por dos semanas se puede observar en la Figura . Se puede observar que las estimaciones de la media para el grupo tratamiento es inferior que para el grupo control.
Estimación de la media de ataques epilépticos por dos semanas en el grupo placebo (línea negra) y el grupo tratamiento (línea roja).