10.6 Regressão Linear

A classe de Modelos Lineares atende a uma ampla gama de problemas aplicados, apresentada por (Neter et al. 2005), podendo incluir modelos lineares de séries temporais (Enders 2014). Há generalizações tais como os Modelos Lineares Generalizados (MLG/GLM) (McCullagh and Nelder 1989), Modelos Aditivos Generalizados (MAG/GAM) (Hastie and Tibshirani 1986), dentre outras. (J. H. Friedman and Stuetzle 1981) sugerem a Regressão com Perseguição da Projeção, procedimento que modela a superfície de regressão como uma soma de funções de suavização de combinações lineares das variáveis preditoras de forma iterativa. (O’Hagan and Forster 1994, 244–76) e (Paulino, Turkman, and Murteira 2003, 194–208) fornecem uma abordagem bayesiana para a regressão.

Quando trabalha-se com dados multivariados (i.e., quase sempre), recomenda-se a exploração criteoriosa dos dados. Não raro ocorrem valores faltantes (missings), rotulados por NA na linguagem R conforme Seção NA no Estatística Básica do mesmo autor. Tais valores faltantes podem dificultar os procedimentos de modelagem, e sua ocorrência se soma a outras questões como a multicolinearidade, discutida na Seção 10.6.3.

10.6.1 Clássico

10.6.1.1 Modelo

O modelo de Regressão Linear Múltipla (RLM/MLR) populacional/universal é construído, na abordagem clássica, com todos as \(N\)-uplas da população ou universo, i.e., \((x_{1i},x_{2i},\ldots,x_{pi},y_i), \; i \in \{1,\ldots,N\}\). Pode ser descrito pela Eq. (10.24).

\[\begin{equation} y_i = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_p x_{pi} + \varepsilon_i, \tag{10.24} \end{equation}\] onde \(\varepsilon_i \sim \mathcal{N}(0,\sigma_{\varepsilon})\).

Por ter dimensionalidade \(p\) usualmente utiliza-se notação matricial na forma

\[\begin{equation} Y = X \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \tag{10.25} \end{equation}\]

onde \(Y = \begin{bmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{bmatrix}\), \(X = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{bmatrix}\), \(\boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{bmatrix}\) e \(\boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix}\).

Estimativa dos parâmetros

