Ajuste e interpretación de Variables categóricas

Comenzamos por cargar nuestros datos:

Datos <- read.csv("ejemplo3RLM.csv", header=TRUE )
summary(Datos)
##        y            X1c                  X2      
##  Min.   :-2.1   Length:300         Min.   :-1.6  
##  1st Qu.: 3.8   Class :character   1st Qu.: 1.8  
##  Median : 5.5   Mode  :character   Median : 2.9  
##  Mean   : 5.7                      Mean   : 3.1  
##  3rd Qu.: 7.5                      3rd Qu.: 4.3  
##  Max.   :16.0                      Max.   : 9.5
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 ...

notamos que la variable x1c no es una variable continua. Para darnos una mejor idea sobre nuestros datos veamos el siguiente gráfico

library(GGally)
ggpairs(data=Datos, title="Datos", aes(color=X1c))

donde vemos que nuestra variable x1c es categórica de tres niveles; vemos también un cierto comportamiento lineal entre la variable x2 y y para los distintos niveles de x1c. Asimismo, parece ser que las medias de los distintos grupos son distintas.

R siempre considera como nivel de referencia el primer nivel guardado en el metadato de etiquetas, no obstante, la variable categórica no ha sido especificada como del tipo factor:

levels(Datos$X1c)
## NULL

es decir, R no reconoce que x1c sea una variable categórica. Para declarar que dicha variable sea de tipo factor escribimos:

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

Graficamos el comportamiento de x2 con y distinguiendo cada uno de los grupos:

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

El modelo que se ajusta por default en R incluyendo sólo las dos variables X1c y X2 es: \(E(y;x)=\beta_0+\beta_1A_2+\beta_2A_3+\beta_3 x_2\). Procedemos a ajustar el modelo

fit <- lm(y ~ X1c + X2, data=Datos)
summary(fit)
## 
## Call:
## lm(formula = y ~ X1c + X2, data = Datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27765 -0.06051  0.00333  0.06529  0.26148 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.50534    0.01385    36.5   <2e-16 ***
## X1cA2        1.31309    0.01405    93.4   <2e-16 ***
## X1cA3        0.49956    0.01400    35.7   <2e-16 ***
## X2           1.49718    0.00304   491.8   <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.999,  Adjusted R-squared:  0.999 
## F-statistic: 8.16e+04 on 3 and 296 DF,  p-value: <2e-16

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

