Diseño de Experimentos
GINNA VALENTINA ROMERO RINCON
23/10/2020
Ejemplo de Avena
En el siguiente ejemplo se van a evaluar tres tipos de avena, estos se encuentran en diez y ocho (18) parcelas completas, donde se encuentran subdivididas en cuatro (4).
Se importan los datos desde Excel y se crea la columna nitroF
library(readxl)
oats <- read_excel("oats.xlsx",
col_types = c("numeric", "text", "text",
"text", "numeric"))
oats <- within(oats, nitroF <- factor(nitro))
#View(oats)
head(oats)## # A tibble: 6 x 6
## Observacion replicate variedad nitro yield nitroF
## <dbl> <chr> <chr> <chr> <dbl> <fct>
## 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.0 117 0.0
## 6 6 I Golden.rain 0.2 114 0.2Se instalan las librerias para poder desarrollar los diversos calculos
library(lattice)
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.3library(agricolae)## Warning: package 'agricolae' was built under R version 4.0.3with(oats, xyplot(yield ~ nitroF | variedad))
Se cran diversos graficos para poder observar la dispercion existente segun las variedades
with(oats, xyplot(yield ~ nitroF | variedad, groups = replicate))
with(oats, xyplot(yield ~ nitroF | variedad, groups = replicate, aspect = "xy"))
with(oats, xyplot(yield ~ nitroF | variedad, groups = replicate, aspect = "xy",
type = "o"))
with(oats, xyplot(yield ~ nitroF | variedad, groups = replicate, aspect = "xy",
type = "a"))
Si no se tiene en cuenta el diseño experimental,los resultados incorrectos serian los siguientes:
res.bad <- lm(yield ~ variedad * nitroF, data = oats)
anova(res.bad)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## variedad 2 1786.4 893.2 1.7949 0.1750
## nitroF 3 20020.5 6673.5 13.4108 8.367e-07 ***
## variedad: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 ' ' 1res.good <- aov(yield ~ variedad * nitroF + Error(replicate:variedad), data = oats)## Warning in aov(yield ~ variedad * nitroF + Error(replicate:variedad), data =
## oats): Error() model is singularsummary(res.good)##
## Error: replicate: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)
## nitroF 3 20020 6673 37.686 2.46e-12 ***
## variedad: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(yield ~ variedad * nitroF + replicate:variedad, data = oats)
summary(res.good2)## Df Sum Sq Mean Sq F value Pr(>F)
## variedad 2 1786 893 5.044 0.0106 *
## nitroF 3 20020 6673 37.686 2.46e-12 ***
## variedad:nitroF 6 322 54 0.303 0.9322
## variedad:replicate 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, 3)
boxCox(res.good2)
Box-Cox plot no necesita ninguna transformacion
with(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':
##
## lowessplot.stuff <- with(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 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.0 79.38889 cbarplot2(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)
# Is this correct?
plot.stuff$groups## yield groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0.0 79.38889 corder(plot.stuff$groups[, 1])## [1] 4 3 2 1plot.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)
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.
library(readr)## Warning: package 'readr' was built under R version 4.0.3spp.rat <- read_csv("split_split_rat (1).csv")##
## -- Column specification ------
## cols(
## glycogen = col_double(),
## food = col_character(),
## rat = col_character(),
## prep = col_character(),
## method = col_character()
## )spp.rat## # A tibble: 48 x 5
## glycogen food rat prep method
## <dbl> <chr> <chr> <chr> <chr>
## 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
## # ... with 38 more rows# View(spp.rat)pr <- with(spp.rat, as.vector(prep))
meth <- with(spp.rat, as.vector(method))
with(spp.rat, xyplot(glycogen~(factor(meth):factor(pr))| food, groups = rat, aspect = "xy"))
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 ' ' 1sp.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