Caso cuando alguna \(x_{i}\) es categórica

Cargamos y vemos los datos

Datos=read.csv("ejemplo3RLM.csv", header=TRUE )

summary(Datos)
##        y              X1c                  X2        
##  Min.   :-2.093   Length:300         Min.   :-1.618  
##  1st Qu.: 3.793   Class :character   1st Qu.: 1.848  
##  Median : 5.528   Mode  :character   Median : 2.912  
##  Mean   : 5.704                      Mean   : 3.069  
##  3rd Qu.: 7.538                      3rd Qu.: 4.264  
##  Max.   :15.984                      Max.   : 9.482
str(Datos)
## 'data.frame':    300 obs. of  3 variables:
##  $ y  : num  3.25 4.23 9.58 5.11 5.34 ...
##  $ X1c: chr  "A1" "A1" "A1" "A1" ...
##  $ X2 : num  1.88 2.54 6.12 3.14 3.26 ...

vemos que la variable x1c viene como cadena. Realizamos las transformaciones correspondientes

Datos$X1c=factor(Datos$X1c)
str(Datos)
## 'data.frame':    300 obs. of  3 variables:
##  $ y  : num  3.25 4.23 9.58 5.11 5.34 ...
##  $ X1c: Factor w/ 3 levels "A1","A2","A3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ X2 : num  1.88 2.54 6.12 3.14 3.26 ...

Después graficamos

par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(Datos)

La gráfica de nuestro interés es la siguiente, que además hemos agregado color (lo que siempre es recomendable)

with(Datos, plot(X2, y, col=c("red", "brown", "blue")[Datos[,2]] ))
legend("topleft",levels(Datos[,2]), col=c("red", "brown", "blue"), pch = c(0,0,0),  pt.cex=1.5,cex = .9, y.intersp = 0.3 , bty="n")

donde comparamos la relación de la variable \(X_{2}\) con la \(y\) y adicionalmente la posible información que podría brindar la variable categórica. Después vemos los niveles de x1c

levels(Datos$X1c)
## [1] "A1" "A2" "A3"

de modo que el modelo con el que estamos trabajando es

\[ \mathbb{E}(y|x)=\beta_{0}+\beta_{1}A_2 + \beta_{2}A_{3}+ \beta_{3}X_{2} \]

donde \(A_{2}\) y \(A_{3}\) son variables dicotómicas y \(X_{2}\) es la variable continua.

fit=lm(y~X1c+X2, data=Datos)
summary(fit)
## 
## Call:
## lm(formula = y ~ X1c + X2, data = Datos)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.277652 -0.060507  0.003333  0.065288  0.261476 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.505342   0.013849   36.49   <2e-16 ***
## X1cA2       1.313092   0.014053   93.44   <2e-16 ***
## X1cA3       0.499562   0.014002   35.68   <2e-16 ***
## X2          1.497182   0.003044  491.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.099 on 296 degrees of freedom
## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9988 
## F-statistic: 8.162e+04 on 3 and 296 DF,  p-value: < 2.2e-16

Analicemos la Interpretación del modelo con estos parámetros. Después regresaremos a la prueba asociada a la tabla ANOVA.

Aquí \(\mathbb{E}(y|x1c, x2)=\beta_0+\beta_1 A_2+\beta_2 A_3+\beta_3 x_2\) con el nivel de referencia igual a A1, además notar que en este caso, si consideramos que x1c toma el valor del nivel A1,entonces el modelo para \(\mathbb{E}(Y|x1c=A_1,x_2)\) es

\[ \mathbb{E}(Y|x1c=A_1,x_2)=\beta_0+\beta_1(0)+\beta_2(0)+\beta_3x_2 =\beta_0+\beta_3x_2 \]

para nivel A2: \(\mathbb{E}(Y|x1c=A_2,x_2)\)

\[ \mathbb{E}(Y|x1c=A_1,x_2)=\beta_0+\beta_1(1)+\beta_2(0)+\beta_3x_2 =(\beta_0+\beta_{1})+\beta_3x_2 \]

para nivel A3: \(\mathbb{E}(Y|x1c=A_3,x_2)\)

\[ \mathbb{E}(Y|x1c=A_1,x_2)=\beta_0+\beta_1(0)+\beta_2(1)+\beta_3x_2 =(\beta_0+\beta_{2})+\beta_3x_2 \]

