Previsão Numérica Usando Regressão Linear
Heitor Gabriel S. Monteiro
28/09/2021
1 Prelúdio
Carregando os pacotes necessários.
library(tidyverse)
library(readxl) # importar arquivos do excel
library(performance) # para ver os vesultados do MQO
library(GGally) # visualizações
library(ggExtra) # plots marginais
library(kableExtra) # tabelas2 Importação e Tratamento
Importo os dados a partir de um arquivo em excel .xlxs, mudo os títulos das variáveis e transformo as categóricas em fatores.
dt <- read_excel("~/Área de Trabalho/R Projects/Análise Macro/Aula 8 - Previsão numérica usando métodos de regressão/insurance3.xlsx",
col_names = FALSE) %>% as_tibble()
attach(dt)
dt <- dt %>% rename('age' = ...1,
'sex' = ...2,
'bmi' = ...3,
'children' = ...4,
'smoker' = ...5,
'region' = ...6,
'expenses' = ...7) %>%
mutate( sex = as_factor(sex),
smoker = as_factor(smoker),
region = as_factor(region))3 Estatísticas Descritivas
Agora, vemos a estrutura dos dados tratados com o comando glimpse. Pelo comando summary, vemos que os dados estão bem distribuídos entre sexo e região, com prevalência de não fumantes. Quando dividimos entre grupos, segundo sexo, tabagismo e região, percebemos uma idade média bem concentrada entre 37 e 40 anos; um IMC também bem concentrado entre 27 e 34. A variável sobre o tabagismo parece influir muito no grupo de médias de gastos médicos.
glimpse(dt)## Rows: 1,338
## Columns: 7
## $ age <dbl> 19, 18, 28, 33, 32, 31, 46, 37, 37, 60, 25, 62, 23, 56, 27, 1…
## $ sex <fct> female, male, male, male, male, female, female, female, male,…
## $ bmi <dbl> 27.9, 33.8, 33.0, 22.7, 28.9, 25.7, 33.4, 27.7, 29.8, 25.8, 2…
## $ children <dbl> 0, 1, 3, 0, 0, 0, 1, 3, 2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0…
## $ smoker <fct> yes, no, no, no, no, no, no, no, no, no, no, yes, no, no, yes…
## $ region <fct> southwest, southeast, southeast, northwest, northwest, southe…
## $ expenses <dbl> 16884.92, 1725.55, 4449.46, 21984.47, 3866.86, 3756.62, 8240.…summary(dt)## age sex bmi children smoker
## Min. :18.00 female:662 Min. :16.00 Min. :0.000 yes: 274
## 1st Qu.:27.00 male :676 1st Qu.:26.30 1st Qu.:0.000 no :1064
## Median :39.00 Median :30.40 Median :1.000
## Mean :39.21 Mean :30.67 Mean :1.095
## 3rd Qu.:51.00 3rd Qu.:34.70 3rd Qu.:2.000
## Max. :64.00 Max. :53.10 Max. :5.000
## region expenses
## southwest:325 Min. : 1122
## southeast:364 1st Qu.: 4740
## northwest:325 Median : 9382
## northeast:324 Mean :13270
## 3rd Qu.:16640
## Max. :63770dt %>%
group_by(sex, smoker, region) %>%
summarise(Age = mean(age),
Children = mean(children),
BMI = mean(bmi),
Expenses = mean(expenses)) %>%
kbl(
#format = "latex",
booktabs = TRUE,
escape = FALSE,
digits = 1,
caption = "Média de Idade, Número de Crianças, IMC e Despesas por Grupos de Características") %>%
kable_styling(latex_options = c("striped", "hold_position"),
position = "center",
full_width = F,
bootstrap_options = c("striped", "hover", "condensed", "responsive"))| sex | smoker | region | Age | Children | BMI | Expenses |
|---|---|---|---|---|---|---|
| female | yes | southwest | 37.0 | 1.0 | 30.1 | 31688.0 |
| female | yes | southeast | 39.2 | 0.9 | 32.3 | 33034.8 |
| female | yes | northwest | 38.8 | 0.8 | 28.3 | 29670.8 |
| female | yes | northeast | 38.7 | 1.2 | 27.3 | 28032.0 |
| female | no | southwest | 40.1 | 1.1 | 30.1 | 8234.1 |
| female | no | southeast | 39.1 | 1.1 | 32.8 | 8440.2 |
| female | no | northwest | 39.8 | 1.2 | 29.5 | 8787.0 |
| female | no | northeast | 39.8 | 1.0 | 29.8 | 9640.4 |
| male | yes | southwest | 35.6 | 1.3 | 31.5 | 32598.9 |
| male | yes | southeast | 40.1 | 1.0 | 33.7 | 36029.8 |
| male | yes | northwest | 39.8 | 1.7 | 30.0 | 30713.2 |
| male | yes | northeast | 37.9 | 0.9 | 29.6 | 30926.3 |
| male | no | southwest | 40.3 | 1.1 | 31.0 | 7778.9 |
| male | no | southeast | 38.3 | 1.1 | 34.1 | 7609.0 |
| male | no | northwest | 38.6 | 1.1 | 28.9 | 8320.7 |
| male | no | northeast | 39.2 | 1.1 | 28.9 | 8664.0 |
Aseguir, fazemos três distribuições de densidade de gastos médicos para cada variável categórica. Constatamos que a variável de tabagismo influi muito na diferença de distribuição dos gastos.
Relação Entre Sexo e Gastos Médicos
Relação Entre Tabagismo e Gastos Médicos
Relação Entre Região e Gastos Médicos
Agora, fiz um correlograma estilizado que mostra uma correlação significante de idade e IMC para a nossa variável de interesse.
Relação Entre as Demais Variáveis Quantitativas e os Gastos Médicos
4 Regressão Linear Normal
Antes dos testes entre dois possíveis modelos, já que queremos somente prever os gastos, e não explicá-lo, pouco nos importa parcimônia ou causalidade. Qualquer informação que eleve nosso \(R^2 \ \ ajustado\) entrará no modelo. Portanto, é esperado usarmos todo o conjunto de variáveis disponívels.
reg1 <- lm(expenses ~ ., data = dt)
summary(reg1)##
## Call:
## lm(formula = expenses ~ ., data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11302.7 -2850.9 -979.6 1383.9 29981.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10946.61 1069.95 10.231 < 2e-16 ***
## age 256.84 11.90 21.586 < 2e-16 ***
## sexmale -131.35 332.94 -0.395 0.693255
## bmi 339.29 28.60 11.864 < 2e-16 ***
## children 475.69 137.80 3.452 0.000574 ***
## smokerno -23847.48 413.14 -57.723 < 2e-16 ***
## regionsoutheast -76.29 470.64 -0.162 0.871253
## regionnorthwest 606.52 477.18 1.271 0.203940
## regionnortheast 959.31 477.91 2.007 0.044921 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6062 on 1329 degrees of freedom
## Multiple R-squared: 0.7509, Adjusted R-squared: 0.7494
## F-statistic: 500.9 on 8 and 1329 DF, p-value: < 2.2e-16Vemos que eliminar as variáveis não significantes diminui marginalmente o \(R^2\), então, manteremos as variáveis não-significantes, para fins de previsão.
reg2 <- lm(expenses ~ age + bmi + smoker, data = dt)
summary(reg2)##
## Call:
## lm(formula = expenses ~ age + bmi + smoker, data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12413 -2970 -977 1476 28960
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12143.56 987.12 12.30 <2e-16 ***
## age 259.53 11.93 21.75 <2e-16 ***
## bmi 322.69 27.49 11.74 <2e-16 ***
## smokerno -23822.61 412.86 -57.70 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6092 on 1334 degrees of freedom
## Multiple R-squared: 0.7475, Adjusted R-squared: 0.7469
## F-statistic: 1316 on 3 and 1334 DF, p-value: < 2.2e-16Ainda sim, mesmo com o objetivo de prever, vamos fazer um diagnóstico dos modelos aqui construídos, por exercício mesmo.
compare_performance(reg1, reg2) %>% kbl()| Name | Model | AIC | BIC | R2 | R2_adjusted | RMSE | Sigma |
|---|---|---|---|---|---|---|---|
| reg1 | lm | 27115.42 | 27167.41 | 0.7509285 | 0.7494292 | 6041.493 | 6061.915 |
| reg2 | lm | 27123.76 | 27149.76 | 0.7474907 | 0.7469228 | 6083.043 | 6092.156 |
check_model(reg1)
check_model(reg2)
check_autocorrelation(reg1)## OK: Residuals appear to be independent and not autocorrelated (p = 0.104).check_autocorrelation(reg2)## OK: Residuals appear to be independent and not autocorrelated (p = 0.170).check_heteroscedasticity(reg1)## Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).check_heteroscedasticity(reg2)## Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).check_normality(reg1)## Warning: Non-normality of residuals detected (p < .001).check_normality(reg2)## Warning: Non-normality of residuals detected (p < .001).test_performance(reg1, reg2) %>% kbl()| Name | Model | log_BF | BF |
|---|---|---|---|
| reg1 | lm | NA | NA |
| reg2 | lm | 8.826684 | 6813.654 |
5 Treinamento
Divido, usando uma amostra aleatória, entre amostra de treino (com 80%) e de teste (20%). Como as quantidades entre sexo, região, e até mesmo a média das variáveis quantitativas são parecidas entre os grupos, exceto entre fumantes e não-fumantes, esperamos duas amostras bem parecidas entre si.
set.seed(777)
train <- dt %>% sample_frac(.,0.8)
sid <- as.numeric(rownames(train))
test <- dt[-sid,]
remove(sid)Vamos treinar o modelo, chmando-o de reg3. Percebemos pelo diagnóstico que ele é realmente parecido em distribuição com o reg1, feito anteriormente.
reg3 <- lm(expenses ~ ., data = train)
summary(reg3)##
## Call:
## lm(formula = expenses ~ ., data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11057.9 -2895.2 -904.1 1535.4 30002.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10467.26 1200.49 8.719 <2e-16 ***
## age 252.49 13.58 18.588 <2e-16 ***
## sexmale -83.46 375.54 -0.222 0.8242
## bmi 357.53 31.89 11.212 <2e-16 ***
## children 415.51 155.11 2.679 0.0075 **
## smokerno -23893.66 462.80 -51.628 <2e-16 ***
## regionsoutheast -23.96 538.83 -0.044 0.9645
## regionnorthwest 1010.71 538.30 1.878 0.0607 .
## regionnortheast 1232.37 542.57 2.271 0.0233 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6120 on 1061 degrees of freedom
## Multiple R-squared: 0.7506, Adjusted R-squared: 0.7487
## F-statistic: 399.2 on 8 and 1061 DF, p-value: < 2.2e-16check_model(reg3)
gg4 <- data.frame(Estimado = reg3$fitted.values,
Verdadeiro = train$expenses
) %>%
ggplot(aes(y=Verdadeiro, x=Estimado)) +
geom_point(alpha=.3) +
geom_abline(slope = 1,
intercept = 0,
color = 'red3')
ggExtra::ggMarginal(gg4,
type = 'density',
margins = 'both',
colour = 'firebrick',
fill = 'tomato3')
6 Teste
pred <- predict(reg3, newdata = test)
dt_pred <- data.frame(Estimado = pred,
Verdadeiro = test$expenses)
gg_pred <- dt_pred %>%
ggplot(aes(y = Verdadeiro,
x = Estimado)) +
geom_point(alpha=.3) +
geom_abline(slope = 1,
intercept = 0,
color = 'deepskyblue')
ggExtra::ggMarginal(gg_pred,
type = 'density',
margins = 'both',
colour = 'skyblue2',
fill = 'lightblue1')