class: center, middle # TEMA 1: Descripción de datos ### José R. Berrendero #### Departamento de Matemáticas, Universidad Autónoma de Madrid --- # Temas a tratar **Objetivos generales** - Introducir las técnicas básicas para describir y visualizar un conjunto de datos - Introducir herramientas de software útiles para estos propósitos. **Objetivos específicos** - Manejo básico de los datos: `dplyr` - Visualizacion de los datos: `ggplot2` - Descripción (grafica y numérica) de la distribución de una variable - Descripcion (gráfica y numérica) de la relación entre variables --- # Datos ordenados. Tipos de variables **Para que un conjunto de datos esté ordenado** - Cada columna corresponde a la información de una variable - Cada fila corresponde a la información de un individuo **Nomenclatura** Las filas se llaman **observaciones**, el número de observaciones se denomina **tamaño muestral**, `\(n\)`, y el número de columnas es la **dimensión**, `\(p\)`. **Tipos de variables** - **variables cualitativas**: contienen información sobre cualidades o atributos, que no se expresan numéricamente. - **variables cuantitativas**: contienen información numérica. Pueden ser: * **Discretas**: toman un número pequeño de valores, normalmente enteros. * **Continuas**: pueden tomar cualquier valor de un intervalo. - **variables ordinales**: toman un número pequeño de valores que se pueden ordenar. --- # La estructura fundamental de datos en R La estructura básica para organizar los datos es un *data frame*: ficheros de datos con `\(n\)` observaciones y `\(p\)` variables. Un *data frame* es una matriz. Las columnas pueden ser de distinta naturaleza (números, palabras, fechas, etc.) pero todas tienen la misma longitud. Para crear un *data frame* con R se usa el comando `data.frame`: ```r x <- c(1, 2, 3, 4) y <- c("a", "b", "c", "d") datos <- data.frame(var1 = x, var2 = y) glimpse(datos) # para ver la estructura del fichero ``` ``` ## Rows: 4 ## Columns: 2 ## $ var1 <dbl> 1, 2, 3, 4 ## $ var2 <chr> "a", "b", "c", "d" ``` --- # Concentraciones de mercurio en peces > Se ha realizado un estudio en los ríos Wacamaw y Lumber en Carolina del Norte (EE.UU.), analizando la cantidad de mercurio que contenían 171 ejemplares capturados de una cierta especie de peces. Leemos los datos de un fichero de texto y los convertimos en un *data frame*: ```r enlace <- "http://matematicas.uam.es/~joser.berrendero/datos/mercurio.txt" mercurio <- read.table(enlace, header = TRUE) head(mercurio) # Primeras filas del fichero ``` ``` ## RIO ESTACION LONG PESO CONC ## 1 0 0 47.0 1616 1.60 ## 2 0 0 48.7 1862 1.50 ## 3 0 0 55.7 2855 1.70 ## 4 0 0 45.2 1199 0.73 ## 5 0 0 44.7 1320 0.56 ## 6 0 0 43.8 1225 0.51 ``` ¿De qué tipo es cada una de las variables? ¿Cuántas observaciones hay? ¿De qué dimensión? --- # Manejo básico de los datos Vamos a usar [el paquete `dplyr`](https://dplyr.tidyverse.org/) para manejar los datos. Está estructurado en una serie de *verbos* que corresponden a las manipulaciones más habituales. Incorpora una *sintaxis en cadena (pipes*) que sigue este esquema: ```r nuevo_data_frame <- antiguo_data_frame %>% verbo1("argumentos") %>% verbo2("argumentos") %>% ... ``` Vamos a ver los *verbos* que más se usan para modificar los datos. --- # Mutate **Crear nuevas variables que resultan de transformar las existentes** En el siguiente ejemplo se recodifica la variable `RIO`, y se declara tanto esta variable como la variable `ESTACIÓN` como variables cualitativas (en R, *factores*): ```r mercurio <- mercurio %>% mutate(RIO = recode(RIO, "0" = "Lumber", "1" = "Wacamaw")) %>% mutate(RIO = factor(RIO)) %>% mutate(ESTACION = factor(ESTACION)) head(mercurio) ``` ``` ## RIO ESTACION LONG PESO CONC ## 1 Lumber 0 47.0 1616 1.60 ## 2 Lumber 0 48.7 1862 1.50 ## 3 Lumber 0 55.7 2855 1.70 ## 4 Lumber 0 45.2 1199 0.73 ## 5 Lumber 0 44.7 1320 0.56 ## 6 Lumber 0 43.8 1225 0.51 ``` --- # Count **Contar las veces que se repite cada valor (frecuencias)** ```r mercurio %>% count(RIO, ESTACION) ``` ``` ## RIO ESTACION n ## 1 Lumber 0 10 ## 2 Lumber 1 19 ## 3 Lumber 2 13 ## 4 Lumber 3 10 ## 5 Lumber 4 14 ## 6 Lumber 5 3 ## 7 Lumber 6 4 ## 8 Wacamaw 7 10 ## 9 Wacamaw 8 14 ## 10 Wacamaw 9 3 ## 11 Wacamaw 10 21 ## 12 Wacamaw 11 9 ## 13 Wacamaw 12 13 ## 14 Wacamaw 13 4 ## 15 Wacamaw 14 13 ## 16 Wacamaw 15 11 ``` --- # Filter **Filtrar las observaciones y quedarnos solo con las que cumplen cierta condición** Nos quedamos solo con las observaciones que corresponden al río Lumber: ```r mercurio_lumber <- mercurio %>% filter(RIO == "Lumber") head(mercurio_lumber) ``` ``` ## RIO ESTACION LONG PESO CONC ## 1 Lumber 0 47.0 1616 1.60 ## 2 Lumber 0 48.7 1862 1.50 ## 3 Lumber 0 55.7 2855 1.70 ## 4 Lumber 0 45.2 1199 0.73 ## 5 Lumber 0 44.7 1320 0.56 ## 6 Lumber 0 43.8 1225 0.51 ``` --- # Select **Seleccionar algunas de las variables** Retenemos solo las variables correspondientes a la longitud y la concentración de mercurio: ```r datos <- mercurio_lumber %>% select(LONG, CONC) head(datos) ``` ``` ## LONG CONC ## 1 47.0 1.60 ## 2 48.7 1.50 ## 3 55.7 1.70 ## 4 45.2 0.73 ## 5 44.7 0.56 ## 6 43.8 0.51 ``` [Más información sobre selección de variables](https://tidyselect.r-lib.org/reference/language.html) --- # Arrange **Ordenar los datos** Seleccionamos los peces del río *Lumber* cuya longitud está por encima de la media y los ordenamos de menor a mayor concentración de mercurio: ```r datos <- mercurio_lumber %>% filter(LONG > mean(LONG)) %>% arrange(CONC) head(datos) ``` ``` ## RIO ESTACION LONG PESO CONC ## 1 Lumber 1 45.1 2920 0.34 ## 2 Lumber 0 40.4 1033 0.50 ## 3 Lumber 0 43.8 1225 0.51 ## 4 Lumber 1 43.5 2674 0.54 ## 5 Lumber 0 44.7 1320 0.56 ## 6 Lumber 4 41.0 1042 0.66 ``` --- # Distribución de una variable Distribución: valores y frecuencias con la que se toman. Si los datos `\(x_1,\ldots,x_n\)` corresponden a `\(n\)` realizaciones independientes de una variable aleatoria entonces la distribución de los datos es una aproximación de la distribución de la variable aleatoria. **Frecuencia absoluta** de un valor es el número de observaciones para las que la variable toma ese valor. Se denotarán `\(N_1,\ldots,N_k\)`, donde `\(k\)` es el número de valores diferentes. **Frecuencia relativa** es la frecuencia absoluta dividida por `\(n\)`. Se denotarán `\(f_1,\ldots,f_k\)`. *Datos agrupados* en `\(k\)` intervalos o clases `\(A_1,\dots,A_k\)`. Los valores `\(m_1,\dots,m_k\)` que representan cada clase se llaman **marcas de clase**. De la distribución nos interesa conocer su posición, dispersión y forma. --- # Nociones básicas de `ggplot2` Para visualizar un conjunto de datos vamos a usar [el paquete `ggplot2`](https://ggplot2.tidyverse.org/) de R. ```r ggplot(data = <DATA>) + <GEOM_FUNCTION>(mapping = aes(<MAPPINGS>)) ``` - El comando `ggplot(data = <DATA>)` crea un sistema de coordenadas vacío. El argumento `data` determina el *data frame* donde están los datos. - Al sistema de coordenadas vacío creado se le van añadiendo *capas* usando el signo `+`. - El comando `<GEOM_FUNCTION>()` añade una capa en la que se determina el tipo de elementos que se van a representar. Catálogo de *geoms* en la [documentación del paquete](https://ggplot2.tidyverse.org/reference/index.html#section-layer-geoms). - El argumento `mapping = aes(<MAPPINGS>)` determina los llamados *aesthetics*, una lista en la que se asignan las variables del fichero a distintos aspectos del gráfico. --- # Gráfico de barras Para representar gráficamente la distribución de una variable cualitativa (o de una cuantitativa discreta) se usa un gráfico de barras [`geom_bar()`]: ```r ggplot(mercurio) + geom_bar(aes(x = ESTACION, fill = RIO), width = 0.5) + labs(x = "Estación", y = "Número de capturas") ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-11-1.png" style="display: block; margin: auto;" /> --- # Histogramas Para representar gráficamente la distribución de una variable cuantitativa continua se usa un histograma [`geom_histogram()`]: - Se divide el rango de los datos en un número adecuado de intervalos. - Sobre cada intervalo se dibuja un rectángulo cuya área es proporcional a la frecuencia (relativa o absoluta) de datos en el intervalo. --- # Histogramas ```r ggplot(mercurio) + geom_histogram(aes(x = CONC, y = after_stat(density)), # y = after_stat(density) para área = 1 bins = 10, # número de intervalos col = 'black', # color del borde fill = 'olivedrab4') + # color del interior labs(x = "Concentración (ppm)", y = NULL) # cambia etiquetas ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-12-1.png" style="display: block; margin: auto;" /> --- # Histogramas A tener en cuenta al mirar un histograma: - Una **moda** es un máximo relativo del histograma. La presencia de varias modas suele indicar que hay varias poblaciones mezcladas. - Si hay datos atípicos o no en la muestra. - Si la distribución es o no simétrica y qué tipo de simetría tiene. - Cuál es la magnitud en torno a la cual se distribuyen los datos. - Si hay mucha dispersión en torno a este centro. --- # Histogramas La forma de una histograma depende del número de clases: <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/nclases-1.png" style="display: block; margin: auto;" /> --- # Estimadores de núcleo Supongamos que `\(x_1,\ldots,x_n\)` son realizaciones independientes de una v.a. `\(X\)` con densidad `\(f\)`. El objetivo es calcular a partir de los datos una función `\(\hat{f}\)` tal que `\(\hat{f}\approx f\)` (un estimador de `\(f\)`). Un histograma es una aproximación de `\(f\)` pero las funciones de densidad no son constantes a trozos. ¿Podemos hacerlo mejor? Las áreas bajo la función de densidad corresponden a las probabilidades de que la variable tome valores en el correspondiente intervalo. Si `\(h\approx 0\)` se tiene: $$ \mbox{P}(x-h\leq X\leq x+h) = \int_{x-h}^{x+h} f(t) dt \approx 2h f(x). $$ Es decir, $$ f(x)\approx \frac{1}{2h} \mbox{P}(x-h\leq X\leq x+h). $$ --- # Estimadores de núcleo Sustituyendo la probabilidad por una proporción: `$$\hat{f}(x) = \frac{1}{2h} \frac{\#\{i:\, |x-x_i| < h\}}{n}.$$` Si `\(K(x) = (1/2) \mathbb{I}_{\{|x|\leq 1\}}\)`, donde `\(\mathbb{I}_A\)` denota la función indicatriz sobre `\(A\)`, `$$\hat{f}(x) = \frac{1}{nh} \sum_{i=1}^nK\left(\frac{x-x_i}{h}\right).$$` Sustituyendo la indicatriz por otra función más suave obtendremos aproximaciones que se asemejan más a una densidad típica: $$ \hat{f}(x) = \frac{1}{nh} \sum_{i=1}^nK\left(\frac{x-x_i}{h}\right) $$ - `\(h>0\)` es el **parámetro de suavizado** - `\(K\)` es el núcleo (que suele verificar `\(K\geq 0\)` y `\(\int K = 1\)`, es decir, que es una función de densidad). --- # Estimadores de núcleo Habitualmente, núcleo gaussiano: `\(K(x)=(2\pi)^{-1/2} e^{-x^2/2}\)`. <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-13-1.png" style="display: block; margin: auto;" /> --- # Estimadores de núcleo Usamos el comando `geom_density` para representar los estimadores de núcleo. - Usa por defecto un núcleo gaussiano. - La correcta selección del parámetro de suavizado `\(h\)` es muy importante. ```r ggplot(mercurio) + geom_density(aes(x = CONC), linewidth = 1.2, fill = 'olivedrab4') + labs(x = "Concentración (ppm)", y = " ") ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-14-1.png" style="display: block; margin: auto;" /> --- # Estimadores de núcleo Si usamos un núcleo rectangular el resultado es menos suave, más parecido al histograma: ```r ggplot(mercurio) + geom_density(aes(x = CONC), linewidth = 1.2, fill = 'olivedrab4', kernel = 'rectangular') + labs(x = "Concentración (ppm)", y = " ") ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-15-1.png" style="display: block; margin: auto;" /> --- # Estimadores de núcleo <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-16-1.png" style="display: block; margin: auto;" /> --- # Estimadores de núcleo Una ventaja de los estimadores del núcleo respecto a los histogramas es que son más apropiados para superponerlos en el mismo gráfico. ```r ggplot(mercurio) + geom_density(aes(x = CONC, linetype = RIO), linewidth = 1.2) + labs(x = "Concentración (ppm)", y = " ") ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-17-1.png" style="display: block; margin: auto;" /> --- # Descripción numérica: media La **media** de un conjunto de números `\(x_1,\ldots,x_n\)` se define como `$$\bar{x}=\frac{x_1+\cdots + x_n}{n}.$$` La suma de las desviaciones de los datos a la media es igual a cero: $$ (x_1-\bar{x}) + \cdots + (x_n -\bar{x}) = 0. $$ Es el centro de gravedad de un histograma (el lugar en el que apoyarlo para que se mantuviera en equilibrio) Es difícil de interpretar si - La distribución es muy asimétrica - Hay datos atípicos --- # Descripción numérica: mediana Una medida de posición que se comporta mejor respecto de estos problemas es la **mediana**. Para calcularla: - Se ordenan los datos de menor a mayor. - Si el número de datos es impar, la mediana es el dato que ocupa la posición central. Si es par, la mediana es la media de los dos datos centrales. La mediana es más **robusta** que la media pero hace un uso **menos eficiente** de la información contenida en los datos. <img src="data:image/png;base64,#figs/salarios2018.jpg" width="65%" style="display: block; margin: auto;" /> (Fuente: [INE](https://t.co/kXf6yRLXyD?amp=1)) --- # Descripción numérica: percentiles Los percentiles dan una idea más completa de la posición de un conjunto de datos. El cuantil `\(p\)` (o percentil `\(100p\)`) es el valor que deja a su **izquierda** una proporción `\(p\)` (o porcentaje `\(100p\)`) de los datos. Al percentil 25 se le llama también **primer cuartil** (notación `\(Q_1\)`). El segundo cuartil es la mediana (que también es el percentil 50). Al percentil 75 se le suele llamar **tercer cuartil** (notación `\(Q_3\)`) --- # Descripción numérica: percentiles <img src="data:image/png;base64,#figs/ine_salarios_2017.jpg" width="110%" style="display: block; margin: auto;" /> (Fuente: [INE, encuesta de estructura salarial](https://www.ine.es/jaxiT3/Tabla.htm?t=28191)) Calcula el valor de `\(x\)` y de `\(y\)` en cada frase: - La media y la mediana de los salarios en España en 2017 fueron `\(x\)` e `\(y\)`, respectivamente. - Un 10% de las mujeres ganaba más de `\(x\)` euros. - Un `\(x\%\)` de lo hombres ganaba entre `\(11393.30\)` y `\(45813.46\)` euros. - Un 90% de las mujeres ganaba más de `\(x\)` euros. --- # Cálculo de medidas numéricas con `dplyr` Se usa el verbo `summarise` de `dplyr`, usualmente en combinación con `group_by`: .center[![gb](data:image/png;base64,#figs/groupby.jpg)] --- # Cálculo de medidas numéricas con `dplyr` Calcula la media y la mediana de la concentración de mercurio: ```r mercurio %>% group_by(RIO) %>% summarise(media_hg = mean(CONC), Q1_hg = quantile(CONC, probs = 0.25), mediana_hg = median(CONC), Q3_hg = quantile(CONC, probs = 0.75)) ``` ``` ## # A tibble: 2 x 5 ## RIO media_hg Q1_hg mediana_hg Q3_hg ## <fct> <dbl> <dbl> <dbl> <dbl> ## 1 Lumber 1.08 0.57 0.91 1.5 ## 2 Wacamaw 1.28 0.618 1.05 1.8 ``` --- # Cálculo de medidas numéricas con `dplyr` El comando `across()` sirve para trabajar con muchas variables/medidas a la vez ```r mercurio %>% group_by(RIO) %>% summarise(across(LONG:CONC, list(media = mean, mediana = median))) ``` ``` ## # A tibble: 2 x 7 ## RIO LONG_media LONG_mediana PESO_media PESO_mediana CONC_media CONC_medi~1 ## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 Lumber 39.4 39 1197. 898 1.08 0.91 ## 2 Wacamaw 40.4 39.2 1111. 860 1.28 1.05 ## # ... with abbreviated variable name 1: CONC_mediana ``` --- # Descripción numérica: varianza La **varianza** es *básicamente* el promedio de las desviaciones al cuadrado de los datos a su media. ```r x <- c(1, 2, 3, 5, 7, 0) x - mean(x) # desviaciones ``` ``` ## [1] -2 -1 0 2 4 -3 ``` ```r (x - mean(x))^2 # desviaciones al cuadrado ``` ``` ## [1] 4 1 0 4 16 9 ``` ```r mean((x-mean(x))^2) # varianza ``` ``` ## [1] 5.666667 ``` --- # Descripción numérica: varianza Solo disponemos de `\(n-1\)` desviaciones independientes por lo que es más correcto dividir por `\(n-1\)` que por `\(n\)`. `$$S^2 = \frac{(x_1-\bar{x})^2+\cdots +(x_n-\bar{x})^2}{n-1}.$$` **Ejercicio:** pensar en el caso `\(n=2\)`. Fórmula alternativa para `\(S^2\)` `$$S^2 = \frac{n}{n-1}\left(\frac{x_1^2+\cdots + x_n^2}{n} - \bar{x}^2 \right).$$` --- # Descripción numérica: desviación típica y RI La **desviación típica** es la raíz cuadrada de `\(S^2\)`. Mide la dispersión en la misma escala que los datos: `$$S =\sqrt{\frac{(x_1-\bar{x})^2+\cdots +(x_n-\bar{x})^2}{n-1}}.$$` Una medida adimensional relacionada es el **coeficiente de variación**: `$$\mbox{CV} = \frac{S}{|\bar{X}|}.$$` Una medida de dispersión más robusta que las anteriores es el **rango intercuartílico**, `\(\mbox{RI}= Q_3-Q_1\)`. --- # Ejemplo Para los datos de los peces, las medidas de dispersión se pueden calcular de la forma siguiente: ```r mercurio %>% group_by(RIO) %>% summarise(var_hg = var(CONC), sd_hg = sd(CONC), RI_hg = quantile(CONC, probs = 0.75) - quantile(CONC, probs = 0.25)) ``` ``` ## # A tibble: 2 x 4 ## RIO var_hg sd_hg RI_hg ## <fct> <dbl> <dbl> <dbl> ## 1 Lumber 0.421 0.649 0.93 ## 2 Wacamaw 0.687 0.829 1.18 ``` La regla por defecto de **selección del parámetro de suavizado en los estimadores de núcleo** que usa `geom_density` es `$$h_n = 0.9 \hat{\sigma} n^{-1/5},$$` donde `\(\hat{\sigma}=\min\{s,\mbox{RI}/1.34\}\)` (regla de Silverman). --- # Cuestiones .small[ 1. ¿Por qué crees que se divide el rango intercuartílico por 1.34 en la regla de Silverman para seleccionar el parámetro de suavizado? (Indicación: haz la hipótesis de que los datos tienen distribución normal). 2. Da un ejemplo de un conjunto de datos tal que `\(S^2=0\)`. 4. Dado un conjunto de observaciones medidas en kg, supongamos que cambiamos las unidades y las pasamos a gramos (es decir, multiplicamos por mil). Determina si son verdaderas o falsas las siguientes afirmaciones: * Tanto la media como la mediana de los nuevos datos se multiplican también por mil. * La varianza se multiplica también por mil. 5. ¿Cómo cambiaría en cada caso la desviación típica? 6. Ahora sumamos 100 a todos los datos. Determina si son verdaderas o falsas las siguientes afirmaciones: * Los cuartiles no cambian. * El rango intercuartílico no cambia. * La desviación típica no cambia. ] --- # Diagramas de cajas Las partes de un diagrama de cajas son las siguientes: <img src="data:image/png;base64,#figs/boxplot.jpg" width="90%" style="display: block; margin: auto;" /> (Fuente de este gráfico: [towardsdatascience](https://towardsdatascience.com/understanding-boxplots-5e2df7bcbd51)) Es útil sobre todo para comparar las características de varios conjuntos de datos. --- # Diagramas de cajas Se usa `geom_boxplot`. ```r ggplot(mercurio) + geom_boxplot(aes(x = ESTACION, y = CONC, fill = RIO)) + labs(title = "Concentración Hg en los peces capturados en varias estaciones", x = "Estación", y = "Concentración Hg (ppm)") ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-25-1.png" style="display: block; margin: auto;" /> --- # Descripción numérica: asimetría Para cuantificar la asimetría de la distribución se usa el siguiente coeficiente: `$$\hat{\gamma} = \frac{\sum_{i=1}^n (x_i - \bar{x})^3}{(n-1)s^3}=\frac{1}{n-1}\sum_{i=1}^n\left(\frac{x_i-\bar{x}}{s}\right)^3.$$` Si la distribución es simétrica se espera que `\(\gamma \approx 0\)`. Si la distribución es asimétrica a la derecha predominarán los sumandos muy positivos por lo que `\(\gamma>0\)`. Lo contrario sucede si la distribución es asimétrica a la izquierda. Se divide por `\(S^3\)` para que la medida no dependa de las unidades. --- # Descripción numérica: asimetría Para los datos de mercurio en los peces obtenemos los valores siguientes (previamente escribimos una función muy sencilla que lo calcula): ```r asimetria <- function(x){ n <- length(x) sum((x - mean(x))^3) / ((n-1) * sd(x)^3) } mercurio %>% group_by(RIO) %>% summarise(asim_hg = asimetria(CONC), asim_longitud = asimetria(LONG), asim_peso = asimetria(PESO)) ``` ``` ## # A tibble: 2 x 4 ## RIO asim_hg asim_longitud asim_peso ## <fct> <dbl> <dbl> <dbl> ## 1 Lumber 1.48 0.559 1.50 ## 2 Wacamaw 0.939 0.363 1.21 ``` --- # Descripción numérica: curtosis La curtosis está relacionada con el momento de cuarto orden: `$$\hat{\kappa} = \frac{\sum_{i=1}^n (x_i - \bar{x})^4}{(n-1)s^4}=\frac{1}{n-1}\sum_{i=1}^n\left(\frac{x_i-\bar{x}}{s}\right)^4.$$` <br> - En general resulta difícil de interpretar - Si `\(Z\equiv \mbox{N}(0,1)\)`, entonces `\(\mbox{E}(Z^4)=3\)` - Cuando hay simetría, a veces se compara `\(\hat{\kappa}\)` con 3 para determinar si la distribución de los datos se puede considerar o no normal --- # Diagramas de dispersión Para visualizar la relación entre dos variables `\(x\)` e `\(y\)` se representa cada observación `\((x_i, y_i)\)` como un punto. La nube de puntos se suele llamar **diagrama de dispersión** (*scatter plot*). ```r ggplot(mercurio) + geom_point(aes(x = LONG, y = CONC)) ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-27-1.png" style="display: block; margin: auto;" /> --- # Diagramas de dispersión Para diferenciar las observaciones de cada uno de los ríos podemos usar distinto color: ```r ggplot(mercurio) + geom_point(aes(x = LONG, y = CONC, col = RIO)) ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-28-1.png" style="display: block; margin: auto;" /> --- # Diagramas de dispersión Algunos aspectos en los que fijarse: - ¿Hay asociación entre `\(X\)` e `\(Y\)`? - ¿Es la asociación lineal o no lineal? ¿Es positiva o negativa? - ¿Hay heterocedasticidad? (Esto significa que la varianza de la `\(Y\)` depende de los valores de `\(X\)`. Lo contrario de la heterocedasticidad es la homocedasticidad.) - ¿Hay datos atípicos? En un diagrama de dispersión hay atípicos que no lo son respecto a ninguna de las variables. - ¿Hay grupos homogéneos (clústeres)? --- # Matriz de diagramas de dispersión Para un número moderado de variables puede ser útil dibujar los diagramas de dispersión para cada par de variables: ```r # requiere la extensión GGally ggpairs(mercurio, columns = 3:5, aes(col = RIO)) ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-29-1.png" style="display: block; margin: auto;" /> --- # Covarianza Se usa para cuantificar el grado de relación **lineal** entre dos variables `\(x\)` e `\(y\)`. A partir de una muestra de datos bidimensionales `\((x_1,y_1),\ldots,(x_n,y_n)\)` la **covarianza** entre `\(x\)` e `\(y\)` es `$$S_{x,y} = \frac{1}{n-1} \sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})$$` Una fórmula alternativa de la covarianza es: `$$S_{x,y} = \frac{n}{n-1} \left[\frac{\sum_{i=1}^n x_iy_i}{n} - \bar{x}\bar{y} \right]$$` --- # Signo de la covarianza <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-30-1.png" style="display: block; margin: auto;" /> --- # Propiedades de la covarianza - La covarianza es simétrica: `\(S_{x,y} = S_{y,x}\)` - La covarianza de una variable consigo misma es la varianza de esa variable: `\(S_{x,x} = S^2_x\)` - Si definimos los vectores `\(x = (x_1-\bar{x},\ldots,x_n-\bar{x})\)` e `\(y = (y_1-\bar{y},\ldots,y_n-\bar{y})\)`, entonces `\((n-1)S_{x,y} = \langle x, y\rangle\)`, donde `\(\langle \cdot, \cdot\rangle\)` denota el producto escalar euclídeo. - De igual forma, `\((n-1)S^2_x = \|x\|^2\)`, donde `\(\|\cdot\|\)` es la norma euclídea. - La covarianza depende de las unidades en las que medimos las variables. De hecho, si `\(\lambda\in\mathbb{R}\)`, entonces `\(S_{\lambda x,y}=?\)` --- # Correlación Para eliminar esta dependencia se *normaliza* la covarianza dividiendo por el producto de desviaciones típicas. Esto da lugar al **coeficiente de correlación** entre `\(x\)` e `\(y\)`: `$$r_{x,y} = \frac{S_{x,y}}{S_xS_y}$$` - Si `\(\lambda\in\mathbb{R}\)`, entonces `\(r_{\lambda x,y}=?\)` - ¿Cuánto vale `\(r_{x,x}\)`? - Si `\(x = (x_1-\bar{x},\ldots,x_n-\bar{x})\)` e `\(y = (y_1-\bar{y},\ldots,y_n-\bar{y})\)`, entonces `$$r_{x,y} = \frac{\langle x, y\rangle}{\|x\|\cdot \|y\|}$$` - Por la desigualdad de Cauchy-Schwarz, `\(|r_{x,y}|\leq 1\)`. Además, `\(|r_{x,y}|=1\)` si y solo si existen `\(a,b\in\mathbb{R}\)` tales que `\(y_i = ax_i + b\)`, para todo `\(i\)`. --- # Cálculo de la covarianza y la correlación ```r mercurio %>% group_by(RIO) %>% summarise(cor_hg_long = cor(CONC, LONG), cor_hg_peso = cor(CONC, PESO), cor_long_peso = cor(LONG, PESO)) ``` ``` ## # A tibble: 2 x 4 ## RIO cor_hg_long cor_hg_peso cor_long_peso ## <fct> <dbl> <dbl> <dbl> ## 1 Lumber 0.554 0.381 0.872 ## 2 Wacamaw 0.707 0.704 0.939 ``` --- # Regresión lineal Observamos dos variables, `\(x\)` e `\(y\)`. El objetivo es predecir o aproximar el valor de la variable `\(y\)` a partir del valor de la variable `\(x\)`. - A la variable `\(y\)` se llama **variable respuesta**. - A La variable `\(x\)` se llama **variable regresora o explicativa**. En un problema de regresión (a diferencia de cuando calculamos el coeficiente de correlación) el papel de las dos variables no es simétrico. --- # Recta de mínimos cuadrados Una recta `\(y = \beta_0+\beta_1 x\)` es adecuada para representar la relación entre `\(y\)` y `\(x\)` si las diferencias `\(y_i - (\beta_0+\beta_1 x_i)\)` tienden a ser pequeñas. La **recta de regresión de mínimos cuadrados** viene dada por los valores `\(\hat{\beta}_0\)` y `\(\hat{\beta}_1\)` para los que se minimiza la función `$$\phi(\beta_0,\beta_1) = \sum_{i=1}^n [y_i - (\beta_0+\beta_1 x_i)]^2$$` Para encontrar la recta basta determinar los valores para los que se anula el gradiente de `\(\phi\)`: `\begin{align*} \sum_{i=1}^n [y_i - (\hat{\beta}_0+\hat{\beta}_1 x_i)] &= 0 \\ \sum_{i=1}^n x_i[y_i - (\hat{\beta}_0+\hat{\beta}_1 x_i)] &= 0. \end{align*}` --- # Recta de mínimos cuadrados De la primera ecuación se deduce fácilmente que la recta pasa por `\((\bar{x},\bar{y})\)`. De la segunda ecuación deducimos el valor de la pendiente. La recta de mínimos cuadrados es (ejercicio): `$$y - \bar{y} = \hat{\beta}_1 (x -\bar{x}),$$` donde la pendiente viene dada por `$$\hat{\beta}_1 = \frac{S_{x,y}}{S^2_x} = r_{x,y} \frac{S_y}{S_x} = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i-\bar{x})^2}$$` --- # Residuos y valores ajustados El valor pronosticado o previsto (también llamado **valor ajustado**) de `\(y_i\)` dado `\(x = x_i\)` es `\(\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i\)`. Las diferencias entre los valores reales de la variable respuesta y los valores ajustados, `\(e_i =y_i - \hat{y}_i\)`, se denominan **residuos**. La recta de mínimos cuadrados es aquella para la que el valor de `\(\sum_{i=1}^n e^2_i\)` es el mínimo posible. Además, las dos ecuaciones que hemos tenido que resolver para calcularla equivalen a `\(\sum_{i=1}^n e_i = 0\)` y `\(\sum_{i=1}^n x_i e_i = 0\)`. Los residuos tienen media cero y están incorrelados con la variable regresora. --- # Varianza residual Para medir el error cometido en la aproximación dada por la recta de regresión se puede usar la **varianza residual** `$$S^2_R = \frac{\sum_{i=1}^n e_i^2}{n-2}$$` Puede demostrarse (ejercicio) que `$$\frac{\sum_{i=1}^n e_i^2}{n-1} = S^2_y (1-r^2_{x,y})$$` Equivalentemente, `$$\frac{n-2}{n-1} S^2_R = S_y^2 (1-r_{x,y}^2)$$` Para `\(n\)` moderado o grande, `\(r_{x,y}^2 \approx 1 - S^2_R/S_y^2\)`. Al valor del coeficiente de correlación al cuadrado también se le llama **coeficiente de determinación**. Es la fracción de la varianza de `\(y\)` que explica la variable `\(x\)`. --- # Cálculo y representación de la recta El comando básico para calcular la recta de mínimos cuadrados es `lm`. En combinación con `geom_smooth` permite también representarla: ```r lm(CONC ~ LONG, data = mercurio) ``` ``` ## ## Call: ## lm(formula = CONC ~ LONG, data = mercurio) ## ## Coefficients: ## (Intercept) LONG ## -1.13165 0.05813 ``` ```r lm(CONC ~ LONG, data = mercurio_lumber) ``` ``` ## ## Call: ## lm(formula = CONC ~ LONG, data = mercurio_lumber) ## ## Coefficients: ## (Intercept) LONG ## -0.62387 0.04318 ``` --- # Cálculo y representación de la recta ```r ggplot(mercurio, aes(x = LONG, y = CONC, col = RIO)) + geom_point() + geom_smooth(method = 'lm', se = FALSE) + scale_colour_brewer(palette="Set1") ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-33-1.png" style="display: block; margin: auto;" /> --- # Cálculo y representación de la recta ```r ggplot(mercurio, aes(x = LONG, y = CONC)) + geom_point() + geom_smooth(method = 'lm', se = FALSE) + facet_wrap(~RIO, ncol = 1) ``` <img src="data:image/png;base64,#est1-descripcion-23_files/figure-html/unnamed-chunk-34-1.png" style="display: block; margin: auto;" />