Parcela Dividida
Kevin Quiroga
21/10/2020
Ejemplo
Este ejemplo se hace con 3 tipos de avena
datatable(sp_oats,filter = 'top',class = 'cell-border stripe', options = list(
pageLength = 5, autoWidth = TRUE))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 se ignora el diseño experimental, se incurriran en resultados erroneos. El error en df es muy alto
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 2.5289 0.08882 .
## nitroF 7 30184.7 4312.1 12.2088 2.348e-09 ***
## variety:nitroF 6 235.9 39.3 0.1113 0.99476
## Residuals 56 19779.0 353.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res.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 1.042 0.37872
## nitroF 1 9883 9883 11.524 0.00436 **
## Residuals 14 12006 858
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## nitroF 6 20302 3384 18.283 2.66e-10 ***
## variety:nitroF 6 236 39 0.212 0.971
## Residuals 42 7773 185
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res.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 4.826 0.013 *
## nitroF 7 30185 4312 23.300 1.48e-12 ***
## variety:nitroF 6 236 39 0.212 0.971
## variety:replicate 14 12006 858 4.634 5.41e-05 ***
## Residuals 42 7773 185
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(res.good2, 1)
plot(res.good2, 2)## Warning: not plotting observations with leverage one:
## 1, 2, 3, 4
plot(res.good2,)## Warning: not plotting observations with leverage one:
## 1, 2, 3, 4
boxCox(res.good2)
HSD <- HSD.test(sp_oats$yield, sp_oats$nitroF, DFerror = 45, MSerror = 117)
with(sp_oats, HSD)## $statistics
## MSerror Df Mean CV
## 117 45 103.9722 10.40341
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey sp_oats$nitroF 8 4.493894 0.05
##
## $means
## sp_oats$yield std r Min Max Q25 Q50 Q75
## 0.0 111.00000 NA 1 111 111 111 111 111
## 0.00 77.52941 18.26238 17 53 117 63 70 89
## 0.2 130.00000 NA 1 130 130 130 130 130
## 0.20 97.05882 21.04599 17 64 140 82 91 108
## 0.4 157.00000 NA 1 157 157 157 157 157
## 0.40 111.70588 20.20138 17 81 161 97 112 121
## 0.6 174.00000 NA 1 174 174 174 174 174
## 0.60 120.41176 19.81180 17 86 156 104 121 133
##
## $comparison
## NULL
##
## $groups
## sp_oats$yield groups
## 0.6 174.00000 a
## 0.4 157.00000 a
## 0.2 130.00000 ab
## 0.60 120.41176 b
## 0.40 111.70588 b
## 0.0 111.00000 bc
## 0.20 97.05882 bc
## 0.00 77.52941 c
##
## attr(,"class")
## [1] "group"plot.stuff <- with(sp_oats, HSD.test(yield, nitroF, DFerror = 45, MSerror = 117))
plot.stuff## $statistics
## MSerror Df Mean CV
## 117 45 103.9722 10.40341
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey nitroF 8 4.493894 0.05
##
## $means
## yield std r Min Max Q25 Q50 Q75
## 0.0 111.00000 NA 1 111 111 111 111 111
## 0.00 77.52941 18.26238 17 53 117 63 70 89
## 0.2 130.00000 NA 1 130 130 130 130 130
## 0.20 97.05882 21.04599 17 64 140 82 91 108
## 0.4 157.00000 NA 1 157 157 157 157 157
## 0.40 111.70588 20.20138 17 81 161 97 112 121
## 0.6 174.00000 NA 1 174 174 174 174 174
## 0.60 120.41176 19.81180 17 86 156 104 121 133
##
## $comparison
## NULL
##
## $groups
## yield groups
## 0.6 174.00000 a
## 0.4 157.00000 a
## 0.2 130.00000 ab
## 0.60 120.41176 b
## 0.40 111.70588 b
## 0.0 111.00000 bc
## 0.20 97.05882 bc
## 0.00 77.52941 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.0 111.00000 NA 1 111 111 111 111 111
## 0.00 77.52941 18.26238 17 53 117 63 70 89
## 0.2 130.00000 NA 1 130 130 130 130 130
## 0.20 97.05882 21.04599 17 64 140 82 91 108
## 0.4 157.00000 NA 1 157 157 157 157 157
## 0.40 111.70588 20.20138 17 81 161 97 112 121
## 0.6 174.00000 NA 1 174 174 174 174 174
## 0.60 120.41176 19.81180 17 86 156 104 121 133plot.stuff$groups## yield groups
## 0.6 174.00000 a
## 0.4 157.00000 a
## 0.2 130.00000 ab
## 0.60 120.41176 b
## 0.40 111.70588 b
## 0.0 111.00000 bc
## 0.20 97.05882 bc
## 0.00 77.52941 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[,2], cex = 1, pos = 3)
#f <- plot.stuff$groups[,];f# Is this correct?
plot.stuff$groups## yield groups
## 0.6 174.00000 a
## 0.4 157.00000 a
## 0.2 130.00000 ab
## 0.60 120.41176 b
## 0.40 111.70588 b
## 0.0 111.00000 bc
## 0.20 97.05882 bc
## 0.00 77.52941 corder(plot.stuff$groups[, 1])## [1] 8 7 6 5 4 3 2 1plot.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)
Split-Split-Split Plot Design
spp.rat <- read_csv("D:/Kevin/Trabajos/Diseno de Experimentos/split_split_rat.csv")## Parsed with column specification:
## cols(
## glycogen = col_double(),
## food = col_character(),
## rat = col_character(),
## prep = col_character(),
## method = col_character()
## )spp.rat## # A tibble: 48 x 5
## glycogen food rat prep method
## <dbl> <chr> <chr> <chr> <chr>
## 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 B
## 7 132 T1 Remy P2 A
## 8 138 T1 Remy P2 B
## 9 146 T1 Templeton P1 A
## 10 144 T1 Templeton P1 B
## # ... with 38 more rowsrat <- spp.rat$rat
food <- spp.rat$food
glycogen <- spp.rat$glycogen
prep <- spp.rat$prep
method <- spp.rat$methodp <- as.vector(prep)
m <- as.vector(method);m## [1] "A" "B" "A" "B" "A" "B" "A" "B" "A" "B" "A" "B" "A" "B" "A" "B" "A" "B" "A"
## [20] "B" "A" "B" "A" "B" "A" "B" "A" "B" "A" "B" "A" "B" "A" "B" "A" "B" "A" "B"
## [39] "A" "B" "A" "B" "A" "B" "A" "B" "A" "B"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(glycogen~(factor(m):factor(p))| food, groups = rat, aspect = "xy");X_Y
with(spp.rat, X_Y)
fact.bad <- lm(glycogen ~ food * prep * method, data = spp.rat)
anova(fact.bad)## Analysis of Variance Table
##
## Response: glycogen
## Df Sum Sq Mean Sq F value Pr(>F)
## food 2 3530.0 1765.02 32.2583 9.401e-09 ***
## prep 1 35.0 35.02 0.6401 0.4289
## method 1 54.2 54.19 0.9904 0.3263
## food:prep 2 133.3 66.65 1.2180 0.3077
## food:method 2 5.4 2.69 0.0491 0.9521
## prep:method 1 11.0 11.02 0.2014 0.6563
## food:prep:method 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 ' ' 1sp.res <- aov(glycogen ~ food * prep + Error(rat/food), data = spp.rat)## Warning in aov(glycogen ~ food * prep + Error(rat/food), data = spp.rat):
## Error() model is singularsummary(sp.res)##
## Error: rat
## Df Sum Sq Mean Sq F value Pr(>F)
## food 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
## food:prep 2 133.3 66.65 2.176 0.127
## Residuals 39 1194.3 30.62sp.res <- aov(glycogen ~ food * prep * method + Error(rat/food:prep), data = spp.rat)## Warning in aov(glycogen ~ food * prep * method + Error(rat/food:prep), data =
## spp.rat): Error() model is singularsummary(sp.res)##
## Error: rat
## Df Sum Sq Mean Sq F value Pr(>F)
## food 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: rat:food:prep
## Df Sum Sq Mean Sq F value Pr(>F)
## prep 1 35.0 35.02 0.269 0.640
## food: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)
## method 1 54.2 54.19 2.232 0.146
## food:method 2 5.4 2.69 0.111 0.896
## prep:method 1 11.0 11.02 0.454 0.506
## food:prep:method 2 5.3 2.65 0.109 0.897
## Residuals 30 728.4 24.28