Modelos de regresión y técnicas de análisis de datos multivariables con el programa R

<br>
<p style="font-size:40px;text-align:center;font-weight:bold">Curso modelos de regresión y técnicas de análisis de datos multivariables con el programa R</p>

<p style="font-size:30px;text-align:center;margin-top:0px">AZTI</p>
<p style="font-size:25px;text-align:center;margin-top:-10px">diciembre 2020</p>
<p style="font-size:25px;text-align:center;margin-top:-10px">Isaac Subirana</p>
<br><br>

<p style="font-size:40px;text-align:center;padding-top:0%;color:white">Bloque 1: Modelos de regresión</p>

---

## Guión del curso

**Primer bloque:** Modelos de regresión

- Modelos de regresión clásicos: modelos lineales y regresión logística.
- Modelos de regresión con muchas variables independientes: elastic net, penalized regression models.

**Segundo bloque:** Medidas repetidas o datos longitudinales

- Respuesta normal: modelos lineales mixtos LMM
- Respuesta binaria: modelos lineales mixtos generalizados GLMM,
- Otros métodos: generalized estimation equation GEE.

---

**Tercer bloque:** Métodos multivariantes para relacionar variables

- Para un conjunto de variables: análisis de componentes principales PCA
- Para dos conjuntos de variables: análisis de correlaciones canónicas CCA
- Para más de dos conjuntos de variables: análisis de correlaciones canónicas generalizado (GCCA)
- Métodos con muchas variables: penalized methods.
- Métodos para variables continuas con distribución no normal o categóricas -> MOFA.

**Cuarto bloque:** Modelos para respuesta ordinal:

- Correlaciones para variables ordinales
- Modelos de regresión ordinal

**Ejemplos y análisis de datos reales**

---
class: inverse, center, middle

# Modelos de regresión

---

## Introducción

**Variables**

- Una variable respuesta

- Una o varias variables explicativas

**Objetivos**

1. Predicción: **modelos predictivos**

2. Estimar el efecto de una variable explicativa: **modelos de causalidad o asociación**

---

## Variable respuesta

**Cuantitativas**

- Continuas: **Normal**, log-Normal, ...
- Conteos: Poisson, Binomial negativa, ...
- Supervivencia: Weibull, ...

**Cualitativas**

- **Binaria** (sí/no)
- **Multinomial** (varias categorías)
- **Ordinal** (categorías ordenadas)

**Otros**
  
- Censuradas
- Truncadas
- Infladas en el cero

---

## Variables explicativas

- Numéricas: efecto lineal, polinómico, splines, ...

- Factores: definir categoría de referencia (variables dummies).

- Interacciones: entre dos o más variables explicativas.

---

## Regresión lineal

- La variable respuesta tiene que ser normal

- Transformaciones (logaritmica, etc.)

- Modelo

`$$y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots \beta_p x_{ip} + e_{i}$$`

- donde

- `$y_i$`: valor observado de la variable respuesta para el individuo `$i$`.

- `$\beta_k$`: coeficiente de la variable `$x_k$`. Cuando `$x_{ik}$` se incrementa una unidad el valor esperado de la variable respuesta aumenta `$\beta_k$` unidades.

- `$x_{ik}$`: valor observado de la variable `$k$` del individuo `$i$`. Puede ser el valor, una interacción, dummy variable,... (**matriz de diseño**)

- `$e_i \sim N(0, \sigma^2)$` independientes (premisas del modelo).

---

## Regresión logística

- La variable respuesta es binaria (0/1), (sí/no),...

- Modelo.

`$$\text{log}\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots \beta_p x_{ip}$$`

- donde

- `$p_i$`: es la probabilidad que ocurra el evento para el individuo `$i$`.

- `$e^{\beta_k}$`: es el Odds Ratio de la variable `$x_k$`. Cuando `$x_{ik}$` se incrementa una unidad, el odds se incrementa `$e^{\beta_k}$`.

- `$x_{ik}$`: lo mismo que para el modelo de regresión lineal.

- Nota que aquí no hay errores.

---

## Matriz de diseño

- Matriz con tantas filas como individuos

- La primera columna continene unos (para la constante)

- A partir de la segunda columna contiene cadad término `$x_k$`, `$k=1,\ldots,p$`.

- Las columnas no siempre corresponden a cada variable explicativa:

- Si es numérica y se supone un efecto lineal `$\rightarrow$` `$x_k$`.
  
  - Si es numérica y se supone un efecto polinómico `$\rightarrow$` tantos términos `$x_k$` como polinomios.
  
  - Si es numérica y se supone un efecto splines `$\rightarrow$` tantos términos `$x_k$` como términos de splines.
  
  - Si es categórica con `$c$` categorías `$\rightarrow$` `$c-1$` dummy variables.
  
---

**Categóricas binarias**

**Categóricas binarias > 2 categorías**

**Interacciones**

---

**Numéricas (polinomio)**

<table class="table" style="margin-left: auto; margin-right: auto;">
 <thead>
<tr>
<th style="border-bottom:hidden; padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="1"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Variable</div></th>
<th style="border-bottom:hidden; padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Polinomios ortogonales</div></th>
<th style="border-bottom:hidden; padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Polinomios</div></th>
</tr>
  <tr>
   <th style="text-align:left;"> tiempo </th>
   <th style="text-align:center;"> poly(tiempo, 2, raw = TRUE)1 </th>
   <th style="text-align:center;"> poly(tiempo, 2, raw = TRUE)2 </th>
   <th style="text-align:center;"> poly(tiempo, 2, raw = FALSE)1 </th>
   <th style="text-align:center;"> poly(tiempo, 2, raw = FALSE)2 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 8.0 </td>
   <td style="text-align:center;"> 8.0 </td>
   <td style="text-align:center;"> 64.00 </td>
   <td style="text-align:center;"> 0.38 </td>
   <td style="text-align:center;"> -0.07 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 7.5 </td>
   <td style="text-align:center;"> 7.5 </td>
   <td style="text-align:center;"> 56.25 </td>
   <td style="text-align:center;"> 0.32 </td>
   <td style="text-align:center;"> -0.18 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3.9 </td>
   <td style="text-align:center;"> 3.9 </td>
   <td style="text-align:center;"> 15.21 </td>
   <td style="text-align:center;"> -0.08 </td>
   <td style="text-align:center;"> -0.29 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3.4 </td>
   <td style="text-align:center;"> 3.4 </td>
   <td style="text-align:center;"> 11.56 </td>
   <td style="text-align:center;"> -0.14 </td>
   <td style="text-align:center;"> -0.22 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3.6 </td>
   <td style="text-align:center;"> 3.6 </td>
   <td style="text-align:center;"> 12.96 </td>
   <td style="text-align:center;"> -0.12 </td>
   <td style="text-align:center;"> -0.25 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 2.0 </td>
   <td style="text-align:center;"> 2.0 </td>
   <td style="text-align:center;"> 4.00 </td>
   <td style="text-align:center;"> -0.29 </td>
   <td style="text-align:center;"> 0.13 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 5.3 </td>
   <td style="text-align:center;"> 5.3 </td>
   <td style="text-align:center;"> 28.09 </td>
   <td style="text-align:center;"> 0.08 </td>
   <td style="text-align:center;"> -0.39 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1.0 </td>
   <td style="text-align:center;"> 1.0 </td>
   <td style="text-align:center;"> 1.00 </td>
   <td style="text-align:center;"> -0.41 </td>
   <td style="text-align:center;"> 0.48 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 9.9 </td>
   <td style="text-align:center;"> 9.9 </td>
   <td style="text-align:center;"> 98.01 </td>
   <td style="text-align:center;"> 0.59 </td>
   <td style="text-align:center;"> 0.57 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1.7 </td>
   <td style="text-align:center;"> 1.7 </td>
   <td style="text-align:center;"> 2.89 </td>
   <td style="text-align:center;"> -0.33 </td>
   <td style="text-align:center;"> 0.22 </td>
  </tr>