Es decir, con este modelo se considera que la \(\mathbb{E}(Y|x1c,x_2)\) es una recta en función del valor de \(x_2\) para cada nivel, con diferente intercepto pero misma pendiente (lo cual se veía en la gráfica de las tres rectas coloreadas).

Otras formas de definir el mismo modelo en R (cuestiones computacionales)

  • Indicando las 2 variebles dicótomicas (niveles) que deben entrar al modelo

    fitAlt=lm(y~I(X1c=="A2")+I(X1c=="A3")+X2, data=Datos)
summary(fitAlt)
    ## 
## Call:
## lm(formula = y ~ I(X1c == "A2") + I(X1c == "A3") + X2, data = Datos)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.277652 -0.060507  0.003333  0.065288  0.261476 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.505342   0.013849   36.49   <2e-16 ***
## I(X1c == "A2")TRUE 1.313092   0.014053   93.44   <2e-16 ***
## I(X1c == "A3")TRUE 0.499562   0.014002   35.68   <2e-16 ***
## X2                 1.497182   0.003044  491.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.099 on 296 degrees of freedom
## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9988 
## F-statistic: 8.162e+04 on 3 and 296 DF,  p-value: < 2.2e-16

  • Cambiando el nivel de referencia que tiene por default R. Aquí hay que obtener de nuevo las expresiones de \(\mathbb{E}(Y|x1c=A1,x_2)\), \(\mathbb{E}(Y|x1c=A_2,x_2)\) y \(\mathbb{E}(Y|x1c=A_3,x_2)\)

    fitAlt2=lm(y~relevel(X1c, "A3")+X2, data=Datos)
summary(fitAlt2)
    ## 
## Call:
## lm(formula = y ~ relevel(X1c, "A3") + X2, data = Datos)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.277652 -0.060507  0.003333  0.065288  0.261476 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1.004904   0.013977   71.89   <2e-16 ***
## relevel(X1c, "A3")A1 -0.499562   0.014002  -35.68   <2e-16 ***
## relevel(X1c, "A3")A2  0.813530   0.014070   57.82   <2e-16 ***
## X2                    1.497182   0.003044  491.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.099 on 296 degrees of freedom
## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9988 
## F-statistic: 8.162e+04 on 3 and 296 DF,  p-value: < 2.2e-16
    fitAlt3=lm(y~I(X1c=="A1")+I(X1c=="A2")+X2, data=Datos)
summary(fitAlt3)
    ## 
## Call:
## lm(formula = y ~ I(X1c == "A1") + I(X1c == "A2") + X2, data = Datos)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.277652 -0.060507  0.003333  0.065288  0.261476 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.004904   0.013977   71.89   <2e-16 ***
## I(X1c == "A1")TRUE -0.499562   0.014002  -35.68   <2e-16 ***
## I(X1c == "A2")TRUE  0.813530   0.014070   57.82   <2e-16 ***
## X2                  1.497182   0.003044  491.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.099 on 296 degrees of freedom
## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9988 
## F-statistic: 8.162e+04 on 3 and 296 DF,  p-value: < 2.2e-16

    El efecto de la función relevel:

    levels(Datos$X1c)
    ## [1] "A1" "A2" "A3"
    Datos$X1c=relevel(Datos$X1c, "A3")
levels(Datos$X1c)
    ## [1] "A3" "A1" "A2"

    El efecto de la función factor

    Datos$X1c=factor(Datos$X1c, levels=c("A1","A2","A3"))
levels(Datos$X1c)
    ## [1] "A1" "A2" "A3"

Regresemos al modelo 

summary(fit)
## 
## Call:
## lm(formula = y ~ X1c + X2, data = Datos)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.277652 -0.060507  0.003333  0.065288  0.261476 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.505342   0.013849   36.49   <2e-16 ***
## X1cA2       1.313092   0.014053   93.44   <2e-16 ***
## X1cA3       0.499562   0.014002   35.68   <2e-16 ***
## X2          1.497182   0.003044  491.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.099 on 296 degrees of freedom
## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9988 
## F-statistic: 8.162e+04 on 3 and 296 DF,  p-value: < 2.2e-16

Consideramos la prueba F asociada a la tabla anova.

Observemos que con la prueba asociada a la tabla ANOVA, se puede concluir que con una significancia de .05, se rechaza \(H_0\) en el contraste

\[ H_0: \beta_1=0\ \ y\ \ \beta_2=0 \ \ y \ \ \beta_3=0 \ \ vs\ \ H_a: \beta_1\neq0\ \ o \ \ \beta_2\neq0\ \ o\ \ \beta_3\neq0 \]

