03-aed-bi.Rmd

```{r}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, comment = "")
library(tidyverse)
library(knitr)
library(flextable)
library(gt)
library(gridExtra)
library(equatiomatic)
```

# Análisis exploratorio bivariante {#aed-bi}





## Datos bivariantes y multivariantes

El análisis univariante es muy útil para describir una única característica de la 
población en estudio. Pero rara vez estudiamos una característica aislada, y lo
habitual es tener conjuntos de datos con varias variables (cuantitativas y cualitativas)
que podemos estudiar por separado (análisis univariante) o conjuntamente (análisis multivariante).

El caso especial del **análisis bivariante** es cuando estudiamos dos características a la vez: $X, Y$.
Nos interesa la **relación** entre ellas, para lo que realizaremos resúmenes numéricos y gráficos.
Los datos bivariantes se encontrarán como pares de valores $(x_i, y_i)$ para cada
observación $i = 1, \ldots, n$.

Cuando estudiamos más de dos variables, tenemos datos multivariantes. En este caso,
estudiamos las relaciones "dos a dos" entre las variables (como en el caso bivariante) y 
la estructura conjunta. Hay algunas técnicas multivariantes específicas para este último caso.
En este capítulo nos vamos a centrar solo en el primer caso.

## Frecuencias conjuntas, marginales y condicionadas

El primer resumen que podemos hacer de datos bivariante es la tabla de frecuencias conjunta.
Igual que en el caso univariante $n$ es número total de datos, es decir, la frecuencia total.
La frecuencia absoluta conjunta, $n_{ij}$, es el número de observaciones en la clase $i$ de $X$ **y** en la clase $j$ de $Y$.
La frecuencia relativa conjunta es $f_{ij}= \frac{n_{ij}}{n}$.


### Distribución conjunta de frecuencias

Las frecuencias conjuntas se representan en una tabla de doble entrada, 
con los valores de una variable en filas y los de la otra en columnas.
En el interior, se ponen las frecuencias conjuntas (absolutas, relativas o ambas).
Si las dos variables son cualitativas, la tabla se denomina **Tabla de contingencia**.

```{r, echo=FALSE}
library(dplyr)
lab <- readxl::read_excel("_data/lab.xlsx") |> 
  mutate(fecha = as.Date(fecha))
```


:::{.rmdejemplo data-latex=""}
La Tabla \@ref(tab:tabs) muestra la tabla de contingencia de los analistas y el tipo de queso en el ejemplo del análisis de la producción de quesos. Asignamos la variable $X$ al Analista (en filas) y la variable $Y$ al Tipo (en columnas).
La Tabla \@ref(tab:trel) muestra la tabla de frecuencias relativas de los mismos datos. El número total de datos es
$n= `r nrow(lab)`$.
:::


```{r tabs, echo=FALSE}
kable(table(lab$analista, lab$tipo),
      caption = "Tabla de contingencia (frecuencias absolutas) de los analistas y el tipo de queso.")
```


```{r trel, echo=FALSE}
kable(prop.table(table(lab$analista, lab$tipo)), digits = 2, 
      caption = "Tabla de frecuencias relativas de los analistas y el tipo de queso.")
```




En el caso de variables continuas, debemos tener los datos agrupados en intervalos (clases).


:::{.rmdejemplo data-latex=""}

La tabla \@ref(tab:tcont) contiene las frecuencias absolutas conjuntas
de las variables:

* $X$ = ph (filas); 

* $Y$ = est (columnas)

:::



```{r, echo = FALSE}
histo <- hist(lab$ph, plot = FALSE)
histo2 <- hist(lab$est, plot = FALSE, breaks = 4)
lab <- lab |> 
  mutate(clase_ph = cut(ph, breaks = histo$breaks)) |> 
  mutate(clase_est = cut(est, breaks = histo2$breaks))
levels(lab$clase_ph) <- c(levels(lab$clase_ph), "(6.4,6.45]")
```



```{r tcont, echo=FALSE}
kable(table(lab$clase_ph, lab$clase_est), 
      caption = "Tabla de frecuencias conjunta del pH y el extracto seco total agrupadas en intervalos.")
```




### Distribución marginal

Si partimos de la distribución conjunta, podemos obtener la de cada
una de las variables (marginal) y estudiarla como datos univariantes.
Basta con hacer las sumas por filas (para la marginal de $X$), o por columnas, (para la marginal de $Y$):

* Frecuencias marginales de $X$: 
  + Absolutas: $n_{i\cdot} = \sum\limits_{j = 1}^{n_j}n_{ij}$
  + Relativas: $f_{i\cdot} = \sum\limits_{j = 1}^{n_j}f_{ij}$
  
