11 Correlación y regresión lineal

Autor/a

Esta sección está en fase de preparación

11.1 Medidas de influencia

11.2 Leverage o apalancamiento

En un modelo de regresión lineal simple, el estimador del coeficiente de regresión es la media ponderada de las pendientes de las rectas que unen el centro de gravedad de la distribución $(\bar{x}, \bar{y})$ con cada uno de los puntos $(x_i, y_i)$ (Peña, 2002). El factor de ponderación es en cada caso \[w_i=\frac{(x_i-\bar{x})^2}{\sum_{j=1}^{n} {(x_j-\bar{x})^2}}\] observando que $0\le w_i \le (n-1)/n$. De aquí se define el leverage o efecto palanca de cada observación como
\[h_i = \frac{1}{n}+w_i \] valor que estará comprendido en el intervalo $(1/n, 1)$. La suma de todos los $h_i$ es

\[\sum{h_i} = 1 + \frac{\sum_{i=1}^{n} {(x_j-\bar{x})^2}}{\sum_{j=1}^{n} {(x_j-\bar{x})^2}} = 2,\]

de manera que el valor promedio es $\bar{h}=2/n$.

Se considera que una observación tiene un efecto palanca importante si $h_i\ge 3\bar{h} =6/n$, aunque A menudo se establece como umbral $2p/n$ con $p$ el número de parámetros del modelo. Sin embargo, un efecto palanca importante no indica que la observación sea necesariamente perjudicial para el modelo (el valor residual puede ser bajo). La influencia real se establece con las distancias de Cook, que combinan la información dada por el valor residual con la del efecto palanca.

Obtención en R: hatvalues(modelo). En la literatura en inglés, es habitual encontrar la denominación hat value para aludir a este indicador del efecto palanca, de aquí el nombre de la función.

11.3 Distancias de Cook

La distancia de Cook es una medida que cuantifica el impacto global que tiene eliminar una observación sobre todos los coeficientes del modelo. Mientras que los puntos palanca indican la potencialidad de ser influyentes la distancia de Cook mide la influencia real de los puntos sobre los coeficientes del modelo. \[D_i=\frac{t_i^2}{k s^2_R h_i}\]

