Breve introducción

El siguiente material pretende ser una guía más práctica que teórica, escrita en un lenguaje amigable. En ella, aplicaremos las herramientas básicas de R y sus paquetes para realizar análisis estadísticos, en su mayoría con bases de datos reales. La primera sección aborda los conceptos esenciales de la estadística descriptiva y muestra cómo implementar sus herramientas en el software R. Posteriormente, utilizamos conocimientos de inferencia estadística para introducir la regresión lineal y aplicarla en la creación de un modelo predictivo.

1 Estadística descriptiva

La estadística descriptiva es una útil herramienta para comprender y resumir conjuntos de datos. Nos permite explorar la estructura, las tendencias y las características de los datos de manera sistemática y cuantitativa. En lugar de simplemente mirar cifras, la estadística descriptiva nos brinda un marco analítico para interpretar y extraer información significativa de los datos.

En esta guía, exploraremos la estadística descriptiva utilizando RStudio. Para hacerlo, descargaremos bases de datos que contienen información sobre enfermedades cardíacas en humanos y zarigüeyas. A través de ejemplos prácticos, aprenderemos a aplicar técnicas estadísticas descriptivas utilizando el lenguaje de programación R y obtener conclusiones de estas.

Al utilizar estas bases de datos, por un lado, podremos analizar cómo se relacionan diferentes variables, como la edad, el sexo, los síntomas y los resultados de las pruebas médicas respecto a enfermedades cardiacas. Para el caso de las zarigüeyas, además de su edad y sexo, analizaremos las medidas de partes de su cuerpo, como el tamaño de sus patas, orejas, pecho y más.

1.1 Importación de bases de datos

El siguiente enlace nos llevará a la plataforma Kaggle, en donde está alojada base de datos sobre enfermedades cardiacas. Una vez abierto, solo debemos pinchar en “Download” y guardar el documento en donde se prefiera. Posteriormente, ejecutaremos RStudio y procederemos a importar la base de datos.

La extensión de la base de datos descargada es .csv, lo cual significa que está escrita en un documento de texto, y cada está dato separado por comas. RStudio nos hace bastante simple el trabajo, solo debemos seleccionar las siguientes opciones:

\[ \textbf{File} \rightarrow \textbf{Import Dataset} \rightarrow \textbf{From Text (base)} \rightarrow \textbf{Seleccionar base de datos} \rightarrow \textbf{Import} \]

Importación de la base de datos desde \(\textsf{RStudio}\)


Ya tenemos lista para usar nuestra base de datos sobre enfermedades cardiacas. Para la base de datos de zarigüeyas, debemos abrir el siguiente enlace y repetir los pasos anteriores.

La siguiente tabla es un resumen de nuestra base de datos sobre enfermedades cardiacas.


Mientras que la siguiente tabla es un resumen de la base de datos sobre zarigüeyas

1.2 Medidas descriptivas

Las medidas descriptivas nos proporcionan información sobre diferentes aspectos de los datos, como su tendencia central, dispersión, posición y correlación. Entre las medidas descriptivas más comunes se incluyen la media, la mediana, la moda, el rango, la varianza, desviación estándar, el rango intercuartílico y la correlación de Pearson. Si no te suena alguna de las medidas mencionadas, no te preocupes porque en las siguientes subsecciones daremos más detalles de cada una.

Gracias a R y sus paquetes de base, podemos calcular fácilmente medidas descriptivas utilizando funciones integradas. Estas herramientas nos permitirán analizar y visualizar rápidamente la distribución y las características de los datos.

Antes de comenzar a calcular estadísticas, calentemos los motores en la consola de RStudio. Utilizando la función head(), e introduciéndole como argumento el nombre de nuestra base de datos de enfermedades cardiacas, heart, podremos ver las primeras filas de datos.

head(heart)


También podemos chequear el número de observaciones que presenta nuestra base de datos mediante la función length() poniéndole como argumento cualquier variable.

length(heart$Sexo)
## [1] 918

Hagamos exactamente lo mismo para la base de datos de las zarigüeyas

head(possum)
length(heart$sex)
## [1] 0

1.2.1 Medidas de tendencia central

Podemos entender las medidas de tendencia central como estadísticas que representan la ubicación central o típica de un conjunto de datos. Las tres medidas de tendencia central más comunes son la media, la mediana y la moda.

  • Media: La media es el promedio aritmético de todos los valores en un conjunto de datos. Se calcula sumando todos los valores y dividiendo por el número total de valores.

  • Mediana: La mediana es el valor que divide al conjunto de datos en dos partes iguales, es decir, el valor medio cuando los datos están ordenados de manera ascendente o descendente.

  • Moda: La moda es el valor que ocurre con mayor frecuencia en un conjunto de datos.

Estas medidas son útiles para resumir y comprender la distribución de los datos. En R, puedes calcular estas medidas fácilmente utilizando funciones como mean() y median(), para la media y la mediana respectivamente.

Podemos, por ejemplo, obtener la edad media de nuestra muestra:

mean(heart$Edad) # Escogemos la variable Edad de la base de datos heart con "$"
## [1] 53.51089

Usando el paquete dplyr podemos obtener el mismo resultado

library(dplyr) # Cargamos el paquete

heart %>% # Escogemos la base de datos y luego usamos el pipe %>% 
  pull(Edad) %>% # Obtenemos la variable Edad como vector
  mean() # Calculamos el promedio de la edad
## [1] 53.51089

O se puede ser más específicos en nuestro cálculo, y obtener la edad media de las personas que tienen alguna enfermedad cardiaca y las que no.

Primero calculemos la edad media de los que tienen enfermedad cardiaca usando las funciones base de R. Usando las herramientas lógicas, definiremos las variables de las personas que tienen una enfermedad cardiaca

tiene_enfermedad_cardiaca <- heart$EnfCard == 1 # Guardamos como variable todos los que tengan como observación de enfermedad cardiaca = 1
mean(heart[tiene_enfermedad_cardiaca, 1])
## [1] 55.89961

Por otro lado, usando el paquete dplyr podemos calcular el promedio de la edad los que no tienen enfermedad cardiaca.

heart %>% 
  filter(EnfCard == 0) %>% # Filtramos las filas que tengan la variable EnfCard = 0
  pull(Edad) %>% 
  mean()
## [1] 50.55122

Podemos decir entonces que en promedio, los pacientes de nuestra muestra que presentan alguna enfermedad cardiaca son cinco años mayor a las que no.

Otro ejercicio sería calcular la mediana de alguna variable que nos interese, por ejemplo, del colesterol.

median(heart$Col)
## [1] 223

O usando dplyr

heart %>% 
  pull(Col) %>% 
  median()
## [1] 223

Nos dio como resultado una mediana de 223 mm/dl para la variable del colesterol. Ahora, notemos que podemos usar la función que sirve crear histogramas,hist(), para crear una tabla de frecuencia agrupada en intervalos. Para ello, dentro de los argumentos opcionales que nos brinda la función, debemos desactivar plot con plot = FALSE. Una vez creada nuestra tabla de frecuencia, obtendremos intervalos de la forma [a,b) para nuestros datos. Para efectos prácticos, le pediremos a R que solo nos lea las primeras dos líneas del output de la función hist que usaremos, finalizando con [1:2].

hist(heart$Col, plot = FALSE)[1:2] # Graficamos una tabla de frecuencia, y solo pedimos las primeras dos líneas del output del comando.
## $breaks
##  [1]   0  50 100 150 200 250 300 350 400 450 500 550 600 650
## 
## $counts
##  [1] 172   3  17 130 291 202  76  14   5   4   2   1   1

Tenemos un llamativo intervalo [0,50) con 172 datos. Significa que hay 172 pacientes con niveles de colesterol sérico entre 0 y 50 mm/dl, números demasiado bajos para niveles de colesterol sérico en seres humanos. ¿Qué está ocurriendo? Hagamos una filtración de nuestra variable colesterol, pidiéndole a R que nos muestre específicamente los valores de cada paciente que tenga los niveles de colesterol sérico por debajo de 50 mm/dl.

heart %>% 
  filter(Col < 50) %>%  
# Filtramos los datos para tener una columna de valores de colesterol bajo a 50. 
  
# Solicitamos a R los valores específicos de cada observación que haya -
# cumplido el filtro anterior.  
  pull(Col) 
