Regresión lineal general

La prueba lineal general asociada a un modelo de regresión lineal múltiple es una extensión de la pruea \(F\) asociada a la tabla ANOVA. En esta prueba \(H_{0}\) corresponde a la igualdad forma simultánea de \(r\) combinaciones lineales de los parámetros con constantes. Es decir

\[ H_{0}: \textbf{K}\beta=\textbf{m} \ \ \ \ \ vs \ \ \ \ \ H_{a}: \textbf{k}\beta\neq \textbf{m} \]

donde \(\textbf{K}\) es una matriz de costantes de \(r\times (p+1)\) y \(rank(\textbf{K})=r\leq p+1\).

Ejemplo:

Comenzamos por cargar y ver los datos

Datos=read.csv("ejemplo1RLM.csv", header=TRUE )
summary(Datos)
##        y                X1              X2       
##  Min.   : 5.810   Min.   :10.38   Min.   :15.40  
##  1st Qu.: 7.989   1st Qu.:14.01   1st Qu.:18.95  
##  Median : 8.651   Median :15.12   Median :20.07  
##  Mean   : 8.711   Mean   :15.18   Mean   :20.15  
##  3rd Qu.: 9.589   3rd Qu.:16.38   3rd Qu.:21.18  
##  Max.   :12.554   Max.   :19.37   Max.   :24.55

de modo que tenemos en este caso dos variables explicativas.

El scatterplot es más díficil de usar en la regresión lineal múltiple pues sólo se muestran relaciones entre pares de variables.

par(mfrow=c(1,2),mar=c(4.5,4.5,1,1))
pairs(Datos)

de modo que debemos tener extremo cuidado al extraer información de dichos scatterplots. Se recalca, no son tan utilizado en el caso de la RLM.

En la regresión lineal múltiple nos preguntamos, por ejemplo, en cómo es el efecto de \(X_2\) en \(\mathbb{E}(Y|X_1,X_2)\), así que debemos considerar que en \(\mathbb{E}(Y|X_1,X_2)\) aparecen tanto \(X_1\) y \(X_2\) variables.

Un posible modelo es:

\[ \mathbb{E}(y|X_1, X_2)=\beta_0+\beta_{1}X_{1}+\beta_{2}X_{2} \]

o en R

fit=lm(y~X1+X2, data=Datos)
summary(fit)
## 
## Call:
## lm(formula = y ~ X1 + X2, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8730 -0.6607 -0.1245  0.6214  2.0798 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.73700    1.85987   0.934    0.353
## X1           0.35404    0.33155   1.068    0.288
## X2           0.07937    0.32998   0.241    0.810
## 
## Residual standard error: 0.9513 on 97 degrees of freedom
## Multiple R-squared:  0.4134, Adjusted R-squared:  0.4013 
## F-statistic: 34.18 on 2 and 97 DF,  p-value: 5.829e-12

donde para agregar más variables al ajuste del modelo sólo debemos agregarlas sumando.

Además tenemos algunas estimaciones puntuales

#\hat{\beta}
coef(fit)
## (Intercept)          X1          X2 
##  1.73699710  0.35404328  0.07937096
#\hat{V(\hat{\beta})}
vcov(fit)
##             (Intercept)         X1         X2
## (Intercept)   3.4591232  0.5079432 -0.5539388
## X1            0.5079432  0.1099237 -0.1080313
## X2           -0.5539388 -0.1080313  0.1088883
#\hat{\sigma}=\sqrt{SCE/(n-p-1)}
sigma(fit)
## [1] 0.9512979

Como bien sabemos, debemos verificar que nuestro modelo tiene sentido, para ello emplearemos la prueba asociada a la tabla ANOVA, donde contrastamos

\[ H_{0}: \beta_{1}=0 \ \ \textrm{y}\ \ \beta_{2}=0 \ \ \ vs\ \ \ H_{a}: \beta_{1}\neq0 \ \ \textrm{o}\ \ \beta_{2}\neq0 \] Del summary tenemos que \(p-value= 5.829e-12\) por lo que se rechaza \(H_0\). Entonces al menos uno de los coeficientes asociados a las variables es diferente de cero y nos ayuda a modelar \(\mathbb{E}(y|x)\). Sin embargo, dicha prueba no nos dice cuál variable es la que se asegura distinta de cero.

Obtención de la Prueba de la Tabla Anova usando la Prueba Lineal General.

Definiremos la matriz k que nos permite recuperar a \(\beta_{1}\) y \(\beta_{2}\), y utilizaremos el vector m=c(0,0) para realizar el contraste en la prueba F.

# Cargamos la librería 
library(multcomp)

# Definimos la matriz correspondiente
K=matrix(c(0,1,0,
           0,0,1), ncol=3, nrow=2, byrow=TRUE)
m=c(0,0)
summary(glht(fit, linfct=K, rhs=m), test=Ftest())
## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0  0.35404
## 2 == 0  0.07937
## 
## Global Test:
##       F DF1 DF2    Pr(>F)
## 1 34.18   2  97 5.829e-12

