Ajuste de modelos lineares e mistos no ambiente R |
|
| 09 e 10 de Outubro de 2014 - Piracicaba - SP |
Prof. Dr. Walmes M. Zeviani |
| Escola Superior de Agricultura “Luiz de Queiroz” - USP |
Lab. de Estatística e Geoinformação - LEG |
| Pós Graduação em Genética e Melhoramento de Plantas | Departamento de Estatística - UFPR |
##-----------------------------------------------------------------------------
## 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
## TRUE
source("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] splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.1 plyr_1.8.1 reshape_0.8.5 multcomp_1.3-6
## [5] TH.data_1.0-3 mvtnorm_1.0-0 doBy_4.5-10 MASS_7.3-34
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.1 lme4_1.1-7 Matrix_1.1-4
## [6] methods_3.1.1 minqa_1.2.3 nlme_3.1-117 nloptr_1.0.4 Rcpp_0.11.2
## [11] sandwich_2.3-1 stringr_0.6.2 tools_3.1.1 zoo_1.7-11
## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)
##-----------------------------------------------------------------------------
## 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 2 0.97 0.15 0.86
## ener 2 1 0.42 0.07 0.94
## sexo:ener 4 2 0.60 0.10 0.98
## Residuals 63 395 6.27
## 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 44.7 1.48 0.24
## ener 2 38 19.1 0.63 0.54
## sexo:ener 4 87 21.7 0.72 0.58
## Residuals 63 1905 30.2
## Não significativo é o resultado esperado se os tratamentos foram
## casualizados aos animais.
\[ \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_{ijk}, \sigma^2) \newline \mu_i &= \mu+\beta_i+\alpha_j+\tau_k \end{aligned} \]
em que \(\beta_i\) é o efeito dos níveis de sexo \(i\), \(\alpha_j\) o efeito dos níveis de energia (Kcal) na dieta \(j\). \(\mu\) é a média de \(Y\) na ausência do efeito dos termos mencioandos. \(\mu_{ijk}\) é 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.
##-----------------------------------------------------------------------------
## 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 200 99.9 3.62 0.033 *
## ener 2 56 27.9 1.01 0.370
## sexo:ener 4 35 8.8 0.32 0.863
## Residuals 63 1741 27.6
## ---
## 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 596 596 32.11 4.2e-07 ***
## id 1 47 47 2.55 0.116
## sexo 2 157 78 4.23 0.019 *
## ener 2 66 33 1.79 0.175
## sexo:ener 4 34 9 0.46 0.762
## Residuals 61 1131 19
## ---
## 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 5 0.28 0.599
## pi 1 638 638 34.38 2e-07 ***
## sexo 2 157 78 4.23 0.019 *
## ener 2 66 33 1.79 0.175
## sexo:ener 4 34 9 0.46 0.762
## Residuals 61 1131 19
## ---
## 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 1741
## 2 61 1131 2 610 16.4 2e-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.15
r <- 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 0 0.03 0.8542
## pi 1 459 459 39.36 4.8e-08 ***
## sexo 2 178 89 7.62 0.0011 **
## ener 2 145 72 6.21 0.0036 **
## sexo:ener 4 67 17 1.43 0.2350
## Residuals 58 677 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow=c(2,2)); plot(m2); layout(1)
##-----------------------------------------------------------------------------
## 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 pi
## 1 127.0 1.211 58 104.82 8.067e-68 124.5 129.4 F alto 138 92.01
## 2 124.9 1.209 58 103.31 1.861e-67 122.5 127.3 MC alto 138 92.01
## 3 129.8 1.291 58 100.50 9.132e-67 127.2 132.3 MI alto 138 92.01
## 4 122.8 1.212 58 101.32 5.713e-67 120.3 125.2 F baixo 138 92.01
## 5 124.0 1.209 58 102.59 2.783e-67 121.6 126.4 MC baixo 138 92.01
## 6 124.6 1.294 58 96.26 1.097e-65 122.0 127.1 MI baixo 138 92.01
## 7 125.7 1.219 58 103.07 2.131e-67 123.2 128.1 F medio 138 92.01
## 8 123.9 1.251 58 99.01 2.160e-66 121.4 126.4 MC medio 138 92.01
## 9 130.0 1.292 58 100.68 8.233e-67 127.5 132.6 MI medio 138 92.01
LSmeans(m2, effect=c("sexo"))
## estimate se df t.stat p.value lwr upr sexo id pi
## 1 125.1 0.7072 58 177.0 5.740e-81 123.7 126.6 F 138 92.01
## 2 124.3 0.7048 58 176.3 7.068e-81 122.9 125.7 MC 138 92.01
## 3 128.1 0.7456 58 171.8 3.150e-80 126.6 129.6 MI 138 92.01
LSmeans(m2, effect=c("ener"))
## estimate se df t.stat p.value lwr upr ener id pi
## 1 127.2 0.7141 58 178.1 3.904e-81 125.8 128.6 alto 138 92.01
## 2 123.8 0.7138 58 173.4 1.863e-80 122.3 125.2 baixo 138 92.01
## 3 126.5 0.7150 58 177.0 5.680e-81 125.1 128.0 medio 138 92.01
X <- 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.1 123.7 126.6 b
## 2 MC 124.3 122.9 125.7 b
## 3 MI 128.1 126.6 129.6 a
X <- 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.2 125.8 128.6 a
## 2 Médio 126.5 125.1 128.0 a
## 3 Baixo 123.8 122.3 125.2 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.12
## [1] 137.9
X <- 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.1 123.0 127.3 ab
## 2 MC 124.3 122.2 126.5 b
## 3 MI 128.2 126.0 130.3 a
X <- 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.1 123.4 126.9 ab
## 2 MC 124.3 122.5 126.0 b
## 3 MI 127.7 126.0 129.5 a
X <- 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.1 123.7 126.5 b
## 2 MC 124.3 122.9 125.7 b
## 3 MI 128.1 126.6 129.6 a
ps0$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 3
l <- 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))