##   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Por las observaciones anteriores, podemos entender que nuestra base de datos nos está diciendo que a 172 pacientes no se les midió el colesterol total en su sangre al momento de examinarlos (o sus niveles de colesterol en la sangre realmente eran 0 mm/dl, pero omitiremos esta posibilidad ya que es imposible que ocurra). Por tanto, nuestros cálculos anteriores para obtener la mediana están incluyendo una cantidad considerable de valores iguales a 0, los cuales más nos vale ignorar o podríamos realizar un análisis sesgado.

Teniendo en cuenta lo anterior, presentamos los valores de la mediana solo considerando los valores de colesterol mayores que 0.

# Calculamos la mediana de la variable colesterol, aunque necesitamos-
# introducirle como argumento a la función median() un vector. -
# Para ello, filtramos de la variable colesterol el vector que contiene -
# los valores de colesterol mayores que 0.
median(heart$Col[heart$Col > 0]) 
## [1] 237

O usando dplyr

heart %>% 
  filter(Col > 0) %>% # Primero filtramos los valores de colesterol mayores a 0.
  pull(Col) %>% # De aquel filtro, solicitamos a R el vector con las observaciones
  median() # Introducimos a la función median() el vector anterior.
## [1] 237

La mediana varió 14 mg/dL respecto a la anterior con la muestra sesgada. ¿Cuánto variará el promedio, en un análisis no sesgado?

Calculando el promedio del colesterol sérico incluyendo las observaciones iguales a 0, es decir, en un cálculo sesgado, obtenemos

mean(heart$Col)
## [1] 198.7996

Hacemos lo mismo usando dplyr para los datos mayores que 0

heart %>% 
  filter(Col > 0) %>% 
  pull(Col) %>% 
  mean()
## [1] 244.6354

Varió 46 mg/dL, lo que es un gran diferencia. Valió la pena ser cuidadosos.

Para calcular la moda no tenemos una función que venga de base en \(\textsf{R}\), pero en internet podemos encontrar un montón de soluciones a este problema. Para no descargar un paquete, usaremos una función creada por un usuario del foro stackoverflow. En este enlace encontraremos la función creada y más opciones programadas por otros usuarios.

Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

Usémosla para comprobar cuál es el tipo de dolor de pecho más frecuente

Mode(heart$TDP)
## [1] "ASY"

La observación que más se repite para la variable tipo de dolor de pecho es que la mayoría de los pacientes no presentan dolor de pecho. Podremos ver a detalle los valores de las otras observaciones en la sección Tablas de frecuencia.

1.2.2 Medidas de posición

Los cuartiles, deciles y percentiles son medidas estadísticas utilizadas para dividir un conjunto de datos ordenados en partes iguales o proporcionales. Estas medidas son importantes en la estadística descriptiva porque proporcionan información sobre la distribución de todos los datos.

  • Cuartiles: Los cuartiles dividen un conjunto de datos ordenados en cuatro partes iguales, es decir, en tres puntos que dividen los datos en cuatro grupos de igual tamaño. Los cuartiles más comúnmente utilizados son el primer cuartil (Q1), que representa el 25% de los datos, el segundo cuartil (Q2) que es igual a la mediana y divide los datos en dos partes iguales, y el tercer cuartil (Q3) que representa el 75% de los datos.

  • Deciles: Los deciles dividen un conjunto de datos en diez partes iguales, es decir, en nueve puntos que dividen los datos en diez grupos de igual tamaño. El primer decil (D1) representa el 10% de los datos, el segundo decil (D2) representa el 20% de los datos, y así sucesivamente.

  • Percentiles: Los percentiles dividen un conjunto de datos en cien partes iguales, es decir, en 99 puntos que dividen los datos en cien grupos de igual tamaño. El percentil P25 representa el 25% de los datos, el percentil P50 es igual a la mediana y representa el 50% de los datos, y el percentil P75 representa el 75% de los datos.

En R, podemos calcular cuartiles, deciles y percentiles utilizando las funciones base quantile() y summary()

Creemos un vector de datos de edad y con él calculemos cuartiles, deciles y percentiles.

datos_edad <- c(10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60)

# Calcular cuartiles
cuartiles <- quantile(datos_edad, probs = c(0.25, 0.5, 0.75)) 
# Con la opción props le pedimos que nos calcule el 25%, 50% y 75%

print(cuartiles)
##  25%  50%  75% 
## 22.5 35.0 47.5
# Calcular deciles
deciles <- quantile(datos_edad, probs = seq(0.1, 0.9, by = 0.1)) 
# Esta vez le pedimos que nos calcule desde el 10% al 90% de los datos, -
# cada uno separado por un 10%

print(deciles)
## 10% 20% 30% 40% 50% 60% 70% 80% 90% 
##  15  20  25  30  35  40  45  50  55
# Calcular percentiles
percentiles <- quantile(datos_edad, probs = seq(0.01, 0.99, by = 0.01)) 
print(percentiles)
##   1%   2%   3%   4%   5%   6%   7%   8%   9%  10%  11%  12%  13%  14%  15%  16% 
## 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0 15.5 16.0 16.5 17.0 17.5 18.0 
##  17%  18%  19%  20%  21%  22%  23%  24%  25%  26%  27%  28%  29%  30%  31%  32% 
## 18.5 19.0 19.5 20.0 20.5 21.0 21.5 22.0 22.5 23.0 23.5 24.0 24.5 25.0 25.5 26.0 
##  33%  34%  35%  36%  37%  38%  39%  40%  41%  42%  43%  44%  45%  46%  47%  48% 
## 26.5 27.0 27.5 28.0 28.5 29.0 29.5 30.0 30.5 31.0 31.5 32.0 32.5 33.0 33.5 34.0 
##  49%  50%  51%  52%  53%  54%  55%  56%  57%  58%  59%  60%  61%  62%  63%  64% 
## 34.5 35.0 35.5 36.0 36.5 37.0 37.5 38.0 38.5 39.0 39.5 40.0 40.5 41.0 41.5 42.0 
##  65%  66%  67%  68%  69%  70%  71%  72%  73%  74%  75%  76%  77%  78%  79%  80% 
## 42.5 43.0 43.5 44.0 44.5 45.0 45.5 46.0 46.5 47.0 47.5 48.0 48.5 49.0 49.5 50.0 
##  81%  82%  83%  84%  85%  86%  87%  88%  89%  90%  91%  92%  93%  94%  95%  96% 
## 50.5 51.0 51.5 52.0 52.5 53.0 53.5 54.0 54.5 55.0 55.5 56.0 56.5 57.0 57.5 58.0 
##  97%  98%  99% 
## 58.5 59.0 59.5

También podemos obtener un resumen rápido usando summary()

summary(datos_edad)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    10.0    22.5    35.0    35.0    47.5    60.0

Podemos sacar algunas conclusiones con el resumen:

  • La edad mínima observada es de 10 años, y la máxima de 60 años.

  • La edad promedio entre los datos es de 35 años, que coincide con la mediana.

  • El 25% de las personas tiene 22.5 años o menos, la mitad tiene 35 años o menos, y el 75% tiene 47.5 años o menos.

1.2.3 Medidas de dispersión

En el análisis estadístico, las medidas de dispersión o variabilidad son fundamentales para comprender qué tanto varían o cuanto se diferencian respecto de su media. Mientras que las medidas de tendencia central, como la media, la mediana y la moda, proporcionan información sobre el valor típico o central de un conjunto de datos, las medidas de dispersión revelan la amplitud y la variabilidad de esos datos.

  • Rango: El rango es la diferencia entre el valor máximo y el valor mínimo en un conjunto de datos. Es la medida más simple de dispersión y proporciona una idea de la amplitud total de los datos, representando la máxima variabilidad presente en los datos.

  • Rango intercuartil (IQR): El rango intercuartil es la diferencia entre el tercer cuartil (Q3) y el primer cuartil (Q1). Representa la dispersión de los datos dentro del 50% central de la distribución y es una medida robusta que no se ve afectada por valores atípicos.

  • Varianza: La varianza es una medida de dispersión que describe cuánto varían los datos respecto a la media. Se calcula como la media de los cuadrados de las diferencias entre cada valor y la media del conjunto de datos.

  • Desviación estándar: La desviación estándar es la raíz cuadrada de la varianza y es una medida de dispersión más interpretable que la varianza. Indica cuánto se alejan, en promedio, los datos de la media y se expresa en las mismas unidades que los datos originales.

