Arquivo

Posts Tagged ‘lapply’

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

    Paralelização de Processos II

    Dando continuidade ao que foi apresentado no post sobre paralelização de processos e motivado pelo comentário do colega J. Franky e fundamentalmente com a ajuda de Benilton Carvalho (aos quais deixo meu agradecimento), nesse “fast-post” vou mostrar uma outra implementação em paralelo.

    O cenário é exatamente o mesmo do apresentado no post anterior (veja para entender). Agora, entretanto o processamento paralelo se dará em dois “níveis”, a saber: As nsmc simulações correm como um fork, ou seja, cópias do mesmo processo. Elas são implementadas como uma função e “correm” com a função mclapply(), o argumento mc.set.seed = TRUE garante que cada uma das nsmc terá uma semente geradora de números diferente.

    Dentro de cada simualação temos ainda dois processos independentes dentro de cada ciclo (as duas estratégias de amostragem). Para um caso em que cada uma dessas estratégias demore um tempo grande, ao invés de esperar que a primeira termine para então começar a segunda, usamos uma versão de thread, com as funções parallel() e collect() e paralelizamos as duas estratégias.

    O argumento mc.set.seed = TRUE tem o mesmo propósito que anteriormente, o argmento name = ‘foo’ é um grande facilitador… Vocês vão entender! Ao passar uma tarefa (a estratégia de amostragem) pela função parallel() é mesmo que dar um comando no terminal com um ‘nohup &’, ou seja, o computador executa a tarefa mas não “bloqueia” tal terminal para outro comando.

    Depois de paralelizados, os processos são coletados pela função collect() na estrutura de uma lista, temos como argumentos da função (que não estão sendo usados, veja a documentação) o wait, ele é muito útil quando queremos que se espera ‘X’ tempo pelo fim das tarefas, após esse tempo, o processo continua.

    A vantagem do uso do argumento name = ‘foo’ é que por default os nomes na estrutura da lista do collect recebem o número PID do processo, ao atribuir um nome fica muito mais fácil distribuir o que foi paralelizado.

    Espero que tenham gostado… Lembrando que este é um exemplo puramente didático, que provavelmente não se aplica ao uso desses procedimentos. Segue portando código da terceira implementação, agora com o uso de paralelização em dois níveis.

    #-----------------------------------------------------------------------------
    ## PARAMETROS DA SIMULAÇÃO -- usar os mesmo para as duas situações
    n <- 5000 # numero de individuos
    ciclos <- 30 # numero de ciclos
    nsmc <- 100 # numero de simulacoes
    
    ## SEGUNDO CASO -- paralelizada
    #-----------------------------------------------------------------------------
    ## IDEM (realiza um simulacao de 'nsmc')
    simula.ii <- function(x) { # inicio da funcao -- sem argumentos
      ## mesmos comentarios do caso acima
      resultados <- matrix(NA, ciclos, 4) 
      p0 <- sample(c(0, 1, 2), n, replace = TRUE, prob = c(.25, .5, .25))
      resultados[1, ] <- rep(c(mean(p0), var(p0)), times = 2)
      estrA <- sample(p0, n/5, replace = FALSE)
      estrB <- p0[seq(1, 5000, 5)]
      resultados[2, ] <- c(mean(estrA), var(estrA), mean(estrB), var(estrB))
      for(k in 3:ciclos) {
        ## pA e pB paralelizados em 'thread'
        pA <- parallel(tipoA(estrA), name = 'A', mc.set.seed = TRUE)
        pB <- parallel(tipoB(estrB), name = 'B', mc.set.seed = TRUE)
        estrG <- collect(list(pA, pB)) # coleta dos processos paralelizados
        estrA <- estrG$A; estrB <- estrG$B # redistribui os processos
        resultados[k, ] <- c(mean(estrA), var(estrA), mean(estrB), var(estrB))
      }
      return(resultados) # retorna uma simulacao
    }
    
    tempoC <- system.time({ # armazena o tempo de processamento
    saida3 <- mclapply(1:nsmc, # numero de nsmc
    # aplica a funcao -- faz a simulacao
                       simula.ii,
                       mc.preschedule = FALSE,
                       # 'expande' quantos processos forem possiveis
                       mc.set.seed = TRUE,
                       # uma semente para cada processo
                       mc.cores = getOption('cores'))
                       # usa quantos processadores forem possiveis
    })
    
    ## TEMPO
    tempoC[3]
    

    Até a próxima!

    Ajuste de muitos modelos de regressão linear

    Pessoas que estão migrando de outro aplicativo para o R normalmente vivem perguntando “o R faz isso?, o R faz aquilo?…”. A resposta é “Sim, o R faz!”. Circula entre usuários o seguinte ditado: “não se pergunta o que o R faz, e sim como ele faz”.

    Acontece que o R não é um software específico de ajuste de modelos de regressão, análise multivariada, análise de experimentos ou qualquer outra técnica. O R é geral. Os pacotes do R, por sua vez, são capazes de concentrar funções específicas. Um aplicativo que só faz ajuste de modelos de regressão tem espaço de sobra para incluir diversas opções, contem lista modelos para ajustar aos dados, enfim, mas tudo isso porque é destinado apenas para esse fim.

    Por outro lado, no R é possível desenvolver procedimentos para qualquer tipo de tarefa. Se você quer ajustar uma lista de modelos e depois escolher o de maior R², ótimo! No você também consegue isso. Veja no CMR abaixo.

    # gera dados
    da <- data.frame(x=runif(100), z=5*rpois(100, lambda=7), w=runif(100, 50, 100))
    da$y <- with(da, 12+0.1*x+0.05*z+0.34*w+0.2*sqrt(z)+0.1*x*w)+rnorm(100,0,0.1)
     
    # vetor com as fórmulas específicando diferentes modelos lineares
    form <- c(mod1=y~x, mod2=y~x+z, mod3=y~x+I(x^2), mod4=y~x+z+w)
     
    # ajuste dos modelos
    ajustes <- lapply(form, function(f){ m0 <- lm(f, data=da); m0 })
     
    lapply(ajustes, summary) # quadro geral de estimativas e qualidade
    lapply(ajustes, anova)   # quadro de anova sequencial
    lapply(ajustes, coef)    # vetor de estimativas
    sapply(ajustes, function(a){ summary(a)$r.squared})     # R²
    sapply(ajustes, function(a){ summary(a)$adj.r.squared}) # R² ajustado
    sapply(ajustes, function(a){ summary(a)$sigma})         # QMR
    sapply(ajustes, deviance)                               # SQR
    sapply(ajustes, df.residual)                            # GLR
    lapply(ajustes, function(a){ summary(a)$coeff})         # tabela de estimativas
    do.call(rbind, lapply(ajustes, function(a){ summary(a)$coeff})) # junta das tabelas
    sapply(ajustes, fitted)    # valores ajustados
    sapply(ajustes, residuals) # resíduos da análise
    sapply(ajustes, vcov)      # matriz de covariância das estimativas
    apply(sapply(ajustes, residuals), 2, shapiro.test) # normalidade dos resíduos
    

    Sempre faça a avaliação das pressuposições do modelo antes de aplicar inferências. A escolha de modelo apenas pelo valor do R² não é aconselhada.

    Realçador de código R para blogs

    Quem trabalha com programação reconhece que a organização do código facilita muito a releitura e checagem. O R, uma vez que é uma linguagem, possui editores que apresentam recursos para organização do código, como Emacs, RStudio, Tinn-R, Rkward entre diversos. A maioria dos editores apresenta um padrão de identação do código e realçador (highlighter) de texto. A necessidade de divulgação de código R, de forma organizada e atraente, fez com que os realçadores aparecessem nas páginas de blogs, wikis, etc. O wordpress.com possui um realçador de código R nativo. Para quem não tem wordpress.com, pode usar o Pretty R syntax highlighter que possui, além dos realces usuais, links para as páginas de documentação das funções. Abaixo está código R com o realçador do wordpress.com e com o pretty-R. Perceba que as funções (negrito azul) possuem links para página de ajuda.

    # importa dados
    soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                       header=TRUE, sep="\t", dec=",")
    str(soja)
     
    # ajusta um modelo e pede anova
    m1 <- aov(rengrao~bloco+agua*potassio, soja)
    anova(m1)
     
    # cria uma lista com as variáveis resposta
    respostas <- do.call(c, apply(soja[,4:7], 2, list))
    do.call(c, respostas)
     
    # faz o ajuste para todas as respostas
    ajustes <- lapply(respostas,
                      function(r){
                        m0 <- aov(r~bloco+agua*potassio, data=soja)
                        m0
                      })
     
    # pede todas as anovas
    lapply(ajustes, anova)
    
    # importa dados
    soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                       header=TRUE, sep="\t", dec=",")
    str(soja)
     
    # ajusta um modelo e pede anova
    m1 <- aov(rengrao~bloco+agua*potassio, soja)
    anova(m1)
     
    # cria uma lista com as variáveis resposta
    respostas <- do.call(c, apply(soja[,4:7], 2, list))
    do.call(c, respostas)
     
    # faz o ajuste para todas as respostas
    ajustes <- lapply(respostas,
                      function(r){
                        m0 <- aov(r~bloco+agua*potassio, data=soja)
                        m0
                      })
     
    # pede todas as anovas
    lapply(ajustes, anova)

    Created by Pretty R at inside-R.org

    Além de no R você fazer praticamente tudo, você pode apresentar de forma elegante. Até a próxima ridícula.