##----------------------------------------------------------------------------- ## Definções do knitr. Não rodar. source("knitr_setup.R") opts_chunk$set(fig.path="fig/anc1") ##----------------------------------------------------------------------------- ## 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/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) ## 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)) ## Testa se id é separada por sexo*energia. anova(lm(id~sexo*ener, data=da)) ## Não significativo é o resultado esperado se os tratamentos foram ## casualizados aos animais. ##----------------------------------------------------------------------------- ## 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) ## 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) ## 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) ## Testa o poder de explicação dos "blocos" contínuos, termos extras. anova(m0, m1) ##----------------------------------------------------------------------------- ## 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) r <- which(im$is.inf[,"dffit"]) m2 <- lm(peso~id+pi+sexo*ener, data=da[-r,]) anova(m2) par(mfrow=c(2,2)); plot(m2); layout(1) ##----------------------------------------------------------------------------- ## Comparações múltiplas para o desdobramento. LSmeans(m2, effect=c("sexo","ener")) LSmeans(m2, effect=c("sexo")) LSmeans(m2, effect=c("ener")) X <- LSmatrix(m2, effect=c("sexo")) rownames(X) <- levels(da$sexo) ps <- apmc(X, model=m2, focus="sexo", test="fdr") ps 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 ## 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) 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 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 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 ps0$model <- "0"; ps1$model <- "1"; ps2$model <- "2" ps <- rbind(ps0, ps1, ps2) ps$model <- factor(ps$model) str(ps) 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))