Para comenzar a aplicar los conceptos en R, comencemos haciendo un cálculo sencillo. Veamos el rango de la edad de la muestra de pacientes.

rango_edad = max(heart$Edad) - min(heart$Edad) 
print(rango_edad)
## [1] 49

El rango de la edad por tanto es 49 años. Parece un número grande, y se podría pensar que la edad varía bastante en los pacientes de la muestra. Sin embargo, y valiéndonos nuevamente del comando hist(heart$Edad, plot=FALSE)[1:2] usado y explicado anteriormente, podremos crear una tabla de frecuencias agrupadas en intervalos

hist(heart$Edad, plot=FALSE)[1:2]
## $breaks
##  [1] 25 30 35 40 45 50 55 60 65 70 75 80
## 
## $counts
##  [1]   5  27  61 103 120 196 185 139  58  20   4

Así, nos damos cuenta que para el grueso de las observaciones la edad varía poco, pues la gran mayoría de las edades se concentran entre los 45 y 65 años, un rango de 20 años. El rango por tanto sí cumplió con su función de representar la máxima variabilidad de los datos, pero tenemos que tener cuidado con su interpretación.

Por otro lado, el rango intercuartil por lo general sí nos entrega una dispersión de datos más centralizada, a diferencia del rango. Podemos llamarla en R con la función IQR(), y nos calculará directamente qué tanto varían el 50% central de los datos.

IQR(heart$Edad)
## [1] 13

Según el rango intercuartil entonces, la dispersión de las edades dentro del 50% central de la distribución de datos es de 13 años.

Otra medida de variabilidad, usada por sus propiedades matemáticas y utilidad en comparación a las anteriores es la varianza, o mejor conocida por R por la función var(), la cual nos permite chequear rápidamente propiedades básicas de esta medida de dispersión. Observemos el siguiente comando, que involucra la variable de la presión arterial en reposo y la constante numérica 2.

var(2*heart$RepPA) == (2^2)*var(heart$RepPA)
## [1] TRUE

El output de la función, TRUE, nos confirma que si multiplicamos la constante numérica 2 por cada uno de los datos de la variable presión arterial en reposo, y calculamos su varianza, es lo mismo que multiplicar el cuadrado de la constante numérica 2 por la varianza de cada uno de los datos. Hay muchas otras propiedades que se pueden revisar fácilmente usando .

Un problema que le podríamos encontrar a la varianza es la poca interpretación que podemos obtener de los calculos realizados con ella. Por ejemplo, Dr. House difícilmente entenderá si le decimos que la varianza de la presión arterial en reposo de sus pacientes es de 342 \((mmHg)^2\). Un forma de evitar estas confusiones es usar la desviación estándar, que saca raíz cuadrada a la varianza, y por tanto, las unidades de nuestros resultados serán iguales a la unidades con la que se realizó la medida inicialmente. Basta con llamar desde R la función sd()

Como la desviación estándar trabaja de modo directo con la media de los datos, a modo de ejemplo, y para entender mejor su funcionamiento, primero tendremos que calcular la media de alguna variable a interés. Elegiremos la variable ritmo cardíaco máximo, aunque filtraremos los datos para los pacientes que sí presentan angina al realizar ejercicio usando funciones base de R

si_anginaejercicio_y_maxrc  <- subset(heart, subset = EjeAng == "S", select = MaxRC)[[1]] 
# Filtramos los pacientes que presentan angina con la función subset() y como argumentos -
# subset para filtramiento con herramientas lógicas, y select para escoger columnas, y -
# dentro de esta filtración solicitamos a R el vector numérico MaxRC con [1]

promedio_si_ang_ejercicio <- mean(si_anginaejercicio_y_maxrc) 
# Calculamos el promedio de la variable creada con la función mean(), -
# que necesita un vector numérico como argumento.

promedio_si_ang_ejercicio # Imprimimos el resultado
## [1] 125.3639

Y ahora calculemos el mismo promedio pero para los que no presentan angina, usando el paquete dplyr

heart %>% 
  filter(EjeAng == "N") %>% # Filtramos las filas que no presentan angina al hacer ejercicio
  pull(MaxRC) %>% # Le pedimos al dataframe que nos suelte la variable MaxRC como vector
  mean() 
## [1] 144.5722

Hay una diferencia considerable entre cada promedio. ¿Y qué nos dice la desviación estándar para cada caso? Veamos para el caso que sí presenta angina usando las funciones base

desv_estandar_si_angina <- sd(subset(heart, subset= EjeAng == "S", select = MaxRC)[[1]])
# Esta vez realizamos el cálculo de la desviación estándar directamente.
desv_estandar_si_angina # Imprimimos el resultado
## [1] 20.45099

Y ahora para el caso que no presenta angina usando dplyr, usando la misma lógica anterior.

heart %>% 
  filter(EjeAng == "N") %>% 
  pull(MaxRC) %>% 
  sd()
## [1] 25.6102

Es posible concluir entonces que las personas que no presentan angina al ejercitarse tienen 5 latidos más por minuto en promedio que las que si presentan angina.

Con los cálculos anteriores de promedios y desviaciones estándar, procedemos a mostrar un par de histogramas para visualizar mejor los datos obtenidos.

La barra donde se encuentra el promedio de cada uno de los casos la representamos con el color rosa, mientras que los datos que están siendo abarcados por una desviación estandar (es decir, si \(\overline{x}\) representa la media de los datos y \(\sigma\) el valor de la desviación estándar, entonces dentro de las barras de colores están los datos entre los valores \(\overline{x} - \sigma\) y \(\overline{x} + \sigma\)) están representados con las barras moradas.

Una observación a hacer, y a modo de adelanto, podemos afirmar que entre las barras coloreadas de cada caso, se debiese concentrar aproximadamente el 68% de los datos. Podremos conocer más detalles estudiando la distribución normal en las siguientes secciones.

1.2.4 Medidas de correlación

Las medidas de correlación son herramientas estadísticas utilizadas para cuantificar la relación entre dos variables. La correlación describe la dirección y la fuerza de la asociación entre las variables, lo que ayuda a comprender cómo cambian entre ellas.

Existen varias medidas de correlación, siendo la más común (y con la única que trabajaremos) la correlación de Pearson. Sin embargo, también se utilizan otras como la correlación de Spearman y la correlación de Kendall, especialmente cuando los datos no siguen una distribución normal o están en una escala ordinal.

  • Correlación de Pearson: La correlación de Pearson mide la relación lineal entre dos variables continuas. Toma valores entre -1 y 1, donde -1 indica una correlación negativa perfecta, 0 indica ninguna correlación y 1 indica una correlación positiva perfecta.
## [1] 1

En R, podemos calcular la correlación de Pearson utilizando la función cor(). Con la base de datos sobre zarigüeyas probaremos qué tanta correlación hay entre sus características físicas.

En una conversación interna con las zarigüeyas, nos dieron una pista de cuáles variables podrían llamarnos la atención para introducir como argumento a la función cor().

# Usaremos el argumento use = "complete.obs" para que cor() haga los cálculos
# a pesar de existir observaciones NA en ambas variables.
cor(possum$earconch, possum$footlgth, use = "complete.obs")
## [1] 0.7830499

El output anterior nos dice que hay una correlación positiva entre el tamaño de la oreja de la zarigüeya y el tamaño de su pata. No podemos interpretar directamente el resultado como que en promedio mientras más grande la oreja de la zarigüeya, también su pata, puesto que podrían haber otras variables que estén afectando, o podría ser simplemente casualidad. En la siguiente sección veremos ahondaremos en estas cuestiones.

¿Qué tan fuerte podría ser la correlación entre la edad y el tamaño de las orejas?

cor(possum$earconch, possum$age, use = "complete.obs")
## [1] 0.05340463

Una correlación casi nula. Quizás sea buena idea considerar el crecimiento de la oreja hasta cierta edad, pues según Google, las zarigüeyas crecen físicamente hasta aproximadamente los 2 años. Intentemos obtener esta correlación usando dplyr

possum %>% 
  select(earconch, age) %>% 
  filter(age <= 2) %>% 
  cor() %>% 
  .[1, 2]
## [1] 0.2331032

La correlación aumentó, pero poco. Puede ser que no sean suficientes datos, o en verdad no hay una relación estrecha entre ambas variables.