</tbody>
</table>

---

**Numéricas (splines)**

<table class="table" style="margin-left: auto; margin-right: auto;">
 <thead>
<tr>
<th style="border-bottom:hidden; padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="1"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Variable</div></th>
<th style="border-bottom:hidden; padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="3"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Splines</div></th>
</tr>
  <tr>
   <th style="text-align:left;"> tiempo </th>
   <th style="text-align:center;"> ns(tiempo, 3)1 </th>
   <th style="text-align:center;"> ns(tiempo, 3)2 </th>
   <th style="text-align:center;"> ns(tiempo, 3)3 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 8.0 </td>
   <td style="text-align:center;"> 0.32 </td>
   <td style="text-align:center;"> 0.34 </td>
   <td style="text-align:center;"> 0.31 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 7.5 </td>
   <td style="text-align:center;"> 0.41 </td>
   <td style="text-align:center;"> 0.34 </td>
   <td style="text-align:center;"> 0.19 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3.9 </td>
   <td style="text-align:center;"> 0.01 </td>
   <td style="text-align:center;"> 0.60 </td>
   <td style="text-align:center;"> -0.35 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3.4 </td>
   <td style="text-align:center;"> -0.10 </td>
   <td style="text-align:center;"> 0.60 </td>
   <td style="text-align:center;"> -0.35 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3.6 </td>
   <td style="text-align:center;"> -0.06 </td>
   <td style="text-align:center;"> 0.60 </td>
   <td style="text-align:center;"> -0.35 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 2.0 </td>
   <td style="text-align:center;"> -0.12 </td>
   <td style="text-align:center;"> 0.33 </td>
   <td style="text-align:center;"> -0.19 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 5.3 </td>
   <td style="text-align:center;"> 0.39 </td>
   <td style="text-align:center;"> 0.46 </td>
   <td style="text-align:center;"> -0.22 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1.0 </td>
   <td style="text-align:center;"> 0.00 </td>
   <td style="text-align:center;"> 0.00 </td>
   <td style="text-align:center;"> 0.00 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 9.9 </td>
   <td style="text-align:center;"> -0.16 </td>
   <td style="text-align:center;"> 0.39 </td>
   <td style="text-align:center;"> 0.77 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1.7 </td>
   <td style="text-align:center;"> -0.09 </td>
   <td style="text-align:center;"> 0.23 </td>
   <td style="text-align:center;"> -0.14 </td>
  </tr>
</tbody>
</table>

---

### Ejemplos con R

```r
library(compareGroups)
data(predimed)
summary(predimed)
```

```
            group          sex            age            smoke     
 Control       :2042   Male  :2679   Min.   :49.00   Never  :3892  
 MedDiet + Nuts:2100   Female:3645   1st Qu.:62.00   Current: 858  
 MedDiet + VOO :2182                 Median :67.00   Former :1574  
                                     Mean   :67.01                 
                                     3rd Qu.:72.00                 
                                     Max.   :87.00                 
      bmi            waist            wth          htn        diab     
 Min.   :19.64   Min.   : 50.0   Min.   :0.3012   No :1089   No :3322  
 1st Qu.:27.23   1st Qu.: 93.0   1st Qu.:0.5839   Yes:5235   Yes:3002  
 Median :29.76   Median :100.0   Median :0.6258                        
 Mean   :29.97   Mean   :100.4   Mean   :0.6283                        
 3rd Qu.:32.46   3rd Qu.:107.0   3rd Qu.:0.6687                        
 Max.   :51.94   Max.   :177.0   Max.   :1.0000                        
 hyperchol  famhist     hormo           p14            toevent       
 No :1746   No :4895   No  :5564   Min.   : 0.000   Min.   :0.01643  
 Yes:4578   Yes:1429   Yes :  97   1st Qu.: 7.000   1st Qu.:2.85832  
                       NA's: 663   Median : 9.000   Median :4.78850  
                                   Mean   : 8.678   Mean   :4.35517  
                                   3rd Qu.:10.000   3rd Qu.:5.79056  
                                   Max.   :14.000   Max.   :6.99795  
 event     
 No :6072  
 Yes: 252  
           
           
           
           
```

---

```r
X <- model.matrix( ~ group + sex + age, data=predimed)
head(X)
```

```
  (Intercept) groupMedDiet + Nuts groupMedDiet + VOO sexFemale age
1           1                   0                  0         0  58
2           1                   0                  0         0  77
4           1                   0                  1         1  72
5           1                   1                  0         0  71
6           1                   0                  1         1  79
8           1                   0                  0         0  63
```

```r
X <- model.matrix( ~ poly(age, 2) + group, data=predimed)
head(X)
```

```
  (Intercept) poly(age, 2)1 poly(age, 2)2 groupMedDiet + Nuts
1           1  -0.018353081   0.014965639                   0
2           1   0.020342002   0.016202915                   0
4           1   0.010159085  -0.005104269                   0
5           1   0.008122502  -0.007567284                   1
6           1   0.024415168   0.028922106                   0
8           1  -0.008170164  -0.005690347                   0
  groupMedDiet + VOO
1                  0
2                  0
4                  1
5                  0
6                  1
8                  0
```

