Ejemplo de estimación. Problema tipo Anova

Datos:

Los datos corresponden a un experimento sobre el crecimiento de plantas. El objetivo es analizar si el crecimiento de las plantas se diferencía de acuerdo con la condición de crecimiento que se establece con dos diferentes tratamientos y el caso donde no se aplica nada (comúnmente llamado tratamiento control). La selección de las plantas a las que se aplicaron los tratamientos se hace de forma aleatoria.

Las preguntas de investigación son:

  1. ¿Existe alguna condición (tratamiento) que favorece el crecimiento de las plantas?
  2. ¿Se puede concluir que existe una condición que favorece el crecimiento de las plantas en comparación con el resto?

Manejo de los datos

Comenzamos por cargar y revisar los datos

data(PlantGrowth)

summary(PlantGrowth)
##      weight       group   
##  Min.   :3.590   ctrl:10  
##  1st Qu.:4.550   trt1:10  
##  Median :5.155   trt2:10  
##  Mean   :5.073            
##  3rd Qu.:5.530            
##  Max.   :6.310

veamos que la variable group ya viene como del tipo factor con tres niveles

is.factor(PlantGrowth$group)
## [1] TRUE
levels(PlantGrowth$group)
## [1] "ctrl" "trt1" "trt2"

Visualización 

library(GGally)
ggpairs(PlantGrowth[,c(1,2)])

del cual puede verse de la segunda gráfica de la primer fila que el tratamiento 2 parece ser el efectivo versus el tratamiento 1. Recordemos que en estos casos es importante centrarnos en el análisis del boxplot

boxplot(weight ~ group, data = PlantGrowth, col = "white", outline=FALSE)
stripchart(weight ~ group, data = PlantGrowth,
           method = "jitter",
           pch = 19,
           col = 2:4,
           vertical = TRUE,
           add = TRUE)

Donde visualmente podemos ver que los supuestos, como se analizaron en regresión simple, se están cumpliendo. Por lo general, en este tipo de problemas se debe revisar si la asignación de los tratamientos se hizo de forma aleatorizada (lo cual por lo general se debe revisar con los investigadores).

Además, para ser tratado como un problema de regresión lineal múltiple el supuesto de homocedasticidad se debe cumplir para cada grupo, lo mismo que la normalidad.

Ejemplo para verificación de supuestos (similar a caso de reg. lineal simple) 

# Homocedasticidad
bartlett.test(weight ~ group, data = PlantGrowth)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  weight by group
## Bartlett's K-squared = 2.8786, df = 2, p-value = 0.2371
library(car)
leveneTest(weight ~ group, data = PlantGrowth)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.1192 0.3412
##       27
# Normalidad por grupos

