00-biostat2019-slides.Rmd

---
title: "Estadística Inferencial en R"
subtitle: "Día 1 del curso:<br/>Análisis de RNA-seq en R"
author: "Andree Valle Campos<br/> @avallecam"
institute: "CDC - Perú<br/>Centro Nacional de Epidemiología,<br/>Prevención y Control de Enfermedades"
date: "2019-02-15"
output:
  xaringan::moon_reader:
    lib_dir: libs
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE, message=FALSE}
library(tidyverse)
library(janitor)
library(corrr)
library(broom)
library(haven)
# library(ggdag)
library(rlang)
library(skimr)
library(compareGroups)
library(knitr)
library(kableExtra)
library(qvalue)
# library(survival)
# library(survminer)
# library(epiR)

set.seed(33)

knitr::opts_chunk$set(echo = TRUE,message=FALSE,warning=FALSE, error=FALSE)
theme_set(theme_bw())

options(htmltools.dir.version = FALSE)

aidsdb <- read_dta("data-raw/aidsdb.dta") %>% as_factor() %>% 
  mutate_if(.predicate = is.factor,.funs = fct_relevel,... = "No")
espir <- read_dta("data-raw/espirometria.dta") %>% as_factor()
multpv <- read_rds("data-raw/createpval.rds")

# to solve skimr histograms charts
# https://docs.ropensci.org/skimr/index.html#support-for-spark-histograms
# fix_windows_histograms()
Sys.setlocale("LC_CTYPE", "Chinese")
```

# Temario

1. __(breve) Introduction a R__

  * Conceptos clave

2. __Variables:__ 

  * Concepto y clasificación
  * Explorar distribuciones
  * __Prueba de Hipótesis:__ limitantes y alternativas

5. __Modelos lineales:__ 

  * __Regresión lineal y logística__
  * Relación con PH
  * __Comparación múltiple:__ corrección del valor p y FDR
  * Aplicación en _microarrays_

---

# Metodología

- Teórico y __práctico__

--

  * Definiciones y procedimientos
  
  * Responder a: _¿y cómo lo hago en R?_

--

- Material

  * 2 archivos de teoría en PDF
  
  * 5 practicas + solucionario 

--

- Contenido disponible en __https://avallecam.github.io/biostat2019/__

---

class: center, middle

# ¡Pregunta 1!

## Del 1 al 5, ¿Qué tan familiarizado estás con cada tema?

---
class: inverse, center, middle

# (breve) introducción a R

---

class: center, middle

__Ir a__ 

__`parte 02`__

[click aquí](https://speakerdeck.com/avallecam/01-biostat2019-slides)

(luego retornar)

---

# Práctica 1

```{r,eval=FALSE}
#aritmética
2+2
x <- 2
x+2

#funciones
seq(1,10)
rep(1,5)

#vectores
heights <- c(147.2, 153.5, 152.5, 162.0, 153.9)
mean(heights)

#subconjuntos
y <- LETTERS[1:10]
y[3:7]
```

Si necesitas ayuda:
  1. Pregúntale a R con la función `help()` o `?`, p.e.: `help(mean)` o `?mean`
  2. Nos pasas la voz :)

---
class: inverse, center, middle

# Variables

<!--
- Enfatizar en los supuestos
- Probabilidades, distribución normal, definición del estadístico z, valor p, 
- ¿teorema del límite central? -> error estandar -> intervalo de confianza
- distribución T student
-->

---

# Variables: Ejemplo

**VEF = Volumen Expiratorio Forzado (mL/s)**

**¿Existe una asociación entre VEF y edad, talla, etc?**

**¿Cómo se ve afectado el VEF en las personas que fuman?**

- base: `espirometria.dta`
- exposición: `edad`, `talla`
- desenlace: `vef`

```{r,echo=FALSE}
knitr::kable(head(espir), format = 'html')
```

---

# Variables numéricas

```{r,eval=FALSE}
espir %>% 
  select(edad,vef,talla) %>% 
  skim() #<<
```

```{r,echo=FALSE,out.width = '50%'}
# render = knitr::normal_print
# fix_windows_histograms()
# Sys.setlocale("LC_CTYPE", "Chinese")
espir %>% 
  select(edad,vef,talla) %>% 
  skim() %>% 
  yank(skim_type = "numeric") %>% 
  select(-complete_rate) %>% 
  kable(format = "html") %>% 
  kable_styling(full_width = F,font_size = 15)
```

---

# Variables numéricas: PH

.pull-left[

```{r,fig.height=4,fig.width=4}
espir %>% 
  ggplot(aes(x=talla,y=vef)) + #<<
  geom_point() +
  geom_smooth(method = "lm")
```

]

--

.pull-right[

```{r,fig.height=4,fig.width=4}
espir %>% 
  select(edad,vef,talla) %>% 
  correlate() %>%  #<<
  rearrange() %>% 
  shave()
```

]

---

# Variables categóricas

```{r}
espir %>% 
  tabyl(sexo) %>%  #<<
  adorn_totals("row") %>% 
  adorn_pct_formatting()
