Introdução

A regressão linear simples é utilizada para analisar relações entre variáveis contínuas (KENNEY & KEEPING, 1962). No R, para especificar um modelo é preciso usar a notação de fórmulas. Por exemplo, para um modelo de regressão com uma variável resposta Y (ou dependente) e uma variável preditora X (ou independente), a notação de fórmula a ser usada no R é: Y ~ X. O argumento ~ (til) significa em relação à. . . ou modelado por. ..(ou seja, Y em relação à X, ou Y modelado por X). Caso haja duas variáveis preditoras no modelo, a fórmula será: Y~X1+X2.

# Carrega o pacote data.table para usar a função fread
library(data.table)
# Lê o arquivo CSV "Tectona.csv" e armazena em um data.frame chamado 'teca'
teca <- fread("Tectona.csv")
# Imprime o data.frame 'teca' para visualização
print(teca)
##     Arvore  DAP     H Volume
##  1:      1  6.6  8.90  0.017
##  2:      2  6.8  7.95  0.017
##  3:      3  7.0  6.87  0.014
##  4:      4  7.1 12.70  0.027
##  5:      5  7.5 11.74  0.028
##  6:      6  9.9 13.84  0.066
##  7:      7 10.8 15.25  0.049
##  8:      8 11.2 13.00  0.058
##  9:      9 11.3 14.60  0.066
## 10:     10 11.7 15.60  0.083
## 11:     11 12.6 17.20  0.100
## 12:     12 14.0 16.60  0.127
## 13:     13 14.3 17.87  0.158
## 14:     14 14.6 17.90  0.141
## 15:     15 14.9 16.50  0.132
## 16:     16 15.0 18.30  0.153
## 17:     17 15.3 17.40  0.162
## 18:     18 15.4 17.55  0.139
## 19:     19 16.2 16.99  0.192
## 20:     20 17.1 17.72  0.158
## 21:     21 17.4 16.20  0.196
## 22:     22 18.0 16.10  0.197
## 23:     23 18.1 15.34  0.250
## 24:     24 19.2 16.15  0.203
## 25:     25 19.8 14.38  0.262
## 26:     26 20.6 22.00  0.317
## 27:     27 22.5 20.60  0.377
## 28:     28 23.5 16.70  0.263
## 29:     29 23.8 22.80  0.395
## 30:     30 32.6 21.55  0.603
##     Arvore  DAP     H Volume

Gráficos de dispersão

# Anexa o data.frame 'teca' para acessar suas colunas diretamente
attach(teca)
# Plota gráfico de dispersão para DAP e Volume
plot(DAP, Volume, type = "p", xlab="DAP (cm)", ylab=expression(Volume~(m^3)))

# Plota gráfico de dispersão para Altura e Volume
plot(H, Volume, type = "p", xlab="Altura (m)", ylab=expression(Volume~(m^3)))

# Desanexa o data.frame 'teca'
detach(teca)

Ajustar os seguintes modelos volumétricos:  Berkhout: V=β0+β1 (DAP)V=β0+β1 (DAP)  Kopezky-Gerhardt: V=β0+β1 (DAP2)  Schumacher-Hall: lnV=β0+β1 (lnDAP)+β2 (lnH)

# Ajuste de modelos
# Ajusta o modelo de regressão linear para Berkhout
Berkhout <- lm(Volume ~ DAP, data=teca)
# Ajusta o modelo de regressão linear para Kopezky-Gerhardt
KGerhardt <- lm(Volume ~ I(DAP^2), data=teca)
# Ajusta o modelo de regressão linear para Schumacher-Hall
SHall <- lm(log(Volume) ~ log(DAP) + log(H), data=teca)

Para obter informações mais detalhadas do ajuste pode-se usar a função summary(). A partir dessa função, os seguintes valores podem ser observados: Erro-Padrão (𝑠𝑥̅), Estimativas dos Coeficientes de Regressão (Valores de 𝛽), Erro-Padrão de Estimativa (𝑠𝑦𝑥), Coeficiente de Determinação (𝑅²) e Coeficiente de Determinação Ajustado (𝑅²𝑎𝑗.). Além disso, pode-se visualizar o teste F de significância global e o do(s) teste(s) t-Student para significância do(s) coeficiente(s) da regressão.

# Resumo dos modelos
# Mostra o resumo estatístico do modelo de Berkhout
summary(Berkhout)
## 
## Call:
## lm(formula = Volume ~ DAP, data = teca)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.081315 -0.020684 -0.006809  0.026882  0.063030 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.160948   0.017105   -9.41 3.62e-10 ***
## DAP          0.021501   0.001053   20.43  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03375 on 28 degrees of freedom
## Multiple R-squared:  0.9371, Adjusted R-squared:  0.9349 
## F-statistic: 417.3 on 1 and 28 DF,  p-value: < 2.2e-16
# Mostra o resumo estatístico do modelo de Kopezky-Gerhardt
summary(KGerhardt)
## 
## Call:
## lm(formula = Volume ~ I(DAP^2), data = teca)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.077604 -0.015718 -0.004496  0.014418  0.064429 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.062e-03  9.014e-03   0.451    0.656    
## I(DAP^2)    6.094e-04  2.682e-05  22.718   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03053 on 28 degrees of freedom
## Multiple R-squared:  0.9485, Adjusted R-squared:  0.9467 
## F-statistic: 516.1 on 1 and 28 DF,  p-value: < 2.2e-16
# Mostra o resumo estatístico do modelo de Schumacher-Hall
summary(SHall)
## 
## Call:
## lm(formula = log(Volume) ~ log(DAP) + log(H), data = teca)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33962 -0.08448  0.00915  0.05905  0.30724 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -9.4394     0.2801 -33.697  < 2e-16 ***
## log(DAP)      1.8938     0.1186  15.968 2.81e-15 ***
## log(H)        0.8283     0.1765   4.693 6.95e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1432 on 27 degrees of freedom
## Multiple R-squared:  0.9803, Adjusted R-squared:  0.9788 
## F-statistic:   672 on 2 and 27 DF,  p-value: < 2.2e-16

