Parcelas divididas

Diseño de parcela dividida en R

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 de la semana pasada. 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 (en bushels / acre) de avena. Son de interés las diferencias en la variedad de avena y los niveles de nitrógeno.

sp.oats <- read.csv2("C:/Users/CK0001/Desktop/ING.AGRONOMICA.UNAL/semestre 6/D.E/Tarea parcelas divididas/oats.csv", sep=";")
sp.oats <- within(sp.oats, nitroF <- factor(nitro))
head(sp.oats)

##   observation replicate     variety nitro yield nitroF
## 1           1         I     Victory   0.0   111    0.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.rain  0.00   117   0.00
## 6           6         I Golden.rain  0.20   114   0.20

library(lattice)  # Can only list one package at a time
library(car)

## Warning: package 'car' was built under R version 4.0.3

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.0.3

library(agricolae)

## Warning: package 'agricolae' was built under R version 4.0.3

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 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.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 ' ' 1

El análisis en R sigue el mismo patrón que la semana pasada. El 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.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 singular

summary(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 ' ' 1

Para verificar las suposiciones, no es necesario utilizar el término de error. Puede agregar el término sin error, pero las pruebas F son incorrectas. Sin embargo, la verificación de suposiciones está bien.

res.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 ' ' 1

plot(res.good2, 1)

plot(res.good2, 2)

## Warning: not plotting observations with leverage one:
##   1, 2, 3, 4

plot(res.good2, 5)

## Warning: not plotting observations with leverage one:
##   1, 2, 3, 4

boxCox(res.good2)

Dado que la gráfica de Box-Cox es buena, no es necesario realizar ninguna transformación. Aún necesita identificar qué niveles de nitrógeno son diferentes.

with(sp.oats, HSD.test(yield, nitroF, DFerror = 45, MSerror = 117))

library(gplots)

## Warning: package 'gplots' was built under R version 4.0.3

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

plot.stuff <- with(sp.oats, HSD.test(yield, nitroF, DFerror = 45, MSerror = 117))

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 133

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      c

barplot2(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)

# 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      c

order(plot.stuff$groups[, 1])

## [1] 8 7 6 5 4 3 2 1

plot.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)

Diseño de parcela 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.

spp.rat <- read.csv("C:/Users/CK0001/Desktop/ING.AGRONOMICA.UNAL/semestre 6/D.E/Tarea parcelas divididas/split_split_rat.csv")
spp.rat

##    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      B

pre <- as.vector(spp.rat$prep)
met <- as.vector(spp.rat$method)
factor(met):factor(pre)

##  [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.rat$glycogen~(factor(met):factor(pre))| spp.rat$food, groups = spp.rat$rat, aspect = "xy");X_Y

Este es un diseño loco, pero sucede. El tratamiento, la preparación y el método tienen EM denominador diferente para su prueba f. ¿Y si lo hiciste 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 ' ' 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(glycogen ~ food * prep + Error(rat/food), data = spp.rat)

## Warning in aov(glycogen ~ food * prep + Error(rat/food), data = spp.rat):
## Error() model is singular

summary(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.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(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 singular

summary(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

Parcelas divididas

Marly Daniela Quiroga Buitrago

22/10/2020

Diseño de parcela dividida en R

Ejemplo de avena

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.

Para verificar las suposiciones, no es necesario utilizar el término de error. Puede agregar el término sin error, pero las pruebas F son incorrectas. Sin embargo, la verificación de suposiciones está bien.

Dado que la gráfica de Box-Cox es buena, no es necesario realizar ninguna transformación. Aún necesita identificar qué niveles de nitrógeno son diferentes.

Diseño de parcela 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.

Este es un diseño loco, pero sucede. El tratamiento, la preparación y el método tienen EM denominador diferente para su prueba f. ¿Y si lo hiciste mal?

Primero, supongamos que se ignora el factor del método. Eso hace que este sea un diseño de parcela dividida.

Ahora agregue el método como una segunda división, por lo que ahora es un diseño de parcela dividida.