library(readr)
library(GGally)
library(ggplot2)
library(sqldf)
library(scatterplot3d)
facebook <- read_delim('dataset_Facebook.csv', delim=';')Análisis de likes en Facebook
Introducción
En este análisis, trabajaremos con una base de datos de Facebook que contiene información sobre publicaciones, incluyendo el número de likes, comentarios, compartidas, alcance total, impresiones y más. El objetivo es explorar las relaciones entre estas variables y desarrollar un modelo de regresión lineal para predecir el número de likes en función de otras métricas.
Carga de datos y preprocesamiento
Echemos un vistazo a la estructura de la base de datos:
str(facebook)spc_tbl_ [500 × 19] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ Page total likes : num [1:500] 139441 139441 139441 139441 139441 ...
$ Type : chr [1:500] "Photo" "Status" "Photo" "Photo" ...
$ Category : num [1:500] 2 2 3 2 2 2 3 3 2 3 ...
$ Post Month : num [1:500] 12 12 12 12 12 12 12 12 12 12 ...
$ Post Weekday : num [1:500] 4 3 3 2 2 1 1 7 7 6 ...
$ Post Hour : num [1:500] 3 10 3 10 3 9 3 9 3 10 ...
$ Paid : num [1:500] 0 0 0 1 0 0 1 1 0 0 ...
$ Lifetime Post Total Reach : num [1:500] 2752 10460 2413 50128 7244 ...
$ Lifetime Post Total Impressions : num [1:500] 5091 19057 4373 87991 13594 ...
$ Lifetime Engaged Users : num [1:500] 178 1457 177 2211 671 ...
$ Lifetime Post Consumers : num [1:500] 109 1361 113 790 410 ...
$ Lifetime Post Consumptions : num [1:500] 159 1674 154 1119 580 ...
$ Lifetime Post Impressions by people who have liked your Page : num [1:500] 3078 11710 2812 61027 6228 ...
$ Lifetime Post reach by people who like your Page : num [1:500] 1640 6112 1503 32048 3200 ...
$ Lifetime People who have liked your Page and engaged with your post: num [1:500] 119 1108 132 1386 396 ...
$ comment : num [1:500] 4 5 0 58 19 1 3 0 0 3 ...
$ like : num [1:500] 79 130 66 1572 325 ...
$ share : num [1:500] 17 29 14 147 49 33 27 14 31 26 ...
$ Total Interactions : num [1:500] 100 164 80 1777 393 ...
- attr(*, "spec")=
.. cols(
.. `Page total likes` = col_double(),
.. Type = col_character(),
.. Category = col_double(),
.. `Post Month` = col_double(),
.. `Post Weekday` = col_double(),
.. `Post Hour` = col_double(),
.. Paid = col_double(),
.. `Lifetime Post Total Reach` = col_double(),
.. `Lifetime Post Total Impressions` = col_double(),
.. `Lifetime Engaged Users` = col_double(),
.. `Lifetime Post Consumers` = col_double(),
.. `Lifetime Post Consumptions` = col_double(),
.. `Lifetime Post Impressions by people who have liked your Page` = col_double(),
.. `Lifetime Post reach by people who like your Page` = col_double(),
.. `Lifetime People who have liked your Page and engaged with your post` = col_double(),
.. comment = col_double(),
.. like = col_double(),
.. share = col_double(),
.. `Total Interactions` = col_double()
.. )
- attr(*, "problems")=<externalptr> La base de datos original contiene 500 observaciones y 19 variables. Renombremos algunas variables para facilitar su manejo:
names(facebook)[1] <- 'Page_total_likes'
names(facebook)[5] <- 'Post_Weekday'Ahora, filtremos las publicaciones con menos de 1500 likes para enfocarnos en el rango más común:
facebook <- sqldf('SELECT like, Page_total_likes, Post_Weekday, comment, share
FROM facebook WHERE like < 1500')Visualización de relaciones entre variables
Antes de construir un modelo, es útil visualizar las relaciones entre las variables de interés. Usemos ggpairs para crear un gráfico de pares:
ggpairs(facebook, lower = list(continuous = wrap('points', colour = "blue")),
diag = list(continuous = wrap("barDiag", colour = "red")))
El gráfico de pares muestra las relaciones entre las variables seleccionadas. Podemos observar una correlación positiva entre el número de likes y el número de comentarios y compartidas. También hay una ligera relación positiva con el total de likes de la página.
Basándonos en estas observaciones, procedamos a construir un modelo de regresión lineal para predecir el número de likes.
Modelo de regresión lineal múltiple
Ajustemos un modelo de regresión lineal múltiple considerando el total de likes de la página, el día de la semana de la publicación, y el número de comentarios y compartidas como predictores:
modelo <- lm(like ~ Page_total_likes + Post_Weekday + comment + share, data = facebook)
summary(modelo)
Call:
lm(formula = like ~ Page_total_likes + Post_Weekday + comment +
share, data = facebook)
Residuals:
Min 1Q Median 3Q Max
-326.80 -42.06 -14.07 17.55 1050.36
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.119e+02 3.897e+01 -2.873 0.004250 **
Page_total_likes 1.038e-03 2.979e-04 3.486 0.000535 ***
Post_Weekday -2.435e+00 2.382e+00 -1.022 0.307275
comment 3.188e+00 5.726e-01 5.567 4.3e-08 ***
share 5.202e+00 2.610e-01 19.929 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 106.7 on 484 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.5946, Adjusted R-squared: 0.5913
F-statistic: 177.5 on 4 and 484 DF, p-value: < 2.2e-16El resumen del modelo nos proporciona varias piezas clave de información:
- Los coeficientes estimados para cada predictor, junto con sus errores estándar, valores t y valores p.
- El R-cuadrado y el R-cuadrado ajustado, que indican la proporción de varianza en los likes explicada por el modelo.
- El estadístico F y su valor p, que prueban la significancia global del modelo.
Interpretemos estos resultados:
- El modelo explica alrededor del 59.5% de la variabilidad en el número de likes (R-cuadrado ajustado = 0.5913).
- Las variables más significativas son el número de comentarios y compartidas (p < 0.001).
- El total de likes de la página también es significativo (p < 0.001), pero su coeficiente es muy pequeño, lo que sugiere un efecto limitado.
- El día de la semana de la publicación no parece tener un efecto significativo (p = 0.307).
Estos hallazgos nos dan una buena idea de qué factores influyen en el número de likes. Sin embargo, antes de finalizar, realicemos algunos chequeos adicionales en el modelo.
Modelo de regresión lineal sin intercepto
Ajustemos ahora un modelo de regresión lineal sin intercepto, considerando solo el número de comentarios y compartidas como predictores:
modelo <- lm(like ~ comment + share - 1, data = facebook)
summary(modelo)
Call:
lm(formula = like ~ comment + share - 1, data = facebook)
Residuals:
Min 1Q Median 3Q Max
-353.61 -36.54 -5.86 16.71 1055.82
Coefficients:
Estimate Std. Error t value Pr(>|t|)
comment 3.4322 0.5751 5.968 4.63e-09 ***
share 5.2668 0.2026 25.995 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 107.9 on 487 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.7705, Adjusted R-squared: 0.7696
F-statistic: 817.7 on 2 and 487 DF, p-value: < 2.2e-16El resumen de este modelo muestra:
- Ambas variables, comentarios y compartidas, son altamente significativas (p < 0.001).
- Los coeficientes estimados son positivos, lo que indica una relación directa con el número de likes.
- El modelo explica alrededor del 77% de la variabilidad en los likes (R-cuadrado ajustado = 0.7696), una mejora sustancial respecto al modelo anterior.
Calculemos los intervalos de confianza para los coeficientes:
confint(modelo) 2.5 % 97.5 %
comment 2.302226 4.562144
share 4.868679 5.664873Y realicemos una prueba ANOVA para evaluar la contribución de cada predictor:
anova(modelo)Analysis of Variance Table
Response: like
Df Sum Sq Mean Sq F value Pr(>F)
comment 1 11167718 11167718 959.60 < 2.2e-16 ***
share 1 7864024 7864024 675.72 < 2.2e-16 ***
Residuals 487 5667667 11638
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La tabla ANOVA descompone la variabilidad en los likes en componentes atribuibles a cada predictor. Ambos, comentarios y compartidas, contribuyen significativamente (p < 0.001).
Gráficos de diagnóstico del modelo
Por último, verifiquemos los supuestos del modelo examinando los gráficos de diagnóstico:
plot(modelo)
Los gráficos de diagnóstico sugieren que:
- Los residuos parecen distribuirse aleatoriamente alrededor de 0, sin patrones claros, lo que indica linealidad.
- El gráfico Q-Q muestra que los residuos se alinean razonablemente bien con la línea, sugiriendo normalidad.
- El gráfico Scale-Location muestra una dispersión relativamente constante de los residuos, indicando homocedasticidad.
- No hay puntos de alto apalancamiento que sean altamente influyentes.
En general, el modelo parece cumplir adecuadamente con los supuestos de regresión lineal.
Gráfico 3D de la relación entre likes, comentarios y compartidas
Para visualizar la relación entre el número de likes, comentarios y compartidas en las publicaciones de Facebook, creemos un gráfico 3D utilizando el paquete scatterplot3d:
facebook2 <- sqldf('SELECT *
FROM facebook
WHERE like < 1500 AND comment < 100')
fb <- data.frame(facebook2$comment, facebook2$share, facebook2$like)
library(scatterplot3d)
s3d <- scatterplot3d(fb, type = "h", color = "steelblue",
angle = 65, scale.y = 0.7, pch = 16,
main = "Relación del número de likes vs comentarios y compartidas",
xlab = "Comment", ylab = "Share", zlab = "Likes")
s3d$plane3d(modelo, lty = "dotted")
El gráfico 3D muestra cada publicación como un punto en el espacio, con el número de comentarios en el eje x, el número de compartidas en el eje y, y el número de likes en el eje z.
Algunas observaciones clave:
La nube de puntos tiene una forma aproximadamente plana y se inclina hacia arriba, indicando una relación positiva entre las tres variables. A medida que aumentan los comentarios y las compartidas, también tienden a aumentar los likes.
La densidad de puntos es mayor en la región con valores bajos de comentarios y compartidas, sugiriendo que la mayoría de las publicaciones tienen un engagement relativamente bajo.
Hay algunos puntos atípicos con un número muy alto de likes en relación con los comentarios y las compartidas. Estas publicaciones podrían ser contenido viral o publicaciones promocionadas.
El plano de regresión ajustado (representado por la línea punteada) pasa por el centro de la nube de puntos, capturando la tendencia general de los datos. La inclinación del plano confirma la relación positiva entre las variables.
La dispersión de los puntos alrededor del plano de regresión sugiere que, aunque el modelo lineal captura la tendencia general, hay una variabilidad sustancial que no se explica por los comentarios y las compartidas solamente. Otros factores no incluidos en el modelo podrían influir en el número de likes.
Este gráfico 3D es una poderosa herramienta de visualización que nos permite comprender intuitivamente la relación entre estas tres métricas clave de engagement. Confirma visualmente los hallazgos del análisis de regresión lineal, mostrando que las publicaciones con más comentarios y compartidas tienden a tener más likes.
Sin embargo, también destaca la complejidad de predecir el engagement en las redes sociales. Mientras que el modelo de regresión lineal captura la tendencia general, hay una considerable dispersión alrededor de esta tendencia, indicando que otros factores también juegan un papel.
En conclusión, este gráfico 3D refuerza la idea de que fomentar las interacciones de los usuarios, como comentarios y compartidas, puede ser una estrategia efectiva para aumentar los likes y el engagement general. Sin embargo, también subraya la necesidad de considerar un rango más amplio de factores, como el tipo de contenido, el tiempo de publicación y las características de la audiencia, para desarrollar una comprensión completa de lo que impulsa el engagement en Facebook.
Predicción con nuevos datos
Ahora que tenemos un modelo de regresión lineal ajustado, podemos usarlo para predecir el número de likes para nuevas publicaciones. Creemos un nuevo marco de datos con valores hipotéticos para el número de comentarios y compartidas:
newpage <- data.frame(comment = 265, share = 153)Luego, usemos la función predict() para obtener una predicción del número de likes para esta nueva publicación, junto con un intervalo de predicción del 90%:
predict(modelo, newdata = newpage, interval = 'prediction', level = 0.90) fit lwr upr
1 1715.346 1431.998 1998.694El resultado muestra que, para una publicación con 265 comentarios y 153 compartidas, el modelo predice aproximadamente 1715 likes, con un intervalo de predicción del 90% de [1432, 1999].
Este intervalo de predicción es más amplio que un intervalo de confianza para la media, ya que tiene en cuenta tanto la incertidumbre en los coeficientes del modelo como la variabilidad residual. Indica que hay un 90% de probabilidad de que el número real de likes para una publicación con estas características caiga dentro de este rango.
Es importante tener en cuenta que esta predicción está sujeta a todas las limitaciones y supuestos del modelo de regresión lineal. Asume que la relación entre las variables es lineal y que el nuevo punto de datos sigue el mismo patrón que los datos utilizados para ajustar el modelo.
Visualización de los valores reales y estimados
Para comparar visualmente los valores reales de likes con los valores estimados por el modelo, podemos crear un gráfico de líneas. Primero, creemos un nuevo marco de datos con los valores reales y estimados:
like <- sort(facebook$like)[1:489]
periodo <- rep(1:length(like), 2)
valores <- c(like, sort(modelo$fitted.values)[1:489])
origen <- c(rep("Reales", length(like)), rep("Estimados", length(like)))
facebook_gplot <- data.frame(periodo, valores, origen)Luego, utilicemos ggplot2 para crear el gráfico de líneas:
library(ggplot2)
ggplot(facebook_gplot, aes(periodo, valores, color = origen)) +
geom_point() +
geom_line() +
labs(x = "Periodo", y = "Likes", color = "Origen") +
theme_minimal()
Este gráfico muestra los valores reales de likes (en azul) y los valores estimados por el modelo (en rojo) a lo largo del tiempo (representado por el período). Permite una comparación visual de qué tan bien el modelo se ajusta a los datos reales.
Podemos ver que el modelo captura la tendencia general de los datos, pero hay discrepancias para publicaciones individuales. Esto es esperado, ya que el modelo es una simplificación de la realidad y no puede explicar toda la variabilidad en los datos.
En resumen, este gráfico proporciona una forma intuitiva de evaluar el rendimiento del modelo y entender sus limitaciones. Muestra que, mientras que el modelo de regresión lineal es una herramienta útil para comprender la relación entre likes, comentarios y compartidas, no es perfecto y debe usarse con precaución al hacer predicciones para publicaciones individuales.
Ejemplos
Ahora que tenemos un modelo de regresión lineal ajustado, podemos usarlo para predecir el número de likes para nuevas publicaciones de Facebook con diferentes características.
Ejemplo 1: Tenemos la página de Facebook de la selección méxicana de futbol con 11,000,000 likes totales, y estamos considerando hacer una publicación que esperamos genere 148 comentarios y 59 compartidas. Creemos un nuevo marco de datos con estos valores:
newpage <- data.frame(Page_total_likes = 11000000, comment = 148, share = 59)Ahora, usemos la función predict() para obtener una predicción del número de likes para esta publicación:
predict(modelo, newdata = newpage) 1
818.7032 El modelo predice que esta publicación recibirá aproximadamente 818 likes. Este es un número sustancial de likes, lo que sugiere que una publicación con estos niveles de comentarios y compartidas en una página con una gran base de seguidores podría tener un alto engagement.
Ejemplo 2: Ahora, consideremos una página de Facebook más pequeña como el de UNAM con 3,500,000 likes totales. Supongamos que estamos planeando una publicación que esperamos genere 1 comentarios y 1 compartidas. Nuevamente, creemos un marco de datos con estos valores:
newpage <- data.frame(Page_total_likes = 3500000, comment = 1, share = 1)Y usemos predict() para obtener una predicción:
predict(modelo, newdata = newpage) 1
8.698961 En este caso, el modelo predice aproximadamente 8.69 likes. Este es un número más bajo de likes en comparación con el primer ejemplo, lo que refleja el menor número de comentarios y compartidas, así como la base de seguidores más pequeña de la página.
Estos ejemplos ilustran cómo el modelo de regresión lineal puede usarse para obtener información valiosa sobre el potencial de engagement de las publicaciones de Facebook antes de que se publiquen. Al proporcionar valores hipotéticos para el número de comentarios y compartidas, así como el total de likes de la página, podemos obtener una estimación del número de likes que una publicación podría recibir.
Sin embargo, es crucial recordar que estas son solo predicciones basadas en un modelo simplificado. Hay muchos otros factores que podrían influir en el engagement real de una publicación, como el contenido de la publicación, la hora del día en que se publica, y las características demográficas de la audiencia de la página. Además, el modelo asume que la relación entre las variables es lineal, lo cual puede no ser siempre el caso.
A pesar de estas limitaciones, el modelo de regresión lineal sigue siendo una herramienta poderosa para comprender los impulsores del engagement en Facebook. Puede ayudar a los administradores de páginas a tomar decisiones informadas sobre su estrategia de contenido, al proporcionar una estimación cuantitativa del potencial de engagement de diferentes tipos de publicaciones.
En conclusión, estos ejemplos demuestran cómo el modelo de regresión lineal desarrollado en este análisis puede aplicarse a situaciones prácticas para predecir el número de likes en publicaciones de Facebook. Al combinar este modelo con otros conocimientos y experiencia en la gestión de redes sociales, los administradores de páginas pueden optimizar su estrategia de contenido para maximizar el engagement y alcanzar sus objetivos de marketing en Facebook.
Validación de supuestos
El modelo de regresión lineal es una buena herramienta de estimación; sin embargo, en dicho proceso se hace uso de diversos supuestos que deben cumplirse para que los resultados obtenidos sean acordes a la teoría desarrollada en los capítulos anteriores, a los que hasta el momento no se les había prestado gran atención en los ejercicios mostrados. Estos supuestos pueden verse en la definición para el caso de regresión simple y en la definición para regresión múltiple. A manera de resumen y de forma general, los supuestos que deben cumplirse son: la esperanza de los errores \(\varepsilon\) tiene media cero (\(\mathbb{E}[\varepsilon]=0\)), varianza constante sobre los errores (\(\mathbb{E}[\varepsilon]=\sigma^2\)), los errores no están correlacionados entre sí (\(\text{Cov}(\varepsilon_i, \varepsilon_j) \forall i\neq j\)) y, por último, los errores tienen distribución normal con media cero y varianza \(\sigma^2\), es decir, \(\varepsilon \sim \mathcal{N} (0,\sigma^2)\).
Es por ello que en este capítulo se analiza y verifica el cumplimiento de los anteriores supuestos. La mayoría de estas validaciones se basan en el análisis de residuales. Por lo tanto, se prestará atención a las propiedades y tipos de residuales.
Análisis de residuales
Anteriormente se han definido los residuales como \(e_i= y_i-\hat{y}_i\), cuyo nombre deriva de la obtención del residual o diferencia que existe entre la línea de regresión ajustada y los valores observados de la variable respuesta \(y_i\). Esta cantidad residual en un buen ajuste debería ser cercana a cero, pues cuando esto sucede se tiene que \(y_i \approx \hat{y}_i\). Debido a ello, se opta por trabajar con los residuales para verificar el cumplimiento de los supuestos que se detallarán a continuación.
Muchas pruebas usan como base modificaciones o tipos específicos de residuos, es por esta razón que se indagará sobre los diversos tipos más conocidos.
Residuales ordinarios
Los residuales ordinarios, o simplemente residuales, miden la diferencia entre la línea de regresión y la variable respuesta \(y_i\), y se definen como:
\[e_i= y_i-\hat{y}_i \qquad i=1,2, \ldots, n.\]
La desventaja que presentan estos residuales es que dependen de cada \(y_i\), por lo que observaciones atípicas pueden generar grandes residuales, ocasionando que pueda existir gran variabilidad al considerarse los errores en forma conjunta.
# Gráfica sencilla
e = modelo$residuals
plot(e)
# Crear un data frame con los residuales ordinarios
residuales_ordinarios <- data.frame(
Observacion = 1:length(e),
Residual = e
)
# Graficar los residuales ordinarios
ggplot(residuales_ordinarios, aes(x = Observacion, y = Residual)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(title = "Gráfico de Residuales Ordinarios",
x = "Observación",
y = "Residual") +
theme_minimal()
Residuales estandarizados
Una forma de disminuir esta variabilidad consiste en dividir a los errores \(e_i\) entre la varianza global, es decir:
\[d_i=\frac{e_i}{\sqrt{\sigma^2}}.\]
Una propiedad importante es que si \(e_i\) sigue una distribución normal, entonces al dividirlo entre la desviación estándar, se tiene que los residuales estandarizados siguen una distribución normal estándar.
# Calcular los residuales estandarizados
residuales_estandarizados <- rstandard(modelo)
# Crear un data frame con los residuales estandarizados
residuales_std <- data.frame(
Observacion = 1:length(residuales_estandarizados),
Residual = residuales_estandarizados
)
# Graficar los residuales estandarizados
ggplot(residuales_std, aes(x = Observacion, y = Residual)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(title = "Gráfico de Residuales Estandarizados",
x = "Observación",
y = "Residual Estandarizado") +
theme_minimal()
Residuales estudentizados
Los residuales estudentizados se basan en la idea de involucrar la varianza de cada observación en el cálculo de los residuales, ya que teóricamente los residuales no tienen varianza constante.
Como se demostró en el teorema de MCO se sabe que \(e=(I-H)\underline{Y}\). De esta manera, calculando la varianza de los errores se tiene el siguiente corolario:
Corolario 6: Sea \(\underline{e}\) los residuales del modelo, entonces la varianza de \(\underline{e}\) está dada por:
\[Var(\underline{e})=\sigma^2(I-H).\]
La varianza de cada residual no es constante para todas las observaciones, es por ello que el resultado depende de la siguiente forma:
\[Var(e_i)=\sigma^2(1-h_{ii}).\]
donde \(h_{ii}\) corresponde al i-ésimo elemento de la diagonal de la matriz sombrero \(H\).
Cada residual se divide entre su varianza, obteniendo de esta forma lo que se conoce como residuales estudentizados, los cuales se definen como:
\[r_i=\frac{e_i}{\sqrt{\hat{\sigma}^2 (1-h_{ii}) } }.\]
Debido a la construcción de los residuales estudentizados, se logra estandarizar de mejor manera a los residuales del modelo. Si se cumple que cada residual \(e_i\) sigue una distribución normal, entonces los residuales estudentizados siguen una distribución aproximada a una t-student con \(n-k-1\) grados de libertad.
# Calcular los residuales estudentizados
residuales_estudentizados <- rstudent(modelo)
# Crear un data frame con los residuales estudentizados
residuales_est <- data.frame(
Observacion = 1:length(residuales_estudentizados),
Residual = residuales_estudentizados
)
# Graficar los residuales estudentizados
ggplot(residuales_est, aes(x = Observacion, y = Residual)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(title = "Gráfico de Residuales Estudentizados",
x = "Observación",
y = "Residual Estudentizado") +
theme_minimal()
Si analizamos todos en conjunto se tiene
# Crear un data frame con los 3 tipos de residuales
residuales <- data.frame(
Observacion = 1:length(modelo$residuals),
Ordinarios = modelo$residuals,
Estandarizados = rstandard(modelo),
Estudentizados = rstudent(modelo)
)
# Cambiar el data frame a formato 'long'
residuales_long <- reshape2::melt(residuales, id.vars = "Observacion")
# Graficar los 3 tipos de residuales en un grid
ggplot(residuales_long, aes(x = Observacion, y = value)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
facet_wrap(~variable, ncol = 3, scales = "free_y") +
labs(title = "Comparación de Residuales",
x = "Observación",
y = "Residual") +
theme_minimal()
Supuesto de multicolinealidad
Se dice que el ajuste del modelo lineal sobre una muestra tiene problemas de multicolinealidad cuando hay una correlación alta entre dos o más variables explicativas, consecuencia de que una variable regresora es linealmente dependiente a alguna o algunas de las demás variables de la matriz diseño \(X\).
Un modelo de regresión lineal con esta característica provoca un error en la estimación de los parámetros de \(\underline{\beta}\), ya que al tener variables linealmente dependientes el determinante \(|X'X|\rightarrow 0\), lo cual provoca que la inversa \((X'X)^{-1}\rightarrow \infty\), ya que \((X'X)^{-1}=\frac{adj((X'X)')}{|X'X|}\). Esto provoca el cálculo erróneo de los parámetros, ya que el producto matricial de \(\hat{ \underline{\beta} }=(X'X)^{-1}X'\underline{Y}\) no pueda ser aproximado de la mejor manera, además de que numéricamente al programar los resultados se puede cometer el llamado underflow'' uoverflow’’ que son las cantidades extremas que pueden ser representadas computacionalmente.
Detección a través de correlaciones
El primer método para detectar la dependencia lineal es analizar el coeficiente de correlación de las variables, ya que si el modelo posee variables altamente correlacionadas entre si, es probable que el ajuste lineal presente problemas de multicolinealidad debido a relación estrecha de éstas variables, una desventaja de este método es que analiza a través de la matriz de correlaciones de una variable comparada con otra, sin embargo, puede existir una baja correlación entre variables analizándolas una a una pero al analizar de dos en dos, o en tres en tres, etcétera, puede encontrase una alta correlación.
Para analizar de mejor forma la información, se obtiene el determinante de la matriz de correlaciones, si este es muy cercano a cero podría proporcionar indicios de problemas de multicolinealidad, debido a que si hay dos variables que son linealmente dependientes el determinante de la matriz de correlaciones tendera a \(0\). Sin embargo, depende del prejuicio del investigador a partir de que cifra debe ser considerado una cantidad mínima.
library(corrplot)corrplot 0.92 loadedvars <- c("like", "comment", "share")
cor_mat <- cor(facebook2[, vars])
print(cor_mat) like comment share
like 1.0000000 0.8377624 0.9040278
comment 0.8377624 1.0000000 0.8683366
share 0.9040278 0.8683366 1.0000000corrplot(cor_mat, method = "color", type = "lower", order = "hclust", tl.coll = "black", addCoef.col = "black")Warning in text.default(pos.xlabel[, 1], pos.xlabel[, 2], newcolnames, srt =
tl.srt, : "tl.coll" is not a graphical parameterWarning in text.default(pos.ylabel[, 1], pos.ylabel[, 2], newrownames, col =
tl.col, : "tl.coll" is not a graphical parameterWarning in title(title, ...): "tl.coll" is not a graphical parameter
A través del Factor de inflación de la varianza
El factor de inflación de la varianza, usualmente abreviado como VIF, (por sus siglas en inglés, Variance inflation factor), se basa en el hecho de que cuando una matriz presenta problemas de multicolinealidad entonces el factor VIF presentará varianzas muy altas.
Se sabe que \(Var(\underline{\hat{ \beta }})=\sigma^2 I(X'X)^{-1}\) por lo que cada marginal de la matriz \(Var(\underline{\beta})\) se comporta de manera \(Var(\underline{\beta}_{j})=\sigma^2 I(X'X)_{jj}^{-1}\) debido a que \(\underline{\beta}\) se distribuye como una normal multivariada. Se puede observar de manera general que cada elemento de la diagonal de \((X'X)_{jj}^{-1}\) se comporta de la forma:
\[C_{jj}=\frac{1}{1-R_j^2}, \quad j=1, 2 \ldots ,k.\]
Donde \(R_j^2\) es el coeficiente de determinación que se obtiene al realizar la regresión (\(R^2=\frac{SC_{reg}}{SC_{TC}}\)), de esta forma, la anterior medida es de gran utilidad para detectar multicolinealidad, pues si \(X_j\) es independiente al resto de las variables regresoras, entonces \(R_j^2\) es cercano a 0 por lo que el valor de la diagonal \(C_{jj}\) se aproximaría a 1, indicando que no hay problema de multicolinealidad en esa variable. Por el contrario si \(X_j\) es dependiente al resto de las variables regresoras, entonces \(R_j^2\) es cercano a 1 por lo que el valor de la diagonal tendería a \(\infty\) teniendo evidencia de multicolinealidad.
De forma general se plantea que cuando el \(VIF\) sea menor o igual a diez (\(VIF\le 10\)), la variable que se analiza no presenta problemas serios de multicolinealidad, en caso contrario existe evidencia de correlación alta entre las variables.
library(car)Loading required package: carDatavif_valores <- vif(modelo)
vif_valores comment share
1.297975 1.297975 A través del índice \(\kappa\)
En 1960, la estadística Jacob Cohen, presentó el método del coeficiente kappa. El coeficiente kappa o simplemente \(\kappa\) se basa en el análisis de los eigenvalores \(\lambda_1, \lambda_2, \ldots, \lambda_k\) de la matriz \((X'X)\), este procedimiento es llamado analisis del eigensistema.
Debido a que \((X'X)\) es una matriz simétrica, el producto de la matriz diseño es igual a sus valores propios o eigenvalores. De esta manera el método propone observar la proporción del eigenvalor más grande respecto al más pequeño, si la proporción es pequeña, no hay problemas de multicolinealidad pues todos los valores propios son similares, sin embargo, proporciones altas indican gran variabilidad en los valores propios por lo que se tiene indicios de multicolinealidad. De esta manera el coeficiente \(\kappa\) es calculado como:
\[\kappa=\frac{\lambda_{max}}{\lambda_{min}}.\]
Cabe destacar que \(\lambda_{max}\) y \(\lambda_{min}\) es el valor propio \(\lambda\) más grande y más pequeño, respectivamente. El coeficiente de kappa tiene un rango permisible en el cual puede variar para ser clasificado de la siguiente manera:
- Si el coeficiente kappa de \((X'X) \le 100\) no se tiene problemas serios de multicolinealidad.
- Si Kappa de \((X'X) \in (100, 1,000)\) se tienen indicios de problemas serios de multicolinealidad.
- Si Kappa de \((X'X) \ge 1,000\) se tiene evidencia de un severo problema de multicolinealidad, ya que hay mucha variabilidad entre los valores. propios por lo que implica mayor inestabilidad al invertir la matriz \((X'X)\).
Este método proporciona mayor confiabilidad para determinar si se posee problemas graves de multicolinealidad pues en el anterior método del determinante de \((X'X)\) otorga una idea al problema de dependencia lineal entre las variables, sin embargo, el método presenta una mayor sensibilidad a los cambios de unidades. En el indice \(\kappa\) busca ser menos sensible al tomar los valores propios como medida de construcción de la estadística.
También se puede obtener el coeficiente \(\kappa\) de manera individual para cada variable explicativa \(j\) con \(j \in 1, \ldots, k\) para comprobar si se posee problemas o no. El procedimiento es similar al de manera global, con la diferencia de que se observa la proporción del valor propio con mayor valor respecto a cualquier otro eigenvalor. De esta manera el coeficiente \(\kappa\) de manera puntual es calculado como:
\[ \kappa_j=\frac{\lambda_{max}}{\lambda_j} \quad \forall j \in 1, \ldots, k. \]
Puede observarse que para toda \(j\) se cumple que \(\kappa\le \kappa_j\); al realizar este análisis de manera puntual la variabilidad del indice reduce de tal forma que:
- Si el coeficiente kappa en la j-ésima variable de \((X'X) \le 10\) no se tiene problemas serios de multicolinealidad.
- Si kappa en la j-ésima variable de \((X'X) \in (10, 30)\) se tienen indicios de problemas serios de multicolinealidad.
- Si kappa en la j-ésima variable de \((X'X) \ge 30\) se tiene evidencia de un severo problema de multicolinealidad.
Una consideración importante, es que algunos paquetes y específicamente en el software \(R\), los indices kappa son mostrados aplicando la raíz cuadrada. Es decir: \[\kappa=\sqrt{\frac{\lambda_{max}}{\lambda_{min}}} \qquad \qquad \kappa_j=\sqrt{\frac{\lambda_{max}}{\lambda_j}} \quad \forall j \in 1, \ldots, k.\]
Al considerar estos indices, entonces los criterios se reajustan obteniendo la raíz cuadrada de estos últimos.
X <- model.matrix(modelo)
XtX <- t(X) %*% X
eigenvalores <- eigen(XtX)$values
# Calcular el índice kappa global
kappa_global <- max(eigenvalores) / min(eigenvalores)
print(kappa_global)[1] 2543.121# Calcular el índice kappa por variable
kappa_vars <- numeric(ncol(X))
for (j in 1:ncol(X)) {
kappa_vars[j] <- max(eigenvalores) / eigenvalores[j]
}
print(kappa_vars)[1] 1.00000 16.10201 2543.12132Mejoras al problema de multicolinealidad
En caso de que se detecte problemas de multicolinealidad entre dos o más variables regresoras, el procedimiento más usual es eliminar del modelo, la variable o variables que presenten correlación muy alta, pues al estar muy correlacionadas entre si, la información que le aporta al modelo es muy parecida, así que una primera solución es eliminar alguna de aquellas variables y volver a realizar el ajuste mediante el nuevo conjunto de variables, con la finalidad de presentar un modelo sin problemas serios de multicolinealidad.
Un problema a la eliminación de variables, es elegir la variable que se desea quitar, puesto que en ocasiones son variables pausibles al modelo, es por ello que se prefiere realizar una variable ficticia, la cual es función de las variables que presentan correlaciones altas entre si, usualmente se realiza la suma, la multiplicación o la división de estas variables, es decir, supongamos que se tienen \(r\) variables que presentan problemas de multicolinealidad entre si, entonces se define a la variable \(X_\alpha\) como la variable ficticia función de las \(r\) variables, los cuales pueden ser expresadas de la forma:
\[X_\alpha=X_1+X_2+\ldots+X_r \qquad \mbox{ó} \qquad X_{\alpha}=X_1X_2\ldots X_r \qquad \mbox{ó} \qquad X_\alpha=\frac{X_1+X_2+\ldots}{\ldots+X_{r-1}+X_r} \qquad etc.\]
Una última mejora al problema de multicolinealidad, aunque más compleja, es el método de componentes principales el cual realiza una otorgonalización entre las variables más correlacionadas obteniendo un vector como nueva variable función de las variables con problemas de multicolinealidad. El lector interesado en el tema puede consultar.
modelo2_share <- lm(like ~ share - 1, data = facebook)
summary(modelo2_share)
Call:
lm(formula = like ~ share - 1, data = facebook)
Residuals:
Min 1Q Median 3Q Max
-389.25 -38.30 -7.14 18.79 1086.79
Coefficients:
Estimate Std. Error t value Pr(>|t|)
share 6.068 0.157 38.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 111.6 on 488 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.7538, Adjusted R-squared: 0.7532
F-statistic: 1494 on 1 and 488 DF, p-value: < 2.2e-16Validación de Supuestos
En el modelo de regresión lineal, es importante verificar que se cumplan ciertos supuestos para garantizar la validez de los resultados obtenidos. A continuación, se describen los principales supuestos y cómo evaluarlos.
Media igual a 0
El supuesto de media igual a cero establece que la media de los residuales debe ser aproximadamente cero. Esto implica que no hay un sesgo sistemático en los residuales. Puedes verificar este supuesto calculando la media de los residuales y comprobando que se aproxime a cero.
mean(residuales$Ordinarios)[1] 2.734154Si la media de los residuales se desvía significativamente de cero, puede indicar que el modelo no está capturando adecuadamente la relación entre las variables.
Varianza Constante (Homocedasticidad)
Se dice que una muestra es homocedástica cuando la varianza es constate a lo largo de todas las observaciones, es decir, no varia conforme se presentan nuevas observaciones. Mientras una muestra heterocedástica se presenta cuando hay variaciones de la varianza conforme se presentan nuevas observaciones.
La desviaciones en el supuesto de homocedasticidad pueden observarse mediante gráficas, la más optima para el análisis es realizar una gráfica de dispersión en el que se muestre la relación entre los valores ajustados \(\hat{ \underline{Y} }\) contra los residuales estandarizados \(d_i\). Si la varianza es constante entonces la gráfica fluctuaría entre el eje horizontal de manera simétrica, asemejando a una distribución uniforme, y sin seguir algún tipo de patrón, ya que típicamente se considera que la mayor parte de los errores deben estar contenidos en franjas horizontales delimitados por el eje vertical entre \(y=-2\) y \(y=2\).
El desarrollo de las gráficas mencionadas, puede ser analizada mediante la gráfica el cual fue replicado mediante el siguiente código de R.
n <- 1000 # tamaño de muestra
set.seed(2019)
x <- rnorm(76, 20, 6) # Se generan las X
y1 <- 28 + 5 * x + rnorm(76, 0, 1) # Se generan las Y con errores homocedásticos
y2 <- 28 + 5 * x + 100*x * rnorm(76, 0, 30) # Se generan las Y con errores heterocedásticos
y3 <- 28 + 5 * x + sin(sort(pi*x)) * rnorm(76, 0, 30) # Se generan las Y con errores heterocedásticos
y4 <- 28 + 5 * x + 1500 * rnorm(76, 0, 30) # Se generan las Y con errores heterocedásticos
modelo1 <- lm(y1 ~ x) # Modelo con errores homocedásticos
summary(modelo1)
Call:
lm(formula = y1 ~ x)
Residuals:
Min 1Q Median 3Q Max
-2.60412 -0.54142 -0.09057 0.70900 2.17184
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.56200 0.37660 75.84 <2e-16 ***
x 4.96474 0.01825 272.10 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9321 on 74 degrees of freedom
Multiple R-squared: 0.999, Adjusted R-squared: 0.999
F-statistic: 7.404e+04 on 1 and 74 DF, p-value: < 2.2e-16modelo2 <- lm(y2 ~ x) # Modelo con errores heterocedásticos
summary(modelo2)
Call:
lm(formula = y2 ~ x)
Residuals:
Min 1Q Median 3Q Max
-157569 -34333 -1719 39704 106435
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2050.2 22620.0 0.091 0.928
x -766.4 1095.9 -0.699 0.487
Residual standard error: 55980 on 74 degrees of freedom
Multiple R-squared: 0.006564, Adjusted R-squared: -0.00686
F-statistic: 0.489 on 1 and 74 DF, p-value: 0.4866modelo3 <- lm (y3~x)
summary(modelo3)
Call:
lm(formula = y3 ~ x)
Residuals:
Min 1Q Median 3Q Max
-45.646 -10.190 0.683 9.231 66.298
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.8346 7.8412 3.295 0.00151 **
x 5.0524 0.3799 13.299 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.41 on 74 degrees of freedom
Multiple R-squared: 0.705, Adjusted R-squared: 0.701
F-statistic: 176.9 on 1 and 74 DF, p-value: < 2.2e-16modelo4 <- lm (y4~x)
summary(modelo4)
Call:
lm(formula = y4 ~ x)
Residuals:
Min 1Q Median 3Q Max
-119375 -24054 -3232 34966 95062
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15036.9 18717.8 0.803 0.424
x -1097.4 906.9 -1.210 0.230
Residual standard error: 46330 on 74 degrees of freedom
Multiple R-squared: 0.01941, Adjusted R-squared: 0.006154
F-statistic: 1.464 on 1 and 74 DF, p-value: 0.2301ajustados1 <- fitted(modelo1)
ajustados2 <- fitted(modelo2)
ajustados3 <- fitted(modelo3)
ajustados4 <- fitted(modelo4)
library(MASS)
residuos1 <- stdres(modelo1)
residuos2 <- stdres(modelo2)
residuos3 <- stdres(modelo3)
residuos4 <- stdres(modelo4)
par(mfrow = c(2, 2))
plot(ajustados1, residuos1, pch = 16, col = 'darkcyan',
xlab = 'V. Ajustados', ylab = 'residuales', main = 'Errores homocedásticos')
lines(loess.smooth(ajustados1, residuos1,span=1), col = 'red3', lwd = 3, lty = 3)
plot(ajustados2, residuos2, pch = 16, col = 'darkcyan',
xlab = 'V. Ajustados', ylab = 'residuales', main = 'Errores heterocedásticos')
lines(loess.smooth(ajustados2, residuos2, span=1), col = 'red3', lwd = 3, lty = 3)
plot(ajustados3, residuos3, pch = 16, col = 'darkcyan',
xlab = 'V. Ajustados', ylab = 'residuales', main = 'Errores heterocedásticos')
lines(loess.smooth(ajustados3, residuos3, span=1), col = 'red3', lwd = 3, lty = 3)
plot(ajustados4, residuos4, pch = 16, col = 'darkcyan',
xlab = 'V. Ajustados', ylab = 'residuales', main = 'Errores heterocedásticos')
lines(loess.smooth(ajustados4, residuos4, span=1), col = 'red3', lwd = 3, lty = 3)
Si la dispersión de los residuales no es constante, sino que aumenta o disminuye sistemáticamente a lo largo de los valores ajustados, se viola el supuesto de homocedasticidad. En este caso, puede ser necesario aplicar transformaciones a las variables o utilizar modelos alternativos.
Prueba de Breusch-Pagan
La prueba de Breusch-Pagan es una prueba estadística utilizada para detectar la presencia de heterocedasticidad (varianza no constante) en un modelo de regresión lineal. La prueba contrasta la hipótesis nula de homocedasticidad (varianza constante) contra la hipótesis alternativa de heterocedasticidad.
La prueba se basa en calcular los residuales estandarizados al cuadrado \(\left(\tilde{e}_j^2=\frac{e_j^2}{\sigma^2}\right)\) y realizar una regresión lineal auxiliar, donde la variable dependiente son los residuales estandarizados al cuadrado, y las variables independientes son las variables explicativas del modelo original.
La estadística de prueba es \(\frac{SC_{reg}}{2}\), donde \(SC_{reg}\) es la suma de cuadrados de la regresión auxiliar. Bajo la hipótesis nula de homocedasticidad, esta estadística sigue una distribución chi-cuadrada con \(k\) grados de libertad, siendo \(k\) el número de variables explicativas.
En R, la prueba de Breusch-Pagan se puede implementar utilizando la función bptest() del paquete lmtest:
library(lmtest)Loading required package: zoo
Attaching package: 'zoo'The following objects are masked from 'package:base':
as.Date, as.Date.numeric# Realizar la prueba de Breusch-Pagan
bptest(modelo)
studentized Breusch-Pagan test
data: modelo
BP = 32.756, df = 2, p-value = 7.712e-08La salida proporcionará la estadística de prueba, los grados de libertad y el valor p. Si el valor p es menor que el nivel de significancia establecido (por ejemplo, 0.05), se rechaza la hipótesis nula de homocedasticidad, indicando la presencia de heterocedasticidad en el modelo.
Prueba de White
La prueba de White es otra prueba estadística utilizada para detectar la presencia de heterocedasticidad en un modelo de regresión lineal. A diferencia de la prueba de Breusch-Pagan, la prueba de White no asume que los errores sigan una distribución normal.
La prueba se basa en calcular los residuales estandarizados \(\tilde{e}_j=\frac{e_j^2}{\sigma^2}\) y realizar una regresión auxiliar, donde la variable dependiente son los residuales estandarizados al cuadrado, y las variables independientes son las variables explicativas del modelo original, sus cuadrados y sus productos cruzados.
La estadística de prueba es \(nR^2\), donde \(n\) es el tamaño de la muestra y \(R^2\) es el coeficiente de determinación de la regresión auxiliar. Bajo la hipótesis nula de homocedasticidad, esta estadística sigue una distribución chi-cuadrada con grados de libertad iguales al número de variables explicativas en la regresión auxiliar.
En R, la prueba de White se puede implementar utilizando la función whites.htest() del paquete het.test:
library(whitestrap)
Please cite as: Lopez, J. (2020), White's test and Bootstrapped White's test under the methodology of Jeong, J., Lee, K. (1999) package version 0.0.1# Realizar la prueba de White
white_test(modelo2_share)White's test results
Null hypothesis: Homoskedasticity of the residuals
Alternative hypothesis: Heteroskedasticity of the residuals
Test Statistic: 31.94
P-value: 0La salida proporcionará la estadística de prueba, los grados de libertad y el valor p. Si el valor p es menor que el nivel de significancia establecido, se rechaza la hipótesis nula de homocedasticidad, indicando la presencia de heterocedasticidad en el modelo.
Tanto la prueba de Breusch-Pagan como la prueba de White son herramientas útiles para detectar la presencia de heterocedasticidad en un modelo de regresión lineal. La elección entre estas dos pruebas puede depender de los supuestos subyacentes y las características específicas del conjunto de datos.
Normalidad
El supuesto de normalidad establece que los residuales deben seguir una distribución normal. Esto es importante para la validez de los intervalos de confianza y las pruebas de hipótesis basadas en el modelo.
Puedes evaluar visualmente este supuesto mediante un gráfico de cuantiles normales (Q-Q plot) de los residuales.
qqnorm(residuales$Ordinarios)
qqline(residuales$Ordinarios)
Si los puntos en el gráfico Q-Q se desvían considerablemente de la línea recta, puede indicar que los residuales no siguen una distribución normal. En este caso, puede ser necesario aplicar transformaciones a las variables o utilizar modelos no paramétricos.
Es importante tener en cuenta que los supuestos mencionados son especialmente relevantes cuando se realizan inferencias basadas en el modelo de regresión lineal. Si los supuestos no se cumplen, los resultados y las conclusiones del modelo pueden ser cuestionables.
Además de evaluar visualmente la normalidad de los residuales mediante un gráfico Q-Q, es recomendable complementar el análisis con una prueba estadística formal, como la prueba de Anderson-Darling.
La prueba de Anderson-Darling es una prueba no paramétrica que se utiliza para verificar si una muestra de datos proviene de una distribución específica, en este caso, la distribución normal. Esta prueba es una variante de la prueba de Kolmogorov-Smirnov y es más sensible a las desviaciones en las colas de la distribución.
La prueba de Anderson-Darling contrasta la hipótesis nula de que los residuales siguen una distribución normal, frente a la hipótesis alternativa de que no siguen una distribución normal.
En R, puedes realizar la prueba de Anderson-Darling utilizando la función ad.test() del paquete nortest:
library(nortest)
# Obtener los residuales del modelo
residuales <- resid(modelo2_share)
# Realizar la prueba de Anderson-Darling
resultado_ad <- ad.test(residuales)
# Imprimir los resultados
print(resultado_ad)
Anderson-Darling normality test
data: residuales
A = 35.701, p-value < 2.2e-16En este caso, se rechaza la prueba de hipotesis, así que los datos no se distribuyen con normalidad
Supuesto de Linealidad
En la construcción del modelo de regresión lineal, se asume que la relación entre cada variable independiente \(X_j\) y la variable dependiente \(Y\) es lineal, para \(j = 1, \ldots, k\), donde \(k\) es el número de variables regresoras con las que se ajustó el modelo.
Sin embargo, esta afirmación no siempre se cumple, por lo que es necesario validar el supuesto de linealidad de manera gráfica. Debido a la dificultad de representar gráficamente más de tres dimensiones, se recuerda que en el modelo de regresión múltiple con \(k\) regresores, se ajusta un hiperplano de dimensión \(k\). Cuando \(k \geq 4\), se recomienda realizar gráficas individuales para comprobar la linealidad entre cada variable explicativa \(X_j\) y la variable de interés \(Y\).
Una forma de evaluar la linealidad es graficar cada variable independiente \(X_j\) contra la variable dependiente \(Y\). Si la relación es lineal, se esperaría observar una tendencia lineal en la gráfica. Sin embargo, este método puede proporcionar conclusiones erróneas cuando los coeficientes tienen magnitudes distintas, ya que se analiza la relación marginal de la variable respuesta con cada variable explicativa.
Por lo tanto, es preferible realizar un análisis de residuales. En este análisis, se grafican los residuales estandarizados contra los valores observados de cada variable explicativa \(X_j\). Si se cumple el supuesto de linealidad, se esperaría observar un patrón de ruido blanco con media cero y varianza \(\sigma^2\).
# Cargar conjunto de datos mtcars
data(mtcars)
# Ajustar el modelo de regresión lineal
modelo <- lm(mpg ~ cyl + disp + hp, data = mtcars)
# Obtener los residuales estandarizados
residuales_estandarizados <- rstandard(modelo)
# Crear una función para graficar los residuales contra cada variable independiente
graficar_residuales <- function(variable) {
plot(mtcars[, variable], residuales_estandarizados, xlab = names(mtcars)[variable],
ylab = "Residuales Estandarizados", main = "Análisis de Linealidad")
lines(loess.smooth(mtcars[, variable], residuales_estandarizados, span = 0.75),
col = "red", lwd = 2)
}
# Graficar los residuales contra cada variable independiente
par(mfrow = c(1, 3))
graficar_residuales("cyl")Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
pseudoinverse used at 6Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
neighborhood radius 2Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
reciprocal condition number -0Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
pseudoinverse used at 6Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
neighborhood radius 2Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
reciprocal condition number -0Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
pseudoinverse used at 6Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
neighborhood radius 2Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
reciprocal condition number -0Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
pseudoinverse used at 6Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
neighborhood radius 2Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
reciprocal condition number -0Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
pseudoinverse used at 6Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
neighborhood radius 2Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
reciprocal condition number -0graficar_residuales("disp")
graficar_residuales("hp")
Si se detectan problemas de linealidad entre las variables explicativas y la variable de interés, el ajuste del modelo puede ser deficiente, ya que la varianza puede presentar problemas en su estimación, lo que afectará a los estadísticos que utilicen \(\sigma^2\) en su desarrollo.
Una posible solución al problema de linealidad es agregar o eliminar variables que causen problemas de linealidad. Otra solución es ajustar una regresión polinomial, generalmente como resultado de la prueba de falta de ajuste (lack of fit). Esta prueba consiste en agregar al modelo cada variable explicativa elevada al cuadrado de manera independiente y observar si la prueba de hipótesis proporciona evidencia para rechazar o no la variable explicativa propuesta mediante la prueba \(\gamma = 0\), donde \(\gamma\) es el coeficiente de la variable explicativa \(j\) al cuadrado. Es decir, se propone el modelo:
\[Y = \hat{\beta}_0 + \hat{\beta}_1X_1 + \hat{\beta}_2X_2 + \ldots + \hat{\beta}_kX_k + \varepsilon + \hat{\gamma}X_j^2\]
Si hay evidencia de que \(H: \gamma \neq 0\), entonces se elimina la variable \(X_j\) y se realiza el ajuste mediante el cuadrado \(X_j^2\), ya que aporta información y se obtiene una posible linealidad en el modelo.
Es importante validar el supuesto de linealidad, ya que su incumplimiento puede afectar la validez de las inferencias realizadas a partir del modelo de regresión lineal.
Valores Atípicos (Outliers)
Los valores atípicos, también conocidos como outliers, son observaciones de la muestra aleatoria que no se comportan de manera similar al resto de los elementos que conforman el conjunto de datos. Gráficamente, una observación con un valor atípico no sigue la tendencia general que presenta la muestra aleatoria.
Una forma de detectar posibles valores atípicos es mediante un análisis de residuales. Este análisis consiste en obtener los residuales, ya sean estandarizados o estudentizados, y observar si estos son significativamente mayores o menores en comparación con un punto crítico con un nivel de significancia \(\alpha\).
Si se trabaja con residuales estandarizados, se recuerda que si \(e_i\) sigue una distribución normal, los residuales estandarizados siguen una distribución normal con media 0 y varianza \(1\). Por lo tanto, los residuales superiores o inferiores al punto crítico \(\pm Z_{1-\alpha/2}\) se consideran posibles outliers con un nivel de significancia \(\alpha\).
Por otro lado, si se trabaja con residuales estudentizados, el punto crítico está determinado por los residuales que se encuentren por encima o por debajo de la banda determinada por el cuantil \(\pm t_{1-\alpha/2, n-k-1}\). En otras palabras, hay evidencia de un valor atípico con un nivel de significancia \(\alpha\) cuando se cumple alguna de las siguientes desigualdades:
\[|d_i| \geq Z_{1-\alpha/2}\] \[|r_i| \geq t_{1-\alpha/2, n-k-1}\]
donde \(d_i\) son los residuales estandarizados, \(r_i\) son los residuales estudentizados, \(n\) es el tamaño de la muestra y \(k\) es el número de variables independientes.
Aquí hay un ejemplo de código en R para detectar valores atípicos utilizando residuales estandarizados con un nivel de confianza del 99%:
data(mtcars)
# Ajustar el modelo de regresión lineal
modelo <- lm(mpg ~ cyl + disp + hp, data = mtcars)
# Detección de valores atípicos
# Obtener el punto crítico (cuantil de la distribución normal)
z_critico <- qnorm(1 - 0.01/2)
# Obtener los residuales estandarizados
residuales <- rstandard(modelo)
# Identificar valores atípicos
outliers <- which(abs(residuales) > z_critico)
# Imprimir los índices de los valores atípicos
print(outliers)named integer(0)boxplot(residuales)
En este código, ajustamos un modelo de regresión lineal, obtenemos el punto crítico z_critico de la distribución normal para un nivel de confianza del 99%, calculamos los residuales estandarizados y luego identificamos los índices de las observaciones cuyos residuales estandarizados (en valor absoluto) exceden el punto crítico z_critico.
Es importante realizar un análisis de influencia de las observaciones que presentan problemas de outliers, ya que aunque se traten de puntos atípicos, pueden ser significativos para el modelo. Eliminar estas observaciones sin un análisis adecuado puede ocasionar desviaciones en la estimación.
Valores Influyentes
Los valores influyentes son observaciones que tienen un gran impacto en el ajuste del modelo de regresión lineal. Es decir, si se eliminan estas observaciones, podría provocar un cambio significativo en los estimadores de los parámetros o en las predicciones del modelo.
Por esta razón, es importante analizar estos puntos para medir su influencia en el modelo y determinar si deben ser eliminados o no de la muestra.
Un método para identificar valores influyentes es a través de los puntos de apalancamiento (leverage points). Este método consiste en examinar la distancia entre cada punto y el punto medio (centroide) de los datos. Para ello, se analiza la matriz sombrero \(H\), y se presta especial atención a los elementos de la diagonal principal, denotados como \(h_{ii}\), conocidos como puntos de apalancamiento.
Se considera que una observación tiene un punto de apalancamiento grande cuando:
\[h_{ii} \geq \frac{\sum_{i=1}^n h_{ii}}{n}\]
Es decir, cuando el punto de apalancamiento de una observación es mayor que el doble del promedio de los puntos de apalancamiento. Sin embargo, esto no es suficiente para afirmar que la observación es influyente, por lo que se requiere el uso de otros métodos estadísticos, como la distancia de Cook.
La distancia de Cook mide la distancia cuadrática entre el modelo ajustado y el modelo ajustado sin la observación \(i\)-ésima. Se calcula de la siguiente manera:
\[C_i = \frac{\left(\hat{\underline{Y}} - \hat{\underline{Y}}_{(i)}\right)'\left(\hat{\underline{Y}} - \hat{\underline{Y}}_{(i)}\right)}{CM_{error}\sum_{i=1}^n h_{ii}}, \quad \forall i \in [1, n]\]
Donde \(\hat{\underline{Y}}\) es el vector de valores ajustados del modelo completo, \(\hat{\underline{Y}}_{(i)}\) es el vector de valores ajustados del modelo sin la observación \(i\)-ésima, y \(CM_{error}\) es la cuadrado medio del error.
Se considera que una observación es influyente si su distancia de Cook es mayor o igual a 1, según Hair et al. (1998), o mayor que \(\frac{4}{n}\), según Bollen y Jackman (1985).
Aquí hay un ejemplo de código en R para calcular la distancia de Cook:
# Detección de valores influyentes
# Calcular la distancia de Cook
cook_dist <- cooks.distance(modelo)
# Identificar valores influyentes
influyentes <- which(cook_dist >= 1)
# Imprimir los índices de los valores influyentes
print(influyentes)named integer(0)En este código, ajustamos un modelo de regresión lineal, calculamos la distancia de Cook utilizando la función cooks.distance(), e identificamos los índices de las observaciones cuya distancia de Cook es mayor o igual a 1.
Es importante analizar cuidadosamente los valores influyentes antes de decidir si se deben eliminar o no de la muestra, ya que, aunque sean observaciones atípicas, pueden ser significativas para el modelo y su eliminación podría introducir sesgos en la estimación.
plot(modelo)