con $t_i=(\hat{y}_i-\hat{y}_{(i)}$ el residuo estudentizado, que es la diferencia entre las estimaciones de $y_i$ a partir del modelo y a partir de un modelo alternativo en el que se ha excluido esa observación. $k$ es el número de parámetros del modelo (incluida la constante) y $h_i$ es el leverage.

Se considera que un punto es influyente si $D_i>1$, aunque, en la práctica, se suele considerar el umbral $4/(n-p)$

Obtención en R: cooks.distance(modelo)

11.4 DFFITS

Es un indicador que mide el cambio del pronóstico para la observación $i$ cuando se elimina dicha observación. Al igual que la distancia de Cook, agrega la información de los residuos estudentizados con la del leverage, pero ahora se trata de evaluar la influencia en el pronóstico individual.

\[DFFITS= t_i \sqrt{\frac{h_i}{1-h_i}}\]

con $t_i=(\hat{y}_i-\hat{y}_{(i)}$ el residuo estudentizado, y $h_i$ el leverage de la observación $i$

El umbral habitual es $2 \sqrt{(p/n)}$

Obtención en R: dffits(modelo)

11.5 DFBETAS

Para cada coeficiente del modelo $\beta_j$ se define un indicador $DFBETAS_{j(i)}$ que indica cuánto cambia dicho coeficiente cuando se elimina la observación $i$, lo que permite detectar el efecto de dicha observación sobre coeficientes específicos. La magnitud de los $DFBETAS$ está estandarizada, de manera que expresan el cambio en unidades de desviación estándar. \[ DFBETAS_{ji} = \frac{\hat{\beta}_j-\hat{\beta}_{j(i)}}{SE(\hat\beta_{j(i)})} \]

siendo $\hat\beta_{j(i)}$ la estimación del coeficiente $\beta_j$ eliminando la observación $i$ y $SE(\hat\beta_{j(i)})$ su error estándar.

Como criterio de evaluación, con muestras pequeñas se puede considerar influyente un valor con $|DFBETAS|_{ji}>1$ mientras que para muestras grandes suele considerarse el umbral $|DFBETAS|_{ji}>2/\sqrt{n}$

Obtención en R: dfbetas(modelo)

11.6 COVRATIO

El indicador covratio evalúa cuánto cambia el volumen de la elipse de confianza para los parámetros del modelo cuando se elimina una observación. Es decir, cuantifica el cambio en la precisión de las estimaciones de los coeficientes. Se obtiene como

\[\mathrm{COVRATIO}_i = \frac{\det(\widehat{\mathrm{Cov}}_{(i)})}{\det(\widehat{\mathrm{Cov}})}\]

siendo $\widehat{Cov}$ el estimador de la matriz de covarianzas de los coeficientes, con todas las observaciones en el numerador y eliminando la observación $i$-ésima en el denominador. Si $CVR \approx 1$ el caso no altera significativamente la matriz de covarianza. El umbral frecuente para este indicador es $1+(3p/n)$.

Obtención en R: covratio(modelo)

11.7 Índice de Hadi

Es una medida que combina información sobre la distancia de Mahalanobis en el espacio de las variables explicativas y el residuo estandarizado. Su objetivo es identificar observaciones que son atípicas en el predictor, en la respuesta o en ambas. La expresión general simplificada para la observación $i$ es:

\[H_i=\frac{h_{ii}}{1-h_{ii}}+\frac{e_i^2}{\hat\sigma^2_\epsilon (1-h_{ii})}\]

con $e_i=y_i-\hat{y_i}$ el residuo crudo ordinario, $h_i$ el leverage y $\sigma^2_\epsilon$ la varianza residual.

La magnitud de la medida está condicionada por el tamaño muestral y por el número de parámetros del modelo. Una regla práctica es considerar como umbral el percentil 95 de la distribución de los $H_i$

11.8 Obtención y evaluación de medidas de influencia con R de manera simultánea

En R, la función influence.measures(modelo) genera una lista con dos objetos de tipo data.frame:

$influence.measures(modelo)$infmat contiene los valores de las medidas de influencia presentadas aquí.
influence.measures(modelo)$is.inf es un data.frame de valores lógicos (TRUE/FALSE) indicadores de si el caso al que alude cada fila ha superado el umbral de la medida presentada en cada columna.

El objeto creado con influence.measures() implementa una función summary() que proporciona un resumen de las medidas de influencia.

11.9 Resumen de funciones

11.10 Notación general

En las expresiones que se presentan a continuación, los términos que aparecen en cursiva aluden a qué tipo de término debe de aparecer en dicha expresión pero su nombre debe sustituirse por el que tenga sentido en el contexto del código. Las funciones, cuyo nombre se debe escribir de forma literal, se indicaran con un sombreado, por ejemplo plot(). En ellas, los argumentos que no estén en cursiva o que estén también sombreados son literales, se deben escribir tal y como se presentan. En general, se hará uso de los siguientes elementos

datos representa el nombre del data.frame con los datos
y es la variable explicada por el modelo (de forma natural, será una columna de datos)
x es la variable que actúa como explicativa (igual que y será una columna de datos)
x₁,…, x_p alude al conjunto de variables explicativas (columnas en el data.frame datos)
modelo alude al nombre del modelo resultante de una expresión. Si es necesario distinguir varios modelos se usarán sufijos, por ejemplo modeloLoess.

Cuando una función no este definida en el lenguaje base de R, se utilizará la sintaxis que hace explícito el paquete en el que está incluida. Por ejemplo, car::scatterplot() quiere decir que la función scatterplot() está definida en el paquete car. Esta sintaxis es totalmente funcional y, si el paquete en cuestión está instalado, evita cargar el paquete entero a través de library(). Su uso resulta muy conveniente para facilitar el seguimiento del código.

11.11 Regresión lineal simple

11.11.1 El diagrama de dispersión

Son equivalentes las expresiones

plot(x, y, data=datos)
plot(y~x, data=datos)

El paquete car presenta la función scatterplot() que incluye diagramas de caja para representar la distribución marginal de cada variable

car::scatterplot(y~x, data=datos, main=“título”)

11.11.2 Modelo de regresión lineal simple

modelo <- lm(y~x, data=datos)

La función lm() (linear model) crea un objeto de clase lm. El primer argumento de esta función es un objeto de clase formula. El símbolo ~ que debe de aparecer en ella es una función que la define como tal. La expresión que se indique a la izquierda de ~ ocupará el rol de variable respuesta y la expresión que aparezca a la derecha será la variable explicativa o regresora. Estos términos de la fórmula no están restringidos a nombres de variable, pueden ser cualquier expresión válida de R, que al ser interpretada actúa como una variable. Por ejemplo

log(y)~x

permite modelizar el logaritmo de y como una función lineal de x.

Algunos modelos con un predictor

Manually changed
Modelo	Fórmula	Observación
Nulo	`y ~ 1`	1 es la intersección en los modelos de regresión lineal. Aquí representa a la media $\bar{y}$
Regresión lineal simple	`y ~ x`	x continua
Regresión por el origen	`y ~ x-1`	Se elimina la constante del modelo
ANOVA (1 vía)	`y ~ A`	A es un factor con k niveles
Polinomio cuadrático	`y ~ I(x^2)` `y ~ poly(x,2)`	las dos expresiones son equivalentes. I() fuerza a la evaluación del argumento.

Componentes de un objeto de clase lm

Los componentes más relevantes son:

modelo$coefficients - vector de coeficientes
modelo$residuals - diferencia entre los valores observados y pronosticados de y para los valores observados de x
modelo$fitted.values - valores pronosticados de y para los valores observados de x
modelo$model - permite acceder a los datos usados para el ajuste (por ejemplo modelo $model$ x)

Funciones aplicables al objeto lm

Un objeto de clase lm() puede constituir el argumento de un buen número de funciones

Listado de funciones

summary(modelo) - Estimación de los coeficientes, su error estándar y el test individual de significación.
anova(modelo) -
aov(modelo) -
coefficients(modelo) - Valores de los coeficientes del modelo
confint(modelo) - Intervalos de confianza para los coeficientes
residuals(modelo) - Proporciona los residuos del modelo
abline(modelo) - Incorpora la recta de regresión a un diagrama de dispersión generado por plot()
predict() - Proporciona los pronósticos del modelo para los valores deseados del regresor
plot(modelo) - Proporciona hasta seis diagramas para la validación del modelo
hatvalues(modelo)

Significación de los coeficientes del modelo lm

summary(modelo)

Descomposición de la varianza

anova(modelo) -
aov(modelo) -

Residuos

residuals(modelo) -

Medidas de influencia

influence.measures(modelo)

el objeto influence.measures() presenta dos campos

influence.measures(modelo)$infmat - es una matriz con los valores de las medidas de influencia
influence.measures(modelo)$is.inf - es una matriz imagen de la anterior pero compuesta por valores lógicos (TRUE/FALSE) que indican si cada valor de cada medida supera o no el umbral.

es útil summary(influence.measures(modelo))

Se puede acceder a cada una de las medidas de influencia de forma individual - hatvalues(modelo)

cooks.distance(modelo)
dfbetamodelo) - sin estandarizar
dfbetas(modelo) - estandarizadas
dffits(modelo) -
covratio(modelo) -

Diagramas de influencia de otros paquetes

car::influencePlot(modelo)
olsrr::ols_plot_cooksd_bar(modelo)
olsrr::ols_plot_cooksd_chart(modelo)

Actualización del modelo

La función update() permite actualizar el modelo para que se ajuste a un subconjunto de datos. Dicho subconjunto se puede definir como un filtro, por ejemplo, el pseudocódigo

filtro<-((x>valor₁) & (x<valor₂))

modelo2 <- update(modelo,subset=filtro)

adaptará el ajuste del modelo a los ****

Pronósticos con el modelo lm()

Pronósticos para el valor medio de la respuesta dada una serie de valores del regresor (los valores de x a pronosticar deben tener la estructura de data.frame)

predict(modelo,newdata=data.frame(x=c(valor₁,…,valor_k)),interval="confidence")

Pronósticos para nuevas observaciones de la respuesta dada una serie de valores del regresor (como antes, los valores de x a pronosticar deben tener la estructura de data.frame)

predict(modelo,newdata=data.frame(x=c(valor₁,…,valor_k)),interval="prediction")

Se puede modificar el nivel de confianza incluyendo el argumento level=nivel_de_confianza (como proporción)

Diagramas de diagnóstico del modelo lm()

La función

plot(modelo)

proporciona cinco diagramas orientados al diagnóstico del modelo. Se puede indicar cuál de los diagramas se desea incluyendo el argumento which=valor

11.11.3 Correlación

Coeficiente de correlación de Pearson

r<-cor(x, y)
cor.test(x, y)

Proporcionan, respectivamente el coeficiente de correlación de Pearson y el test de correlación lineal. No es válida la especificación de una fórmula como argumento.

Coeficiente de correlación de Spearman

r<-cor(x, y, method="spearman")

En este caso, el test no es realizable cuando hay valores empatados.

11.11.4 Modelo Loess

Modelo de regresión no lineal. Se trata de un modelo no paramétrico. Permite diagnosticar visualmente la linealidad de la relación entre las variable implicada

modelo <- loess(y~x, data=datos)

Por defecto, el tamaño de la ventana es span=0.75. Se puede variar este tamaño entre 0 (ventana muy pequeña) y 1 (ventana muy grande)

modeloLoess <- loess(y~x, data=datos, span=valor)

Representación del modelo loess

Representación directa:

plot(x, modeloLoess$fitted, type="l",col="color", xlab="x",ylab="y")

Inclusión del modelo en un diagrama de dispersión ya elaborado:

points(x, modeloLoess$fitted, type="l"``)

La función car::scatterplot() representa el modelo loess con sus bandas de confianza cuando se indica el argumento smooth=TRUE

car::scatterplot(y~x,smooth=T)

11.12 Modelo lineal múltiple

modelo <- lm(y~x₁+…+x_k, data=datos)

El objeto de tipo formula

11.12.1 Diagrama de dispersión

Diagramas para cada variable

car::scatterplot(y~x_i, data=datos, smooth=T)

Matriz de diagramas de dispersión

pairs(y~x₁+…+x_k,data=datos,panel=panel.smooth)

el argumento panel=panel.smoothgenera la representación del modelo loess

car::scatterplotMatrix(y~x₁+…+x_k, data=datos)

11.12.2 Exploración de la linealidad con modelos GAM

library(mgcv)

modgam<-gam(y~s(x₁)+s(x₂)+…+s(x_p),data=datos)

plot(modgam)

11.12.3 Correlación entre las variables

cor(datos)

Puede ser conveniente eliminar alguna variable del data.frame (por ejemplo la primera, cuando esta representa el número de caso) y limitar el número de decimales.

round(cor(datos[,-1]),3)

Otros paquetes dan cierta versatilidad a la matriz a la matriz de correlaciones

Hmisc::rcorr(as.matrix(datos)) - muestra la significación de las correlaciones
PerformanceAnalytics::chart.Correlation(R=datos, histogram=TRUE)

11.12.4 Regresión por pasos

step(lm(y~x₁+…+x_k, data=datos), direction="both")

El argumento direction puede tomar valores "both", "forward" y "backward"

11.12.5 Validación del modelo

gvlma::gvlma(modelo)

11.12.6 Medidas y diagramas de influencia

influence.measures(modelo)

Diagramas de influencia

car::influencePlot(modelo, main=titulo)
olsrr::ols_plot_resid_stand(modelo) # Residuos estandarizados
olsrr::ols_plot_resid_stud(modelo) # Residuos estudentizados
olsrr::ols_plot_cooksd_bar(modelo) # Distancias de Cook:
olsrr::ols_plot_dfbetas(modelo) # DFBetas
olsrr::ols_plot_dffits(modelo) # DFFits
olsrr::ols_plot_resid_lev(modelo) # residuos estudentizados vs. leverage

11.12.7 Colinealidad

car::vif(modelo)

11.12.8 Modelo estandarizado

modelo_std<-lm(data.frame(scale(modelo$model)))

lm.beta::lm.beta(modelo)

11.12.9 Pronósticos

predict(modelo,newdata=data.frame(x₁=valor₁,x₂=valor₂),interval="confidence")
predict(modelo,newdata=data.frame(x₁=valor₁,x₂=valor₂),interval="prediction")

11.12.10 Incorporación de términos no lineales

Modelos polinómicos

Polinomio de tercer grado en x:

modelopol<-lm(y ~ x+I(x^2)+I(x^3), data=datos)