```

--


```{r}
espir %>% 
  tabyl(sexo,fumar) %>%  #<<
  adorn_totals(c("col")) %>% 
  adorn_percentages() %>% 
  adorn_pct_formatting() %>% 
  adorn_ns(position = "front") %>% 
  adorn_title()
```

---

# Variables categóricas: PH

```{r}
espir %>% 
  tabyl(sexo,fumar) %>% 
  chisq.test() %>%  #<<
  tidy()
```

---

# Variables numéricas y categóricas

.pull-left[

```{r,fig.height=4,fig.width=4}
espir %>% 
  ggplot(aes(x=fumar,y=vef)) +
  geom_point(
    aes(color=fumar), #<<
    position = "jitter") + #<<
  geom_boxplot(alpha=0) #<<
```

]

--

.pull-right[

```{r,fig.height=4,fig.width=4}
espir %>% 
  ggplot(aes(x=fumar,y=vef)) +
  geom_violin(aes(color=fumar)) #<<
```

]

---

# Variables numéricas y categóricas: PH

```{r,eval=FALSE}
espir %>% 
  select(fumar,vef) %>% 
  group_by(fumar) %>%  #<<
  skim()
```

```{r,echo=FALSE}
espir %>% 
  select(fumar,vef) %>% 
  group_by(fumar) %>%  
  skim() %>%
  yank(skim_type = "numeric") %>% 
  select(-complete_rate) %>% 
  kable(format = "html") %>% 
  kable_styling(full_width = F,font_size = 15)
```

--

```{r,eval=FALSE}
t.test(vef ~ fumar, data = espir, var.equal=FALSE) %>% 
  tidy() #<<
```

```{r,echo=FALSE}
t.test(vef ~ fumar, data = espir, var.equal=FALSE) %>% 
  tidy() %>% 
  kable(format = "html",digits = 2) %>% 
  kable_styling(full_width = F,font_size = 15)
```

---

# Práctica 2

```{r,eval=FALSE}
#importar base en formato DTA
espir <- read_dta("data-raw/espirometria.dta") %>% as_factor()

#generar resumen descriptivo
espir %>% skim()

#grafica distribución
espir %>% 
  ggplot(aes(vef)) +
  geom_histogram()

#realizar pruebas de hipótesis
#experimenta y describe qué cambios genera cada uno de los argumentos?
compareGroups(formula = fumar ~ edad + vef + talla + sexo,
              data = espir, 
              byrow=T,
              method=c(vef=2)
              ) %>% 
  createTable(show.all = T) %>% 
  export2xls("table/tab1.xls")
```

```{r,eval=FALSE,echo=FALSE}

espir %>% 
  select(edad,vef,talla) %>% 
  correlate() %>% 
  stretch()

espir %>% 
  select(edad,vef,talla) %>% 
  as.matrix() %>% 
  Hmisc::rcorr() %>% 
  broom::tidy()
```

```{r,eval=FALSE,echo=FALSE}
#regresar luego de finalizado el capítulo de regresiones
t.test(vef ~ fumar, data = espir,var.equal=TRUE) %>% 
  tidy()
lm(vef ~ 1+ fumar, data = espir) %>% tidy()
lm(vef ~ 1+ fumar, data = espir) %>% tidy(conf.int=TRUE)
```

---

# Comentario: Pruebas de hipótesis

- __Limitante:__ 
  - dicotomización de resultados en significativos y no significativos
  - valor p da indicativo de _tamaño de efecto_, pero no lo cuantifica

- __Alternativa<sup>1</sup>:__

```{r, echo=FALSE, fig.cap="", out.width = '100%'}
knitr::include_graphics("notrack/estimation_plot.PNG")
```

.footnote[
[1] [Moving beyond P values: data analysis with estimation graphics](https://www.nature.com/articles/s41592-019-0470-3.epdf?author_access_token=Euy6APITxsYA3huBKOFBvNRgN0jAjWel9jnR3ZoTv0Pr6zJiJ3AA5aH4989gOJS_dajtNr1Wt17D0fh-t4GFcvqwMYN03qb8C33na_UrCUcGrt-Z0J9aPL6TPSbOxIC-pbHWKUDo2XsUOr3hQmlRew%3D%3D)
]

---
class: inverse, center, middle

# Modelos lineales


---

# ¿**Por qué** usar una **regresión**?

--

```{r, echo=FALSE, fig.cap="", out.width = '90%',fig.align='center'}
knitr::include_graphics("notrack/nmeth.3627-F1.jpg")
```

--


- Porque:

  - Fija una **variable independiente** o **exposición** y observa una respuesta en la **variable dependiente** o **desenlace**.

  - Permite **explicar** el cambio promedio de un evento $Y$ en base a cambios en $X$, usando coeficientes o medidas de asociación.

  - Permite **predecir** la probabilidad asociada a un evento.
  
  - Permite **cuantificar** el tamaño del efecto de la comparación.


---

# Regresión Lineal Simple

.pull-left[

## características

- **una** variable independiente (simple)
- **una** variable dependiente (univariada)
- ambas variables deben ser **numéricas**

## objetivos

- ajustar datos a una recta
- interpretar medida de bondad de ajuste ( $R^2$ ) y coeficientes
- evaluar supuestos
- visualizar el modelo


]


.pull-right[

```{r, echo=FALSE,fig.height=6,fig.width=6}
women %>% ggplot(aes(weight,height)) + geom_point() + geom_smooth(method = "lm")
```

$$Y = \beta_0 + \beta_1 X_1 + \epsilon$$

]

---

# RLinS: Ejemplo 

.pull-left[

```{r,fig.height=5.5,fig.width=5.5}
espir %>% 
  ggplot(aes(x=edad,y=vef)) + #<<
  geom_point()
