##============================================================================= ## Ajuste de modelos lineares e mistos no ambiente R ## 08 e 09 de Novembro de 2014 - Esalq/USP ## Pós Graduação de Genética e Melhoramento de Plantas ## ## Prof. Walmes Zeviani ## LEG - DEST - UFPR ##============================================================================= ##----------------------------------------------------------------------------- ## Definções do knitr. Não rodar. opts_chunk$set( cache=FALSE, tidy=FALSE, fig.width=7, fig.height=6, fig.align="center", dpi=150, dev="png", dev.args=list(png=list(family="Ubuntu Light", type="cairo"))) opts_chunk$set(fig.path="fig/script06") options(width=90) ##----------------------------------------------------------------------------- ## 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/ramalho_feijao1.txt" ## Importa os dados a partir do arquivo na web. da <- read.table(file=url, header=TRUE, sep="\t", colClasses=c("factor","factor",NA,NA)) ## Para facilitar, abreviar os nomes. names(da) <- abbreviate(names(da)) ## Mostra a estrutura do objeto. str(da) da$faml <- with(da, { neword <- levels(faml)[order(as.integer(levels(faml)))] factor(faml, levels=neword) }) ## Tabela de frequência. xtabs(~rptc+faml, data=da) ## Não são esses blocos grandes demais para serem considerados ## homogêneos? ## Tem registros perdidos? sum(is.na(da$prdc)) ##----------------------------------------------------------------------------- ## Análise exploratória. xyplot(prdc~faml|rptc, data=da, type=c("p","a")) ##----------------------------------------------------------------------------- ## Estimação. ## Esse modelo assume que Var(ue) == 0. Pode ser uma suposição ## forte. Na prática esse modelo não é útil pois não representa ## adequadamente às fontes de variação presentes. Será considerado aqui ## por razões didáticas. m0 <- lm(prdc~rptc+faml, data=da) ## Declarando o modelo especificando termo de efeito aleatório que ## corresponde à um estrato de variação. da$ue <- with(da, interaction(rptc, faml)) m0 <- aov(prdc~rptc+faml+Error(ue), data=da) ##----------------------------------------------------------------------------- ## Quadro de análise de variância estratificado. summary(m0) ##----------------------------------------------------------------------------- ## A classe 'aovlist' possuí poucos metodos disponíveis, o que a torna menos ## atraente. class(m0) methods(class="aovlist") ## Análise dos resíduos e comparações múltiplas são coisas complicadas ## de fazer com objetos de classe 'aovlist'. Nesse caso, pode-se reduzir ## os dados à totais ou médias e usar o número de observações como ## ponderador na análise. ## Por curiosidade, pode-se declarar o modelo com os termos todos de ## efeito fixo. O quadro de anova será igual exceto que não será ## estratificado em com valores de F obtidos usando o QMR com ## denominador (o que está errado). Isso permite acessar os resíduos ## porque o modelo será de classe 'lm' e é só para isso que se deve ## usá-lo. m1 <- lm(prdc~rptc+faml+ue, data=da) ## Diganóstico dos resíduos. par(mfrow=c(2,2)); plot(m1, which=1:3) MASS::boxcox(m1) abline(v=1/3, col=2); layout(1) ## Aponta clara relação esperança-variância. ##----------------------------------------------------------------------------- ## Agregando os dados para médias e pesos. Não esquecer que existem NA. db <- ddply(da, .(rptc, faml), summarise, n=sum(!is.na(prdc)), prdcm=mean(prdc, na.rm=TRUE)) str(db) ##----------------------------------------------------------------------------- ## Ajuste do modelo. m2 <- lm(prdcm~rptc+faml, db, weight=n) ## Diagnóstico. par(mfrow=c(2,2)); plot(m2, which=1:3) MASS::boxcox(m2) abline(v=1/3, col=2); layout(1) ## Cadê aquele afastamento todo? ## Lembre-se de uma propriedade envolvendo ordem de funções lineares e ## não lineares. sqrt(mean(x)) é diferente de mean(sqrt(x)). A ordem ## importa. ## Qual a correspondência entre os quadros de anova? summary(m0) anova(m2) ## O primeiro estrato é igual ao quadro de anova usando as médias e ## pesos. A vantagem é ter um objeto de classe 'lm' que tem muito mais ## métodos disponíveis. methods(class="lm") ##----------------------------------------------------------------------------- ## Comparações múltiplas. ## Haja visto que o F da análise de variância não rejeitou a hipótese do ## efeito de família ser nulo, não há necessidade de conduzir ## comparações múltiplas. Por razões didáticas, o código para isso será ## mantido. ## Médias ajustadas (LSmeans). p0 <- LSmeans(m2, effect="faml", 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(m2, linfct=mcp(faml="Tukey")), test=adjusted(type="fdr")) g0 ## Resumo compacto com letras. l0 <- cld(g0, decreasing=TRUE) l0 ## Formatando média seguida de letras. p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters) ## p0$let <- sprintf("%5.2f %s", p0$estimate, l0$mcletters$Letters) p0$let ##----------------------------------------------------------------------------- ## Representação gráfica dos resultados. xlab <- "Produção" ylab <- "Família de feijão" sub <- list( "Médias seguidas de mesma letra indicam diferença nula à 5%.", font=1, cex=0.8) p0 <- transform(p0, faml=reorder(faml, estimate)) levels(p0$faml) p0 <- p0[order(p0$faml),] xlim <- c(-1, 15) segplot(faml~lwr+upr, centers=estimate, data=p0, draw=FALSE, xlim=xlim, 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.integer(z), x=xlim[1], label=letras, pos=4, cex=0.8) }) ## A amplitude dos intervalos reflete o número diferente de observações ## para cada nível de família. ##----------------------------------------------------------------------------- ## Experimento com mudas de ?? onde está se avaliando os substratos. da <- read.table("http://www.leg.ufpr.br/~walmes/data/mudas.txt", header=TRUE, sep=";") da$bloc <- factor(da$bloc) str(da) ##----------------------------------------------------------------------------- ## Seleciona só uma variável para não tomar espaço com o print dos ## resultados. Em situação real, não precisa fazer isso. Omitir linhas ## com NA. Cuidado ao usar na.omit() com muitas colunas de resposta, ## pois no caso as linhas serão eliminadas pela união de linhas com NA, ## ou seja, se houver NA para apenas uma variável na linha, a linha toda ## é removida, o que descarta valores úteis das demais. da <- subset(da, select=c("trat","bloc","rep","psr")) da <- droplevels(na.omit(da)) str(da) ## Alto nível de desbalanceamento. xtabs(~trat+bloc, da) ## Este é um caso em que começou-se com 20 plantas por parcela e algumas ## morreram. Para-se fazer a análise deve-se assumir que as perdas ## foram ao acaso (MAR: missing at random) e não podem de forma alguma ## estarem associadas aos níveis dos fatores pois enviesaria conclusões. ## psr: peso seco de raizes xyplot(psr~trat|bloc, groups=rep, data=da, layout=c(NA,1), scales=list(x=list(rot=90, alternating=1))) ##----------------------------------------------------------------------------- ## Fazendo a análise com as médias e os pesos. db <- ddply(da, .(bloc, trat), summarise, psrm=mean(psr), n=length(psr)) str(db) m0 <- lm(psrm~bloc+trat, db, weight=n) par(mfrow=c(2,2)); plot(m0); layout(1) anova(m0) ##----------------------------------------------------------------------------- ## Comparações múltiplas (embora não tenha necessidade, mas é usual ## reportar o valor das médias ajustadas). X <- LSmatrix(m0, effect="trat") rownames(X) <- levels(da$trat) ## A função apmc() engloba todos os comandos antes usados para chegar a ## tabela de comparações múltiplas. ## grid <- apmc(X=X, model=m0, focus="trat", test="bonferroni") grid <- apmc(X=X, model=m0, focus="trat", test="fdr") grid <- transform(grid, trat=reorder(trat, estimate)) grid <- grid[order(grid$trat),] ##----------------------------------------------------------------------------- ## Gráfico. xlab <- "Peso seco de raízes" ylab <- "Substrato" sub <- list( "Médias seguidas de mesma letra indicam diferença nula à 5%.", font=1, cex=0.8) xlim <- c(0.1, 1) segplot(trat~lwr+upr, centers=estimate, data=grid, draw=FALSE, xlim=xlim, xlab=xlab, ylab=ylab, sub=sub, letras=sprintf("%.2f %s", grid$estimate, grid$cld), panel=function(x, y, z, centers, letras, ...){ panel.segplot(x=x, y=y, z=z, centers=centers, ...) panel.text(y=as.integer(z), x=xlim[1], label=letras, pos=4, cex=0.8) }) ## Porque o F e as comparações múltiplas às vezes discordam?