| This post was kindly contributed by Ridículas - go there to comment and to read the full post. |
Valores observados e ajustados acompanhados do intervalo de 95% em função da dose de potássio e nível de água para as variáveis resposta estudadas.
Na realização de experimento é comum a observação de várias variáveis respostas em uma mesma unidade amostral. Em termos estatísticos, o que observamos é um vetor aleatório. Por exemplo, em experimento de com soja para fins de ensaio de competição de cultivares o pesquisador pode observar: altura da planta, altura da primeira vagem, massa seca da planta, número de vagens por planta, teores foliares de N, P, K, Mg, Ca, peso de 100 grãos, proporção de sementes viáveis, dias para fechar o ciclo, todas elas observadas em cada unidade experimental. Nesse caso então, observamos um vetor aleatório de 12 variáveis aleatórias. O ideal seria empregarmos um modelo multivariado, porém como as respostas podem ser contínuas não limitadas, contagens, contínuas limitadas, contínuas com censura, proporções, etc, ainda não dispomos de um modelo de distribuição multivariado que contemple essa situação. O mais comum é fazer a análise univariada e como todas essas variáveis estão sob o efeito das mesmas fontes de variação, a análise estatística vai empregar o mesmo modelo (pelo menos no tocante a parte determinística). Sendo assim, tem como fazer o estudo de todas as variáveis reposta de uma vez só? Sim. Para saber como continue lendo.
Existem diversas formas de aplicar tarefas em grupos ou colunas. A família apply do stats oferece recurso para um número bem grande de situações. A família ply do plyr acrescenta mais opções. Além destes temos o doBy com mais algumas funcionalidades.
Bem, vamos ao que interessa. No CMR eu faço ajustes de regressão com polinômios para dados de um experimento com a cultura da soja (vasos em casa de vegetação). Foram estudados o nível de água (umidade do solo em relação à capacidade de campo) e potássio aplicado ao solo. O objetivo desse experimento foi verificar se o fornecimento de potássio pode minimizar impactos do deficit hídrico. Para acessar o artigo completo clique aqui: Umidade do solo e doses de potássio na cultura da soja. Não deixe de conferir a lista de links abaixo. Em próximos artigos do Rídiculas vamos aprofundar o estudo nesses dados comparando curvas. Até lá.
#-----------------------------------------------------------------------------
# importa dados
soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
header=TRUE, sep="\t", dec=",")
str(soja)
soja <- transform(soja, a=agua, A=factor(agua),
k=potassio, K=factor(potassio))
resp <- 4:7 # índice das variáveis resposta
str(soja)
#-----------------------------------------------------------------------------
# ver
soja2 <- cbind(stack(soja[,resp]), soja[,-resp])
require(lattice)
xyplot(values~k|ind, groups=A, data=soja2,
scales="free", jitter.x=TRUE, type=c("p","a"))
#-----------------------------------------------------------------------------
# ajuste do modelo saturado de efeitos categóricos (modelo de médias)
m0 <- apply(soja[,resp], 2,
function(r){ aov(r~bloco+A*K, soja) })
lapply(m0, anova)
#-----------------------------------------------------------------------------
# ajuste de polinômio para potassio por nível de água
m1 <- apply(soja[,resp], 2,
function(r){ aov(r~bloco+A/poly(k, 4), soja) })
lapply(m1, anova)
names(coef(m1[[1]]))
inter <- m1[[1]]$assign==3 # indices das interações
eff <- names(coef(m1[[1]]))[inter]
togrep <- outer(c("A37.5","A50","A62.5"), c("1$","2$","(3|4)$"),
paste, sep=":.*") # nomes de busca
or <- order(togrep)
tolist <- outer(c("A37.5","A50.0","A62.5"), c("lin","qua","lof"),
paste, sep="|") # nomes de rótulo
list.desd <- sapply(togrep, grep, eff)
names(list.desd) <- tolist
list.desd <- list.desd[or]
list.desd # lista com os índices dos desdobramentos
# anovas com desdobramento de efeito de potássio em cada nível de água
lapply(m1, summary, split=list("A:poly(k, 4)"=list.desd))
#-----------------------------------------------------------------------------
# com base nas anovas acima selecionamos os modelos
form <- list( rengrao~bloco+A/poly(k, 2, raw=TRUE),
pesograo~bloco+A/poly(k, 2, raw=TRUE),
kgrao~bloco+A/poly(k, 2, raw=TRUE),
pgrao~bloco+A/poly(k, 3, raw=TRUE))
names(form) <- names(soja)[resp]
m2 <- lapply(form, aov, data=soja, contrasts=list(bloco=contr.sum))
lapply(m2, coef)
#-----------------------------------------------------------------------------
# predição a partir do modelo
str(soja)
pred <- expand.grid(bloco="I", A=levels(soja$A), k=seq(0,180,5))
#aux <- lapply(m2, predict, newdata=pred, interval="confidence")
aux <- lapply(m2, function(r){ # remove o efeito do bloco na predição
predict(r, newdata=pred, interval="confidence")-coef(r)["bloco1"]})
aux <- lapply(aux, as.data.frame)
str(aux)
require(plyr)
pred <- cbind(ldply(aux), pred)
pred$.id <- factor(pred$.id, levels=levels(soja2$ind))
#-----------------------------------------------------------------------------
# gráfico final
fl <- c("K nos grãos (mg)", "Peso 100 grãos (g)",
"P nos grãos (mg)", "Rendimento de grãos (g)")
cbind(levels(pred$.id), levels(soja2$ind), fl)
require(latticeExtra)
display.brewer.all()
mycol <- brewer.pal(6, "Set1")
panel.ciH <- function(x, y, ly, uy, subscripts, ...){
y <- as.numeric(y)
x <- as.numeric(x)
or <- order(x)
ly <- as.numeric(ly[subscripts])
uy <- as.numeric(uy[subscripts])
panel.polygon(c(x[or], x[rev(or)]), c(ly[or], uy[rev(or)]),
col=1, alpha=0.05, border=NA)
panel.lines(x[or], ly[or], lty=3, lwd=0.5, col=1)
panel.lines(x[or], uy[or], lty=3, lwd=0.5, col=1)
panel.xyplot(x, y, subscripts=subscripts, ...)
}
#png(file="f031.png", width=500, height=500)
trellis.par.set(superpose.symbol=list(col=mycol),
superpose.line=list(col=mycol),
strip.background=list(col="gray90"))
xyplot(values~k|ind, groups=A, data=soja2,
scales="free", jitter.x=TRUE,
xlab="Potássio no solo", ylab=NULL,
strip=strip.custom(factor.levels=fl),
auto.key=list(title="Níveis de água (%)", cex.title=1.2,
columns=3, points=FALSE, lines=TRUE))+
as.layer(xyplot(fit~k|.id, groups=A, pred, type="l", lwd=1.5,
ly=pred$lwr, uy=pred$upr,
panel.groups=panel.ciH,
panel=panel.superpose))
#dev.off()
#-----------------------------------------------------------------------------
# pode-se tentar comparar as curvas ajustando modelos apropriados sob H0
# exemplo, será que as curvas vermelha e verde para kgrao e pesograo não
# podem ser representados pela mesma curva?
# as curvas azul e ver para pgrao diferem de uma reta horizontal?
#-----------------------------------------------------------------------------
