Curso de R no IPÊ
Script do dia 5.2: regressão múltipla
Marcos V. C. Vital e participantes da turma de 2014
Para mais material, acesse http://marcosvital.wordpress.com/
Este documento é a 2a parte do quinto (e último!) de uma série que nasceu durante as atividades do “Curso de introdução ao uso do software R e suas aplicações estatísticas nas ciências biológicas”, realizado no IPÊ em dezembro de 2014.
Nesta parte do guia nós vamos realizar uma regressão múltipla. No caminho, vamos ver um jeito bacana de visualizar a correlação entre variáveis.
Como sempre, todo o documento foi criado utilizando o software R, o programa RStudio e a linguagem Markdown. Para saber mais, acesse:
Antes de começar, carregue os dados e os pacotes necessários:
setwd("C:/R/Curso IPÊ")
dados<-read.table("DensPraga.txt", h=T)
attach(dados)
summary(dados)## Parcelas Temperatura Altitude Pluviosidade
## Min. : 1 Min. :22.70 Min. :160.9 Min. : 519.3
## 1st Qu.:14 1st Qu.:23.81 1st Qu.:364.5 1st Qu.: 810.7
## Median :27 Median :24.86 Median :559.9 Median :1159.9
## Mean :27 Mean :25.00 Mean :519.5 Mean :1174.1
## 3rd Qu.:40 3rd Qu.:26.23 3rd Qu.:649.6 3rd Qu.:1465.2
## Max. :53 Max. :27.63 Max. :771.3 Max. :1985.1
## Mata Diversidade Dens_Praga
## Min. :0.0039 Min. :0.4434 Min. : 5.703
## 1st Qu.:0.2631 1st Qu.:1.2776 1st Qu.: 88.777
## Median :0.3841 Median :1.7430 Median :158.635
## Mean :0.3731 Mean :1.6921 Mean :144.019
## 3rd Qu.:0.5134 3rd Qu.:2.1638 3rd Qu.:194.387
## Max. :0.6881 Max. :2.8500 Max. :308.207
Explorando a correlação entre as variáveis
Uma coisa fundamental em uma regressão múltipla é utilizar variáveis explicativas que não estejam correlaciomnadas entre si. Podemos avaliar se isto acontece com nossos dados calculando um coeficiente de correlação entre as variáveis de nosso interesse.
Calculando uma matriz de correlação de spearman:
cor(dados[2:6])## Temperatura Altitude Pluviosidade Mata Diversidade
## Temperatura 1.00000000 -0.96273540 0.07872156 -0.03100336 -0.02726357
## Altitude -0.96273540 1.00000000 -0.07251465 0.07052293 0.03335541
## Pluviosidade 0.07872156 -0.07251465 1.00000000 0.09636866 0.36479850
## Mata -0.03100336 0.07052293 0.09636866 1.00000000 0.63927535
## Diversidade -0.02726357 0.03335541 0.36479850 0.63927535 1.00000000Simples e prático, basta interpretar os valores.
Mas como sou chato, quero algo que me permita olhar para as correlações ao mesmo tempo que vejo em gráficos como os dados se relacionam.
O nosso ponto de partida é a função pairs:
pairs(dados[2:6])
Que já é bem legal, mas pode ser melhor aproveitada. Que tal colocar coisas úteis na diagonal principal e no painel superior?
Nós vamos copiar e colar um trecho do exemplo do help da função, fazer uma modificação bem pequena, e executar. Se não quiser mexer com isso diretamente, basta executar o código abaixo:
panel.hist <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = 1)
}Beleza, agora podemos re-escrever nossa função, e criar um gráfico bem útil para se explorar os padrões de correlação entre as variáveis (com um histograma de cada uma de brinde no meio!).
Vamos lá:
pairs(dados[2:6],
diag.panel=panel.hist, upper.panel=panel.cor)
Vejamos:
- o argumento
diag.paneldefine o que vai aparecer na diagonal principal da matriz de gráficos; - o argumento
upper.paneldefine o que aparece no painel superior; panel.histepanel.corforam definidos naquele blocão de código logo acima.
Ficou bonito, né? Agora vamos seguir e realizar a nossa regressão. O importante aqui é que agora sabemos que algumas das variáveis estão correlacionadas, e podemos escolher quais variáveis entram e quais não entram no nosso modelo estatístico.
Fazendo a regressão múltipla
É facinho! Se você já fez uma regressão múltipla, verá que não tem nada demais aqui:
modelo<-lm(Dens_Praga~Altitude+Pluviosidade+Mata)Vamos dar uma olhada nos pressupostos?
par(mfrow=c(1,2))
#Homocedasticidade:
plot(modelo$fitted, modelo$residuals,
abline(h=0,lty=2))
#Normalidade:
qqnorm(modelo$residuals)
qqline(modelo$residuals, lty=2)
Não está fantástico, mas considero como aceitável.
Então vamos ver o resultado do modelo:
summary(modelo)##
## Call:
## lm(formula = Dens_Praga ~ Altitude + Pluviosidade + Mata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.83 -13.73 -0.05 11.97 59.47
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.378e+02 1.551e+01 21.779 < 2e-16 ***
## Altitude -1.197e-01 1.989e-02 -6.016 2.21e-07 ***
## Pluviosidade 7.005e-03 7.703e-03 0.909 0.368
## Mata -3.748e+02 1.898e+01 -19.744 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.57 on 49 degrees of freedom
## Multiple R-squared: 0.9013, Adjusted R-squared: 0.8953
## F-statistic: 149.2 on 3 and 49 DF, p-value: < 2.2e-16Beleza! A interpretação não tem mistério, é similar a de uma regressão simples.
Se quisermos, podemos fazer os gráficos para enxergar melhor o efeito das variáveis que fora estatisticamente significantes:
par(mfrow=c(1,2))
plot(Dens_Praga~Altitude)
plot(Dens_Praga~Mata)
Beleza, pessoal, estamos quase terminando, e o próximo é o último da nossa compilação do curso.