Parcelas Divididas
Angie Alejandra Juyo Buitrago
22/10/2020
Diseño en parcelas divididas
El rendimiento de la avena (Variable respuesta) fue determinado a través de un ensayo de campo con diseño en parcelas divididas utilizando tres variedades y cuatro niveles de Nitrogeno.
Hay 3 parcelas principales, cada una dividida en 4 subparcelas;las variedades se aplicaron a las parcelas principales y los valores de nitrogeno a las subparcelas.
sp.avena <- read.csv("/cloud/project/avena.csv")
sp.avena <- within(sp.avena,nitroF <- factor(nitro))
head(sp.avena)## observation replicate variety nitro yield nitroF
## 1 1 I Victory 0.0 111 0
## 2 2 I Victory 0.2 130 0.2
## 3 3 I Victory 0.4 157 0.4
## 4 4 I Victory 0.6 174 0.6
## 5 5 I Golden 0.0 117 0
## 6 6 I Golden 0.2 114 0.2library(lattice)
library(car)## Loading required package: carDatalibrary(agricolae)
with(sp.avena, xyplot(yield ~ nitroF | variety))
with(sp.avena, xyplot(yield ~ nitroF | variety, groups = replicate))
with(sp.avena, xyplot(yield ~ nitroF | variety, groups = replicate, aspect = "xy"))
with(sp.avena, xyplot(yield ~ nitroF | variety, groups = replicate, aspect = "xy",
type = "o"))
with(sp.avena, xyplot(yield ~ nitroF | variety, groups = replicate, aspect = "xy",
type = "a"))
Luego de visualizar el comportamiento de los factores , se realiza un análisis de varianza.
Primero , se realiza un “Anova malo” que no toma en cuenta los errores y por tanto los grados de libertad de los residuales serán mayores y estos pueden realizar suposiciones o detecciones erradas
res.bad <- lm(yield ~ variety * nitroF, data = sp.avena)
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 factor “variedad” sólo tiene (18-1) grados de libertad total para utilizarlo al examinar las diferencias en la variedad. Por lo tanto, necesitamos especificar correctamente el término de error para la variedad.
res.buena <- aov(yield ~ variety * nitroF + Error(replicate:variety), data = sp.avena)## Warning in aov(yield ~ variety * nitroF + Error(replicate:variety), data =
## sp.avena): Error() model is singularsummary(res.buena)##
## 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 ' ' 1Como el P-valor es mayor de 5% , asumimos que no hay interacción entre los factores “variedad” y “nitro F”.Además, se rechaza la hipotesis de que las medias del factor nitrogeno son iguales pues , como no hay interacción con la variedad se puede analizar como un independiente y su p valor es menor al 5%.
Para comprobar las suposiciones, no es necesario utilizar el término de error.
res.buena2 <- aov(yield ~ variety * nitroF + replicate:variety, data = sp.avena)
summary(res.buena2)## 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 ' ' 1Se comprueban los supuestos para el análisis de varianza
plot(res.buena2, 1)
plot(res.buena2, 2)
boxCox(res.buena2)
Como en los valores de nitrogeno se rechazó la hipotesis de que las medias eran iguales , con el test HSD de Tukey para comparar las medias de los niveles de nitrógeno
HSD = with(sp.avena, HSD.test(yield, nitroF, DFerror = 45, MSerror = 117));HSD## $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 79.38889 19.39417 18 53 117 63.25 72.0 94.25
## 0.2 98.88889 21.84407 18 64 140 83.75 95.0 112.50
## 0.4 114.22222 22.31738 18 81 161 97.75 115.0 124.75
## 0.6 123.38889 22.99908 18 86 174 106.25 121.5 139.00
##
## $comparison
## NULL
##
## $groups
## yield groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0 79.38889 c
##
## attr(,"class")
## [1] "group"En la salida , podemos ver que las medias con la misma letra no son significativamente diferentes.(Niveles de 0.6 y 0.4 de nitrógeno)
library(gplots)##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowessplot.stuff <- with(sp.avena, HSD.test(yield, nitroF, 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 nitroF 4 3.772697 0.05
##
## $means
## yield std r Min Max Q25 Q50 Q75
## 0 79.38889 19.39417 18 53 117 63.25 72.0 94.25
## 0.2 98.88889 21.84407 18 64 140 83.75 95.0 112.50
## 0.4 114.22222 22.31738 18 81 161 97.75 115.0 124.75
## 0.6 123.38889 22.99908 18 86 174 106.25 121.5 139.00
##
## $comparison
## NULL
##
## $groups
## yield groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0 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 79.38889 19.39417 18 53 117 63.25 72.0 94.25
## 0.2 98.88889 21.84407 18 64 140 83.75 95.0 112.50
## 0.4 114.22222 22.31738 18 81 161 97.75 115.0 124.75
## 0.6 123.38889 22.99908 18 86 174 106.25 121.5 139.00plot.stuff$groups## yield groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0 79.38889 cbarplot2(plot.stuff$means[, 1])
barplot2(plot.stuff$means[, 1], names.arg = rownames(plot.stuff$means), xlab = "Nitrogeno",
ylab = "Rendimiento por Acre")
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 = "Rendimeinto por Acre", plot.ci = TRUE, ci.l = mu.i - se.i, ci.u = mu.i +
se.i)
text(bp, 0, plot.stuff$groups[, 2], cex = 1, pos = 3)
plot.stuff$groups## yield groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0 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"El resultado que obtenemos finalmente es un gráfico de barras que nos representa el rendimiento por acre de la avena en función de la cantidad de nitrógeno, con el nombre de los grupos y las respectivas barras de error
bp <- barplot2(mu.i, names.arg = rownames(plot.stuff$means), xlab = "Nitrogeno",
ylab = "Rendimienti 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)
Diseño en subparcelas divididas
En este experimento se desea medir los efectos de tres factores en la cantidad de glicógeno en el hígado.
En el experimento, hay 6 ratas (parcelas enteras). A cada rata se le asignó aleatoriamente una de las tres dietas alimenticias (T1, T2 y T3). A cada rata se le quitó el hígado y se dividió en cuatro segmentos. Cada segmento se preparó utilizando uno de dos químicos diferentes (P1 y P2). Finalmente, a cada segmento de hígado se le midió el nivel de glucógeno utilizando dos técnicas analíticas diferentes (A y B).
La unidad experimental para la dieta es la rata. La unidad experimental para el producto químico de preparación del hígado es la tira de hígado. La unidad experimental para la técnica analítica es un pedazo de hígado. Todos son diferentes.
spp.rata <- read.csv("/cloud/project/split_split_rat.csv"); spp.rata## glycogen food rat prep method
## 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
## 11 136 T1 Templeton P2 A
## 12 139 T1 Templeton P2 B
## 13 156 T1 Templeton P1 A
## 14 146 T1 Templeton P1 B
## 15 134 T1 Templeton P2 A
## 16 133 T1 Templeton P2 B
## 17 157 T2 Scabbers P1 A
## 18 145 T2 Scabbers P1 B
## 19 154 T2 Scabbers P2 A
## 20 142 T2 Scabbers P2 B
## 21 147 T2 Scabbers P1 A
## 22 153 T2 Scabbers P1 B
## 23 155 T2 Scabbers P2 A
## 24 157 T2 Scabbers P2 B
## 25 151 T2 Splinter P1 A
## 26 155 T2 Splinter P1 B
## 27 147 T2 Splinter P2 A
## 28 147 T2 Splinter P2 B
## 29 162 T2 Splinter P1 A
## 30 152 T2 Splinter P1 B
## 31 145 T2 Splinter P2 A
## 32 144 T2 Splinter P2 B
## 33 130 T3 Nicodemus P1 A
## 34 121 T3 Nicodemus P1 B
## 35 134 T3 Nicodemus P2 A
## 36 134 T3 Nicodemus P2 B
## 37 131 T3 Nicodemus P1 A
## 38 132 T3 Nicodemus P1 B
## 39 128 T3 Nicodemus P2 A
## 40 127 T3 Nicodemus P2 B
## 41 134 T3 Rizzo P1 A
## 42 136 T3 Rizzo P1 B
## 43 135 T3 Rizzo P2 A
## 44 134 T3 Rizzo P2 B
## 45 130 T3 Rizzo P1 A
## 46 123 T3 Rizzo P1 B
## 47 136 T3 Rizzo P2 A
## 48 133 T3 Rizzo P2 Bp <- as.vector(spp.rata$prep)
m <- as.vector(spp.rata$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:P2 X_Y <- xyplot(spp.rata$glycogen~(factor(m):factor(p))| spp.rata$food, groups = spp.rata$rat, aspect = "xy");X_Y
with(spp.rata, X_Y)
sp.res <- aov(glycogen ~ food * prep + Error(rat/food), data = spp.rata)## Warning in aov(glycogen ~ food * prep + Error(rat/food), data = spp.rata):
## 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.rata)## Warning in aov(glycogen ~ food * prep * method + Error(rat/food:prep), data =
## spp.rata): 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