```

]

--

.pull-right[

- __Pregunta__

  + **¿Existe una relación lineal entre edad y VEF?**
  
  + **¿En cuánto incrementa el VEF, por cada incremento en un año de edad?**

]

---

# RLinS: Ejemplo 

.pull-left[

```{r,fig.height=5.5,fig.width=5.5}
espir %>% 
  ggplot(aes(x=edad,y=vef)) +
  geom_point()
```

]

.pull-right[

- __Pregunta__

  + **¿Existe una relación lineal entre edad y VEF?**
  
  + **¿En cuánto incrementa el VEF, por cada incremento en un año de edad?**

- __Evaluación de supuestos:__

  + **#1** independencia de observaciones
  
  + **#2** linealidad

]


---

# RLinS: residuales

- Definición: Diferencia entre el _valor observado_ y el _valor predicho_ en el eje vertical.

.pull-left[

```{r, echo=FALSE, fig.cap="", out.width = '60%',fig.align='center'}
knitr::include_graphics("notrack/Inkedlss01_LI.jpg")
```

```{r, echo=FALSE, fig.cap="", out.width = '60%',fig.align='center'}
knitr::include_graphics("notrack/Inkedlss02_LI.jpg")
```

]

--

.pull-right[

```{r, echo=FALSE, fig.cap="", out.width = '60%',fig.align='center'}
knitr::include_graphics("notrack/Inkedlss02_LI.jpg")
```

```{r, echo=FALSE, fig.cap="", out.width = '60%',fig.align='center'}
knitr::include_graphics("notrack/Inkedlss03_LI.jpg")
```

]

---

# RLinS: suma de mínimos cuadrados

.pull-left[

<!--![fuente: `https://youtu.be/nk2CQITm_eo`](notrack/lss.png)-->

```{r, echo=FALSE, fig.cap="", out.width = '90%'}
knitr::include_graphics("notrack/lss.png")
```

]

--

.pull-right[

Cálculo de la **sumatoria del cuadrado de los residuales** hacia la media y la recta:

$$ SSE(mean) = \sum (data - mean)^2 $$
$$ SSE(fit) = \sum (data - fit)^2 $$

]

---

# RLinS: suma de mínimos cuadrados

.pull-left[

<!--![fuente: `https://youtu.be/nk2CQITm_eo`](notrack/lss.png)-->

```{r, echo=FALSE, fig.cap="", out.width = '90%'}
knitr::include_graphics("notrack/lss.png")
```

]

.pull-right[

Cálculo de la **sumatoria del cuadrado de los residuales** hacia la media y la recta:

$$ SSE(mean) = \sum (data - mean)^2 $$
$$ SSE(fit) = \sum (data - fit)^2 $$

$$ Var(x) = \frac{SSE(x)^2}{n} $$

Medida de **bondad de ajuste**:

$$ R^2 = \frac{Var(mean) - Var(fit)}{Var(mean)}  $$

]

---

# RLinS: bondad de ajuste

```{r,eval=FALSE}
wm1 <- lm(vef ~ edad, data = espir)

wm1 %>% 
  glance() #<<
```

```{r,echo=FALSE}
wm1 <- lm(vef ~ edad, data = espir)
wm1 %>% glance() %>% dplyr::select(r.squared:df)
```

**INTERPRETACIÓN**

--

- **edad** _explica_ el **57%** de la variabilidad de **VEF**

- existe un **57%** _de reducción en_ la variabilidad de **VEF** al tomar en cuenta la **edad**


---

# RLinS: coeficientes

```{r}
wm1 %>% tidy()
```

**INTERPRETACIÓN**

--

```{r, echo=FALSE, fig.cap="", out.width = '90%'}
knitr::include_graphics("notrack/logan-06.PNG")
```

---

# RLinS: coeficientes

```{r}
wm1 %>% tidy()
wm1 %>% tidy(conf.int=TRUE) #<<
```

- $\beta_{edad}$: 
- En la población, por cada incremento de **edad** en _una unidad_, el **VEF** en promedio _incrementa_ en 0.22 mL/s,
- con un intervalo de confianza al 95% de 0.21 a 0.24 mL/s.
- Este resultado es estadísticamente significativo con un valor **p < 0.001**