1.3 Tablas de frecuencia

Las tablas de frecuencia son herramientas fundamentales en estadística que nos permiten presentar la información resumida y ordenadamente por variable. Estas tablas muestran la frecuencia de ocurrencia de cada categoría o valor en una variable, así como las frecuencias acumuladas, porcentajes y otras medidas descriptivas importantes. Son especialmente útiles para analizar y comprender la distribución de datos categóricos y para realizar comparaciones entre diferentes grupos.

En R, podemos generar fácilmente tablas de frecuencia utilizando la función base table(). Esta función cuenta la frecuencia de ocurrencia de cada valor único en una variable y crea una tabla de frecuencia. Además, podemos calcular frecuencias acumuladas, porcentajes y otras medidas descriptivas utilizando funciones como cumsum(), prop.table() y otras funciones de manipulación de datos en R, que vienen incluidas en el paquete dplyr.

Creemos nuestra primera tabla de frecuencia usando table() y la base de datos possum.

table(possum$sex, useNA = "always") #S eleccionamos la variable sexo y queremos que contabilice los "NA"
## 
##    f    m <NA> 
##   43   61    0

Podemos usar más de una variable como muestra la siguiente tabla

edad_y_sexo <- table(possum$age, possum$sex)
edad_y_sexo
##    
##      f  m
##   1  3  7
##   2  7  9
##   3 11 16
##   4  6  8
##   5  6  7
##   6  7  5
##   7  1  6
##   8  0  1
##   9  2  0

O pedirle a nuestra tabla que nos muestre los datos en términos porcentuales usando la función prop.table(), y luego (usando un pipe), redondeamos los decimales a 4 mediante la función round()

prop.table(edad_y_sexo) %>% round(2) # El argumento de prop.table siempre es una tabla.
##    
##        f    m
##   1 0.03 0.07
##   2 0.07 0.09
##   3 0.11 0.16
##   4 0.06 0.08
##   5 0.06 0.07
##   6 0.07 0.05
##   7 0.01 0.06
##   8 0.00 0.01
##   9 0.02 0.00

Como en la base de datos hay aproximadamente observaciones de 100 zarigüeyas, no era tan difícil predecir los porcentajes anteriores. En su mayoría, tenemos zarigüeyas de 3 años, principalmente machos.

Necesitaremos el paquete scales para lo siguiente.

library(scales)

Crear tablas más avanzadas usando dplyr puede parecer más laborioso, pero es solamente listar los pasos mediante pipes. Creemos una tabla (en estricto rigor, un dataframe) de frecuencia para los pacientes de sexo masculino que presentan colesterol alto. La siguiente puede parecer intimidante por su extensión, pero no es más que una serie de funciones que ya hemos usado antes, y alguna extra.

heart %>%  
   filter(Sexo == "M", Col >= 200) %>% # Filtramos las filas mediante argumentos lógicos
   count(TDP) %>% # Contamos los datos de la variable TDP y nos crea una nueva columna "n"
   arrange(n) %>% # Ordenamos los valores del conteo anterior
   mutate( # Creamos nuevas columnas
     frec_absoluta_p = scales::percent(n/sum(n)), # Calculamos el porcentaje de los datos contados por count()
     frec_acumulada_p = scales::percent(cumsum(n/sum(n))) # Usamos cumsum() para la frecuencia acumulada
   ) %>% 
   select( # Cambiamos el nombre a las variables para mayor estética
    "Tipo dolor pecho" = TDP,
    "Frecuencia absoluta" = n,
    "Frecuencia absoluta porcentual"  = frec_absoluta_p,
    "Frecuencia acumulada porcentual" = frec_acumulada_p,
    )

Básicamente, filtramos a los pacientes de sexo masculino y que tenían un colesterol alto, y de ellos, contabilizamos qué tipo de dolor de pecho presenta cada uno. En base a ese conteo, calculamos las frecuencias relativas y acumuladas.

Una conclusión interesante de la tabla anterior es que entre los pacientes de sexo masculino que poseen el colesterol alto, un porcentaje no menor de ellos (40% aproximadamente), presenta algún tipo de dolor en el pecho atípico, o no presenta dolor en el pecho pero sí en otra zona del cuerpo.

1.4 Gráficos

Los gráficos nos permiten representar de forma efectiva la información de manera visual, facilitando la comprensión y la identificación de patrones, tendencias y relaciones en los datos. Hay una gran variedad de gráficos que podemos crear, presentamos algunos de ellos a continuación.

  • Gráficos de Barra : Representan datos cualitativos o comparaciones de cantidades mediante barras verticales u horizontales.

  • Histograma: Muestra la distribución de frecuencia de datos numéricos mediante barras contiguas.

  • Boxplot (Diagrama de caja y bigotes): Representa la distribución de un conjunto de datos numéricos mediante cuartiles, mediana y valores atípicos.

  • Gráfico Circular: Muestra la proporción de diferentes categorías como segmentos de un círculo.

  • Ojiva: Muestra la acumulación o frecuencia acumulativa de datos ordenados.

  • Gráfico de Dispersión: Representa la relación entre dos variables numéricas mediante puntos en un plano cartesiano.

En R, existen múltiples paquetes y funciones para crear gráficos que se adaptan a diferentes tipos de datos y objetivos de visualización. Para esta sección, usaremos exclusivamente el paquete ggplot2, el cual ofrece una sintaxis amigable para crear gráficos estéticamente agradables, personalizables y la capacidad de trabajar con datos de forma más intuitiva que las funciones base de R.

La razón de por qué el paquete ggplot2 se considera más intuitivo a comparación de las funciones base de R, es porque construimos los gráficos sobreponiendo en orden diferentes “capas” creadas por funciones del paquete, y especificando cómo se deben mapear los datos a elementos visuales.

Mostramos un resumen de cómo graficar con ggplot2:

  1. Inicialización: Para comenzar a usar ggplot2, primero debemos cargar el paquete con library(ggplot2).

  2. Sintaxis básica: La función ggplot() se usa para iniciar la construcción de un gráfico. La sintaxis básica es ggplot(datos = <datos>), donde <datos> es el conjunto de datos que deseamos visualizar.

  3. Conector de capas: Si consideramos a la misma función ggplot() como una capa, para continuar con la siguiente capa debemos conectar con un + al final de nuestra última función usada.

  4. Capas geométricas: Luego del conector, debemos agregar capas geométricas para representar la forma de nuestros datos. Tenemos varias capas a elegir:

  • geom_bar() para un gráfico de barras.
  • geom_histogram() para un gráfico de histograma.
  • geom_boxplot() para un gráfico boxplot.
  • geom_bar() con coord_polar() para un gráfico circular.
  • geom_line() para un gráfico ojiva.
  • geom_point() para un gráfico de dispersion.
  1. Estética y atributos: Podemos personalizar la apariencia del gráfico escogido especificando atributos como color, tamaño, forma, etc., dentro de las capas geométricas.

Dentro de esta capa, se suele iniciar con la función aes(), que solo en caso de ser necesario, contiene las variables que se utilizarán para definir la estructura básica del gráfico escogido. Además, se utiliza para mapear cada variable seleccionada del conjunto de datos a propiedades visuales de los elementos gráficos, como color, tamaño, forma, posición, etc. Si las personalizaciones estéticas no están dentro de esta función, significa que están siendo aplicadas a todas las variables.