siendo lo anterior un camino alternativo al del summary para verificar que nuestro modelo tiene sentido.

Una vez que se realiza la prueba asociada a la tabla Anova la pregunta es si ambas variables aportan información al modelado para eso podemos hacer pruebas individuales sobre cada coeficiente usando las pruebas t o bien la prueba lineal general

Las pruebas t que se obtienen del summary están asociadas a cada coeficiente de forma “no simultánea”. Para su interpretación se sugiere identificar qué modelo resulta de considerar el parámetro asociado igual a cero.

En este caso tenemos \(\beta_0\), \(\beta_1\) y \(\beta_2\). Luego, recordemos del summary

summary(fit)
## 
## Call:
## lm(formula = y ~ X1 + X2, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8730 -0.6607 -0.1245  0.6214  2.0798 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.73700    1.85987   0.934    0.353
## X1           0.35404    0.33155   1.068    0.288
## X2           0.07937    0.32998   0.241    0.810
## 
## Residual standard error: 0.9513 on 97 degrees of freedom
## Multiple R-squared:  0.4134, Adjusted R-squared:  0.4013 
## F-statistic: 34.18 on 2 and 97 DF,  p-value: 5.829e-12

donde Pr(<|t|) representa el \(p-value\) obtenido de la prueba donde se contrasta que \(\beta_{i}=0\) versus \(\beta_{i}\neq 0\). Son estos \(p-values\) los obtenidos de las pruebas no simultáneas.

Veamos un ejemplo de cómo interpretaríamos esta información:

La prueba \(H_0:\) \(\beta_1=0\) vs \(H_a: \beta_1 \neq 0\), nos lleva a concluir que no hay evidencia para rechazar \(H_0\) con una significancia de \(\alpha=.05\). Se toma el \(p-value\) del summary, renglón asociado a \(X_1\): \(p-value=0.29 > .05\).

El modelo que considera que \(\beta_1=0\), sería uno en el que sólo está \(X_2\) es decir, el modelo reducido \(\mathbb{E}(y|X_1, X_2)=\beta_0+\beta_2 X_2\). En este caso, la prueba está elaborada para responder: ¿La inclusión de la variable \(X_1\) una vez que se tiene el modelo \(y=\beta_0+\beta_2 X_2+e\) nos está o no agregando información adicional; en otras palabras, condicionar en \(X_2\), \(X_1\) nos agrega información adicional para modelar \(\mathbb{E}(Y|X)\)?

De tal manera, en nuestro caso la variable \(X_{1}\) no está agregando mayor información al modelo reducido previamente considerado \(\mathbb{E}(y|X_1, X_2)=\beta_0+\beta_2 X_2\). De froma análogo, dado que \(Pr(<|t|)=0.810\) para la variable \(X_{2}\) en el summary, entonces agregar a la variable \(X_{2}\) al modelo reducido \(\mathbb{E}(y|X_1, X_2)=\beta_0+\beta_1 X_1\) no nos estaría brindando mayor información.

Las pruebas t se enfocan en una sóla combinación de los parámetros. También se pueden obtener de forma directa con el paquete multcomp de dos maneras:

#H0: b1=0 vs Ha: b1 != 0
library(multcomp)
K=matrix(c(0,1,0), ncol=3, nrow=1, byrow=TRUE)
m=c(0)
#usando prueba F equivalente (no permite pruebas de una cola con alternativa > o <)
summary(glht(fit, linfct=K, rhs=m), test=Ftest())
## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0    0.354
## 
## Global Test:
##      F DF1 DF2 Pr(>F)
## 1 1.14   1  97 0.2882
#usando prueba t (permite prueba de una cola)
#Notar que sólo eliminamos el argumento test=
summary(glht(fit, linfct=K, rhs=m))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ X1 + X2, data = Datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)
## 1 == 0   0.3540     0.3315   1.068    0.288
## (Adjusted p values reported -- single-step method)
#Ejemplo de una prueba de una cola
#H0: b1<=0 vs Ha: b1 > 0
summary(glht(fit, linfct=K, rhs=m, alternative = "greater") )
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ X1 + X2, data = Datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)
## 1 <= 0   0.3540     0.3315   1.068  0.144
## (Adjusted p values reported -- single-step method)

Nota: Se puede ver que las pruebas t de dos colas, así como las pruebas que se pueden incluir en la prueba lineal general son casos particulares en donde se comparan modelos anidados el reducido vs el completo. R tiene otras formas de obtener las pruebas usando el lenguaje de modelos anidados, por ejemplo, la función anova nos sirve para comparar dos modelos anidados sin necesidad de usar la definición matricial (\(K\) y \(m\)).

Por ejemplo si \(\beta_1=0\), entonces el modelo reducido es \(y=\beta_0+\beta_2 X_2+e\). Lo ajustamos en R