Para a obtenção das estimativas dos parâmetros utiliza-se a Eq. (10.26). Para ajustar um modelo RPO, basta eliminar a coluna unitária de \(X\). \[\begin{equation} \boldsymbol{\hat{\beta}} = (X'X)^{-1} X'Y \tag{10.26} \end{equation}\]

Exemplo 10.8 Em http://archive.ics.uci.edu/ml/datasets/Energy+efficiency está disponível uma análise de energia feita por (Tsanas and Xifara 2012) usando 12 formas diferentes de construção simuladas no Ecotect. Os edifícios diferem em relação à área envidraçada, à distribuição da área envidraçada e à orientação, entre outros parâmetros. Foram simuladas várias configurações como funções das características acima mencionadas para obter 768 formas de construção. O conjunto de dados detalhado a seguir compreende 768 amostras e 8 características (X1 a X8), com o objetivo de prever duas respostas reais (Y1 e Y2).

X1: Compactação Relativa
X2: Superfície
X3: Área da parede
X4: Área do telhado
X5: Altura total
X6: Orientação
X7: Área de Envidraçamento
X8: Distribuição da Área de Envidraçamento
Y1: Carga de aquecimento
Y2: Carga de resfriamento

# arquivo
url1 = 'http://archive.ics.uci.edu/ml/machine-learning-databases/00242/ENB2012_data.xlsx'
download.file(url1, 'temp.xlsx', mode = 'wb')
energy = readxl::read_excel('temp.xlsx')

# exploratória
skimr::skim(energy)
Tabela 10.1: Data summary
Name energy
Number of rows 768
Number of columns 10
_______________________
Column type frequency:
numeric 10
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
X1 0 1 0.76 0.11 0.62 0.68 0.75 0.83 0.98 ▇▆▃▃▂
X2 0 1 671.71 88.09 514.50 606.38 673.75 741.12 808.50 ▅▅▇▅▇
X3 0 1 318.50 43.63 245.00 294.00 318.50 343.00 416.50 ▃▅▇▂▂
X4 0 1 176.60 45.17 110.25 140.88 183.75 220.50 220.50 ▃▃▁▁▇
X5 0 1 5.25 1.75 3.50 3.50 5.25 7.00 7.00 ▇▁▁▁▇
X6 0 1 3.50 1.12 2.00 2.75 3.50 4.25 5.00 ▇▇▁▇▇
X7 0 1 0.23 0.13 0.00 0.10 0.25 0.40 0.40 ▂▇▁▇▇
X8 0 1 2.81 1.55 0.00 1.75 3.00 4.00 5.00 ▇▆▆▆▆
Y1 0 1 22.31 10.09 6.01 12.99 18.95 31.67 43.10 ▇▇▃▆▅
Y2 0 1 24.59 9.51 10.90 15.62 22.08 33.13 48.03 ▇▂▅▃▁
# matriz de dispersão
pairs(energy, panel = panel.smooth, main = "Energy")

# histograma de Y1
hist(energy$Y1, 30)

# teste de normalidade
shapiro.test(energy$Y1)
## 
##  Shapiro-Wilk normality test
## 
## data:  energy$Y1
## W = 0.9121, p-value < 2.2e-16
# modelo saturado
fit0 <- lm(Y1 ~ . -Y2, data = energy)
summary(fit0)
## 
## Call:
## lm(formula = Y1 ~ . - Y2, data = energy)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.8965 -1.3196 -0.0252  1.3532  7.7052 
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  84.013418  19.033613   4.414 1.16e-05 ***
## X1          -64.773432  10.289448  -6.295 5.19e-10 ***
## X2           -0.087289   0.017075  -5.112 4.04e-07 ***
## X3            0.060813   0.006648   9.148  < 2e-16 ***
## X4                  NA         NA      NA       NA    
## X5            4.169954   0.337990  12.338  < 2e-16 ***
## X6           -0.023330   0.094705  -0.246  0.80548    
## X7           19.932736   0.813986  24.488  < 2e-16 ***
## X8            0.203777   0.069918   2.915  0.00367 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.934 on 760 degrees of freedom
## Multiple R-squared:  0.9162, Adjusted R-squared:  0.9154 
## F-statistic:  1187 on 7 and 760 DF,  p-value: < 2.2e-16

No modelo atribuído a fit0, é possível notar que não é possível calcular o coeficiente da variável X4 devido a singularidades. Para mais detalhes, veja esta discussão do Stackoverflow.

Stepwise

O método de regressão stepwise foi proposto por (Efroymson 1960) com o intuito de selecionar variáveis em regressões múltiplas conforme Eq. (10.24). De acordo com o autor11, no procedimento de stepwise os resultados intermediários dos métodos para a obtenção das estimativas de regressões múltiplas “são usados para fornecer informações estatísticas valiosas em cada etapa do cálculo”. O autor também ressalta12 que o método “é baseado nos fatos que (a) uma variável pode ser indicada como significativa em qualquer estágio inicial e, assim, entrar na equação e (b) após várias outras variáveis serem adicionadas à equação de regressão, o variável inicial pode ser indicada como insignificante. A variável insignificante será removida da equação de regressão antes de adicionar uma nova variável. Portanto, apenas variáveis significativas são incluídas na regressão final”.

Este método busca automaticamente o melhor conjunto de variáveis de maneira a minimizar alguma medida. Na implementação de R documentada em stats::step(), o modelo selecionado será aquele com o menor valor de AIC (Akaike 1974) (k = 2) ou BIC (Schwarz 1978) (k = log(n)). Note que ‘seleção de modelos’ e ‘seleção de variáveis’ são sinônimos.

Exemplo 10.9 A partir do modelo saturado do Exemplo 10.8 pode-se ajustar uma regressão stepwise.

# filtrando as variáveis com stepwise
fit1 <- step(fit0, trace = 0)
summary(fit1)
## 
## Call:
## lm(formula = Y1 ~ X1 + X2 + X3 + X5 + X7 + X8, data = energy)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.9315 -1.3189 -0.0262  1.3587  7.7169 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  83.931762  19.018978   4.413 1.17e-05 ***
## X1          -64.773432  10.283096  -6.299 5.06e-10 ***
## X2           -0.087289   0.017065  -5.115 3.97e-07 ***
## X3            0.060813   0.006644   9.153  < 2e-16 ***
## X5            4.169954   0.337781  12.345  < 2e-16 ***
## X7           19.932736   0.813484  24.503  < 2e-16 ***
## X8            0.203777   0.069875   2.916  0.00365 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.933 on 761 degrees of freedom
## Multiple R-squared:  0.9162, Adjusted R-squared:  0.9155 
## F-statistic:  1387 on 6 and 761 DF,  p-value: < 2.2e-16

Pelos resultados obtidos pode-se verificar que o modelo ajustado em fit1 possui todas as variávies significativas para um \(\alpha\) inferior a 0.01, estatística F de 1387 para 6 e 761 graus de liberdade, levando a um p-value geral do modelo menor que \(2.2 \times 10^{-16}\), o que indica boa aderência aos dados. O valor do Multiple R-squared é o coeficiente de determinação de valor 0.9162, indicando que o modelo explica em torno de 92% da variação de Y1.

Com este diagnóstico pode-se considerar um modelo aceitável, que possui coeficientes de X1 e X2 negativos, indicando que um aumento destas variáveis (respectivamente compactação relativa e superfície) deve reduzir a carga de aquecimento. Mais especificamente, um aumento de uma unidade na compactação relativa (X1) gera uma redução esperada de 64.77 unidades na carga de aquecimento, mantidas constantes as demais variáveis. As variáveis X3, X5, X7 e X8 possuem coeficientes positivos, levando a um impacto esperado positivo em Y1. Como exemplo, para cada aumento de uma unidade na altura total (X5) espera-se um aumento aproximando de 4.17 unidades na carga de aquecimento. As outras variáveis possuem interpretação análoga, devendo-se sempre observar o sinal dos coeficientes.

10.6.1.2 LASSO

O modelo LASSO (Least Absolute Shrinkage and Selection Operator) também é definido pela Eq. (10.24), porém com coeficientes obtidos ao \[\begin{equation} \text{minimizar} \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_{i1} - \beta_2 X_{i2} - \cdots - \beta_p X_{ip})^{2} \\ \text{ sujeito a} \sum_{j=1}^p |\beta_j| \le t \tag{4} \end{equation}\] onde \(t\) é o parâmetro livre ou de regularização.

