Modelos lineales
E. Martinez y E. Maure
2022-10-05
Introducción
Un modelo es una representación formal de un fenómeno, una reducción de dimensionalidad que posee utilidad práctica. Dicha representación normalmente puede ser condensada en una expresión matemática, una fórmula, que indica cómo una variable se relaciona con otra(s). Empíricamente, el paradigma se basa en estudiar la relación matemática entre variables aleatorias respuesta, con una distribución de probabilidades dada, y aquellas variables que la predicen, con el fin de explicar asociaciones entre variables y realizar inferencia. De modo muy general, podemos escribir:
\(y = f(x)\)
Donde \(y\) es la variable que modelamos, \(f\) es una función de una o múltiples variable(s) explicatoria(s).
¿Qué es un modelo lineal?
Un modelo lineal suele escribirse como: \[y=\beta\ _0 + \beta\ _1 x_1 + ... + \beta\ _m x_m \] donde el efecto de cada \(\beta\ \) debe interpretarse como el cambio en \(y\) dado por un cambio unitario en la variable predictora asociada a ese coeficiente, siempre que las demás variables \(x\) se mantengan constantes. Además, la ordenada al origen es la media general del modelo (cuando ninguna de las variables está asociada a la respuesta).
El modelo que presentamos es determinista ya que el valor de \(y\) para cada uno de los valores de \(x\) es fijo. Pero en la realidad cada uno de los valores de \(y\) se ve afectado por un error aleatorio que debe agregarse al modelo: \[y_i = \beta\ _0 + \beta\ _1 x_1i + ... + \beta\ _m x_mi + \varepsilon_i\] debido a que existen diferencias entre las unidades observacionales o experimentales.
Análisis de la correlación
Gráfico de dispersión
Antes del plantear una relación funcional entre las variables en estudio es necesario un análisis descriptivo del comportamiento de las variables. En nuestro caso, seguiremos trabajando con la base penguns_size.csv. Para este análisis es necesario tomar sólo aquellas variables que sean cuantitativas.
Para observar la asociación entre las variables dela base en estudio, haremos un plot:
penguins <- read.csv("penguins_size.csv")
head(penguins, n=8) species island culmen_length_mm culmen_depth_mm flipper_length_mm
1 Adelie Torgersen 39.1 18.7 181
2 Adelie Torgersen 39.5 17.4 186
3 Adelie Torgersen 40.3 18.0 195
4 Adelie Torgersen NA NA NA
5 Adelie Torgersen 36.7 19.3 193
6 Adelie Torgersen 39.3 20.6 190
7 Adelie Torgersen 38.9 17.8 181
8 Adelie Torgersen 39.2 19.6 195
body_mass_g sex
1 3750 MALE
2 3800 FEMALE
3 3250 FEMALE
4 NA <NA>
5 3450 FEMALE
6 3650 MALE
7 3625 FEMALE
8 4675 MALEpenguins = na.omit(penguins) # Para eliminar las filas con NA
penguins$sex[penguins$sex=="."]="MALE"
plot(penguins)
Tomaremos dos variables: largo de aleta y masa corporal. En principio no tendremos en cuenta ninguna clasificación de especie, hábitat o sexo.
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.4.1
✔ readr 2.1.2 ✔ forcats 0.5.2
✔ purrr 0.3.4
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
Loading required package: carData
Attaching package: 'car'
The following object is masked from 'package:dplyr':
recode
The following object is masked from 'package:purrr':
some
Attaching package: 'rstatix'
The following object is masked from 'package:stats':
filterggplot(data = penguins) +
geom_point(mapping = aes(x = flipper_length_mm, y = body_mass_g))+
ggtitle("Diagrama de dispersión del largo de aleta (mm) y la masa corporal (g) de los pinguinos en estudio")
ggplot(data = penguins) +
geom_point(mapping = aes(x = flipper_length_mm, y = body_mass_g,color=species))
ggplot(data = penguins) +
geom_point(mapping = aes(x = flipper_length_mm, y = body_mass_g,
shape=species))
ggplot(data = penguins) +
geom_point(mapping = aes(x = flipper_length_mm,
y = body_mass_g,
color=species))+
facet_wrap(~species, nrow = 3)
ggplot(data = penguins) +
geom_point(mapping = aes(x = flipper_length_mm,
y = body_mass_g,
color=species))+
facet_wrap(species~sex, nrow = 2)
ggplot(data = penguins) +
geom_point(mapping = aes(x = flipper_length_mm,
y = body_mass_g,
color=species))+
facet_grid(sex~species)
Coeficiente de correlación simple lineal
Para el cálculo del coeficiente de correlación simple lineal muestral, se utilizan los comandos base de R:
r1=cor(penguins$body_mass_g, penguins$flipper_length_mm)
r1[1] 0.873211La función cor() tiene varios argumentos, entre ellos el método para el cálculo del coeficiente de CSL muestral que por defecto utiliza Pearson. Recordar que siempre puede utilizarse la ayuda para cada función, librería o paquete.
Al interpretar el valor de \(r_1\) podemos decir que exite una correlación simple lineal fuerte entre el largo de la aleta y la masa corporal de los pingüinos.
Prueba parael coeficiente de CSL
\(H_c : \quad\text{"Existe correlación simple lineal entre la masa corporal (g) y el largo de aleta (mm) de los pingüinos en estudio"}\)
\(H_0 : \rho =0\)
\(H_1 : \rho \neq 0\)
\(\alpha = 0,05\)
Con la prueba de hipótesis para el Coeficiente de Correlación Simple Lineal \(\rho\) también sale el cálculo de la estimación:
r1_test = cor.test(penguins$body_mass_g,penguins$flipper_length_mm)
r1_test
Pearson's product-moment correlation
data: penguins$body_mass_g and penguins$flipper_length_mm
t = 32.648, df = 332, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.8450881 0.8965147
sample estimates:
cor
0.873211 Por lo obtenido, se rechaza la hipótesis nula bajo un \(\alpha = 0,05\). Existe evidencia muestral suficiente para decir que existe asociación entre el largo de aleta (mm) y la masa corporal (g) de los pingüinos en estudio.
Análisis de Regresión simple lineal
Tenemos evidencia que puede existir una relación funcional entre las variables en estudio. Por ello se plantea el modelo: \[Y_i = \beta _0 + \beta _1 x_i + \varepsilon_i i=1,...,n\] Para la estimación de los parámetros y las pruebas de hipótesis correspondientes, se utiliza:
modelo1 = lm(body_mass_g~flipper_length_mm,
data=penguins)
modelo1
Call:
lm(formula = body_mass_g ~ flipper_length_mm, data = penguins)
Coefficients:
(Intercept) flipper_length_mm
-5865.82 50.12 En la salida principal del modelo tenemos las estimaciones de los parámetros. (Interpretar) Con ellos podemos graficar la recta que ajusta a los datos:
g1 = ggplot(penguins, aes(flipper_length_mm,body_mass_g))+
geom_smooth(method="lm", se=FALSE, color="blue")+
geom_point(alpha = 0.5) +
theme(plot.background = element_rect(colour = NA))+
xlab("Largo de aleta (mm)")+
ylab("Masa corporal (g)")
g1`geom_smooth()` using formula 'y ~ x'
Para las pruebas debemos llamar a los diferentes elementos que se generaron con la función lm() en el objeto lista “modelo1”.
names(modelo1) [1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "xlevels" "call" "terms" "model" Para realizar la inferencia de los coeficientes:
summary(modelo1)
Call:
lm(formula = body_mass_g ~ flipper_length_mm, data = penguins)
Residuals:
Min 1Q Median 3Q Max
-1057.23 -259.15 -11.13 242.86 1293.73
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5865.818 309.340 -18.96 <2e-16 ***
flipper_length_mm 50.120 1.535 32.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 392.8 on 332 degrees of freedom
Multiple R-squared: 0.7625, Adjusted R-squared: 0.7618
F-statistic: 1066 on 1 and 332 DF, p-value: < 2.2e-16Para la del modelo, es decir teniendo en cuenta la partición de la variabilidad del modelo:
anova(modelo1)Analysis of Variance Table
Response: body_mass_g
Df Sum Sq Mean Sq F value Pr(>F)
flipper_length_mm 1 164474102 164474102 1065.9 < 2.2e-16 ***
Residuals 332 51230376 154308
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Asimismo, la función confint() permite obtener de manera muy simple intervalos de confianza para los parámetros de la regresión:
confint(modelo1) 2.5 % 97.5 %
(Intercept) -6474.33214 -5257.30369
flipper_length_mm 47.10013 53.13991Supuestos
Cuando construimos modelos de este tipo, asumimos ciertas cosas. Los supuestos principales en este caso son:
-Independencia entre observaciones. -Homocedasticidad. -Los residuos son normales.
Vimos anteriormente cómo evaluar el cumplimiento de los supuestos. Otra herramienta muy utilizada es con el análisis de los residuos.
Residuos vs Predichos (Residuals vs Fitted)
En este gráfico buscamos ver la dispersión respecto de la recta para cada valor predicho (cuánto se alejan de nuestro ajuste). Esperamos no observar ningún tipo de patrón en los residuos. Esperamos no ver datos atípicos (datos con residuos muy grandes).
predichos1 = modelo1$fitted.values
residuos1 = modelo1$residuals
res_stud1 = rstandard(modelo1)Grafiquemos
ggplot(data=modelo1)+
geom_point(mapping=aes(x=modelo1$fitted.values, y=modelo1$residuals))+
xlab("Predichos")+
ylab("Residuos")
Q-Q Plot
Este gráfico muestra cómo se acumulan los residuos respecto de los cuantiles teóricos de una distribución normal. Si la distribución de residuos es normal, los veremos cercanos a la recta. Desviaciones de la recta indican que la distribución de los residuos no es normal.
ggplot(modelo1, aes(sample = modelo1$residuals))+
stat_qq()+
stat_qq_line()+
xlab("Cuartiles normales")+
ylab("Cuartiles residuos")
#### Residuos estandarizados vs Predichos
Este gráfico es similar al primero, pero los residuos están estandarizados. Esperamos no ver ningún patrón, con los residuos distribuidos normalmente alrededor de cero.
ggplot(data=modelo1)+
geom_point(mapping=aes(x=modelo1$fitted.values, y=rstandard(modelo1)))+
geom_hline(yintercept = c(-3,0,3),col="red")+
xlab("Predichos")+
ylab("Residuos estandarizados")
Con todos los análisis realizados se puede decir que el modelo1 es un buen ajuste a los datos observados en los pingüinos
Ajustamiento de modelos intrínsecamente lineales
De la misma manera que se evalúa el modelo lineal lo realizamos de manera explícita en la función.
Modelo polinómico de 2° grado
Teniendo en cuenta que el modelo es. \[Y_i = \beta _0 + \beta _1 x_i + \beta _2 x_i^2 + \varepsilon _i\] La función se escribe:
modelo2 = lm(penguins$body_mass_g~penguins$flipper_length_mm + I(penguins$flipper_length_mm^2))Modelo potencial
Teniendo en cuenta que el modelo es. \[Y_i = \beta _0 * x_i^{\beta_1}*\varepsilon _i\] Aplicamos logarirmo en ambos miembros del modelo y obtenemos la linealización del modelo exponencial: \[ln(Y_i) = ln(\beta _0) +ln(x_i)\beta_1+ln(\varepsilon _i)\] \[Y_i^* = \beta _0^* +\beta_1 x_i^*+\varepsilon _i^*\] La función se escribe:
modelo3 = lm(log(penguins$body_mass_g)~log(penguins$flipper_length_mm))Modelo exponencial
Teniendo en cuenta que el modelo es. \[Y_i = \beta _0 * e^{\beta_1x_i}*\varepsilon _i\] Aplicamos logarirmo en ambos miembros del modelo y obtenemos la linealización del modelo exponencial: \[ln(Y_i) = ln(\beta _0) +\beta_1x_i+ln(\varepsilon _i)\] \[Y_i^* = \beta _0^* +\beta_1x_i+\varepsilon _i^*\] La función se escribe:
modelo4 = lm(log(penguins$body_mass_g)~penguins$flipper_length_mm)