##----------------------------------------------------------------------------- ## Definções do knitr. Não rodar. opts_chunk$set( cache=FALSE, tidy=FALSE, fig.width=6, fig.height=4.5, fig.align="center", dpi=100, dev="png", dev.args=list(png=list(family="Ubuntu Light", type="cairo"))) options(width=90) ##============================================================================= ## Aplicação de modelos de regressão linear e ## não linear em ciências agrárias ## ## 09 à 11 de Dezembro de 2014 - Goiânia/GO ## Embrapa Arroz e Feijão ## ## Prof. Dr. Walmes M. Zeviani ## LEG/DEST/UFPR ##============================================================================= ##----------------------------------------------------------------------------- ## Definições da sessão. pkg <- c("lattice", "latticeExtra", "gridExtra", "car", "alr3", "plyr", "reshape", "doBy", "multcomp", "ellipse", "wzRfun") sapply(pkg, require, character.only=TRUE) trellis.device(color=FALSE) extendseq <- function(x, f=0.1, length.out=20){ er <- extendrange(x, f=f) seq(er[1], er[2], length.out=length.out) } ##----------------------------------------------------------------------------- ## Tamanho de tartarugas em função do sexo. url <- "http://westfall.ba.ttu.edu/isqs5349/Rdata/turtles.txt" tur <- read.table(url, header=TRUE, sep="\t") str(tur) xtabs(~gender, tur) ## xyplot(width+height~length, groups=gender, data=tur, type=c("p","smooth"), ## scales=list(y=list(relation="free"))) splom(tur[,1:3], groups=tur$gender, type=c("p","smooth"), auto.key=TRUE) ##----------------------------------------------------------------------------- ## Ajustar reta separada por gender. São só 2 níveis, então codificada ## como numérica (dummy) ou fator implica no mesmo modelo. ## Possíveis modelos: ## fórmula | n. par. | equação ## ~1 | 1 | b0 ## ~length | 2 | b0+b1*x ## ~gender+length | 3 | b0+a0*(sexo==1)+b1*x ## ~gender*length | 4 | b0+a0*(sexo==1)+b1*x+a1*(sexo==1) ## Modelo com interação (modelo maior). m0 <- lm(height~gender*length, data=tur) par(mfrow=c(2,2)); plot(m0); layout(1) summary(m0) ## Usando gender como fator de dois níveis. tur$Gender <- factor(tur$gender) str(tur) m1 <- lm(height~Gender*length, data=tur) m1 <- lm(height~Gender*length, data=tur, contrasts=list(Gender=contr.SAS)) summary(m1) ## contr.treatment ## contr.SAS ## contr.sum ## contr.helmert ## contr.poly ## Diagnóstico par(mfrow=c(2,2)); plot(m1); layout(1) ## Essa parametrização de nível de referência é útil para testar ## hipóteses. A parametrização de estimativas por categoria é mais ## interessante para interpretação. m2 <- update(m1, .~0+Gender/length) summary(m2) ## Veja que o R² aumentou pelo fato de o modelo declarado não ter o ## modelo nulo como modelo mínimo por considerar efeitos iguais a zero. ## Mesmo nessa parametrização é possível comparar as inclinações. Basta ## escrever a matriz de funções lineares corretamente. L <- rbind(c(0,0,-1,1)) summary(glht(m2, linfct=L)) ##----------------------------------------------------------------------------- ## Predição. pred <- ddply(tur, .(Gender), summarise, length=extendseq(length, l=12)) pred <- cbind(pred, predict(m2, newdata=pred, interval="confidence")) p1 <- xyplot(height~length, groups=Gender, data=tur, xlab="Altura", ylab="Comprimento") p2 <- xyplot(fit~length, groups=Gender, data=pred, type="l", ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.5, prepanel=prepanel.cbH, panel=panel.superpose, panel.groups=panel.cbH) p1+as.layer(p2, under=TRUE) ## ##----------------------------------------------------------------------------- ## ## require(rpanel) ## ## form <- c("null"=.~1, ## "gend"=.~Gender, ## "leng"=.~length, ## "gend+leng"=.~Gender+length, ## "gend*leng"=.~Gender*length) ## ## pred <- ddply(tur, .(Gender), summarise, length=extendseq(length, l=12)) ## ## draw.model <- function(panel){ ## m3 <- update(m1, as.formula(form[[panel$f]])) ## pred <- cbind(pred, predict(m3, newdata=pred, interval="confidence")) ## p2 <- xyplot(fit~length, groups=Gender, data=pred, type="l", ## ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.5, ## prepanel=prepanel.cbH, panel=panel.superpose, ## panel.groups=panel.cbH) ## print(p1+as.layer(p2, under=TRUE)) ## panel ## } ## ## panel <- rp.control() ## rp.listbox(panel, variable=f, vals=names(form), action=draw.model) ## ##----------------------------------------------------------------------------- ## Total de tempo dormindo em função do peso corporal e nível de ## agressividade de mamíferos. str(sleep1) ## help(sleep1, h="html") ## Danger Index. sleep1$D <- factor(sleep1$D) sleep1$lbw <- log(sleep1$BodyWt) xyplot(TS~lbw|D, data=sleep1) xyplot(TS~lbw|D, data=sleep1, type=c("p","r")) ##----------------------------------------------------------------------------- ## Ajuste do modelo. m0 <- lm(TS~D*lbw, data=sleep1) ## Diagnóstico. par(mfrow=c(2,2)); plot(m0); layout(1) ## Relação média variância. ## Transformar a resposta? MASS::boxcox(m0) m1 <- update(m0, log(TS)~.) par(mfrow=c(2,2)); plot(m1); layout(1) ## Quadro de anova e de estimativas. anova(m1) summary(m1) ## summary(glht(m1)) summary(glht(m1), test=adjusted(type="fdr")) ##----------------------------------------------------------------------------- ## Estimativas por grupo. ## R² com SQT corrigida para média. summary(m1) m2 <- update(m1, .~0+D/lbw) ## R² com SQT **não** corrigida para a média. summary(m2) deviance(m1) deviance(m2) m3 <- update(m1, .~0+D+lbw) summary(m3) ## ##----------------------------------------------------------------------------- ## ## require(rpanel) ## ## form <- c("null"=.~1, ## "D"=.~D, ## "lbw"=.~lbw, ## "D+lbw"=.~D+lbw, ## "D*lbw"=.~D*lbw) ## ## pred <- ddply(sleep1, .(D), summarise, lbw=extendseq(lbw, l=12)) ## ## draw.model <- function(panel){ ## m3 <- update(m1, as.formula(form[[panel$f]])) ## pred <- cbind(pred, ## predict(m3, newdata=pred, interval="confidence")) ## p1 <- xyplot(log(TS)~lbw|D, data=sleep1) ## p2 <- xyplot(fit~lbw|D, data=pred, type="l", ## ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.25, ## prepanel=prepanel.cbH, ## panel=panel.cbH) ## print(p1+as.layer(p2, under=TRUE)) ## panel ## } ## ## panel <- rp.control() ## rp.listbox(panel, variable=f, vals=names(form), action=draw.model, ## title="Select a model") ## ##----------------------------------------------------------------------------- ## Preço dos carros em função da categoria (simples, sedan ou cross). url <- "http://www.leg.ufpr.br/~walmes/data/hb20_venda_webmotors_280314.txt" db <- read.table(url, header=TRUE, sep="\t") str(db) db <- subset(db, select=c("carro", "km", "preco")) db <- transform(db, km=km/1000, preco=preco/1000) ## Variável indicadora de carro zero km. ## db$novo <- ifelse(db$km==0, 1, 0) db$novo <- as.integer(db$km==0) xyplot(preco~km|carro, data=db) xyplot(preco~km|carro, groups=novo, data=db, type=c("p", "r")) ##----------------------------------------------------------------------------- ## Ajuste do modelo. m0 <- lm(preco~carro*(km+novo), data=db) ## Diagnóstico. par(mfrow=c(2,2)); plot(m0); layout(1) summary(m0) anova(m0) m1 <- update(m0, .~ carro*km) summary(m1) anova(m1) ## m2 <- update(m1, .~. -carro:km) m2 <- update(m1, .~carro+km) summary(m2) ##----------------------------------------------------------------------------- ## Predição. pred <- expand.grid(carro=levels(db$carro), km=seq(0,50,l=10)) aux <- predict(m2, newdata=pred, interval="confidence") pred <- cbind(pred, aux) p1 <- xyplot(preco~km, groups=carro, data=db, auto.key=TRUE) p2 <- xyplot(fit~km, groups=carro, data=pred, type="l", ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.5, panel=panel.superpose, panel.groups=panel.cbH) p1+as.layer(p2, under=TRUE) ##----------------------------------------------------------------------------- ## Dados de índice agronômico de milho. da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt", header=TRUE, sep="\t") str(da) xyplot(indice~dose|cultivar, data=da, type=c("p","a"), jitter.x=TRUE) ## xyplot(indice~dose|cultivar, data=da, groups=bloco, type="b") ## Por que sempre dois pontos são iguais? ##----------------------------------------------------------------------------- ## Ajuste do modelo saturado (equivalente ao poli grau 5). da$Dose <- factor(da$dose) m0 <- lm(indice~bloco+cultivar*Dose, data=da) ## Diagnóstico. par(mfrow=c(2,2)); plot(m0); layout(1) anova(m0) ## Ajuste dos modelos de trabalho. m1 <- update(m0, .~bloco+cultivar*(dose+I(dose^2)+I(dose^3))) m2 <- update(m0, .~bloco+cultivar*(dose+I(dose^2))) ## Teste da falta de ajuste. anova(m0, m1) anova(m0, m2) ##----------------------------------------------------------------------------- ## Predição (considerar o efeito de um bloco qualquer, o 1 primeiro). pred <- with(da, expand.grid(bloco=levels(bloco)[1], cultivar=levels(cultivar), dose=extendseq(dose))) pred <- cbind(pred, predict(m2, newdata=pred, interval="confidence")) p1 <- xyplot(indice~dose|cultivar, data=da) p2 <- xyplot(fit~dose|cultivar, data=pred, type="l", ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.25, prepanel=prepanel.cbH, panel=panel.cbH) p1+as.layer(p2) ## Não passa no meio dos pontos porque essa predição está ## particularizada para o efeito do bloco 1 e não na média dos blocos. ##----------------------------------------------------------------------------- ## Caso alguém tenha necessidade, é possível conduzir testes para saber ## se as inclinações são iguais na origem, ou seja, comparar o b1 das ## curvas. Só é necessário montar a matriz de contrastes. levels(da$cultivar) cbind(coef(m2)) L <- rbind("Ag vs BR"=c(0, 0,0,0, 0,0, 0, 0, 1,0, 0,0), "Ag vs Pi"=c(0, 0,0,0, 0,0, 0, 0, 0,1, 0,0), "Br vs Pi"=c(0, 0,0,0, 0,0, 0, 0,-1,1, 0,0)) ## Estimativa das diferenças em inclinação na origem. L%*%coef(m2) summary(glht(m2, linfct=L)) ## A produtividade na dose 240 é a mesma nessas cultivares? Como fazer ## esse teste de médias? LSmeans(m2, effect="cultivar", at=list("dose"=240)) L <- LSmatrix(m2, effect="cultivar", at=list("dose"=240)) levels(da$cultivar) L <- rbind( "Ag"=LSmatrix(m2, at=list("cultivar"="Ag-1002", "dose"=320)), "Br"=LSmatrix(m2, at=list("cultivar"="BR-300", "dose"=250)), "Pi"=LSmatrix(m2, at=list("cultivar"="Pioneer-B815", "dose"=255))) ## Estimativas das produtividades nessa dose. summary(glht(m2, linfct=L)) ## Estimativas dos contrastes. summary(glht(m2, linfct=apc(L))) summary(glht(m2, linfct=apc(L)), test=Ftest()) ##----------------------------------------------------------------------------- ## Como obter as equações (se é que fato são necessárias). m3 <- update(m2, .~0+cultivar/(dose+I(dose^2))+bloco, contrasts=list(bloco=contr.sum)) summary(m3)$coeff split(coef(m3), m3$assign) ##----------------------------------------------------------------------------- ## É possível também, caso considerar necessário, repartir as somas de ## quadrados. m4 <- aov(indice~bloco+cultivar/(dose+I(dose^2)), data=da) m4$assign anova(m4) coef(m4) summary(m4, split=list( "cultivar:dose"=list(Ag=1, BR=2, Pi=3), "cultivar:I(dose^2)"=list(BR=2, Pi=3, Ag=1) )) m5 <- aov(indice~bloco+cultivar/(dose+I(dose^2)+I(dose^3)), data=da) summary(m5, split=list( "cultivar:dose"=list(Ag=1, BR=2, Pi=3), "cultivar:I(dose^2)"=list(Ag=1, BR=2, Pi=3), "cultivar:I(dose^3)"=list(Ag=1, BR=2, Pi=3) )) levels(da$cultivar) da$ind <- as.integer(da$cultivar=="Pioneer-B815") m6 <- lm(indice~bloco+cultivar*(dose+I(dose^2))+ind:I(dose^3), data=da) summary(m6) ## As mesmas conclusões são obtidas se olhar para o quadro de ## estimativas. summary.lm(m4)$coeff ##----------------------------------------------------------------------------- ## Obter os valores preditos na média dos blocos. pred <- with(da, expand.grid(bloco=levels(bloco)[1], cultivar=levels(cultivar), dose=extendseq(dose))) pred$ind <- as.integer(pred$cultivar=="Pioneer-B815") str(pred) pred <- equallevels(pred, da) m2 <- m6 L <- model.matrix(formula(m2)[-2], data=pred) head(L) ## Nas colunas de correspondem ao efeito dos blocos, colocar 0.25 (1/4). m2$assign ## Coeficientes de bloco são indicados por 1. L[,m2$assign==1] <- L[,m2$assign==1]+0.25 head(L) ## Obter os ic com a glht. ic <- confint(glht(m2, linfct=L), calpha=univariate_calpha())$confint pred <- cbind(pred, ic) p3 <- xyplot(Estimate~dose|cultivar, data=pred, type="l", ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.25, col=2, fill=2, prepanel=prepanel.cbH, panel=panel.cbH) p1+as.layer(p2)+as.layer(p3) ## As bandas de confiança no efeito ponderado dos blocos são visualmente ## apelativas pois são mais precisas e passam entre os pontos. ##----------------------------------------------------------------------------- ## Dados de experimento com soja sob 3 níveis de umidade do solo (massa ## de água/massa de solo) e níveis de potássio fornecidos na adubação. A ## hipótese sob estudo é que presença de potássio dá suporte para a ## cultura superar a ocorrência de períodos sem chuva. Experimento ## conduzido em casa de vegetação com 5 blocos, 2 plantas por ## u.e. (vaso). ## Para acessar o artigo (pdf online). ## browseURL("http://www.scielo.br/pdf/rca/v43n2/a03v43n2.pdf") da <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt", header=TRUE, sep="\t", dec=",") names(da) <- substr(names(da), 1, 4) str(da) xyplot(reng~pota, groups=agua, data=da, type=c("p","a")) xyplot(reng~pota|agua, data=da, type=c("p","a")) ##----------------------------------------------------------------------------- ## Ajuste do modelo. da$A <- factor(da$agua) da$K <- factor(da$pota) ## Ajuste. m0 <- lm(reng~bloc+A*K, data=da) par(mfrow=c(2,2)); plot(m0); layout(1) ## Remove obs 74 (uma planta que definou). da <- da[-74,] m0 <- update(m0, data=da) par(mfrow=c(2,2)); plot(m0); layout(1) ## Quadro de anova. anova(m0) ## Estimativas dos parâmentros (efeitos). summary(m0) ##----------------------------------------------------------------------------- ## Polinômio de segundo grau para o efeito de potássio. ## m1 <- lm(reng~bloc+A*(pota+I(pota^2)), data=da) m1 <- update(m0, .~bloc+A*(pota+I(pota^2))) par(mfrow=c(2,2)); plot(m1); layout(1) ## Testa o abandono dos termos (falta de ajuste). anova(m0, m1) anova(m1) ## Com b2 comum aos níveis de água. ## m2 <- lm(reng~bloc+A*pota+I(pota^2), data=da) m2 <- update(m1, .~bloc+A*pota+I(pota^2)) anova(m2, m0) ##----------------------------------------------------------------------------- ## Predição com bandas (considerando efeito de bloco I). pred <- expand.grid(bloc="I", A=levels(da$A), pota=seq(0, 180, by=3)) aux <- predict(m2, newdata=pred, interval="confidence") pred <- cbind(pred, aux) str(pred) ## Melhorar a legenda. tx <- paste("Umidade de ", levels(da$A), "%", sep="") lgd <- simpleKey(columns=1, text=tx, points=FALSE, lines=TRUE, corner=c(0.03,0.97)) p1 <- xyplot(reng~pota, groups=A, data=da) p2 <- xyplot(fit~pota, groups=A, data=pred, type="l", key=lgd, xlab=expression("Potássio no solo"~(mg~dm^{-3})), ylab=expression("Rendimento de grãos"~(g~vaso^{-1})), ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.1, fill=1, prepanel=prepanel.cbH, panel.groups=panel.cbH, panel=panel.superpose) p1+as.layer(p2) ## Como calcular os pontos de máximo? São dados por x_max = -0.5*b1/b2. ##----------------------------------------------------------------------------- ## Testar se o rendimento sem potássio (0) é o mesmo para as águas. X <- LSmatrix(m2, effect="A", at=list("pota"=0, "I(pota^2)"=0)) rownames(X) <- levels(da$A) X g <- glht(m2, linfct=X) confint(g, calpha=univariate_calpha()) ## Contrastes par a par. Xc <- apc(X) summary(glht(m2, linfct=Xc)) summary(glht(m2, linfct=Xc), test=Ftest()) ## O teste simultâneo, intercepto comum. formula(m1) m3 <- update(m2, formula=.~.-A) coef(m3) anova(m3, m2) anova(m2) ##----------------------------------------------------------------------------- ## Coeficientes das equações. Mudar a formula do modelo para ter uma ## interpretação de parâmetros por estrato. ## Usar contraste tipo soma para blocos. formula(m2) m3 <- lm(reng~0+A/pota+I(pota^2)+bloc, data=da, contrasts=list(bloc=contr.sum)) summary(m3) ## Prova de que são o mesmo modelo. deviance(m3); deviance(m2) ## O R² é o quadrado da correlação entre observado e predito. cor(da$reng, fitted(m2))^2 ## R² separado por água. ddply(data.frame(o=da$reng, f=fitted(m2)), .(da$A), summarise, R2=cor(o, f)^2) ##----------------------------------------------------------------------------- ## Os pontos de máximo. ## Buscas essas ocorrências de texto. oc <- paste0("^A", levels(da$A), ":") ## Estimativas e suas posições de ordem de efeito. c3 <- coef(m3) a <- m3$assign cbind(c3, a) ## Estimativas do ponto de máximo. ptmax <- -0.5*c3[a==4]/c3[a==2] xyplot(fit~pota, groups=A, data=pred, type="l", key=lgd, xlab=expression("Potássio no solo"~(mg~dm^{-3})), ylab=expression("Rendimento de grãos"~(g~vaso^{-1})), ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.1, fill=1, prepanel=prepanel.cbH, panel.groups=panel.cbH, panel=panel.superpose)+ layer( panel.abline(v=ptmax, lty=2)) ##----------------------------------------------------------------------------- ## E se ajustar um modelo linear-platô? E se for um quadrático-platô? ##----------------------------------------------------------------------------- ## Dados. da <- read.table("http://www.leg.ufpr.br/~walmes/data/algodão.txt", header=TRUE, sep="\t", encoding="latin1") da$desf <- da$desf/100 levels(da$estag) <- c("Vegetativo", "Botão floral", "Florescimento", "Maçã", "Capulho") str(da) xyplot(pcapu~desf|estag, data=da, type=c("p","smooth"), xlab="Nível de desfolha artificial", ylab="Peso de capulhos produzidos (g)") ##----------------------------------------------------------------------------- ## 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. ##----------------------------------------------------------------------------- ## Especificação dos modelos. ## 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) ## 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")) Xc <- 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)) ##----------------------------------------------------------------------------- ## Ler tabela de dados. da <- read.table("http://www.leg.ufpr.br/~walmes/data/maracuja_planta.txt", header=TRUE, sep="\t") da$bloc <- factor(da$bloc) str(da) u <- unique(da$cinz) cbind(u, u/1.5, sqrt(u/1.5), log(u/1.5, base=2)) ## Doses de cinza em uma escala com distância mais regular entre níveis. da$cin <- sqrt(da$cinz/1.5) ##----------------------------------------------------------------------------- ## Ver. p1 <- xyplot(mfpa~cin|fam, groups=agreg, data=da, type=c("p","smooth"), auto.key=list(columns=2)) p2 <- xyplot(mspa~cin|fam, groups=agreg, data=da, type=c("p","smooth"), auto.key=list(columns=2)) p3 <- xyplot(ds~cin|fam, groups=agreg, data=da, type=c("p","smooth"), auto.key=list(columns=2)) p4 <- xyplot(cav~cin|fam, groups=agreg, data=da, type=c("p","smooth"), auto.key=list(columns=2)) grid.arrange(p1, p2, nrow=2) grid.arrange(p3, p4, nrow=2) ##----------------------------------------------------------------------------- ## Especificação e ajuste dos modelos em batelada. resp <- c("mfpa", "mspa", "ds", "cav") form0 <- lapply(paste(resp, "~bloc+fam*agreg*(cinz+I(cinz^2))"), as.formula) names(form0) <- resp ## Ajustes. m0 <- lapply(form0, lm, data=da) ##----------------------------------------------------------------------------- ## Quadros de anova. lapply(m0, anova) ##----------------------------------------------------------------------------- ## Quadro de estimativas dos efeitos. lapply(m0, summary) ##----------------------------------------------------------------------------- ## Avaliação dos pressupostos. par(mfrow=c(2,2)); plot(m0[["mfpa"]]); layout(1) ## ok! par(mfrow=c(2,2)); plot(m0[["mspa"]]); layout(1) ## ok! par(mfrow=c(2,2)); plot(m0[["ds"]]); layout(1) ## bom. par(mfrow=c(2,2)); plot(m0[["cav"]]); layout(1) ## ok! ## par(mfrow=c(2,2)) ## lapply(m0, plot, which=1) ## layout(1) ## par(mfrow=c(2,2)) ## lapply(m0, plot, which=2) ## layout(1) ## par(mfrow=c(2,2)) ## lapply(m0, plot, which=3) ## layout(1) ##----------------------------------------------------------------------------- ## Modelos após abandono de termos não significativos. form1 <- list(mfpa=mfpa~bloc+agreg*cinz, mspa=mspa~bloc+agreg*cinz, ds=ds~bloc+(fam+agreg+cinz)^3+(fam+agreg)*I(cinz^2), cav=cav~bloc+fam*agreg*(cinz+I(cinz^2))) m1 <- lapply(form1, lm, data=da, contrast=list(bloc=contr.sum)) lapply(m1, anova) ## Anova entre modelos sequenciais. sapply(names(m1), simplify=FALSE, function(i){ anova(m0[[i]],m1[[i]]) }) ##----------------------------------------------------------------------------- ## summary(m0[["mfpa"]]) ##----------------------------------------------------------------------------- ## Valores preditos. gridlist <- list(bloc="1", fam=levels(da$fam), agreg=levels(da$agreg), cin=seq(0,6,l=100), cinz=seq(0,50,l=100)) m1 <- lapply(m1, function(i){ predvars <- attr(terms(i), "term.labels") pred <- do.call( expand.grid, gridlist[predvars[!sapply(gridlist[predvars], is.null)]]) i$newdata <- pred i }) mypredict <- function(m){ cbind(m$newdata, predict(m, newdata=m$newdata, interval="confidence")- coef(m)["bloc1"]) } all.pred <- lapply(m1, mypredict) ## str(all.pred) ## lapply(all.pred, names) ##----------------------------------------------------------------------------- ## Gráficos. xlab <- expression(Cinza~(t~ha^{-1})) ylab <- list(expression(Massa~fresca~de~parte~aérea~(g)), expression(Massa~seca~de~parte~aérea~(g)), expression(Densidade~do~solo~(Mg~t^{-1})), expression(Água~consumida~no~ciclo~(mL))) names(ylab) <- names(m1) ##----------------------------------------------------------------------------- ## Mfpa. p1 <- xyplot(mfpa~cinz, groups=agreg, data=da, ylab=ylab[["mfpa"]], xlab=xlab, strip=strip.custom(bg="gray90"), auto.key=list(columns=2, points=FALSE, lines=TRUE, title="Classe de agregado (mm)", cex.title=1.2))+ as.layer(with(all.pred[["mfpa"]], xyplot(fit~cinz, groups=agreg, type="l", ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2, prepanel=prepanel.cbH, panel=panel.superpose, panel.groups=panel.cbH))) ##----------------------------------------------------------------------------- ## Mspa. p2 <- xyplot(mspa~cinz, groups=agreg, data=da, ylab=ylab[["mspa"]], xlab=xlab, strip=strip.custom(bg="gray90"), auto.key=list(columns=2, points=FALSE, lines=TRUE, title="Classe de agregado (mm)", cex.title=1.2))+ as.layer(with(all.pred[["mspa"]], xyplot(fit~cinz, groups=agreg, type="l", ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2, prepanel=prepanel.cbH, panel=panel.superpose, panel.groups=panel.cbH))) ##----------------------------------------------------------------------------- ## Ds. p3 <- xyplot(ds~cinz|fam, groups=agreg, data=da, layout=c(NA,1), ylab=ylab[["ds"]], xlab=xlab, strip=strip.custom(bg="gray90"), auto.key=list(columns=2, points=FALSE, lines=TRUE, title="Classe de agregado (mm)", cex.title=1.2))+ as.layer(with(all.pred[["ds"]], xyplot(fit~cinz|fam, groups=agreg, type="l", ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2, prepanel=prepanel.cbH, panel=panel.superpose, panel.groups=panel.cbH))) ##----------------------------------------------------------------------------- ## Cav. p4 <- xyplot(cav~cinz|fam, groups=agreg, data=da, layout=c(NA,1), ylab=ylab[["cav"]], xlab=xlab, strip=strip.custom(bg="gray90"), auto.key=list(columns=2, points=FALSE, lines=TRUE, title="Classe de agregado (mm)", cex.title=1.2))+ as.layer(with(all.pred[["cav"]], xyplot(fit~cinz|fam, groups=agreg, type="l", ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2, prepanel=prepanel.cbH, panel=panel.superpose, panel.groups=panel.cbH))) ##----------------------------------------------------------------------------- ## Arranjo 1. grid.arrange(p1, p2, ncol=2) grid.text(label="A", x=0.025, y=0.975) grid.text(label="B", x=0.5+0.025, y=0.975) ##----------------------------------------------------------------------------- ## Arranjo 2. grid.arrange(p3, p4, ncol=2) grid.text(label="A", x=0.025, y=0.975) grid.text(label="B", x=0.5+0.025, y=0.975) ##----------------------------------------------------------------------------- ## Informações da sessão. Sys.time() sessionInfo()