10.6.1.3 Diagnóstico inferencial

Teste para \(\boldsymbol{\beta}\)

Exemplo 10.10 Considere os dados do Exemplo 10.9. Pode-se gerar um relatório do modelo desejado com base::summary(). Um diagnóstico visual pode ser obtido via base::plot(), onde which=1:4 indica que estão sendo solicitados os 4 primeiros gráficos disponíves de um total de 6.

summary(fit1)
## 
## Call:
## lm(formula = Y1 ~ X1 + X2 + X3 + X5 + X7 + X8, data = energy)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.9315 -1.3189 -0.0262  1.3587  7.7169 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  83.931762  19.018978   4.413 1.17e-05 ***
## X1          -64.773432  10.283096  -6.299 5.06e-10 ***
## X2           -0.087289   0.017065  -5.115 3.97e-07 ***
## X3            0.060813   0.006644   9.153  < 2e-16 ***
## X5            4.169954   0.337781  12.345  < 2e-16 ***
## X7           19.932736   0.813484  24.503  < 2e-16 ***
## X8            0.203777   0.069875   2.916  0.00365 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.933 on 761 degrees of freedom
## Multiple R-squared:  0.9162, Adjusted R-squared:  0.9155 
## F-statistic:  1387 on 6 and 761 DF,  p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(fit1, which=1:4)

10.6.1.4 Diagnóstico preditivo

Exemplo 10.11 Considere novamente os dados do Exemplo ??.

# arquivo
url1 = 'http://archive.ics.uci.edu/ml/machine-learning-databases/00242/ENB2012_data.xlsx'
download.file(url1, 'temp.xlsx', mode = 'wb')
energy = readxl::read_excel('temp.xlsx')

# treino-teste
set.seed(999)
itrain <- sort(sample(1:nrow(energy), ceiling(0.7*nrow(energy))))
train <- energy[itrain,]
test <- energy[-itrain,]

# modelo
fit0 <- lm(Y1 ~ . -Y2, data = train)
fit1 <- step(fit0, trace = 0)

# predição
pred <- predict(fit1, test)

# resíduos de predição
r <- pred-test$Y1 # predito - realizado
(mse <- mean(r^2)) # https://en.wikipedia.org/wiki/Mean_squared_error
## [1] 9.330395
(rmse <- sqrt(mse)) # https://en.wikipedia.org/wiki/Root_mean_square_deviation
## [1] 3.05457
(mae <- mean(abs(r))) # https://en.wikipedia.org/wiki/Mean_absolute_error
## [1] 2.219659
(mape <- mean(abs(r/test$Y1))) # https://en.wikipedia.org/wiki/Mean_absolute_percentage_error
## [1] 0.1052415
hist(r, 30)

# 100 repetições
M <- 100
residuos <- data.frame(modelo=1:M, mse=NA, rmse=NA, mae=NA, mape=NA)
for(i in 1:M){
  # treino-teste
  set.seed(i)
  itrain <- sort(sample(1:nrow(energy), ceiling(0.7*nrow(energy))))
  train <- energy[itrain,]
  test <- energy[-itrain,]

  # modelo
  fit0 <- lm(Y1 ~ . -Y2, data = train)
  fit1 <- step(fit0, trace=0)

  # predição
  pred <- predict(fit1, test)
  
  # resíduos de predição
  r <- pred-test$Y1
  residuos[i,2] <- mean(e^2)
  residuos[i,3] <- sqrt(residuos[i,2])
  residuos[i,4] <- mean(abs(r))
  residuos[i,5] <- mean(abs(r/test$Y1))
}
par(mfrow=c(2,2))
hist(residuos[,2], main = 'MSE', 30)
hist(residuos[,3], main = 'RMSE', 30)
hist(residuos[,4], main = 'MAE', 30)
hist(residuos[,5], main = 'MAPE', 30)

skimr::skim(residuos)
Tabela 10.2: Data summary
Name residuos
Number of rows 100
Number of columns 5
_______________________
Column type frequency:
numeric 5
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
modelo 0 1 50.50 29.01 1.00 25.75 50.50 75.25 100.00 ▇▇▇▇▇
mse 0 1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 ▁▁▇▁▁
rmse 0 1 0.01 0.00 0.01 0.01 0.01 0.01 0.01 ▁▁▇▁▁
mae 0 1 2.08 0.10 1.84 2.02 2.07 2.15 2.42 ▂▇▆▃▁
mape 0 1 0.10 0.01 0.09 0.09 0.10 0.10 0.12 ▃▇▅▁▁

