Experimento em arranjo fatorial e com covariáveis
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "reshape",
"plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)## lattice latticeExtra doBy multcomp reshape plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUEsource("lattice_setup.R")
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()## R version 3.1.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.3 plyr_1.8.1 reshape_0.8.5 multcomp_1.3-7
## [5] TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-11 MASS_7.3-34
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] rmarkdown_0.3.3 knitr_1.7
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 grid_3.1.1 htmltools_0.2.6
## [6] Matrix_1.1-4 Rcpp_0.11.3 sandwich_2.3-0 stringr_0.6.2 tools_3.1.1
## [11] yaml_2.1.13 zoo_1.7-10## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)Importação dos dados
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/ancova.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t")
## Trunca os nomes para 4 digitos.
names(da) <- substr(names(da), 1, 4)
## Mostra a estrutura do objeto.
str(da)## 'data.frame': 72 obs. of 8 variables:
## $ ener: Factor w/ 3 levels "alto","baixo",..: 2 2 2 2 2 2 2 2 3 3 ...
## $ sexo: Factor w/ 3 levels "F","MC","MI": 2 2 2 2 2 2 2 2 2 2 ...
## $ rep : int 1 2 3 4 5 6 7 8 1 2 ...
## $ pi : num 93.3 94 91.6 89.3 91.5 93 88.8 93.7 93.9 88.5 ...
## $ id : int 134 137 133 133 144 147 138 142 134 137 ...
## $ pviv: num 127 125 126 119 122 ...
## $ rend: num 82.3 83.1 80.7 82.7 82.9 ...
## $ peso: num 127 125 126 119 122 ...## Dados de experimento com nutrição de suínos. Animais foram pesados
## antes do experimento e tinham idade conhecida. Essas variáveis
## contínuas foram usadas para explicar/corrigir parte da variação
## presente e melhor comparar os níveis de energia na ração fornecidos.
##-----------------------------------------------------------------------------
## Análise exploratória.
## Gráficos, distribuição das id e pi nos grupos formados pelos fatores
## categóricos.
## xyplot(pi~id|sexo, groups=energia, data=da, auto.key=TRUE)
## xyplot(pi~id|energia, groups=sexo, data=da, auto.key=TRUE)
##-----------------------------------------------------------------------------
## Pode-se fazer uma anova dessas covariáveis para verificar se elas são
## separadas ou confundidas com os fatores experimentais. Como as
## variáveis são medidas antes do inicio do experimento e assumindo que
## os animais são aleatóriamente designinados aos níveis de energia,
## espera-se resultado não significativo.
## Testa se pi é separada por sexo*energia.
anova(lm(pi~sexo*ener, data=da))## Analysis of Variance Table
##
## Response: pi
## Df Sum Sq Mean Sq F value Pr(>F)
## sexo 2 1.93 0.9654 0.1541 0.8575
## ener 2 0.84 0.4204 0.0671 0.9352
## sexo:ener 4 2.38 0.5958 0.0951 0.9837
## Residuals 63 394.75 6.2658## Testa se id é separada por sexo*energia.
anova(lm(id~sexo*ener, data=da))## Analysis of Variance Table
##
## Response: id
## Df Sum Sq Mean Sq F value Pr(>F)
## sexo 2 89.36 44.681 1.4778 0.2359
## ener 2 38.11 19.056 0.6303 0.5358
## sexo:ener 4 86.89 21.722 0.7185 0.5825
## Residuals 63 1904.75 30.234## Não significativo é o resultado esperado se os tratamentos foram
## casualizados aos animais.Especificação e estimação do modelo
\[ \begin{aligned} Y|\text{fontes de variação} &\,\sim \text{Normal}(\mu_{ij}, \sigma^2) \newline \mu_{ij} &= \mu+\alpha_i+\tau_j+\gamma_{ij} \end{aligned} \]
em que \(\alpha\) é o efeito dos níveis de sexo \(i\), \(\tau_j\) o efeito dos níveis de energia (Kcal) na dieta \(j\), \(gamma_{ij}\) representa a interação . \(\mu\) é a média de \(Y\) na ausência do efeito dos termos mencioandos. \(\mu_{ij}\) é a média para as observações na combinação \(ij\) e \(\sigma^2\) é a variância das observações ao redor desse valor médio. No modelo de análise de covariância, os valores de peso inicial e idade serão usados para explicar parte da variação ao final do experimento. Sendo assim, o modelo considerando o efeito dessas variáveis é representado por
\[ \begin{aligned} Y|\text{fontes de variação} &\,\sim \text{Normal}(\mu+\beta_1 x+\beta_2 z+\mu_{ij}, \sigma^2) \newline \mu_{ij} &= \alpha_i+\tau_j+\gamma_{ij} \end{aligned} \]
em que \(\beta_1\) representa o efeito da covariável peso (\(x\)) expresso por um coeficiente angular e \(\beta_2\) o mesmo para o efeito da covariável idade inicial (\(z\)).
##-----------------------------------------------------------------------------
## Especificação dos modelos.
## A análise dos resíduos será conduzida depois. Será verificado
## primeiro o efeito da ordem dos termos.
## Só com as fontes de variação de interesse.
m0 <- lm(peso~sexo*ener, data=da)
anova(m0)## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## sexo 2 199.81 99.903 3.6153 0.03263 *
## ener 2 55.84 27.921 1.0104 0.36990
## sexo:ener 4 35.38 8.846 0.3201 0.86348
## Residuals 63 1740.92 27.634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Análise de variância (em experimentos não ortogonais a ordem dos
## termos é importante!).
m1 <- lm(peso~pi+id+sexo*ener, data=da)
anova(m1)## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## pi 1 595.58 595.58 32.1128 4.209e-07 ***
## id 1 47.23 47.23 2.5464 0.11571
## sexo 2 156.97 78.48 4.2318 0.01901 *
## ener 2 66.46 33.23 1.7918 0.17531
## sexo:ener 4 34.38 8.59 0.4634 0.76228
## Residuals 61 1131.33 18.55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Ocorre redução na SQ pois parte da varição é explicada por pi e id.
## Troca ordem de entrada apenas por curiosidade.
m1 <- lm(peso~id+pi+sexo*ener, data=da)
anova(m1)## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## id 1 5.17 5.17 0.2789 0.59937
## pi 1 637.63 637.63 34.3803 1.98e-07 ***
## sexo 2 156.97 78.48 4.2318 0.01901 *
## ener 2 66.46 33.23 1.7918 0.17531
## sexo:ener 4 34.38 8.59 0.4634 0.76228
## Residuals 61 1131.33 18.55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Testa o poder de explicação dos "blocos" contínuos, termos extras.
anova(m0, m1)## Analysis of Variance Table
##
## Model 1: peso ~ sexo * ener
## Model 2: peso ~ id + pi + sexo * ener
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 63 1740.9
## 2 61 1131.3 2 609.59 16.434 1.953e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m1); layout(1)
## Medidas de influência para as observações.
im <- influence.measures(m1)
summary(im)## Potentially influential observations of
## lm(formula = peso ~ id + pi + sexo * ener, data = da) :
##
## dfb.1_ dfb.id dfb.pi dfb.sxMC dfb.sxMI dfb.enrb dfb.enrm dfb.sxMC:nrb dfb.sxMI:nrb
## 25 -0.60 -0.42 0.95 -0.04 0.00 -0.02 0.02 0.05 0.68
## 37 -0.20 0.96 -0.43 0.07 0.00 0.01 -0.04 -0.04 0.09
## 46 0.40 -0.57 -0.05 -0.04 -0.82 0.00 0.02 0.01 0.53
## 48 0.28 0.30 -0.52 0.03 0.81 0.01 -0.01 -0.03 -0.54
## dfb.sxMC:nrm dfb.sxMI:nrm dffit cov.r cook.d hat
## 25 -0.01 -0.06 1.76_* 0.15_* 0.23 0.18
## 37 0.12 -0.43 -1.48_* 0.44_* 0.18 0.23
## 46 -0.09 0.51 -1.30_* 0.30_* 0.14 0.16
## 48 0.02 -0.53 1.27 0.31_* 0.13 0.15r <- which(im$is.inf[,"dffit"])
m2 <- lm(peso~id+pi+sexo*ener, data=da[-r,])
anova(m2)## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## id 1 0.40 0.40 0.0341 0.854171
## pi 1 459.25 459.25 39.3607 4.8e-08 ***
## sexo 2 177.89 88.94 7.6230 0.001150 **
## ener 2 144.97 72.49 6.2125 0.003592 **
## sexo:ener 4 66.83 16.71 1.4319 0.234964
## Residuals 58 676.73 11.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow=c(2,2)); plot(m2); layout(1)
Comparações múltiplas
##-----------------------------------------------------------------------------
## Comparações múltiplas para o desdobramento.
LSmeans(m2, effect=c("sexo","ener"))## estimate se df t.stat p.value lwr upr sexo ener id
## 1 126.9699 1.211308 58 104.82046 8.066598e-68 124.5452 129.3946 F alto 138.0145
## 2 124.9146 1.209095 58 103.31250 1.861161e-67 122.4944 127.3349 MC alto 138.0145
## 3 129.7600 1.291111 58 100.50258 9.132251e-67 127.1755 132.3444 MI alto 138.0145
## 4 122.7648 1.211616 58 101.32318 5.713283e-67 120.3395 125.1901 F baixo 138.0145
## 5 124.0100 1.208737 58 102.59473 2.782637e-67 121.5905 126.4296 MC baixo 138.0145
## 6 124.5507 1.293884 58 96.26108 1.097349e-65 121.9607 127.1406 MI baixo 138.0145
## 7 125.6840 1.219403 58 103.07010 2.131268e-67 123.2431 128.1249 F medio 138.0145
## 8 123.9038 1.251383 58 99.01350 2.159913e-66 121.3989 126.4087 MC medio 138.0145
## 9 130.0355 1.291530 58 100.68331 8.233295e-67 127.4503 132.6208 MI medio 138.0145
## pi
## 1 92.0087
## 2 92.0087
## 3 92.0087
## 4 92.0087
## 5 92.0087
## 6 92.0087
## 7 92.0087
## 8 92.0087
## 9 92.0087LSmeans(m2, effect=c("sexo"))## estimate se df t.stat p.value lwr upr sexo id pi
## 1 125.1396 0.7071868 58 176.9540 5.740376e-81 123.7240 126.5552 F 138.0145 92.0087
## 2 124.2761 0.7048357 58 176.3193 7.067888e-81 122.8653 125.6870 MC 138.0145 92.0087
## 3 128.1154 0.7456101 58 171.8262 3.149698e-80 126.6229 129.6079 MI 138.0145 92.0087LSmeans(m2, effect=c("ener"))## estimate se df t.stat p.value lwr upr ener id pi
## 1 127.2148 0.7141424 58 178.1365 3.903753e-81 125.7853 128.6443 alto 138.0145 92.0087
## 2 123.7752 0.7138463 58 173.3919 1.863060e-80 122.3462 125.2041 baixo 138.0145 92.0087
## 3 126.5411 0.7149777 58 176.9861 5.680494e-81 125.1099 127.9723 medio 138.0145 92.0087X <- LSmatrix(m2, effect=c("sexo"))
rownames(X) <- levels(da$sexo)
ps <- apmc(X, model=m2, focus="sexo", test="fdr")
ps## sexo estimate lwr upr cld
## 1 F 125.1396 123.7240 126.5552 b
## 2 MC 124.2761 122.8653 125.6870 b
## 3 MI 128.1154 126.6229 129.6079 aX <- LSmatrix(m2, effect=c("ener"))
rownames(X) <- levels(da$ener)
pe <- apmc(X, model=m2, focus="ener", test="fdr")
pe$ener <- factor(pe$ener, levels=c("alto","medio","baixo"),
labels=c("Alto","Médio","Baixo"))
pe <- arrange(pe, ener)
pe## ener estimate lwr upr cld
## 1 Alto 127.2148 125.7853 128.6443 a
## 2 Médio 126.5411 125.1099 127.9723 a
## 3 Baixo 123.7752 122.3462 125.2041 b## names(ps)[1] <- names(pe)[1] <- "nivel"
## p0 <- cbind(rbind(ps, pe), fator=rep(c("sexo","ener"), each=3))
## dput(levels(p0$nivel))
## str(p0)##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- expression("Peso aos 28 dias"~(kg))
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
p1 <-
segplot(sexo~lwr+upr, centers=estimate,
data=ps, draw=FALSE,
ylab="Nível de sexo", xlab=xlab, sub=sub,
letras=sprintf("%0.2f %s", ps$estimate, ps$cld),
panel=function(x, y, z, subscripts, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers,
subscripts=subscripts, ...)
panel.text(y=as.numeric(z)[subscripts],
x=centers[subscripts],
label=letras[subscripts], pos=3)
})
p2 <-
segplot(ener~lwr+upr, centers=estimate,
data=pe, draw=FALSE,
ylab="Nível de energia da dieta", xlab=xlab, sub=sub,
letras=sprintf("%0.2f %s", pe$estimate, pe$cld),
panel=function(x, y, z, subscripts, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers,
subscripts=subscripts, ...)
panel.text(y=as.numeric(z)[subscripts],
x=centers[subscripts],
label=letras[subscripts], pos=3)
})
## Gráfico final.
plot(p1, split=c(1,1,2,1), more=TRUE)
plot(p2, split=c(2,1,2,1), more=FALSE)
##-----------------------------------------------------------------------------
## Como fica o resultado sem usar as covariáveis?
mean(da$pi); mean(da$id)## [1] 92.11667## [1] 137.8889X <- LSmatrix(m0, effect=c("sexo"), at=list(pi=92, id=138))
rownames(X) <- levels(da$sexo)
ps0 <- apmc(X, model=m0, focus="sexo", test="fdr")
ps0## sexo estimate lwr upr cld
## 1 F 125.1167 122.9724 127.2610 ab
## 2 MC 124.3250 122.1807 126.4693 b
## 3 MI 128.1875 126.0432 130.3318 aX <- LSmatrix(m1, effect=c("sexo"), at=list(pi=92, id=138))
rownames(X) <- levels(da$sexo)
ps1 <- apmc(X, model=m1, focus="sexo", test="fdr")
ps1## sexo estimate lwr upr cld
## 1 F 125.1306 123.3506 126.9106 ab
## 2 MC 124.2653 122.4914 126.0392 b
## 3 MI 127.7405 125.9733 129.5076 aX <- LSmatrix(m2, effect=c("sexo"), at=list(pi=92, id=138))
rownames(X) <- levels(da$sexo)
ps2 <- apmc(X, model=m2, focus="sexo", test="fdr")
ps2## sexo estimate lwr upr cld
## 1 F 125.1316 123.7156 126.5476 b
## 2 MC 124.2682 122.8576 125.6787 b
## 3 MI 128.1074 126.6150 129.5999 aps0$model <- "0"; ps1$model <- "1"; ps2$model <- "2"
ps <- rbind(ps0, ps1, ps2)
ps$model <- factor(ps$model)
str(ps)## 'data.frame': 9 obs. of 6 variables:
## $ sexo : Factor w/ 3 levels "F","MC","MI": 1 2 3 1 2 3 1 2 3
## $ estimate: num 125 124 128 125 124 ...
## $ lwr : num 123 122 126 123 122 ...
## $ upr : num 127 126 130 127 126 ...
## $ cld : chr "ab" "b" "a" "ab" ...
## $ model : Factor w/ 3 levels "0","1","2": 1 1 1 2 2 2 3 3 3l <- c("Fatorial", "Ancova", "Limpo")
key <- list(title="Modelo", cex.title=1.1,
corner=c(0.05,0.95), type="o", divide=1,
text=list(l), lines=list(pch=1:length(l)))
segplot(sexo~lwr+upr, centers=estimate,
data=ps, groups=model, draw=FALSE,
ylab="Nível de sexo", xlab=xlab, sub=sub, key=key,
panel=panel.segplotBy, f=0.075, pch=as.integer(ps$model))