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)
source("lattice_setup.R")
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()
## 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/volume.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t", nrow=27)
## Mostra a estrutura do objeto.
str(da)
## Os dados no formato amplo (não empilhado).
unstack(da, form=volu~gen)
## Os dados são de um experimento para verificar o ganho de peso em
## função de tipos de rações. Foram usandos 4 níveis do fator ração e 5
## repetições.
##-----------------------------------------------------------------------------
## Análise exploratória.
## xyplot(volu~gen, data=da, type=c("p","a"))
xyplot(volu~reorder(gen, volu), data=da, type=c("p","a"))
\[ \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_i, \sigma^2) \newline \mu_i &= \mu+\alpha_i \end{aligned} \]
em que \(\alpha_i\) é o efeito do genótipo \(i\) sobre a variável resposta \(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito dos genótipos, \(\mu_i\) é a média para as observações em cada nível de genótipo e \(\sigma^2\) é a variância das observações ao redor desse valor médio.
##-----------------------------------------------------------------------------
## Estimação.
m0 <- lm(volu~gen, data=da)
## Estimativas e erros padrões.
summary(m0)
## Restrição paramétrica de zerer o primeiro nível.
contrasts(da$gen)
## model.matrix(m0)
##-----------------------------------------------------------------------------
## Existe efeito de genótipos? Ou seja, todos os \alpha_i podem ser
## considerados iguais à zero?
## Quadro de análise de variância.
anova(m0)
##-----------------------------------------------------------------------------
## Análise de diagnóstico dos pressupostos via resíduos.
## As inferências acima só passam ser valiadas após confirmação de não
## haver afastamento dos pressupostos.
par(mfrow=c(2,2))
plot(m0); layout(1)
##-----------------------------------------------------------------------------
## Transformar os dados pode diminuir a fuga dos pressupostos?
MASS::boxcox(m0)
##-----------------------------------------------------------------------------
## Ajuste com a variável transformada.
m1 <- lm(log(volu)~gen, data=da)
par(mfrow=c(2,2))
plot(m1); layout(1)
anova(m1)
## Passa m1 para m0.
m0 <- m1
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## Médias ajustadas (LSmeans).
p0 <- LSmeans(m0, effect="gen", level=0.95)
## Criando a tabela com as estimativas.
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, da)
p0
## Comparações múltiplas, contrastes de Tukey.
## Método de correção de p-valor: single-step.
g0 <- summary(glht(m0, linfct=mcp(gen="Tukey")),
test=adjusted(type="fdr"))
g0
## Resumo compacto com letras.
l0 <- cld(g0)
l0
## Formatando média seguida de letras.
p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
p0$let
##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- "Volume das raízes (log)"
ylab <- "Genótipo de sorgo"
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
p0 <- transform(p0, gen=reorder(gen, estimate))
levels(p0$gen)
p0 <- p0[order(p0$gen),]
segplot(gen~lwr+upr, centers=estimate,
data=p0, draw=FALSE,
xlab=xlab, ylab=ylab, sub=sub, letras=p0$let,
panel=function(x, y, z, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers, ...)
panel.text(y=as.numeric(z), x=centers,
label=letras, pos=3)
})
##-----------------------------------------------------------------------------
## Representando na escala original.
sel <- c("estimate","lwr","upr")
p0[,sel] <- exp(p0[,sel])
xlab <- "Volume das raízes"
p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
segplot(gen~lwr+upr, centers=estimate,
data=p0, draw=FALSE,
xlab=xlab, ylab=ylab, sub=sub, letras=p0$let,
panel=function(x, y, z, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers, ...)
panel.text(y=as.numeric(z), x=centers,
label=letras, pos=3)
})