shapiro.test( PlantGrowth$weight[PlantGrowth$group==levels(PlantGrowth$group)[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  PlantGrowth$weight[PlantGrowth$group == levels(PlantGrowth$group)[1]]
## W = 0.95668, p-value = 0.7475
shapiro.test( PlantGrowth$weight[PlantGrowth$group==levels(PlantGrowth$group)[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  PlantGrowth$weight[PlantGrowth$group == levels(PlantGrowth$group)[2]]
## W = 0.93041, p-value = 0.4519
shapiro.test( PlantGrowth$weight[PlantGrowth$group==levels(PlantGrowth$group)[3]])
## 
##  Shapiro-Wilk normality test
## 
## data:  PlantGrowth$weight[PlantGrowth$group == levels(PlantGrowth$group)[3]]
## W = 0.94101, p-value = 0.5643

Si no se cumplieran los supuestos habría que realizar algún tipo de transformación de los datos.

Uso de una regresión lineal múltiple 

fit1=lm(weight ~ group, data = PlantGrowth)
summary(fit1)
## 
## Call:
## lm(formula = weight ~ group, data = PlantGrowth)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0710 -0.4180 -0.0060  0.2627  1.3690 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.0320     0.1971  25.527   <2e-16 ***
## grouptrt1    -0.3710     0.2788  -1.331   0.1944    
## grouptrt2     0.4940     0.2788   1.772   0.0877 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6234 on 27 degrees of freedom
## Multiple R-squared:  0.2641, Adjusted R-squared:  0.2096 
## F-statistic: 4.846 on 2 and 27 DF,  p-value: 0.01591

De modo que el modelo quedaría como \(\mathbb{E}(Y;grupo)=\beta_0 + \beta_1 (\textrm{grouptrt}1) + \beta_2 (\textrm{grouptrt}2)\), donde

\[\begin{align*} &\mathbb{E}(Y|\textrm{grupo}=\textrm{control})=\beta_{0}\\ &\mathbb{E}(Y|\textrm{grupo}=\textrm{trt1})=\beta_{0}+\beta_{1}\\ &\mathbb{E}(Y|\textrm{grupo}=\textrm{trt2})=\beta_{0}+\beta_{2} \end{align*}\]

¿Dichas medias son iguales?

Son iguales si \(\beta_1=0\) y \(\beta_2=0\). Prueba de la tabla ANOVA justo sirve para hacer ese test

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

#Otra forma de hacer la prueba
drop1(fit1, test="F")
## Single term deletions
## 
## Model:
## weight ~ group
##        Df Sum of Sq    RSS     AIC F value  Pr(>F)  
## <none>              10.492 -25.517                  
## group   2    3.7663 14.258 -20.316  4.8461 0.01591 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

donde se rechaza \(H_{0}\).

#Prueba ANOVA usando la prueba lineal general (la global sólo acepta igualdad)
library(multcomp)
K=matrix(c(0,1,0,
           0,0,1), ncol=3, nrow=2, byrow=TRUE)
m=c(0,0)
summary(glht(fit1, linfct=K, rhs=m),test=Ftest())
## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0   -0.371
## 2 == 0    0.494
## 
## Global Test:
##       F DF1 DF2  Pr(>F)
## 1 4.846   2  27 0.01591

nuevamente se rechaza \(H_{0}\).

La prueba asociada a la tabla anova nos indica que se rechaza \(H_0\), es decir, al menos uno de los grupos tiene una peso promedio diferente en el peso de las plantas.

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ñoo diferente (en papers: no se encontraron diferencias significativas entre los grupos)

Pruebas para identificar diferencias y dirección

La dirección estricta y diferencias se ponen en la hipótesis alternativa.

Trat 1 mejor que trat 2

Lo que se quiere se pone en la hipótesis alternativa

\[ \mathbb{E}(Y;\textrm{grupo}=trt1)=\beta_0+\beta_1 >\mathbb{E}(Y;\textrm{grupo}=trt2)=\beta_0+\beta_2 \]

equivalentemente \(\beta_1>\beta_2\) ó \(\beta_1-\beta_2>0\). Aquí \(H_0: \beta_1-\beta_2\leq 0\).

K=matrix(c(0,1,-1), ncol=3, nrow=1, byrow=TRUE)
m=c(0)
summary(glht(fit1, linfct=K, rhs=m, alternative="greater"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = weight ~ group, data = PlantGrowth)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)
## 1 <= 0  -0.8650     0.2788  -3.103  0.998
## (Adjusted p values reported -- single-step method)

donde no se rechaza \(H_{0}\), es deicr, es plausible asumir que \(\beta_{1}\leq \beta_{2}\).

Trat 2 mejor que trat 1

En este caso se considera \(H_{a}: \beta_{1}-\beta_{2}<0\) y \(H_{0}: \beta_{1}-\beta_{2}\geq 0\)

K=matrix(c(0,1,-1), ncol=3, nrow=1, byrow=TRUE)
m=c(0)
summary(glht(fit1, linfct=K, rhs=m, alternative="less"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = weight ~ group, data = PlantGrowth)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value  Pr(<t)   
## 1 >= 0  -0.8650     0.2788  -3.103 0.00223 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
# -----

#equivalentemente podría ser b2-b1>0
#Aquí Ho: b2-b1<=0

# K=matrix(c(0,-1,1), ncol=3, nrow=1, byrow=TRUE)
# m=c(0)
# summary(glht(fit1, linfct=K, rhs=m, alternative="greater"))

rechazándose \(H_{0}\) y por ende se considera \(H_{a}: \beta_{1}-\beta_{2}<0\) ó \(\beta_{1}<\beta_{2}\).

Lo natural después de la conclusión anterior es lo siguiente

Trat 2 mejor que control

\(\mathbb{E}(Y;\textrm{grupo}=\textrm{control})=\beta_0<\mathbb{E}(Y;\textrm{grupo}=trt2)= \beta_0+\beta_2\), de donde \(\beta_0<\beta_0+\beta_2\), o equivalentemente \(0<\beta_2\). Aquí Ho: \(0\geq \beta_2\).

K=matrix(c(0,0,1), ncol=3, nrow=1, byrow=TRUE)
m=c(0)
summary(glht(fit1, linfct=K, rhs=m, alternative="greater"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = weight ~ group, data = PlantGrowth)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)  
## 1 <= 0   0.4940     0.2788   1.772 0.0438 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

de modo que el tratamiento 2 es mejor que el control. Finalmente, para ver si el trat2 es el mejor trataiento pasamos a lo siguiente

El trat 2 mejor que control y mejor que trat 1. Aquí hacemos prueba simultánea

  1. Trat 2 mejor que control: \(0<\beta_2\)

  2. Trat 2 mejor que trat 1: \(\beta_2-\beta_1>0\)

Revisar que la alternativa tenga la misma dirección. Aquí \(H_0: 0\geq \beta_2\) y \(\beta_2-\beta_1\leq0\).

K=matrix(c(0,0,1,
           0,-1,1 ), ncol=3, nrow=2, byrow=TRUE)
m=c(0,0)
summary(glht(fit1, linfct=K, rhs=m, alternative="greater"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = weight ~ group, data = PlantGrowth)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value  Pr(>t)   
## 1 <= 0   0.4940     0.2788   1.772 0.07684 . 
## 2 <= 0   0.8650     0.2788   3.103 0.00423 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Para concluir que se rechaza \(H_0\), se debe rechazar con cierta significancia al menos una de las pruebas, aunque la interpretación deseada es sobre el cumplimiento de las dos condiciones en la alternativa entonces buscamos que se rechacen ambas. Con base en lo último no se puede conluir con \(\alpha=.05\) que el tratamiento 2 es mejor que ambos, el tratamiento 1 y el control.

Se puede concluir que el tratamiento 2 es mejor que 1, pero no necesariamente que el control. Sin embargo con \(\alpha=.1\) sí es posible concluir eso.

En la práctica cuando no se conoce nada de antemano, o el objetivo sólo es identificar qué grupo es donde se observan las diferencias, se procede a comparar todos los pares con \(H_{0i}: \textrm{igualdad de dos grupos vs} \ \ H_{ai}: \textrm{diferencia de dos grupos}\)

# todas las comparaciones por pares de tipo Tukey
summary(glht(fit1, linfct = mcp(group = "Tukey")))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = weight ~ group, data = PlantGrowth)
## 
## Linear Hypotheses:
##                  Estimate Std. Error t value Pr(>|t|)  
## trt1 - ctrl == 0  -0.3710     0.2788  -1.331   0.3909  
## trt2 - ctrl == 0   0.4940     0.2788   1.772   0.1979  
## trt2 - trt1 == 0   0.8650     0.2788   3.103   0.0119 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se concluye que sólo se podría considerar la diferencia de forma simultánea entre tratamiento 1 y tratamiento 2, sin embargo, este procedimiento no incluye dirección además que incluye tres combinaciones lineales al mismo tiempo, lo que hace que sólo se identifiquen los casos más evidentes en donde hay diferencias.

Recomendación

Hacer primero la prueba \(F\) asociada a la tabla ANova, después hacer la prueba por pares tukey.

Parte asociada a la sección 2.10.3.2: Versión A) sección 2.10.3.2

Sin redundancias:

K=matrix(c(0,1,0,
           0,0,1), ncol=3, nrow=2, byrow=TRUE)
m=c(0,0)
summary(glht(fit1, linfct=K, rhs=m),test=Ftest())
## 
##   General Linear Hypotheses
## 
## Linear Hypotheses:
##        Estimate
## 1 == 0   -0.371
## 2 == 0    0.494
## 
## Global Test:
##       F DF1 DF2  Pr(>F)
## 1 4.846   2  27 0.01591

Para este ejemplo es la prueba F asociada a la tabla anova

Versión B) sección 2.10.3.2

Con redundancias:

K=matrix(c(0,1,0,
           0,0,1,
           0,-1,1), ncol=3, nrow=3, byrow=TRUE)
m=c(0,0,0)
summary(glht(fit1, linfct=K, rhs=m))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = weight ~ group, data = PlantGrowth)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)  
## 1 == 0  -0.3710     0.2788  -1.331    0.391  
## 2 == 0   0.4940     0.2788   1.772    0.198  
## 3 == 0   0.8650     0.2788   3.103    0.012 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Para este ejemplo, son todas las pruebas por pares tipo “Tukey”