TRABAJO viernes
Juan Narvaez
20/10/2020
library(readxl)
oats <- read_excel("oats.xlsx")
head(oats)## # A tibble: 6 x 6
## Observacion Repliaca Variedad Nitrogeno Rendimiento NitroF
## <dbl> <chr> <chr> <chr> <dbl> <chr>
## 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.2library(lattice) # Can only list one package at a time
library(car)## Loading required package: carDatalibrary(agricolae)
with(oats, xyplot(oats$Rendimiento ~ oats$NitroF | oats$Variedad))
with(oats, xyplot(oats$Rendimiento ~ oats$NitroF | oats$Variedad, groups = oats$Repliaca))
with(oats,xyplot(oats$Rendimiento ~ oats$NitroF | oats$Variedad, groups = oats$Repliaca, aspect = "xy"))
with(oats,xyplot(oats$Rendimiento ~ oats$NitroF | oats$Variedad, groups = oats$Repliaca, aspect = "xy", type = "o"))
res.bad <- lm(oats$Rendimiento ~ oats$Variedad * oats$NitroF, data = oats)
anova(res.bad)## Analysis of Variance Table
##
## Response: oats$Rendimiento
## Df Sum Sq Mean Sq F value Pr(>F)
## oats$Variedad 2 1786.4 893.2 1.7949 0.1750
## oats$NitroF 3 20020.5 6673.5 13.4108 8.367e-07 ***
## oats$Variedad:oats$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(oats$Rendimiento ~ oats$Variedad * oats$NitroF + Error(oats$Repliaca:oats$Variedad), data = oats)## Warning in aov(oats$Rendimiento ~ oats$Variedad * oats$NitroF +
## Error(oats$Repliaca:oats$Variedad), : Error() model is singularsummary(res.good)##
## Error: oats$Repliaca:oats$Variedad
## Df Sum Sq Mean Sq F value Pr(>F)
## oats$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)
## oats$NitroF 3 20020 6673 37.686 2.46e-12 ***
## oats$Variedad:oats$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(oats$Rendimiento ~ oats$Variedad * oats$NitroF + oats$Repliaca:oats$Variedad, data = oats)
summary(res.good2)## Df Sum Sq Mean Sq F value Pr(>F)
## oats$Variedad 2 1786 893 5.044 0.0106 *
## oats$NitroF 3 20020 6673 37.686 2.46e-12 ***
## oats$Variedad:oats$NitroF 6 322 54 0.303 0.9322
## oats$Variedad:oats$Repliaca 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, 5)## hat values (leverages) are all = 0.375
## and there are no factor predictors; no plot no. 5boxCox(res.good2)
with(oats, HSD.test(oats$Rendimiento, oats$NitroF, DFerror = 45, MSerror = 117))library(gplots)##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowessplot.stuff <- with(oats, HSD.test(oats$Rendimiento, oats$NitroF, 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 oats$NitroF 4 3.772697 0.05
##
## $means
## oats$Rendimiento 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.00
##
## $comparison
## NULL
##
## $groups
## oats$Rendimiento groups
## 0.6 123.38889 a
## 0.4 114.22222 a
## 0.2 98.88889 b
## 0.0 79.38889 c
##
## attr(,"class")
## [1] "group"names(plot.stuff)## [1] "statistics" "parameters" "means" "comparison" "groups"plot.stuff$means## oats$Rendimiento 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## oats$Rendimiento 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])
barplot2(plot.stuff$means[, 1], names.arg = rownames(plot.stuff$means), xlab = "Nitrogen",
ylab = "Yield per Acre")
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## oats$Rendimiento 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]## NULLThis 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)
library(readxl)
split_split_rat <- read_excel("split_split_rat.xlsx")
method<- split_split_rat$method
method<- as.vector(method)
prep<- split_split_rat$prep
prep<- as.vector(prep)
with(split_split_rat, xyplot(glycogen ~ (factor(method):factor(prep)) | food, groups = rat, aspect = "xy"))
fact.bad <- lm(glycogen ~ food * prep * method, data = split_split_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 =split_split_rat )## Warning in aov(glycogen ~ food * prep + Error(rat/food), data =
## split_split_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 = split_split_rat )## Warning in aov(glycogen ~ food * prep * method + Error(rat/food:prep), data =
## split_split_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