Ejemplo Abeto

Abeto blanco (Abes alba), es un árbol de gran belleza por la elegancia de sus formas y el exquisito perfume balsámico que destilan sus hojas y cortezas. Destilando hojas y madera se obtiene aceite de trementina muy utilizado en medicina contra torceduras y contusiones. En estos últimos años se ha observado que la producción de semillas ha descendido y con objeto de conseguir buenas producciones se proponen tres tratamientos. Se observa que árboles diferentes tienen distintas características naturales de reproducción, este efecto de las diferencias entre los árboles se debe de controlar y este control se realiza mediante bloques. En el experimento se utilizan 10 abetos, dentro de cada abeto se seleccionan tres ramas semejantes. Cada rama recibe exactamente uno de los tres tratamientos que son asignados aleatoriamente. Constituyendo cada árbol un bloque completo. Los datos obtenidos se presentan en archivo PracticaDBCA.

Nos interesa saber si los distintos tratamientos influyen en la producción de semillas, para ello debemos contrastar la hipótesis de que no hay diferencia en las medias de los tres tratamientos frente a la alternativa de que al menos una media difiere de otra

Tabla ANOVA

Para probar las hipótesis dada utilizaremoas ANOVA para un diseño en bloques completos al azar (DBCA), en el cual SE descompone la variabilidad total de los datos en sus tres componentes: la variabilidad debida a tratamientos, la debida al factor de bloque y la que corresponde al error aleatorio.

library(readxl)
PracticaDBCA <- read_excel("Practica_diseño _Bloques.xlsx", 
    sheet = "datos") # leer los datos
PracticaDBCA

continuación debemos transformar tanto la columna de los tratamientos como la de los bloques en un factor para podemos realizar los cálculos posteriores adecuadamente.

PracticaDBCA$Tratamiento = factor(PracticaDBCA$Tratamiento)
PracticaDBCA$Tratamiento
##  [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## Levels: 1 2 3
PracticaDBCA$Abeto = factor(PracticaDBCA$Abeto)
PracticaDBCA$Abeto
##  [1] 1  1  1  2  2  2  3  3  3  4  4  4  5  5  5  6  6  6  7  7  7  8  8  8  9 
## [26] 9  9  10 10 10
## Levels: 1 2 3 4 5 6 7 8 9 10

Para calcular la tabla ANOVA primero hacemos uso de la función aov de la siguiente forma:

# obtencion de la tabla ANOVA
fit <- aov(y~Tratamiento+Abeto, data = PracticaDBCA)
summary(fit)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Tratamiento  2   16.2   8.100   9.228 0.00174 ** 
## Abeto        9   54.8   6.089   6.937 0.00026 ***
## Residuals   18   15.8   0.878                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

y es el nombre de la columna de las observaciones (variable respuesta)

Tratamiento es el nombre de la columna en la que están representados los tratamientos (Factor principal)

Abeto es el nombre de la columna en la que están representados los bloques data = data.frame en el que están guardados los datos

La Tabla ANOVA, muestra que:

El valor del estadístico de contraste de igualdad de tratamiento, F = 9.228 deja a su derecha un p-valor de 0.00174, menor que el nivel de significación del 1%.

En cuanto a los bloques bloques, F = 6.937 deja a su derecha un p-valor menor que 0.00026, menor que el nivel de significación del 1%, por lo que se podria decir que la eficacia de este diseño depende de los efectos de los bloques. Un valor grande de F de los bloques (6.937) implica que el factor bloque tiene un efecto grande. En este caso el diseño es más eficaz que el diseño completamente aleatorizado ya que si el cuadrado medio entre bloques es grande (6.089), el término residual será mucho menor (0.878) y el contraste principal de las medias de los tratamientos será más sensible a las diferencias entre tratamientos. Por lo tanto la inclusión del factor bloque en el modelo es acertada. Así, la producción de semillas depende del abeto.

El modelo que hemos propuesto hay que validarlo, para ello hay que comprobar si se verifican los supuestos

Verificación de supuestos