---

# RLinS: supuesto **#3** normalidad residuales

.pull-left[

```{r, echo=FALSE, fig.cap="", out.width = '75%'}
knitr::include_graphics("notrack/Inkedlss03_LI.jpg")
```
]

.pull-right[

```{r, echo=FALSE, fig.cap="", out.width = '100%'}
knitr::include_graphics("notrack/nmeth.3854-F2.jpg")
```

]

---

# RLinS: supuesto **#3** normalidad residuales

- generación de _dataframe_ con:
  * `.fitted` : valor predicho de Y para valor de X ( $\hat Y$ )
  * `.resid` : valor crudo del residual ( $Y - \hat Y$ )
  * `.std.resid` : valor _estudiantizado_ del residual

```{r,fig.height=4,fig.width=4,eval=TRUE, echo=TRUE}
wm1 %>% 
  augment() #<<
```


---

# RLinS: supuesto **#3** normalidad residuales

.pull-left[

```{r,fig.height=4,fig.width=4,eval=TRUE, echo=TRUE}
wm1 %>% 
  augment() %>% 
  ggplot(aes(.std.resid)) + #<<
  geom_histogram()
```

]

--

.pull-right[

- **META:** puntos sobre la línea

```{r,fig.height=4,fig.width=4}
wm1 %>% 
  augment() %>% 
  ggplot(
    aes(sample=.std.resid) #<<
    ) + 
  geom_qq() + #<<
  geom_qq_line() #<<
```

]

---

# RLinS: supuesto **#4** homoscedasticidad


```{r, echo=FALSE, fig.cap="", out.width = '70%',fig.align='center'}
knitr::include_graphics("notrack/logan-07.PNG")
```

---

# RLinS: supuesto **#4** homoscedasticidad

```{r, echo=FALSE, fig.cap="", out.width = '80%',fig.align='center'}
knitr::include_graphics("notrack/nmeth.3854-F1.jpg")
```

---

# RLinS: supuesto **#4** homoscedasticidad

- **META:** distribución idéntica a ambos lados de la línea

```{r,fig.height=4,fig.width=4,message=FALSE,error=FALSE,fig.align='center'}
wm1 %>% 
  augment() %>% 
  ggplot(aes(.fitted,.std.resid)) +
  geom_point() +
  geom_smooth() + #<<
  geom_hline(yintercept = c(0))
```

---

# RLinS: ¿cómo se ve el modelo?

$$ Y = \beta_0 + \beta_1 X_1 + \epsilon $$
$$ VEF = 0.60 + 0.22 (edad) + \epsilon $$

```{r,fig.height=4,fig.width=4,fig.align='center'}
espir %>% 
  ggplot(aes(edad,vef)) +
  geom_point() + 
  geom_smooth(method = "lm") #<<
```

<!-- - ¿qué pasa si borramos `method = "lm"`? -->

---

# RLinS: retroalimentación

- $R^2$ indica el porcentaje de _variabilidad_ del desenlace (var. dependiente) explicada por la exposición (var. independiente).

- los **coeficientes** permiten __cuantificar__ el _tamaño del efecto_ de la exposición en el desenlace en base a un modelo estadístico.

- los **supuestos** permiten evaluar qué tan adecuado es el ajuste de los datos al modelo.

---

# Práctica 3

- **¿Existe una relación lineal entre talla y VEF?**

```{r,fig.height=3,fig.width=3,eval=FALSE}
espir %>% 
  ggplot(aes(x = talla,y = vef)) +
  geom____()
```

- identificar coeficientes y $R^2$

```{r,eval=FALSE}
# recordar: y ~ x
wm1 <- lm(_____ ~ _____, data = _____)
wm1 %>% g_____
wm1 %>% t_____
wm1 %>% c_____
```

- evaluar supuestos: normalidad y homoscedasticidad

```{r,fig.height=3,fig.width=3,eval=FALSE}
wm1 %>% augment() %>% 
  ggplot() + _____

wm1 %>% augment() %>% 
  ggplot(aes(.fitted,.std.resid)) + __________
```


---

# Regresión Logística Simple

.pull-left[

## características

- **una** variable independiente (simple)
- **una** variable dependiente (univariada)
- la variable _dependiente_ debe ser **categórica dicotómica**

## objetivos

- ajustar datos a una recta
- interpretar coeficientes


]


.pull-right[

```{r, echo=FALSE, fig.cap="", out.width = '100%'}
knitr::include_graphics("notrack/nmeth.3904-F23.jpg")
```

$$Y = \beta_0 + \beta_1 X_1 + \epsilon$$

]

---

# RlogS: Ecuación y coeficiente

- $p$ posee distribución _binomial_, linealizada $[-\infty,\infty]$ por la función _logit_

$$logit(p) = log(odds)  = log\left(\frac{p}{1-p}\right) \\ logit(p) = \beta_0 + \beta_1 X_1 + . . . + \beta_n X_n + \epsilon \enspace ; \enspace y \sim Binomial$$