fitred1=lm(y~X2, data=Datos)
summary(fitred1)
## 
## Call:
## lm(formula = y ~ X2, data = Datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.91862 -0.70204 -0.08324  0.59596  2.22091 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.10101    1.05526   0.096    0.924    
## X2           0.42732    0.05216   8.192 9.84e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.952 on 98 degrees of freedom
## Multiple R-squared:  0.4065, Adjusted R-squared:  0.4004 
## F-statistic: 67.12 on 1 and 98 DF,  p-value: 9.842e-13

Ahora comparamos el modelo completo y el reducido

anova(fitred1, fit)
## Analysis of Variance Table
## 
## Model 1: y ~ X2
## Model 2: y ~ X1 + X2
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     98 88.814                           
## 2     97 87.782  1    1.0319 1.1403 0.2882
# no importa el orden de los parámetros
anova(fit,fitred1)
## Analysis of Variance Table
## 
## Model 1: y ~ X1 + X2
## Model 2: y ~ X2
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     97 87.782                           
## 2     98 88.814 -1   -1.0319 1.1403 0.2882

donde \(p-value=0.28\), es decir que no se rechaza \(H_{0}\), i.e. es plausible asumir el modelo reducido.

Cálculos a mano 

#(SCEred-SCEc)/SCEc  * (n-p-1)/r
SCEred=sigma(fitred1)^2*fitred1$df.residual 
SCEcom=sigma(fit)^2*fit$df.residual 
#n-p-1 son los grados de libertad del modelo completo
#r las combinaciones lineales necesarias para obtener el reducido
( (SCEred-SCEcom)/SCEcom *(fit$df.residual/1) )
## [1] 1.140306

Algo similar para obtener la prueba de la tabla anova: en ese caso el modelo reducido es \(y=\beta_0+e\) cuando \(H_0: \beta_1=0\) y \(\beta_2=0\)

fitred2=lm(y~1, data=Datos)
anova(fitred2, fit)
## Analysis of Variance Table
## 
## Model 1: y ~ 1
## Model 2: y ~ X1 + X2
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     99 149.638                                  
## 2     97  87.782  2    61.856 34.176 5.829e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rechazándose \(H_{0}\).

¿Qué modelo elegir?

Selección entre los modelos. Notemos que

summary(fit)
## 
## Call:
## lm(formula = y ~ X1 + X2, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8730 -0.6607 -0.1245  0.6214  2.0798 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.73700    1.85987   0.934    0.353
## X1           0.35404    0.33155   1.068    0.288
## X2           0.07937    0.32998   0.241    0.810
## 
## Residual standard error: 0.9513 on 97 degrees of freedom
## Multiple R-squared:  0.4134, Adjusted R-squared:  0.4013 
## F-statistic: 34.18 on 2 and 97 DF,  p-value: 5.829e-12
fitred1=lm(y~X2, data=Datos)
summary(fitred1)
## 
## Call:
## lm(formula = y ~ X2, data = Datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.91862 -0.70204 -0.08324  0.59596  2.22091 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.10101    1.05526   0.096    0.924    
## X2           0.42732    0.05216   8.192 9.84e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.952 on 98 degrees of freedom
## Multiple R-squared:  0.4065, Adjusted R-squared:  0.4004 
## F-statistic: 67.12 on 1 and 98 DF,  p-value: 9.842e-13
fitred2=lm(y~X1, data=Datos)
summary(fitred2)
## 
## Call:
## lm(formula = y ~ X1, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8820 -0.6696 -0.1450  0.6054  2.0528 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.14077    0.79684   2.687  0.00848 ** 
## X1           0.43279    0.05212   8.304 5.67e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9467 on 98 degrees of freedom
## Multiple R-squared:  0.413,  Adjusted R-squared:  0.407 
## F-statistic: 68.96 on 1 and 98 DF,  p-value: 5.672e-13

Aquí casi los tres modelos indican un mismo coeficiente de determinación, se seleccionaría en general el de menor número de parámetros y con el mayor coeficiente de determinación aunque más adelante se considerarán otros criterios (BIC y AIC) o bien se realiza un análisis del poder predictivo del modelo. De tal manera, el modelo con el cual nos quedaremos puede ser aquel que sólo contempla a \(X_1\) o a \(X_{2}\).

Alerta de error común

La lectura de las pruebas t NO debe hacerse de forma simultánea, ver por ejemplo

summary(fit)
## 
## Call:
## lm(formula = y ~ X1 + X2, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8730 -0.6607 -0.1245  0.6214  2.0798 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.73700    1.85987   0.934    0.353
## X1           0.35404    0.33155   1.068    0.288
## X2           0.07937    0.32998   0.241    0.810
## 
## Residual standard error: 0.9513 on 97 degrees of freedom
## Multiple R-squared:  0.4134, Adjusted R-squared:  0.4013 
## F-statistic: 34.18 on 2 and 97 DF,  p-value: 5.829e-12

muchas personas cometen el error de leer la salida de forma simultánea concluyendo en este caso que ni \(X_1\) ni \(X_2\) aportan información al modelo lo que en este caso no es cierto basado en el análisis realizado de forma correcta con la prueba asociada a la tabla ANOVA.