Regresión lineal
Lina Angélica Buitrago Reyes
1 Introducción
Los modelos de regresión lineal pertenecen a la familia de modelos lineales, es decir de aquellos en los cuales se explicar una variable respuesta (\(y\)) a partir de una combinación lineal de variables explicativas (\(x_j, j=1,...,p\)), así:
\[y=\beta_0+\beta_1x_1+...+\beta_px_p+e\] donde los \(\beta_j\)’s son los parámetros del modelo y \(e\) es el error entre lo ajustado por el modelo y la realidad.
Permite cumplir con dos objetivos, predecir el valor de la variable respuesta o determinar si existe relación lineal entre las \(x_j\) y \(y\).
2 Objetivos
- Comprender el proceso de estimación de los parámetros del modelo de regresión lineal múltiple, así como su interpretación
- Comprender el proceso inferencial en un modelo de regresión lineal múltiple
- Comprender y aprender a validar cada uno de los supuestos del modelo de regresión lineal múltiple
3 Competencias
Al finalizar este módulo el estudiante estará en la capacidad de ajustar, interpretar y validar un modelo de regresión múltiple.
4 Estimación de los parámetros
Existen dos caminos para estimar los parámetros del modelo, vía mínimos cuadrados o por máxima verosimilitud. El primero no exige supuesto alguno, mientras que en el segundo se debe asumir que los errores tienen distribución normal.
4.1 Mínimos cuadrados
Por construcción, la recta de mínimos cuadrados es aquella que pasa justo por el “centro” de la nube de puntos, por lo tanto la suma de los errores es siempre cero. Luego, con el objetivo de encontrar una recta que represente de la mejor forma el comportamiento de los datos, se estiman los parámetros minimizando:
\[\sum_{i=1}^{n}e_i^2=\sum_{i=1}^{n}(y_i-(\beta_0+\beta_1x_{i1}+...+\beta_px_{ip}))^2=\boldsymbol{e^te}\]
Siendo \(n\) el tamaño de la muestra,\(y_i\) es el valor de la variable respuesta para el \(i\)-ésimo individuo, \(x_{ij}\) es el valor de la j-ésima variable independiente para el individuo \(i\), \(e_i\) el error del modelo para el \(i\)-ésimo individuo y \(\boldsymbol{e^t}=(e_1,e_2,...,e_n)\).
Minimizando \(L\) se obtiene que el estimador de mínimos cuadrados está dado por: \[\boldsymbol{\hat\beta}=(\boldsymbol{X}^t\boldsymbol{X})^{-1}\boldsymbol{X}^t\boldsymbol{Y}\] Donde:\(\boldsymbol{\beta}^t=(\beta_0,\beta_1,...,\beta_p)\),\(\boldsymbol{Y}^t=(y_1,y_2,...,y_n)\) y \[\boldsymbol{X}= \left( \begin{array}{ll} 1 & x_{11} &...& x_{1p}\\ 1 & x_{21} &...& x_{2p}\\ \vdots & \vdots &...&\vdots \\ 1 & x_{n1} &...& x_{np} \end{array} \right)\]
4.2 Ejemplo-aplicación en R
#Antes de iniciar, descargue el archivo "capacidad_aerobica.txt" y copie la dirección de la carpeta
#en dónde lo guardó. Recuerde que debe cambiar el directorio de R a dicha carpeta, en mi caso:
#setwd("C:\\Users\\lange\\Google Drive\\Andes\\R")
#La siguiente es una base de datos de 31 pacientes.
datos<-read.table("capacidad_aerobica.csv",h=T, sep=";")
#Se busca determinar el estado físico de estos pacientes, razón por la cuál se midieron las siguientes
#variables
names(datos)## [1] "edad" "peso" "Oxig" "Tiempo"
## [5] "Pulso_reposo" "Pulso_corriend" "Pulso_maximo"#edad: edad del paciente en años
#peso: peso en Kg
#Oxig: Nivel de oxigenación
#Tiempo: tiempo de la prueba de esfuerzo (min)
#Pulso_reposo: Pulsaciones por minuto con el paciente en reposo
#Pulso_corriendo: Promedio de las pulsaciones por minuto durante la prueba de esfuerzo
#Pulso_maximo:
#Siempre antes de ajustar un modelo lineal se debe hacer un análisis descriptivo de los datos
summary(datos)## edad peso Oxig Tiempo
## Min. :38.00 Min. :59.08 Min. :37.39 Min. : 8.17
## 1st Qu.:44.00 1st Qu.:73.20 1st Qu.:44.96 1st Qu.: 9.78
## Median :48.00 Median :77.45 Median :46.77 Median :10.47
## Mean :47.68 Mean :77.44 Mean :47.38 Mean :10.59
## 3rd Qu.:51.00 3rd Qu.:82.33 3rd Qu.:50.13 3rd Qu.:11.27
## Max. :57.00 Max. :91.63 Max. :60.05 Max. :14.03
## Pulso_reposo Pulso_corriend Pulso_maximo
## Min. :40.00 Min. :146.0 Min. :155.0
## 1st Qu.:48.00 1st Qu.:163.0 1st Qu.:168.0
## Median :52.00 Median :168.0 Median :172.0
## Mean :53.45 Mean :169.1 Mean :173.8
## 3rd Qu.:58.50 3rd Qu.:176.0 3rd Qu.:180.0
## Max. :70.00 Max. :186.0 Max. :192.0#¿Qué tipo de población tiene en este estudio?, ¿son pacientes jóvenes?, ¿cómo está su condición física?
#¿cómo es la distribución de estas variables?, ¿simétricas?, ¿sesgadas?, ¿hay datos atípicos?
par(mfrow=c(2,4))
boxplot(datos$edad, main="Edad")
boxplot(datos$peso,main="Peso")
boxplot(datos$Oxig,main="Oxígeno")
boxplot(datos$Tiempo,main="Tiempo")
boxplot(datos$Pulso_reposo, main="Pulso en reposo")
boxplot(datos$Pulso_corriend, main="Pulso corriendo")
boxplot(datos$Pulso_maximo, main="Pulso máximo")
#Se busca determinar si el nivel de oxígeno está relacionado con las demás variables
#Inicialmente exploremos esta relación gráficamente: la siguiente gráfica nos muestra
#los gráficos de dispersión entre pares de variables de la base de datos
plot(datos)
#Se puede observar una relación lineal negativa clara entre el nivel de oxígeno y el tiempo (cor=-0.8621949), una menos
#fuerte con el pulso en reposo (cor=-0.39935611). También existe una relación lineal fuerte entre el pulso corriendo y
#el pulso máximo (cor=0.8513621), por lo que a la hora de pensar en un modelo de regresión lineal múltiple, no
#sería apropiado incluir las dos, ya que en primer lugar brindarían información redudante, además
#tendríamos problemas con el supuesto de multicolinealidad (La decisión sobre cuál de las dos incluir
#puede tomarse en dos sentidos: por la relevancia clínica o porque produce mejor ajuste).
#Lo anterior se puede corroborar a partir de la matriz de correlación de los datos
cor(datos)## edad peso Oxig Tiempo Pulso_reposo
## edad 1.0000000 -0.23353903 -0.3045924 0.1887453 -0.16409995
## peso -0.2335390 1.00000000 -0.1627528 0.1435076 0.04397417
## Oxig -0.3045924 -0.16275285 1.0000000 -0.8621949 -0.39935611
## Tiempo 0.1887453 0.14350758 -0.8621949 1.0000000 0.45038260
## Pulso_reposo -0.1640999 0.04397417 -0.3993561 0.4503826 1.00000000
## Pulso_corriend -0.3221848 0.11485895 -0.3017162 0.1762425 0.34844077
## Pulso_maximo -0.4329159 0.24938123 -0.2367402 0.2261030 0.30512400
## Pulso_corriend Pulso_maximo
## edad -0.3221848 -0.4329159
## peso 0.1148589 0.2493812
## Oxig -0.3017162 -0.2367402
## Tiempo 0.1762425 0.2261030
## Pulso_reposo 0.3484408 0.3051240
## Pulso_corriend 1.0000000 0.8513621
## Pulso_maximo 0.8513621 1.0000000#Ajustando un modelo con todas las variables, excepto el pulso máximo:
#y=beta0+beta1*edad+beta2*peso+...+beta5*Pulso_corriendo+e
m1<-lm(datos$Oxig~datos$edad+datos$peso+datos$Tiempo+datos$Pulso_reposo+datos$Pulso_corriend)
#se obtiene:
summary(m1)##
## Call:
## lm(formula = datos$Oxig ~ datos$edad + datos$peso + datos$Tiempo +
## datos$Pulso_reposo + datos$Pulso_corriend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1332 -1.0986 0.3669 1.1069 5.3256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 117.63416 11.70136 10.053 2.88e-10 ***
## datos$edad -0.26928 0.10051 -2.679 0.0129 *
## datos$peso -0.05665 0.05726 -0.989 0.3320
## datos$Tiempo -2.87843 0.39462 -7.294 1.21e-07 ***
## datos$Pulso_reposo -0.01268 0.07156 -0.177 0.8608
## datos$Pulso_corriend -0.12944 0.05176 -2.501 0.0193 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.476 on 25 degrees of freedom
## Multiple R-squared: 0.82, Adjusted R-squared: 0.7841
## F-statistic: 22.78 on 5 and 25 DF, p-value: 1.427e-08###Interpretación de las estimaciones
#beta0=117.63416: Se estima que el nivel de oxígeno en promedio tiene un valor de 117.63416 cuando todas las variables
# explicativas son cero
#beta1=-0.26928:Se estima que, en promedio, por cada año de edad que se incremente, el nivel de oxígeno disminuye
# 0.26928
#beta2=-0.05665:Se estima que, en promedio, por cada kg de peso que se incremente, el nivel de oxígeno disminuye
# 0.05665
#beta3=-2.87843:Se estima que, en promedio, por cada min que se incremente la prueba, el nivel de oxígeno disminuye
# 2.87843
#¿Cómo se interpretan las estimaciones de los dos parámetros que restan en el modelo?
#Nuestra base de datos no tiene variables categóricas, a modo de ejercicio vamos a crear
#una variable dicotómica que represente el sexo:
datos$sexo<-factor(c(rep("F",20),rep("M",dim(datos)[1]-20)))
#y la incluimos en el modelo:#y=beta0+beta1*edad+beta2*peso+...+beta5*Pulso_corriendo+beta*sexo(M)+e
m2<-lm(datos$Oxig~datos$edad+datos$peso+datos$Tiempo+datos$Pulso_reposo+datos$Pulso_corriend+datos$sexo)
summary(m2)##
## Call:
## lm(formula = datos$Oxig ~ datos$edad + datos$peso + datos$Tiempo +
## datos$Pulso_reposo + datos$Pulso_corriend + datos$sexo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6890 -0.9647 0.0826 0.9978 5.7985
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 119.39013 12.31550 9.694 8.98e-10 ***
## datos$edad -0.31266 0.13018 -2.402 0.0244 *
## datos$peso -0.04922 0.05973 -0.824 0.4180
## datos$Tiempo -2.87016 0.40066 -7.164 2.10e-07 ***
## datos$Pulso_reposo -0.01255 0.07261 -0.173 0.8642
## datos$Pulso_corriend -0.13303 0.05294 -2.513 0.0191 *
## datos$sexoM 0.70255 1.31067 0.536 0.5969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.512 on 24 degrees of freedom
## Multiple R-squared: 0.8222, Adjusted R-squared: 0.7777
## F-statistic: 18.49 on 6 and 24 DF, p-value: 6.358e-08#si no construímos nosotros mismos las variables dummys para las variables categóricas, R tomará
#uno de los niveles como base (el que no aparece en la tabla de resultados) y calculará las dummys para las restantes.
#En este caso tomó como base a la categoría femenino y muestra la estimación para la categoría masculino
###Interpretación de las estimaciones
#beta0=119.39013: Se estima que el nivel de oxígeno en promedio tiene un valor de 119.39013 cuando todas las variables
# explicativas son cero (en este punto hay que tener en cuenta que para este modelo
# la variable que incluímos de sexo es cero para las mujeres)
#beta1=-0.31266:Se estima que, en promedio, por cada año de edad que se incremente, el nivel de oxígeno disminuye
# 0.31266
#beta2=-0.04922:Se estima que por cada kg de peso que se incremente, el nivel de oxígeno disminuye
# 0.04922
#beta3=-2.87016:Se estima que, en promedio, por cada min que se incremente la prueba, el nivel de oxígeno disminuye
# 2.87016
#beta6=0.70255: Se estima que, en promedio, un hombre tiene un nivel de oxígeno de 0.70255 más que el de una mujer
###Estimación por intervalo
#La estimación por intervalo de los parámetros del modelo se obtiene así:
confint(m2)## 2.5 % 97.5 %
## (Intercept) 93.9721948 144.80807379
## datos$edad -0.5813300 -0.04398234
## datos$peso -0.1724891 0.07405061
## datos$Tiempo -3.6970858 -2.04322923
## datos$Pulso_reposo -0.1624056 0.13729862
## datos$Pulso_corriend -0.2422929 -0.02376358
## datos$sexoM -2.0025489 3.407647684.3 Máxima verosimilitud
Como vimos, para encontrar las estimaciones de mínimos cuadrados, no es necesario hacer supuesto alguno, simplemente éstas nos brindarán los parámetros de tal forma que se minimice la suma de cuadrados de los errores. En cambio, para encontrar los estimadores vía máxima verosimilitud se deben cumplir el numeral 5 de los siguientes supuestos.
4.3.1 Supuestos del modelo
Los siguientes supuestos son necesarios para verificar algunas propiedades de los estimadores, construir intervalos de confianza y hacer pruebas de hipótesis sobre los parámetros del modelo.
Independencia: Las mediciones de cada uno de los individuos son independientes (\(y_1,y_2,...y_n\) son estadísticamente independientes)
Linealidad: \(E(Y|x_1,x_2,...,x_n)=\beta_0+\beta_1x_1+...+\beta_px_p\)
\(E(e_i)=0\)
\(V(e_i)=\sigma^2\) (homocedasticidad)
\(e_i\sim N(0,\sigma^2)\) lo cual es equivalente a \(y_i|x_{i1},x_{i2},...,x_{in}\sim N(\beta_0+\beta_1x_1+...+\beta_px_p, \sigma^2)\)
Los errores son independientes, \(Cov(e_i,e_s)=0, i\neq s\)
No hay relación lineal entre las variables explicativas y el error, \(cov(e,x_j)=0, j=1,...,p\)
No hay multicolinealidad entre las variables explicativas, \(cov(x_j,x_k)=0, j,k=1,...,p, j\neq k\)
Para encontrar los estimadores máximo verosímiles, es necesario el supuesto 5, pues de allí parte la construcción de la función de verosimilitud:
\[L=\prod_{i=1}^n \frac{1}{2\pi\sigma}\exp\left\lbrace-\frac{1}{2}\frac{(y_i-(\beta_0+\beta_1x_1+...+\beta_px_p))^2}{\sigma^2}\right\rbrace\] Maximizar \(L\) es equivalente a minimizar la suma de cuadrados de los errores, por lo que el estimador encontrado es igual al encontrado por el método de mínimos cuadrados:
\[\boldsymbol{\hat\beta}=(\boldsymbol{X}^t\boldsymbol{X})^{-1}\boldsymbol{X}^t\boldsymbol{Y}\]
5 Inferencia-Pruebas de hipótesis para el modelo de RLM
Como ya sabemos, hay dos formas de hacer inferencia sobre parámetros, a partir de la estimación (puntual o por intervalo) o con pruebas de hipótesis. En la sección anterior vimos la forma de hacer las estimaciones, en esta nos concentraremos en las pruebas de hipótesis, las cuáles pueden ser generales, \(H_0: \beta_1=\beta_2=...=\beta_p=0\) o más específicas sobre un parámetro o un conjunto de éstos. En primer lugar debemos recordar que para que la inferencia que vamos a hacer sobre el modelo sea válida, es necesario que se cumplan los supuestos del modelo, es decir que las conclusiones que saquemos sobre nuestro problema sólo serán sólidas en el momento en que verifiquemos que los supuestos se cumplen.
5.1 Pruebas de Wald
Las pruebas de Wald se utilizan para probar si existe relación lineal entre cada una de las variables explicativas y la variable respuesta, en presencia de las demás variables independientes. El sistema de hipótesis está dado por:
\[ H_0: \beta_j=0\text{ vs }H_1: \beta_j\neq 0\] El estadístico de prueba es: \[t_c=\frac{\hat\beta_j}{\hat\sigma_{\hat\beta_j}}\]
y el test \(\tau:\)“Rechazar \(H_0\) si \(|t_c|>t_{1-\alpha/2(n-p-1)}\)”.
De igual forma se puede obtener un intervalo de confianza del \(100(1-\alpha)\%\), así:
\(L_{inf}=\hat\beta_j-t_{1-\alpha/2(n-p-1)}\hat\sigma_{\hat\beta_j}\) y \(L_{sup}=\hat\beta_j+t_{1-\alpha/2(n-p-1)}\hat\sigma_{\hat\beta_j}\)
5.1.1 Aplicación en R
## Ajuste del modelo
## y=beta0+beta1*edad+beta2*peso+...+beta5*Pulso_corriendo+e
m1<-lm(datos$Oxig~datos$edad+datos$peso+datos$Tiempo+datos$Pulso_reposo+datos$Pulso_corriend)
#Podemos obtener los resultados de las pruebas de Wald (t student) para los parámetros
#H0:betaj=0 a partir de (valor p=Pr(>|t|)):
summary(m1)##
## Call:
## lm(formula = datos$Oxig ~ datos$edad + datos$peso + datos$Tiempo +
## datos$Pulso_reposo + datos$Pulso_corriend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1332 -1.0986 0.3669 1.1069 5.3256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 117.63416 11.70136 10.053 2.88e-10 ***
## datos$edad -0.26928 0.10051 -2.679 0.0129 *
## datos$peso -0.05665 0.05726 -0.989 0.3320
## datos$Tiempo -2.87843 0.39462 -7.294 1.21e-07 ***
## datos$Pulso_reposo -0.01268 0.07156 -0.177 0.8608
## datos$Pulso_corriend -0.12944 0.05176 -2.501 0.0193 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.476 on 25 degrees of freedom
## Multiple R-squared: 0.82, Adjusted R-squared: 0.7841
## F-statistic: 22.78 on 5 and 25 DF, p-value: 1.427e-085.2 Análisis de varianza (ANOVA)
La tabla de análisis de varianza general permite determinar qué parte de la variabilidad de la variable respuesta es explicada por el modelo, de igual forma, probar si alguna de las variables explicativas tiene relación lineal con la respuesta y por lo tanto debería ser la primera prueba de hipótesis sobre los parámetros:
\[H_0:\beta_1=\beta_2=...=\beta_p=0\text{ vs }H_1:\text{Al menos un }\beta_j\neq 0\]
| Fuente de var | gl | SC | CM | F |
|---|---|---|---|---|
| Regresión | \(p\) | \(SCR=\sum_{i=1}^{n}(\hat y_i-\bar y)^2\) | \(CMR=SCR/p\) | \(CMR/CME\) |
| Error | \(n-p-1\) | \(SCE=\sum_{i=1}^{n}(\hat y_i-y_i)^2\) | \(CME=SCE/(n-p-1)\) | |
| Total | \(n-1\) | \(SCT=\sum_{i=1}^{n}(y_i-\bar y)^2\) |
El estadístico de prueba es \(F_c=CMR/CME\) y el test: "Rechazar \(H_0\) si \(F_c>f_{1-\alpha(p,n-p-1)}\).
5.2.1 Aplicación en R
#Matriz de diseño: contiene por columnas los valores de las variables explicativas para todos
#los individuos
X<-model.matrix(m1)
#tamaño de muestra
n<-dim(X)[1]
#número de variables
p<-dim(X)[2]-1 #la matriz de diseño incluye la columna del intercepto
#Podemos obtener las pruebas F que se calculan a partir de la descomposición de la
#variabilidad
anova(m1)## Analysis of Variance Table
##
## Response: datos$Oxig
## Df Sum Sq Mean Sq F value Pr(>F)
## datos$edad 1 78.99 78.99 12.8889 0.001407 **
## datos$peso 1 49.26 49.26 8.0380 0.008940 **
## datos$Tiempo 1 528.02 528.02 86.1598 1.41e-09 ***
## datos$Pulso_reposo 1 3.58 3.58 0.5838 0.451982
## datos$Pulso_corriend 1 38.32 38.32 6.2533 0.019315 *
## Residuals 25 153.21 6.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#en primer lugar, la suma de cuadrados total corregida por la media está dada por:
var(datos$Oxig)*(dim(datos)[1]-1)## [1] 851.3815#que es equivalente a sumar las sumas de cuadrados del anova
sct<-sum(anova(m1)$"Sum Sq")
#la suma de cuadrados del error está dada por
sce<-anova(m1)$"Sum Sq"[p+1]
#por lo tanto la suma de cuadrados de la regresión:
scr<-sct-sce
scr## [1] 698.1713## [1] 698.1713#Para probar H0: todos los betaj son cero vs H1: al menos un betaj es diferente de cero (esta debe ser el primer sistema a probar),
#utilizamos el estadístico de prueba F=CMR/CME (CMR:cuadrado medio de la regresión corregido por
#la media, CME: cuadrado medio del error)
#Como ya calculamos las sumas de cuadrados, ahora podemos calcular los cuadrados medios:
#CMR
cmr<-scr/sum(anova(m1)$"Df"[1:p])
cmr## [1] 139.6343## [1] 6.12841## [1] 22.78475## [1] 2.6029875.2.2 Sumas de cuadrados secuenciales
Las sumas de cuadrados secuenciales, es decir que va probando si incluir las variables en orden resulta significativo el parámetro que las acompaña,mientras que con la prueba de Wald vemos si la variable resulta significativa en presencia de todas las demás que componen el modelo. Es decir, las pruebas que se presenta en una anova con sumas de cuadrados secuenciales se presentan así:
- \(H_0:\beta_1=0\) vs \(H_1:\beta_1\neq0\), en el modelo \(y=\beta_0+\beta_1x_1+e\)
- \(H_0:\beta_2=0\) vs \(H_1:\beta_2\neq0\), en el modelo \(y=\beta_0+\beta_1x_1+\beta_2x_2+e\)
- \(H_0:\beta_3=0\) vs \(H_1:\beta_3\neq0\), en el modelo \(y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_3+e\)
- Así sucesivamente…
5.2.2.1 Aplicación en R
## Analysis of Variance Table
##
## Response: datos$Oxig
## Df Sum Sq Mean Sq F value Pr(>F)
## datos$edad 1 78.99 78.99 12.8889 0.001407 **
## datos$peso 1 49.26 49.26 8.0380 0.008940 **
## datos$Tiempo 1 528.02 528.02 86.1598 1.41e-09 ***
## datos$Pulso_reposo 1 3.58 3.58 0.5838 0.451982
## datos$Pulso_corriend 1 38.32 38.32 6.2533 0.019315 *
## Residuals 25 153.21 6.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# y notemos que los resultados en términos de significancia son diferentes a los obtenidos con:
summary(m1)##
## Call:
## lm(formula = datos$Oxig ~ datos$edad + datos$peso + datos$Tiempo +
## datos$Pulso_reposo + datos$Pulso_corriend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1332 -1.0986 0.3669 1.1069 5.3256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 117.63416 11.70136 10.053 2.88e-10 ***
## datos$edad -0.26928 0.10051 -2.679 0.0129 *
## datos$peso -0.05665 0.05726 -0.989 0.3320
## datos$Tiempo -2.87843 0.39462 -7.294 1.21e-07 ***
## datos$Pulso_reposo -0.01268 0.07156 -0.177 0.8608
## datos$Pulso_corriend -0.12944 0.05176 -2.501 0.0193 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.476 on 25 degrees of freedom
## Multiple R-squared: 0.82, Adjusted R-squared: 0.7841
## F-statistic: 22.78 on 5 and 25 DF, p-value: 1.427e-08#¿Por qué ocurre ésto?
#Recordemos que existen diferentes sumas de cuadrados, por defecto R calcula las sumas de cuadrados
#secuenciales, es decir que va probando si incluir las variables en orden resulta significativo,
#mientras que con la prueba de Wald vemos si la variable resulta significativa en presencia de todas las
#demás que componen el modelo, veamos:
#la suma de cuadrados en la tabla anova correspondiente a la variable edad es equivalente a la
#que se obtiene del modelo:
anova(lm(datos$Oxig~datos$edad))## Analysis of Variance Table
##
## Response: datos$Oxig
## Df Sum Sq Mean Sq F value Pr(>F)
## datos$edad 1 78.99 78.988 2.9657 0.0957 .
## Residuals 29 772.39 26.634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#es decir que estamos probando si incluir la variable edad a un modelo que tiene únicamente el intercepto
#resulta significativo. La suma de cuadrados para la variable peso se obtiene de:
anova(lm(datos$Oxig~datos$edad+datos$peso))## Analysis of Variance Table
##
## Response: datos$Oxig
## Df Sum Sq Mean Sq F value Pr(>F)
## datos$edad 1 78.99 78.988 3.0585 0.09127 .
## datos$peso 1 49.26 49.260 1.9074 0.17818
## Residuals 28 723.13 25.826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#luego, estamos probando si resulta significativo que se incluya el peso al modelo que ya contiene la edad,
#diferente de lo que probamos con la prueba de Wald, pues en ésta probamos si la variable es significativa
#en presencia de todas las demás. Así sucesivamente se obtienen las restantes sumas de cuadrados secuenciales:
anova(lm(datos$Oxig~datos$edad+datos$peso+datos$Tiempo))## Analysis of Variance Table
##
## Response: datos$Oxig
## Df Sum Sq Mean Sq F value Pr(>F)
## datos$edad 1 78.99 78.99 10.9306 0.002679 **
## datos$peso 1 49.26 49.26 6.8167 0.014561 *
## datos$Tiempo 1 528.02 528.02 73.0694 3.668e-09 ***
## Residuals 27 195.11 7.23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Variance Table
##
## Response: datos$Oxig
## Df Sum Sq Mean Sq F value Pr(>F)
## datos$edad 1 78.99 78.99 10.7224 0.002993 **
## datos$peso 1 49.26 49.26 6.6869 0.015669 *
## datos$Tiempo 1 528.02 528.02 71.6775 6.016e-09 ***
## datos$Pulso_reposo 1 3.58 3.58 0.4857 0.492053
## Residuals 26 191.53 7.37
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Variance Table
##
## Response: datos$Oxig
## Df Sum Sq Mean Sq F value Pr(>F)
## datos$edad 1 78.99 78.99 12.8889 0.001407 **
## datos$peso 1 49.26 49.26 8.0380 0.008940 **
## datos$Tiempo 1 528.02 528.02 86.1598 1.41e-09 ***
## datos$Pulso_reposo 1 3.58 3.58 0.5838 0.451982
## datos$Pulso_corriend 1 38.32 38.32 6.2533 0.019315 *
## Residuals 25 153.21 6.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16 Residuales
6.1 Residual clásico
El residual clásico es la diferencia entre el valor real de la variable dependiente y su estimación por medio del modelo de regresión:
\[\hat{e_i}=y_i-\hat{y_i}\]
Sin embargo, estos residuales no tienen varianza constante, por lo que no serían los residuales apropiados para evaluar los supuestos de homocedasticidad y normalidad, ya que:
\[V(\hat{\boldsymbol{e}})=\sigma^2(I-H)\]
Con \(H=X(X^tX)^{-1}X^t\).
6.2 Residuales estandarizados
Como su nombre lo indica, corresponden a la estandarización de los residuales clásicos:
\[r_i=\frac{\hat{e_i}}{\sqrt{(1-h_{ii})CME}}\] ## Residuales studentizados
Los residuales estudentizados permiten encontrar datos influyentes, se definen como:
\[r_i^*=\frac{\hat{e_i}}{\sqrt{(1-h_{ii})\hat{\sigma_{(i)}^2}}}\] donde \(\hat{\sigma_{(i)}^2}\) es el CME de la regresión obtenida sin el i-ésimo elemento.
7 Supuestos y su verificación
7.1 Supuestos del modelo de regresión lineal múltiple
Para que se cumpla que la distribución de los estadísticos de prueba sean las que vimos en la sección anterior, es necesario que el modelo ajustado cumpla con algunos supuestos
7.1.1 Linealidad
El supuesto de linealidad significa que efectivamente existe relación lineal entre la varibale independiente y las variables explicativas, es decir:
\[E(Y|x_1,...,x_p)=\beta_0+\beta_1x_1+...+\beta_px_p\] Lo cuál puede verificarse de la siguiente así:
#En la sección anterior ajustamos el modelo
m1<-lm(datos$Oxig~datos$edad+datos$peso+datos$Tiempo+datos$Pulso_reposo+datos$Pulso_corriend)
summary(m1)##
## Call:
## lm(formula = datos$Oxig ~ datos$edad + datos$peso + datos$Tiempo +
## datos$Pulso_reposo + datos$Pulso_corriend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1332 -1.0986 0.3669 1.1069 5.3256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 117.63416 11.70136 10.053 2.88e-10 ***
## datos$edad -0.26928 0.10051 -2.679 0.0129 *
## datos$peso -0.05665 0.05726 -0.989 0.3320
## datos$Tiempo -2.87843 0.39462 -7.294 1.21e-07 ***
## datos$Pulso_reposo -0.01268 0.07156 -0.177 0.8608
## datos$Pulso_corriend -0.12944 0.05176 -2.501 0.0193 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.476 on 25 degrees of freedom
## Multiple R-squared: 0.82, Adjusted R-squared: 0.7841
## F-statistic: 22.78 on 5 and 25 DF, p-value: 1.427e-08# y vimos que las variables peso y Pulso_reposo no son significativas, sacando el Pulso_reposo, que es
#la variable con mayor valor p, se obtiene:
m2<-lm(datos$Oxig~datos$edad+datos$peso+datos$Tiempo+datos$Pulso_corriend)
summary(m2)##
## Call:
## lm(formula = datos$Oxig ~ datos$edad + datos$peso + datos$Tiempo +
## datos$Pulso_corriend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1727 -1.0245 0.3605 1.0230 5.3242
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 117.38330 11.39702 10.299 1.14e-10 ***
## datos$edad -0.26546 0.09632 -2.756 0.0105 *
## datos$peso -0.05554 0.05585 -0.994 0.3292
## datos$Tiempo -2.91088 0.34300 -8.487 5.74e-09 ***
## datos$Pulso_corriend -0.13152 0.04947 -2.659 0.0132 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.429 on 26 degrees of freedom
## Multiple R-squared: 0.8198, Adjusted R-squared: 0.7921
## F-statistic: 29.57 on 4 and 26 DF, p-value: 2.459e-09#Como la edad sigue siendo no significativa, el modelo queda:
m2<-lm(datos$Oxig~datos$edad+datos$Tiempo+datos$Pulso_corriend)
summary(m2)##
## Call:
## lm(formula = datos$Oxig ~ datos$edad + datos$Tiempo + datos$Pulso_corriend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7949 -1.2439 0.2537 1.0508 4.8400
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 112.55751 10.31016 10.917 2.09e-11 ***
## datos$edad -0.24121 0.09316 -2.589 0.0153 *
## datos$Tiempo -2.97637 0.33655 -8.844 1.85e-09 ***
## datos$Pulso_corriend -0.13115 0.04945 -2.652 0.0132 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.429 on 27 degrees of freedom
## Multiple R-squared: 0.813, Adjusted R-squared: 0.7922
## F-statistic: 39.12 on 3 and 27 DF, p-value: 5.738e-10#Dado que el modelo tiene más de una variable y que el R2 se ve afectado por el número de variables
#en el modelo, se analiza la variabilidad explicada por el modelo a partir del R2 ajustado. En este
#caso R2 aj=0.813, es decir que el modelo explica el 81.3% de la variabilidad del nivel de oxígeno
####Linealidad
#estimaciones del modelo para los datos incluidos
yest<-predict(m2)
#gráfica de y (nivel de oxígeno) vs yest (estimaciones del nivel de oxígeno con el modelo estimado)
plot(datos$Oxig,yest)
##un indicativo de un buen ajuste del modelo es que la nube de puntos de la gráfica anterior
#esté lo más cerca posible de la recta y=x
abline(a=0, b=1, col="red")
7.1.2 Independencia
El supuesto de independencia implica que los errores del modelo no deben estar correlacionados, es decir: \[cor(e_i, e_j)=0\] Este supuesto no se cumpliría por ejemplo si tenemos datos de un estudio antes-después o cuando, por ejemplo, se tienen pacientes de un mismo hospital o de una misma familia, que haga pensar que tanto las variables dependientes como independientes se comporten de forma similar.
Este supuesto se puede verificar a partir de los residuales estandarizados, así:
####Residuales
#para obtener los residuales de un modelo de regresión lineal utilizaremos el comando residuals
?residuals
#residuales ordinarios
res<-residuals(m2)
plot(res)
#No se observa un patrón definido, lo cual es deseable, por otro lado se puede ver que toman valores
#en un rango entre -4 (subestimación del nivel de oxígeno) y 4 (sobrestimación del nivel de oxígeno)
#Recordemos que los residuales ordinarios no tienen varianza constante, por lo que se utilizan los
#residuales estandarizados para verificar los supuestos del modelo:
restan<-rstandard(m2)
####Para verificar de forma exploratoria el supuesto de no correlación:
#i vs ri
plot(restan)
#Los residuales se pueden ubicar dentro de una franja con amplitud constante (-2,2), es decir
#que no se presentan datos atípicos.Tampoco se observa alguna tendencia, adicionalmente, los datos por el diseño del estudio
#son independientes por el diseño del estudio, ya que provienen de pacientes diferentesEn general, tampoco se cumple cuando se tienen mediciones en el tiempo, en este caso se puede verificar si existe o no dicha autocorrelación con la prueba de Durbin Watson.
7.1.3 Varianza constante
También llamada homocedasticidad y significa que la variabilidad del error debe ser constante: \[V(e_i)=\sigma^2\] Se pueden verificar gráficamente y comprobar inferencialmente a partir de la prueba de Bartlett:
#En general, los residuales están dentro de una franja constante, excepto aquellos que corresponden
#a estimaciones menores a 44. Con la prueba de Bartlett veremos si esta diferencia es significativa.
#H0: Hay homocedasticidad vs H1: No hay homocedasticidad
#Para hacer esta prueba, en primer lugar es necesario construir grupos en los cuáles veamos que el
#comportamiento de la varibilidad de los residuales es similar, en este caso los construiremos
#observando la gráfica anterior:
grupos<-numeric()
grupos[yest<41]<-1
grupos[yest>=41 & yest<48]<-2
grupos[yest>=48 & yest<53]<-3
grupos[yest>=53]<-4
table(grupos)## grupos
## 1 2 3 4
## 4 11 12 4##
## Bartlett test of homogeneity of variances
##
## data: restan by grupos
## Bartlett's K-squared = 3.9665, df = 3, p-value = 0.2651#Como p-value = 0.2651>0.05 no se rechaza la hipótesis nula, por lo tanto no hay evidencia de
#heterocedasticidad y se verifica el supuesto.Otras pruebas de homocedasticidad son: - Razón de verosimilitud - Hartley - Bartlett - Hartley - Levene - O’Brien
7.1.4 Normalidad
Los errores deben tener distribución normal: \[e_i\sim N(0,\sigma^2)\]
#Prueba de normalidad
#H0: Los residuales tienen distribución normal vs H1: Los residuales no tienen distribución normal
shapiro.test(restan)##
## Shapiro-Wilk normality test
##
## data: restan
## W = 0.97037, p-value = 0.5293#Como p-value = 0.5293>0.05 no se rechaza la hipótesis, luego se verifica el supuesto de normalidad.Otras pruebas de normalidad son:
- Jarque-Bera
- Kolmogorov-Smirnov
7.1.5 No multicolinealidad
Significa que las variables independientes no deben estar correlacionadas entre sí, puesto esto causa que no se encuentren estimaciones únicas para los parámetros del modelo.
Se puede explorar inicialmente a partir de las correlaciones entre las variables independientes, o calculando el factor de inflación de la varianza, el cuá para la variable independiente \(j\) se define como:
\[VIF_j=\frac{1}{1-R_j^2}\] Donde \(R_j^2\) es el coeficiente de determinación de \(x_j\) vs las demás variables explicativas.
Si no hay multicolinealidad \(VIF_j=1\), si la multicolinealidad es moderada. \(5<VIF_j<10\), si hay alta multicolinealidad \(VIF_j\geq10\).
###Multicolinealidad
#Exploración descriptiva
plot(data.frame(datos$edad,datos$peso,datos$Tiempo,datos$Pulso_corriend))
cor(data.frame(datos$edad,datos$peso,datos$Tiempo,datos$Pulso_corriend))## datos.edad datos.peso datos.Tiempo datos.Pulso_corriend
## datos.edad 1.0000000 -0.2335390 0.1887453 -0.3221848
## datos.peso -0.2335390 1.0000000 0.1435076 0.1148589
## datos.Tiempo 0.1887453 0.1435076 1.0000000 0.1762425
## datos.Pulso_corriend -0.3221848 0.1148589 0.1762425 1.0000000#Tanto en la gráfica como en la matriz de correlación no se evidencia una fuerte relación
#lineal entre las variables explicativas, por lo que el supuesto de multicolinealidad también
#se verifica, sin embargo calcularemos los VIFs:
#la función qye utilizaremos se encuentra en el paquete faraway
library(faraway)
## datos$edad datos$Tiempo datos$Pulso_corriend
## 1.199069 1.109051 1.1934227.1.6 No relación entre las variables explicativas y el error
En caso de existir relación entre alguna variable explicativa y el error, significa que el modelo aún no ha logrado explicar completamente la relación entre dicha variable y la variable respuesta, para lo cual se hace un análisis exploratorio gráfico, que puede dar luz sobre lo que está ocurriendo.
#Es importante graficar cada una de las variables explicativas vs los residuales, si el modelo
#ajustado es el correcto, estas gráficas deben ser nubes de puntos sin tendencia alguna, si por
#el contrario muestran algún comportamiento, por ejemplo cuadrático, significa que debemos incluir
#dicho efecto en el modelo
par(mfrow=c(2,2))
plot(datos$edad,restan)
plot(datos$peso,restan)
plot(datos$Tiempo,restan)
plot(datos$Pulso_corriend,restan)
7.2 Datos influyentes
Los datos influyentes son aquellos que potencialmente pueden cambiar la pendiente de la recta ajustada. Algunos de los métodos para identificarlos son:
- Puntos leverage: Mide la influencia de cada dato con base en resultados de la matriz de diseño del modelo
- Distancia de Cook: Compara las estimaciones de los parámetros del modelo con todos los datos vs el modelo sin cada uno de ellos.
- Residuos studentizados
##Datos influyentes
##Cálculo de Leverage
#Matriz de diseño
X<-model.matrix(m2)
#tamaño de muestra
n<-dim(X)[1]
#número de variables en el modelo
p<-dim(X)[2]-1
#Matriz H
H<-X%*%solve(t(X)%*%X)%*%t(X)
#puntos de palanca
#leverage:elementos de la diagonal de la matriz H
h<-diag(H)
plot(h)
#identificación de valores con leverage alto
abline(h=2*(p+1)/n,col="red")
## [1] 12##Distancia de Cook
cook<-cooks.distance(m2)
plot(cook)
#identificación de elementos con distancia de Cook alta
abline(h=qf(0.5,p+1,n-p-1), col="red")
#ninguna distancia de cook está por encima de qf(0.5,p,n-p)=0.8082103
##Residuos studentizados
stud<-rstudent(m2)
plot(stud)
## [1] 3.504931#todos los residuales studentizados tienen valores inferiores al percentil de la distribución t
#Algunas de las gráficas de diagnóstico pueden obtenerse a partir de:
par(mfrow=c(2,2))
plot(m2)
8 Soluciones al no cumplimiento de los supuestos
8.1 Transformación de Box-Cox
La transformación de Box-Cox es utilizada cuandono se cumplen los supuestos de normalidad o varianza constante. Se define así:
\[z=\frac{y_i^\lambda-1}{\lambda y°^{\lambda-1}}, \lambda\neq 0\] \[y° \ln y_i, \lambda=0\] Donde \(y°\) es la media gemoétrica de los valores de \(y\). La estimación de lambda se hace utilizando el método de máxima verosimilitud, que en R se implementa así:
lambda <- t$x # valores de lambda
l<- t$y # valores de la log verosimilitud
bc <- cbind(lambda, l)
bc <- bc[order(-l),] # ordena los valores de lambda de mayor a menor verosimilitud
head(bc, n = 10)## lambda l
## [1,] -0.5454545 42.41810
## [2,] -0.5858586 42.41734
## [3,] -0.5050505 42.41667
## [4,] -0.6262626 42.41441
## [5,] -0.4646465 42.41306
## [6,] -0.6666667 42.40930
## [7,] -0.4242424 42.40726
## [8,] -0.7070707 42.40203
## [9,] -0.3838384 42.39926
## [10,] -0.7474747 42.39259