--


```{r,fig.height=4,fig.width=8,eval=TRUE,echo=FALSE,fig.align='center'}
data_frame(p=seq(0,1,0.01)) %>% 
  mutate(odd=p/(1-p),
         logit=log(odd)) %>% 
  pivot_longer(cols = -p,names_to = "key",values_to = "value") %>% 
  mutate(key=fct_rev(key)) %>% 
  ggplot(aes(p,value)) +
  geom_point() +
  facet_wrap(~key,scales = "free") +
  coord_flip()
```

---

# RlogS: Ecuación y coeficiente


- $p$ posee distribución _binomial_, linealizada $[-\infty,\infty]$ por la función _logit_

$$logit(p) = log(odds)  = log\left(\frac{p}{1-p}\right) \\ logit(p) = \beta_0 + \beta_1 X_1 + . . . + \beta_n X_n + \epsilon \enspace ; \enspace y \sim Binomial$$

- El valor exponenciado de los coeficientes se pueden interpretar como **Odds Ratio (OR)**

$$Y = \beta_0 + \beta_1 X_1 \\ \begin{cases} Y_{x=1} = log(odds_{x=1}) = \beta_0 + \beta_1(1) \\ Y_{x=0} = log(odds_{x=0}) = \beta_0 + \beta_1(0) \end{cases} \\ Y_{x=1} - Y_{x=0} = \beta_1 \\ log(odds_{x=1}) - log(odds_{x=0}) = \beta_1 \\ log \left(\frac{odds_{x=1}}{odds_{x=0}}\right)  = \beta_1 \\ OR = exp(\beta_1)$$

---

# RlogS: Modelos Lineales Generalizados

GLM ajusta modelos lineales de $g(y)$ con covaribles $x$

$$g(Y) = \beta_0 + \sum_{n=1}^{n} \beta_n X_n \enspace , \enspace y \sim F$$

Donde: 

- $F$ es la familia de distribución
- $g(\space)$ es la función de enlace

--

```{r, echo=FALSE, fig.cap="", out.width = '100%'}
knitr::include_graphics("notrack/logan-05.PNG")
```


---

# GLM: máxima verosimilitud

.pull-left[

```{r, echo=FALSE, fig.cap="", out.width = '100%',fig.align='center'}
knitr::include_graphics("notrack/nmeth.3904-F3_A.jpg")
```

]

.pull-right[

- La estimación de parámetros se da por un proceso de optimización llamado **máxima verosimilitud** (o _likelihood_).

- un estimado por _máxima verosimilitud_ es aquel que maximiza la verosimilitud de obtener las actuales observaciones dado el modelo elegido.

- estimación numérica mediante _proceso iterativo_ hasta la convergencia.


]

---

# RlogS: Ejemplo

**En personas VIH+, ¿el polimorfismo CCR5 está *asociado* con el desarrollo del SIDA?**

- base: `aidsdb.dta`
- exposición: `ccr5`
- desenlace: `aids`

```{r,echo=FALSE}
knitr::kable(head(aidsdb), format = 'html')
```

---

# RlogS: Interpretar variables categóricas 

```{r}
wm1 <- glm(aids ~ ccr5, data = aidsdb, 
           family = binomial(link = "logit")) #<<
```

```{r,echo=FALSE}
#wm1 %>% glance()
m1 <- epitidy::epi_tidymodel_or(wm1) %>% print()
```

**INTERPRETACIÓN**

- $\beta_0$
  - En la población, el **odds** de tener sida dado que **no poseen CCR5** es **`r round(m1$or[1],2)`**,
  - con un intervalo de confianza al 95% de `r round(m1$conf.low[1],2)` a `r round(m1$conf.high[1],2)`.

---

# RlogS: Interpretar variables categóricas 

```{r}
wm1 <- glm(aids ~ ccr5, data = aidsdb, 
           family = binomial(link = "logit")) #<<
```

```{r,echo=FALSE}
#wm1 %>% glance()
m1 <- epitidy::epi_tidymodel_or(wm1) %>% print()
```

**INTERPRETACIÓN**

- $\beta_{CCR2:Yes}$
  - En la población, el **odds** de tener sida dado que **sí poseen CCR5** es
  - **`r round(m1$or[2],2)`** veces, 
  - el **odds** de tener depresión dado que **no poseen CCR5**,
  - con un intervalo de confianza al 95% de `r round(m1$conf.low[2],2)` a `r round(m1$conf.high[2],2)`.
  - Este resultado no es estadísticamente significativo con un valor **p = `r round(m1$p.value[2],3)`**

```{r,fig.height=3,fig.width=3, echo=FALSE, eval=FALSE}
#fitted values, residuals, etc
wm1 %>% augment()
```

---

# RlogS: Interpretar variables contínuas 

```{r}
wm1 <- glm(aids ~ loghivrna, data = aidsdb, 
           family = binomial(link = "logit"))
```

