Ejemplo de Diseño en Parcelas Divididas
Sandra Katteryne Rodríguez Hurtado
22/10/2020
Diseño de parcela dividida en R
El siguiente ejercicio fue tomado de la página: http://www.personal.psu.edu/mar36/stat_461/split_plot/split_plot.html.
El diseño consta de bloques (o trazados completos) en los que se aplica un factor (todo el factor de trazado) al azar. Dentro de cada trazado/bloque completo, se divide en unidades más pequeñas y los niveles de segundo factor se aplican aleatoriamente a las piezas más pequeñas de toda la gráfica. La clave es que la unidad experimental es diferente para cada factor.
Ejemplo de avena
En este ejemplo, campos enteros se plantan con uno de los tres tipos de avena. Estas son todas las parcelas, de las cuales hay 18. Cada parcela completa se divide en cuatro parcelas divididas, cada una de las cuales se le asignó aleatoriamente uno de los cuatro niveles de nitrógeno. Se determinó el rendimiento (en bushels/acre) de avena. De interés son las diferencias en la variedad de avena y los niveles de nitrógeno.
sp.oats <- read.csv("D:/Escritorio/oats.csv")
sp.oats <- within (sp.oats, nitroF <- factor(nitro))
head(sp.oats)## observation replicate variety nitro yield nitroF
## 1 1 I Victory 0.0cwt 111 0.0cwt
## 2 2 I Victory 0.2cwt 130 0.2cwt
## 3 3 I Victory 0.4cwt 157 0.4cwt
## 4 4 I Victory 0.6cwt 174 0.6cwt
## 5 5 I Golden.rain 0.0cwt 117 0.0cwt
## 6 6 I Golden.rain 0.2cwt 114 0.2cwtV: Variedad, B: Replicate, X: Observación
Cargar las librerías
library(lattice)
library(car)## Warning: package 'car' was built under R version 4.0.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.0.3library(agricolae)## Warning: package 'agricolae' was built under R version 4.0.3library(gplots)## Warning: package 'gplots' was built under R version 4.0.3##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowesslibrary(readr)with(sp.oats, xyplot(yield ~ nitroF | variety))
with(sp.oats, xyplot(yield ~ nitroF | variety, groups = replicate))
with(sp.oats, xyplot(yield ~ nitroF | variety, groups = replicate, aspect = "xy"))
with(sp.oats, xyplot(yield ~ nitroF | variety, groups = replicate, aspect = "xy",
type = "o"))
with(sp.oats, xyplot(yield ~ nitroF | variety, groups = replicate, aspect = "xy",
type = "a"))
Si ignora el diseño experimental, obtendrá los siguientes resultados incorrectos. Observe, el error df es mayor, por lo que es más fácil detectar diferencias que no están realmente allí.
res.bad <- lm(yield ~ variety * nitroF, data = sp.oats)
anova(res.bad)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## variety 2 1786.4 893.2 1.7949 0.1750
## nitroF 3 20020.5 6673.5 13.4108 8.367e-07 ***
## variety:nitroF 6 321.7 53.6 0.1078 0.9952
## Residuals 60 29857.3 497.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El análisis en R sigue el mismo patrón que la semana pasada. El factor “variedad” sólo tiene (18-1) df total para usar cuando se prueban las diferencias en la variedad. Por lo tanto, necesitamos especificar correctamente el término de error para la variedad. Cabe señalar que la forma de la “Error(A:wholeplot)” puede cambiar dependiendo de los datos que se establecen. La clave es conocer el df correcto para que sepas tener los resultados correctos.
res.good <- aov(yield ~ variety * nitroF + Error(replicate:variety), data = sp.oats)## Warning in aov(yield ~ variety * nitroF + Error(replicate:variety), data =
## sp.oats): Error() model is singularsummary(res.good)##
## Error: replicate:variety
## Df Sum Sq Mean Sq F value Pr(>F)
## variety 2 1786 893.2 0.612 0.555
## Residuals 15 21889 1459.2
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## nitroF 3 20020 6673 37.686 2.46e-12 ***
## variety:nitroF 6 322 54 0.303 0.932
## Residuals 45 7969 177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Para marcar las suposiciones, usted necesita no utilizar el término de error. Puede agregar el término sin error, pero las pruebas F son incorrectas. Sin embargo, la comprobación de la suposición está bien.
res.good2 <- aov(yield ~ variety * nitroF + replicate:variety, data = sp.oats)
summary(res.good2)## Df Sum Sq Mean Sq F value Pr(>F)
## variety 2 1786 893 5.044 0.0106 *
## nitroF 3 20020 6673 37.686 2.46e-12 ***
## variety:nitroF 6 322 54 0.303 0.9322
## variety:replicate 15 21889 1459 8.240 1.61e-08 ***
## Residuals 45 7969 177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(res.good2, 1)
plot(res.good2, 2)
boxCox(res.good2)
Dado que la trama de Box-Cox es buena, no hay necesidad de hacer ninguna transformación. Todavía necesita identificar qué niveles de nitrógeno son diferentes.
(with(sp.oats, HSD.test(yield, nitroF, DFerror = 45, MSerror = 117)))## $statistics
## MSerror Df Mean CV MSD
## 117 45 103.9722 10.40341 9.618527
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey nitroF 4 3.772697 0.05
##
## $means
## yield std r Min Max Q25 Q50 Q75
## 0.0cwt 79.38889 19.39417 18 53 117 63.25 72.0 94.25
## 0.2cwt 98.88889 21.84407 18 64 140 83.75 95.0 112.50
## 0.4cwt 114.22222 22.31738 18 81 161 97.75 115.0 124.75
## 0.6cwt 123.38889 22.99908 18 86 174 106.25 121.5 139.00
##
## $comparison
## NULL
##
## $groups
## yield groups
## 0.6cwt 123.38889 a
## 0.4cwt 114.22222 a
## 0.2cwt 98.88889 b
## 0.0cwt 79.38889 c
##
## attr(,"class")
## [1] "group"library(gplots)
(plot.stuff <- with(sp.oats, HSD.test(yield, nitroF, DFerror = 45, MSerror = 117)))## $statistics
## MSerror Df Mean CV MSD
## 117 45 103.9722 10.40341 9.618527
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey nitroF 4 3.772697 0.05
##
## $means
## yield std r Min Max Q25 Q50 Q75
## 0.0cwt 79.38889 19.39417 18 53 117 63.25 72.0 94.25
## 0.2cwt 98.88889 21.84407 18 64 140 83.75 95.0 112.50
## 0.4cwt 114.22222 22.31738 18 81 161 97.75 115.0 124.75
## 0.6cwt 123.38889 22.99908 18 86 174 106.25 121.5 139.00
##
## $comparison
## NULL
##
## $groups
## yield groups
## 0.6cwt 123.38889 a
## 0.4cwt 114.22222 a
## 0.2cwt 98.88889 b
## 0.0cwt 79.38889 c
##
## attr(,"class")
## [1] "group"names(plot.stuff)## [1] "statistics" "parameters" "means" "comparison" "groups"plot.stuff$means## yield std r Min Max Q25 Q50 Q75
## 0.0cwt 79.38889 19.39417 18 53 117 63.25 72.0 94.25
## 0.2cwt 98.88889 21.84407 18 64 140 83.75 95.0 112.50
## 0.4cwt 114.22222 22.31738 18 81 161 97.75 115.0 124.75
## 0.6cwt 123.38889 22.99908 18 86 174 106.25 121.5 139.00Me dan diferentes los errores estándar???
plot.stuff$groups## yield groups
## 0.6cwt 123.38889 a
## 0.4cwt 114.22222 a
## 0.2cwt 98.88889 b
## 0.0cwt 79.38889 cbarplot2(plot.stuff$means[, 1])
# Add some lables
barplot2(plot.stuff$means[, 1], names.arg = rownames(plot.stuff$means), xlab = "Nitrogen",
ylab = "Yield per Acre")
# Add some error bars
mu.i <- plot.stuff$means[, 1]
se.i <- qt(1 - 0.05/2, 45) * plot.stuff$means[, 2]
bp <- barplot2(mu.i, names.arg = rownames(plot.stuff$means), xlab = "Nitrogen",
ylab = "Yield per Acre", plot.ci = TRUE, ci.l = mu.i - se.i, ci.u = mu.i +
se.i)
# Group them together
text(bp, 0, plot.stuff$groups[0, 3], cex = 1, pos = 3)
# No me lee, me toca añadirle un númer pero, qué es ese número?Desde aquí está mal.
# Esto es correcto?
plot.stuff$groups## yield groups
## 0.6cwt 123.38889 a
## 0.4cwt 114.22222 a
## 0.2cwt 98.88889 b
## 0.0cwt 79.38889 corder(plot.stuff$groups[, 1])## [1] 4 3 2 1#Por qué sale null????
plot.stuff$groups[order(plot.stuff$groups[, 1]), 3]## NULL# This is better
bp <- barplot2(mu.i, names.arg = rownames(plot.stuff$means), xlab = "Nitrogen",
ylab = "Yield per Acre", plot.ci = TRUE, ci.l = mu.i - se.i, ci.u = mu.i +
se.i)
text(bp, 0, plot.stuff$groups[order(plot.stuff$groups[, 1]), 3], cex = 1, pos = 3)
# Sometimes it is easier to do it by hand
bp <- barplot2(mu.i, names.arg = rownames(plot.stuff$means), xlab = "Nitrogen",
ylab = "Yield per Acre", plot.ci = TRUE, ci.l = mu.i - se.i, ci.u = mu.i +
se.i)
text(bp, 0, c("A", "B", "C", "C"), cex = 1, pos = 3)