Exercício 10.1 Acesse o endereço http://archive.ics.uci.edu/ml/datasets.php e escolha um banco de dados com mais de 100 atributos (variáveis/colunas). Faça análises considerando modelos lineares e, se pertinente, utilizando componentes principais. \(\\\)

Exercício 10.2 Refaça o Exemplo 10.8 utilizando a variável Y2 como resposta, ajustando o modelo saturado, filtrando as variáveis com a função step e interpretando os resultados.

Exemplo 10.12 (Belsley, Kuh, and Welsch 2004, 39) analisam a hipótese desenvolvida por (Modigliani and Papademos 1975), de que a a razão de poupança ao longo da vida (poupança pessoal agregada dividido pelo rendimento disponível) é explicada pelo rendimento disponível per capita, a taxa de variação percentual do rendimento disponível per capita e duas variáveis demográficas: a proporção da população com menos de 15 anos e a proporção da população com mais de 75 anos. Os dados são calculados ao longo da década de 1960-1970 para remover o ciclo de negócios ou outras flutuações de curto prazo.

par(mfrow=c(3,1))
pairs(LifeCycleSavings, panel = panel.smooth, main = "LifeCycleSavings")

hist(LifeCycleSavings$sr, 30)
qqnorm(LifeCycleSavings$sr)
qqline(LifeCycleSavings$sr, col = 'red')
shapiro.test(LifeCycleSavings$sr)
## 
##  Shapiro-Wilk normality test
## 
## data:  LifeCycleSavings$sr
## W = 0.98074, p-value = 0.5836
fm0 <- lm(sr ~ ., data = LifeCycleSavings)
fm1 <- step(fm0, trace = 0)
summary(fm1)
## 
## Call:
## lm(formula = sr ~ pop15 + pop75 + ddpi, data = LifeCycleSavings)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2539 -2.6159 -0.3913  2.3344  9.7070 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.1247     7.1838   3.915 0.000297 ***
## pop15        -0.4518     0.1409  -3.206 0.002452 ** 
## pop75        -1.8354     0.9984  -1.838 0.072473 .  
## ddpi          0.4278     0.1879   2.277 0.027478 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.767 on 46 degrees of freedom
## Multiple R-squared:  0.3365, Adjusted R-squared:  0.2933 
## F-statistic: 7.778 on 3 and 46 DF,  p-value: 0.0002646
mctest::mctest(fm1)
## 
## Call:
## omcdiag(mod = mod, Inter = TRUE, detr = detr, red = red, conf = conf, 
##     theil = theil, cn = cn)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.1739         0
## Farrar Chi-Square:        82.4971         1
## Red Indicator:             0.5254         1
## Sum of Lambda Inverse:    12.4857         0
## Theil's Method:            0.9827         1
## Condition Number:         31.7169         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test
sort(car::vif(fm1), decreasing = TRUE) 
##    pop15    pop75     ddpi 
## 5.745478 5.736014 1.004186

Exercício 10.3 Considere o Exemplo 10.12.
a. A partir da saída de car::vif, você recomenda retirar alguma variável do modelo? Em caso afirmativo, execute sua recomendação.
b. Teste jtools::summ(fm1).

Exemplo 10.13 Considere o banco de dados disponível em https://archive.ics.uci.edu/dataset/162/forest+fires baseado em (Cortez and Morais 2007).

  1. X - x-axis spatial coordinate within the Montesinho park map: 1 to 9
  2. Y - y-axis spatial coordinate within the Montesinho park map: 2 to 9
  3. month - month of the year: ‘jan’ to ‘dec’
  4. day - day of the week: ‘mon’ to ‘sun’
  5. FFMC - FFMC index from the FWI system: 18.7 to 96.20
  6. DMC - DMC index from the FWI system: 1.1 to 291.3
  7. DC - DC index from the FWI system: 7.9 to 860.6
  8. ISI - ISI index from the FWI system: 0.0 to 56.10
  9. temp - temperature in Celsius degrees: 2.2 to 33.30
  10. RH - relative humidity in %: 15.0 to 100
  11. wind - wind speed in km/h: 0.40 to 9.40
  12. rain - outside rain in mm/m2 : 0.0 to 6.4
  13. area - the burned area of the forest (in ha): 0.00 to 1090.84 (this output variable is very skewed towards 0.0, thus it may make sense to model with the logarithm transform).
dat <- read.table('https://filipezabala.com/data/forestfires.csv',
                 sep = ',', head = TRUE)