```{r,echo=FALSE}
#wm1 %>% glance()
m1 <- epitidy::epi_tidymodel_or(wm1) %>% print()
```

**INTERPRETACIÓN**

- $\beta_{loghivrna}$
  - En la población, por cada incremento en una unidad del **log de RNA viral**,
  - el **odds** de sida _cambia_ en **`r round(m1$or[2],2)`**, con un intervalo de confianza al 95% de `r round(m1$conf.low[2],2)` a `r round(m1$conf.high[2],2)`.
  - Este resultado es estadísticamente significativo con un valor **p < 0.001**
 
```{r,fig.height=3,fig.width=3, echo=FALSE, eval=FALSE}
#fitted values, residuals, etc
wm1 %>% augment()
```

---

# Práctica 4

- **¿Qué otros marcadores están asociados con el desarrollo de SIDA?**

- **Intenta ajustar por otras covariables como edad, raza o sexo**

```{r,eval=FALSE}
# recordar: y ~ x
wm1 <- glm(_____ ~ _____, data=_____, family= _____(link="_____"))

wm1 %>% g_____
wm1 %>% t_____
wm1 %>% c_____
wm1 %>% a_____

wm1 %>% epitidy::epi_tidymodel_or()
```

---

# Relación con Prueba de Hipótesis

```{r, echo=FALSE, fig.cap="", out.width = '100%',fig.align='center'}
knitr::include_graphics("notrack/linear_tests_cheat_sheet.png")
```