* Frecuencias marginales de $Y$: 
  + $n_{\cdot j} = \sum\limits_{i = 1}^{n_i}n_{ij}$
  + $f_{\cdot j} = \sum\limits_{i = 1}^{n_i}f_{ij}$

donde $n_i$ es el número de clases de la variable $X$ y $n_j$ es el 
número de clases de la variable $Y$. Análogamente, para frecuencias marginales relativas sumamos las frecuencias relativas conjuntas,
o bien dividimos las frecuencias absolutas marginales entre el número total de datos $n$.




:::{.rmdejemplo data-latex=""}
La tabla \@ref(tab:tmar) contiene las frecuencias marginales como suma de filas y
columnas de la distribución conjunta en la tabla \@ref(tab:tcont). Las distribuciones
marginales de $X$ y de $Y$ por separado se muestran en las tablas \@ref(tab:tmarx) y \@ref(tab:tmary)
respectivamente.
:::


```{r tmar, echo=FALSE}
tabla <- table(lab$clase_ph, lab$clase_est) |> addmargins() 
tabla <- tabla |> 
  as.data.frame.matrix(tabla, row.names = rownames(tabla)) |> 
  rownames_to_column("X\\Y") 
tabla |> 
  gt(caption = "Frecuencias marginales como suma de filas y columnas de la distribución conjunta") |> 
  tab_style(style = list(cell_text(weight = "bold")),
            locations = cells_body(columns = Sum)) |> 
  tab_style(style = list(cell_text(weight = "bold")),
            locations = cells_body(rows = nrow(tabla)))
```



```{r tmarx, echo=FALSE}
tabla |> select(1, ncol(tabla)) |> 
  mutate(rel = round(Sum/nrow(lab), 3)) |> 
  gt(caption = "Frecuencias marginales de X (pH)") |> 
  cols_label(Sum = "\\(n_{i·}\\)",
             rel = "\\(f_{i·}\\)",
             `X\\Y` = "\\(x_i\\)")
```


```{r tmary, echo=FALSE}
tabla |> slice_tail() |> select(-1) |> t() |> as_tibble(rownames = "Y") |> 
  mutate(rel = round(V1/nrow(lab), 3)) |> 
  gt(caption = "Frecuencias marginales de Y (est)") |> 
  cols_label(Y = "\\(y_j\\)",
             V1 = "\\(n_{·j}\\)",
             rel = "\\(f_{·j}\\)")
```




### Distribución condicionada

La distribución de la variable $Y$ condicionada a un valor $x_i$ de la variable $X$ 
se representa por $Y | X = x_i$. Análogamente se puede definir para varios valores
de $X$, o de $X$ condicionada a $Y$.
Se lee "Distribución de $Y$ condicionada a que $X$ es igual a $x_i$".
Estas distribuciones condicionadas son variables univariantes que se pueden estudiar 
con análisis univariante.

A partir de la distribución conjunta, tomaríamos la fila o columna que se corresponde
con el valor "conocido", es decir, el de la condición. Las frecuencias relativas se calcularán 
dividiendo entre la frecuencia marginal, y no entre la frecuencia total. Por ejemplo:

$$f_{x_i|y=y_j}=\frac{n_{ij}}{n_{·j}}.$$



:::{.rmdejemplo data-latex=""}
La tabla \@ref(tab:tcondx) muestra la distribución De $X$ condicionada a que 
$Y = y_2$. Es decir, como es una variable numérica que hemos dividido en intervalos,
condicionada a que $Y \in `r colnames(tabla)[3]`$. La tabla \@ref(tab:tcondy)
muestra la distribución de $Y$ condicionada a que
$X = x_5$, es decir, de $Y | X \in `r tabla[5,1]`$.
:::



```{r tcondx, echo=FALSE}
tabla |> select(1, 3) |> rename(X = 1, n = 2) |> 
  mutate(f = round(n/602, 3)) |>
  gt(caption = "Frecuencias condicionadas de X") |>
  cols_label(X = "\\(x_i\\)",
             n = "\\(n_{i|j=2}\\)",
             f = "\\(f_{i|j=2}\\)")
```

```{r tcondy, echo=FALSE}
tt <- tabla |> slice(5) |> t() |>  as.data.frame() |> rownames_to_column("y")
# colnames(tt) <- tt[1,] 
tt <- tt |> slice(-1) |> mutate(V1 = as.numeric(V1))
tt |> 
  mutate(f = round(V1/184, 3)) |> 
  gt(caption = "Frecuencias condicionadas de Y") |> 
  cols_label(y = "\\(y_{j}\\)",
             V1 = "\\(n_{j|i=5}\\)",
             f = "\\(f_{j|i=5}\\)")

```