Hay muchas opciones para configurar para cada capa, las cuales pueden ser encontradas en manuales específicos de ggplot2. Presentamos un breve resumen con algunas de ellas.

  • Gráfico de barras:
    • Opciones de personalización:
      • fill: Define la variable asociada al color de las barras.
      • color: Establece el color del contorno de las barras.
      • width: Controla el ancho de las barras.
    • Ejemplo de personalización:
      ggplot(data = datos) +
      geom_bar(
      aes(x = variable, y = conteo),
      fill = "categoría escogida",
      color = "black",
      width = 0.7)
  • Gráfico de histograma:
    • Opciones de personalización:
      • fill:
      • color: Color del contorno de las barras.
      • bins: Número de intervalos/bins en el histograma.
    • Ejemplo de personalización:
      ggplot(data = datos) +
      geom_histogram(
      fill = ,
      color = "black",
      bins = 30)
  • Gráfico boxplot:
    • Opciones de personalización:
      • fill: Color de las cajas y los bigotes.
      • color: Color del contorno de las cajas y los bigotes.
      • notch: Si se muestra un indicador de la mediana.
    • Ejemplo de personalización:
      ggplot(data = datos) +
      geom_boxplot(
      fill = "lightpink",
      color = "darkred",
      notch = TRUE)
  • Gráfico circular:
    • Opciones de personalización:
      • fill: Color de las secciones del gráfico circular.
      • color: Color del contorno de las secciones.
      • width: Ancho de las secciones.
    • Ejemplo de personalización:
      ggplot(data = datos) +
      geom_bar(aes(x = variable, y = valores),
      fill = "lightblue",
      color = "black", width = 0.8) +
      coord_polar(theta = "y")
  • Gráfico de ojiva:
    • Opciones de personalización:
      • color: Color de la línea.
      • size: Grosor de la línea.
    • Ejemplo de personalización: ggplot(data = datos) + geom_line(aes(x = variable, y = frecuencia_acumulada), color = "purple", size = 1.5)

  • Ejes y etiquetas: Podemos agregar ejes X e Y con xlab("Etiqueta eje X") y ylab("Etiqueta eje Y"), respectivamente. Además, podemos cambiar el título del gráfico con ggtitle("Título del gráfico").

  • Facetas: Para dividir el gráfico en subgráficos según una variable, podemos usar facetas con facet_wrap(~variable) o facet_grid(variable1 ~ variable2).

  • Escalas y temas: Podemos ajustar las escalas de los ejes con scale_x_continuous() o scale_y_continuous() y aplicar temas visuales con theme() u otros temas predefinidos.

A cada función ggplot() al paquete ggplot2, el argumento <data> debe ser un dataframe. A causa de lo anterior, el paquete dplyr será un buen complemento a ggplot2, pues los outputs de sus funciones son precisamente dataframes.

library(ggplot2) #Cargamos el paquete ggplot2

Para comenzar a graficar, diseñaremos nuestro primer gráfico de barras usando ggplot2 y la base de datos heart. Tenemos que escoger una variable cualitativa que queramos representar. Para ello, crearemos un dataframe, por medio del paquete dplyr que nos contabilice qué tipo de electrocardiograma en reposo muestran los pacientes, por lo que la variable a utilizar será RepECG. Posteriormente, introduciremos el dataframe creado a la función ggplot(),

# Creamos el dataframe que contabiliza la variable RepECG
conteo_RepECG <- heart %>% 
  count(RepECG)

# Imprimimos el dataframe
conteo_RepECG
# Creamos el gráfico de barras con ggplot2
ggplot(conteo_RepECG) +
  geom_bar(aes(x = RepECG, y = n, fill = RepECG), stat = "identity") + # Agregamos la capa geométrica de gráfico de barras. Agregamos el argumento stat = "identity"
  labs(x = "Tipo de electrocardiograma en reposo", y = "Número de observaciones", title = "Cantidad de observaciones para cada tipo de electrocardiograma en reposo") + # Añadimos etiquetas y título
  theme(axis.text.x = element_text(angle = 45, hjust = 1))  # Rotar etiquetas del eje x

Notemos que agregamos el argumento stat = "identity". la razón de por qué no lo mencionamos antes, es que la función geom_bar en ggplot2 no siempre debe llevar la función stat. La función geom_bar por defecto utiliza stat = "count", lo que significa que automáticamente cuenta el número de ocurrencias de cada categoría en el eje X y representa estas frecuencias como alturas de barras en el gráfico. También podemos especificar explícitamente una función stat diferente en geom_bar si necesitamos realizar otro tipo de manipulación de los datos antes de la representación visual. Por ejemplo, podemos usar stat = "sum" para sumar los valores de una variable dentro de cada categoría en el eje X y mostrar estas sumas como alturas de barras en el gráfico. O como acabamos de hacer, podemos usar stat = "identity" cuando ya tenemos los valores precalculados de las alturas de las barras y no necesitamos que ggplot2 realice ningún cálculo adicional.

Para el caso de la función theme(), axis.text.x = element_text(angle = 45, hjust = 1) rota las etiquetas del eje X (axis.text.x) en un ángulo de 45 grados y las alinea horizontalmente hacia la derecha (hjust = 1).

Ahora, volviendo al análisis del gráfico, es claro que el número de observaciones para el tipo de electrocardiograma en reposo “Normal” presentado por los pacientes es ampliamente superior a “LVH” y “ST”, los cuales son similares en su frecuencia.

Continuemos nuestro estudio de gráficos en ggplot2 con un histograma. Para ello, necesitamos un conjunto de datos cuantitativos, los cuales obtendremos de la variable colesterol de la base de datos heart. Realizaremos una comparativa entre la frecuencia para ciertos rangos de colesterol entre hombres y mujeres, por lo que necesitamos modelar dos histogramas. Además, ya mencionamos en la sección anterior que habían 172 observaciones para la variable colesterol iguales a 0, así que tendremos que realizar un filtrado de datos.

# Realizamos el filtrado, creando un dataframe llamado colesterol_pos
# que no contenga ningun valor de colesterol igual a 0.
colesterol_pos <- heart %>% 
  filter(Col > 0) 

# Comenzamos introduciendo como argumento a nuestra función ggplot
# el dataframe creado.
ggplot(colesterol_pos) + 
  
# Usamos las personalizaciones explicadas anteriormente.
  geom_histogram(aes(x = Col), bins = 10, fill="pink", colour="black") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  
# Especificamos los puntos de quiebre en el eje x con la opción "breaks" 
# Dentro de c va cada quiebre.
  scale_x_continuous(breaks = c(0, 90, 140, 200, 260, 320, 370, 430, 490, 550,  max(colesterol_pos$Col))) +
  
  labs(
    title = "Frecuencia de niveles de colesterol entre hombres y mujeres ",
    x = "Niveles de colesterol [mg/dL]",
    y = "Frecuencia"
  ) +

  
# Dividimos el histograma en dos, en base a la variable Sexo.
  facet_grid(Sexo ~.)

Claramente hay un intervalo donde que agrupa la mayor frecuencia para las observaciones. La gran mayoría de los pacientes, tanto sexo masculino como femenino, tenían entre 220 a 260 mg/dL de colesterol sérico cuando se les tomó la muestra. La comparación entre las frecuencias tanto para sexo masculino y femenino parece seguir una misma tendencia, pese a que tengamos un número menor de observaciones para las pacientes de sexo femenino.
Según valores de la literatura médica, los niveles óptimos de colesterol sérico son por debajo de 200 mg/dL, por lo que sería una buena idea aumentar la cantidad de barras, y así tener una mayor claridad, al tener un histograma más detallado, de los valores que se presentan en los intervalos cercanos a niveles no recomendados. Bastaría con cambiar el número de bins.

Proseguimos el estudio de gráficos con un boxplot, esta vez usando la base de datos possum. Esta base de datos posee bastantes variables numéricas que podemos escoger como nuestro conjunto de datos para graficar. Observaremos visualmente, un hecho que difícilmente podríamos haber asimilado si hubiésemos tenido todos los datos dispersos frente a nuestros ojos. Grafiquemos un boxplot entre el tamaño del pie de las zarigüeya y su ubicación.

ggplot(possum, aes(x = Pop, y = footlgth)) +
  
  # Color de llenado como fill, color de bordes como color, 
  # alpha indica la transparencia de las cajas y bigote, va de 0 a 1,
  # width controla el ancho de las cajas y el bigote.
  geom_boxplot(fill = "lightblue", color = "darkblue", alpha = 0.7, width = 0.5, notch = FALSE) +
  
  # theme_minimal: diseño de gráfico minimalista
  theme_minimal() + 
  
  # Ajustamos las etiquetas y el título.
  labs(
    title = "Distribución de longitud de pies de zarigüeya por locación",
    x = "Locación",
    y = "Longitud del pie [cm]"
  )

De este gráfico, se pueden sacar varias conclusiones útiles, por ejemplo:

