Arquivo

Posts Tagged ‘do.call’

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.