---
class: inverse, center, middle

# Regresión lineal con R

---

Ejemplo de `predimed` del package `compareGroups`

- Variable respuesta: `bmi`

- Variables predictoras: `group`, `age`, `sex`, `smoke`, `diab` y `p14`

---

## Descriptiva univariante de las variables

```r
tabla <- descrTable(~ bmi + group + age + sex + smoke + diab + p14, 
                    data = predimed)
tabla
```

```

--------Summary descriptives table ---------

________________________________________ 
                          [ALL]      N   
                          N=6324         
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
Body mass index        30.0 (3.82)  6324 
Intervention group:                 6324 
    Control            2042 (32.3%)      
    MedDiet + Nuts     2100 (33.2%)      
    MedDiet + VOO      2182 (34.5%)      
Age                    67.0 (6.17)  6324 
Sex:                                6324 
    Male               2679 (42.4%)      
    Female             3645 (57.6%)      
Smoking:                            6324 
    Never              3892 (61.5%)      
    Current            858 (13.6%)       
    Former             1574 (24.9%)      
Type-2 diabetes:                    6324 
    No                 3322 (52.5%)      
    Yes                3002 (47.5%)      
MeDiet Adherence score 8.68 (1.94)  6324 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
```

Comprobar también normalidad de las variables

```r
plot(tabla)
```

---

## Asociación de cada variable los terciles de BMI

```r
library(Hmisc)
predimed$terciles <- cut2(predimed$bmi, g=3) 
descrTable(terciles ~ group + age + sex + smoke + diab + p14, 
           data = predimed, 
           show.p.trend=TRUE,
           show.p.overall=FALSE)
```

```

--------Summary descriptives table by 'Body mass index'---------

_____________________________________________________________________ 
                       [19.6,28.1)  [28.1,31.5)  [31.5,51.9]  p.trend 
                          N=2118       N=2105       N=2101            
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
Intervention group:                                            0.047  
    Control            625 (29.5%)  680 (32.3%)  737 (35.1%)          
    MedDiet + Nuts     765 (36.1%)  705 (33.5%)  630 (30.0%)          
    MedDiet + VOO      728 (34.4%)  720 (34.2%)  734 (34.9%)          
Age                    67.1 (6.28)  67.1 (6.15)  66.8 (6.08)   0.082  
Sex:                                                           0.000  
    Male               1016 (48.0%) 984 (46.7%)  679 (32.3%)          
    Female             1102 (52.0%) 1121 (53.3%) 1422 (67.7%)         
Smoking:                                                      <0.001  
    Never              1212 (57.2%) 1227 (58.3%) 1453 (69.2%)         
    Current            337 (15.9%)  302 (14.3%)  219 (10.4%)          
    Former             569 (26.9%)  576 (27.4%)  429 (20.4%)          
Type-2 diabetes:                                               0.442  
    No                 1087 (51.3%) 1132 (53.8%) 1103 (52.5%)         
    Yes                1031 (48.7%) 973 (46.2%)  998 (47.5%)          
MeDiet Adherence score 8.91 (1.91)  8.69 (1.94)  8.43 (1.96)  <0.001  
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
```

---

## Asociación de cada variable con el BMI

- `bmi` vs `group`: ANOVA

```r
oneway.test(bmi ~ group, data=predimed)
```

```

One-way analysis of means (not assuming equal variances)

data:  bmi and group
F = 12.203, num df = 2.0, denom df = 4194.8, p-value = 5.196e-06
```

- `bmi` vs `age`: correlación de Pearson

```r
with(predimed, cor.test(bmi, age, method="pearson"))
```

```

Pearson's product-moment correlation

data:  bmi and age
t = -3.1839, df = 6322, p-value = 0.00146
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.06459451 -0.01537902
sample estimates:
        cor 
-0.04001103 
```

---

- `bmi` vs `sex`: t-student

```r
t.test(bmi ~ sex, data=predimed)
```

```

Welch Two Sample t-test

data:  bmi by sex
t = -12.504, df = 6247.2, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.3484451 -0.9829292
sample estimates:
  mean in group Male mean in group Female 
            29.29430             30.45998 
```

- `bmi` vs `smoking`: ANOVA

```r
oneway.test(bmi ~ smoke, data=predimed)
```

```

One-way analysis of means (not assuming equal variances)

data:  bmi and smoke
F = 48.084, num df = 2.0, denom df = 2189.5, p-value < 2.2e-16
```

---

- `bmi` vs `diab`: t-test

```r
t.test(bmi ~ diab, data=predimed)
```

```

Welch Two Sample t-test

data:  bmi by diab
t = 1.2256, df = 6078.6, p-value = 0.2204
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.07099949  0.30786502
sample estimates:
 mean in group No mean in group Yes 
         30.02239          29.90396 
```

- `bmi` vs `p14`: correlación e Pearson

```r
with(predimed, cor.test(bmi, p14, method="pearson"))
```

```

Pearson's product-moment correlation

data:  bmi and p14
t = -6.9872, df = 6322, p-value = 3.091e-12
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.11194566 -0.06302883
sample estimates:
        cor 
-0.08754002 
```

---

- Todas las variables parecen significativas de forma bivariada.

- Pero, como entre ellas hay correlación, en el modelo puede que algunas no lo sean.

- Para ajustar un modelo de regresión lineal se usa la función **`lm`**

- Las variables se especifican con la fórmula (`y ~ x1 + x2 + ...`)

- Interacciones: `:`

- Términos principales e interacciones: `*` = `var1 + var2 + var1:var2`

- Transformaciones: `I(.)`

---

## Ajuste del modelo

```r
modelo <- lm(bmi ~ group + age + sex + smoke + diab + p14, data=predimed)
```

```r
modelo
```

```

Call:
lm(formula = bmi ~ group + age + sex + smoke + diab + p14, data = predimed)

Coefficients:
        (Intercept)  groupMedDiet + Nuts   groupMedDiet + VOO  
           34.14550             -0.49463             -0.28499  
                age            sexFemale         smokeCurrent  
           -0.04394              0.95750             -0.73847  
        smokeFormer              diabYes                  p14  
           -0.19282             -0.04584             -0.15607  
```

---

## Extracción de los resultados

```r
summary(modelo)
```