El bigote inferior de cada cada (las líneas verticales debajo de cada caja) son más cortas que los bigotes superiores, lo cual significa el rango de la longitud del tamaño de los pies de las zarigüeyas, para el 25% menor de la muestra, independiente de su población, está más concentrado que el 25% mayor. La concentración de datos también la podemos apreciar dentro de la caja, ya que para el boxplot de la derecha, la linea del medio (la mediana de los datos para las zarigüeyas de la población Victoria, es decir, cuartil 2) está más cerca hacia la línea superior de la caja (cuartil 3) que para la línea inferior de la caja (cuartil 1), por lo que la longitud de los pies de las zarigüeyas de la población Victoria presenta valores más dispersos para los datos que están entre el 25% y 50%, en comparación a los datos que longitud de pies que están entre el 50% y 75%.

Además, se puede apreciar por los gráficos que claramente las zarigüeyas de nuestra muestra para la población Victoria tienen un por lo general la longitud de sus pies más grande que las zarigüeyas de New South Wales or Queensland. De igual modo, podemos observar hay casos extremos para la longitud de los pies, que son los puntos aislados en ambos boxplot.

Para finalizar nuestro estudio de gráficos, presentamos un gráfico de dispersión de datos entre el tamaño de los pies y orejas de las zarigueyas.

ggplot(data = possum, aes(x = footlgth, y = earconch)) + geom_point(alpha = 0.4, color = "blue", aes(shape = sex))

Observamos una relación lineal positiva entre ambas variables: mientras más crece una, la otra también. Si intentamos trazar una recta que mejor aproximara esta relación entre los datos, ¿qué pendiente tendría esta recta? ¿qué metodos usaríamos para que la recta este suficientemente cerca de todos los datos, y así logre reflejar una aproximación de ellos? Esas interrogantes responderemos en la sección siguiente.

2 Regresión lineal

2.1 Primeros conceptos

La regresión lineal es una técnica utilizada para modelar la relación entre una variable dependiente (o respuesta) y y una o más variables independientes X (o predictores). En nuestro caso, solamente trabajaremos con una variable independiente X, es decir, un modelo regresión lineal simple. Su objetivo es encontrar la mejor línea recta que se ajuste a los datos observados, o sea, una línea recta que minimice la distancia entre la mayoría de los puntos pertenecientes a un gráfico de dispersión, para así, predecir el valor de la variable de respuesta para valores que no sean parte de la muestra de la variable predictora.

El nombre del modelo nos entrega pistas de qué necesitamos (en esencia) para poder construir el nuestro: una relación lineal entre dos variables. Es sencillo comprobar si existe tal relación entre dos variables mediante un gráfico scatterplot.

No es difícil imaginarse que una recta podría encajar entre los puntos del gráfico de alta relación lineal entre dos variables. Sería una recta de la forma \(y=ax+b\) y algunos puntos quedarían por arriba de la línea, otros por debajo, y también algunos que sean parte de la línea. Podemos preguntarnos entonces, ¿cuál sería la pendiente de aquella línea que estuviera cercana (en términos de distancia vertical) a la mayoría de los puntos? ¿Y su intercepto? ¿cómo podríamos expresar matemáticamente esa distancia aparentemente aleatoria entre la recta y los puntos \((x_i,y_i)\)?

Existen varios métodos para conocer la ecuación de aquella recta, el más conocido quizá es el método de los mínimos cuadrados. Así, si conocemos los datos \((x_1, y_1), (x_2, y_2),..., (x_n, y_n)\) de dos variables, podemos calcular la pendiente e intercepto de una recta que mejor aproxime los datos, mientras que la aleatoridad entre la distancia de los puntos a la recta podemos expresarla mediamente una variable aleatoria \(\epsilon\), que representaría un “error” de aproximación en la distancia.

Ahora, para entrar en terreno estadístico, necesitamos una población de interés, a la cual le podemos generar un modelo de regresión lineal poblacional que sigue básicamente la misma idea anterior, aunque es más refinado en algunos detalles. Definimos el modelo poblacional de regresión simple como la ecuación \[\begin{equation*} Y = \beta_{0} + \beta_{1}X + \epsilon \ \text{,donde} \end{equation*}\]

  • \(\epsilon\) es una variable aleatoria normalmente distribuida con \(E(\epsilon) = 0\) y \(V(\epsilon) = \sigma^2\).
  • \(\beta_{0}\) y \(\beta_{1}\) son las variables de la población.

Tenemos que tener en cuenta que el modelo poblacional de regresión es una descripción idealizada entre las variables en toda la población de interés. Para hacer un ajuste de este modelo conforme a nuestros datos, necesitamos estimar los parámetros del modelo poblacional utilizando la muestra de datos que tenemos disponible para construir una aproximación de la verdadera relación entre las variables de la población, \(\beta_{0}\) y \(\beta_{1}\), y obtendríamos así el modelo muestral de regresión simple, que se representa por la ecuación \[\begin{equation*} \hat{Y}=\hat{\beta_{0}}+\hat{\beta_{1}}X+e, \end{equation*}\] donde la variable \(e\) es llamada residuo.

Vale recalcar que los errores \(\epsilon_i\) de un modelo poblacional de regresión simple son una variable aleatoria, por lo que no podemos simplemente estimarlos usando herramientas de inferencia estadística (como sí haremos para los términos \(\beta_0, \beta_1\)). Sin embargo, una aproximación a \(\epsilon_i\) con la información que tendremos disponible para nuestras investigaciones son los residuos \(e_i\), que expresan la distancia vertical entre la observación \((\hat{x}_i,\hat{y}_i)\) dada en nuestra muestra, y el punto \((x_i,y_i)\) que pertenece a la recta que mejor aproxima a todas las observaciones, por lo que \(e_i = y_i-\hat{y}_i\) sería el i-ésimo residuo.

Para el método de mínimos cuadrados, se define la suma de los residuos al cuadrado para cada observación, lo cual se conoce como \(RSS\), y posteriormente cálculo diferencial para minimizar \(RSS\), y así encontrar los coeficientes \(\hat{\beta_0}, \hat{\beta_1}\), obteniendo las siguientes ecuaciones

\[\begin{equation} \hat{\beta}_1=\frac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}\\ \hat{\beta}_0=\bar{y}-\hat{\beta}_1\bar{x} \end{equation}\]

En R, la función lm() nos creará un modelo de regresión lineal en base a nuestros datos. El output de la función summary(lm()) será una tabla que contendrá toda la información que necesitemos.

Argumentos de la función lm()

  1. formula

Este argumento especifica la relación entre la variable dependiente y las variables independientes en forma de una fórmula.

  • La fórmula generalmente tiene la forma y ~ x_1 + x_2 + …, donde y es la variable dependiente y x_1, x_2, etc., son las variables independientes.
  1. data

Se encarga de especificar el conjunto de datos que contiene las variables utilizadas en la fórmula.

  • Los nombres de las variables en la fórmula deben corresponder a las columnas en este conjunto de datos.
  1. subset

Permite especificar una condición para incluir un subconjunto de las observaciones en el modelo

  1. na.action

Determina cómo manejar los valores NA en los datos. Las opciones comunes son na.omit (eliminar las filas con NAs) o na.exclude (excluir las filas con NAs pero mantener el número de filas para la predicción posterior).

Procederemos con un ejemplo del uso de la función lm(), que se extenderá hasta la siguiente seccion Creación de un modelo. Continuaremos utilizando la base de datos de las zarigüeyas, pues sabemos, por la sección de Gráficos, que existe una relación lineal entre las variables footlgth = X y earconch = Y.

M1 <- lm(formula=footlgth ~ earconch, data=possum)
summary(M1)
## 
## Call:
## lm(formula = footlgth ~ earconch, data = possum)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2996 -1.8673 -0.0468  1.6837  7.8867 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 27.54609    3.24482   8.489 1.87e-13 ***
## earconch     0.84870    0.06708  12.653  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.747 on 101 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.6132, Adjusted R-squared:  0.6093 
## F-statistic: 160.1 on 1 and 101 DF,  p-value: < 2.2e-16

  • Intercept: El valor estimado del intercepto es 27.54609. Esto significa que se espera que la longitud del pie sea de aproximadamente 27.55 centímetros cuando la longitud de la oreja es cero. Así, \(\hat{\beta_{0}} = 27.55\)

  • earconch: El coeficiente estimado para earconch es 0.84870. Esto indica que, en promedio, por cada centímetro de aumento en la longitud de la oreja, se espera un aumento de aproximadamente 0.85 centímetros en la longitud del pie. Así \(\hat{\beta_{1}} = 0.85\)

Entonces, la recta del modelo es \(\hat{Y}= 27.55+0.85X\). Grafiquemos.