### Independencia de variables

Cuando la distribución de una variable $X$ no varía en función
de los valores de otra variable $Y$, entonces las variables $X$ e
$Y$ son estadísticamente independientes. Equivalentemente, la distribución de $X | Y=y_j$ es igual para
cualquier valor de $y_j$. Es decir, si son independientes se cumple:

$$f_{ij} = f_{i.}\cdot f_{.j} \;\forall i, j.$$


## Representación gráfica conjunta

### Gráficos de barras

Al igual que en el caso univariante, la mejor representación de tablas de
frecuencias para variables cualitativas y para variables numéricas discretas
sigue siendo el gráfico de barras. Hay varias variantes para representar las
dos variables:

1. El eje horizontal del gráfico representa las clases de ambas variables
$X$ e $Y$ de las variables, y se representan barras que representan las frecuencias absolutas conjuntas.

2. El eje horizontal del gráfico representa las clases de una de las variables,
la altura total de las barras representa la frecuencia absoluta de esa variable, y cada barra
se divide en trozos cuya altura es la frecuencia absoluta conjunta.

3.  El eje horizontal del gráfico representa las clases de una de las variables,
la altura total de las barras representa el 100% (todas la misma altura), y cada barra
se divide en trozos cuya altura es la frecuencia relativa condicionada a la clase de la barra.



:::{.rmdejemplo data-latex=""}
La figura \@ref(fig:barras2) muestra la representación de la tabla de
frecuencias \@ref(tab:tabs) con las barras adyacentes, mientras que la figura
\@ref(fig:barras2api) la representa con las barras apiladas. 
La figura \@ref(fig:barras2cond) muestra la representación de la tabla de
frecuencias de la variable "analista" condicionadas a cada tipo de queso.
:::



```{r barras2, fig.width=7, fig.height=4, fig.align='center', out.width="80%", fig.cap="Gráfico de barras con frecuencias conjuntas", echo=FALSE}
lab |> 
  ggplot(aes(x = tipo, fill = analista)) +
  geom_bar(position = position_dodge2()) +
  theme_bw() +
  labs(x = "Tipo de queso",
       y = "Frecuencia absoluta",
       title = "Gráfico de barras tabla de frecuencias",
       fill = "Analista") 
  
```

```{r barras2api, fig.width=7, fig.height=4, fig.align='center', out.width="80%", fig.cap="Gráfico de barras con frecuencias conjuntas apiladas", echo=FALSE}
lab |> 
  ggplot(aes(x = tipo, fill = analista)) +
  geom_bar(position = position_stack()) +
  theme_bw() +
  labs(x = "Tipo de queso",
       y = "Frecuencia absoluta",
       title = "Gráfico de barras tabla de frecuencias",
       fill = "Analista") 
  
```

```{r barras2cond, fig.width=7, fig.height=4, fig.align='center', out.width="80%", fig.cap="Gráfico de barras con frecuencias condicionadas", echo=FALSE}
lab |> 
  ggplot(aes(x = tipo, fill = analista)) +
  geom_bar(position = position_fill()) +
  theme_bw() +
  labs(x = "Tipo de queso",
       y = "Frecuencia relativa",
       title = "Gráfico de barras tabla de frecuencias condicionadas",
       fill = "Analista") 
  
```

Los gráficos de sectores son una representación muy popular para representar
tablas de frecuencias, aunque en general no están recomendados ya que el ojo humano
no es muy bueno para apreciar las diferencias entre ángulos, en comparación
con la capacidad para distinguir alturas.

Las tablas de frecuencias conjuntas de variables continuas agrupadas en intervalos
se pueden representar también mediante gráficos de barras igual que las 
cualitativas y las discretas. No obstante, si disponemos de todos los datos,
es más adecuado utilizar las representaciones que se explican a continuación.


### El gráfico de dispersión

Para dos variables continuas, la representación más adecuada es el gráfico de 
dispersión. El gráfico de dispersión es una "nube  de puntos" representada
en unos ejes cartesianos $X$, $Y$, que representan las escalas de cada una
de las variables. Cada punto representa un par de valores $(x_i, y_i)$.
Este gráfico permite apreciar a simple vista si hay relación lineal o de otro tipo
entre ambas variables.