Aquí \(E(y;x1c, x2)=\beta_0+\beta_1A_2+\beta_2A_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 \(E(Y;x1c=A1,x2)\) es

\[ E(Y;x1c=A1,x_2)=\beta_0+\beta_1(0)+\beta_2(0)+\beta_3x_2=\beta_0+\beta_3x_2 \] para nivel A2: \[ E(Y;x1c=A2,x_2)=\beta_0+\beta_1(1)+\beta_2(0)+\beta_3x_2=(\beta_0+\beta_1)+\beta_3x_2 \] para nivel A3:

\[ E(Y;x1c=A3,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 \(E(Y;x1c,x2)\) es una recta en función del valor de x2 para cada nivel, con diferente intercepto pero misma pendiente.

Otras formas de definir el mismo modelo en R (cuestiones computacionales) indicando las 2 variables 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.27765 -0.06051  0.00333  0.06529  0.26148 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.50534    0.01385    36.5   <2e-16 ***
## I(X1c == "A2")TRUE  1.31309    0.01405    93.4   <2e-16 ***
## I(X1c == "A3")TRUE  0.49956    0.01400    35.7   <2e-16 ***
## X2                  1.49718    0.00304   491.8   <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.999,  Adjusted R-squared:  0.999 
## F-statistic: 8.16e+04 on 3 and 296 DF,  p-value: <2e-16

donde I() es una función indicadora.

Podemos cambiar el nivel de referencia que tiene por default R. Aquí hay que obtener de nuevo las expresiones de \(E(Y;x1c=A1,x_2), E(Y;x1c=A2,x_2)\) y \(E(Y|;x1c=A3,x_2)\)

# Indicamos que el nivel de referencia para el modelo
# sera A3 en vez del A1:
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.27765 -0.06051  0.00333  0.06529  0.26148 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1.00490    0.01398    71.9   <2e-16 ***
## relevel(X1c, "A3")A1 -0.49956    0.01400   -35.7   <2e-16 ***
## relevel(X1c, "A3")A2  0.81353    0.01407    57.8   <2e-16 ***
## X2                    1.49718    0.00304   491.8   <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.999,  Adjusted R-squared:  0.999 
## F-statistic: 8.16e+04 on 3 and 296 DF,  p-value: <2e-16

o alternativamente

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.27765 -0.06051  0.00333  0.06529  0.26148 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.00490    0.01398    71.9   <2e-16 ***
## I(X1c == "A1")TRUE -0.49956    0.01400   -35.7   <2e-16 ***
## I(X1c == "A2")TRUE  0.81353    0.01407    57.8   <2e-16 ***
## X2                  1.49718    0.00304   491.8   <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.999,  Adjusted R-squared:  0.999 
## F-statistic: 8.16e+04 on 3 and 296 DF,  p-value: <2e-16

podemos trabajar de manera directa sobre los niveles:

# El efecto de la función relevel:
# Niveles dados por default
levels(Datos$X1c)
## [1] "A1" "A2" "A3"
# Cambiamos que A3 sea el nivel de referencia
Datos$X1c <- relevel(Datos$X1c, "A3")
levels(Datos$X1c)
## [1] "A3" "A1" "A2"
# Podemos especificar explicitamente los niveles y el orden
Datos$X1c <- factor(Datos$X1c, levels=c("A1","A2","A3"))
levels(Datos$X1c)
## [1] "A1" "A2" "A3"

Regresamos al análisis del modelo.

summary(fit)
## 
## Call:
## lm(formula = y ~ X1c + X2, data = Datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27765 -0.06051  0.00333  0.06529  0.26148 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.50534    0.01385    36.5   <2e-16 ***
## X1cA2        1.31309    0.01405    93.4   <2e-16 ***
## X1cA3        0.49956    0.01400    35.7   <2e-16 ***
## X2           1.49718    0.00304   491.8   <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.999,  Adjusted R-squared:  0.999 
## F-statistic: 8.16e+04 on 3 and 296 DF,  p-value: <2e-16

Observemos que con la prueba asociada a la tabla ANOVA, se puede concluir que con una significancia de .05, que 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 X2, podemos observar que se rechaza \(H_0\) en el contraste

\[ H_0: \beta_3=0 \ \ vs\ \ H_a: \beta_3 \neq 0 \]

Esto nos indicaría que una vez incluida la variable categórica X1c en el modelo, la variable X2 está agregando información. Para analizar si la variable X1c agrega información dado que el modelo incluye a X2, debemos contrastar

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

para lo cual:

# Opcion 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.31
## 2 == 0     0.50
## 
## Global Test:
##      F DF1 DF2    Pr(>F)
## 1 4442   2 296 1.76e-221

Se observa que se rechaza \(H_0\), entonces la variable X1c nos proporciona información para modelar \(E(Y;x)\) aún cuando en el modelo ya está la variable X2.

# Opcion 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 90.0                             
## 2    296  2.9  2      87.1 4442 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Opcion 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   2    4442 <2e-16 ***
## X2          2371   1  241851 <2e-16 ***
## Residuals      3 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).

En cualquier caso se está rechazando la hipótesis nula.

summary(fit)
## 
## Call:
## lm(formula = y ~ X1c + X2, data = Datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27765 -0.06051  0.00333  0.06529  0.26148 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.50534    0.01385    36.5   <2e-16 ***
## X1cA2        1.31309    0.01405    93.4   <2e-16 ***
## X1cA3        0.49956    0.01400    35.7   <2e-16 ***
## X2           1.49718    0.00304   491.8   <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.999,  Adjusted R-squared:  0.999 
## F-statistic: 8.16e+04 on 3 and 296 DF,  p-value: <2e-16

Retomando la información sobre nuestro modelo, 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 condicionar en un valor fijo de X2, la esperanza de Y es diferente entre el nivel A2 y el A1 (de referencia). De manera análoga para el coeficiente \(\beta_{2}\), donde al condicionar en un valor fijo de \(X_2\), la esperanza de \(Y\) es diferente entre el nivel A3 y el A1 (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.27765 -0.06051  0.00333  0.06529  0.26148 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.50534    0.01385    36.5   <2e-16 ***
## X1cA2        1.31309    0.01405    93.4   <2e-16 ***
## X1cA3        0.49956    0.01400    35.7   <2e-16 ***
## X2           1.49718    0.00304   491.8   <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.999,  Adjusted R-squared:  0.999 
## F-statistic: 8.16e+04 on 3 and 296 DF,  p-value: <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 X2: \(E(y;X1c, X2)=\beta_0+\beta_1A_2+\beta_2A_3+\beta_3x_2\).

Además, \(\beta_1\) 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 A2 contra el A1. (mayor esperanza en nivel A2 contra nivel A1) Recordar que nivel A1 es el de referencia en este modelo. Esta interpretación se obtiene al comparar

\[ E(Y;A2, X2)-E(Y;A1, X2)=\beta_0+\beta_1+\beta_3x_2-(\beta_0+\beta_3x_2)=\beta_1 \]

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

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

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

\(E(Y;A2, X2)=\beta_o+\beta_1+\beta_3X_2\)

\(E(Y;A3, X2)=\beta_o+\beta_2+\beta_3X_2\)

\(E(Y;A3, X2)-E(Y;A2, X2)=\beta_2-\beta_1\)

Por ejemplo si se tiene interés en contrastar

\[ H_0:E(Y;A2, X2) =E(Y;A3, X2)\ \ vs \ \ H_a:E(Y;A2, X2)\neq E(Y;A3, X2) \]

Las hipótesis se pueden escribir como \(H_0:\beta_2-\beta_1=0 \ \ vs \ \ H_a:\beta_2-\beta_1\neq0\), donde 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.814
## 
## Global Test:
##      F DF1 DF2    Pr(>F)
## 1 3343   1 296 2.55e-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.8135     0.0141   -57.8   <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 X2, la esperanza de Y es diferente entre los niveles A3 y A2 De hecho, la esperanza disminuye en .8135 unidades al comparar el nivel A3 contra el A2.

Para visualizar resultados del modelo \(E(y;X1c,X2)=\beta_0+\beta_1A_2+\beta_2A_3+\beta_3x_2\), recordar que en este caso dicho modelo se relaciona con tres rectas, una por nivel:

  • Nivel A1: \(E(Y;A1,x2)=\beta_0+\beta_3x_2\).
  • Nivel A2: \(E(Y;A2,x2)=(\beta_0+\beta_1)+\beta_3x_2\).
  • Nivel A3: \(E(Y;A3,x2)=(\beta_0+\beta_2)+\beta_3x_2\).

Es decir, como función de X2, son tres rectas que tienen diferente intercepto pero misma pendiente b3. Para ello

# Accedemos a los coeficiente betas
fit$coefficients
## (Intercept)       X1cA2       X1cA3          X2 
##        0.51        1.31        0.50        1.50

Código para graficar las rectas

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=0.9, y.intersp=1.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)