Tarea Diseño en parcelas
Catherine O
21/10/2020
Diseño en parcelas divididas
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 midio el rendimiento en en (bushels / acre) (es decir 14,515kg de avena en EUU/0.404686 ha) de avena. Son de interés las diferencias en la variedad de avena y los niveles de nitrógeno. Los niveles de nitrogeno expresados en cwt corresponden a centenas (50.8023kg).
- Resumiendo:
- 18 Parcelas principales dividia en 4, y estas cuatro se les asigna aletoriamente los 4 niveles de nitrogeno
- Rendimiento (bushels / acre)= 14,515kg/0.404686 ha
- Nitrogeno (ctw) centenas = 50.8023kg
library(readxl)
oats_tarea <- read_excel("C:/Users/Catis/Desktop/Disenno de experimentos/oats tarea.xlsx")
View(oats_tarea)
#Se importan los datosavena <- (oats_tarea)
avena <- within(avena, nitrogenoF <- factor(nitrogeno))
head(avena)## # A tibble: 6 x 6
## Observacion replica variedad nitrogeno rendimiento nitrogenoF
## <dbl> <chr> <chr> <dbl> <dbl> <fct>
## 1 1 I Victory 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.rain 0 117 0
## 6 6 I Golden.rain 0.2 114 0.2library(lattice) # Se va a graficar
library(car) #Calcula tablas de análisis de varianza de tipo II o tipo III para objetos modelo producidos por lm, glm,## Loading required package: carDatalibrary(agricolae)## Warning: package 'agricolae' was built under R version 4.0.3with(avena, xyplot(rendimiento ~ nitrogenoF | variedad))
with(avena, xyplot(rendimiento ~ nitrogenoF | variedad, groups = replica))
with(avena, xyplot(rendimiento ~ nitrogenoF | variedad, groups = replica, aspect = "xy"))
with(avena, xyplot(rendimiento ~ nitrogenoF | variedad, groups = replica, aspect = "xy",
type = "o"))
with(avena, xyplot(rendimiento ~ nitrogenoF | variedad, groups = replica, aspect = "xy",
type = "a"))
## Modelo
El factor de la parcela principal es la variedad. Una parcela dividida está dada por una subparcela de terreno.
El factor de la subparcela viene dado por las dosis de nitrógeno.
El factor Variedad se bloquea, por ende, la unidad experimental de la parcela principal vienedada por la combinación de bloque y variedad.
\[\require{cancel} y_{ijk}=\mu+\tau_i+\gamma_k+(\eta)_{ik}+(\beta)_j+(\alpha\beta)_{ij}+\epsilon_{ijk}\] - Donde:
- \(y_{ijk}\) Es el rendimiento de la avena
- \(\tau_{i}\) Es el efecto de la variedad
- \(\gamma_{k}\) Es el efecto del bloqueo
- \((\eta)_{ik}\) Es el error de la parcela principal
- \((\beta)_j\) Es el efecto del nitrogeno
- \((\alpha\beta)_{ij}\) Es el efecto de la interaccion variedad: nitrogeno
- \(\epsilon_{ijk}\) Es el error de las subparcelas
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.malos <- lm(rendimiento ~ variedad * nitrogenoF, data = avena)
anova(res.malos)## Analysis of Variance Table
##
## Response: rendimiento
## Df Sum Sq Mean Sq F value Pr(>F)
## variedad 2 1786.4 893.2 1.7949 0.1750
## nitrogenoF 3 20020.5 6673.5 13.4108 8.367e-07 ***
## variedad:nitrogenoF 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 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 df correcto para que sepa tener los resultados correctos.
- la clave son los grados de libertad correctos (df)
res.buenos <- aov(rendimiento ~ variedad * nitrogenoF + Error(replica:variedad), data = avena)## Warning in aov(rendimiento ~ variedad * nitrogenoF + Error(replica:variedad), :
## Error() model is singularsummary(res.buenos)##
## Error: replica:variedad
## Df Sum Sq Mean Sq F value Pr(>F)
## variedad 2 1786 893.2 0.612 0.555
## Residuals 15 21889 1459.2
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## nitrogenoF 3 20020 6673 37.686 2.46e-12 ***
## variedad:nitrogenoF 6 322 54 0.303 0.932
## Residuals 45 7969 177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Para 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.buenos2 <- aov(rendimiento ~ variedad * nitrogenoF + replica:variedad, data = avena)
summary(res.buenos2)## Df Sum Sq Mean Sq F value Pr(>F)
## variedad 2 1786 893 5.044 0.0106 *
## nitrogenoF 3 20020 6673 37.686 2.46e-12 ***
## variedad:nitrogenoF 6 322 54 0.303 0.9322
## variedad:replica 15 21889 1459 8.240 1.61e-08 ***
## Residuals 45 7969 177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Comparando res.buenos y res.buenos2:
- PAra variedad, cambia el valor F, y el p valor, pero y en res.buenos2 la variedad si es significativa.
- Nitrogeno no cambia entre las dos, sigue siendo significativo
- en res.buenos2 aparece variedad:replica (corresponde al error), que en res.buenos se llama residuals en la parte de error: replica:variedad, y tiene ** lo que quiere decir que es significativa
Revision de supuestos
- Independencia
plot(res.buenos2$residuals) #el que yo supongo que es 
plot(res.buenos2, 1) #el de la pagina
- Normalidad
qqnorm(res.buenos2$residuals) # los que creo que son
qqline(res.buenos2$residuals)
plot(res.buenos2, 2) #el de la pagina
#plot(res.buenos2, 5)boxCox(res.buenos2)
Una transformación de Box Cox es una transformación de una variable dependiente no normal en una forma normal. La normalidad es un supuesto importante para muchas técnicas estadísticas; Si sus datos no son normales, la aplicación de Box-Cox significa que puede ejecutar una mayor cantidad de pruebas
Se utiliza para modificar la forma de distribución de un conjunto de datos para que se distribuyan de manera más normal, de modo que las pruebas y los límites de confianza que requieren normalidad se puedan utilizar de forma adecuada.
- 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
prueba<-with(avena, HSD.test(rendimiento, nitrogenoF, DFerror = 45, MSerror = 117))
prueba## $statistics
## MSerror Df Mean CV MSD
## 117 45 103.9722 10.40341 9.618527
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey nitrogenoF 4 3.772697 0.05
##
## $means
## rendimiento 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
## rendimiento groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0 79.38889 c
##
## attr(,"class")
## [1] "group"#no sale como en la pagina
#Comparaciones Múltiples: Tukey, Realiza múltiples comparaciones de tratamientos mediante Tukey. El nivel por defecto alfa es 0.05.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(avena, HSD.test(rendimiento, nitrogenoF, DFerror = 45, MSerror = 117)) #no quedo como en la pagina
plot.stuff## $statistics
## MSerror Df Mean CV MSD
## 117 45 103.9722 10.40341 9.618527
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey nitrogenoF 4 3.772697 0.05
##
## $means
## rendimiento 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
## rendimiento 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 #tampoco sale bien, la columna de std es diferente, y en la pagina tampoco aparecen los Q25, Q50 y Q75## rendimiento 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 #si sale bien, lo unico es que el rendimientoaqui no aparece aproximado## rendimiento groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0 79.38889 cbarplot2(plot.stuff$means[, 1])
# añadimos algunas tablas
barplot2(plot.stuff$means[, 1], names.arg = rownames(plot.stuff$means), xlab = "Nitrogeno",
ylab = "Rendimiento per 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 = "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[1, 3], cex = 1, pos = 3)
# ES esto correcto?
plot.stuff$groups## rendimiento 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#ESta mal, deberia de aparecer otra cosa# esto es mejor
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
Split split
En este experimento, desea medir los efectos de tres factores sobre la cantidad de glucógeno en el hígado.
- Hay 6 ratas (parcelas completas)
- 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("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 Bpre <- 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:P2X_Y <- xyplot(spp.rat$glycogen~(factor(met):factor(pre))| spp.rat$food, groups = spp.rat$rat, aspect = "xy");X_Y
with(spp.rat, X_Y) El tratamiento, la preparación y el método tienen EM denominador diferente para su prueba F.
Haciendolo 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 ' ' 1Supongamos 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 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 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 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