Ajuste simultâneo de regressão/anova para várias respostas

April 23, 2012
By

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á.

  • A brief introduction to “apply” in R;
  • plyr – The split-apply-combine strategy for R;
  • plyr: Tools for splitting, applying and combining data;
  • R Grouping functions: sapply vs. lapply vs. apply. vs. tapply vs. by vs. aggregate vs.;
  • Say it in R with “by”, “apply” and friends;
  • The doBy package;
  • doBy: doBy – Groupwise summary statistics, general linear contrasts, population means (least-squares-means), and other utilities;
  • #-----------------------------------------------------------------------------
    # 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?
    
    #-----------------------------------------------------------------------------
    


    Tags: , , , ,

    Comments are closed.