Parcelas divididas y subdivididas
Laura Sofia Hernandez Sanchez
20/10/2020
Diseño de parcela dividida (DPD)
El diseño recibe el nombre de parcelas divididas debido a que generalmente se asocia uno de los factores a unidades experimentales de mayor tamaño (parcela principal) y dentro de cada parcela principal se identifican “subparcelas” sobre las cuales se asigna al azar el segundo factor (UNLP, sf).
Ejemplo de cultivo de avena
Se realiza un experimento de parcelas divididas con tres variedades de avena en 6 parcelas principales, a las cuales se les mide su rendimiento apartir de cuatro dosis diferentes de Nitrogeno. Cabe açadir que para tener una mejor compresión, se realizara de una manera incorrecta para hacer un contraste con la forma adecuada.
Modelo que describe al diseño
\[Y_{ijk}=\mu+ \alpha_i+ \gamma_k+ \eta_{ ik} + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ijk}\]
\(Y_{ijk}\) = Rendimiento
\(\alpha_i\) = Efecto fijo de la variedad
\(\gamma_k\) = Efecto fijo del bloqueo
\(\eta_{ik}\) = Error en parcela principal
\(\beta_j\) = Efecto fijo Nitrogeno \(\alpha\beta_{ij}\) = Efecto interacción \(\epsilon_{ijk}\) = Error en subparcela
sp.oats_ <- read.csv("C:/Users/Sofia Hernandez/Documents/2020-2/Diseño de experimentos/oats.csv")
sp.oats<- within(sp.oats_, NF <- factor(N))
library(DT)
sp.oats1<-datatable(sp.oats, class='cell-border stripe', filter='top', options = list(pageLength=10,autoWidth=T), rownames=FALSE, colnames = c("X", "Repeticion (B)", "Variedad (V)", "Nitrogeno (N)", "Rendimiento (Y)", "NitroF"));sp.oats1table(sp.oats$V, sp.oats$NF, sp.oats$B) ## , , = I
##
##
## 0.0cwt 0.2cwt 0.4cwt 0.6cwt
## Golden.rain 1 1 1 1
## Marvellous 1 1 1 1
## Victory 1 1 1 1
##
## , , = II
##
##
## 0.0cwt 0.2cwt 0.4cwt 0.6cwt
## Golden.rain 1 1 1 1
## Marvellous 1 1 1 1
## Victory 1 1 1 1
##
## , , = III
##
##
## 0.0cwt 0.2cwt 0.4cwt 0.6cwt
## Golden.rain 1 1 1 1
## Marvellous 1 1 1 1
## Victory 1 1 1 1
##
## , , = IV
##
##
## 0.0cwt 0.2cwt 0.4cwt 0.6cwt
## Golden.rain 1 1 1 1
## Marvellous 1 1 1 1
## Victory 1 1 1 1
##
## , , = V
##
##
## 0.0cwt 0.2cwt 0.4cwt 0.6cwt
## Golden.rain 1 1 1 1
## Marvellous 1 1 1 1
## Victory 1 1 1 1
##
## , , = VI
##
##
## 0.0cwt 0.2cwt 0.4cwt 0.6cwt
## Golden.rain 1 1 1 1
## Marvellous 1 1 1 1
## Victory 1 1 1 1table(sp.oats$NF, sp.oats$B, sp.oats$V)## , , = Golden.rain
##
##
## I II III IV V VI
## 0.0cwt 1 1 1 1 1 1
## 0.2cwt 1 1 1 1 1 1
## 0.4cwt 1 1 1 1 1 1
## 0.6cwt 1 1 1 1 1 1
##
## , , = Marvellous
##
##
## I II III IV V VI
## 0.0cwt 1 1 1 1 1 1
## 0.2cwt 1 1 1 1 1 1
## 0.4cwt 1 1 1 1 1 1
## 0.6cwt 1 1 1 1 1 1
##
## , , = Victory
##
##
## I II III IV V VI
## 0.0cwt 1 1 1 1 1 1
## 0.2cwt 1 1 1 1 1 1
## 0.4cwt 1 1 1 1 1 1
## 0.6cwt 1 1 1 1 1 1table(sp.oats$B, sp.oats$V, sp.oats$NF)## , , = 0.0cwt
##
##
## Golden.rain Marvellous Victory
## I 1 1 1
## II 1 1 1
## III 1 1 1
## IV 1 1 1
## V 1 1 1
## VI 1 1 1
##
## , , = 0.2cwt
##
##
## Golden.rain Marvellous Victory
## I 1 1 1
## II 1 1 1
## III 1 1 1
## IV 1 1 1
## V 1 1 1
## VI 1 1 1
##
## , , = 0.4cwt
##
##
## Golden.rain Marvellous Victory
## I 1 1 1
## II 1 1 1
## III 1 1 1
## IV 1 1 1
## V 1 1 1
## VI 1 1 1
##
## , , = 0.6cwt
##
##
## Golden.rain Marvellous Victory
## I 1 1 1
## II 1 1 1
## III 1 1 1
## IV 1 1 1
## V 1 1 1
## VI 1 1 1# Balanceado, con 12 tratamientos
sp.oats$trat= interaction(sp.oats$V, sp.oats$NF)
med_trat = tapply(sp.oats$Y, sp.oats$trat, mean); med_trat## Golden.rain.0.0cwt Marvellous.0.0cwt Victory.0.0cwt Golden.rain.0.2cwt
## 80.00000 86.66667 71.50000 98.50000
## Marvellous.0.2cwt Victory.0.2cwt Golden.rain.0.4cwt Marvellous.0.4cwt
## 108.50000 89.66667 114.66667 117.16667
## Victory.0.4cwt Golden.rain.0.6cwt Marvellous.0.6cwt Victory.0.6cwt
## 110.83333 124.83333 126.83333 118.50000des_trat = tapply(sp.oats$Y, sp.oats$trat, sd);des_trat## Golden.rain.0.0cwt Marvellous.0.0cwt Victory.0.0cwt Golden.rain.0.2cwt
## 21.004761 16.573071 20.598544 13.472194
## Marvellous.0.2cwt Victory.0.2cwt Golden.rain.0.4cwt Marvellous.0.4cwt
## 26.853305 22.509257 29.944393 9.786044
## Victory.0.4cwt Golden.rain.0.6cwt Marvellous.0.6cwt Victory.0.6cwt
## 26.010895 20.875025 20.292035 30.091527cv_trat= 100*des_trat/med_trat;cv_trat ### Problemas de deviacion >20 %## Golden.rain.0.0cwt Marvellous.0.0cwt Victory.0.0cwt Golden.rain.0.2cwt
## 26.255952 19.122774 28.809152 13.677354
## Marvellous.0.2cwt Victory.0.2cwt Golden.rain.0.4cwt Marvellous.0.4cwt
## 24.749590 25.103261 26.114296 8.352243
## Victory.0.4cwt Golden.rain.0.6cwt Marvellous.0.6cwt Victory.0.6cwt
## 23.468477 16.722316 15.998976 25.393694CF_V=100*tapply(sp.oats$Y, sp.oats$V, sd)/ tapply(sp.oats$Y, sp.oats$V, mean)
CF_NF=100*tapply(sp.oats$Y, sp.oats$NF, sd)/ tapply(sp.oats$Y, sp.oats$NF, mean)
CF_TF=100*sd(sp.oats$Y)/ mean(sp.oats$Y)Acontinueción se presenta un diagrama de arbol para ejemplificar la distribucion de los factores
library(collapsibleTree)
collapsibleTree(sp.oats, c("V", "NF"))Se exponen diferentes graficos con sus pasos para evidenciar la disposición de los factores
library(lattice)
library(car)
library(agricolae)
with(sp.oats, xyplot(Y~NF|V))
with(sp.oats, xyplot(Y~NF|V, groups = B)) #replicate
# El nitrogeno en las diferentes replicaciones por color que esta anidada en la variedadwith(sp.oats, xyplot(Y~NF|V, groups = B, aspect = "xy"))
with(sp.oats, xyplot(Y~NF|V, groups = B, aspect = "xy", type="o"))
with(sp.oats, xyplot(Y~NF|V, groups = B, aspect = "xy", type="a"))
res.mal<-lm(Y~V*NF, data=sp.oats) # No tiene en cuenta el error
anova(res.mal)## 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
## NF 3 20020.5 6673.5 13.4108 8.367e-07 ***
## V:NF 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.buena<-aov(Y~V*NF+Error(B:V), data=sp.oats) # Mas error
summary(res.buena)##
## 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)
## NF 3 20020 6673 37.686 2.46e-12 ***
## V:NF 6 322 54 0.303 0.932
## Residuals 45 7969 177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se analiza la tabla comenzando de abajo hacia arriba; se observa "***" en el NF lo que indica que existen diferencias estadísticas entre los niveles (grupos) de nitrogeno, para observar estas diferencias se debe realizar la prueba de comparación multiple denominada TuckeyHSD o en este caso, HSD.test correspondiente a la libreria agricolae
0 ‘***’ = Tienen diferencia entre las dosis NF; No tienen interacción
0.001 ‘**’ = Interacción, no deben ser interpretados por separado
0.01 ‘*’ = Hay diferencia
0.05 ‘.’ = Cercano al 5%,
library(lme4)
res.buena.pp<-lmer(Y~B+V*NF+(1|B:V), data=sp.oats)
res.buena.pp## Linear mixed model fit by REML ['lmerMod']
## Formula: Y ~ B + V * NF + (1 | B:V)
## Data: sp.oats
## REML criterion at convergence: 485.1565
## Random effects:
## Groups Name Std.Dev.
## B:V (Intercept) 10.30
## Residual 13.31
## Number of obs: 72, groups: B:V, 18
## Fixed Effects:
## (Intercept) BII BIII
## 111.3611 -28.0833 -39.4167
## BIV BV BVI
## -37.1667 -44.4167 -39.0833
## VMarvellous VVictory NF0.2cwt
## 6.6667 -8.5000 18.5000
## NF0.4cwt NF0.6cwt VMarvellous:NF0.2cwt
## 34.6667 44.8333 3.3333
## VVictory:NF0.2cwt VMarvellous:NF0.4cwt VVictory:NF0.4cwt
## -0.3333 -4.1667 4.6667
## VMarvellous:NF0.6cwt VVictory:NF0.6cwt
## -4.6667 2.1667Verificación de las suposiciones, sin usar el error
res.buena2<-aov(Y~V*NF+B:V, data=sp.oats) # Sin error, pero pruebas F incorrectas
summary(res.buena2)## Df Sum Sq Mean Sq F value Pr(>F)
## V 2 1786 893 5.044 0.0106 *
## NF 3 20020 6673 37.686 2.46e-12 ***
## V:NF 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.buena2, 1) #Selecciona solo el primer grafico
plot(res.buena2, 2)
boxCox(res.buena2) # Estandarización 
Se necesita identificar qué niveles de nitrógeno son diferentes como se concluyo anteriormente al usar el anova, por lo que se usara HSD.test de la libreria agricolae
library(agricolae)
HSD<-with(sp.oats,(HSD.test(Y, NF, DFerror = 45, MSerror = 117))); HSD # RESIDUALS (error within); Df: error, MSerror= media## $statistics
## MSerror Df Mean CV MSD
## 117 45 103.9722 10.40341 9.618527
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey NF 4 3.772697 0.05
##
## $means
## Y 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
## Y 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"Se observo la clasificación de cuatro niveles de nitrogeno en tres grupos, “a” “b” y “c”, al realizar la prueba HSD.test se percibe que los niveles correspondientes a 0.4 y 0.6 son cercanos entre si, ya que se clasifican en un mismo grupo “a”, esto indica que para tomar desiciones agronomicas para un optimo rendimiento y menores costos se debe aplicar la dosificación de 0.4 de Nitrogeno.
library(gplots)
plot.stuff<-with(sp.oats, HSD.test(Y, NF, 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 NF 4 3.772697 0.05
##
## $means
## Y 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
## Y 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"medias_N=plot.stuff$means
names(medias_N)## [1] "Y" "std" "r" "Min" "Max" "Q25" "Q50" "Q75"gp<-plot.stuff$groups # Cuatro dosis de Nitrogeno library(gplots)
barplot2(medias_N [, 1], names.arg = rownames(medias_N), xlab= "Nitrogeno", ylab= "Rendimiento por Acre")
## Tres variedades de avena y cuantro dosis diferentes de Nitrogeno en 18 parcelas divididas en cuatro
mu.i <-medias_N[,1]
se.i <- qt(1-0.05/2, 45) * medias_N[,2]
bp <- barplot2(mu.i, names.arg = rownames(medias_N), 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)
http://www.personal.psu.edu/mar36/stat_461/split_plot/split_plot.html
bp2 <- barplot2(mu.i, names.arg = rownames(medias_N), xlab = "Nitrogeno", ylab = "Rendimiento por Acre", plot.ci = TRUE, ci.l = mu.i - se.i, ci.u = mu.i + se.i)
text(bp2, 0, gp[order(gp[, 1]), 3], cex = 1, pos = 3)
Modelo de parcelas divididas
Se desea saber los factores sobre la cantidad de glucogeno en el higado de 6 ratas, se les asigno una de tres dietas (T1, T2, T3) a cada una. Posteriormente se extrajo el higado y se dividio en cuatro partes, este se preparo usando dos productos quimicos diferentes P1 Y P2, finalmente se midio el nivel de glucogeno usando dos tecnicas A y B.
library(collapsibleTree)
library(readxl)
sp_rat <- read_excel("C:/Users/Sofia Hernandez/Documents/2020-2/Computacion estadistica/split_split_rat.xlsx")
datatable(sp_rat, class='cell-border stripe', filter='top', options = list(pageLength=10,autoWidth=T), rownames=FALSE, colnames = c("Glucogeno", "Dieta", "Rat", "Prep", "Metodo"))Se facilita la observaión de los datos por medio de un diagrama de arbol
collapsibleTree(sp_rat, c("food", "rat", "prep", "method"))pp <- as.vector(sp_rat$prep)
mt <- as.vector(sp_rat$method)
factor(mt):factor(pp)## [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:P2with(sp_rat, xyplot(glycogen ~ factor(factor(mt):factor(pp)) | food, groups = rat, aspect = "xy"))
# Mal
ft.malo<-anova(lm(glycogen ~ food*prep*method, data=sp_rat)) # No se tiene en cuenta la interaccion del error en las "parcelas"
ft.malo## 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=sp_rat)## Warning in aov(glycogen ~ food * prep + Error(rat/food), data = sp_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 si, un diseño de parcela dividida correcto
sp.res.b<- aov(glycogen~food*prep*method+Error(rat/food), data=sp_rat) # Se tiene en cuenta la interaccion del error## Warning in aov(glycogen ~ food * prep * method + Error(rat/food), data =
## sp_rat): Error() model is singularsummary(sp.res.b)##
## 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.033 0.317
## method 1 54.2 54.19 1.599 0.215
## food:prep 2 133.3 66.65 1.966 0.156
## food:method 2 5.4 2.69 0.079 0.924
## prep:method 1 11.0 11.02 0.325 0.572
## food:prep:method 2 5.3 2.65 0.078 0.925
## Residuals 33 1118.4 33.89Se realiza la prueba de varianza (aov) con la interacción de las variables, en donde se observa que no hay diferencia en “food” (dieta). Se debe tener en cuenta el “.” que designa un p.valor cercano al 5%.
https://aulavirtual.agro.unlp.edu.ar/pluginfile.php/4850/mod_resource/content/0/TP_-_Diseno_de_Parcela_Dividida.pdf http://pages.stat.wisc.edu/~larget/