```

Call:
lm(formula = bmi ~ group + age + sex + smoke + diab + p14, data = predimed)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.2504  -2.6500  -0.1203   2.4652  21.2437

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         34.14550    0.58540  58.329  < 2e-16 ***
groupMedDiet + Nuts -0.49463    0.11701  -4.227 2.40e-05 ***
groupMedDiet + VOO  -0.28499    0.11567  -2.464   0.0138 *  
age                 -0.04394    0.00781  -5.626 1.93e-08 ***
sexFemale            0.95750    0.12263   7.808 6.75e-15 ***
smokeCurrent        -0.73847    0.16134  -4.577 4.80e-06 ***
smokeFormer         -0.19282    0.13711  -1.406   0.1597    
diabYes             -0.04584    0.09555  -0.480   0.6314    
p14                 -0.15607    0.02441  -6.394 1.73e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.746 on 6315 degrees of freedom
Multiple R-squared:  0.03877,	Adjusted R-squared:  0.03755 
F-statistic: 31.84 on 8 and 6315 DF,  p-value: < 2.2e-16
```

**Fíjate en la `$R^2$` y la `$R^2$` ajustada.**

---

### Coeficientes

```r
confint(modelo)
```

```
                          2.5 %      97.5 %
(Intercept)         32.99792935 35.29308045
groupMedDiet + Nuts -0.72401711 -0.26525234
groupMedDiet + VOO  -0.51174501 -0.05824197
age                 -0.05924769 -0.02862609
sexFemale            0.71710324  1.19789484
smokeCurrent        -1.05474568 -0.42219744
smokeFormer         -0.46160733  0.07596329
diabYes             -0.23315389  0.14147235
p14                 -0.20391137 -0.10822088
```

```r
coef(modelo)
```

```
        (Intercept) groupMedDiet + Nuts  groupMedDiet + VOO                 age 
        34.14550490         -0.49463472         -0.28499349         -0.04393689 
          sexFemale        smokeCurrent         smokeFormer             diabYes 
         0.95749904         -0.73847156         -0.19282202         -0.04584077 
                p14 
        -0.15606613 
```

---

### Residuos

```r
resid(modelo)[1:10]
```

```
         1          2          4          5          6          8          9 
 3.6863179  2.6326091  0.6926374 -1.4100000  5.9975989 11.9257770 -4.3283740 
        12         13         14 
-3.5602408  2.2817620 -4.1273939 
```

---

### Valores predichos

```r
predict(modelo)[1:10]
```

```
       1        2        4        5        6        8        9       12 
29.84368 28.41739 30.16736 29.09000 29.94240 29.73422 30.22837 29.51024 
      13       14 
28.64824 28.61739 
```

---

Predicción del BMI para una mujer de 55 años asignada al aceite de oliva (MedDiet+VOO), fumadora, no diabética y con p14=7

```r
nuevo <- data.frame(group = "MedDiet + VOO",
                    age = 55,
                    sex="Female",
                    smoke = "Current",
                    diab="No",
                    p14 = 7)
predict(modelo, newdata=nuevo, interval="confidence")
```

```
       fit     lwr      upr
1 30.57055 30.1929 30.94819
```

Y para un grupo de mujeres de las mismas caracteísticas

```r
predict(modelo, newdata=nuevo, interval="prediction")
```

```
       fit      lwr      upr
1 30.57055 23.21771 37.92339
```

---

## Validación del modelo

---

### Normalidad de los residuos

```r
plot(modelo, which=2)
```

![](bloque1_files/figure-html/unnamed-chunk-28-1.png)

---

### Igualdad de varianzas

```r
plot(modelo, which=1)
```

![](bloque1_files/figure-html/unnamed-chunk-29-1.png)

---

### Colinealidad entre las variables independientes