Hipótesis de aditividad entre los bloques y tratamientos

La interacción entre el factor bloque y los tratamientos vamos a estudiarla analíticamente mediante el Test de Interacción de un grado de Tukey

Para realizar este test en R tenemos que utilizar la library “daewr” y dentro de ella la función Tukey1df. De la siguiente forma:

library(daewr)
## Registered S3 method overwritten by 'DoE.base':
##   method           from       
##   factorize.factor conf.design
Tukey1df(data.frame(PracticaDBCA))
## Source           df     SS        MS        F     Pr>F 
## A                 2   16.2    8.1 
## B                 9   54.8    6.0889 
## Error            18   15.8    71.1 
## NonAdditivity     1   3.5573    3.5573    4.94    0.0401 
## Residual         17   12.2427    0.7202

Analizando el tukey test podemos observar que no hay interacción entre los tratamientos aplicados y los abetos.

Supuesto de normalidad de residuos

Para verificar la normalidad de los residuos utilizaremos la prueba de Shapiro-Wilks cuyo script es el siguiente:

shapiro.test(residuals(fit))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(fit)
## W = 0.96415, p-value = 0.3935

Con el resultado anterior con una confianza del 95 % no podemos rechazar la idea de que los residuos provengan de una distribución normal \((W=0.9688, p-value > 0.05)\)

de foma gráfica este mismo supuesto

qqnorm(fit$residuals)

Homocedasticidad de Varianzas

Para verificar el supuesto de homocedasticidad de las varianzas utilizaremos la prueba de Bartlett script es el siguiente:

bartlett.test(PracticaDBCA$y, PracticaDBCA$Tratamiento)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  PracticaDBCA$y and PracticaDBCA$Tratamiento
## Bartlett's K-squared = 4.1729, df = 2, p-value = 0.1241

Con una confianza del 95 % no podemos rechazar la idea de que las varianzas sean homogeneas \((chi^2= 4.1729, df = 2, p-value = 0.1241)\)

De la misma manera procedemos para el factor bloque:

bartlett.test(PracticaDBCA$y, PracticaDBCA$Abeto)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  PracticaDBCA$y and PracticaDBCA$Abeto
## Bartlett's K-squared = 4.0723, df = 9, p-value = 0.9066

Con una confianza del 95 % no podemos rechazar la idea de que las varianzas sean homogeneas \((chi^2= 4.0723, df = 9, p-value = 0.9066)\)

plot(fit)

Independencia

plot(fit)

Nos fijamos en el primer gráfico que representa los residuos frente a los valores ajustados y observamos que no hay ninguna tendencia sistemática. El grafico no muestra evidencia de falta de independencia Concluimos que no hay sospechas para que se incumpla la hipótesis de independencia

Dado que se han verificado los supuestos del modelo ahora podriamos concluir que con una confianza del 99% rechamos la hipótesis nula de igualdad de tratamientos. Así, los tratamientos influyen en el número de semillas. Es decir, existen diferencias significativas en el número de semillas entre los tres tratamientos.

Para averiguar cuál de los tratamientos es el que difiere procedemos a realizar un test de rangos multiples.

Test de rangos multiples

Test LSD Tratamiento

library(agricolae)
 LSD.test(fit, "Tratamiento", console = TRUE, group= TRUE, alpha = 0.05)
## 
## Study: fit ~ "Tratamiento"
## 
## LSD t Test for y 
## 
## Mean Square Error:  0.8777778 
## 
## Tratamiento,  means and individual ( 95 %) CI
## 
##      y      std  r      LCL       UCL Min Max
## 1  8.3 1.337494 10 7.677553  8.922447   7  11
## 2  9.2 1.135292 10 8.577553  9.822447   8  12
## 3 10.1 2.183270 10 9.477553 10.722447   7  14
## 
## Alpha: 0.05 ; DF Error: 18
## Critical Value of t: 2.100922 
## 
## least Significant Difference: 0.8802727 
## 
## Treatments with the same letter are not significantly different.
## 
##      y groups
## 3 10.1      a
## 2  9.2      b
## 1  8.3      c
 LSD.test(fit, "Tratamiento", console = TRUE, group= FALSE, alpha = 0.05)
