---
title: "Respuesta conteos"
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
---  

# Regresión Poisson  

## La distribución Poisson 
- La función de probabilidad de una distribución de Poisson es
$$
f(y ; \mu) = \frac{e^{-\mu} \mu^y}{y!},
$$
con $y = 0, 1, \ldots $.
- La podemos expresar como un elemento de la familia exponencial natural  
$$
f(y; \mu) = \frac{e^{-\mu}\mu^{y}}{y!} =  \exp\{y \ln \mu - \mu - \ln(y!)\},
$$
tomando $\theta_i = \ln \mu_i$, $b(\theta_i) = \exp{\theta_i}$,
  $a(\phi) =1$, $c(y_i,\phi) = - \ln(y_i!)$. 


## Ecuaciones de verosimilitud  
- Son las siguientes  
$$
\sum_{i=1}^n  (y_i - \mu_i) x_{ij} = 0.
$$
- La desviación (escalada y sin escalar) y con término constante es 
$$
D(\mathbf y; \hat{\mathbf \mu}) =  2 \sum_{i=1}^n y_i \ln \frac{y_i}{\hat{\mu}_i}.
$$
- Adopta la expresión
$$
D(\mathbf y; \hat{\mathbf \mu}) = 2 \sum Observado \ln \left (\frac{Observado}{Esperado}
\right ) .
$$

# `tamidata2:PRJNA218851` 

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

```{r}
table(colData(PRJNA218851)[,"Stage"])
```
- Nos fijamos en el gen que ocupa la fila 1000.
- Construimos un `data.frame` en donde incluimos la variable
`Stage` y el conteo correspondiente a este gen. 

```{r}
df = data.frame(count = assay(PRJNA218851)[1000,],
                Stage=colData(PRJNA218851)[,"Stage"])
head(df)
```
## Regresión Poisson  

- Ajustamos un modelo loglineal de Poisson. 
```{r}
fit = glm(count ~ Stage, family = poisson(link = log),data = df)
summary(fit)
```
- No es adecuado por varias razones:
    - Distintas profundidades de secuenciación. 
	- Posiblemente sobre dispersión.  
	- No lo podemos usar.  

- Vemos la diferencia de las desviaciones.  
```{r}
fit$null.deviance - fit$deviance
```

## Sobre dispersión  

- Vamos a estimar la sobre dispersión.  

```{r}
fit1 = glm(count ~ Stage, family = quasipoisson(link = log),data = df)
summary(fit1)$dispersion
```


# Regresión binomial negativa  

## Distribución binomial negativa  
- La función de probabilidad es 
$$
f(y; \phi, \mu) = \frac{\Gamma(y+\phi)}{\Gamma(\phi)\Gamma(y+1)} \bigg (
\frac{\phi}{\mu+\phi} \bigg)^\phi \bigg ( 1-  \frac{\phi}{\mu+\phi} \bigg)^y
$$
con $y= 0, 1, 2 , \ldots$ donde $\phi$ y $\mu$ son los parámetros.  

- Se tiene que 
$$
E(Y) = \mu,
$$
$$
var(Y) = \mu + \frac{\mu^2}{\phi}.
$$
- El parámetro $1/\phi$ es un **parámetro de dispersión**. 


## Ejemplo  

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