Diseño por parcelas divididas en R
Paula Natalia Paez Monroy
22/10/2020
El diseño tradicional de parcela dividida es desde el punto de vista del análisis estadístico, similar al diseño de medidas repetidas de dos factores. El diseño consta de bloques (o parcelas enteras) en los que se aplica un factor (el factor de parcela completa) al azar. Dentro de cada parcela / bloque completo, se divide en unidades más pequeñas y los niveles del segundo factor se aplican aleatoriamente a las partes más pequeñas de toda la parcela. La clave es que la unidad experimental es diferente para cada factor.
Ejemplo de avena
En este ejemplo, se plantan campos enteros con uno de los tres tipos de avena. Estas son las parcelas completas, de las cuales hay 18. Cada parcela completa se divide en cuatro parcelas divididas, a cada una de las cuales se le asignó aleatoriamente uno de los cuatro niveles de nitrógeno. Se determinó el rendimiento de avena. Son de interés las diferencias en la variedad de avena y los niveles de nitrógeno.
datos.avena <- read.csv("/cloud/project/oats (1).csv")
head(datos.avena)## X Repetición Variedad Nitrogeno Rendimiento
## 1 1 I Victory 0.0cwt 111
## 2 2 I Victory 0.2cwt 130
## 3 3 I Victory 0.4cwt 157
## 4 4 I Victory 0.6cwt 174
## 5 5 I Golden.rain 0.0cwt 117
## 6 6 I Golden.rain 0.2cwt 114datos.avena_1 = within(datos.avena, NitrogenoF <-factor(Nitrogeno))
head(datos.avena_1)## X Repetición Variedad Nitrogeno Rendimiento NitrogenoF
## 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.2cwtlibrary(lattice)
library(car)## Loading required package: carDatalibrary(agricolae)
with(datos.avena_1, xyplot(Rendimiento ~ NitrogenoF | Variedad))
with(datos.avena_1, xyplot(Rendimiento ~ NitrogenoF | Variedad, groups = Repetición))
with(datos.avena_1, xyplot(Rendimiento ~ NitrogenoF | Variedad, groups = Repetición, aspect = "xy"))
with(datos.avena_1, xyplot(Rendimiento ~ NitrogenoF | Variedad, groups = Repetición, aspect = "xy",type = "o"))
with(datos.avena_1, xyplot(Rendimiento ~ NitrogenoF | Variedad, groups = Repetición, aspect = "xy",type = "a"))
Si ignora el diseño experimental, obtendrá los siguientes resultados incorrectos. Tenga en cuenta que el error df es mayor, lo que facilita la detección de diferencias que realmente no existen.
res.malo <- lm(Rendimiento ~ Variedad * NitrogenoF, data = datos.avena_1)
anova(res.malo)## Analysis of Variance Table
##
## Response: Rendimiento
## Df Sum Sq Mean Sq F value Pr(>F)
## Variedad 2 1786.4 893.2 1.7949 0.1750
## NitrogenoF 3 20020.5 6673.5 13.4108 8.367e-07 ***
## Variedad:NitrogenoF 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 ' ' 1res.bueno <- aov(Rendimiento ~ Variedad* NitrogenoF + Error(Repetición:Variedad), data = datos.avena_1)## Warning in aov(Rendimiento ~ Variedad * NitrogenoF +
## Error(Repetición:Variedad), : Error() model is singularsummary(res.bueno)##
## Error: Repetición:Variedad
## Df Sum Sq Mean Sq F value Pr(>F)
## Variedad 2 1786 893.2 0.612 0.555
## Residuals 15 21889 1459.2
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## NitrogenoF 3 20020 6673 37.686 2.46e-12 ***
## Variedad:NitrogenoF 6 322 54 0.303 0.932
## Residuals 45 7969 177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El factor de “variedad” solo tiene (18-1) gl total para usar cuando se prueban las diferencias en la variedad. Por lo tanto, necesitamos especificar correctamente el término de error para variedad. Cabe señalar que la forma del “Error (A: parcela completa)” puede cambiar dependiendo de la disposición de los datos. La clave es saber el gl correcto para que sepa tener los resultados correctos.
res.bueno2 = aov(Rendimiento ~ Variedad * NitrogenoF + Repetición:Variedad, data = datos.avena_1)
summary(res.bueno2)## Df Sum Sq Mean Sq F value Pr(>F)
## Variedad 2 1786 893 5.044 0.0106 *
## NitrogenoF 3 20020 6673 37.686 2.46e-12 ***
## Variedad:NitrogenoF 6 322 54 0.303 0.9322
## Variedad:Repetición 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.bueno2, 1)
plot(res.bueno2, 2)
boxCox(res.bueno2)
library(agricolae)
HSD_test = with(datos.avena_1, HSD.test(Rendimiento, NitrogenoF, DFerror = 45, MSerror = 117))
HSD_test## $statistics
## MSerror Df Mean CV MSD
## 117 45 103.9722 10.40341 9.618527
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey NitrogenoF 4 3.772697 0.05
##
## $means
## Rendimiento 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
## Rendimiento 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)##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowessplot.stuff <- with(datos.avena_1, HSD.test(Rendimiento, NitrogenoF, DFerror = 45, MSerror = 117))
plot.stuff## $statistics
## MSerror Df Mean CV MSD
## 117 45 103.9722 10.40341 9.618527
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey NitrogenoF 4 3.772697 0.05
##
## $means
## Rendimiento 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
## Rendimiento 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## Rendimiento 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.00plot.stuff$groups## Rendimiento 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])
# Añadir algunas etiquetas
barplot2(plot.stuff$means[, 1], names.arg = rownames(plot.stuff$means), xlab = "Nitrogeno", ylab = "Rendimiento por acre")
# Añadir las barras de error
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 = "Nitrogeno",
ylab = "Rendimiento por acre", plot.ci = TRUE, ci.l = mu.i - se.i, ci.u = mu.i +
se.i)
text(bp, 0, c("a", "a", "b", "c"), cex = 1, pos = 3)
# ¿Es correcto?
plot.stuff$groups## Rendimiento 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 1plot.stuff$groups[order(plot.stuff$groups[, 1]), 2]## [1] "c" "b" "a" "a"# Esto es mejor
bp <- barplot2(mu.i, names.arg = rownames(plot.stuff$means), xlab = "Nitrogeno",
ylab = "Rendimiento por 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]), 2], cex = 1, pos = 3)
# A veces es más fácil hacerlo a mano
bp <- barplot2(mu.i, names.arg = rownames(plot.stuff$means), xlab = "Nitrogeno",
ylab = "Rendimiento por 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)
Diseño de parcelas dividida-dividida-dividida
En este experimento, desea medir los efectos de tres factores sobre la cantidad de glucógeno en el hígado. En su experimento, hay 6 ratas (parcelas completas).
A cada rata, se le asignó al azar una de las tres dietas alimentarias (T1, T2 y T3).
De cada rata, se extrajo el hígado y se dividió en cuatro segmentos. Cada segmento se preparó utilizando uno de dos productos químicos diferentes (P1 y P2).
Finalmente, se midió el nivel de glucógeno de cada pedazo de hígado usando dos técnicas analíticas diferentes (A y B).
La unidad experimental de la dieta es la rata. La unidad experimental para la preparación química del hígado es una tira de hígado. La unidad experimental de la técnica analítica es un trozo de hígado. Todos son diferentes.
datos.ratas <- read.csv("datos.rat.csv")
head(datos.ratas)## Glucogeno Comida Rata Prep Metodo
## 1 127 T1 Remy P1 A
## 2 126 T1 Remy P1 B
## 3 127 T1 Remy P2 A
## 4 121 T1 Remy P2 B
## 5 124 T1 Remy P1 A
## 6 125 T1 Remy P1 Bp =as.vector(datos.ratas$Prep)
m = as.vector(datos.ratas$Metodo)
factor(m):factor(p)## [1] A:P1 B:P1 A:P2 B:P2 A:P1 B:P1 A:P2 B:P2 A:P1 B:P1 A:P2 B:P2 A:P1 B:P1 A:P2
## [16] B:P2 A:P1 B:P1 A:P2 B:P2 A:P1 B:P1 A:P2 B:P2 A:P1 B:P1 A:P2 B:P2 A:P1 B:P1
## [31] A:P2 B:P2 A:P1 B:P1 A:P2 B:P2 A:P1 B:P1 A:P2 B:P2 A:P1 B:P1 A:P2 B:P2 A:P1
## [46] B:P1 A:P2 B:P2
## Levels: A:P1 A:P2 B:P1 B:P2X_Y = xyplot(datos.ratas$Glucogeno~(factor(m):factor(p))| datos.ratas$Comida, groups = datos.ratas$Rata, aspect = "xy")
X_Y
with(datos.ratas, X_Y)
#¿Algo esta mal?
fact.malo <- lm(Glucogeno ~ Comida * Prep * Metodo, data = datos.ratas)
anova(fact.malo)## Analysis of Variance Table
##
## Response: Glucogeno
## Df Sum Sq Mean Sq F value Pr(>F)
## Comida 2 3530.0 1765.02 32.2583 9.401e-09 ***
## Prep 1 35.0 35.02 0.6401 0.4289
## Metodo 1 54.2 54.19 0.9904 0.3263
## Comida:Prep 2 133.3 66.65 1.2180 0.3077
## Comida:Metodo 2 5.4 2.69 0.0491 0.9521
## Prep:Metodo 1 11.0 11.02 0.2014 0.6563
## Comida:Prep:Metodo 2 5.3 2.65 0.0484 0.9529
## Residuals 36 1969.7 54.72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Primero, supongamos que se ignora el factor del método. Eso hace que este sea un diseño de parcela dividida.
sp.res <- aov(Glucogeno ~ Comida * Prep + Error(Rata/Comida), data = datos.ratas)## Warning in aov(Glucogeno ~ Comida * Prep + Error(Rata/Comida), data =
## datos.ratas): Error() model is singularsp.res##
## Call:
## aov(formula = Glucogeno ~ Comida * Prep + Error(Rata/Comida),
## data = datos.ratas)
##
## Grand Mean: 138.8542
##
## Stratum 1: Rata
##
## Terms:
## Comida Residuals
## Sum of Squares 3530.042 851.312
## Deg. of Freedom 2 3
##
## Residual standard error: 16.8455
## 2 out of 4 effects not estimable
## Estimated effects may be unbalanced
##
## Stratum 2: Within
##
## Terms:
## Prep Comida:Prep Residuals
## Sum of Squares 35.0208 133.2917 1194.3125
## Deg. of Freedom 1 2 39
##
## Residual standard error: 5.533841
## Estimated effects may be unbalancedsummary(sp.res)##
## Error: Rata
## Df Sum Sq Mean Sq F value Pr(>F)
## Comida 2 3530 1765.0 6.22 0.0856 .
## Residuals 3 851 283.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Prep 1 35.0 35.02 1.144 0.291
## Comida:Prep 2 133.3 66.65 2.176 0.127
## Residuals 39 1194.3 30.62# Ahora agregue el método como una segunda división, por lo que ahora es un diseño de parcela dividida.
sp.res <- aov(Glucogeno ~ Comida * Prep * Metodo + Error(Rata/Comida:Prep), data = datos.ratas)## Warning in aov(Glucogeno ~ Comida * Prep * Metodo + Error(Rata/Comida:Prep), :
## Error() model is singularsp.res##
## Call:
## aov(formula = Glucogeno ~ Comida * Prep * Metodo + Error(Rata/Comida:Prep),
## data = datos.ratas)
##
## Grand Mean: 138.8542
##
## Stratum 1: Rata
##
## Terms:
## Comida Residuals
## Sum of Squares 3530.042 851.312
## Deg. of Freedom 2 3
##
## Residual standard error: 16.8455
## 6 out of 8 effects not estimable
## Estimated effects may be unbalanced
##
## Stratum 2: Rata:Comida:Prep
##
## Terms:
## Prep Comida:Prep Residuals
## Sum of Squares 35.0208 133.2917 390.0625
## Deg. of Freedom 1 2 3
##
## Residual standard error: 11.40267
## 3 out of 6 effects not estimable
## Estimated effects may be unbalanced
##
## Stratum 3: Within
##
## Terms:
## Metodo Comida:Metodo Prep:Metodo Comida:Prep:Metodo Residuals
## Sum of Squares 54.1875 5.3750 11.0208 5.2917 728.3750
## Deg. of Freedom 1 2 1 2 30
##
## Residual standard error: 4.927389
## Estimated effects may be unbalancedsummary(sp.res)##
## Error: Rata
## Df Sum Sq Mean Sq F value Pr(>F)
## Comida 2 3530 1765.0 6.22 0.0856 .
## Residuals 3 851 283.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Rata:Comida:Prep
## Df Sum Sq Mean Sq F value Pr(>F)
## Prep 1 35.0 35.02 0.269 0.640
## Comida:Prep 2 133.3 66.65 0.513 0.643
## Residuals 3 390.1 130.02
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Metodo 1 54.2 54.19 2.232 0.146
## Comida:Metodo 2 5.4 2.69 0.111 0.896
## Prep:Metodo 1 11.0 11.02 0.454 0.506
## Comida:Prep:Metodo 2 5.3 2.65 0.109 0.897
## Residuals 30 728.4 24.28