VIF ("Variance Inflation Factor) para cada variable

- Entre 1 y 5: colinealidad pequeña-moderada
- Entre 5 y 10: colinealidad moderada-grande
- Superior a 10: colinealidad exceciva `$\rightarrow$` quitar la variable

```r
library(car)
vif(modelo)
```

```
          GVIF Df GVIF^(1/(2*Df))
group 1.013132  2        1.003267
age   1.048162  1        1.023798
sex   1.654889  1        1.286425
smoke 1.681441  2        1.138729
diab  1.026104  1        1.012968
p14   1.014632  1        1.007289
```

---

### Valores anómalos

```r
plot(modelo, which=3)
```

![](bloque1_files/figure-html/unnamed-chunk-31-1.png)

---

### Valores influyentes

```r
plot(modelo, which=4)
```

![](bloque1_files/figure-html/unnamed-chunk-32-1.png)

---

### Linealidad

```r
plot(modelo, which=1)
```

![](bloque1_files/figure-html/unnamed-chunk-33-1.png)

---

## Variables independientes

---

### Test de una variable categórica > 2 categorías

```r
anova(modelo, update(modelo, . ~ . - smoke))
```

```
Analysis of Variance Table

Model 1: bmi ~ group + age + sex + smoke + diab + p14
Model 2: bmi ~ group + age + sex + diab + p14
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1   6315 88608                                  
2   6317 88908 -2   -299.53 10.674 2.357e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

### Interacción

```r
modeloInt <- update(modelo, . ~ . + sex:group) 
coef(summary(modelo))
```

```
                       Estimate  Std. Error    t value     Pr(>|t|)
(Intercept)         34.14550490 0.585396266 58.3288737 0.000000e+00
groupMedDiet + Nuts -0.49463472 0.117011549 -4.2272299 2.399221e-05
groupMedDiet + VOO  -0.28499349 0.115669503 -2.4638602 1.377143e-02
age                 -0.04393689 0.007810281 -5.6255198 1.928804e-08
sexFemale            0.95749904 0.122629663  7.8080541 6.754218e-15
smokeCurrent        -0.73847156 0.161336386 -4.5772165 4.802260e-06
smokeFormer         -0.19282202 0.137111595 -1.4063145 1.596800e-01
diabYes             -0.04584077 0.095551358 -0.4797501 6.314217e-01
p14                 -0.15606613 0.024406608 -6.3944209 1.726860e-10
```

```r
anova(modelo, modeloInt)
```

```
Analysis of Variance Table

Model 1: bmi ~ group + age + sex + smoke + diab + p14
Model 2: bmi ~ group + age + sex + smoke + diab + p14 + group:sex
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1   6315 88608                           
2   6313 88604  2    4.5092 0.1606 0.8516
```

---

### Cambio de grupo de referencia

```r
predimed$sex <- factor(predimed$sex, c("Female","Male"))
modeloInt2 <- update(modeloInt)
coef(summary(modeloInt2))
```

```
                               Estimate  Std. Error     t value     Pr(>|t|)
(Intercept)                 35.12082823 0.581024784 60.44635135 0.000000e+00
groupMedDiet + Nuts         -0.54867343 0.154718765 -3.54626298 3.935698e-04
groupMedDiet + VOO          -0.29295543 0.149753783 -1.95624726 5.048007e-02
age                         -0.04393249 0.007811339 -5.62419383 1.943627e-08
sexMale                     -1.00698596 0.186742226 -5.39238491 7.205187e-08
smokeCurrent                -0.73846291 0.161362182 -4.57643110 4.820284e-06
smokeFormer                 -0.19192102 0.137141457 -1.39943836 1.617307e-01
diabYes                     -0.04647612 0.095574014 -0.48628402 6.267827e-01
p14                         -0.15587159 0.024412517 -6.38490442 1.837044e-10
groupMedDiet + Nuts:sexMale  0.12380927 0.235826771  0.52500091 5.996010e-01
groupMedDiet + VOO:sexMale   0.02091249 0.235060669  0.08896634 9.291115e-01
```

---

### Elección de las variables bajo un modelo predictivo

- Criterio stepwise (AIC)

- No puede haber ningun missing en ninguna variable (`na.omit(datos)`)

```r
modeloFin <- step(modelo, direction = "backward")
```

```
Start:  AIC=16712.55
bmi ~ group + age + sex + smoke + diab + p14

Df Sum of Sq   RSS   AIC
- diab   1      3.23 88611 16711
<none>               88608 16713
- group  2    252.35 88861 16727
- smoke  2    299.53 88908 16730
- age    1    444.04 89052 16742
- p14    1    573.72 89182 16751
- sex    1    855.43 89464 16771

Step:  AIC=16710.79
bmi ~ group + age + sex + smoke + p14

Df Sum of Sq   RSS   AIC
<none>               88611 16711
- group  2    251.57 88863 16725
- smoke  2    296.43 88908 16728
- age    1    451.66 89063 16741
- p14    1    570.90 89182 16749
- sex    1    878.78 89490 16771
```

---

```r
coef(summary(modeloFin))
```

```
                       Estimate  Std. Error   t value     Pr(>|t|)
(Intercept)         35.09568415 0.578095892 60.709105 0.000000e+00
groupMedDiet + Nuts -0.49364546 0.116986248 -4.219688 2.480665e-05
groupMedDiet + VOO  -0.28614888 0.115637381 -2.474536 1.336700e-02
age                 -0.04420148 0.007790309 -5.673906 1.457683e-08
sexMale             -0.96418741 0.121827188 -7.914386 2.913497e-15
smokeCurrent        -0.73207838 0.160775291 -4.553426 5.376955e-06
smokeFormer         -0.19121398 0.137062266 -1.395088 1.630382e-01
p14                 -0.15549854 0.024376431 -6.379053 1.908088e-10
```

---

Podemos pensar en incluir todas las interacciones con el sexo. Y luego ver qué variables y/o interacciones son significativas

```r
modeloSex <- update(modelo, . ~ (.)*sex)
coef(summary(modeloSex))
```

```
                                Estimate Std. Error     t value     Pr(>|t|)
(Intercept)                 35.668097186 0.80320925 44.40697996 0.000000e+00
groupMedDiet + Nuts         -0.523817484 0.15443789 -3.39176802 6.986999e-04
groupMedDiet + VOO          -0.276493701 0.14934989 -1.85131511 6.417094e-02
age                         -0.045231949 0.01084694 -4.17002057 3.086536e-05
sexMale                     -2.310769693 1.15400099 -2.00239837 4.528463e-02
smokeCurrent                -1.826568169 0.28253102 -6.46501807 1.088700e-10
smokeFormer                 -0.890614922 0.23812963 -3.74004240 1.856203e-04
diabYes                      0.151222291 0.12565335  1.20348791 2.288327e-01
p14                         -0.208262432 0.03248473 -6.41108710 1.549432e-10
groupMedDiet + Nuts:sexMale  0.084389047 0.23570210  0.35803265 7.203309e-01
groupMedDiet + VOO:sexMale  -0.011881472 0.23477129 -0.05060871 9.596389e-01
age:sexMale                 -0.002155661 0.01559959 -0.13818702 8.900970e-01
sexMale:smokeCurrent         1.857431256 0.34867740  5.32707670 1.032613e-07
sexMale:smokeFormer          1.314324692 0.29517069  4.45276146 8.623400e-06
sexMale:diabYes             -0.477571362 0.19262464 -2.47928487 1.319054e-02
sexMale:p14                  0.126551430 0.04902233  2.58150561 9.859397e-03
```

---

```r
modeloSexFin <- step(modeloSex)
```

```
Start:  AIC=16670.77
bmi ~ group + age + sex + smoke + diab + p14 + group:sex + age:sex + 
    sex:smoke + sex:diab + sex:p14

Df Sum of Sq   RSS   AIC
- group:sex  2      2.85 87833 16667
- age:sex    1      0.27 87830 16669
<none>                   87830 16671
- sex:diab   1     85.59 87916 16675
- sex:p14    1     92.79 87923 16675
- sex:smoke  2    541.75 88372 16706

Step:  AIC=16666.97
bmi ~ group + age + sex + smoke + diab + p14 + age:sex + sex:smoke + 
    sex:diab + sex:p14

Df Sum of Sq   RSS   AIC
- age:sex    1      0.32 87833 16665
<none>                   87833 16667
- sex:diab   1     86.69 87920 16671
- sex:p14    1     94.44 87927 16672
- group      2    244.04 88077 16681
- sex:smoke  2    540.49 88373 16702

Step:  AIC=16664.99
bmi ~ group + age + sex + smoke + diab + p14 + sex:smoke + sex:diab + 
    sex:p14

Df Sum of Sq   RSS   AIC
<none>                   87833 16665
- sex:diab   1     87.98 87921 16669
- sex:p14    1     94.45 87928 16670
- group      2    244.10 88077 16679
- age        1    490.66 88324 16698
- sex:smoke  2    553.35 88387 16701
```

---

```r
coef(summary(modeloSexFin))
```

```
                        Estimate  Std. Error   t value     Pr(>|t|)
(Intercept)          35.73306627 0.612425834 58.346765 0.000000e+00
groupMedDiet + Nuts  -0.48641854 0.116557016 -4.173224 3.043537e-05
groupMedDiet + VOO   -0.28176462 0.115203671 -2.445795 1.448042e-02
age                  -0.04627639 0.007793827 -5.937570 3.046596e-09
sexMale              -2.43450454 0.471121796 -5.167463 2.445797e-07
smokeCurrent         -1.82913724 0.280912425 -6.511415 8.018103e-11
smokeFormer          -0.89265768 0.237582947 -3.757246 1.733421e-04
diabYes               0.15365556 0.125481160  1.224531 2.207977e-01
p14                  -0.20882094 0.032417929 -6.441526 1.270000e-10
sexMale:smokeCurrent  1.86180704 0.345321115  5.391524 7.239725e-08
sexMale:smokeFormer   1.31553169 0.294547839  4.466275 8.098049e-06
sexMale:diabYes      -0.48252989 0.191914501 -2.514296 1.195170e-02
sexMale:p14           0.12726712 0.048852217  2.605145 9.205081e-03
```

Han quedado una pocas interacciones

---

### Conclusiones elección de variables

- En general es muy difícil decidir qué variables y/o interacciones incluir en el modelo porque las posibilidades son muy grandes.

- No existe un método infalible y es importante "acotar" las variables/interacciones candidatas a partir de un criterio clínico.

- También se puede testar el efecto lineal de las variables cuantitativas.

---

### Selección de variables bajo un modelo causal

- Cuando ajustamos un modelo causal, no es tan importante que el modelo tenga buenas predicciones sino elegir las covariables (confusoras) adecuadas.

- Las variables confusoras son aquellas que están asociadas a la variable de interés y también a la variable respuesta.

- Aunque a veces se eligen un set de confusoras a priori (o varios sets) y se va jugando cómo cambia el coeficiente de la variable de interés.

---

- **Efecto del grupo sin ajustar**

```r
coef(summary(lm(bmi ~ group, predimed)))
```

```
                      Estimate Std. Error    t value     Pr(>|t|)
(Intercept)         30.2804407 0.08434088 359.024493 0.000000e+00
groupMedDiet + Nuts -0.5931884 0.11844958  -5.007940 5.650786e-07
groupMedDiet + VOO  -0.3399412 0.11734719  -2.896884 3.781828e-03
```

- **Ajustado por edad y sexo**

```r
coef(summary(lm(bmi ~ group + age + sex, predimed)))
```

```
                       Estimate  Std. Error    t value     Pr(>|t|)
(Intercept)         33.36688604 0.533482073  62.545468 0.000000e+00
groupMedDiet + Nuts -0.54241085 0.117165758  -4.629432 3.740016e-06
groupMedDiet + VOO  -0.33511405 0.115871435  -2.892120 3.839584e-03
age                 -0.03871097 0.007731921  -5.006643 5.688870e-07
sexMale             -1.20598508 0.096668819 -12.475430 2.641717e-35
```

---

- **Ajustado por edad, sexo y diabetes**

```r
coef(summary(lm(bmi ~ group + age + sex + diab, predimed)))
```

```
                       Estimate  Std. Error     t value     Pr(>|t|)
(Intercept)         33.36732146 0.533528330  62.5408617 0.000000e+00
groupMedDiet + Nuts -0.54198385 0.117198469  -4.6244960 3.829880e-06
groupMedDiet + VOO  -0.33548702 0.115898622  -2.8946592 3.808703e-03
age                 -0.03882887 0.007759897  -5.0037865 5.773549e-07
sexMale             -1.20777395 0.097180034 -12.4282108 4.707292e-35
diabYes              0.01729506 0.095547384   0.1810103 8.563653e-01
```

---

## Efecto dosis-respuesta

- Evaluación del efecto de una variable cuantitativa sobre la respuesta

- Los tests polinómicos pueden ser poco flexibles `$\Rightarrow$` métodos no paramétricos de splines.

- Para introducir términos splines hay que usar la función `gam`.

---

Efecto dosis-respuesta de la edad ajustando por el grupo de intervención y el sexo:

```r
library(gam)
modelo.gam <- gam(bmi ~ s(age) + group + sex, predimed, family="gaussian")
summary(modelo.gam)
```

```

Call: gam(formula = bmi ~ s(age) + group + sex, family = "gaussian", 
    data = predimed)
Deviance Residuals:
     Min       1Q   Median       3Q      Max 
-11.2839  -2.6824  -0.1445   2.4848  21.3781

(Dispersion Parameter for gaussian family taken to be 14.1506)

Null Deviance: 92182.36 on 6323 degrees of freedom
Residual Deviance: 89375 on 6316 degrees of freedom
AIC: 34713.78

Number of Local Scoring Iterations: NA

Anova for Parametric Effects
            Df Sum Sq Mean Sq F value    Pr(>F)    
s(age)       1    148  147.57  10.429  0.001247 ** 
group        2    387  193.56  13.679 1.181e-06 ***
sex          1   2233 2232.56 157.771 < 2.2e-16 ***
Residuals 6316  89375   14.15                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova for Nonparametric Effects
            Npar Df Npar F  Pr(F)
(Intercept)                      
s(age)            3 1.6315 0.1798
group                            
sex                              
```

---

```r
summary.glm(modelo.gam)
```

```

Call:
gam(formula = bmi ~ s(age) + group + sex, family = "gaussian", 
    data = predimed)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-11.2839   -2.6824   -0.1445    2.4848   21.3781

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         33.374150   0.533402  62.568  < 2e-16 ***
s(age)              -0.038793   0.007731  -5.018 5.36e-07 ***
groupMedDiet + Nuts -0.542011   0.117148  -4.627 3.79e-06 ***
groupMedDiet + VOO  -0.330666   0.115854  -2.854  0.00433 ** 
sexMale             -1.214047   0.096654 -12.561  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 14.15057)

