Tres sucedaneos com diferente concentracao de proteina (tratamentos) foram testados em bezerros da raca holandes:
O experimento consistiu de um delineamento em blocos (9) completos casualizados.
pkg <- c("agricolae", "MASS", "ggplot2", "Rmisc", "plyr")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)
## agricolae MASS ggplot2 Rmisc plyr
## TRUE TRUE TRUE TRUE TRUE
#setwd("/media/DATA/Dropbox/Varios-BRA/Vachi/vachi_R/exp2") #PC Juanchi
#source('/media/DATA/Dropbox/MyR/Sources/boxplot.with.outlier.label.r', encoding='UTF-8') # pc juanchi
setwd("~/Dropbox/MyR/Análise otros/vachi_R/exp2")#lab juanchi
source('~/Dropbox/MyR/Sources/boxplot.with.outlier.label.r', encoding='UTF-8') # lab
#setwd("C:/Users/evangelina/Dropbox/Varios-BRA/Vachi/vachi_R/exp2") #pc vachi
#source('C:/Users/evangelina/Dropbox/Varios-BRA/Vachi/vachi_R/exp2/boxplot.with.outlier.label.r', encoding='UTF-8')
#list.files()
dat = read.csv("d_diarreia e vida.csv", h=T, s=";", d=",")
## ani trat bloc dd dv
## 1 2 B 2 37 56
## 2 3 C 1 23 56
## 3 4 C 3 4 17
## 4 5 B 3 40 56
## 5 6 A 1 35 56
## 6 7 A 3 24 56
## 'data.frame': 33 obs. of 5 variables:
## $ ani : int 2 3 4 5 6 7 8 9 10 11 ...
## $ trat: Factor w/ 3 levels "A","B","C": 2 3 3 2 1 1 2 3 1 2 ...
## $ bloc: int 2 1 3 3 1 3 1 2 2 6 ...
## $ dd : int 37 23 4 40 35 24 39 48 41 39 ...
## $ dv : int 56 56 17 56 56 56 56 56 56 56 ...
## bloc
## trat 1 2 3 4 5 6 7 8 9
## A 1 1 1 1 2 1 1 2 1
## B 1 1 1 1 2 1 1 1 1
## C 1 1 1 1 3 1 1 2 1
## Loading required package: TeachingDemos
# Ajuste de modelo
mod = aov(dd~ trat + bloc, data=dat[-12,])
# Diagnostico
par(mfrow=c(2,2))
plot(mod, which=1:2); boxcox(mod)
car::influencePlot(mod); layout(1)
## StudRes Hat CookD
## 14 2.2221839 0.3936189 0.4953667
## 15 -2.4034049 0.1871486 0.3138464
## 17 0.8999226 0.5437751 0.2975802
anova(mod) # cuadro anova
## Analysis of Variance Table
##
## Response: dd
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 2 1828.9 914.46 7.0200 0.004628 **
## bloc 8 1121.7 140.22 1.0764 0.415926
## Residuals 21 2735.6 130.26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# comparaciones multiples
TukeyHSD(x=mod, 'trat', conf.level=0.95)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = dd ~ trat + bloc, data = dat[-12, ])
##
## $trat
## diff lwr upr p adj
## B-A 7.484848 -5.445494 20.4151912 0.3301994
## C-A -11.015152 -23.023674 0.9933713 0.0760201
## C-B -18.500000 -31.185582 -5.8144178 0.0038523
(tuk = HSD.test(mod, 'trat', group=TRUE))
## $statistics
## Mean CV MSerror HSD r.harmonic
## 35.15625 32.46469 130.265 12.54754 10.51327
##
## $parameters
## Df ntr StudentizedRange
## 21 3 3.564625
##
## $means
## dd std r Min Max
## A 37.18182 8.244006 11 24 53
## B 44.66667 6.204837 9 37 54
## C 26.16667 16.151743 12 3 50
##
## $comparison
## NULL
##
## $groups
## trt means M
## 1 B 44.66667 a
## 2 A 37.18182 ab
## 3 C 26.16667 b
(sumdat <- summarySEwithin(dat[-12,], measurevar="dd", withinvars="trat",
na.rm=FALSE, conf.interval=.95))
## trat N dd sd se ci
## 1 A 11 37.18182 10.096804 3.044301 6.783125
## 2 B 9 44.66667 7.599342 2.533114 5.841371
## 3 C 12 26.16667 19.781764 5.710503 12.568733
df = data.frame(
let = tuk$groups[with(tuk$groups, order(trt)), ],
me = tuk$means$dd,
se = tuk$means$std/sqrt(tuk$means$r))
df
## let.trt let.means let.M me se
## 2 A 37.18182 ab 37.18182 2.485661
## 1 B 44.66667 a 44.66667 2.068279
## 3 C 26.16667 b 26.16667 4.662606
(p1<-ggplot(df, aes(let.trt, me)) + theme_bw() +
geom_point(size=3) +
geom_errorbar(aes(x = let.trt, ymin = me-se,
ymax = me+se),colour="black", width=.1) +
geom_line() +
xlab("Sucedáneos") + ylab("Dias com diarreia") +
annotate("text", x = (1:3)-0.1, y = sumdat$dd,
label = df$let.M)
)
## geom_path: Each group consist of only one observation. Do you need to adjust the group aesthetic?