ggplot(data = possum, mapping = aes(x = footlgth, y = earconch)) +
  geom_point(color = "grey50", size = 2)  + #Gráfico de puntos para las observaciones
  geom_smooth(method = "lm", se = FALSE, color = "black")  + #Gráfico para la recta del modelo muestral
  labs(title  =  'Diagrama de dispersión junto con recta de regresión lineal', 
        x  =  'longitud de pies de  zarigüeyas [mm]', 
        y  =  'longitud de orejas de  zarigüeyas [mm]') +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")

Dentro del código del gráfico, para la opción geom_smooth(), el argumento method="lm" significa que le estamos pidiendo que grafique la mejor recta usando un modelo lineal en base a la información que introdujimos en ggplot(), y el argumento se crea un una banda de error estándar alrededor de la recta, pero aún no la usaremos.

2.2 Inferencia en parámetros y ajuste del modelo.

Test de Hipótesis para los Coeficientes.

Para evaluar la significancia de los coeficientes, se llevan a cabo pruebas t para cada parámetro del modelo:

  • Hipótesis nula (\(H_0\)): El coeficiente es igual a cero (\(\beta_i=0\)).
  • Hipótesis alternativa (\(H_a\)): El coeficiente es diferente de cero (\(\beta_i \neq 0\)).

El valor \(t\) se compara con una distribución t Student para determinar el valor p. Un valor \(p\) pequeño (generalmente < 0.05) indica que podemos rechazar la hipótesis nula y concluir que el coeficiente es significativo, lo cual para \(\beta_0\) no tiene una interpretación clara, pero para \(\beta_1\) implica que el coeficiente es significativamente diferente de cero, lo que indica que la variable independiente \(X\) tiene un efecto considerable sobre la variable dependiente \(Y\). La prueba \(t\) se calcula como:

\[\begin{equation} t= \frac{\hat{\beta_i-0}}{SE(\hat{\beta_i)}} \end{equation}\]

donde \(SE(\hat{\beta_i})\) es el error estándar de los coeficientes, y se puede interpretar como una cantidad que indica cuánto difiere nuestra estimación, en promedio, del valor verdadero del coeficiente (el valor poblacional), y se calcula como

\[\begin{equation} SE(\hat{\beta_0})^2=\sigma^2[\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}],\\ SE(\hat{\beta_1})^2=\frac{\sigma^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}, \end{equation}\]

con \(\sigma^2=Var(\epsilon)\) (asumiendo que \(Var(\epsilon)\) es igual para todas las observaciones y cada \(\epsilon_i\) no está correlacionado con \(\epsilon_j\), con \(i,j\) indexaciones para distintas observaciones.). Sin embargo, \(Var(\epsilon)\) es un valor que generalmente no es conocido, por lo que se estima su raíz cuadrada, es decir la desviación estándar de \(\epsilon\), mediante el error estándar residual (RSE) con la siguiente ecuación

\[\begin{equation} RSE=\sqrt(\frac{SRR}{n-2})=\sqrt(\frac{(e_1)^2+(e_2)^2+...+(e_n)^2}{n-2})= \end{equation}\]

Con el valor \(RSE\) podemos conocer la desviación promedio de los residuos, es decir, qué tanta es la distancia entre las observaciones de la muestra y la recta generada por el modelo. Bajo esa interpretación, podemos deducir que este valor es importante para ajustar nuestro modelo.

Ajuste del modelo.

Para el modelo de regresión lineal simple, tenemos dos herramientas que nos ayudarán a saber qué tan bueno es nuestro modelo: \(RSE\) y \(R^2\)

\(RSE\)

El error estándar residual nos dirá qué tan ajustado está nuestro modelo. Si el \(RSE\) es pequeño, significa que los valores \(\hat{y_i}\) observados, son muy próximos a los valores \(y_i\) generados por la recta de nuestro modelo, por lo que podemos asegurar que la recta del modelo (creada en base a la muestra) estaría aproximando bien a la población.

¿Qué tan pequeño debería ser el \(RSE\)? Va a depender del contexto del estudio y las unidades de la variable dependiente. Por ejemplo, en estudios médicos donde los riesgos son muy altos (la vida de un paciente), un \(RSE\) muy pequeño puede ser vital, mientras que en un estudio de comportamiento del consumidor, se puede tolerar un \(RSE\) mayor. También dependerá del tipo del objetivo del modelo a generar:

  • Modelos de Predicción: Destinados a hacer predicciones precisas, se suele necesitar un \(RSE\) pequeño.

  • Modelos Explicativos: Diseñados para entender relaciones entre variables, el \(RSE\) puede ser mayor si las relaciones encontradas son evidentes.

Por último, vale la pena mencionar que comparar el \(RSE\) del modelo actual con los \(RSE\) de modelos anteriores o de otros estudios similares puede ser una buena referencia.

\(R^2\)

Mientras que el \(RSE\) depende de las unidades de \(Y\), \(R^2\) es una alternativa que es independiente de unidades de medida.

\[\begin{equation} R^2=1-\frac{RSS}{TSS}, \end{equation}\] donde \(TSS=\sum_{i=1}^{n}(y_i-\bar{y})^2\) representa la suma total de los cuadrados.

Para un modelo de regresión lineal simple, se cumple que \(R^2=r^2\), con \(r^2=Cor(X,Y)\)

2.3 Supuestos del modelo

Las siguientes suposiciones sustentan un modelo de regresión lineal válido. Algunas ya se han mencionado anteriormente, pero la importancia de ellas hace que tengamos que puntualizarlas y revisar cada una.Muchas de ellas dependen de los residuos, por lo que son de suma importancia, pues los residuos tienen un papel fundamental a la hora de hacer inferencia. Al crear nuestro modelo, deberemos verificar que se cumplen.

  1. Linealidad entre las variables: La relación entre la variable dependiente y las variables independientes debe ser lineal. Si la relación no es lineal, las predicciones y las inferencias basadas en el modelo serán incorrectas. La linealidad la podemos verificar mediante gráficos de dispersión.

  2. Distribución normal de los residuos: La normalidad de los residuos es importante para la validación de los intervalos de confianza y las pruebas de hipótesis. Aunque la regresión lineal puede ser robusta a desviaciones de la normalidad, grandes desviaciones pueden afectar la precisión de las inferencias. Esto se puede verificar con histogramas de residuos y gráficos Q-Q.

  3. Homocedasticidad: La varianza de los errores debe ser constante a lo largo de todas las observaciones. La homoscedasticidad se puede verificar mediante gráficos de dispersión de losresiduos.

  4. Valores atípicos: Los valores atípicos (outliers) son observaciones que se desvían significativamente de la tendencia general de los datos. Estos pueden influir fuertemente en los resultados del modelo de regresión lineal, afectando la estimación de los coeficientes, los errores estándar y, en última instancia, las inferencias que hacemos a partir del modelo.Debemos identificarlos y tenerlos en cuenta a la hora de realizar nuestro modelo.

  5. Independencia de los residuos: Los errores de la predicción (residuos) deben ser independientes entre sí. La independencia de los errores nos asegura que no hay patrones sistemáticos en los residuos que puedan sesgar los resultados.

2.4 Creación de un modelo

El ejemplo siguiente, será construido con la base de datos de las zarigüeyas. Procedemos a verificar qué tan bien está ajustado nuestro modelo y si se cumplen los supuestos anteriores

  1. Verificar si un modelo lineal tiene sentido.

Para el punto 2.1 de este apartado, finalizamos graficando un scatterplot para ambas variables, notando que existe una relación lineal entre ambas. Notemos que

Además, podemos realizar el siguiente test de correlación

