Archive

Posts Tagged ‘plyr’

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

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?
    
    #-----------------------------------------------------------------------------
    
    Categorias:regressão Tags:, , ,

    Média e desvio-padrão de muitas variáveis separado por grupos

    Gráfico de barras com as médias e desvios-padrões para variáveis resposta de um experimento avaliando cultivares de banana e meios de cultura em cultura de tecidos.

    Em estatística descritiva é comum a síntese dos dados em medidas de posição e dispersão. As medidas mais frequentemente usadas são a média e o desvio-padrão. Não raramente, existe a necessidade de obter essas medidas separando as observações por alguma variável categórica do conjunto de dados. Nesse post apresento algumas formas de obter isso. Estou certo de que as soluções que apresento abaixo não cobrem 5% das disponíveis considerando a extensa quantidade de pacotes/funções do R ligadas à análise descritiva, portanto, considere isso apenas como um ponto de partida. Por fim, sintetizo em um gráfico os resultamos aplicados em um conjunto de dados reais. Até a próxima ridícula.

    #-----------------------------------------------------------------------------
    # gerando um conjunto de dados artificial
    
    da <- expand.grid(A=1:3, B=1:2, r=1:15)
    da <- cbind(da, as.data.frame(matrix(rnorm(nrow(da)*3), nrow(da))))
    str(da)
    
    #-----------------------------------------------------------------------------
    # usando funções do reshape e plyr
    
    require(reshape)
    require(plyr)
    
    ddply(da[,-c(1:3)], .(da$A, da$B), mean) # as médias
    ddply(da[,-c(1:3)], .(da$A, da$B), sd)   # os desvios-padrões
    
    db <- melt(da, id.vars=1:3, id.measures=4:6)
    str(db)
    
    #-----------------------------------------------------------------------------
    # fazendo empilhamento das observações
    
    ddply(db[,-3], .(A, B, variable), summarise, m=mean(value), sd=sd(value))
    
    #-----------------------------------------------------------------------------
    # usando funções do pacote doBy
    
    library(doBy)
    
    summaryBy(V1+V2+V3~A+B, data=da,
              FUN=function(x){ c(m=mean(x), sd=sd(x)) })
    
    #-----------------------------------------------------------------------------
    # usando funções do Hmisc
    
    require(Hmisc)
    
    with(da, bystats(cbind(V1, V2, V3), A, B,
                     fun=function(x){ c(m=colMeans(x), s=sd(x)) }))
    
    #-----------------------------------------------------------------------------
    # dados reais
    
    bnn <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/cnpaf/banana.txt",
                      header=TRUE, sep="\t")
    str(bnn)
    bnn <- melt(bnn, id.measures=4:12)
    bnn <- ddply(bnn, .(cult, meio, variable),
                 summarise, m=mean(value), sd=sd(value))
    str(bnn)
    
    #-----------------------------------------------------------------------------
    # gráfico de barras com desvio-padrão
    
    require(lattice)
    require(RColorBrewer)
    display.brewer.all()
    cols <- brewer.pal(nlevels(bnn$meio), "Set3")
    
    png("f028.png", 500, 500)
    trellis.par.set(superpose.polygon=list(col=cols),
                    strip.background=list(col="gray90"))
    barchart(m~cult|variable, groups=meio, data=bnn,
             sdv=bnn$sd, scales="free", layout=c(3,3),
             xlab="Cultivar", ylab=expression(bar(x)+sd),
             key=list(title="Meio de cultura", cex.title=1.2,
               columns=3, text=list(levels(bnn$meio)),
               rectangle=list(col=trellis.par.get("superpose.polygon")$col)),
             prepanel=function(y, sdv, subscripts=subscripts, ...){
               ly <- as.numeric(y+sdv[subscripts])
               list(ylim=range(y, ly, finite=TRUE))
             },
             panel=function(x, y, subscripts, groups, sdv, box.ratio, ...){
               panel.barchart(x, y, subscripts=subscripts,
                              groups=groups, box.ratio=box.ratio, ...)
               d <- 1/(nlevels(groups)+nlevels(groups)/box.ratio)
               g <- (as.numeric(groups[subscripts])-1); g <- (g-median(g))*d
               panel.arrows(as.numeric(x)+g, y, as.numeric(x)+g, y+sdv[subscripts],
                            angle=90, length=0.025)
             })
    dev.off()
    
    #-----------------------------------------------------------------------------
    
    Categorias:gráficos Tags:, , ,

    Converter dados em formato amplo e formato longo

    Respostas observadas em função dos fatores em experimento com a cultura da soja.

    A coisa mais comum para quem analisa dados e ter que editar o formato dos dados. Dificilmente você recebe os dados no formato correto. Muitas pessoas fazem de suas planilhas de dados uma arte, com cédulas coloridas, linhas e colunas com textos em diversos formatos.

    Dados de medidas repetidas tem a característica de cada nova medida repetida ser uma nova coluna na planilha (formato amplo). A maioria dos métodos de análise requer que os dados estejam no formato longo, ou seja, cada linha é um registro, cada coluna é uma variável. Portanto, temos que converter o formato dos dados. Isso vai se tornar claro/simples com os comandos que apresento nessa ridícula.

    Os pacotes reshape e plyr contém funções úteis para a manipulação/restruturação dos dados. Os comandos abaixo importam dados no formato amplo, passam para o formato longo e retornam para o formato amplo.

    > #-----------------------------------------------------------------------------
    > # dados de espessura de fibra de Schinus
    > 
    > require(reshape)
    > 
    > # dados estão no formato amplo, colunas são secções no tronco da árvore
    > fibra <- read.table("http://www.leg.ufpr.br/~walmes/ridiculas/fibra.txt",
    +                     header=TRUE, sep="\t")
    > str(fibra)
    'data.frame':	150 obs. of  7 variables:
     $ arvore : int  1 1 1 1 1 1 1 1 1 1 ...
     $ amostra: int  1 2 3 4 5 6 7 8 9 10 ...
     $ A      : num  2.61 2.43 2.09 2.5 2.06 ...
     $ B      : num  3.17 2.68 2.43 2 1.85 ...
     $ C      : num  2.17 4.19 4.82 5.54 3.26 ...
     $ D      : num  3.7 2.2 3.5 2.61 3.43 ...
     $ E      : num  2.4 3.5 3.11 2.4 3.26 ...
    > 
    > # passamos os dados para o formato longo, criou uma coluna para secção
    > fibra <- melt(fibra, id=c("arvore", "amostra"))
    > str(fibra)
    'data.frame':	750 obs. of  4 variables:
     $ arvore  : int  1 1 1 1 1 1 1 1 1 1 ...
     $ amostra : int  1 2 3 4 5 6 7 8 9 10 ...
     $ variable: Factor w/ 5 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
     $ value   : num  2.61 2.43 2.09 2.5 2.06 ...
    > 
    > # retornamos os dados para o formato amplo
    > fibra <- adply(cast(fibra, amostra~arvore~variable,
    +                     value="value", fun=NULL), c(1,2))
    > str(fibra)
    'data.frame':	150 obs. of  7 variables:
     $ amostra: chr  "1" "2" "3" "4" ...
     $ arvore : chr  "1" "1" "1" "1" ...
     $ A      : num  2.61 2.43 2.09 2.5 2.06 ...
     $ B      : num  3.17 2.68 2.43 2 1.85 ...
     $ C      : num  2.17 4.19 4.82 5.54 3.26 ...
     $ D      : num  3.7 2.2 3.5 2.61 3.43 ...
     $ E      : num  2.4 3.5 3.11 2.4 3.26 ...
    > 
    > #-----------------------------------------------------------------------------
    

    Os dados do experimento em soja já estão no formato longo, exigido para fazer ajuste de modelos lineares. Nesse experimento foram observadas quatro variáveis resposta. Para fazer uma análise exploratória conjunta, eu passei os dados para o formato longo. Assim, tenho uma coluna que identifica a resposta observada e outra o valor. Com os recursos do pacote lattice eu fiz um gráfico de dispersão para tirar informações iniciais dos resultados do experimento.

    #-----------------------------------------------------------------------------
    # dados
    
    soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                       header=TRUE, sep="\t", dec=",")
    names(soja)[4:7] <- c("Rendimento de grãos (g/vaso)",
                          "Massa de 100 grãos (g)",
                          "K no grão (g/kg)", "P no grão (g/kg)")
    str(soja)
    
    #-----------------------------------------------------------------------------
    # criando uma coluna que identifica as variáveis resposta
    
    soja2 <- melt(soja, id=1:3)
    str(soja2)
    
    #-----------------------------------------------------------------------------
    # gráfico para análise exploratória
    
    require(lattice)
    
    #png("f006.png", w=500, h=400)
    mycols <- colors()[c(47,81,614)]
    trellis.par.set(superpose.symbol=list(col=mycols, pch=c(15,17,19)),
                    superpose.line=list(col=mycols),
                    strip.background=list(col="gray90"))
    xyplot(value~potassio|variable, groups=agua, data=soja2,
           xlab="Potássio no solo (mg/dm³)", ylab="Valor observado",
           key=list(title="Níveis de umidade do solo (%)", cex.title=1,
             columns=3, type="o", divide=1,
             lines=list(pch=trellis.par.get("superpose.symbol")$pch,
               col=trellis.par.get("superpose.symbol")$col),
             text=list(c(format(unique(soja$agua))))),
           type=c("p","a","g"), jitter.x=TRUE,
           scales=list(y=list(relation="free")), layout=c(2,2))
    #dev.off()
    
    #-----------------------------------------------------------------------------
    # para ver as configurações gráficas da lattice
    
    library(TeachingDemos)
    TkListView(trellis.par.get())
    
    #-----------------------------------------------------------------------------
    
    Categorias:dados, gráficos Tags:, , , , ,

    Fatorial com desdobramento da interação, gráfico e tabela de médias

    Rendimento de grãos em função dos níveis de potássio e umidade do solo. Letras minúsculas comparam níveis de potassio em um mesmo nível de umidade, as letras maiúsculas fazem o inverso.

    O R não é um aplicativo específico de análise de experimentos. Mesmo assim, ele oferece funções que permitem a análise dos mais diversos tipos de experimentos, ficando a cargo do usuário saber usar corretamente as funções/pacotes que o R oferece.

    Pode-se dizer que os experimentos fatoriais são os mais comuns na ciência. Nas ciências agrárias são amplamente usados. Têm-se uma tradição de representar a interação entre os fatores (qualitativos) empregando testes de médias para níveis de um fator fixando níveis dos demais.

    Nessa ridícula, apresento procedimentos para obter o desdobramento de interação em um fatorial duplo. Os dados são do rendimento de grãos de soja em função dos fatores umidade do solo e dose de fertilizante potássico. O experimento foi instalado em casa de vegetação com as unidades experimentais (vasos) agrupadas em blocos que controlaram variações de luminosidade/umidade/ventilação.

    Os procedimentos incluem o ajuste do modelo, o teste sobre os efeitos dos fatores pela análise de variância, a checagem das pressuposições do modelo, o fatiamento das somas de quadrados, o teste de médias. No final são apresentados procedimentos para representar os resultados em gráfico e tabela de dupla entrada. É importante dizer que os fatores deveriam ser analisados como quantitativos, onde espera-se uma interpretação mais interessante sobre o fenômeno. Entretanto, o método abaixo é ilustrativo do procedimento para fatores qualitativos. Em uma futura oportunidade analisarei dos dados considerando os fatores como quantitativos.

    #-----------------------------------------------------------------------------
    # importa e prepara os dados
     
    rend <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/rendimento.txt",
                       header=TRUE, sep="\t")
    rend <- transform(rend, K=factor(K), A=factor(A), bloc=factor(bloc))
    str(rend)
     
    #-----------------------------------------------------------------------------
    # análise gráfica
     
    require(lattice)
    xyplot(rg~K|A, groups=bloc, data=rend, type="b", auto.key=TRUE)
     
    #-----------------------------------------------------------------------------
    # ajuste do modelo e faz checagem
     
    m0 <- aov(rg~bloc+A*K, data=rend)
    summary(m0)
    par(mfrow=c(2,2)); plot(m0); layout(1)
     
    #-----------------------------------------------------------------------------
    # desdobrando somas de quadrados para a variação de K dentro de A e vice-versa
    
    m1 <- aov(rg~bloc+K/A, data=rend)
    desAinK <- sapply(paste("K", levels(rend$K), sep=""), simplify=FALSE,
                      grep, x=names(coef(m1)[m1$assign==3]))
    summary(m1, split=list("K:A"=desAinK))
    
    m2 <- aov(rg~bloc+A/K, data=rend)
    desKinA <- sapply(paste("A", levels(rend$A), sep=""), simplify=FALSE,
                      grep, x=names(coef(m2)[m2$assign==3]))
    summary(m2, split=list("A:K"=desKinA))
     
    #-----------------------------------------------------------------------------
    # desdobrando a interação em testes de Tukey
     
    require(agricolae)
    require(plyr)
    
    KinA <- sapply(levels(rend$A), simplify=FALSE,
                   function(a){
                     with(subset(rend, A==a),
                          HSD.test(rg, K,
                                   DFerror=df.residual(m0),
                                   MSerror=deviance(m0)/df.residual(m0)))
                   })
    
    KinA <- ldply(KinA, NULL)
    KinA$M <- gsub(" ", "", KinA$M, fixed=TRUE)
    KinA$trt <- as.factor(as.numeric(as.character(KinA$trt)))
    str(KinA)
    
    AinK <- sapply(levels(rend$K), simplify=FALSE,
                   function(k){
                     with(subset(rend, K==k),
                          HSD.test(rg, A,
                                   DFerror=df.residual(m0),
                                   MSerror=deviance(m0)/df.residual(m0)))
                   })
    
    AinK <- ldply(AinK, NULL)
    AinK$M <- toupper(gsub(" ", "", AinK$M, fixed=TRUE))
    AinK$trt <- as.factor(as.numeric(as.character(AinK$trt)))
    str(AinK)
    
    #-----------------------------------------------------------------------------
    # combinando os data.frames para obter o gráfico e a tabela
    
    KinA$comb <- paste(KinA$trt, KinA$.id)
    AinK$comb <- paste(AinK$.id, AinK$trt)
    desd <- merge(KinA, AinK, by=c("comb","comb"))
    desd$let <- with(desd, paste(M.y, M.x, sep=""))
    desd <- desd[,c("trt.x","trt.y","means.x","let")]
    names(desd) <- c("K","A","means","let")
    str(desd)
    
    #-----------------------------------------------------------------------------
    # gráfico das médias com as letras sobre as barras
    
    #png("f005.png", w=500, h=250)
    barchart(means~K|A, data=desd, horiz=FALSE, layout=c(3,1),
             ylab="Rendimento de grãos (g/vaso)", ylim=c(NA,45),
             xlab="Dose de potássio (mg/dm³)",
             key=list(text=list("Níveis de umidade do solo (%)")),
             sub=list("*letras minúsculas comparam médias em um mesmo nível de umidade.", cex=0.8, font=3),
             strip=strip.custom(bg="gray90"), col=colors()[51],
             panel=function(x, y, subscripts, ...){
               panel.barchart(x, y, subscripts=subscripts, ...)
               panel.text(x, y, label=desd[subscripts,"let"], pos=3)
             })
    #dev.off()
    
    #-----------------------------------------------------------------------------
    # tabela de dupla entrada com as médias das combinações A e K
    
    require(reshape)
    desd$means.let <- with(desd, paste(means, let))
    cast(desd, K~A, value="means.let")
    
    #-----------------------------------------------------------------------------