---
title: "Modelos lineales generalizados"
author: "Guillermo Ayala Gallego"  
date: "`r Sys.Date()`"  
address: Departamento de Estadística e I.O. Universidad de Valencia  
format: 
   revealjs:  
      scrollable: true
      echo: true
      embed-resources: true  
   pdf:
       toc: true
       number-sections: false
       colorlinks: true
---

# Modelo lineal generalizado  

## Componentes de  un modelo lineal generalizado  

- **Componente aleatoria**
- **Componente sistemática**
- **Función de enlace**

##   La componente aleatoria  

-  Consiste de una variable aleatoria $Y$ con observaciones
independientes $(Y_1,\ldots,Y_n)$.

- La respuesta $Y$ sigue una distribución en la **familia de dispersión 
exponencial** 
$$
f(y_i; \theta_i,\phi) = e^{\frac{y_i \theta_i - b(\theta_i)}{a(\phi)} + 
c(y_i,\phi)}.
$$

- El parámetro $\theta_i$ el  **parámetro natural**.

- El parámetro $\phi$ es el **parámetro de dispersión**.

-  Familia exponencial natural ss un caso particular de la 
**familia de dispersión exponencial**.
   
- Ocurre cuando
	$$
		a(\phi) = 1, 
		$$
		$$
c(y_i,\phi) = c(y_i).
$$

- La densidad tiene la forma
$$
f(y_i; \theta_i) = h(y_i) \exp{[y_i \theta_i - b(\theta_i)]}.
$$

- Verosimilitud de una sola observación
$$
f(y_i; \theta_i,\phi) = e^{\frac{y_i \theta_i - b(\theta_i)}{a(\phi)} + 
	c(y_i,\phi)}.
$$
- Logverosimilitud de una sola observación
$$
\ell_i = \frac{y_i \theta_i - b(\theta_i)}{a(\phi)} + c(y_i,\phi).
$$

- Logverosimitud 
$$
\ell = \sum_{i=1}^n \ell_i.
$$

## Distribución binomial  
  
- $Y_i$ es la **proporción muestral   de éxitos**.

- $n_i Y_i \sim Bi(n_i,\pi_i)$.

-  $\mu_i = EY_i = \pi_i$.

- Consideramos 
    $$
    \theta_i = \log \frac{\pi_i}{1-\pi_i}
    $$
    esto es, definimos $\theta_i$ como el logit de la probabilidad de éxito.
 
 - La transformación inversa es 
    $$
    \pi_i = \frac{e^{\theta_i}}{1+e^{\theta_i}},
    $$
y que $\log(1-\pi_i) = - \log\bigg (1+e^{\theta_i}\bigg )$.
     
- 
 $$
 f(y_i; \pi_i,n_i) = \binom{n_i}{n_i y_i} \pi_i^{n_i y_i}
      (1-\pi_i)^{n_i - n_i y_i} = 
$$
$$
\exp \bigg [ \frac{y_i \theta_i - \log[1+\exp(\theta_i)]}{1/n_i} 
      + \log \binom{n_i}{n_i y_i} \bigg ],
$$
    siendo
$$
      b(\theta_i) = \log [1+\exp(\theta_i)], \  a(\phi) = 1/n_i,
    \ c(y_i,\phi) = \log \binom{n_i}{n_i y_i}.
$$

-   El parámetro natural es
$$
    \theta_i =  \log \frac{\pi_i}{1-\pi_i},
$$
el logit de $\pi_i$. 


## Distribución Poisson

- La función de probabilidad es 
 $$
    f(y_i; \mu_i) = \frac{e^{-\mu_i}\mu_i^{y_i}}{y_i!} = 
    \exp (y_i \log \mu_i - \mu_i - \log(y_i!)).
    $$
	
 -  Tomamos
    $$\theta_i = \log \mu_i,$$
	$$b(\theta_i) = \exp{\theta_i},$$
    $$a(\phi) =1,$$
	$$c(y_i,\phi) = - \log(y_i!).$$


##  La componente sistemática  
  
- La componente sistemática de un modelo lineal generalizado
    es el vector $(\eta_1,\ldots,\eta_n)$.  
	
-  Cada uno  de los $\eta_i$ es la combinación lineal de los predictores
    correspondientes a la $i$-ésima observación, es decir,
    $$
    \eta_i = \sum_{j=1}^p \beta_j x_{ij},
    $$
    con  $i=1, \ldots, n$ y $x_{ij}$ es $j$-ésimo predictor 
    en el $i$-ésimo individuo.
  
- La combinación lineal $\sum_j \beta_j x_{ij}$  es el
    **predictor lineal**.
 
- Se suele asumir que $x_{i1} = 1$. 
 
## Función de enlace   
  
