tarea
Nombre: Erika Milena Aparicio Aponte
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, 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
Se plantaron 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 (72 observaciones) 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.csv("C:/Users/amor/Desktop/Rstudio/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
## 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.0 117 0
## 6 6 I Golden.rain 0.2 114 0.2library(lattice) # Can only list one package at a time## Warning: package 'lattice' was built under R version 4.0.3library(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.3library(agricolae)## Warning: package 'agricolae' was built under R version 4.0.3with(sp.oats, xyplot(yield ~ nitroF | variety))
En este grafico se observa el rendimiento segun la variedad de avena y el nivel de nitrogeno, en general se observa que al aumentar el nivel de nitrogeno aumenta el rendimiento en las tres variedades, sin embargo, tambien se observa a primera vista que aunque la variedad Victory presenta una variabilidad de los datos mas marcada que las demas variedades.
with(sp.oats, xyplot(yield ~ nitroF | variety, groups = replicate))
En el grafico anterior se le asigno un color diferente a los valores de obtenidos de rendimiento segun la replicacion (son 6 replicaciones para cada variedad), entonces, por ejemplo, los datos de la replicacion 1 tiene el color azul.
with(sp.oats, xyplot(yield ~ nitroF | variety, groups = replicate, aspect = "xy"))
with(sp.oats, xyplot(yield ~ nitroF | variety, groups = replicate, aspect = "xy",
type = "o"))
En el grafico anterior se observa la union de los valores de rendimiento segun la replicacion, se observa que la replicacion 1 de la variedad Victory se aleja bastante de los resultados de las otras 5 replicaciones, mientras que en la variedad Golden.rain y Marvellous no se observa esto.
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(grados de libertad) es mayor, lo que facilita la detección de diferencias que realmente no existen.
res.bad <- lm(yield ~ variety * nitroF, data = sp.oats) # lm:modelo lineal, su formula esta diciendo que la variable respuesta es el rendimiento y las covariables son la variable * el nitrogeno
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 ' ' 1No son 60 grados de libertad sino menos en residuales Segun lo anterior, el analisis tabla de la varianza, daria a entender que el nivel de nitrogeno es la principal fuente de variabilidad de los datos de rendimiento
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 singularsummary(res.good)##
## 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 ' ' 1La variabilidad de los datos es por los datos de nitrogeno
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 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 ' ' 1Disminuye los grados de libertad de los residuales, hay interaccion entre la variedad y las replicaciones. lA variabilidad es por los datos de nitrogeno.
Én este, la variabilidad esta en los niveles de nitrogeno
plot(res.good2, 1)
Esta grafica significaria que los datos atipicos son 65, 39, 30 la relacion entre residuales y valores predichos.
plot(res.good2, 2)
Significa como los residuales se ajustan a la distribucion normal
# plot(res.good2, 5)Grafico de los residuales vs el valor predicho dependiendo de la variedad, se observa que los atipicos esta en la variedad victory (65) y en la variedad marvellous (30 y 39)
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)) # (Honestly-significant-difference) de Tukey es un test de comparaciones múltiples. Permite comparar las medias de los t niveles del vactor nitrogeno despues de haber rechazado la Hipótesis nula de igualdad de medias mediante la técnica ANOVA.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(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 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])
# Add some lables
barplot2(plot.stuff$means[, 1], names.arg = rownames(plot.stuff$means), xlab = "Nitrogeno",
ylab = "Rendimiento por 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 = "Nitrogeno",
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)
Barras de error
# Is this correct?
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 1mirar esto
plot.stuff$groups[order(plot.stuff$groups[, 1]), 3]## NULL# This is better
bp <- barplot2(mu.i, names.arg = rownames(plot.stuff$means), xlab = "Nitrogeno",
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 = "Nitrogeno",
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/amor/Desktop/Rstudio/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 Bp <- as.vector(spp.rat$prep)
m <- as.vector(spp.rat$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:P2X_Y <- xyplot(spp.rat$glycogen~(factor(m):factor(p))| spp.rat$food, groups = spp.rat$rat, aspect = "xy");X_Y
with(spp.rat, X_Y)
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 ' ' 1La variabilidad se da en la dietas y no hay interaccciones
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.62sp.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