glicose
setwd("~/Dropbox/MyR/Análise otros/vachi_R")pkg <- c("ggplot2", "knitr","reshape2", "dplyr", "nlme")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)## ggplot2 knitr reshape2 dplyr nlme
## TRUE TRUE TRUE TRUE TRUElibrary(data.table)
outlierReplace = function(dataframe, cols, rows, newValue = NA) {
if (any(rows)) {
set(dataframe, rows, cols, newValue)
}
}Carregar dataset
## ani trat bloco 1 2 3 4 5 6 7 8
## 1 21 2 2 NA 97.24 99.56 85.03 112.10 109.50 110.20 104.70
## 2 22 2 2 99.21 89.63 116.20 81.45 96.17 119.50 98.76 100.80
## 3 23 1 1 88.92 89.56 92.11 129.90 76.62 82.13 111.40 106.80
## 4 24 1 1 102.50 106.20 118.10 89.74 135.80 142.10 78.82 115.00
## 5 26 3 2 94.23 99.13 77.84 79.09 103.30 90.43 108.20 111.70
## 6 27 1 2 89.32 78.88 66.20 70.32 104.00 100.80 82.60 93.11Adaptar dataset
# Dados long expandidos
dat = melt(datlab, id.vars = c("ani", "trat", "bloco"))
names(dat)[4:5] = c("semana", "glico") ; head(dat) ## ani trat bloco semana glico
## 1 21 2 2 1 NA
## 2 22 2 2 1 99.21
## 3 23 1 1 1 88.92
## 4 24 1 1 1 102.50
## 5 26 3 2 1 94.23
## 6 27 1 2 1 89.32dat = dat[with(dat, order(trat, ani,semana)), ]
head(dat)## ani trat bloco semana glico
## 3 23 1 1 1 88.92
## 44 23 1 1 2 89.56
## 85 23 1 1 3 92.11
## 126 23 1 1 4 129.90
## 167 23 1 1 5 76.62
## 208 23 1 1 6 82.13for (i in 1:3) dat[,i] <- as.factor(dat[,i])
for (i in 4:5) dat[,i] <- as.numeric(dat[,i])
str(dat)## 'data.frame': 328 obs. of 5 variables:
## $ ani : Factor w/ 41 levels "21","22","23",..: 3 3 3 3 3 3 3 3 4 4 ...
## $ trat : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ bloco : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ semana: num 1 2 3 4 5 6 7 8 1 2 ...
## $ glico : num 88.9 89.6 92.1 129.9 76.6 ...datc <- dat[complete.cases(dat),] # tirar NA's Grafico exploratorio
explot = ggplot(datc, aes(semana, glico)) +
geom_point(shape=20, size=2) +
geom_smooth() +
facet_wrap(~trat, ncol=3, nrow=1) +
scale_y_continuous("Glicose") +
scale_x_continuous(breaks=seq(1, 8, by=1))
explot## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
# outliers?
qplot(data = datc, x = glico) +
scale_x_continuous(breaks=seq(0, 170, by=10))## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
datc = filter(datc, glico > 60 , glico < 140 )
qplot(data = datc, x = glico) +
scale_x_continuous(breaks=seq(0, 170, by=10))## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
explot1 = ggplot(datc, aes(semana, glico)) +
geom_point(shape=20, size=2) +
geom_smooth() +
facet_wrap(~trat, ncol=3, nrow=1) +
scale_y_continuous("Glicose") +
scale_x_continuous(breaks=seq(1, 8, by=1))
explot1## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
Análise da variância
Efeitos fixos
mod_glico = aov(glico ~ trat * semana + bloco, data=datc)
summary(mod_glico)## Df Sum Sq Mean Sq F value Pr(>F)
## trat 2 972 486 1.847 0.159514
## semana 1 3386 3386 12.865 0.000389 ***
## bloco 3 450 150 0.570 0.635033
## trat:semana 2 30 15 0.058 0.943984
## Residuals 307 80791 263
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Modelo misto
mod_glico <- lme(glico ~ trat * semana, random=~1|ani, data=datc)
summary(mod_glico)## Linear mixed-effects model fit by REML
## Data: datc
## AIC BIC logLik
## 2651.955 2681.848 -1317.978
##
## Random effects:
## Formula: ~1 | ani
## (Intercept) Residual
## StdDev: 3.968443 15.72267
##
## Fixed effects: glico ~ trat * semana
## Value Std.Error DF t-value p-value
## (Intercept) 90.52277 3.494466 272 25.904607 0.0000
## trat2 0.48337 5.136924 38 0.094097 0.9255
## trat3 3.35616 4.964298 38 0.676060 0.5031
## semana 1.26796 0.654852 272 1.936260 0.0539
## trat2:semana 0.34260 0.958772 272 0.357331 0.7211
## trat3:semana 0.16913 0.932004 272 0.181474 0.8561
## Correlation:
## (Intr) trat2 trat3 semana trt2:s
## trat2 -0.680
## trat3 -0.704 0.479
## semana -0.849 0.577 0.597
## trat2:semana 0.580 -0.854 -0.408 -0.683
## trat3:semana 0.596 -0.406 -0.851 -0.703 0.480
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.55856036 -0.72323285 -0.02076069 0.66852816 2.49532366
##
## Number of Observations: 316
## Number of Groups: 41