Parcela_Dividida
Sergio Rocha
20/10/2020
Diseño de Parcela dividida en R
El diseño de parcela dividida tradicional es, desde un punto de vista estadístico, similar al diseño de dos factores con medidas repetidas(?). El diseño consiste en bloques (o parcelas completas) en los cuales un factor (El factor de parcela completo(?)) es aplicado aleatoriamente. Dentro cada parcela/bloque completo se divide en unidades más pequeñas (subunidades) y los niveles del segundo factor son aplicadas aleatoriamente a cada pieza pequeña de la parcela completa. La clave es que la unidad experimental es diferente para cada factor.
Ejemplo en Avena
En este ejemplo campos completos son cultivados con una de tres variedades de avena. Éstas son las parcelas completas, de las cuales hay 18. Cada parcela completa está dividida en 4 parcelas divididas, a cada una de las cuales les fue asignado aleatoriamente uno de cuatro niveles de nitrogeno. El rendimiento (en fanegadas/acre(?)) de avena fueron determinados. Son de interés las diferencias en la variedad de avena y los niveles de nitrógeno.
Modelo del diseño
\[y_{ijm}=\mu+\tau_i+\beta_j+(\tau\beta)_{ij}+ \epsilon_{ijm}\]
sp.oats <- read.csv("oats.csv")
N1 <- rep(seq(0.0, 0.6, 0.2), 18)
sp.oats <- cbind(sp.oats, N1)
sp.oats <- within(sp.oats, nitroF <- factor(N1))
library(DT)
datatable(sp.oats, class = 'cell-border stripe',filter = 'top', options = list(
pageLength = 6, autoWidth = TRUE))library(lattice) # Solo se puede listar un paquete a la vez
library(car)## Warning: package 'car' was built under R version 4.0.3## Loading required package: carDatalibrary(agricolae)
library(ggplot2)
library(collapsibleTree)
collapsibleTree(sp.oats, hierarchy = c('V', 'nitroF'), fill = rgb(0.5,0.2,0.5))#with(sp.oats, xyplot(Y ~ nitroF | V))
#with(sp.oats, xyplot(Y ~ nitroF | V, groups = B, pch = 22))
#with(sp.oats, xyplot(Y ~ nitroF | V, groups = B,aspect = 'xy', pch = 24))
#with(sp.oats, xyplot(Y ~ nitroF | V, groups = B,aspect = 'xy', type = "o"))
#with(sp.oats, xyplot(Y ~ nitroF | V, groups = B,aspect = 'xy', type = "a"))
fig1 <- ggplot(sp.oats, type="scatter", aes(x=N1, y=Y)) + geom_point(shape=1)
fig1 <- fig1 + facet_wrap(~V, ncol = 3)
fig1 <- fig1 + theme(plot.background = element_rect(fill = rgb(0.4,0.7,0.4, 0.8)), panel.background = element_rect(fill = rgb(1.0,0.9,0.9, 0.8)))
fig2 <- ggplot(sp.oats, type="scatter", aes(x=N1, y=Y, color=B)) + geom_point(shape=22)
fig2 <- fig2 + facet_wrap(~V, ncol = 3)
fig2 <- fig2 + theme(plot.background = element_rect(fill = rgb(0.4,0.7,0.4, 0.8)), panel.background = element_rect(fill = rgb(1.0,0.9,0.9, 0.8)))
fig3 <- ggplot(sp.oats, type="scatter", aes(x=N1, y=Y, color=B)) + geom_point(shape=16)
fig3 <- fig3 + facet_wrap(~V, ncol = 3)
fig3 <- fig3 + geom_line()
fig3 <- fig3 + theme(plot.background = element_rect(fill = rgb(0.4,0.7,0.4, 0.8)), panel.background = element_rect(fill = rgb(1.0,0.9,0.9, 0.8)))
fig4 <- ggplot(sp.oats, type="scatter", aes(x=N1, y=Y, color=B)) + geom_line()
fig4 <- fig4 + facet_wrap(~V, ncol = 3)
fig4 <- fig4 + theme(plot.background = element_rect(fill = rgb(0.4,0.7,0.4, 0.8)), panel.background = element_rect(fill = rgb(1.0,0.9,0.9, 0.8)))
fig1
fig2
fig3
fig4
Si ignoramos el diseño experimental, obtenemos los siguientes resutlados incorrectos. Note, el error en los grados de libertad es mayor, haciendo facil el detectar diferencias que no están realmente allí.
res.bad <- lm(Y ~ V * nitroF, data = sp.oats)
anova(res.bad)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## V 2 1786.4 893.2 1.7949 0.1750
## nitroF 3 20020.5 6673.5 13.4108 8.367e-07 ***
## V: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 analisis en R sigue el mismos análisis que la semana pasada, el factor ‘Variedad’ solo tiene (18-1) de grados de libertad totales para usar cuando probamos las diferencias en la variedad. Por lo tanto, necesitamos especificar el término del error de la variedad correctamente. Debe notarce que la forma de “Error(A:wholeplot)” puede cambiar dependiendo de la dispocición de los datos.
res.good <- aov(Y ~ V * nitroF + Error(B:V), data = sp.oats)## Warning in aov(Y ~ V * nitroF + Error(B:V), data = sp.oats): Error() model is
## singularsummary(res.good)##
## Error: B:V
## Df Sum Sq Mean Sq F value Pr(>F)
## V 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 ***
## V:nitroF 6 322 54 0.303 0.932
## Residuals 45 7969 177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res.good2 <- aov(Y ~ V * nitroF + B:V, data = sp.oats)
summary(res.good2)## Df Sum Sq Mean Sq F value Pr(>F)
## V 2 1786 893 5.044 0.0106 *
## nitroF 3 20020 6673 37.686 2.46e-12 ***
## V:nitroF 6 322 54 0.303 0.9322
## V:B 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)
## plot(res.good2, 5)
boxCox(res.good2)
Mientras el gráfico the box-cox sea bueno, no hay neceidad de hacer ninguna transformación. Aun asi necesitamos identificar cuales niveles de nitrógeno son diferentes.
## Error
library(agricolae)
hsd <- with(sp.oats, HSD.test(Y, 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
## Y 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
## Y groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0 79.38889 c
##
## attr(,"class")
## [1] "group"library(gplots)## Warning: package 'gplots' was built under R version 4.0.3##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowessplot.stuff <- with(sp.oats, HSD.test(Y, 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
## Y 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
## Y 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## Y 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## Y groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0 79.38889 cbarplot2(plot.stuff$means[, 1])
# Se añaden etiquetas
barplot2(plot.stuff$means[, 1], names.arg = rownames(plot.stuff$means), xlab = "Nitrógeno",
ylab = "Rendimiento por Acre")
# Se añaden 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 = "Nitrógeno",
ylab = "Rendimiento por Acre", plot.ci = TRUE, ci.l = mu.i - se.i, ci.u = mu.i +
se.i)
# Agruparlos
## Error
text(bp, 0, plot.stuff$groups[,2], cex = 1, pos = 3)
# ¿Es esto correcto?
plot.stuff$groups## Y 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]), 3]## NULL# esto es mejor
bp <- barplot2(mu.i, names.arg = rownames(plot.stuff$means), xlab = "Nitrógeno",
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]), 3], cex = 1, pos = 3)
## ErrorDiseño de parcela divida-divida-dividia(sub-subdividida?)
En este experimento se desea medir los efectos de 3 factores en la cantidad de glucógeno en el hígado. En el experimento hay 6 ratas (parcelas completas)
Para cada rata, una de tres dietas de alimento fue asignada aleatoriamente (T1, T2, T3).
De cada rata, el hígado fue extraido y dividido en 4 segmentos. Cada segmento fue preparado usando uno de dos químicos diferentes (P1, P2).
Finalmente, a cada pieza de hígado se le midio el nivel de glucógeno utilizando dos técnicas analíticas diferentes (A, B).
La unidad experimental para la dieta es la rata. La unidad experimental para la preparación química del hígado es tira de hígado. la unidad experimental para la técnica analítica es pedazo de hígado. Las tres son diferentes.
Modelo del diseño
\[y_{ijkm}=\mu+\tau_i+\beta_{j(i)}+\alpha_k+(\tau\alpha)_{ik}+(\beta\alpha)_{j(i)k}+ \epsilon_{ijkm}\]
spp.rat <- read.csv("split_split_rat.csv")
library(DT)
datatable(spp.rat, class = 'cell-border stripe',filter = 'top', options = list(
pageLength = 8, autoWidth = TRUE))p <- with(spp.rat, as.vector(prep))
m <- with(spp.rat, as.vector(method))
with(spp.rat, xyplot(glycogen~(factor(m):factor(p))| food, groups = rat, aspect = "xy"))
este es un diseño chifaldo pero ocurre. El tratamiento, la preparacion y el metodo todos tienen diferente denominador MS para su f-test.
¿y si lo hacemos mal?
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 ' ' 1Primero pretendamos que el factor método es ignorado. Esto lo convierte en un diseño de parcela dividida.
sp.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.62Ahora introduzcamos el metodo como segunda division, ahora es un diseño de parcela divida-dividida (subparcela)
sp.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