:::{.rmdejemplo data-latex=""}
La figura \@ref(fig:disp) representa la variable "est" (extracto seco total) frente 
al pH en el ejemplo de los quesos. Se ha añadido transparencia a los puntos para
apreciar dónde hay mayor densidad de datos. A simple vista no se ve una relación clara.
Si añadimos algunos elementos al gráfico, podemos descubrir algo más. El gráfico
de la izquierda de la figura \@ref(fig:disp2) añade una línea de tendencia suavizada
en la que parece que hay un intervalo del pH en el que el est crece en la misma dirección
(entre 6,5 y 6,6 aproximadamente). En el gráfico de la derecha de la misma figura \@ref(fig:disp2)
hemos añadido una "capa" de color representando el tipo de queso. Aquí se puede ver claramente
como la "masa" inferior de puntos se corresponde con quesos distintos a los de la
"masa" superior.
:::


```{r disp, fig.cap="Gráfico de dispersión del ejemplo de los quesos", echo=FALSE}
lab |> 
  ggplot(aes(x = ph, y = est)) +
  geom_point(alpha = 0.6) +
  theme_bw() +
  labs("Gráfico de dispersión del extracto seco total frente al pH",
       x = "pH",
       y = "Extracto seco total")
```


```{r disp2, echo=FALSE, fig.height = 4, fig.width=9, out.width="100%", fig.cap="Gráficos de dispersión 'enriquecidos'"}

p1 <- lab |> 
  ggplot(aes(x = ph, y = est)) +
  geom_point(alpha = 0.6) +
  geom_smooth() +
  labs(title = "Gráfico de dispersión con ajuste de regresión",
       x = "pH",
       y = "Extracto seco total") +
  theme_bw()
p2 <- lab |> 
  ggplot(aes(x = ph, y = est, col = tipo)) +
  geom_point(alpha = 0.6) +
  labs(title = "Gráfico de dispersión identificando grupos",
       x = "pH",
       y = "Extracto seco total",
       colour = "Tipo") +
  theme_bw()

grid.arrange(p1, p2, nrow = 1, widths = c(0.4, 0.6))
```


### Gráficos de cajas por grupos

Cuando queremos estudiar conjuntamente una variable cuantitativa y una 
variable cualitativa, la representación más adecuada es mediante gráficos 
de cajas representando en el eje horizontal las clases de la variable 
cualitativa. De esta forma podemos comparar los principales estadísticos
en cada uno de los grupos. Los gráficos de cajas pueden ocultar información importante acerca de la distribución, por lo que otra representación muy útil
es el gráfico de densidad por grupos. También se pueden añadir capas de puntos,
o utilizar el gráfico de violín que suaviza la caja a la forma de la distribución.


:::{.rmdejemplo data-latex=""}
En la figura \@ref(fig:bp2) se representa la relación entre el tipo de queso
y el pH. El gráfico de la izquierda muestra el gráficos de cajas, y el de la derecha, un gráfico de violín. La figura \@ref(fig:dens2) mustra gráficos de 
densidad de la variable pH para cada tipo de queso. En ambos casos se aprecia claramente el comportamiento diferente de los quesos de tipo C.
:::




```{r bp2, echo=FALSE, fig.height = 4, fig.width=9, out.width="100%", fig.cap="Representación de la relación entre pH y tipo de queso"}
p1 <- lab |> 
  ggplot(aes(x = tipo, y = ph)) + 
  geom_boxplot() +
  theme_bw() +
  labs(x = NULL, y = NULL)
p2 <- lab |> 
  ggplot(aes(x = tipo, y = ph)) + 
  geom_violin() +
  geom_jitter(aes(col = analista), alpha = 0.2) + 
  theme_bw() +
  labs(x = NULL, y = NULL, colour = "Analista")
grid.arrange(p1, p2, nrow = 1, 
             top = "pH por tipo de queso",
             bottom = "Tipo de queso",
             left = "pH", widths = c(0.4, 0.6))
```

```{r dens2, fig.cap="Gráficos de densidad por grupos", fig.height=4}
lab |> 
  ggplot(aes(x = ph, colour = tipo, fill = tipo)) +
  geom_density(alpha = 0.4) +
  theme_bw() +
  labs(title = "Gráficos de densidad por tipo de queso",
       x = "pH",
       y = "Densidad",
       fill = "Tipo de queso",
       colour = "Tipo de queso") 
  
```





## Medidas conjuntas

Al igual que en el caso univariante, podemos obtener medidas conjuntas
de dos variables que resuman su relación. 

### Covarianza

La medida más importante es 
la **covarianza**. La covarianza mide la dependencia **lineal** entre las variables, y se define como:

$$\sigma_{xy} = \frac{1}{N} \sum\limits_{i=1}^n(X_i-\bar X)(Y_i-\bar Y)$$

Para datos muestrales, dividimos por $n - 1$:

$$s_{xy} = \frac{1}{n-1} \sum\limits_{i=1}^n(x_i-\bar x)(y_i-\bar y).$$



* Si $s_{xy} > 0$: Relación lineal positiva (cuando una crece, la otra también).
* Si $s_{xy} < 0$: Relación lineal negativa (cuando una crece, la otra decrece).
* Si $s_{xy} = 0$: No hay relación lineal.


Decimos que hay relación lineal cuando una variable se puede expresar en función de la otra como una transformación lineal, es decir: $Y=a+bX$. Esta sería una relación perfecta entre dos variables. En la práctica, habrá cierto "error" en esta relación, como veremos más adelante.

Podemos realizar el cálculo abreviado al igual que con la varianza. Según el caso:

* Covarianza poblacional con todos los datos:

$$\sigma_{xy} = \frac{1}{N} \sum\limits_{i=1}^N(X_i \cdot Y_i) - \bar X \cdot \bar Y.$$

* Covarianza poblacional con tablas de frecuencias:

$$\sigma_{xy} = \frac{1}{N} \sum\limits_{i=1}^N(n_i \cdot X_i \cdot Y_i) - \bar x \cdot \bar y.$$

* Covarianza muestral, con todos los datos:

$$s_{xy} = \frac{1}{n-1} \left ( \sum\limits_{i=1}^n(x_i \cdot y_i) - n \cdot \bar x \cdot \bar y \right ).$$


Una propiedad importante de la covarianza es que si $X$ e $Y$ son variables independientes, su covarianza es cero. Lo contrario no tiene por qué ser verdad, ya que puede existir una relación no lineal, como veremos más adelante.

La covarianza es un estadístico muy importante para aplicar muchas técnicas estadísticas. Sin embargo, su interpretación más allá del signo no es muy útil, ya que no está acotada: depende de las unidades de las variables. 

### Correlación

Para interpretar la relación lineal entre dos variables, en vez de la covarianza calcularemos el **coeficiente de correlación lineal**, que sí está acotado, en concreto entre -1 y 1:

$$r_{xy}=\frac{s_{xy}}{s_x \cdot s_y}.$$

* Si $r_{xy}$ próximo a 0, no hay relación lineal
* Si $r_{xy}$ próximo a -1 o 1, relación lineal **fuerte**, negativa o positiva respectivamente
* Si $r_{xy} = 1$ relación lineal perfecta. Aparece también al correlacionar una variable con sí misma.


Un conjunto de datos bivariantes quedará caracterizado por su vector de medias
y su matriz de varianzas-covarianzas. El vector de medias contiene las medias
de las dos variables:

$$\left [
\begin{array}{c}\bar x\\
\bar y\end{array}\right ].$$


La Matriz de varianzas-covarianzas, o simplemente matriz de covarianzas, sería:

$$\mathbf{S} = \left [\begin{array}{cc}
s^2_x & s_{xy}\\
s_{xy} & s_y^2
\end{array}\right ]$$

Mientras que la matriz de correlaciones tendría valores igual a 1 en la diagonal,
y el coeficiente de correlación fuera de la diagonal:

$$\mathbf{R} = \left [\begin{array}{cc}
1 & r_{xy}\\
r_{xy} & 1
\end{array}\right ]$$




:::{.rmdejemplo data-latex=""}
En el ejemplo de la producción de quesos,
si tomamos $X$ = ph; $Y$ = est, tenemos los siguientes datos resumidos:

$$n = `r fnum(nrow(lab))`;\; \sum x_i = `r fnum(sum(lab$ph))`;\; \sum y_i = `r fnum(sum(lab$est))`$$
$$\sum x_i^2 = `r fnum(sum(lab$ph^2))`;\; \sum y_i^2 = `r fnum(sum(lab$est^2))` ;\; \sum x_i y_i = `r fnum(sum(lab$ph*lab$est))`$$


el vector de medias sería:

$$\left [
\begin{array}{c} `r fnum(mean(lab$ph))`\\
`r fnum(mean(lab$est))`\end{array}\right ],$$

y la matriz de covarianzas (muestrales)

$$\mathbf{S} = \left [\begin{array}{cc}
s^2_x & s_{xy}\\
s_{xy} & s_y^2
\end{array}\right ]=\left [\begin{array}{cc}
`r fnum(var(lab$ph), 4)` & `r fnum(var(lab$ph, lab$est), 4)`\\
`r fnum(var(lab$ph, lab$est), 4)` & `r fnum(var(lab$est), 4)`
\end{array}\right ]$$

El coeficiente de correlación lineal será entonces:

$$r_{xy}=\frac{s_{xy}}{s_x \cdot s_y}= \frac{`r fnum(var(lab$ph, lab$est),4)`}{\sqrt{`r fnum(var(lab$ph),4)`\cdot `r fnum(var(lab$est),4)`}} = `r fnum(cor(lab$ph, lab$est),4)`.$$
:::

### Caso multivariante

En el caso de estudiar más de dos características numéricas, se extiende la matriz de 
covarianzas a cualquier número de variables. En este caso las matrices de covarianzas y de correlaciones tendrán tantas filas o columnas como variables. Estas matrices se utilizan en cálculo multivariante y multitud de métodos estadísticdos. Una propiedad importante de la matriz de covarianzas es que si previamente **tipificamos** las variables, la matriz de covarianzas de la variable tipificada es igual a la matriz de correlaciones de las variables originales.

:::{.rmdpractica data-latex=""}
Con el código a continuación obtenemos la matriz de correlaciones de todas
las variables numéricas del conjunto de datos de la producción de queso.
:::


```{r}
lab |> select(where(is.numeric), -codigo) |> 
  drop_na() |> 
  cor() |> 
  round(2)
```



### Patrones en los gráficos

La nube de puntos del gráfico de dispersión nos da una idea de la relación entre 
las variables. También nos da pistas de cómo será el coeficiente de correlación.
Es importante volver a resaltar que Los resúmenes numéricos pueden estar escondiendo información importante, por ejemplo relaciones no lineales. La figura \@ref(fig:patronesr) muestra los gráficos de dispersión de dos variables
y sus correspondientes coeficientes de correlación. Los dos primeros gráficos de la primera columna muestran una relación perfecta, mientras que el tercero muestra ausencia de correlación, y además independencia (no se aprecia otro tipo de relación, solo "ruido"). En la segunda columna, el primer gráfico muestra una alta correlación lineal, es decir hay una relación **fuerte** entre las variables. El segundo sin embargo también muestra relación lineal, pero más moderada. El último gráfico ilustra cómo una correlación muy baja, próxima a cero, indica que no hay relación lineal, como en el gráfico de la izquierda. Si embargo, no podemos decir que las variables sean independientes, porque hay claramente una dependencia funcional entre la $X$ y la $Y$, pero no es lineal.




```{r patronesr, echo=FALSE, fig.width=7, fig.height=6, fig.align='right', out.width="100%", fig.cap="Patrones en la nube de puntos según el coeficiente de correlación"}
set.seed(2)
micolor <- "orange"
x1  <-  runif(100, 6, 7)
x7 <- runif(100, -1, 1)
cdata <- data.frame(
  x1 = x1,
  x2 = 5 + 0.5*x1,
  x3 = 5+2*x1 + rnorm(100, sd = 0.2),
  x4 = sample(rnorm(100, 100, 4)),
  x5 = 5 - 0.2*x1 ,
  x6 = 5 - 0.2*x1 + rnorm(100, mean = 10, sd = 0.06),
  x7,
  x8 = -0.8*x7^2 + 2 + rnorm(100, sd = 0.02))
p1 <- cdata |> ggplot(aes(x1, x2)) + geom_point(col = micolor) + theme_bw() +
  labs(title = paste0("r = ", cor(cdata$x1, cdata$x2)))
p2 <- cdata |> ggplot(aes(x1, x3)) + geom_point(col = micolor) + theme_bw() +
  labs(title = paste0("r = ", round(cor(cdata$x1, cdata$x3), 4)))
p3 <- cdata |> ggplot(aes(x1, x4)) + geom_point(col = micolor) + theme_bw() +
  labs(title = paste0("r = ", round(cor(cdata$x1, cdata$x4), 4)))
p4 <- cdata |> ggplot(aes(x1, x5)) + geom_point(col = micolor) + theme_bw() +
  labs(title = paste0("r = ", round(cor(cdata$x1, cdata$x5), 4)))
p5 <- cdata |> ggplot(aes(x1, x6)) + geom_point(col = micolor) + theme_bw() +
  labs(title = paste0("r = ", round(cor(cdata$x1, cdata$x6), 4)))
p6 <- cdata |> ggplot(aes(x7, x8)) + geom_point(col = "red3") + theme_bw() +
  labs(title = paste0("r = ", round(cor(cdata$x7, cdata$x8), 4)))


library(gridExtra)
grid.arrange(p1, p2,  p4, p5, p3, p6)
```






## Regresión lineal simple

Las técnicas de regresión, de manera general tratan de encontrar una
relación funcional entre una variable dependiente $Y$ y una o más variables
independientes $X$:

$$Y = f(\mathbf{X}) + \varepsilon.$$