head(dat)
##   X Y month day FFMC  DMC    DC  ISI temp RH wind rain area
## 1 7 5   mar fri 86.2 26.2  94.3  5.1  8.2 51  6.7  0.0    0
## 2 7 4   oct tue 90.6 35.4 669.1  6.7 18.0 33  0.9  0.0    0
## 3 7 4   oct sat 90.6 43.7 686.9  6.7 14.6 33  1.3  0.0    0
## 4 8 6   mar fri 91.7 33.3  77.5  9.0  8.3 97  4.0  0.2    0
## 5 8 6   mar sun 89.3 51.3 102.2  9.6 11.4 99  1.8  0.0    0
## 6 8 6   aug sun 92.3 85.3 488.0 14.7 22.2 29  5.4  0.0    0
table(dat$X)
## 
##  1  2  3  4  5  6  7  8  9 
## 48 73 55 91 30 86 60 61 13
dat$X <- factor(dat$X, ordered = T)
dat$Y <- factor(dat$Y, ordered = T)
str(dat)
## 'data.frame':    517 obs. of  13 variables:
##  $ X    : Ord.factor w/ 9 levels "1"<"2"<"3"<"4"<..: 7 7 7 8 8 8 8 8 8 7 ...
##  $ Y    : Ord.factor w/ 7 levels "2"<"3"<"4"<"5"<..: 4 3 3 5 5 5 5 5 5 4 ...
##  $ month: chr  "mar" "oct" "oct" "mar" ...
##  $ day  : chr  "fri" "tue" "sat" "fri" ...
##  $ FFMC : num  86.2 90.6 90.6 91.7 89.3 92.3 92.3 91.5 91 92.5 ...
##  $ DMC  : num  26.2 35.4 43.7 33.3 51.3 ...
##  $ DC   : num  94.3 669.1 686.9 77.5 102.2 ...
##  $ ISI  : num  5.1 6.7 6.7 9 9.6 14.7 8.5 10.7 7 7.1 ...
##  $ temp : num  8.2 18 14.6 8.3 11.4 22.2 24.1 8 13.1 22.8 ...
##  $ RH   : int  51 33 33 97 99 29 27 86 63 40 ...
##  $ wind : num  6.7 0.9 1.3 4 1.8 5.4 3.1 2.2 5.4 4 ...
##  $ rain : num  0 0 0 0.2 0 0 0 0 0 0 ...
##  $ area : num  0 0 0 0 0 0 0 0 0 0 ...
# descritivas
hist(dat$area, 30)

larea <- log(dat$area+1)
hist(larea, 30)

# modelos
fit0 <- lm(area ~ ., data = dat)
fit1 <- step(fit0, trace = 0)
summary(fit1)
## 
## Call:
## lm(formula = area ~ temp, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -27.34  -14.68  -10.39   -3.42 1071.33 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -7.4138     9.4996  -0.780   0.4355  
## temp          1.0726     0.4808   2.231   0.0261 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 63.41 on 515 degrees of freedom
## Multiple R-squared:  0.009573,   Adjusted R-squared:  0.00765 
## F-statistic: 4.978 on 1 and 515 DF,  p-value: 0.0261
par(mfrow=c(2,2))
plot(fit1, which=1:4)

# modelando log area
fit2 <- lm(larea ~ ., data = dat[,-13]) # tirando a coluna 13, area
fit3 <- step(fit2, trace = 0)
summary(fit3)
## 
## Call:
## lm(formula = larea ~ X + Y + month + DMC + DC + temp + wind, 
##     data = dat[, -13])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1879 -0.9764 -0.4079  0.7768  5.3173 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.780632   0.575046   1.358 0.175247    
## X.L          0.984402   0.345491   2.849 0.004567 ** 
## X.Q          1.229688   0.339002   3.627 0.000317 ***
## X.C          0.453981   0.291630   1.557 0.120191    
## X^4          0.406301   0.259677   1.565 0.118317    
## X^5          0.377971   0.224246   1.686 0.092529 .  
## X^6         -0.116953   0.219999  -0.532 0.595241    
## X^7          0.233594   0.174057   1.342 0.180203    
## X^8         -0.564175   0.218151  -2.586 0.009994 ** 
## Y.L          0.188704   0.713436   0.264 0.791507    
## Y.Q         -1.551066   0.453561  -3.420 0.000679 ***
## Y.C         -2.270578   0.674401  -3.367 0.000821 ***
## Y^4         -2.771889   0.808894  -3.427 0.000663 ***
## Y^5         -1.888967   0.609695  -3.098 0.002059 ** 
## Y^6         -0.753255   0.304775  -2.472 0.013795 *  
## monthaug     0.643724   0.788255   0.817 0.414531    
## monthdec     2.213019   0.773107   2.863 0.004384 ** 
## monthfeb     0.365346   0.551503   0.662 0.507993    
## monthjan    -0.682928   1.076240  -0.635 0.526020    
## monthjul     0.436781   0.683625   0.639 0.523176    
## monthjun    -0.394812   0.630963  -0.626 0.531786    
## monthmar    -0.230436   0.493123  -0.467 0.640495    
## monthmay     0.736262   1.056583   0.697 0.486240    
## monthnov    -0.967716   1.428329  -0.678 0.498400    
## monthoct     1.315924   0.957744   1.374 0.170079    
## monthsep     1.360462   0.887515   1.533 0.125952    
## DMC          0.004546   0.001774   2.563 0.010675 *  
## DC          -0.002498   0.001231  -2.029 0.042954 *  
## temp         0.023922   0.014498   1.650 0.099583 .  
## wind         0.058852   0.036372   1.618 0.106294    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.345 on 487 degrees of freedom
## Multiple R-squared:  0.1268, Adjusted R-squared:  0.07477 
## F-statistic: 2.438 on 29 and 487 DF,  p-value: 6e-05
# retirar coordenadas X e Y
fit4 <- lm(larea ~ ., data = dat[,-c(1:3,13)]) # tirando a colunas 1,2, 13
fit5 <- step(fit4, trace = 0)
summary(fit5)
## 
## Call:
## lm(formula = larea ~ DMC + RH + wind, data = dat[, -c(1:3, 13)])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5093 -1.1071 -0.6537  0.9216  5.7712 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.9129363  0.2433065   3.752 0.000195 ***
## DMC          0.0017552  0.0009657   1.817 0.069731 .  
## RH          -0.0055830  0.0037785  -1.478 0.140138    
## wind         0.0624134  0.0345109   1.809 0.071112 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.392 on 513 degrees of freedom
## Multiple R-squared:  0.01425,    Adjusted R-squared:  0.008485 
## F-statistic: 2.472 on 3 and 513 DF,  p-value: 0.06101