Null deviance: 92182  on 6323  degrees of freedom
Residual deviance: 89375  on 6316  degrees of freedom
AIC: 34714

Number of Fisher Scoring iterations: NA
```

---

```r
plot(modelo.gam, terms="s(age)", se = TRUE)
```

![](bloque1_files/figure-html/unnamed-chunk-47-1.png)

---

# Modelos de regresión logística con R

---

Seguimos con la misma base de datos del predimed

- Ahora la variable respuesta es `event` (evento cardiovascular (sí/no))

- Y las variables independientes las mismas que en ejemplo anterior (`group`, `age`, `sex`, `smoke`, `diab` y `p14`)

- La función de R para ajustar modelos de regresión logística es **`glm`**, especificando el argumento `family=binomial` (o también `family=binomial()` o `family="binomial"`)

- La fórmula funciona de la misma manera que para la función `lm`. Pero a la izquierda de `~` tiene que haber una variable codificada como 0/1 ó `TRUE/FALSE`.

---

## Ajuste del modelo

- Tenemos que poner la variable `event` que está codificada como `Yes`, `No` a TRUE/FALSE, mediante la función `I()`

- Fíjate en el uso de `I(age/10)`. Ahora el coeficiente se refiere al cambio de 10 unidades de edad.

```r
modelo <- glm(I(event=='Yes') ~ group + I(age/10) + sex + smoke + diab + p14, 
              predimed, family=binomial)