.footnote[
[1] [Jonas Kristoffer Lindeløv](https://lindeloev.github.io/tests-as-linear/)
]

---
class: inverse, center, middle

# Comparaciones múltiples


---

# Problema: _multiple hypothesis testing_

- Cuando monitoreamos un gran número de resultados experimentales, esperamos observar descubrimientos que ocurren al azar.

- Los métodos desarrollados para ajustar o reinterpretar los valores p son llamados _métodos de corrección por múltiples comparaciones_.

- p.e.: experimentos ómicos:

```{r, echo=FALSE, fig.cap="", out.width = '75%',fig.align='center'}
knitr::include_graphics("notrack/nmeth.2900-F1.jpg")
```

---

background-image: url(notrack/statquest-3.png)
background-size: contain
.footnote[
Stamer J. [StatQuest](https://www.youtube.com/user/joshstarmer/playlists)
]

---

background-image: url(notrack/statquest-6.png)
background-size: contain
.footnote[
Stamer J. [StatQuest](https://www.youtube.com/user/joshstarmer/playlists)
]

---

background-image: url(notrack/statquest-7.png)
background-size: contain
.footnote[
Stamer J. [StatQuest](https://www.youtube.com/user/joshstarmer/playlists)
]

---

background-image: url(notrack/statquest-10.png)
background-size: contain
.footnote[
Stamer J. [StatQuest](https://www.youtube.com/user/joshstarmer/playlists)
]

---

background-image: url(notrack/statquest-12.png)
background-size: contain
.footnote[
Stamer J. [StatQuest](https://www.youtube.com/user/joshstarmer/playlists)
]

---

background-image: url(notrack/statquest-13.png)
background-size: contain
.footnote[
Stamer J. [StatQuest](https://www.youtube.com/user/joshstarmer/playlists)
]

---

# Método Benjamini Hochberg

1. ordenar ascendentemente y rankear valores p
2. el mayor valor p tiene el mismo valor de p ajustado
3. el siguiente es el menor entre: (i) el previo o (ii) la solución:

.pull-left[

$$p_{adjusted} = p_{value} * \left(\frac{N_{p_{value}}}{rank_{p_{value}}}\right)$$

```{r,eval=FALSE}
tibble( #<<
  p=seq(from = 0.01,
        to = 0.91,
        by = 0.1),
  r=1:10)
```

]

--

.pull-right[

```{r,echo=FALSE}
tibble(p=seq(0.01,0.91,0.1),r=1:10)
```

]

---

# Método Benjamini Hochberg

1. ordenar ascendentemente y rankear valores p
2. el mayor valor p tiene el mismo valor de p ajustado
3. el siguiente es el menor entre: (i) el previo o (ii) la solución:

.pull-left[

$$p_{adjusted} = p_{value} * \left(\frac{N_{p_{value}}}{rank_{p_{value}}}\right)$$

```{r,eval=FALSE}
tibble(
  p=seq(from = 0.01,
        to = 0.91,
        by = 0.1),
  r=1:10) %>% 
  mutate(  #<<
    p.adjust=  #<<
      p.adjust(  #<<
        p = p,  #<<
        method = "BH")) #<<
```

]

.pull-right[

```{r,echo=FALSE}
tibble(p=seq(0.01,0.91,0.1),r=1:10)
```

]


---

# Método Benjamini Hochberg

1. ordenar ascendentemente y rankear valores p
2. el mayor valor p tiene el mismo valor de p ajustado
3. el siguiente es el menor entre: (i) el previo o (ii) la solución:

.pull-left[

$$p_{adjusted} = p_{value} * \left(\frac{N_{p_{value}}}{rank_{p_{value}}}\right)$$

```{r,eval=FALSE}
tibble(
  p=seq(from = 0.01,
        to = 0.91,
        by = 0.1),
  r=1:10) %>% 
  mutate(  #<<
    p.adjust=  #<<
      p.adjust(  #<<
        p = p,  #<<
        method = "BH")) #<<
```

]

.pull-right[

```{r,echo=FALSE}
tibble(p=seq(0.01,0.91,0.1),r=1:10) %>% 
  mutate(p.adjust=p.adjust(p = p,method = "BH")) #<<
```

]

---

background-image: url(notrack/statquest-14.png)
background-size: contain
.footnote[
Stamer J. [StatQuest](https://www.youtube.com/user/joshstarmer/playlists)
]

---

background-image: url(notrack/statquest-15.png)
background-size: contain
.footnote[
Stamer J. [StatQuest](https://www.youtube.com/user/joshstarmer/playlists)
]

---

background-image: url(notrack/statquest-16.png)
background-size: contain
.footnote[
Stamer J. [StatQuest](https://www.youtube.com/user/joshstarmer/playlists)
]

---

# Ejemplo: método BH

```{r,fig.height=5,fig.width=7}
multpv %>% 
  ggplot(aes(p.value)) + 
  geom_histogram(binwidth=.05) +
  facet_wrap(~ type, scale="free_y", nrow=2) + #<<
  labs(x="P-values")
```

---

# Ejemplo: método BH

```{r,fig.height=4,fig.width=8}
multpv %>% 
  group_by(type) %>%  #<<
  mutate(p.adjust=p.adjust(p = p.value,method = "BH"), #<<
         p.adjust_pass=if_else(p.adjust<0.05,"TRUE","FALSE")) %>% 
  ungroup() %>% 
  ggplot(aes(p.value,fill=p.adjust_pass)) + geom_histogram(binwidth=.05) +
  facet_wrap(~ type, scale="free_y", nrow=2) +labs(x="P-values",fill="BH")
```

---

# Aplicación: Microarray

- Diseño:
  * Medición de la expresión en inanición (Brauer, 2008)
  * Cultivos de _S. cerevisiae_ expuestos a seis concentraciones de nutrientes, suministrado en seis tasas o _rates_ por un quimiostato.
  * Un _gene expression microarray_ mide qué tanto un gen se expresa en determinada condición.

--

```{r}
cleaned_data <- read_rds("data-raw/microarraydata.rds")
```

.pull-left[

```{r}
cleaned_data %>% 
  count(nutrient) #<<
```

]

.pull-right[

```{r}
cleaned_data %>% 
  count(rate) #<<
```

]

```{r,echo=FALSE}
plot_expression_data <- function(expression_data) {
  ggplot(expression_data, aes(rate, expression, color = nutrient)) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    facet_wrap(~name + systematic_name, scales = "free_y")
}
```

---

# Microarray: visualización

```{r,fig.height=6,fig.width=9}
cleaned_data %>%
  filter(BP == "leucine biosynthesis") %>% #<<
  plot_expression_data() #función detallada en: practica-05.R
```

---

# Microarray: visualización

```{r,fig.height=6,fig.width=9}
cleaned_data %>%
  filter(BP == "cell wall organization and biogenesis") %>%
  plot_expression_data()
```

---

# Microarray: Regresión Lineal Simple

```{r,fig.height=4,fig.width=6}
cleaned_data %>%
  #elegimos 01 gen y 01 nutriente
  filter(name == "LEU1", nutrient == "Leucine") %>% 
  #graficamos la relación Y: expresión ~ X: rate
  ggplot(aes(rate, expression)) +
  #empleamos la geometría punto
  geom_point()
```

---

# Microarray: Regresión Lineal Simple

```{r,fig.height=4,fig.width=6}
cleaned_data %>%
  #elegimos 01 gen y 01 nutriente
  filter(name == "LEU1", nutrient == "Leucine") %>% 
  #graficamos la relación Y: expresión ~ X: rate
  ggplot(aes(rate, expression)) +
  #empleamos la geometría punto
  geom_point() + geom_smooth(method = "lm")
```


---

# Microarray: Regresión Lineal Simple

```{r}
cleaned_data %>%
  #elegimos 01 gen y 01 nutriente
  filter(name == "LEU1", nutrient == "Leucine") %>% 
  #ajustamos una regresión lineal
  #dado que data no es el primer argumento
  #necesitamos especificarlo en data con "."
  lm(expression ~ rate, data = .) %>% 
  #visualizamos tabla con estimados
  tidy()
```

---

# Microarray: a todas los pares gen-nutriente

```{r}
linear_models <- cleaned_data %>%
  filter(nutrient=="Ammonia") %>% #filtramos por nutriente #<<
  group_by(name, systematic_name, nutrient) %>% #agrupamos por gen
  nest() %>% #anidamos los datos en una columna lista de df #<<
  # ajustamos un modelo lineal a cada fila -ver paquete purrr::map-
  mutate(model = map(data, ~ lm(expression ~ rate, data = .x)), #<<
         tidym = map(model,tidy)) #<<
```

--

```{r}
linear_models
```


---

# Microarray: conservamos las pendientes

```{r,fig.height=4,fig.width=10}
slope_terms <- linear_models %>% 
  unnest(cols = c(tidym)) %>% #<<
  ungroup() %>% 
  filter(term=="rate" & !is.na(p.value)) 

slope_terms %>% 
  ggplot(aes(p.value)) +
  geom_histogram(binwidth = .01) +
  facet_wrap(~nutrient)
```

---

# Microarray: corregimos valor p

```{r,fig.height=4,fig.width=10}
slope_terms_adj <- slope_terms %>% 
  mutate(q.value = qvalue(p.value)$qvalues, #<<
         q.value_pass=if_else(q.value < .01,"TRUE","FALSE")) 

slope_terms_adj %>% 
  ggplot(aes(p.value,
             fill=q.value_pass)) + #<<
  geom_histogram(binwidth = .01) +
  facet_wrap(~nutrient)
```

---

# Microarray: Lista de genes de interés

```{r}
slope_terms_adj %>% 
  filter(q.value_pass=="TRUE") %>% 
  arrange(q.value) %>%  #<<
  select(systematic_name,nutrient,term,estimate,p.value,q.value)
```

---

# Práctica 5

- reproducir el proceso para los demás nutrientes

- explorar paquete `limma`

- explorar la documentación del paquete `qvalue` [link](https://github.com/StoreyLab/qvalue)

```{r,eval=FALSE}
library(qvalue)

#datos microarray
qmicro <- qvalue(p = slope_terms$p.value)
summary(qmicro)
hist(qmicro)
plot(qmicro)
```


---

# Referencias

- Logan M. (2011). Biostatistical design and analysis using R: a practical guide. John Wiley & Sons. [Library Genesis](http://gen.lib.rus.ec/book/index.php?md5=9CF463AE32C405E29FCEDD99440225BE)
- Akehurst H. (2016) Bioestadística con R. Curso dictado en UPCH.
- Statistics for Biologist: Points of Significance. [Nature Collection](https://www.nature.com/collections/qghhqm/pointsofsignificance)
- Ho J, et al. (2019). Moving beyond P values: data analysis with estimation graphics. [Nature Methods](https://www.nature.com/articles/s41592-019-0470-3.epdf?author_access_token=Euy6APITxsYA3huBKOFBvNRgN0jAjWel9jnR3ZoTv0Pr6zJiJ3AA5aH4989gOJS_dajtNr1Wt17D0fh-t4GFcvqwMYN03qb8C33na_UrCUcGrt-Z0J9aPL6TPSbOxIC-pbHWKUDo2XsUOr3hQmlRew%3D%3D)
- Stamer J. StatQuest. False Discovery Rates, FDR, clearly explained. [YouTube channel](https://www.youtube.com/user/joshstarmer/playlists)
- Ding J, et al. (2018). Association of gene polymorphism of SDF1(CXCR12) with susceptibility to HIV-1 infection and AIDS disease progression: A meta-analysis. [PLoS One](https://www.ncbi.nlm.nih.gov/pubmed/29420545)
- Estrada-Aguirre JA, et al. (2013). Protective effect of CCR5 Delta-32 allele against HIV-1 in Mexican women. [Curr HIV Res](https://www.ncbi.nlm.nih.gov/pubmed/24382026)
- Lindeløv J. (2019) Common statistical tests are linear models (or: how to teach stats) [Github personal webpage](https://lindeloev.github.io/tests-as-linear/)
- Robinson D. (2014). Modeling gene expression with broom: a case study in tidy analysis. [Variance explained blog](http://varianceexplained.org/r/tidy-genomics-broom/)

---

# Material complementario

- Rstudio Education: [Learn](https://education.rstudio.com/learn/)
  * beginner, intermediate, expert tracks
  * transform table [here](https://rstudio-education.github.io/tidyverse-cookbook/transform-tables.html)
  * tidy explain *_join verbs [here](https://www.garrickadenbuie.com/project/tidyexplain/)

- Rstudio [Cheat sheets](https://rstudio.com/resources/cheatsheets/)
  * data import
  * data transformation
  * data visualization

---

class: center, middle

# ¡Gracias!

## Andree Valle Campos<br/><br/> `r icons::fontawesome$brands$twitter` `r icons::fontawesome$brands$github` __`@avallecam`__ <br/><br/> `r icons::fontawesome$solid$inbox` __avallecam@gmail.com__

.footnote[
Slides created via the R package [**xaringan**](https://github.com/yihui/xaringan).
]

---

class: center, middle

background-image: url(https://upload.wikimedia.org/wikipedia/commons/3/39/Naruto_Shiki_Fujin.svg)
background-size: contain