10.6.2 Bayesiano

Uma prévia pode ser encontrada https://alpopkes.com/posts/machine_learning/bayesian_linear_regression/

Exemplo 10.14 Ajustando os modelos lineares clássico e bayesiano para o banco de dados datasets::cars.

ggplot2::qplot(speed, dist, data = cars)

# Regressão linear clássica por MQO
fit_lm <- lm(dist ~ speed, data = cars)
summary(fit_lm)
## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12
par(mfrow=c(2,2))
plot(fit_lm, which = 1:4)

# Regressão linear bayesiana via `rstanarm::stan_glm`
library(rstan)
library(rstanarm)
options(mc.cores = parallel::detectCores())

fit_stan <- rstanarm::stan_glm(dist ~ speed, data = cars, family = gaussian)
rstan::stan_trace(fit_stan, pars = c('(Intercept)', 'speed', 'sigma'))

summary(fit_stan)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      dist ~ speed
##  algorithm:    sampling
##  sample:       4000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 50
##  predictors:   2
## 
## Estimates:
##               mean   sd    10%   50%   90%
## (Intercept) -17.7    6.9 -26.5 -17.7  -8.9
## speed         3.9    0.4   3.4   3.9   4.5
## sigma        15.7    1.6  13.7  15.6  17.9
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 42.8    3.1 38.8  42.9  46.7 
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.1  1.0  3609 
## speed         0.0  1.0  3466 
## sigma         0.0  1.0  3401 
## mean_PPD      0.0  1.0  3986 
## log-posterior 0.0  1.0  1438 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
rstanarm::pp_check(fit_stan)

rstanarm::posterior_vs_prior(fit_stan, group_by_parameter = TRUE,
                             pars = c('(Intercept)', 'speed', 'sigma'))

Exercício 10.4 Veja as seguintes documentações:
a. ?rstanarm::print.stanreg.
b. ?rstanarm::prior_summary.
c. ?rstanarm::summary.stanreg.
d. ?rstanarm::pp_check.
e. ?rstanarm::launch_shinystan.
f. ?rstanarm::posterior_vs_prior.

10.6.3 Multicolinearidade

Segundo (B. S. Everitt and Skrondal 2006), multicolinearidade é um termo utilizado na análise de regressão para indicar situações em que as variáveis explicativas estão relacionadas por uma função linear, impossibilitando a estimação dos coeficientes de regressão.

Exemplo 10.15 Dalpiaz (2022) - Applied Statistics with R traz alguns exemplos sobre colinearidade.

gen_exact_collin_data <- function(num_samples = 100) {
  x1 = rnorm(n = num_samples, mean = 80, sd = 10)
  x2 = rnorm(n = num_samples, mean = 70, sd = 5)
  x3 = 2 * x1 + 4 * x2 + 3
  y = 3 + x1 + x2 + rnorm(n = num_samples, mean = 0, sd = 1)
  data.frame(y, x1, x2, x3)
}

set.seed(42)
col_data <- gen_exact_collin_data()
head(col_data)
##          y       x1       x2       x3
## 1 170.7135 93.70958 76.00483 494.4385
## 2 152.9106 74.35302 75.22376 452.6011
## 3 152.7866 83.63128 64.98396 430.1984
## 4 170.6306 86.32863 79.24241 492.6269
## 5 152.3320 84.04268 66.66613 437.7499
## 6 151.3155 78.93875 70.52757 442.9878
pairs(col_data, col = "dodgerblue")