Además, si analizamos la prueba t, para el coeficiente \(\beta_3\), asociado a \(X_2\), podemos observar que se rechaza \(H_0\) en el contraste \(H_0: \beta_3=0\ \ vs\ \ Ha: \beta_3 \neq 0\),esto nos indicaría que una vez incluida la variablecategórica X1c en el modelo, la variable \(X_2\) está agregando información.

Para analizar si la variable X1c agrega información dado que el modelo incluye a \(X_2\), debemos constrastar \(H_0: \beta_1=0\ \ y \ \ \beta_2=0\ \ vs\ \ H_a: \beta_1\neq0 \ \ o \ \ \beta_2\neq0\). En este caso se debe tener cuidado puesto que no podemos emplear el mismo análisis anterior dado que x1c es categórica.

#Opción 1, usando multcomp
library(multcomp)
K=matrix(c(0,1,0,0,
           0,0,1,0), ncol=4, 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   1.3131
## 2 == 0   0.4996
## 
## Global Test:
##      F DF1 DF2     Pr(>F)
## 1 4442   2 296 1.763e-221

se rechaza \(H_{0}\) entonces la variable X1c nos proporciona información para modelar \(\mathbb{E}(Y|x)\) aún cuando en el modelo ya está la variable \(X_2\).

#Opción 2, usando anova con modelos reducidos
fitred=lm(y~X2, data=Datos)
anova(fitred, fit)
## Analysis of Variance Table
## 
## Model 1: y ~ X2
## Model 2: y ~ X1c + X2
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    298 89.985                                  
## 2    296  2.901  2    87.084 4442.3 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

de igual manera el \(p-value\) es pequeño.

#Opción 3, usar directamente la función drop1
drop1(fit, test = "F")
## Single term deletions
## 
## Model:
## y ~ X1c + X2
##        Df Sum of Sq     RSS      AIC  F value    Pr(>F)    
## <none>                 2.90 -1383.59                       
## X1c     2     87.08   89.99  -357.24   4442.3 < 2.2e-16 ***
## X2      1   2370.54 2373.44   626.49 241851.0 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Opcón 4, usar directamente la función Anova en librería car
library(car)
Anova(fit, type="II")
## Anova Table (Type II tests)
## 
## Response: y
##            Sum Sq  Df  F value    Pr(>F)    
## X1c         87.08   2   4442.3 < 2.2e-16 ***
## X2        2370.54   1 241851.0 < 2.2e-16 ***
## Residuals    2.90 296                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notar que lo anterior sólo se vale cuando la variable categórica aparece en el modelo sólo a través de sus \(k-1\) variables dicotómicas (no hay interacciones entre variables).

Por otro lado, si analizamos la prueba t, para el coeficiente \(\beta_1\), se rechaza \(H_0\). Aquí el contraste es \(H_0: \beta_1=0\ \ vs\ \ H_a: \beta_1 \neq 0\). Esto nos indica, que condicional en un valor fijo de \(X_2\), la esperanza de \(Y\) es diferente entre el nivel \(A_2\) y el \(A_1\) (de referencia).

Ahora, si analizamos la prueba t, para el coeficiente \(\beta_2\), se rechaza \(H_0\). Aquí el contraste es \(H_0:\ beta_2=0\ \ vs\ \ H_a: \beta_2 \neq 0\). Esto nos indica, que condicional en un valor fijo de \(X_2\), la esperanza de \(Y\) es diferente entre el nivel \(A_3\) y el \(A_1\) (de referencia)

Con base en lo anterior, parece que no podríamos reducir el modelo, todos los coeficientes parecen significativos.

Interpretación del modelo 

summary(fit)
## 
## Call:
## lm(formula = y ~ X1c + X2, data = Datos)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.277652 -0.060507  0.003333  0.065288  0.261476 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.505342   0.013849   36.49   <2e-16 ***
## X1cA2       1.313092   0.014053   93.44   <2e-16 ***
## X1cA3       0.499562   0.014002   35.68   <2e-16 ***
## X2          1.497182   0.003044  491.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.099 on 296 degrees of freedom
## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9988 
## F-statistic: 8.162e+04 on 3 and 296 DF,  p-value: < 2.2e-16

\(R^2\), el coeficiente de determinación, nos indica que se está explicando el \(99.9\%\) de la variabilidad observada en \(Y\) a través del modelo que incluye ambas variables X1c y \(X_2\): \(\mathbb{E}(y|X1c, X_2)=\beta_0+\beta_1A_2+\beta_2A_3+\beta_3X_2\).

Además, b1 se puede interpretar como: condicionado en un valor fijo de \(X_2\), el promedio de la variable. Y aumenta en 1.31 unidades al comparar el nivel \(A_2\) contra el \(A_1\). (mayor esperanza en nivel \(A_2\) contra nivel \(A_1\)). Recordar que nivel \(A_1\) es el de referencia en este modelo. Esta interpretación se obtiene al comparar

\[ \mathbb{E}(Y|A_2, X_2)-\mathbb{E}(Y|A_1, X_2)=\beta_0+\beta_1+\beta_3 X_2-(\beta_0+\beta_3 x_2)=\beta_1 \]

\(\beta_2\) se puede interpretar como: condicionado en un valor fijo de \(X_2\), el promedio de la variable \(Y\) aumenta en .5 unidades al comparar el nivel \(A_3\) contra el \(A_1\).

\(\beta_3\) se puede interpretar como: condicionado en un nivel fijo de la variable X1c, el promedio de la variable \(Y\) aumenta en 1.5 unidades al aumentar en una unidad la variable \(X_2\).

¿Cómo se comparan los niveles A2 y A3, dado un valor fijo de X2?

\[\begin{align*} &\mathbb{E}(Y|A_2, X_2)=\beta_0+\beta_1+\beta_3 X_2\\ &\mathbb{E}(Y|A_3, X_2)=\beta_0+\beta_2+\beta_3 X_2 \end{align*}\]

entonces

\[ \mathbb{E}(Y|A_3, X_2)-\mathbb{E}(Y|A_2, X_2)=\beta_2-\beta_1 \]

Por ejemplo si se tiene interés en contrastar

\[ H_0:\mathbb{E}(Y|A_2, X_2)=\mathbb{E}(Y|A_3, X_2)\ \ vs\ \ H_a:\mathbb{E}(Y|A_2, X_2)\neq \mathbb{E}(Y|A_3, X_2) \]

Se puede usar multcomp

library(multcomp)
K=matrix(c(0,-1,1,0), ncol=4, nrow=1, byrow=TRUE)
m=c(0)
summary(glht(fit, linfct=K, rhs=m), test=Ftest())
## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0  -0.8135
## 
## Global Test:
##      F DF1 DF2     Pr(>F)
## 1 3343   1 296 2.545e-163
summary(glht(fit, linfct=K, rhs=m))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ X1c + X2, data = Datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0 -0.81353    0.01407  -57.82   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza la igualdad de medias, es decir, dado un valor de \(X_2\), la esperanza de \(Y\) es diferente entre los niveles \(A_3\) y \(A_2\). De hecho, la esperanza disminuye en .8135 unidades al comparar el nivel \(A_3\) contra el \(A_2\).

Visualización de los resultados del modelo

\(\mathbb{E}(y|X1c,X_2)=\beta_0+\beta_1A_2+\beta_2A_3+\beta_3 X_2\)

Recordar que en este caso, el modelo se relaciona con tres rectas, una por nivel. Además recordemos que como función de \(X_2\), son tres rectas que tienen diferente intercepto pero misma pendiente \(\beta_3\)

fit$coefficients
## (Intercept)       X1cA2       X1cA3          X2 
##   0.5053418   1.3130916   0.4995617   1.4971818
#coef(fit)
fA1 <- function(X2) {fit$coefficients[1]+ fit$coefficients[4]*X2}
fA2 <- function(X2) {fit$coefficients[1]+fit$coefficients[2]+ fit$coefficients[4]*X2}
fA3 <- function(X2) {fit$coefficients[1]+fit$coefficients[3]+ fit$coefficients[4]*X2}

with(Datos, plot(X2, y, col=c("red", "green", "blue")[Datos[,2]] ))
legend("topleft",levels(Datos[,2]), col=c("red", "green", "blue"), pch = c(0,0,0), pt.cex=1.5,cex = .9, y.intersp = 0.4 , bty="n" )
curve(fA1, from = min(Datos$X2), to = max(Datos$X2),
      col = "red", add = T)
curve(fA2, from = min(Datos$X2), to = max(Datos$X2),
      col = "green", add = T)
curve(fA3, from = min(Datos$X2), to = max(Datos$X2),
      col = "blue", add = T)