```

---

## Extracción de resultados

---

### Sumario

```r
summary(modelo)
```

```

Call:
glm(formula = I(event == "Yes") ~ group + I(age/10) + sex + smoke + 
    diab + p14, family = binomial, data = predimed)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8195  -0.3142  -0.2428  -0.1881   3.0449

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -7.73196    0.80214  -9.639  < 2e-16 ***
groupMedDiet + Nuts -0.32354    0.16268  -1.989 0.046729 *  
groupMedDiet + VOO  -0.17937    0.15414  -1.164 0.244539    
I(age/10)            0.71293    0.10647   6.696 2.14e-11 ***
sexMale              0.48193    0.17304   2.785 0.005351 ** 
smokeCurrent         0.69454    0.20984   3.310 0.000933 ***
smokeFormer          0.47458    0.18348   2.587 0.009693 ** 
diabYes              0.58218    0.13651   4.265 2.00e-05 ***
p14                 -0.11432    0.03255  -3.513 0.000444 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2118.1  on 6323  degrees of freedom
Residual deviance: 1990.6  on 6315  degrees of freedom
AIC: 2008.6

Number of Fisher Scoring iterations: 6
```

---

### Odds Ratio e intervalos de confianza

```r
exp(cbind(coef(modelo), confint(modelo)))[-1,]
```

```
                                  2.5 %    97.5 %
groupMedDiet + Nuts 0.7235866 0.5245047 0.9935529
groupMedDiet + VOO  0.8357964 0.6170549 1.1301512
I(age/10)           2.0399665 1.6574556 2.5165602
sexMale             1.6191884 1.1536162 2.2733098
smokeCurrent        2.0027903 1.3216527 3.0115796
smokeFormer         1.6073357 1.1224756 2.3048405
diabYes             1.7899396 1.3730933 2.3463327
p14                 0.8919697 0.8370634 0.9510152
```

El OR de la edad es el OR de un cambio de 10 años de edad tal y como hemos especificado en la fórmula.

---

### Valores predichos

```r
#?predict.glm
predict(modelo, type="response")[1:10]
```

```
         1          2          4          5          6          8          9 
0.02223350 0.16424832 0.06685520 0.04965404 0.03529569 0.06118599 0.05231548 
        12         13         14 
0.04028074 0.03508895 0.17024069 
```

---

## Validación del modelo

- En regresión logística no hay la `$R^2$`.

- Para validar el modelo se calcula su capacidad de discriminación con el AUC.

- Y también se realiza el test de Hosmer-Lemeshow para ver si calibra bien (observados vs esperados)

---

### Discriminación

```r
library(pROC)
roc.res <- pROC::roc(modelo$y, fitted(modelo))
plot(roc.res)
```

![](bloque1_files/figure-html/unnamed-chunk-52-1.png)

```r
auc(roc.res)
```

```
Area under the curve: 0.7033
```

```r
ci(auc(roc.res))
```

```
95% CI: 0.6724-0.7342 (DeLong)
```

---

### Calibración

```r
library(ResourceSelection)
hl <- hoslem.test(modelo$y, fitted(modelo), g=10)
hl
```

```

Hosmer and Lemeshow goodness of fit (GOF) test

data:  modelo$y, fitted(modelo)
X-squared = 7.0185, df = 8, p-value = 0.5346
```

---

```r
tab <- rbind(hl$observed[,2], hl$expected[,2])
barplot(tab, beside=TRUE, names=1:10, ylab="observados")
```

![](bloque1_files/figure-html/unnamed-chunk-54-1.png)

---

## Otros

- Para regresión logística, la selección de las variables, interacciones, splines, etc. aplica de la misma manera que para los modelos de regresión lineal.

- Cuando se usa la función `gam` hay que especificar `family="binomial"`

---

# _Penalized regression models_

---

## _Penalized regression models_

- Modelos de regresión con una variable respuesta y muchas variables independientes

- Puede haber más variables que individuos

- No se pueden estimar los modelos "clásicos" de regresión

- Los _penalized regression models_ usan técnicas de estimación penalizada: compromiso entre ajuste y complejidad (número de variables).

---

## Modelo

- Criterio estimación de los coeficientes:

$$ \text{max}_{\beta} L(\beta) + \lambda * f(\beta)$$
- donde

- `$L(\beta)$` es la función de log-verosimilitud
  - `$\lambda$` es el parámetro de penalización
  - `$f(\beta)$` es la función de penalización de los parámetros
  
- Métodos según `$f(\beta)$`:
  
  - **_Lasso_** : `$f(\beta) = \sum_{k=1}^p |\beta_k|$`

- **_Ridge regression_**: `$f(\beta)  = \sum_{k=1}^p \beta_k^2 / 2$`

---

- Tanto con el método _lasso_ como _ridge regression_, cuando más grandes los efectos en valor absoluto más grande es la penalización.

- Los _penalized regression models_ se pueden considerar como técnicas de selección de variables y son útiles cuando a priori se sabe que algunas o la mayoría de las variables no son significativas.

- El método _ridge regression_ suele ser más conservador (elige menos variables).

---

## Elastic Net

- El método _Eslastic Net_ engloba tanto el método Lasso como Ridge regression.

`$$f(\beta) = (1-\alpha) \cdot \sum_{k=1}^p \beta_k^2 / 2 + \alpha \cdot \sum_{k=1}^p |\beta_k|$$`

- Cuando `$\alpha=1$` corresponde al método _lasso_ y cuando `$\alpha=0$` corresponde al _ridge regression_

- Función **`glmnet`** del package `glmnet` que puede estimar modelos _Elastic Net_ tanto para respuesta normal como binaria.

```r
library(glmnet)
# ?glmnet
```

---

## Ejemplo respuesta normal con R

```r
data(QuickStartExample)
dim(x)
```

```
[1] 100  20
```

```r
length(y)
```

```
[1] 100
```

---

```r
par(mfrow=c(1,2))
fit <- glmnet(x, y, alpha = 0)
plot(fit, xvar="lambda", label=TRUE); title("lasso")
fit <- glmnet(x, y, alpha = 1)
plot(fit, xvar="lambda", label=TRUE); title("ridge")
```

![](bloque1_files/figure-html/unnamed-chunk-57-1.png)

---

### Elección de `$\lambda$`

```r
set.seed(123456)
cvfit <- cv.glmnet(x, y, alpha=1)
cvfit
```

```