cor(col_data)
##            y         x1         x2        x3
## y  1.0000000 0.91124848 0.43068603 0.9557550
## x1 0.9112485 1.00000000 0.03127984 0.7638609
## x2 0.4306860 0.03127984 1.00000000 0.6689586
## x3 0.9557550 0.76386089 0.66895856 1.0000000
symnum(cor(col_data))
##    y x1 x2 x3
## y  1         
## x1 * 1       
## x2 .    1    
## x3 B ,  ,  1 
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1
col_data_fit <- lm(y ~ x1 + x2 + x3, data = col_data)
summary(col_data_fit)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = col_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.57662 -0.66188 -0.08253  0.63706  2.52057 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.957336   1.735165   1.704   0.0915 .  
## x1          0.985629   0.009788 100.702   <2e-16 ***
## x2          1.017059   0.022545  45.112   <2e-16 ***
## x3                NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.014 on 97 degrees of freedom
## Multiple R-squared:  0.9923, Adjusted R-squared:  0.9921 
## F-statistic:  6236 on 2 and 97 DF,  p-value: < 2.2e-16
X <- cbind(1, as.matrix(col_data[,-1]))

fit1 <- lm(y ~ x1 + x2, data = col_data)
fit2 <- lm(y ~ x1 + x3, data = col_data)
fit3 <- lm(y ~ x2 + x3, data = col_data)
all.equal(fitted(fit1), fitted(fit2))
## [1] TRUE
all.equal(fitted(fit2), fitted(fit3))
## [1] TRUE

10.6.3.1 Pacote mctest

O pacote mctest (Imdad et al. 2019) pode ser utilizado para calcular medidas de diagnóstico de multicolinearidade populares e amplamente utilizadas.

Exemplo 10.16 (Woods, Steinour, and Starke 1932) estudam a associação do calor evoluído durante a fabricação de 13 misturas de cimento de quatro ingredientes básicos. Cada porcentagem de ingrediente parece ser arredondada para um número inteiro. A soma das quatro porcentagens de mistura varia 95% a 99%. Se todas as quatro variáveis \(X\) do regressor sempre somassem 100%, a matriz \(X\) centrada seria então de posto 3. A primeira coluna é a resposta e as quatro colunas restantes são os preditores.

Y: Calor (cals/gm) evoluído na composição.
X1: Porcentagem inteira de \(3Ca O Al_2 O_3\) na mistura.
X2: Porcentagem inteira de \(3Ca O Si O_2\) na mistura.
X3: Porcentagem inteira de \(4Ca O Al_2 O_3 Fe_2 O_3\) na mistura.
X4: Porcentagem inteira de \(2Ca O Si O_2\) na mistura.

library(mctest)
data(Hald) ## Hald Cement data
model <- lm(y ~ ., data = as.data.frame(Hald))
## Medidas de diagnóstico geral e autovalores com intercepto
mctest(model)
## 
## Call:
## omcdiag(mod = mod, Inter = TRUE, detr = detr, red = red, conf = conf, 
##     theil = theil, cn = cn)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0011         1
## Farrar Chi-Square:        67.2825         1
## Red Indicator:             0.5414         1
## Sum of Lambda Inverse:   622.3006         1
## Theil's Method:            0.9981         1
## Condition Number:        249.5783         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test
## Medidas de diagnóstico individuais
mctest(model, type = 'i')
## 
## Call:
## imcdiag(mod = mod, method = method, corr = FALSE, vif = vif, 
##     tol = tol, conf = conf, cvif = cvif, ind1 = ind1, ind2 = ind2, 
##     leamer = leamer, all = all)
## 
## 
## All Individual Multicollinearity Diagnostics Result
## 
##         VIF    TOL       Wi        Fi Leamer    CVIF Klein   IND1   IND2
## X1  38.4962 0.0260 112.4886  187.4811 0.1612 -0.5846     0 0.0087 0.9875
## X2 254.4232 0.0039 760.2695 1267.1158 0.0627 -3.8635     1 0.0013 1.0099
## X3  46.8684 0.0213 137.6052  229.3419 0.1461 -0.7117     0 0.0071 0.9923
## X4 282.5129 0.0035 844.5386 1407.5643 0.0595 -4.2900     1 0.0012 1.0103
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test
## 
## X1 , X2 , X3 , X4 , coefficient(s) are non-significant may be due to multicollinearity
## 
## R-square of y on all x: 0.9824 
## 
## * use method argument to check which regressors may be the reason of collinearity
## ===================================

10.6.3.2 Fator de Inflação da Variância (VIF)

(Chatterjee and Hadi 2012, 248–51) definem o Fator de Inflação da Variância (VIF - Variance Inflation Factor) para a variável \(X_j\) tal que

\[\begin{equation} VIF_j = \frac{1}{1-r_{j}^2}, \; j=1,\ldots,p \tag{10.27} \end{equation}\]

onde \(r_{j}^2\) é o coeficiente de determinação que resulta quando a variável preditora \(X_j\) é regredida contra todas as outras variáveis preditoras. Se \(X_j\) tiver um forte relacionamento linear com as outras variáveis preditoras, \(r_{j}^2\) deve ficar próximo de 1 e \(VIF_j\) seria grande. Como regra de bolso, valores de \(VIF\) maiores que 10 geralmente são tomados como um sinal de que os dados têm problema de colinearidade.