cor.test(x = possum$footlgth, y = possum$earconch, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  possum$footlgth and possum$earconch
## t = 12.653, df = 101, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6948175 0.8480600
## sample estimates:
##       cor 
## 0.7830499

Obteniendo así, una correlación de \(0.78\) con un p-valor muy pequeño.

  1. Chequear qué tan bien ajustado está nuestro modelo

Al finalizar la sección 2.1, en la tabla obtenida por la función summary(M1), aparece el valor de Multiple R-squared siendo este de \(0.6132\), un valor moderadamente alto, lo cual significa el 61.32% de la variabilidad en la longitud de las orejas de las zarigüeyas (\(Y\)) puede ser explicada por la longitud de sus pies (\(X\)). El restante 38.68% de la variabilidad en la longitud de las orejas no es explicada por la longitud de los pies, lo cual hace indicar que hay otros factores que afectan la longitud de las orejas que no están incluidos en nuestro modelo.

  1. Verificación de los supuestos del modelo

Tenemos cinco supuestos que verificar, y lo haremos graficando y haciendo algunos tests. Para comprobar si se cumplen los supuestos de linealidad entre las variables X e Y, la homocedasticidad e identificar los outliers, haremos un gráfico de los residuos para cada observación.

Gráfico de residuos para cada observación

ggplot(data = possum, aes(x = prediccion, y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
  geom_segment(aes(xend = prediccion, yend = 0), alpha = 0.2) +
  geom_smooth(se = FALSE, color = "firebrick") +
  labs(title = "Distribución de los residuos", x = "predicción modelo",
       y = "residuo") +
  geom_hline(yintercept = 0) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Notemos que el gráfico de residuos presenta una curva de tendencia que se mueve en sintonía a los valores de los residuos. En un buen ajuste de regresión lineal, esta línea debería ser aproximadamente recta y horizontal en torno a la línea de cero.

  • Linealidad: Si los residuos muestran un patrón sistemático (por ejemplo, una curva, una forma en U, o aumento o disminución constante), esto indica que la relación entre \(X\) e \(Y\) no es lineal. En nuestro caso, los residuos para cada observación se distribuyen de forma aleatoria a lo largo de la recta de regresión.

  • Homocedasticidad: La homocedasticidad es la condición en la que los residuos tienen una varianza constante. Si los residuos se dispersan de manera desigual (por ejemplo, formando una especie de embudo), esto sugiere heterocedasticidad. En nuestro caso, los resiudos se están dispersando de manera aleatoria en torno al 0.

  • Outliers:. Los outliers son puntos que se desvían significativamente del resto de los datos y pueden distorsionar los resultados. En nuestro caso, se evidencian algunos puntos que se están desviando de los demás valores, por lo que tendremos que observarlos con cuidado.

Hay un leve patrón “negativo” para las longitudes de pies menores a 70, y “positivo” para mayores a 70”

Distribución normal de los residuos

Para chequear la distribución normal de los residuos, haremos un histograma de densidad y posteriormente un gráfico QQPlot: si nuestros residuos están muy cerca de la recta diagonal que genera el gráfico, entonces siguen una distribución aproximadamente normal. ¿Por qué un histograma de densidad y no de frecuencia? La respuesta es que un histograma de densidad es más adecuado para visualizar la distribución subyacente de los datos.

ggplot(data = possum, aes(x = residuos)) +
  geom_histogram(aes(y = after_stat(density))) + # con after_stat(density) pedimos un histograma de densidad. 
  labs(title = "Histograma de los residuos") +
  theme_light()

Si pusiéramos una curva por encima del histograma, su forma sería parecida a una campana.

Podemos hacer otra comprobación para ver si los residuos distribuyen aproximadamente normal, y es con la función qqnorm(), la cual crea un gráfico de probabilidad normal (o también llamado Q-Q plot) que compara los cuantiles de los residuos con los cuantiles de una distribución normal estándar. Si los residuos se distribuyen aproximadamente de forma normal, los puntos deberían alinearse aproximadamente a lo largo de una línea recta. Tal línea recta de referencia también se puede añadir al gráfico con la función qqline(), y se calcula utilizando la media y la desviación estándar de los residuos. Ayuda a visualizar más fácilmente cualquier desviación de la normalidad en la distribución de los residuos.

qqnorm(M1$residuals) 
qqline(M1$residuals)

Como los residuos están alineados a lo largo de la línea de referencia, podemos concluir que distribuyen normal.

Autocorrelación de residuos

Para el siguiente análisis, graficaremos los residuos pero indexados en función del orden en el que aparecen los datos, ya que cualquier patrón secuencial o temporal se vuelve evidente (podrían existir patrones relacionados con el orden de recogida de datos.). En muchas aplicaciones, las observaciones tienen un orden natural (por ejemplo, en datos recolectados secuencialmente). En caso de haber autocorrelación (donde los residuos están correlacionados con los residuos en índices cercanos), se manifestaría como un patrón no aleatorio en el gráfico.

ggplot(data = possum, aes(x = seq_along(residuos), y = residuos)) + 
# seq_along(residuos) genera una secuencia de índices (1, 2, 3, ...) para los residuos.
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
  geom_line(linewidth = 0.3) + 
# geom_line conecta los puntos con una línea delgada (qué tan delgada lo establece linewidth), lo que ayuda a visualizar tendencias a lo largo del tiempo o el índice.
  labs(title = "Distribución de los residuos", x = "index", y = "residuo") +
  geom_hline(yintercept = 0) + # Añade una línea horizontal en y = 0 para referencia.
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")

Los residuos están incorrelacionados pues se distribuyen aleatoriamente por encima y por debajo del 0.

Valores influyentes

Utilicemos el siguiente test para revisar rápidamente valores influyentes que tienen potencial para afectar nuestro modelo.

library(car)
outlierTest(M1, cutoff=Inf)
##     rstudent unadjusted p-value Bonferroni p
## 42 -3.178290          0.0019717      0.20309
## 59  3.010265          0.0033042      0.34033
## 55  2.582486          0.0112570           NA
## 11  2.290125          0.0241140           NA
## 21  2.040958          0.0438900           NA
## 65  1.978108          0.0506680           NA
## 74 -1.854392          0.0666300           NA
## 54  1.744105          0.0842140           NA
## 87  1.636563          0.1048700           NA
## 86 -1.607680          0.1110600           NA

Por los resultados del test, los valores 42 y 59 son sospechosos. ¿Qué tanto cambiarían los coeficientes del modelo si sacamos esas observaciones?

lm(formula = footlgth ~ earconch, data = possum)$coefficients
## (Intercept)    earconch 
##  27.5460934   0.8487005
lm(formula = footlgth ~ earconch, data = possum[c(-42,-59),])$coefficients
## (Intercept)    earconch 
##  25.9746406   0.8813786

Graficamos ambas rectas para comparar

ggplot(data = possum, mapping = aes(x = footlgth, y = earconch)) +
  geom_point(color = "grey50", size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  #se resalta el valor excluido
  geom_point(data = possum[c(42,59),], color = "red", size = 2) +
  #se añade la nueva recta de mínimos cuadrados
  geom_smooth(data = possum[c(-42,-59),], method = "lm", se = FALSE, color = "blue") +
  labs(title  =  'Diagrama de dispersión', x  =  'número  de bateos') +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")

En realidad no cambia mucho la recta. Aquellas observaciones no afectan en mayor medida. De hecho, puede que hagan nuestro modelo aún mejor.

par(mfrow = c(1,2))
plot(M1)

La dispersión de los datos se mide mediante la desviación residual. Aquí, el parámetro de dispersión es aproximadamente 7.55, lo que indica que hay cierta variabilidad en la longitud del pie que no está explicada por la longitud de la oreja. Luego, ϵ es un error que sigue una distribución normal con media cero y varianza constante \((V(ϵ)=7.5)\)

2.5 Predicción

3 Anexo

3.1 Instalación de paquetes

Para cargar un paquete, solo debemos usar la función library(*paquete a instalar*)

library(dplyr)

3.1.1 Uso básico del pipe (%>%)

El pipe se utiliza para encadenar operaciones de manera secuencial, pasando el resultado de una operación como entrada a la siguiente operación. La sintaxis básica del pipe es la siguiente:

objeto_inicial %>% operacion1() %>% operacion2() %>% ... %>% operacionN()

En esta sintaxis:

  • objeto_inicial es el objeto inicial que deseamos transformar o manipular.
  • operacion1(), operacion2(), …, operacionN() son las operaciones que deseamos aplicar al objeto inicial en secuencia.

Ejemplo de uso:

# Creamos un vector de números
numeros <- c(1, 2, 3, 4, 5)

# Utilizamos el pipe para aplicar operaciones secuenciales
resultado <- numeros %>%
  sum() %>%           # Calculamos la suma de los números
  sqrt() %>%          # Calculamos la raíz cuadrada del resultado
  round(digits = 2)   # Redondeamos el resultado a 2 decimales

print(resultado)      # Imprimir el resultado
## [1] 3.87