En el caso más sencillo de una sola variable independiente, cuando existe relación lineal los puntos del gráfico de dispersión se disponen alrededor de una **línea recta**, y podemos expresar la relación mediante la ecuación de la recta que mejor
se aproxima a la nube de puntos. Esta aproximación estará ligada a un error, representado en la expresión anterior por $\varepsilon$, que mediremos con la variabilidad alrededor de la línea, de forma análoga a como la varianza expresa la variabilidad de los datos entorno a la media. Así, la recta de regresión quedaría expresada como:

$$Y = a + bX + \varepsilon,$$
donde $X$ es la variable independiente o explicativa, $Y$ es la variable dependiente o variable respuesta, $a$ y $b$ son los coeficientes de regresión (término independiente y pendiente de la recta respectivamente) y $\varepsilon$ es el error. 

### Estimación de los coeficienetes

Tenemos pares de datos $(x_i, y_i)$ que queremos ajustar al modelo matemático de la recta. Para cualesquiera valores de $a$ y $b$ que elijamos, podemos
aproximar los valores de $y_i$ como:

$$\hat{y}_i = a + bx_i.$$

Esta aproximación tiene un error, que es el **residuo** de esa observación:

$$y_i = \hat y + \varepsilon_i \implies \varepsilon_i=y_i - \hat{y}_i.$$
La figura \@ref(fig:resid) esquematiza la relación entre los valores de la
variable bivariante (puntos naranjas), la recta de regresión (línea azul) y
los residuos (líneas negras).

```{r resid, echo = FALSE, fig.cap="Esquematización de los valores, recta y residuos en regresión lineal"}
set.seed(2)
micolor <- "orange"
x1  <-  runif(100, 6, 7)
x7 <- runif(100, -1, 1)
cdata <- data.frame(
  x1 = x1,
  x2 = 5 + 0.5*x1,
  x3 = 5+2*x1 + rnorm(100, sd = 0.2),
  x4 = sample(rnorm(100, 100, 4)),
  x5 = 5 - 0.2*x1 ,
  x6 = 5 - 0.2*x1 + rnorm(100, mean = 10, sd = 0.06),
  x7,
  x8 = -0.8*x7^2 + 2 + rnorm(100, sd = 0.02))

cdata2 <- lm(x3~x1, cdata |> slice_head(n =10) ) |> broom::augment() 
cdata2 |> 
  ggplot(aes(x1, x3)) + geom_point(col = micolor) + theme_bw() +
  labs(title = "Valores, recta y residuos") +
  geom_smooth(method = "lm", se = FALSE) +
  geom_segment(aes(x= x1, xend = x1, y = x3, yend = .fitted))

```



Una forma de estimar los coeficientes de regresión es minimizando la suma de los
cuadrados de los residuos:

$${a; b}: \min \sum\limits_{i=1}^n \varepsilon_i^2 =\min \sum\limits_{i=1}^n(y_i-(a+bx_i))^2.$$
Derivando e igualando a cero tenemos los siguientes valores que minimizan la función:

$$b = \frac{s_{xy}}{s_x^2},$$
$$a = \bar y - b \bar x.$$
Es inmediato demostrar que la pendiente de la recta también la podemos
calcular como:

$$b = \frac{s_y}{s_x}r_{xy}.$$

### Interpretación de los coeficientes

La recta de regresión nos sirve para dos cometidos: explicar la relación entre
las variables, y hacer predicciones. En cuanto a la interpretación, el coeficiente $a$ es el valor medio esperado de $Y$ cuando $X=0$. 
No siempre tiene interpretación ya que el 0 puede estar fuera del dominio de $X$.
El coeficiente $b$ significa "cuántas unidades aumenta $Y$ por cada unidad 
que aumenta $X$ (si $b<0$, cuánto disminuye).

En cuanto a la predicción, podemos predecir cuánto valdrá, en promedio, la variable $Y$ para un
determinado valor de $X, x_{n+1}$. Para hacer la predicción, sustuimos $x$ por el nuevo valor en la ecuación de la recta:

$$\hat{y}_{n+1} = a + bx_{n+1}.$$
Hay que tener en cuenta que esta predicción solo es válida para el rango de $X$ utilizado para obtener los coeficientes.



### Error y bondad de ajuste

Pare resumir el error cometido, podríamos buscar una medida resumen de los residuos. Por las propiedades del método de mínimos cuadrados, la media de los 
residuos es cero:

$$\bar \varepsilon = \frac{1}{n}\sum\limits_{i=1}^n \varepsilon_i = 0.$$

Entonces, para cuantificar el error, calculamos la media de los residuos al cuadrado, a
la que llamamos "varianza residual":