?mctest::mc.plot
?mctest::imcdiag(fit1)

10.6.3.3 Técnica de Farrar and Glauber

(Farrar and Glauber 1967)

10.6.3.4 Diagnóstico de Belsley, Kuh e Welsch

(Belsley, Kuh, and Welsch 2004, 112–52) sugerem um procedimento de diagnóstico para detectar colinearidade em matrizes \(X\).

Referências

Akaike, Hirotugu. 1974. “A New Look at the Statistical Model Identification.” In Selected Papers of Hirotugu Akaike, 215–22. Springer.
Belsley, David A, Edwin Kuh, and Roy E Welsch. 2004. Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. John Wiley & Sons.
Chatterjee, Samprit, and Ali S Hadi. 2012. Regression Analysis by Example. John Wiley & Sons.
Cortez, Paulo, and Anı́bal de Jesus Raimundo Morais. 2007. “A Data Mining Approach to Predict Forest Fires Using Meteorological Data.” http://www3.dsi.uminho.pt/pcortez/fires.pdf.
Efroymson, MA. 1960. “Multiple Regression Analysis.” Mathematical Methods for Digital Computers, 191–203.
Enders, Walter. 2014. Applied Econometric Time Series. 4th ed. Wiley Series in Probability and Statistics.
Everitt, Brian S, and Anders Skrondal. 2006. “The Cambridge Dictionary of Statistics.” https://www.cambridge.org/9780521766999.
Farrar, Donald E, and Robert R Glauber. 1967. “Multicollinearity in Regression Analysis: The Problem Revisited.” The Review of Economic and Statistics, 92–107. https://www.jstor.org/stable/pdf/1937887.pdf.
Friedman, Jerome H, and Werner Stuetzle. 1981. “Projection Pursuit Regression.” Journal of the American Statistical Association 76 (376): 817–23. https://doi.org/10.2307/2287576.
Hastie, Trevor, and Robert Tibshirani. 1986. “[Generalized Additive Models]: Comment.” Statistical Science 1 (3): 297–318.
Imdad, M. U., M. Aslam, S. Altaf, and A. Munir. 2019. “Some New Diagnostics of Multicollinearity in Linear Regression Model.” Sains Malaysiana 48(9): 2051–60. http://dx.doi.org/10.17576/jsm-2019-4809-26.
McCullagh, Peter, and John Ashworth Nelder. 1989. Generalized Linear Models. Chapman Hall, London. 2nd ed. http://www.utstat.toronto.edu/~brunner/oldclass/2201s11/readings/glmbook.pdf.
Modigliani, Franco, and Lucas Papademos. 1975. “Targets for Monetary Policy in the Coming Year.” Brookings Papers on Economic Activity 1975 (1): 141–65. http://pombo.free.fr/modipapa.pdf.
Neter, John, Michael H Kutner, Christopher J Nachtsheim, and William Wasserman. 2005. Applied Linear Statistical Models. 5th ed. McGraw Hill/Irwin New York. https://www.ime.unicamp.br/~dias/John%20Neter%20Applied%20linear%20regression%20models.pdf.
O’Hagan, Anthony, and Jonathan Forster. 1994. “Kendall’s Advanced Theory of Statistics: Volume 2B.” Bayesian Inference.
Paulino, Carlos Daniel Mimoso, Maria Antónia Amaral Turkman, and Bento Murteira. 2003. Estatı́stica Bayesiana. Fundação Calouste Gulbenkian, Lisboa. http://primo-pmtna01.hosted.exlibrisgroup.com/PUC01:PUC01:puc01000334509.
Schwarz, Gideon. 1978. “Estimating the Dimension of a Model.” The Annals of Statistics, 461–64. https://projecteuclid.org/journals/annals-of-statistics/volume-6/issue-2/Estimating-the-Dimension-of-a-Model/10.1214/aos/1176344136.full.
Tsanas, Athanasios, and Angeliki Xifara. 2012. “Accurate Quantitative Estimation of Energy Performance of Residential Buildings Using Statistical Machine Learning Tools.” Energy and Buildings 49: 560–67. http://people.maths.ox.ac.uk/tsanas/Preprints/ENB2012.pdf.
Woods, Hubert, Harold H Steinour, and Howard R Starke. 1932. “Effect of Composition of Portland Cement on Heat Evolved During Hardening.” Industrial & Engineering Chemistry 24 (11): 1207–14. https://pubs.acs.org/doi/pdf/10.1021/ie50275a002.

  1. In the stepwise procedure, to be described below, intermediate results, which are not even recorded by normal calculation methods, are used to give valuable statistical information at each step in the calculation. These intermediate answers are also used to control the method of calculation.↩︎

  2. An important property of the stepwise procedure is based on the facts that (a) a variable may be indicated to be significant in any early stage and thus enter the equation, and (b) after several other variables are added to the regression equation, the initial variable may be indicated to be insignificant. The insignificant variable will be removed from the regression equation before adding an additional variable. Therefore, only significant variables are included in the final regression.↩︎