- Mediante la **función de enlace** $g$ relacionamos las
    componentes aleatoria y sistemática 
    $$
    g(\mu_i) = \eta_i = \sum_{j=1}^p \beta_j x_{ij}.
    $$
  
  - La función de enlace que nos transforma la media $\mu_i$ en el parámetro
    natural recibe el nombre de **enlace canónico**
    $$
    \theta_i = \sum_{j=1}^p \beta_j x_{ij}.
    $$

## Estimación de los coeficientes  

- La función de verosimilitud es
$$
L(\mathbf  \beta) = \prod_{i=1}^n
\exp\bigg \{\frac{y_i \theta_i - b(\theta_i)}{a(\phi)}+ c(y_i,\phi)
\bigg \}
$$
y la función de logverosimilitud es 
$$
\ell(\mathbf \beta) = \sum_{i=1}^n \frac{y_i \theta_i - b(\theta_i)}{a(\phi)}+ 
\sum_{i=1}^n  c(y_i,\phi).
$$

- Las ecuaciones de estimación son
    $$
    \frac{\partial \ell(\mathbf \beta)}{\partial \beta_j} = 
    \sum_{i=1}^n \frac{\partial \ell_i(\mathbf \beta)}{\partial \beta_j}
    $$
y tienen la expresión 
		
$$
\frac{\partial \ell(\mathbf \beta)}{\partial \beta_j} =
\sum_{i=1}^n \frac{(y_i  - \mu_i) x_{ij}}{var(Y_i)} 
    \frac{\partial \mu_i}{\partial \eta_i} = 0,
$$
con $j=1, \ldots, p$ siendo
$$
\eta_i = \sum_{j=1}^p \beta_j x_{ij} = g(\mu_i)
$$
para la función de enlace $g$.

- Estas ecuaciones han de ser resueltas de un modo iterativo. 

- La distribución de la variable $Y_i$ aparece en las ecuaciones
    de verosimilitud solamente a través de su media $\mu_i$ y su 
    varianza $var(Y_i)$. 
  
## Distribución de los estimadores máximo verosímiles  
  
- La distribución asintótica de los coeficientes $\hat{\mathbf \beta}$
  viene dada por 
  $$
  \hat{\mathbf \beta} \sim N_p(\mathbf \beta, (\mathbf X^T \mathbf W \mathbf X)^{-1}),
  $$
  siendo $\mathbf W$ una matriz diagonal con las entradas 
  $$
  w_i =  \bigg( \frac{\partial \mu_i}{\partial \eta_i} \bigg )^2
  \frac{1}{var(Y_i)}.
  $$

##  Covarianzas de $\hat{\mathbf \beta}$  
  
- La matriz de covarianzas asintótica de $\hat{\mathbf \beta}$ viene 
  dada por
$$
\widehat{var(\hat{\mathbf \beta)}} = (\mathbf X^T \hat{\mathbf W} \mathbf X)^{-1}
$$
donde $\hat{\mathbf W}$ es la matriz $\mathbf W$ evaluada en $\hat{\mathbf \beta}$.

## Comparando modelos anidados  

- Podemos comparar modelos anidados.
- $\mathbf X_0$ y $\mathbf X_1$ matrices de modelo tales que
  $C(\mathbf X_0) \subset C(\mathbf X_1)$.
- El contraste de hipótesis que nos interesa es 
\begin{align*}  \label{align:110123a}
  H_0: \mathbf g(\mathbf \mu ) & = \mathbf X_0 \mathbf \beta_0,\\
  H_1:  \mathbf g(\mathbf \mu ) &  = \mathbf X_1 \mathbf \beta_1
\end{align*}
siendo $\mathbf g(\mathbf \mu) = (g(\mu_1),\ldots,g(\mu_p))^T$.
- El test del cociente
de verosimilitudes se basa en comparar los máximos que podemos alcanzar
en la logverosimilitud bajo cada una de las hipótesis,
$\ell(\hat{\mathbf \beta_0})$ y $\ell(\hat{\mathbf \beta_1})$, y utilizar
el resultado que afirma 
$$
  -2(\ell(\hat{\mathbf \beta_0}) - \ell(\hat{\mathbf \beta_1})) \sim \chi^2_{p_1-p_0},
$$
siendo $p_0$ y $p_1$ las dimensiones de los espacios columna $C(\mathbf X_0)$ y
$C(\mathbf X_1)$. 
 
 
# Regresión logística  

## Datos 

- Son datos anonimizados. 

```{r}
#| label: GLM1
library(tamidata3)
finput = system.file("extdata","SNPs_RM.csv",package="tamidata3")
df = read.table(file = finput,header = TRUE,sep=",")
```
- Realizamos alguna modificación. 