$$s_\varepsilon^2= \frac{1}{n}\sum\limits_{i=1}^n \varepsilon_i^2.$$

Si desarrollamos los términos de la varianza residual: 

$$s_\varepsilon^2= \frac{1}{n}\sum\limits_{i=1}^n (y_i - (a + b x_i))^2,$$
llegamos a la siguiente expresión:

$$\frac{s_\varepsilon^2}{s_y^2}=(1-r_{xy}^2),$$
y de ahí:

$$r^2_{xy} = 1- \frac{s_\varepsilon^2}{s_y^2} \equiv R^2,$$

que definimos como **coeficiente de determinación**, y su interpretación es como sigue. La varianza residual $s_\varepsilon^2$ es la variabilidad no explicada por la recta. Cuanto más se acerquen los puntos a la recta, mejor explicado está el modelo, y más pequeña es esta varianza. La varianza de $Y$, $s_y^2$, es la variabilidad Total. Por tanto, $\frac{s_\varepsilon^2}{s_y^2}$ es la proporción de variabilidad no explicada por la recta de regresión. Y entonces $R^2$ es **la proporción de variabilidad explicada por la regresión**. 



La expresión general $R^2 = 1- \frac{s_\varepsilon^2}{s_y^2}$ es válida para cualquier modelo de regresión en el que podamos calcular los residuos y su varianza. Pero la equivalencia $R^2=r^2$ solo es válida en el caso de la regresión lineal simple.

```{r, echo=FALSE}
ldata <- lab |> slice(1:15) |> select(mg, est) 
mod <- lm(est ~ mg, ldata)
mod.a <- mod |> broom::augment()
mod.m <- mod |> broom::glance()
mod.c <- mod |> broom::tidy()
# summary(mod)
```


:::{.rmdejemplo data-latex=""}
En el ejemplo de la fábrica de quesos, tomemos $X$ la variable "materia grasa" e $Y$ el extracto seco total. Para simplificar el ejemplo, tomamos solo las 15 primeras observaciones, que se muestran en la tabla \@ref(tab:reglabt). La figura \@ref(fig:reglab) muestra la nube de puntos
y la recta de regresión estimada por mínimos cuadrados. La recta de regresión estimada es:

```{r, echo=FALSE}
extract_eq(mod, use_coefs = TRUE, coef_digits = 3)
```

El coeficiente $a$ no tiene interpretación, ya que el cero no es un valor
posible en la variable materia grasa, al menos en esta muestra. La interpretación del coeficiente $b$ es la siguiente: por cada unidad que aumente `mg`, porcentaje de materia grasa, aumenta `r round(coef(mod)[2], 3)` unidades `est` (extracto seco total).

La varianza residual es `r round(mod.m$sigma^2, 3)`, y el coeficiente de determinación:

$$R^2 = `r fnum(mod.m$r.squared)`.$$

Por tanto, la regresión explica el 68,48% de la variabilidad de los datos.


:::

```{r reglabt, echo=FALSE}
ldata |> t() |> kable(digits = 2, caption = "Datos de quesos para la regresión")
```


```{r reglab, echo = FALSE, fig.height=4, fig.align='center', out.width="100%", fig.cap="Recta de regresión del extracto seco frente a la materia grasa"}

ldata |> ggplot(aes(x=mg, y = est)) +
  geom_point() + theme_bw() +
  geom_smooth(method = lm, se = FALSE) +
  labs(title = "Recta de regresión ejemplo quesos",
       x = "Materia grasa",
       y = "Extracto seco total")
```

### Extensiones del modelo lineal

El modelo lineal explicado aquí se puede extender fácilmente para recoger algunas
relaciones no lineales. Una manera sencilla es aplicar logaritmos a la variable explicativa, a la respuesta o a ambas, y estimar los coeficientes con las variables transformadas. Del mismo modo, se pueden transformar usando potencias o exponenciales, y ajustar polinomios añadiendo términos cuadráticos, cúbicos, etc.

Otra extensión es la regresión multivariante, generalizando a $p$ variables independientes, que incluso pueden ser cualitativas. Por último, los modelos
lineales generalizados pueden estimar modelos con respuestas no continuas, por 
ejemplo binarias (regresión logística) o de recuentos (regresión de Poisson).


:::{.rmdcafe data-latex=""}
¡Correlación no implica **causalidad**!

El hecho de que exista relación entre las variables, no significa
que esta relación sea de causa-efecto. Para demostrar este tipo de
relaciones, sería necesario confirmar con diseño de experimentos.

Visita este enlace para ver algunas correlaciones espúreas (¡y curiosas!).

https://tylervigen.com/spurious-correlations


:::



## Intro multivariante