Problema tipo Ancova

Teoría

En este caso, además de considerar diferentes grupos como en el caso Anova, los investigadores sospechan sobre un efecto diferenciado con base en una variable continua observada en la misma muestra.

Ejemplo simple:

\[\begin{align*} x_{1}=\left\{\begin{matrix}\textrm{grupo 1}\\ \textrm{grupo 2}\end{matrix}\right. \ \ \ x_{g1}=\left\{\begin{matrix}1 \ \ \textrm{si grupo 1}\\ 0 \ \ \textrm{c.o.c}\end{matrix}\right. \ \ \ x_{2} \ \ \textrm{continua} \end{align*}\]

con el modelo

\[ y_{i}=\beta_{0}+\beta_{1}\underbrace{x_{i g_{1}}}_{\textrm{efecto principal}}+\beta_{2}x_{i2}+\beta_{3}\underbrace{x_{ig1}x_{i1}}_{\textrm{interacción}}+\varepsilon_{i} \]

donde \(\varepsilon\sim N(\textbf{0},\sigma^{2}\mathbb{I}_{n\times n})\).

Para interpretar el modelo, basta encontrar las expresiones de la esperanza de \(y\) para las diferentes categorías de \(x_{1}\), dejando fijo el valor de \(x_{2}\). De tal manera se tiene

\[\begin{align*} &\mathbb{E}(y; \textrm{grupo 1}, x_{2})=(\beta_{0}+\beta_{1})+(\beta_{2}+\beta_{3})x_{2}\\ &\mathbb{E}(y; \textrm{grupo 2}, x_{2})=\beta_{0}+\beta_{2}x_{2} \end{align*}\]

Notamos de las expresiones anteriores, dos rectas en función de \(x_{2}\), una por grupo con diferente intercepto y diferente pendiente.

Se presentan algunas preguntas de interés:

  • El modelo tiene sentido

\[ H_{0}: \beta_{1}=\beta_{2}=\beta_{3}=0 \ \ vs \ \ H_{a}: \beta_{j}\neq0 \ \ \textrm{para alguna j} \]

  • La variable \(x_{2}\) tiene efecto al menos en algún grupo

\[ H_{0}: \beta_{2}=0 \ \ y \ \ \beta_{3}=0 \ \ vs \ \ H_{a}: \beta_{2}\neq0 \ \ o \ \ \beta_{3}\neq0 \] * La variable \(x_{1}\) tiene efecto para algún valor de \(2_{2}\)

\[ H_{0}: \beta_{1}=0 \ \ y \ \ \beta_{3}=0 \ \ vs \ \ H_{a}: \beta_{1}\neq0 \ \ o \ \ \beta_{3}\neq0 \]

  • El efecto de la variable \(x_{2}\) es el mismo en las diferentes clases de \(x_{1}\) (rectas con la misma pendiente; se pueden quitar las restricciones)

Ejemplo

ANCOVA caso 1. Igualdad de pendientes. Continuación de ejemplo. 

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 ...
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 ...

Cuando tenemos una variable categórica y una varibale continua lo recomendable es realizar el siguiente scatterplot

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 = 1.3 , bty="n")

donde parecen tener una nube de puntos (uno por recta), asimismo, parece que las tres rectas tienen la misma pendiente con diferente intercepto.

Vemos los niveles de nuestra variable

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

Se empieza con el modelo que incluye interacciones

\[ \mathbb{E}(y;x1c, x_2)= \beta_0 + \beta_1 A_2 + \beta_2 A_3 + \beta_3 X_2 + \beta_4(A_2\cdot X_2) + \beta_5(A_3\cdot X_2) \]

fit=lm(y~X1c*X2, data=Datos) otra alternativa para indicar todas las interacciones entre dos variables

fit=lm(y~X1c+X2+X1c:X2, data=Datos)
summary(fit)
## 
## Call:
## lm(formula = y ~ X1c + X2 + X1c:X2, data = Datos)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.278987 -0.060431  0.002723  0.065852  0.261398 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.504355   0.019996  25.223   <2e-16 ***
## X1cA2        1.324839   0.026535  49.928   <2e-16 ***
## X1cA3        0.488498   0.028057  17.411   <2e-16 ***
## X2           1.497492   0.005459 274.322   <2e-16 ***
## X1cA2:X2    -0.004174   0.007507  -0.556    0.579    
## X1cA3:X2     0.003408   0.007571   0.450    0.653    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09916 on 294 degrees of freedom
## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9988 
## F-statistic: 4.882e+04 on 5 and 294 DF,  p-value: < 2.2e-16

de donde

\[\begin{align*} &\mathbb{E}(Y;X1c=A_1, X_2)= \beta_0 + \beta_3 X_2\\ &\mathbb{E}(Y;X1c=A 2, X_2)= \beta_0 + \beta_1 + \beta_3 X_2 + \beta_4 X_2 = (\beta_0 + \beta_1) + (\beta_3 + \beta_4) X_2\\ &\mathbb{E}(Y;X1c=A_3, X_2)= \beta_0 + \beta_2 + \beta_3 X_2 + \beta_5 X_2 = (\beta_0 + \beta_2) + (\beta_3 + \beta_5) X_2 \end{align*}\]

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

\[\begin{align*} &H_0: \beta_1=0\ \ y\ \ \beta_2=0\ \ y \ \ \beta_3=0\ \ y \ \ \beta_4=0 \ \ y \ \ \beta_5=0 \ \ vs\\ &H_a: \beta_1\neq0\ \ o\ \ \beta_2\neq0\ \ o\ \ \beta_3\neq0\ \ o\ \ \beta_4\neq0\ \ o\ \ \beta_5\neq0 \end{align*}\]

En la práctica se realiza de forma inmediata una prueba sobre la igualdad de pendientes (efecto de las interacciones), es decir, si el efecto de la variable \(X_2\) en la media de \(Y\) es la misma para los tres niveles de la variable \(X1c\) (se asume que hay un efecto de la variable \(X_2\)).

Pendientes iguales:

\[ H_0: \beta_4=0\ \ y\ \ \beta_5=0 \ \ vs\ \ H_a: \beta_4\neq0\ \ o\ \ \beta_5\neq0 \]

donde se busca no rechazar \(H_{0}\), es decir, que dea plausible quitar las interacciones.

library(multcomp)
K=matrix(c(0,0,0,0,1,0,
           0,0,0,0,0,1), ncol=6, 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.004174
## 2 == 0  0.003408
## 
## Global Test:
##        F DF1 DF2 Pr(>F)
## 1 0.5339   2 294 0.5869

Observar que la función drop1, por ejemplo, sólo da esta prueba en este caso, al considerarse lo más común a proceder cuando hay interacciones

drop1(fit, test = "F")
## Single term deletions
## 
## Model:
## y ~ X1c + X2 + X1c:X2
##        Df Sum of Sq    RSS     AIC F value Pr(>F)
## <none>              2.8908 -1380.7               
## X1c:X2  2  0.010499 2.9013 -1383.6  0.5339 0.5869

Se observa que se NO rechaza \(H_0\), entonces es plausible realizar el modelado sin considerar interacciones (eso ayudará a interpretar el modelo en particular, cuando el interés es concluir que una de las categorías en X1c es superior o inferior). Entonces el modelo reducido queda como

\[ \mathbb{E}(y;x1c, x_2)= \beta_0 + \beta_1 A_2 + \beta_2 A_3 + \beta_3 X_2 \]

y volvemos a realizar el ajuste

fit2=lm(y~X1c+X2, data=Datos)
summary(fit2)
## 
## 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

\[\begin{align*} &\mathbb{E}(Y;X1c=A_1, X_2)= \beta_0 + \beta_3 X_2\\ &\mathbb{E}(Y;X1c=A 2, X_2)= (\beta_0 + \beta_1) + \beta_3 X_2\\ &\mathbb{E}(Y;X1c=A_3, X_2)= (\beta_0 + \beta_2) + \beta_3 X_2 \end{align*}\]

Ahora se revisa la prueba asociada a la tabla ANOVA, en este caso se rechaza \(H_0\). Después se realiza la prueba sobre el parámetro asociado a la variable continua \(X_2\)

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

Se puede hacer directamente a partir de la salida, se rechaza \(H_0\). Ahora se procede a realizar un prueba global sobre si la media de \(Y\) para alguno de los niveles de la variable X1c presenta diferencias considerando que el modelo se corrige (ajusta, condiciona, toma en cuenta) con \(X_2\)

\[\begin{align*} &H_0: \mathbb{E}(Y;X1c=A_1, X_2) = \mathbb{E}(Y;X1c=A_2, X_2) = \mathbb{E}(Y;X1c=A_3, X_2)\\ &H_0: \beta_0 + \bera_3 X_2 = (\beta_0 + \beta_1) + \beta_3 X_2 = (\beta_0 + \beta_2) + \beta_3 X_2\\ &H_0: 0=\beta_1=\beta_2, \ \ \beta_2-\beta_1=0 \ \ \textrm{es redundante} \end{align*}\]

Notar que aquí fijar un valor de X2 o considerar a todos los valores es similar, esto es porque se tienen las misma pendientes.

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(fit2, 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

Con esto se rechaza H0.

Nota. En caso de no rechazar, hasta aquí terminaría el análisis concluyendo que con los datos observados no se cuenta con evidencia para indicar que algún grupo tiene un desempeño diferente (En papers: no se encontraron diferencias significativas entre los grupos)

Cuando se rechaza Ho.

Se procede a graficar el modelo y vendrán las preguntas de los investigadores. Por ejemplo, si algún nivel de X1c da mayores valores promedio en la variable Y en comparación con el resto de niveles. Se busca identificar y concluir. Esto depende de las hipótesis de los investigadores principalmente.

fit$coefficients
##  (Intercept)        X1cA2        X1cA3           X2     X1cA2:X2     X1cA3:X2 
##  0.504355463  1.324838723  0.488498303  1.497491913 -0.004174038  0.003407903
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 = 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)

Por ejemplo: \(A_2\) es mejor que \(A_3\) y \(A_1\) (notar que sería para todo valor de \(X_2\))

\[\begin{align*} &\mathbb{E}(Y;X1c=A_1, X_2) < \mathbb{E}(Y;X1c=A_2, X_2)\ \ y\ \ \mathbb{E}(Y;X1c=A_3, X_2) < \mathbb{E}(Y;X1c=A_2, X_2)\\ &\beta_0 + \beta_3 X_2 < (\beta_0 + \beta_1) + \beta_3 X_2\ \ y\ \ (\beta_0 + \beta_2) + \beta_3 X_2 < (\beta_0 + \beta_1) + \beta_3 X_2\\ &0 < \beta_1 \ \ y \ \ \beta_2 < \beta_1\\ &0 < \beta_1\ \ y \ \ 0 < \beta_1-\beta_2 \end{align*}\]

Aquí $H_0: 0 _1   y   0 _1-_2.

K=matrix(c(0,1,0,0,
           0,1,-1,0 ), ncol=4, nrow=2, byrow=TRUE)
m=c(0,0)
summary(glht(fit2, linfct=K, rhs=m, alternative="greater"))
## 
##   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  1.31309    0.01405   93.44 <1e-10 ***
## 2 <= 0  0.81353    0.01407   57.82 <1e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza \(H_0\). Y sí se concluye que el nivel \(A_2\) es el mejor,pues ambas \(H_{oi}\) se rechazaron. Nuevamente, si los investigadores no tienen alguna idea sobre qué nivel podría diferir, posterior de la prueba asociada a la tabla ANOVA .se procede a hacer comparaciones por pares

# Lo anterior es equivalente a especificar los contrastes individuales
# H0_1: E(Y;X1c=A1, X2) = E(Y;X1c=A2, X2)   H0_1: b1=0 
# H0_2: E(Y;X1c=A1, X2) = E(Y;X1c=A3, X2)   H0_2: b2=0 
# H0_3: E(Y;X1c=A2, X2) = E(Y;X1c=A3, X2)   H0_3: b2-b1=0 

# con #E(Y|X1c=A1, X2)= b0 + b3X2
#E(Y|X1c=A2, X2)= b0 + b1 + b3X2  = (b0 + b1) + b3X2
#E(Y|X1c=A3, X2)= b0 + b2 + b3X2  = (b0 + b2) + b3X2

Notar que aquí se dejan las tres combinaciones lineales obtenidas para cada hipótesis individual.

#Opción 1 (válida sólo con igualdad de pendientes)
summary(glht(fit2, linfct = mcp(X1c = "Tukey")))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = y ~ X1c + X2, data = Datos)
## 
## Linear Hypotheses:
##              Estimate Std. Error t value Pr(>|t|)    
## A2 - A1 == 0  1.31309    0.01405   93.44   <2e-16 ***
## A3 - A1 == 0  0.49956    0.01400   35.68   <2e-16 ***
## A3 - A2 == 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)
K=matrix(c(0,1,0,0,
           0,0,1,0,
           0,-1,1,0), ncol=4, nrow=3, byrow=TRUE)
m=c(0,0,0)
summary(glht(fit2, 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  1.31309    0.01405   93.44   <2e-16 ***
## 2 == 0  0.49956    0.01400   35.68   <2e-16 ***
## 3 == 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)

Aquí se puede concluir que todos los niveles son diferentes, pero no se podría hablar de dirección, aunque la interpretación se puede complementar con lo observado en la gráfica.

ANCOVA caso 2. Sin igualdad de pendientes.

Limpiamos R y cargamos otros datos

rm(list = ls(all.names = TRUE))
gc()
##           used (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 1678065 89.7    2969075 158.6  2969075 158.6
## Vcells 2982763 22.8    8388608  64.0  4931698  37.7
depression <- read.table("depression.txt", header=TRUE )
str(depression)
## 'data.frame':    36 obs. of  5 variables:
##  $ y  : int  56 41 40 28 55 25 46 71 48 63 ...
##  $ age: int  21 23 30 19 28 23 33 67 42 33 ...
##  $ x2 : int  1 0 0 0 1 0 0 0 0 1 ...
##  $ x3 : int  0 1 1 0 0 0 1 0 1 0 ...
##  $ TRT: chr  "A" "B" "B" "C" ...

La variable \(Y\) es una medida de la efectividad del tratamiento para la depresión entre mayor valor de \(Y\) es mejor la efectividad.

# transformamos la variable a tipo factor
depression$TRT=factor(depression$TRT)

# graficamos
plot(depression$age, depression$y, pch=19, col = c("red", "green", "blue")[depression$TRT])
legend("bottomright", levels(depression$TRT),
       col = c("red", "green", "blue"), pch = 19, inset = 0.01,  pt.cex=1.5,cex = .9, y.intersp = 1.3 , bty="n")

levels(depression$TRT)
## [1] "A" "B" "C"

Aquí la referencia es TRT A. Ajuste del modelo con interacciones:

\[ \mathbb{E}(y;x)= \beta_0 + \beta_1 age + \beta_2 TRTB + \beta_3 TRTC + \beta_4(age\cdot TRTB) + \beta_5(age\cdot TRTC) \]

fit <- lm(y ~ age * TRT, data = depression) 
summary(fit)
## 
## Call:
## lm(formula = y ~ age * TRT, data = depression)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4366 -2.7637  0.1887  2.9075  6.5634 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  47.51559    3.82523  12.422 2.34e-13 ***
## age           0.33051    0.08149   4.056 0.000328 ***
## TRTB        -18.59739    5.41573  -3.434 0.001759 ** 
## TRTC        -41.30421    5.08453  -8.124 4.56e-09 ***
## age:TRTB      0.19318    0.11660   1.657 0.108001    
## age:TRTC      0.70288    0.10896   6.451 3.98e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.925 on 30 degrees of freedom
## Multiple R-squared:  0.9143, Adjusted R-squared:  0.9001 
## F-statistic: 64.04 on 5 and 30 DF,  p-value: 4.264e-15

donde

#E(Y;TRT=A, age)= b0 + b1 age
#E(Y;TRT=B, age)=  b0 + b1 age +b2 + b4 age = (b0 + b2) + (b1 + b4) age
#E(Y;TRT=C, age)=  b0 + b1 age +b3 + b5 age = (b0 + b3) + (b1 + b5) age

Se rechaza \(H_o\) en la prueba asociada a la tabla ANOVA, entonces el modelo tiene sentido. Seguimos con prueba de igualdad de pendientes (coeficientes asociados a las interacciones)

\[ H_0: \beta_4=0\ \ y\ \ \beta_5=0\ \ vs\ \ H_a: \beta_4\neq0\ \ o\ \ \beta_5\neq0 \]

K=matrix(c(0,0,0,0,1,0,
           0,0,0,0,0,1), ncol=6, 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.1932
## 2 == 0   0.7029
## 
## Global Test:
##       F DF1 DF2   Pr(>F)
## 1 22.83   2  30 9.41e-07

Se rechaza \(H_0\), al menos hay una pendiente que difiere. Escribimos una prueba de hipótesis simultánea para analizar qué provocó el rechazo anterior. Notar que aquí se incluyen las tres hipótesis individuales asociadas a

# H0_1: pendientes niveles A y B iguales (b4=0)
# H0_2: pendientes niveles A y C iguales (b5=0)
# H0_3: pendientes niveles B y C iguales (b4-b5=0)
#¿Todas las pendientes difieren?
K=matrix(c(0,0,0,0,1,0,
           0,0,0,0,0,1,
           0,0,0,0,1,-1), ncol=6, nrow=3, byrow=TRUE)
m=c(0,0,0)
summary(glht(fit, linfct=K, rhs=m))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ age * TRT, data = depression)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0   0.1932     0.1166   1.657    0.238    
## 2 == 0   0.7029     0.1090   6.451   <0.001 ***
## 3 == 0  -0.5097     0.1104  -4.617   <0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

De forma simultánea podemos identificar que es plausible identificar que \(\beta_4\) no es diferente de 0, esto implica que tratamiento A y B podrían tener las mismas pendientes pero sí se observa que el tratamiento C tiene diferentes pendientes que A y también que B.

Se procede a reducir, es decir, considerar un modelo con menos parámetros al eliminar la interacción (age*TRTB), parámetro \(\beta_4\)

fitred <- lm(y ~ age + TRT + I(age*(TRT=="C")), data = depression) 
summary(fitred)
## 
## Call:
## lm(formula = y ~ age + TRT + I(age * (TRT == "C")), data = depression)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.0517 -2.2394 -0.1463  3.2348  6.4728 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            43.28524    2.92725  14.787 1.37e-15 ***
## age                     0.42486    0.05990   7.093 5.74e-08 ***
## TRTB                  -10.02721    1.64773  -6.085 9.62e-07 ***
## TRTC                  -37.07386    4.51889  -8.204 2.88e-09 ***
## I(age * (TRT == "C"))   0.60853    0.09547   6.374 4.25e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.034 on 31 degrees of freedom
## Multiple R-squared:  0.9065, Adjusted R-squared:  0.8944 
## F-statistic: 75.14 on 4 and 31 DF,  p-value: 1.679e-15

y se vuelven a repetir los pasos.

Con un modelo donde hay alguna interacción aún es posible realizar la prueba sobre si un tratamiento es mejor condicional en un valor de la edad supongamos que se quiere argumentar que el tratamiento A es mejor que el B y C para personas de edad 20

#E(Y;TRT=A, 20)= b0 + b1 20
#E(Y;TRT=B, 20)= (b0 + b2) + (b1 + b4) 20
#E(Y;TRT=C, 20)= (b0 + b3) + (b1 + b5) 20

#E(Y;TRT=A, 20)>E(Y;TRT=B, 20) y E(Y;TRT=A, 20)>E(Y;TRT=C, 20)
#b0 + b1 20 > (b0 + b2) + (b1 + b4) 20  y   b0 + b1 20 >  (b0 + b3) + (b1 + b5) 20
#0 > b2 + b4 20  y  0 >  b3 +  b5 20

Se especifican

\[\begin{align*} H_{o1}: 0 \leq \beta_2 + \beta_4 \cdot 20\\ H_{o1}: 0 \leq \beta_3 + \beta_5 \cdot 20 \end{align*}\]

Notar que aquí sí hay dependencia sobre el valor de la variable continua y no se puede generalizar la prueba (es especifica).

K=matrix(c(0,0,1,0,20,0,
           0,0,0,1,0,20), ncol=6, nrow=2, byrow=TRUE)
m=c(0,0)
summary(glht(fit, linfct=K, rhs=m, alternative="less"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ age * TRT, data = depression)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value   Pr(<t)    
## 1 >= 0  -14.734      3.262  -4.517 8.85e-05 ***
## 2 >= 0  -27.247      3.094  -8.806 8.04e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza \(H_0\), además se puede concluir que ambas \(H_{oi}\) se rechazan, así que A sí es mejor que B y C para edad 20.

#Para edad 55? 
K=matrix(c(0,0,1,0,55,0,
           0,0,0,1,0,55), ncol=6, nrow=2, byrow=TRUE)
m=c(0,0)
summary(glht(fit, linfct=K, rhs=m, alternative="less"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ age * TRT, data = depression)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value   Pr(<t)    
## 1 >= 0   -7.973      2.027  -3.933 0.000446 ***
## 2 >= 0   -2.646      1.984  -1.334 0.161468    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza \(H_0\), pero ahora sólo se puede concluir que A es mejor que B, pero no que C para edad 55.

Veamos ahora una alternativa simple para interpretar con intervalos de confianza simultáneos, para varios valores de edad. Veamos nuevamente la siguiente gráfica

plot(depression$age, depression$y, pch=19, col = c("red", "green", "blue")[depression$TRT])
legend("bottomright", levels(depression$TRT),
       col = c("red", "green", "blue"), pch = 19, inset = 0.01,  pt.cex=1.5,cex = .9, y.intersp = 1.3 , bty="n")

veamos que la varibe age tiene como mínimo el 19 y cómo máximo el 67, así, crearemos un intervalo entre estas edades

summary(depression)
##        y              age              x2               x3         TRT   
##  Min.   :25.00   Min.   :19.00   Min.   :0.0000   Min.   :0.0000   A:12  
##  1st Qu.:46.75   1st Qu.:32.25   1st Qu.:0.0000   1st Qu.:0.0000   B:12  
##  Median :58.00   Median :44.00   Median :0.0000   Median :0.0000   C:12  
##  Mean   :55.17   Mean   :44.11   Mean   :0.3333   Mean   :0.3333         
##  3rd Qu.:63.25   3rd Qu.:56.50   3rd Qu.:1.0000   3rd Qu.:1.0000         
##  Max.   :73.00   Max.   :67.00   Max.   :1.0000   Max.   :1.0000
# Intervalo
age <- seq(from = 19, to = 67, by = .5)

Después, bajaremos la confianza al 90% pues estamos trabajando con intervalos simultáneos. Creamos pues las bands para las rectas

# Para una banda para la recta del tratamiento A
# E(Y;TRT=A, X2)= b0 + b1 age
KA <- cbind(1, age, 0, 0,0,0)
(KA)
##          age        
##  [1,] 1 19.0 0 0 0 0
##  [2,] 1 19.5 0 0 0 0
##  [3,] 1 20.0 0 0 0 0
##  [4,] 1 20.5 0 0 0 0
##  [5,] 1 21.0 0 0 0 0
##  [6,] 1 21.5 0 0 0 0
##  [7,] 1 22.0 0 0 0 0
##  [8,] 1 22.5 0 0 0 0
##  [9,] 1 23.0 0 0 0 0
## [10,] 1 23.5 0 0 0 0
## [11,] 1 24.0 0 0 0 0
## [12,] 1 24.5 0 0 0 0
## [13,] 1 25.0 0 0 0 0
## [14,] 1 25.5 0 0 0 0
## [15,] 1 26.0 0 0 0 0
## [16,] 1 26.5 0 0 0 0
## [17,] 1 27.0 0 0 0 0
## [18,] 1 27.5 0 0 0 0
## [19,] 1 28.0 0 0 0 0
## [20,] 1 28.5 0 0 0 0
## [21,] 1 29.0 0 0 0 0
## [22,] 1 29.5 0 0 0 0
## [23,] 1 30.0 0 0 0 0
## [24,] 1 30.5 0 0 0 0
## [25,] 1 31.0 0 0 0 0
## [26,] 1 31.5 0 0 0 0
## [27,] 1 32.0 0 0 0 0
## [28,] 1 32.5 0 0 0 0
## [29,] 1 33.0 0 0 0 0
## [30,] 1 33.5 0 0 0 0
## [31,] 1 34.0 0 0 0 0
## [32,] 1 34.5 0 0 0 0
## [33,] 1 35.0 0 0 0 0
## [34,] 1 35.5 0 0 0 0
## [35,] 1 36.0 0 0 0 0
## [36,] 1 36.5 0 0 0 0
## [37,] 1 37.0 0 0 0 0
## [38,] 1 37.5 0 0 0 0
## [39,] 1 38.0 0 0 0 0
## [40,] 1 38.5 0 0 0 0
## [41,] 1 39.0 0 0 0 0
## [42,] 1 39.5 0 0 0 0
## [43,] 1 40.0 0 0 0 0
## [44,] 1 40.5 0 0 0 0
## [45,] 1 41.0 0 0 0 0
## [46,] 1 41.5 0 0 0 0
## [47,] 1 42.0 0 0 0 0
## [48,] 1 42.5 0 0 0 0
## [49,] 1 43.0 0 0 0 0
## [50,] 1 43.5 0 0 0 0
## [51,] 1 44.0 0 0 0 0
## [52,] 1 44.5 0 0 0 0
## [53,] 1 45.0 0 0 0 0
## [54,] 1 45.5 0 0 0 0
## [55,] 1 46.0 0 0 0 0
## [56,] 1 46.5 0 0 0 0
## [57,] 1 47.0 0 0 0 0
## [58,] 1 47.5 0 0 0 0
## [59,] 1 48.0 0 0 0 0
## [60,] 1 48.5 0 0 0 0
## [61,] 1 49.0 0 0 0 0
## [62,] 1 49.5 0 0 0 0
## [63,] 1 50.0 0 0 0 0
## [64,] 1 50.5 0 0 0 0
## [65,] 1 51.0 0 0 0 0
## [66,] 1 51.5 0 0 0 0
## [67,] 1 52.0 0 0 0 0
## [68,] 1 52.5 0 0 0 0
## [69,] 1 53.0 0 0 0 0
## [70,] 1 53.5 0 0 0 0
## [71,] 1 54.0 0 0 0 0
## [72,] 1 54.5 0 0 0 0
## [73,] 1 55.0 0 0 0 0
## [74,] 1 55.5 0 0 0 0
## [75,] 1 56.0 0 0 0 0
## [76,] 1 56.5 0 0 0 0
## [77,] 1 57.0 0 0 0 0
## [78,] 1 57.5 0 0 0 0
## [79,] 1 58.0 0 0 0 0
## [80,] 1 58.5 0 0 0 0
## [81,] 1 59.0 0 0 0 0
## [82,] 1 59.5 0 0 0 0
## [83,] 1 60.0 0 0 0 0
## [84,] 1 60.5 0 0 0 0
## [85,] 1 61.0 0 0 0 0
## [86,] 1 61.5 0 0 0 0
## [87,] 1 62.0 0 0 0 0
## [88,] 1 62.5 0 0 0 0
## [89,] 1 63.0 0 0 0 0
## [90,] 1 63.5 0 0 0 0
## [91,] 1 64.0 0 0 0 0
## [92,] 1 64.5 0 0 0 0
## [93,] 1 65.0 0 0 0 0
## [94,] 1 65.5 0 0 0 0
## [95,] 1 66.0 0 0 0 0
## [96,] 1 66.5 0 0 0 0
## [97,] 1 67.0 0 0 0 0

continuando

# Para una banda para la recta del tratamiento B
# E(Y;TRT=B, X2)=  b0 + b1 age +b2 + b4 age = (b0 + b2) + (b1 + b4) age
KB <- cbind(1, age, 1, 0, age, 0)

#E(Y;TRT=C, X2)=  b0 + b1 age +b3 + b5 age = (b0 + b3) + (b1 + b5) age
KC <- cbind(1, age, 0, 1, 0, age)

# Juntamos las bandas
K = rbind(KA, KB, KC)

Intervalos simultáneos:

library(multcomp)

fitE <- glht(fit, linfct = K)
fitci <- confint(fitE, level = 0.90)

y finalmente graficamos

plot(depression$age, depression$y, pch=19, col = c("red", "green", "blue")[depression$TRT])
legend("bottomright", levels(depression$TRT),
       col = c("red", "green", "blue"), pch = 19, inset = 0.01,  pt.cex=1.5,cex = .9, y.intersp = 1.3 , bty="n")

lines(age, coef(fitE)[1:97], col="red")
lines(age, fitci$confint[1:97,"upr"], col="red")
lines(age, fitci$confint[1:97,"lwr"], col="red")

lines(age, coef(fitE)[98:194], col="green")
lines(age, fitci$confint[98:194,"upr"], col="green")
lines(age, fitci$confint[98:194,"lwr"], col="green")

lines(age, coef(fitE)[195:291], col="blue")
lines(age, fitci$confint[195:291,"upr"], col="blue")
lines(age, fitci$confint[195:291,"lwr"], col="blue")

Estas bandas de confianza o intervalos de confianza simultáneos sirven para identificar las diferencias más evidentes con los datos, pero son muy conservadoras.

Por ejemplo, el tratamiento A es mejor a los otros tratamientos entre 19 y 49.5 años. En efecto:

cbind(age, fitci$confint[1:97,"lwr"], fitci$confint[195:291,"upr"])
##     age                  
## 1  19.0 47.35663 31.52382
## 2  19.5 47.61830 31.95862
## 3  20.0 47.87953 32.39394
## 4  20.5 48.14029 32.82977
## 5  21.0 48.40058 33.26616
## 6  21.5 48.66036 33.70313
## 7  22.0 48.91962 34.14069
## 8  22.5 49.17832 34.57889
## 9  23.0 49.43643 35.01775
## 10 23.5 49.69393 35.45729
## 11 24.0 49.95080 35.89756
## 12 24.5 50.20698 36.33859
## 13 25.0 50.46246 36.78042
## 14 25.5 50.71719 37.22307
## 15 26.0 50.97113 37.66660
## 16 26.5 51.22425 38.11104
## 17 27.0 51.47649 38.55644
## 18 27.5 51.72782 39.00285
## 19 28.0 51.97818 39.45032
## 20 28.5 52.22753 39.89889
## 21 29.0 52.47579 40.34863
## 22 29.5 52.72293 40.79958
## 23 30.0 52.96886 41.25182
## 24 30.5 53.21354 41.70539
## 25 31.0 53.45688 42.16037
## 26 31.5 53.69881 42.61682
## 27 32.0 53.93927 43.07482
## 28 32.5 54.17815 43.53443
## 29 33.0 54.41538 43.99574
## 30 33.5 54.65087 44.45881
## 31 34.0 54.88453 44.92373
## 32 34.5 55.11624 45.39058
## 33 35.0 55.34591 45.85945
## 34 35.5 55.57344 46.33041
## 35 36.0 55.79871 46.80355
## 36 36.5 56.02160 47.27896
## 37 37.0 56.24200 47.75672
## 38 37.5 56.45980 48.23691
## 39 38.0 56.67486 48.71963
## 40 38.5 56.88707 49.20494
## 41 39.0 57.09631 49.69292
## 42 39.5 57.30245 50.18366
## 43 40.0 57.50538 50.67721
## 44 40.5 57.70498 51.17363
## 45 41.0 57.90115 51.67300
## 46 41.5 58.09379 52.17535
## 47 42.0 58.28280 52.68072
## 48 42.5 58.46810 53.18916
## 49 43.0 58.64961 53.70068
## 50 43.5 58.82728 54.21531
## 51 44.0 59.00105 54.73304
## 52 44.5 59.17091 55.25387
## 53 45.0 59.33682 55.77780
## 54 45.5 59.49879 56.30480
## 55 46.0 59.65682 56.83484
## 56 46.5 59.81095 57.36787
## 57 47.0 59.96122 57.90386
## 58 47.5 60.10769 58.44275
## 59 48.0 60.25041 58.98448
## 60 48.5 60.38948 59.52897
## 61 49.0 60.52499 60.07617
## 62 49.5 60.65702 60.62599
## 63 50.0 60.78570 61.17836
## 64 50.5 60.91112 61.73319
## 65 51.0 61.03341 62.29039
## 66 51.5 61.15268 62.84990
## 67 52.0 61.26907 63.41161
## 68 52.5 61.38268 63.97545
## 69 53.0 61.49364 64.54133
## 70 53.5 61.60207 65.10916
## 71 54.0 61.70809 65.67887
## 72 54.5 61.81181 66.25038
## 73 55.0 61.91334 66.82360
## 74 55.5 62.01281 67.39846
## 75 56.0 62.11030 67.97490
## 76 56.5 62.20592 68.55282
## 77 57.0 62.29976 69.13218
## 78 57.5 62.39192 69.71289
## 79 58.0 62.48249 70.29491
## 80 58.5 62.57155 70.87817
## 81 59.0 62.65917 71.46260
## 82 59.5 62.74544 72.04816
## 83 60.0 62.83043 72.63480
## 84 60.5 62.91419 73.22246
## 85 61.0 62.99680 73.81109
## 86 61.5 63.07832 74.40065
## 87 62.0 63.15880 74.99110
## 88 62.5 63.23830 75.58240
## 89 63.0 63.31686 76.17450
## 90 63.5 63.39453 76.76737
## 91 64.0 63.47137 77.36098
## 92 64.5 63.54740 77.95529
## 93 65.0 63.62268 78.55026
## 94 65.5 63.69723 79.14588
## 95 66.0 63.77110 79.74211
## 96 66.5 63.84432 80.33892
## 97 67.0 63.91691 80.93629