## 
## Study: fit ~ "Tratamiento"
## 
## LSD t Test for y 
## 
## Mean Square Error:  0.8777778 
## 
## Tratamiento,  means and individual ( 95 %) CI
## 
##      y      std  r      LCL       UCL Min Max
## 1  8.3 1.337494 10 7.677553  8.922447   7  11
## 2  9.2 1.135292 10 8.577553  9.822447   8  12
## 3 10.1 2.183270 10 9.477553 10.722447   7  14
## 
## Alpha: 0.05 ; DF Error: 18
## Critical Value of t: 2.100922 
## 
## Comparison between treatments means
## 
##       difference pvalue signif.       LCL         UCL
## 1 - 2       -0.9 0.0456       * -1.780273 -0.01972731
## 1 - 3       -1.8 0.0004     *** -2.680273 -0.91972731
## 2 - 3       -0.9 0.0456       * -1.780273 -0.01972731

Test HSD Tratamiento

library(agricolae)
HSD.test(fit, "Tratamiento", console = TRUE, group= TRUE, alpha = 0.05)
## 
## Study: fit ~ "Tratamiento"
## 
## HSD Test for y 
## 
## Mean Square Error:  0.8777778 
## 
## Tratamiento,  means
## 
##      y      std  r Min Max
## 1  8.3 1.337494 10   7  11
## 2  9.2 1.135292 10   8  12
## 3 10.1 2.183270 10   7  14
## 
## Alpha: 0.05 ; DF Error: 18 
## Critical Value of Studentized Range: 3.609304 
## 
## Minimun Significant Difference: 1.06934 
## 
## Treatments with the same letter are not significantly different.
## 
##      y groups
## 3 10.1      a
## 2  9.2     ab
## 1  8.3      b
HSD.test(fit, "Tratamiento", console = TRUE, group= FALSE, alpha = 0.05)
## 
## Study: fit ~ "Tratamiento"
## 
## HSD Test for y 
## 
## Mean Square Error:  0.8777778 
## 
## Tratamiento,  means
## 
##      y      std  r Min Max
## 1  8.3 1.337494 10   7  11
## 2  9.2 1.135292 10   8  12
## 3 10.1 2.183270 10   7  14
## 
## Alpha: 0.05 ; DF Error: 18 
## Critical Value of Studentized Range: 3.609304 
## 
## Comparison between treatments means
## 
##       difference pvalue signif.      LCL        UCL
## 1 - 2       -0.9 0.1081         -1.96934  0.1693398
## 1 - 3       -1.8 0.0012      ** -2.86934 -0.7306602
## 2 - 3       -0.9 0.1081         -1.96934  0.1693398

Test Duncan Tratamiento

duncan <- duncan.test(fit, "Tratamiento" , group = F)
duncan
## $statistics
##     MSerror Df Mean       CV
##   0.8777778 18  9.2 10.18367
## 
## $parameters
##     test      name.t ntr alpha
##   Duncan Tratamiento   3  0.05
## 
## $duncan
##      Table CriticalRange
## 2 2.971152     0.8802727
## 3 3.117384     0.9235973
## 
## $means
##      y      std  r Min Max  Q25 Q50   Q75
## 1  8.3 1.337494 10   7  11 7.25   8  8.75
## 2  9.2 1.135292 10   8  12 9.00   9  9.00
## 3 10.1 2.183270 10   7  14 9.25  10 11.50
## 
## $comparison
##       difference pvalue signif.       LCL         UCL
## 1 - 2       -0.9 0.0456       * -1.780273 -0.01972735
## 1 - 3       -1.8 0.0006     *** -2.723597 -0.87640270
## 2 - 3       -0.9 0.0456       * -1.780273 -0.01972735
## 
## $groups
## NULL
## 
## attr(,"class")
## [1] "group"