Hasta el momento hemos visto el modelo lineal más sencillo, que funciona bajo el supuesto distribucional de normalidad. Los modelos lineales generalizados (MLG), fueron desarrollados por Nelder - Wedderburn (1972), permiten ampliar la gama de distribuciones de la variable respuesta a todas aquellas que pertenezcan a la familia exponencial de densidades.
Al igual que la regresión lineal múltiple, permite cumplir con dos objetivos, determinar si existe relación lineal entre las \(x_j\) y \(y\) y cuál es su magnitud.
Bajo el supuesto \(e_i\sim N(0,\sigma^2)\): \[Y_i=\beta_1X_{i1}+\dots+\beta_pX_{ip}+e_i\qquad i=1,\dots,n\]
\[E(Y_i)=\beta_1X_{i1}+\dots+\beta_pX_{ip}\qquad i=1,\dots,n\]
Es posible:
Encontrar \(\hat{\beta_j}\) via máxima verosimilitud (MV), \(\hat{\boldsymbol{\beta}}=(\boldsymbol{X}^t\boldsymbol{X})^{-1}\boldsymbol{X}^t\boldsymbol{Y}\)
Determinar la distribución de \(\hat{\beta_j}\), \(\hat{\boldsymbol{\beta}}\sim N_{p}(\boldsymbol{\beta},\sigma^2 (\boldsymbol{X}^t\boldsymbol{X})^{-1})\)
\[f(y,\theta,\phi)=\exp\left\lbrace\phi \left( y\theta-b(\theta)\right) +c(y,\phi)\right\rbrace\]
Estimación MV. \(\hat{\beta}_j\).
Distribución asintótica de \(\hat{\beta}_j\).
| Clásico | MLG | |
|---|---|---|
| Estimación | MV-Mín cuadrados | MV (Métodos numéricos) |
| Pruebas | SC-F-Wald | Deviance-\(\chi^2\)- \(F\) |
| ANOVA | Análisis de desvío | |
| Calidad del ajuste | \(R^2\)-AIC | Desvío-AIC |
Los MLGs se tienen tres componentes:
Componente aleatoria: \(\boldsymbol{Y}^t=(Y_1,\dots,Y_n)\), \(f(y,\theta,\phi)\) familia exponencial, \(Y_i's\) independientes.
Componente sistemática: Predictor lineal: \[\boldsymbol{\eta}=\boldsymbol{X}\boldsymbol{\beta}\] \[\eta_i=x_i^t\boldsymbol{\beta}\]
con \(x_i=(x_{i1},\dots,x_{i1})\) y \(\boldsymbol{\beta}\). Función de enlace: \[g(\mu_i)=\eta_i\]
Definición: Una familia de distribuciones es una familia exponencial sii: \[f(y,\theta,\phi)=\exp\left\lbrace \phi \left( y\theta-b(\theta)\right)+c(y,\phi) \right\rbrace\]
Resultado 1: Sea \(Y\sim f(y_i,\theta_i,\phi)\), con \(f(y_i,\theta_i,\phi)\) una familia exponencial de densidades, la función generadora de momentos de \(Y\) está dada por:
\[M_{Y_i}(t)=\exp\left\lbrace \phi\left[ b\left(\frac{t}{\phi}+\theta_i \right)-b(\theta_i)\right] \right\rbrace\] Corolario 1: Sea \(Y\sim f(y_i,\theta_i,\phi)\), con \(f(y_i,\theta_i,\phi)\) una familia exponencial de densidades.
\(E(Y_i)=\mu_i=b'(\theta_i)\)
\(V(Y_i)=\phi^{-1}b''(\theta_i)=\phi^{-1}V(\mu_i)=\phi^{-1}V_i\), con \(V(\mu_i)=\frac{d\mu_i}{d\theta_i}\) y \(\phi^{-1}>0\) el parámetro de dispersión.
Corolario 2: Sea \(Y\sim f(y_i,\theta_i,\phi)\), con \(f(y_i,\theta_i,\phi)\) una familia exponencial de densidades, entonces se verifican las condiciones de regularidad:
\[f(y,\theta,\phi)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left\lbrace-\frac{1}{2\sigma^2}(y-\mu)^2 \right\rbrace\] \[f(y,\theta,\phi)=\exp\left\lbrace \phi \left( y\theta-b(\theta)\right) \right\rbrace+c(y,\phi)\]
Con \(\theta=\mu\), \(b(\theta)=\frac{\theta^2}{2}\), \(\phi=\frac{1}{\sigma^2}\), \(c(y,\phi)=-\frac{1}{2}\left(\frac{y^2}{\sigma^2} +\log(2\pi\sigma^2)\right)\), \(V(\mu)=b''(\theta)=1\).
\[f(y)=\frac{e^{-\mu}\mu^y}{y!}=\exp\left\lbrace y\log\mu-\mu-\log y! \right\rbrace\] \[f(y)=\exp\left\lbrace \phi \left( y\theta-b(\theta)\right) \right\rbrace+c(y,\phi)\] Con \(\theta=\log\mu\), \(b(\theta)=\mu=e^\theta\),\(\phi=1\), \(c(y,\phi)=-\log y!\), \(V(\mu)=b''(\theta)=e^\theta=\mu\).
Se dice que un MLG tiene enlace canónico si: \[\theta_i=\eta_i=\sum_{j=1}^{p}\beta_jx_{ij}\] Cuando trabajamos con el enlace canónico garantizamos dos propiedades importantes, en primer lugar garantizamos que se tenga una única estimación máximo verosímil y en segundo lugar que dichos estimadores sean suficientes. Veamos, la función de verosimilitud está dada por: \[\mathfrak{l}(\boldsymbol{\beta})=\exp\left\lbrace \phi\sum_{i=1}^{n} \left( y_i\theta_i-b(\theta_i)\right)+\sum_{i=1}^{n}c(y_i,\phi)\right\rbrace \] Como \(\theta_i=\eta_i=\sum_{j=1}^{p}\beta_jx_{ij}\) y sea \(s_j=\phi\sum_{i=1}^{n}y_ix_{ij}\):
\[\mathfrak{l}(\boldsymbol{\beta})=\exp\left\lbrace \sum_{j=1}^{p}s_j\beta_j-\phi\sum_{i=1}^{n}b\left(\sum_{j=1}^{p}x_{ij}\beta_j\right)+\sum_{i=1}^{n}c(y_i,\phi)\right\rbrace\] \[\mathfrak{l}(\boldsymbol{\beta})=\exp\left\lbrace \sum_{j=1}^{p}s_j\beta_j-\phi\sum_{i=1}^{n}b\left(\sum_{j=1}^{p}x_{ij}\beta_j\right)+\sum_{i=1}^{n}c(y_i,\phi)\right\rbrace\] \[\mathfrak{l}(\boldsymbol{\beta})=g_{\boldsymbol{\beta}}(S)h(\boldsymbol{y})\] \(\Rightarrow S=(S_1,\dots,S_p)\) es suficiente minimal para \(\boldsymbol{\beta}\). Adicionalmente, el enlace canónico garantiza que \(\mathfrak{l}(\boldsymbol{\beta})\) es estrictamente convexa.
La estimación de los parámetros de un MLG se hace vía máxima verosimilitud, en R:
#Recordemos que al igual que con el modelo de regresion lineal y en general antes de aplicar cualquier metodo
#estadistico, es necesario hacer un analisis descriptivo de las variables a incluir.
#El siguiente es el ejemplo 1.12.2 Proceso infeccioso pulmonar, de Paula (2013)
library(readxl)
#El archivo contiene los datos de 175 pacientes con proceso infeccioso pulmonar, con las siguientes
#variables:
#tipo de tumor (1: maligno, 0: benigno); idade (edad en anos); sexo (1: masculino, 2: femenino);
#hl, intensidad de la celula histiocitos-linfocitos (1:ausente, 2: discreta, 3: moderada, 4: intensa)
#ff, intensidad de la celula fibrosis-floja (1: ausente, 2: discreta, 3: moderada, 4: intensa)
canc3.dat <- read_excel("canc3.xlsx")
head(canc3.dat)## # A tibble: 6 x 5
## tipo idade sexo hl ff
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 26 1 3 1
## 2 0 72 2 4 3
## 3 0 26 2 3 3
## 4 0 27 1 4 1
## 5 0 82 2 3 2
## 6 0 20 2 2 1attach(canc3.dat)
#se deben definir las variables como factores, pues de lo contrario cuando el modelo las tomara como
#variables cuantitativas
sexo = factor(sexo)
hl = factor(hl)
ff = factor(ff)
#El ajuste del modelo se hace con la funcion glm, enla opcion family se especifica la distribucion
#de la variable respuesta, en la opcion link se indica la funcion de enlace a utilizar, en caso de
#no hacerlo, R ajusta el modelo con la funcion de enlace canonica:
fit1.canc3 = glm( tipo ~ sexo + idade + hl + ff, family=binomial)
#El siguiente objeto nos muestra el resultado del ajuste,en la tabla tenemos el resultados de las
#estimaciones de los parametros del modelo (Estimate), su error estandar (Std. Error), el estadistico de prueba calculado
#(z value) y el valor p para la prueba de Wald (Pr(>|z|))
summary(fit1.canc3)##
## Call:
## glm(formula = tipo ~ sexo + idade + hl + ff, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4487 -0.6713 -0.2621 0.7367 1.9773
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.85044 1.06250 -1.742 0.0816 .
## sexo2 0.78392 0.47066 1.666 0.0958 .
## idade 0.06513 0.01328 4.906 9.31e-07 ***
## hl2 -0.86926 0.94750 -0.917 0.3589
## hl3 -2.24903 0.97178 -2.314 0.0206 *
## hl4 -3.29561 1.47997 -2.227 0.0260 *
## ff2 -0.68718 0.50397 -1.364 0.1727
## ff3 -1.02472 0.52721 -1.944 0.0519 .
## ff4 0.43149 1.12452 0.384 0.7012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 236.34 on 174 degrees of freedom
## Residual deviance: 157.40 on 166 degrees of freedom
## AIC: 175.4
##
## Number of Fisher Scoring iterations: 5#Note que como en el caso de regresion lineal multiple, el paquete construye automaticamente las
#variables dummies para los factores, dejando por defecto como grupo de comparacion el nivel mas
#bajo (el que no se muestra en los resultados)
#Estimacion por intervalo
confint(fit1.canc3)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -3.91741623 0.398309188
## sexo2 -0.12563188 1.732478772
## idade 0.04059474 0.093048020
## hl2 -2.99077007 0.867905944
## hl3 -4.41300698 -0.467768524
## hl4 -6.77085234 -0.687364335
## ff2 -1.68733707 0.301717484
## ff3 -2.08753568 -0.005401418
## ff4 -1.90183312 2.747609690Supongamos: \[\log \mathfrak{l}(\boldsymbol{\beta})=L(\boldsymbol{\mu},\boldsymbol{y})=\sum_{i=1}^{n}L(\mu_i,y_i)\] con \(\mu_i=g^{-1}(\eta_i)\), \(\eta_i=x_i^t\boldsymbol{\beta}\).
Para el modelo saturado \((p=n)\), \(L(\mu,y)\) es estimado por: \[L(\boldsymbol{y},\boldsymbol{y})=\sum_{i=1}^{n}L(y_i,y_i)\] Sea \(\hat{\mu}_i\) la estimación MV de \(\mu_i\), la calidad del ajuste se mide por:
\[D^*(\boldsymbol{y},\boldsymbol{\mu})=\phi D(\boldsymbol{y},\boldsymbol{\mu})=2\left\lbrace L(\boldsymbol{y},\boldsymbol{y})-L(\boldsymbol{\mu},\boldsymbol{y}) \right\rbrace\]
Siendo \(D(\boldsymbol{y},\boldsymbol{\mu})\) la función de desvío.
Sean: \(\hat{\theta_i}=\theta_i(\hat{\mu_i})\) la estimación para el modelo con \(p\) parámetros y \(\tilde{\theta_i}=\theta_i(\tilde{\mu_i})\) la estimación para el modelo saturado, entonces:
\[D(\boldsymbol{y},\boldsymbol{\mu})=2\sum_{i=1}^{n}\left\lbrace y_i(\tilde{\theta_i}-\hat{\theta_i})+b(\hat{\theta_i})-b(\tilde{\theta_i}) \right\rbrace\]
\[D(\boldsymbol{y},\boldsymbol{\mu})=\sum_{i=1}^{n}d^2(y_i,\mu_i)\]
\(\theta_i=\mu_i\), \(b(\theta_i)=\frac{\theta_i^2}{2}\) \[\Longrightarrow D(\boldsymbol{y},\boldsymbol{\mu})=\sum_{i=1}^{n}\left(y_i-\hat{\mu_i} \right)^2=SCE\]
\(\theta_i=\log\mu_i\), \(b(\theta_i)=\mu_i\) \[\Longrightarrow d^2(y_i,\mu_i)= \begin{cases} 2\left(y_i\log\frac{y_i}{\mu_i}-\left( y_i-\hat{\mu_i}\right) \right) & y_i>0\\ 2\hat{\mu_i} & y_i=0 \end{cases}\]
\[d^2(y_i,\mu_i)= \begin{cases} 2\left\lbrace y_i\log\frac{y_i}{n_i\hat{\mu_i}}+(n_i-y_i)\log\frac{1-y_i/n_i}{1-\hat{\mu_i}} \right\rbrace& 0<y_i<n_i\\ -2n_i\log(1-\hat{\mu_i}) & y_i=0\\ -2n_i\log\hat{\mu_i} & y_i=n_i \end{cases}\]
Supongamos \(\boldsymbol{\beta}=(\boldsymbol{\beta}_1^t,\boldsymbol{\beta}_2^t)^t\), con \({\boldsymbol{\beta}_1^t}_{(q\times 1)}\), \({\boldsymbol{\beta}_2^t}_{((p-q)\times 1)}\). Se quiere juzgar:
\[H_0: \boldsymbol{\beta}_1=0\qquad H_1: H_0: \boldsymbol{\beta}_1\neq 0\]
Sean:
\(D(\boldsymbol{y},\hat{\boldsymbol{\mu}^0})\) el desvío con la estimación bajo \(H_0\).
\(D(\boldsymbol{y},\hat{\boldsymbol{\mu}})\) el desvío con la estimación MV.
\[\xi_{RV}=\phi\left\lbrace D(\boldsymbol{y},\hat{\boldsymbol{\mu}^0})-D(\boldsymbol{y},\hat{\boldsymbol{\mu}}) \right\rbrace\longrightarrow^d \chi^2_{(q)} \]
bajo \(H_0\), cuando \(\phi\) es conocido.
Modelo logístico para la ocurrencia de cáncer de pulmón en pacientes con un proceso infeccioso pulmonar, con parte sistemática:
\[\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_3+\beta_4x_4\]
dónde: \(x_1=1\) si es mujer, \(x_1=0\) si es hombre, \(x_2\):Edad en años, \(x_3\) y \(x_4\) son factores con 4 niveles que representan la intensidad de dos tipos de célula.
| Modelo | Desvío | Dif | gl | Prueba | \(\chi^2_{0.95}\) | Valor p |
|---|---|---|---|---|---|---|
| Constante | 236.34 | |||||
| \(+x_1\) | 235.20 | 1.14 | 1 | \(x_1\) | 3.84 | 0.286 |
| \(+x_2\) | 188.22 | 46.98 | 1 | \(x_2|x_1\) | 3.84 | \(<0.001\) |
| \(+x_3\) | 162.55 | 25.67 | 3 | \(x_3|x_1+x_2\) | 7.81 | \(<0.001\) |
| \(+x_4\) | 157.40 | 5.15 | 3 | \(x_4|x_1+x_2+x_3\) | 7.81 | 0.161 |
McCullagh,P.; Nelder, J. (1989). Generalized Linear Models. Chapman Hall, London.
Paula, G. (2013). Modelos de Regress~{a}o com apoio computacional. Instituto de Matemática e Estatística. Universidade de S~{a}o Paulo, Brasil.
Shao, J. (2007). Mathematical Statistics. Springer. USA.