Call:  cv.glmnet(x = x, y = y, alpha = 1)

Measure: Mean-Squared Error

Lambda Measure     SE Nonzero
min 0.08307   1.073 0.1583       9
1se 0.15933   1.190 0.1867       8
```

---

```r
plot(cvfit)
```

![](bloque1_files/figure-html/unnamed-chunk-59-1.png)

---

### Coeficientes estimados al valor óptimo de `$\lambda$`

```r
coef(cvfit, s = "lambda.min")
```

```
21 x 1 sparse Matrix of class "dgCMatrix"
                      1
(Intercept)  0.14936467
V1           1.32975267
V2           .         
V3           0.69096092
V4           .         
V5          -0.83122558
V6           0.53669611
V7           0.02005438
V8           0.33193760
V9           .         
V10          .         
V11          0.16239419
V12          .         
V13          .         
V14         -1.07081121
V15          .         
V16          .         
V17          .         
V18          .         
V19          .         
V20         -1.04340741
```

---

### Predicciones

```r
predict(cvfit, newx = x[1:5,], s = "lambda.min")
```

```
              1
[1,] -1.3647490
[2,]  2.5686013
[3,]  0.5705879
[4,]  1.9682289
[5,]  1.4964211
```

***

**Ejercicio:** repite este ejemplo pero cambiando al método _lasso_.

---

## Ejemplo respuesta binaria con R

```r
data(BinomialExample)
dim(x)
```

```
[1] 100  30
```

```r
length(y)
```

```
[1] 100
```

```r
y[1:10]
```

```
 [1] 0 1 1 0 1 0 0 0 1 1
```

---

### Ajuste del modelo

```r
fit <- glmnet(x, y, family = "binomial", alpha=1)
plot(fit, xvar="lambda", label=TRUE)
```

![](bloque1_files/figure-html/unnamed-chunk-63-1.png)

---

### Selección de `$\lambda$`

```r
cvfit <- cv.glmnet(x, y, family = "binomial", type.measure = "class")
plot(cvfit)
```

![](bloque1_files/figure-html/unnamed-chunk-64-1.png)

---

### Coeficientes

```r
coef(cvfit, s = "lambda.min")
```

```
31 x 1 sparse Matrix of class "dgCMatrix"
                      1
(Intercept)  0.24460001
V1           0.02971466
V2           0.59025113
V3          -0.48659340
V4          -1.06018546
V5          -0.16440097
V6          -0.87042836
V7           .         
V8          -0.50010966
V9           0.68643360
V10         -1.32067765
V11         -0.03686667
V12         -0.04540248
V13          .         
V14          .         
V15          .         
V16          0.29231866
V17          .         
V18          .         
V19          .         
V20          .         
V21          .         
V22          0.18958230
V23          0.30256489
V24          .         
V25          0.60351837
V26         -0.31077805
V27         -0.07552772
V28          0.20812559
V29         -0.20239433
V30          0.07212866
```

---

### Predicciones

```r
predict(cvfit, newx = x[1:10,], s = "lambda.min", type = "class")
```

```
      1  
 [1,] "0"
 [2,] "1"
 [3,] "1"
 [4,] "0"
 [5,] "1"
 [6,] "0"
 [7,] "0"
 [8,] "0"
 [9,] "1"
[10,] "1"
```

```r
predict(cvfit, newx = x[1:10,], s = "lambda.min", type = "response")
```

```
               1
 [1,] 0.14974672
 [2,] 0.97467179
 [3,] 0.76584057
 [4,] 0.05882152
 [5,] 0.69260770
 [6,] 0.39180699
 [7,] 0.10592421
 [8,] 0.01240285
 [9,] 0.78809117
[10,] 0.80620636
```

---

## Ejercicio: Ejemplo con los datos del predimed

Ejecute este ejemplo e interpreta los resultados.

```r
y <- as.integer(predimed$event=='Yes')
x <- model.matrix(~ age + sex + group + diab + bmi, predimed)[,-1]

fit <- glmnet(x, y, family="binomial", alpha=1)
plot(fit, xvar="lambda", label=TRUE)

cvfit <- cv.glmnet(x, y, family = "binomial", 
                   lambda=seq(0,1,len=1000), type.measure = "class")
cvfit

coef(cvfit, s = "lambda.min")
```

---

## Conclusiones

- Los _penalized regression models_ ajustan modelos a los datos penalizando los coeficientes `$\beta_k$`.

- Los _elastic net_ engloban las técnicas _ridge regression_ y _lasso_

- No se obtienen p-valores ni intervalos de confianza.

- Sirven para elegir qué variables tienen un coeficiente diferente de cero.

---

# Tu turno

---

## Ejercicio 1

Con los datos `regicor` del package `compareGroups`

```r
data(regicor)
```

- Haz una descriptiva de las variables

- Ajusta un modelo lineal para predecir la variable `chol` a partir de las demás variables excepto: `id`, `cv`, `tocv`, `death`, `todeath`.

- Ajusta un modelo de regresión logística para predecir el evento cardiovascular `cv` a partir de todas las variables excepto: `id`, `tocv`, `death`, `todeath`.

Valida en ambos casos el modelo. Mira qué interacciones son significativas con el sexo.

---

## Ejercicio 2

Lee los datos `partoFin.sav` que hay en la sección datos del moodle

```r
library(haven)
datos <- read_sav("./datos/partoFin.sav")
```

- Lee la base de datos.

- Describe las variables.

- Construye un modelo para predecir el peso del bebé. Valídalo.

- ¿Cuál es el efecto (odds ratio) de la variable `tip_par` sin ajustar y ajustado por la edad de la madre y por la variable `tx`? Da también su intervalo de confianza.

---

## Ejercicio 3

- Analiza la base de datos `SNPs` del package `compareGroups`.

- Escoge mediante _Elastic Net_ los SNPs que predigan mejor el ser caso o control.

- Prueba para diferentes valores de `$\alpha$`.

```r
library(compareGroups)
library(glmnet)
data(SNPs)
names(SNPs)

# variable caso control
y <- SNPs$casco

# variables SNPs sin missings
x <- SNPs[,-(1:5)]
x <- x[,colSums(is.na(x))==0]

# conversión de los SNPs a valores enteros
for (j in 1:ncol(x)) x[,j] <- as.integer(as.factor(x[,j]))

x <- as.matrix(x)
```