```{r}
#| label: GLM2
df$clinica.diabetes.mellitus.type2 =  (df$clinica.diabetes.mellitus.type2 == "YES")
df$clinica.gender = factor(df$clinica.gender)
```
- Podemos ver las primeras filas de los datos.  
```{r}
#| label: GLM3
head(df,n=2)
```
- Los SNPs aparecen en las primeras columnas. 
- Hemos de transformar estas variables para convertirlas en posibles predictores de acuerdo
con distintos modelos genéticos.

```{r}
#| label: GLM4

#' Transformation of the SNP to a genetic model
#' @description
#' Transformation of the SNP to a genetic model
#' @param x SNPs
#' @param type Model to be used
#' @param sep Separator
#' @export
snp2model = function(x,type=c("codominant","dominant","recessive"),
                     sep=""){
  x1 = substr(x,1,1)
  x2 = substr(x,2,2)
  a = table(c(x1,x2))
  recessive  = names(which.min(a))
  dominant = names(which.max(a))
  x1 = (x1 == recessive)*1
  x2 = (x2 == recessive)*1
  if(type == "codominant") rs =  factor(x1+x2)
  if(type == "dominant") rs = factor(x1+x2 == 0) 
  if(type == "recessive") rs = factor(x1+x2 == 2)
  rs
}
```
- Vamos a evaluar cada uno de los modelos genéticos buscando asociación con la respuesta.
```{r}
#| label: GLM5
snpscol = 1:5
df1 = df 
df1[,snpscol] = apply(df1[,snpscol],2,snp2model,type="codominant")
df1$SNP1 = factor(df1$SNP1)
df1$SNP2 = factor(df1$SNP2)
df1$SNP3 = factor(df1$SNP3)
df1$SNP4 = factor(df1$SNP4)
df1$SNP5 = factor(df1$SNP5)
head(df1)
```
Ajustamos el modelo. 

```{r}
#| label: GLM6
fit = glm(clinica.diabetes.mellitus.type2 =="YES" ~ ., family="binomial",data=df1)
```
- Vemos un resumen del ajuste. 
```{r}
#| label: GLM7
summary(fit)
```

- Utilizamos el modelo recesivo.  
```{r}
#| label: GLM8
snpscol = 1:5
df1 = df 
df1[,snpscol] = apply(df1[,snpscol],2,snp2model,type="recessive")
fit = glm(clinica.diabetes.mellitus.type2 =="YES" ~ ., family="binomial",data=df1)
summary(fit)
```
 
# Modelo loglineales de Poisson 

## Datos  
```{r}
#| label: GLM9
pacman::p_load(SummarizedExperiment)
data(PRJNA218851,package="tamidata2")  
```

```{r}
#| label: GLM10
table(colData(PRJNA218851)[,"Stage"])
```

```{r}
#| label: GLM11
df = data.frame(count = assay(PRJNA218851)[1001,],
                Stage=colData(PRJNA218851)[,"Stage"])
head(df)
```

## Modelo loglineal de Poisson  

Ajustamos un modelo loglineal de Poisson.
```{r}
#| label: GLM12
fit = glm(count ~ Stage, family = poisson(link = log),data = df)
summary(fit)
```

- La desviacion nula es la desviacion para el modelo que tiene solo la 
constante.  
- La desviacion residual es la desviacion del modelo que tiene 
la constante y las variables binarias que describen `Stage`.   
- La diferencia entre los valores tiene una distribucion ji-cuadrado con 
dos grados de libertad y nos permite 
contrastar si los coeficientes de  `StageMetastasis`
y `StageNormal` pueden considerarse simultáneamente nulos. 
```{r}
#| label: GLM13
fit$null.deviance - fit$deviance
```
Podemos rechazar confortablemente la hipotesis nula. 

## Sobredispersión  

- En una distribución de Poisson, la media y la varianza son iguales. 
- Cuando trabajamos con conteos reales no suele ser cierta esta hipótesis. 
- Con frecuencia la varianza es mayor que la media. 
- A esto se le llama **sobredispersión**. 
```{r}
#| label: GLM14
fit1 = glm(count ~ Stage, family = quasipoisson(link = log),data = df)
summary(fit1)$dispersion
```

## GLM binomiales negativos  

- La distribución binomial negativa no está en la familia de dispersión exponencial salvo que 
consideremos conocido el parámetro de dispersión.  
- Hemos de estimarlo previamente.  

```{r}
#| label: GLM15
library(MASS)
fit = glm.nb(count~Stage,data=df)
summary(fit)
```

- Podemos ver que los errores estándar de los coeficientes son mucho mayores porque
tenemos en cuenta la sobre dispersión que el modelo de Poisson no lo considera. 