A função anova() permite obter a tabela de ANOVA da Regressão:

# ANOVA da REGRESSÃO
# Obtém a tabela de ANOVA para o modelo de Berkhout
anova(Berkhout)
# Obtém a tabela de ANOVA para o modelo de Kopezky-Gerhardt
anova(KGerhardt)
# Obtém a tabela de ANOVA para o modelo de Schumacher-Hall
anova(SHall)

Pode-se visualizar a reta ajustada para o modelo de Berkhout (V=β0+β1 (DAP)V=β0+β1 (DAP)):

# Reta da Regressão ajustada
# Mostra a reta ajustada para o modelo de Berkhout
plot(teca$Volume~teca$DAP, main="Berkhout", xlab="DAP (cm)", ylab = expression(Volume~(m^3)))
abline(Berkhout, lty=2, col="red")

# Obtém os valores preditos pelo modelo de Berkhout
predict(Berkhout)
##             1             2             3             4             5 
## -0.0190447630 -0.0147446517 -0.0104445404 -0.0082944848  0.0003057378 
##             6             7             8             9            10 
##  0.0519070732  0.0712575740  0.0798577966  0.0820078522  0.0906080748 
##            11            12            13            14            15 
##  0.1099585756  0.1400593545  0.1465095215  0.1529596884  0.1594098553 
##            16            17            18            19            20 
##  0.1615599110  0.1680100779  0.1701601335  0.1873605787  0.2067110795 
##            21            22            23            24            25 
##  0.2131612464  0.2260615802  0.2282116359  0.2518622480  0.2647625818 
##            26            27            28            29            30 
##  0.2819630269  0.3228140842  0.3443146406  0.3507648075  0.5399697040
# Obtém os valores dos resíduos do modelo de Berkhout
residuals(Berkhout)
##            1            2            3            4            5            6 
##  0.036044763  0.031744652  0.024444540  0.035294485  0.027694262  0.014092927 
##            7            8            9           10           11           12 
## -0.022257574 -0.021857797 -0.016007852 -0.007608075 -0.009958576 -0.013059355 
##           13           14           15           16           17           18 
##  0.011490479 -0.011959688 -0.027409855 -0.008559911 -0.006010078 -0.031160134 
##           19           20           21           22           23           24 
##  0.004639421 -0.048711079 -0.017161246 -0.029061580  0.021788364 -0.048862248 
##           25           26           27           28           29           30 
## -0.002762582  0.035036973  0.054185916 -0.081314641  0.044235192  0.063030296

Análise de Resíduos Os resíduos da regressão linear devem ser avaliados quanto às pressuposições de normalidade, autocorrelação e homocedasticidade. O teste para normalidade dos resíduos pode ser realizado usado a função shapiro.test(). A autocorrelação pode ser detectada com uso da função durbinWatsonTest() e a heterocedasticidade dos resíduos com a função Breusch-Pagan().

# Análise de resíduos
# Teste de Shapiro-Wilk para normalidade dos resíduos
shapiro.test(Berkhout$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  Berkhout$residuals
## W = 0.97872, p-value = 0.7905
shapiro.test(KGerhardt$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  KGerhardt$residuals
## W = 0.96479, p-value = 0.4081
shapiro.test(SHall$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  SHall$residuals
## W = 0.98013, p-value = 0.829
# Teste Durbin-Watson para autocorrelação dos resíduos
library(car)
## Loading required package: carData
durbinWatsonTest(Berkhout)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.04712723      1.740416   0.352
##  Alternative hypothesis: rho != 0
durbinWatsonTest(KGerhardt)
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.2838086      2.469606   0.236
##  Alternative hypothesis: rho != 0
durbinWatsonTest(SHall)
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.1307352      2.182513    0.88
##  Alternative hypothesis: rho != 0
# Teste Breusch-Pagan para heterocedasticidade dos resíduos
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(Berkhout)
## 
##  studentized Breusch-Pagan test
## 
## data:  Berkhout
## BP = 9.7973, df = 1, p-value = 0.001748
bptest(KGerhardt)
## 
##  studentized Breusch-Pagan test
## 
## data:  KGerhardt
## BP = 12.264, df = 1, p-value = 0.0004617
bptest(SHall)
## 
##  studentized Breusch-Pagan test
## 
## data:  SHall
## BP = 1.9754, df = 2, p-value = 0.3724
#Plotar gráficos
plot(Berkhout)

plot(KGerhardt)

plot(SHall)

A inexistência de multicolinearidade é outra importante pressuposição que deve ser analisada em modelos de regressão linear múltipla. O grau de confundimento entre variáveis preditoras pode ser verificado por meio da estatística Variance Inflation Factor – (VIF), que detecta problemas de multicolinearidade. É recorrente o uso de VIF > 10 como critério para afirmar a existência de multicolinearidade no modelo de regressão linear múltipla. Para obter o valor desse índice, pode-se usar a função vif() do pacote faraway (FARAWAY, 2016).

# Multicolinearidade
# Verifica a multicolinearidade através do VIF
library(faraway)
## 
## Attaching package: 'faraway'
## The following objects are masked from 'package:car':
## 
##     logit, vif
vif(SHall)
## log(DAP)   log(H) 
## 3.339837 3.339837