Objetivo

Objetivo é executar uma analise afim de buscar relações entre as variáveis e conseguir um melhor ajuste para fazer determinadas predições , para esta análise vamos utilizar o pacote RStanArm e o software Rstúdio

Metodologia

para está análise utilizaremos banco de dados Automobile Data Set contendo informações relevantes sobre características dos carros Variáveis:

Definindo o Problema de Negócio

Nosso objetivo é construir um modelo de regressão generalizado que seja capaz de explicar quais características influenciam na precificação do automovel. A variável a ser prevista é um valor numérico que representa o preço de cada automovel observado. Para cada automóvel temos diversas variáveis explanatórias. Sendo assim, podemos buscar resolver este problema utilizando MLG.

Definindo o Dataset

Este conjunto de dados consiste em dados do Anuário Automotivo da Ward de 1985. Aqui estão as fontes

Fontes:

  1. 1985 Modelo de Importação de Carro e Especificações de Caminhão, Anuário Automotivo de 1985 Ward.
  2. Personal Auto Manuals, Insurance Services Office, 160 Water Street, Nova York, NY 10038
  3. Insurance Collision Report, Insurance Institute for Highway Safety, Watergate 600, Washington, DC 20037

Esse conjunto de dados consiste em três tipos de entidades: (a) a especificação de um automóvel em termos de várias características, (b) sua classificação de risco de seguro atribuída, (c) suas perdas normalizadas em uso em comparação com outros carros. A segunda classificação corresponde ao grau em que o automóvel é mais arriscado do que seu preço indica. Os carros recebem inicialmente um símbolo de fator de risco associado ao seu preço. Então, se for mais arriscado (ou menos), esse símbolo é ajustado movendo-o para cima (ou para baixo) na escala. Os atuários chamam esse processo de “simbolização”. Um valor de +3 indica que o automóvel é arriscado, -3 que provavelmente é bastante seguro.

O terceiro fator é o pagamento de perda média relativa por ano de veículo segurado. Este valor é normalizado para todos os automóveis dentro de uma determinada classificação de tamanho (pequeno de duas portas, carrinhas, desporto/especialidade, etc…), e representa a perda média por automóvel por ano.

Nota: Vários dos atributos no banco de dados podem ser usados como um atributo de “classe”.

  • total de variáveis: 26
  • total de observações :205

Introdução

O primeiro meio de transporte a fazer uso de um motor a gasolina para se locomover foi um automóvel que continha somente três rodas e foi criado no ano de 1885 por um alemão de nome Karl Benz.A partir de então teve início a corrida pela produção e venda de automóveis, iniciada por uma empresa francesa conhecida pelo nome de Panhard et Levassor. No ano de 1892, o conhecido Henry Ford fabricou seu primeiro carro, o Ford, na América do Norte.Estima-se que a quantidade de autóveis no Brasil já ultrapassa a casa dos 11 milhões de veículos, segundo o senado, Dados da Secretaria Nacional de Trânsito do Ministério da Infraestrutura indicam haver mais de 3,5 milhões de caminhões em circulação no Brasil e, desse total, cerca de 26% dos veículos possuem mais de 30 anos de fabricação, a expectativa é que esse número aumente ainda mais nós próximos anos. Neste material buscaremos criar um modelo capaz de explicar quais fatores influenciar no preço do automóvel, para este fim sera utilizado os modelos lineares generalizados.

Os modelos lineares Generalizados são utilizados quando utilizar os modelos clássicos que pressupoem que os dados seguem uma distribuição normal e principalmente onde temos dados que segue uma base de dados discretos . O princípio do GLM é pegar uma família de distribuições ao ínves de uma única distribuição que no nosso caso será a famíla exponencial , ou sejá , ao ínves de trabalhar apenas com a normal , trabalhamos com várias distribuições : gamma,expenencial,poisson ,binomial,binomial negativa,weibell. Um estudo de regressão busca,essencialmente,associar uma variável Y (denominada variável dependente)a uma outra variável X (denominada variável explanatória ou variável independente ).

Partes de um MLG

  1. \(Y_{i} \sim p(y | \theta) \rightarrow E(y_{i}) = \mu_{i}\)

  2. \(n_{i} = g(\mu_{i})\)

  3. \(n_{i} = \beta_{0} + \beta_{1} *x_{I}\)

Conhecendo o Nossos Dados

Análises Descritivas

Banco de Dados

library(readr)
dados <- read_csv("carrosTrat.csv", 
    col_types = cols(stroke = col_double()))
## Warning: One or more parsing issues, see `problems()` for details
dados$price<-as.numeric(dados$price)
head (dados)
## # A tibble: 6 × 26
##   symboling normalized.losses make  fuel.type aspiration num.of.doors body.style
##       <dbl>             <dbl> <chr> <chr>     <chr>      <chr>        <chr>     
## 1         3               115 alfa… gas       std        two          convertib…
## 2         3               115 alfa… gas       std        two          convertib…
## 3         1               115 alfa… gas       std        two          hatchback 
## 4         2               164 audi  gas       std        four         sedan     
## 5         2               164 audi  gas       std        four         sedan     
## 6         2               115 audi  gas       std        two          sedan     
## # … with 19 more variables: drive.wheels <chr>, engine.location <chr>,
## #   wheel.base <dbl>, length <dbl>, width <dbl>, height <dbl>,
## #   curb.weight <dbl>, engine.type <chr>, num.of.cylinders <chr>,
## #   engine.size <dbl>, fuel.system <chr>, bore <dbl>, stroke <dbl>,
## #   compression.ratio <dbl>, horsepower <dbl>, peak.rpm <dbl>, city.mpg <dbl>,
## #   highway.mpg <dbl>, price <dbl>

Dados Faltantes NA’s

round(mean(is.na(dados))*100,10)
## [1] 0.0750469

Tem a presença de 0.075% de dados faltantes, analisando onde ele se encontra:

library(skimr)
## Warning: package 'skimr' was built under R version 4.2.1
skim(dados)
Data summary
Name dados
Number of rows 205
Number of columns 26
_______________________
Column type frequency:
character 10
numeric 16
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
make 0 1 3 13 0 22 0
fuel.type 0 1 3 6 0 2 0
aspiration 0 1 3 5 0 2 0
num.of.doors 0 1 1 4 0 3 0
body.style 0 1 5 11 0 5 0
drive.wheels 0 1 3 3 0 3 0
engine.location 0 1 4 5 0 2 0
engine.type 0 1 1 5 0 7 0
num.of.cylinders 0 1 3 6 0 7 0
fuel.system 0 1 3 4 0 8 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
symboling 0 1.00 0.83 1.25 -2.00 0.00 1.00 2.00 3.00 ▃▇▆▃▃
normalized.losses 0 1.00 120.60 31.81 65.00 101.00 115.00 137.00 256.00 ▆▇▃▁▁
wheel.base 0 1.00 98.76 6.02 86.60 94.50 97.00 102.40 120.90 ▁▇▂▁▁
length 0 1.00 174.05 12.34 141.10 166.30 173.20 183.10 208.10 ▁▅▇▃▁
width 0 1.00 65.91 2.15 60.30 64.10 65.50 66.90 72.30 ▁▇▇▂▁
height 0 1.00 53.72 2.44 47.80 52.00 54.10 55.50 59.80 ▁▆▇▆▂
curb.weight 0 1.00 2555.57 520.68 1488.00 2145.00 2414.00 2935.00 4066.00 ▃▇▅▃▁
engine.size 0 1.00 126.91 41.64 61.00 97.00 120.00 141.00 326.00 ▇▆▂▁▁
bore 0 1.00 3.33 0.27 2.54 3.15 3.31 3.58 3.94 ▁▅▇▇▂
stroke 4 0.98 3.26 0.32 2.07 3.11 3.29 3.41 4.17 ▁▂▇▇▁
compression.ratio 0 1.00 10.14 3.97 7.00 8.60 9.00 9.40 23.00 ▇▁▁▁▁
horsepower 0 1.00 104.17 39.53 48.00 70.00 95.00 116.00 288.00 ▇▅▂▁▁
peak.rpm 0 1.00 5126.10 477.04 4150.00 4800.00 5200.00 5500.00 6600.00 ▂▇▇▂▁
city.mpg 0 1.00 25.22 6.54 13.00 19.00 24.00 30.00 49.00 ▆▇▅▂▁
highway.mpg 0 1.00 30.75 6.89 16.00 25.00 30.00 34.00 54.00 ▂▇▆▁▁
price 0 1.00 13150.31 7879.12 5118.00 7788.00 10295.00 16500.00 45400.00 ▇▃▁▁▁

Observe que os dados faltantes se encontram em STROKE (derramamento de óleo), eles precisam ser tratados:

dados[is.na(dados$stroke),]$stroke = median(dados$stroke,na.rm = T)

Foi atribuido para cada valor faltante a sua mediana da coluna, agora poderá seguir com o restante da análise dos dados.

Histograma e Normalidade

library(ggplot2)
dados %>% ggplot(aes(price))+
  geom_density(fill = "darkblue")+
  labs(x= "Preço", y="Densidade")

Observe que os dados de preço são assimétricos a esquerda.

shapiro.test(dados$price)
## 
##  Shapiro-Wilk normality test
## 
## data:  dados$price
## W = 0.79542, p-value = 1.145e-15

Pelo teste de normalidade univariada de Shapiro-Wilk, o valor deu menor que 0.05, então os dados de preço não seguem a normalidade dos resíduos.

library(skimr)
skim(dados)
Data summary
Name dados
Number of rows 205
Number of columns 26
_______________________
Column type frequency:
character 10
numeric 16
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
make 0 1 3 13 0 22 0
fuel.type 0 1 3 6 0 2 0
aspiration 0 1 3 5 0 2 0
num.of.doors 0 1 1 4 0 3 0
body.style 0 1 5 11 0 5 0
drive.wheels 0 1 3 3 0 3 0
engine.location 0 1 4 5 0 2 0
engine.type 0 1 1 5 0 7 0
num.of.cylinders 0 1 3 6 0 7 0
fuel.system 0 1 3 4 0 8 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
symboling 0 1 0.83 1.25 -2.00 0.00 1.00 2.00 3.00 ▃▇▆▃▃
normalized.losses 0 1 120.60 31.81 65.00 101.00 115.00 137.00 256.00 ▆▇▃▁▁
wheel.base 0 1 98.76 6.02 86.60 94.50 97.00 102.40 120.90 ▁▇▂▁▁
length 0 1 174.05 12.34 141.10 166.30 173.20 183.10 208.10 ▁▅▇▃▁
width 0 1 65.91 2.15 60.30 64.10 65.50 66.90 72.30 ▁▇▇▂▁
height 0 1 53.72 2.44 47.80 52.00 54.10 55.50 59.80 ▁▆▇▆▂
curb.weight 0 1 2555.57 520.68 1488.00 2145.00 2414.00 2935.00 4066.00 ▃▇▅▃▁
engine.size 0 1 126.91 41.64 61.00 97.00 120.00 141.00 326.00 ▇▆▂▁▁
bore 0 1 3.33 0.27 2.54 3.15 3.31 3.58 3.94 ▁▅▇▇▂
stroke 0 1 3.26 0.31 2.07 3.11 3.29 3.41 4.17 ▁▂▇▇▁
compression.ratio 0 1 10.14 3.97 7.00 8.60 9.00 9.40 23.00 ▇▁▁▁▁
horsepower 0 1 104.17 39.53 48.00 70.00 95.00 116.00 288.00 ▇▅▂▁▁
peak.rpm 0 1 5126.10 477.04 4150.00 4800.00 5200.00 5500.00 6600.00 ▂▇▇▂▁
city.mpg 0 1 25.22 6.54 13.00 19.00 24.00 30.00 49.00 ▆▇▅▂▁
highway.mpg 0 1 30.75 6.89 16.00 25.00 30.00 34.00 54.00 ▂▇▆▁▁
price 0 1 13150.31 7879.12 5118.00 7788.00 10295.00 16500.00 45400.00 ▇▃▁▁▁

Análise Exploratória

O ponta pé inicial de qualquer análise é a análise exploratória, com ela começas a entender a grandeza dos nossos dados e os seus comportamentos, dessa forma, vamos carregar os pacotes necessários e executar a análise exploratória:

library(tidyverse)
library(DT)
library(readxl)


dados<-read.table("dados.txt",sep = ",",header = TRUE)

Vamos verificar se nosso conjunto de dados pussui valores ausentes, caso tenhamos, teremos que remover essas observações

sum(dados[!complete.cases(dados),])
## [1] 0

Como podemos ver o nosso conjunto de dados não apresentou nenhum valor faltante, logo, podemos prosseguir da forma que está.

Observando a Grandeza dos Nossos Dados

Para verificar a grandeza dos nossos dados vamos utilizar o seguinte comando:

glimpse(dados)
## Rows: 205
## Columns: 26
## $ symboling         <int> 3, 3, 1, 2, 2, 2, 1, 1, 1, 0, 2, 0, 0, 0, 1, 0, 0, 0…
## $ normalized.losses <chr> "?", "?", "?", "164", "164", "?", "158", "?", "158",…
## $ make              <chr> "alfa-romero", "alfa-romero", "alfa-romero", "audi",…
## $ fuel.type         <chr> "gas", "gas", "gas", "gas", "gas", "gas", "gas", "ga…
## $ aspiration        <chr> "std", "std", "std", "std", "std", "std", "std", "st…
## $ num.of.doors      <chr> "two", "two", "two", "four", "four", "two", "four", …
## $ body.style        <chr> "convertible", "convertible", "hatchback", "sedan", …
## $ drive.wheels      <chr> "rwd", "rwd", "rwd", "fwd", "4wd", "fwd", "fwd", "fw…
## $ engine.location   <chr> "front", "front", "front", "front", "front", "front"…
## $ wheel.base        <dbl> 88.6, 88.6, 94.5, 99.8, 99.4, 99.8, 105.8, 105.8, 10…
## $ length            <dbl> 168.8, 168.8, 171.2, 176.6, 176.6, 177.3, 192.7, 192…
## $ width             <dbl> 64.1, 64.1, 65.5, 66.2, 66.4, 66.3, 71.4, 71.4, 71.4…
## $ height            <dbl> 48.8, 48.8, 52.4, 54.3, 54.3, 53.1, 55.7, 55.7, 55.9…
## $ curb.weight       <int> 2548, 2548, 2823, 2337, 2824, 2507, 2844, 2954, 3086…
## $ engine.type       <chr> "dohc", "dohc", "ohcv", "ohc", "ohc", "ohc", "ohc", …
## $ num.of.cylinders  <chr> "four", "four", "six", "four", "five", "five", "five…
## $ engine.size       <int> 130, 130, 152, 109, 136, 136, 136, 136, 131, 131, 10…
## $ fuel.system       <chr> "mpfi", "mpfi", "mpfi", "mpfi", "mpfi", "mpfi", "mpf…
## $ bore              <chr> "3.47", "3.47", "2.68", "3.19", "3.19", "3.19", "3.1…
## $ stroke            <chr> "2.68", "2.68", "3.47", "3.40", "3.40", "3.40", "3.4…
## $ compression.ratio <dbl> 9.00, 9.00, 9.00, 10.00, 8.00, 8.50, 8.50, 8.50, 8.3…
## $ horsepower        <chr> "111", "111", "154", "102", "115", "110", "110", "11…
## $ peak.rpm          <chr> "5000", "5000", "5000", "5500", "5500", "5500", "550…
## $ city.mpg          <int> 21, 21, 19, 24, 18, 19, 19, 19, 17, 16, 23, 23, 21, …
## $ highway.mpg       <int> 27, 27, 26, 30, 22, 25, 25, 25, 20, 22, 29, 29, 28, …
## $ price             <chr> "13495", "16500", "16500", "13950", "17450", "15250"…

Como podemos ver a maioria das variáveis não estão no formato adequado,vamos ter que converte-las para fatores as que são referentes a fatores e numericos a referente a valores númericos

Transformações de Variáveis

Vamos utilizar a função mutate_if para colocar a condição de transformação das colunas, por exemplo , no primeiro código abaixo as culanas de dados da coluna 3 a coluna 9 que são do tipo character serão convertidas para factor , os demais códigos abaixo seguem as mesmas linhas de raciocínio.

dados[3:9]<-dados[3:9] %>% mutate_if(is.character,as.factor)
  
dados[15:16]<-dados[15:16] %>% mutate_if(is.character,as.factor) 
  
 

dados[18]<-dados[18] %>% mutate_if(is.character,as.factor)  


dados[20]<-dados[20] %>% mutate_if(is.character,as.factor) 



dados<-dados %>% mutate_if(is.character,as.numeric)

Vamos observar os nossos dados

library(gt)
# dados |> datatable(
#           rownames = FALSE,
#           options = list(
#             columnDefs = list(list(className = 'dt-center', targets = 0:4))
#             )
#           )






dados |> head(10) |> gt()
symboling normalized.losses make fuel.type aspiration num.of.doors body.style drive.wheels engine.location wheel.base length width height curb.weight engine.type num.of.cylinders engine.size fuel.system bore stroke compression.ratio horsepower peak.rpm city.mpg highway.mpg price
3 NA alfa-romero gas std two convertible rwd front 88.6 168.8 64.1 48.8 2548 dohc four 130 mpfi 3.47 2.68 9.0 111 5000 21 27 13495
3 NA alfa-romero gas std two convertible rwd front 88.6 168.8 64.1 48.8 2548 dohc four 130 mpfi 3.47 2.68 9.0 111 5000 21 27 16500
1 NA alfa-romero gas std two hatchback rwd front 94.5 171.2 65.5 52.4 2823 ohcv six 152 mpfi 2.68 3.47 9.0 154 5000 19 26 16500
2 164 audi gas std four sedan fwd front 99.8 176.6 66.2 54.3 2337 ohc four 109 mpfi 3.19 3.40 10.0 102 5500 24 30 13950
2 164 audi gas std four sedan 4wd front 99.4 176.6 66.4 54.3 2824 ohc five 136 mpfi 3.19 3.40 8.0 115 5500 18 22 17450
2 NA audi gas std two sedan fwd front 99.8 177.3 66.3 53.1 2507 ohc five 136 mpfi 3.19 3.40 8.5 110 5500 19 25 15250
1 158 audi gas std four sedan fwd front 105.8 192.7 71.4 55.7 2844 ohc five 136 mpfi 3.19 3.40 8.5 110 5500 19 25 17710
1 NA audi gas std four wagon fwd front 105.8 192.7 71.4 55.7 2954 ohc five 136 mpfi 3.19 3.40 8.5 110 5500 19 25 18920
1 158 audi gas turbo four sedan fwd front 105.8 192.7 71.4 55.9 3086 ohc five 131 mpfi 3.13 3.40 8.3 140 5500 17 20 23875
0 NA audi gas turbo two hatchback 4wd front 99.5 178.2 67.9 52.0 3053 ohc five 131 mpfi 3.13 3.40 7.0 160 5500 16 22 NA

Podemos perceber que após as transformações temos a presença de valores ausentes,teremos que tratar esses valores ausentes

Trantando valores ausentes nas variáveis númericas

# symboling é fator ,mas vai permanecer como númerica   


# subiistituindo os valores ausentes pelo valor mediano 

dados[is.na(dados$normalized.losses),]$normalized.losses =median(dados$normalized.losses,na.rm = T)

dados[is.na(dados$price),]$price = median(dados$price,na.rm = T)


dados$stroke<-as.numeric(dados$stroke)

#dados[is.na(dados$stroke),]$stroke = median(dados$stroke,na.rm = T)


dados[is.na(dados$bore),]$bore = median(dados$bore,na.rm = T)
#dados[is.na(dados$compression.ratio),]$compression.ratio = median(dados$compression.ratio,na.rm = T)
dados[is.na(dados$horsepower),]$horsepower= median(dados$horsepower,na.rm = T)

dados[is.na(dados$peak.rpm),]$peak.rpm= median(dados$peak.rpm,na.rm = T)

Verificando se Ainda temos Algum Valor Ausente

colnms <- colnames(dados)

# filter
teste<-dados %>%
  filter_at(vars(all_of(colnms)), any_vars(is.na(.)))

sum(teste)
## [1] 0

Como podemos ver não temos mais a presença de valores ausentes, podemos prosseguir.

Vamos observar novamente nossos dados

library(gt)
# dados |> datatable(
#           rownames = FALSE,
#           options = list(
#             columnDefs = list(list(className = 'dt-center', targets = 0:4))
#             )
#           )






dados |> head(10) |> gt()
symboling normalized.losses make fuel.type aspiration num.of.doors body.style drive.wheels engine.location wheel.base length width height curb.weight engine.type num.of.cylinders engine.size fuel.system bore stroke compression.ratio horsepower peak.rpm city.mpg highway.mpg price
3 115 alfa-romero gas std two convertible rwd front 88.6 168.8 64.1 48.8 2548 dohc four 130 mpfi 3.47 6 9.0 111 5000 21 27 13495
3 115 alfa-romero gas std two convertible rwd front 88.6 168.8 64.1 48.8 2548 dohc four 130 mpfi 3.47 6 9.0 111 5000 21 27 16500
1 115 alfa-romero gas std two hatchback rwd front 94.5 171.2 65.5 52.4 2823 ohcv six 152 mpfi 2.68 29 9.0 154 5000 19 26 16500
2 164 audi gas std four sedan fwd front 99.8 176.6 66.2 54.3 2337 ohc four 109 mpfi 3.19 26 10.0 102 5500 24 30 13950
2 164 audi gas std four sedan 4wd front 99.4 176.6 66.4 54.3 2824 ohc five 136 mpfi 3.19 26 8.0 115 5500 18 22 17450
2 115 audi gas std two sedan fwd front 99.8 177.3 66.3 53.1 2507 ohc five 136 mpfi 3.19 26 8.5 110 5500 19 25 15250
1 158 audi gas std four sedan fwd front 105.8 192.7 71.4 55.7 2844 ohc five 136 mpfi 3.19 26 8.5 110 5500 19 25 17710
1 115 audi gas std four wagon fwd front 105.8 192.7 71.4 55.7 2954 ohc five 136 mpfi 3.19 26 8.5 110 5500 19 25 18920
1 158 audi gas turbo four sedan fwd front 105.8 192.7 71.4 55.9 3086 ohc five 131 mpfi 3.13 26 8.3 140 5500 17 20 23875
0 115 audi gas turbo two hatchback 4wd front 99.5 178.2 67.9 52.0 3053 ohc five 131 mpfi 3.13 26 7.0 160 5500 16 22 10295

Vendo as Correlações das Variáveis Númericas

Problemas de multicolinearidade são muitos comuns em modelos de regressão,isso acaba prejudicando o modelo, devido a este fato vamos verificar a correção entre as variáveis preditoras, e depois vamos verificar quais variaveis remover, para verificar a correlação vamos executar o gráfico de correlação de pearson através do comando abaixo:

library(corrplot)
## corrplot 0.92 loaded
par(bg = '#586573')
dados[-26] |>
  dplyr::select(where(is.numeric))  |>
  cor( ) |>
corrplot::corrplot( type = "upper")

Como podemos ver várias colunas apresentaram correlação,logo teremos que remove-lás

Verificando Presença de Multicolinearidade das Variáveis Númericas

library(caret)

t<-dados[-26] |>
  dplyr::select(where(is.numeric)) 

 t %>%
  cor() %>%
  findCorrelation( cutoff = .50, verbose = FALSE)
## [1]  7  4  5  3 15  8 14  9  6

Logo, essas as variáveis cotadas para remoção ,mas antes disso vamos ver se temos colunas com repetidas

Verificando se Temos Colunas Identicas

t %>%
  cor() %>%  findLinearCombos ()
## $linearCombos
## list()
## 
## $remove
## NULL

Como podemos ver não temos colunas repetidas , vamos remover as colunas sugeridas

t<-dados |>
  dplyr::select(where(is.numeric))  
   

novo<-t %>% 
  dplyr::select(normalized.losses,compression.ratio,horsepower,peak.rpm,price )

nov1<-dados |>
  dplyr::select(where(is.factor))  


#unido os dados 

# novo<-cbind(nov1,novo)
# 
# glimpse(novo)
# novo<-novo |>
#   dplyr::filter(make !="?" , fuel.type !="?",aspiration  !="?" ,
#                 num.of.doors !="?", engine.type !="l")


write.csv(novo,"carros2Trat.csv",row.names = F)

Agora podemos proseguir para criação do nosso modelo

Criando o Modelo MLG

O primeiro modelo irá receber todas as variáveis , o segundo modelo recebera as variaveis selecionadas , os demais modelos serão adicionando as variáveis um a um . fazendo os modelos

# nv<- novo[11:16]
# 
mod<-glm(price~.,family = Gamma(link=log),data = t)



summary(mod)
## 
## Call:
## glm(formula = price ~ ., family = Gamma(link = log), data = t)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.04513  -0.13004  -0.03686   0.09170   0.52441  
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        6.050e+00  9.775e-01   6.190 3.66e-09 ***
## symboling          1.522e-02  1.678e-02   0.907  0.36535    
## normalized.losses  6.144e-04  5.557e-04   1.106  0.27030    
## wheel.base         5.769e-03  6.921e-03   0.834  0.40556    
## length             2.649e-03  3.501e-03   0.757  0.45015    
## width              2.598e-03  1.572e-02   0.165  0.86886    
## height             8.837e-03  9.075e-03   0.974  0.33143    
## curb.weight        2.269e-04  1.101e-04   2.061  0.04066 *  
## engine.size        4.034e-03  8.865e-04   4.550 9.57e-06 ***
## bore              -2.238e-02  7.671e-02  -0.292  0.77082    
## stroke            -5.792e-03  1.752e-03  -3.306  0.00113 ** 
## compression.ratio  2.286e-02  5.260e-03   4.346 2.26e-05 ***
## horsepower         2.261e-03  1.033e-03   2.188  0.02987 *  
## peak.rpm           9.480e-05  4.269e-05   2.221  0.02755 *  
## city.mpg          -2.304e-02  1.139e-02  -2.024  0.04443 *  
## highway.mpg        9.245e-03  1.014e-02   0.912  0.36305    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.04040951)
## 
##     Null deviance: 55.8608  on 204  degrees of freedom
## Residual deviance:  7.8485  on 189  degrees of freedom
## AIC: 3780.9
## 
## Number of Fisher Scoring iterations: 8

Para determinar o melhor modelo vamos utilizar o critério informação de Akaike que é uma métrica que mensura a qualidade de um modelo estatístico visando também a sua simplicidade. Fornece, portanto, uma métrica para comparação e seleção de modelos, em que menores valores de AIC representam uma maior qualidade e simplicidade, segundo este critério.

#mod2<-glm(price~ compression.ratio+horsepower+peak.rpm,family = Gamma(link=log),data = novo)
# modelo com as variaveis selecionadas 
mod2<-glm(price~normalized.losses+compression.ratio+horsepower+peak.rpm
          ,family = Gamma(link=log),data = novo)
summary(mod2)
## 
## Call:
## glm(formula = price ~ normalized.losses + compression.ratio + 
##     horsepower + peak.rpm, family = Gamma(link = log), data = novo)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.60361  -0.18808  -0.05306   0.14854   0.67458  
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        8.475e+00  2.505e-01  33.830  < 2e-16 ***
## normalized.losses  8.700e-04  5.963e-04   1.459  0.14612    
## compression.ratio  2.807e-02  5.170e-03   5.429 1.63e-07 ***
## horsepower         1.153e-02  4.762e-04  24.212  < 2e-16 ***
## peak.rpm          -1.330e-04  4.346e-05  -3.060  0.00252 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.06771908)
## 
##     Null deviance: 55.861  on 204  degrees of freedom
## Residual deviance: 14.625  on 200  degrees of freedom
## AIC: 3887.6
## 
## Number of Fisher Scoring iterations: 6
mod3<-glm(price~horsepower+peak.rpm+compression.ratio ,family = Gamma(link=log),data = novo)
  
summary(mod3)
## 
## Call:
## glm(formula = price ~ horsepower + peak.rpm + compression.ratio, 
##     family = Gamma(link = log), data = novo)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.62086  -0.18467  -0.06676   0.13946   0.64136  
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        8.507e+00  2.510e-01  33.887  < 2e-16 ***
## horsepower         1.163e-02  4.736e-04  24.556  < 2e-16 ***
## peak.rpm          -1.207e-04  4.269e-05  -2.827  0.00517 ** 
## compression.ratio  2.803e-02  5.194e-03   5.398 1.89e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.06834542)
## 
##     Null deviance: 55.861  on 204  degrees of freedom
## Residual deviance: 14.771  on 201  degrees of freedom
## AIC: 3887.7
## 
## Number of Fisher Scoring iterations: 6
#mod4<-glm(price~compression.ratio+horsepower+peak.rpm,family = Gamma(link=log),data = novo)
  
#summary(mod4)

Verificando a Importancia das Variáveis para o Modelo

Vamos verificar quais variáveis são mais importantes para o modelo , para isso vamos utilizar o pacote e os comandos abaixo:

library(randomForest)
importancia  = randomForest(price ~ ., data = t)

col = importance(importancia)
options(scipen=999) 
par(bg = '#586573')
varImpPlot(importancia)

Podemos perceber que uma das variáveis removidas anteriormente engine.size (cilindradas do motor) possui a maior importância para o modelo , logo,vamos tentar ajustar um modelo com esssa covariável

t2<-t %>% 
  dplyr::select(engine.size,curb.weight,horsepower,highway.mpg,city.mpg,length ,price)
cor(t2)
##             engine.size curb.weight horsepower highway.mpg   city.mpg
## engine.size   1.0000000   0.8505941  0.8102160  -0.6774699 -0.6536579
## curb.weight   0.8505941   1.0000000  0.7509275  -0.7974648 -0.7574138
## horsepower    0.8102160   0.7509275  1.0000000  -0.7707804 -0.8021695
## highway.mpg  -0.6774699  -0.7974648 -0.7707804   1.0000000  0.9713370
## city.mpg     -0.6536579  -0.7574138 -0.8021695   0.9713370  1.0000000
## length        0.6833599   0.8777285  0.5533374  -0.7046616 -0.6709087
## price         0.8603427   0.8198167  0.7499191  -0.6930373 -0.6688215
##                 length      price
## engine.size  0.6833599  0.8603427
## curb.weight  0.8777285  0.8198167
## horsepower   0.5533374  0.7499191
## highway.mpg -0.7046616 -0.6930373
## city.mpg    -0.6709087 -0.6688215
## length       1.0000000  0.6865674
## price        0.6865674  1.0000000
# modelo sugerido 1
yu<-glm(price~engine.size+highway.mpg+length , family = Gamma(link=log),data = t2)
summary(yu)
## 
## Call:
## glm(formula = price ~ engine.size + highway.mpg + length, family = Gamma(link = log), 
##     data = t2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.58832  -0.16092  -0.05071   0.11449   0.77670  
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  7.5305555  0.4293509  17.539 < 0.0000000000000002 ***
## engine.size  0.0058535  0.0006155   9.511 < 0.0000000000000002 ***
## highway.mpg -0.0198601  0.0038295  -5.186          0.000000522 ***
## length       0.0098396  0.0021535   4.569          0.000008542 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.06120744)
## 
##     Null deviance: 55.861  on 204  degrees of freedom
## Residual deviance: 11.191  on 201  degrees of freedom
## AIC: 3830.2
## 
## Number of Fisher Scoring iterations: 5
# modelo sugerido 2 
yu2<-glm(price~curb.weight+length , family = Gamma(link=log),data = t2)

summary(yu2)
## 
## Call:
## glm(formula = price ~ curb.weight + length, family = Gamma(link = log), 
##     data = t2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.76966  -0.15834  -0.04829   0.06193   1.05941  
## 
## Coefficients:
##                Estimate  Std. Error t value            Pr(>|t|)    
## (Intercept)  7.50654984  0.40039726  18.748 <0.0000000000000002 ***
## curb.weight  0.00092313  0.00007622  12.112 <0.0000000000000002 ***
## length      -0.00280031  0.00321675  -0.871               0.385    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.07376687)
## 
##     Null deviance: 55.861  on 204  degrees of freedom
## Residual deviance: 12.437  on 202  degrees of freedom
## AIC: 3850
## 
## Number of Fisher Scoring iterations: 4

Modelo Binomial Negativo

require(MASS)
mod.nb <- glm.nb(price~., data = novo)

summary(mod.nb)
## 
## Call:
## glm.nb(formula = price ~ ., data = novo, init.theta = 14.18961524, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.0389  -0.7081  -0.1993   0.5592   2.5399  
## 
## Coefficients:
##                      Estimate  Std. Error z value             Pr(>|z|)    
## (Intercept)        8.47508193  0.25571653  33.142 < 0.0000000000000002 ***
## normalized.losses  0.00086989  0.00060863   1.429              0.15293    
## compression.ratio  0.02806643  0.00527721   5.318          0.000000105 ***
## horsepower         0.01152785  0.00048597  23.721 < 0.0000000000000002 ***
## peak.rpm          -0.00013301  0.00004436  -2.998              0.00272 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(14.1896) family taken to be 1)
## 
##     Null deviance: 791.77  on 204  degrees of freedom
## Residual deviance: 207.32  on 200  degrees of freedom
## AIC: 3887.7
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  14.19 
##           Std. Err.:  1.39 
## 
##  2 x log-likelihood:  -3875.675
# modelo sugerido 2 
yu2<-glm(price~curb.weight+length , family = Gamma(link=log),data = t2)

summary(yu2)
## 
## Call:
## glm(formula = price ~ curb.weight + length, family = Gamma(link = log), 
##     data = t2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.76966  -0.15834  -0.04829   0.06193   1.05941  
## 
## Coefficients:
##                Estimate  Std. Error t value            Pr(>|t|)    
## (Intercept)  7.50654984  0.40039726  18.748 <0.0000000000000002 ***
## curb.weight  0.00092313  0.00007622  12.112 <0.0000000000000002 ***
## length      -0.00280031  0.00321675  -0.871               0.385    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.07376687)
## 
##     Null deviance: 55.861  on 204  degrees of freedom
## Residual deviance: 12.437  on 202  degrees of freedom
## AIC: 3850
## 
## Number of Fisher Scoring iterations: 4

Tabela de modelos

library(hnp)

ajuste = c('mod', 'mod2','mod3','yu','yu2','mod.nb')
aic    = c(AIC(mod),AIC(mod2),AIC(mod3),AIC(yu),AIC(yu2),AIC(mod.nb))
deviance    = c(deviance(mod),deviance(mod2),deviance(mod3),deviance(yu),deviance(yu2),deviance(mod.nb))
#verossimilhanca = c(logLik(glmod1),logLik(glmod2),logLik(glmod3))
df<-data.frame(ajuste, round(aic,2),round(deviance,2)) 
colnames(df)=c("ajuste","AIC","Deveiance")

df%>% datatable()

Como podemos ver pelo O critério de informação de Akaike (AIC) e o deviance, o modelo que apresentou os melhores resultados segundo os critérios foi o modelo mod, que o modelo que contem variáveis removidas anteriormente,entretanto ao verificar o gráfico de envelope simulado normalmente distribuído para os resíduos, notamos que temos vários pontos fora do envelope, logo mesmo o critério de Akaike (AIC) e o deviance sendo o mais baixo, o modelo aparenta não ser o mais indicado

set.seed(1234)
par(bg = '#586573')
hnp(mod$residuals, sim = 99,resid.type ='deviance',how.many.out=T ,conf = 0.95,scale = T)
## Half-normal plot with simulated envelope generated assuming the residuals are 
##         normally distributed under the null hypothesis. 
## Estimated mean: -0.00000000002031851 
## Estimated variance: 0.03743822

## Total points: 205 
## Points out of envelope: 55 ( 26.83 %)

Como podemos observar de fato temos vários pontos fora do envelope, com isso o vamos observar os outros modelos vamos fazer a verificação do envelope para o modelo 2

Verificando o envelope do modelo 3

set.seed(1234)
par(bg = '#586573')
hnp(mod3$residuals, sim = 99,resid.type ='deviance',how.many.out=T ,conf = 0.95,scale = T)
## Half-normal plot with simulated envelope generated assuming the residuals are 
##         normally distributed under the null hypothesis. 
## Estimated mean: -0.0000000002925218 
## Estimated variance: 0.06734035

## Total points: 205 
## Points out of envelope: 1 ( 0.49 %)
# mod<-glm(test$price~.,family = Gamma(link=log),data = test)
 
 summary(mod2)
## 
## Call:
## glm(formula = price ~ normalized.losses + compression.ratio + 
##     horsepower + peak.rpm, family = Gamma(link = log), data = novo)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.60361  -0.18808  -0.05306   0.14854   0.67458  
## 
## Coefficients:
##                      Estimate  Std. Error t value             Pr(>|t|)    
## (Intercept)        8.47472093  0.25050641  33.830 < 0.0000000000000002 ***
## normalized.losses  0.00086997  0.00059626   1.459              0.14612    
## compression.ratio  0.02806842  0.00517010   5.429          0.000000163 ***
## horsepower         0.01152907  0.00047617  24.212 < 0.0000000000000002 ***
## peak.rpm          -0.00013297  0.00004346  -3.060              0.00252 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.06771908)
## 
##     Null deviance: 55.861  on 204  degrees of freedom
## Residual deviance: 14.625  on 200  degrees of freedom
## AIC: 3887.6
## 
## Number of Fisher Scoring iterations: 6
library('hnp')
set.seed(1234)
par(bg = '#586573')
hnp(mod2$residuals, sim = 99,resid.type ='deviance',how.many.out=T ,
    conf = 0.95,scale = T)
## Half-normal plot with simulated envelope generated assuming the residuals are 
##         normally distributed under the null hypothesis. 
## Estimated mean: -0.0000000002211503 
## Estimated variance: 0.06639125

## Total points: 205 
## Points out of envelope: 0 ( 0 %)
# hnp(yu, sim = 99,resid.type ='deviance',how.many.out=T ,
#     conf = 0.95,scale = T)

podemos ver que os resídeuos do segundo modelo ficou bem envelopado,sendo um bom ajuste caso não tenhamos outra escolha.Como nós já ajustamos o modelo binomial negativo vamos comparar seu gráfico de envelope simulado para os resíduos:

set.seed(1234)
par(bg = '#586573')
hnp(mod.nb$residuals, sim = 99,resid.type ='deviance',how.many.out=T ,
    conf = 0.95,scale = T)
## Half-normal plot with simulated envelope generated assuming the residuals are 
##         normally distributed under the null hypothesis. 
## Estimated mean: -0.00001070088 
## Estimated variance: 0.06639428

## Total points: 205 
## Points out of envelope: 0 ( 0 %)

Como podemos ver o modelo binomial negativo também aparenta um bom ajuste pelo gráfico de envelope ,logo, vamos comparar o AIC e o deviance desses modelos para decidir qual será o modelo final

ajuste = c('mod2','mod3','mod.nb')
aic    = c(AIC(mod2),AIC(mod3),AIC(mod.nb))
deviance    = c(deviance(mod2),deviance(mod3),deviance(mod.nb))
#verossimilhanca = c(logLik(glmod1),logLik(glmod2),logLik(glmod3))
df<-data.frame(ajuste, round(aic,2),round(deviance,2)) 
colnames(df)=c("ajuste","AIC","Deveiance")

df%>% datatable()

Podemos observar que o modelo mod2 aparenta ser o melhor deviance,vamos observar o Rsquad

Comparando o R2

#install.packages('DescTools')
library(DescTools)
PseudoR2(mod2)
PseudoR2(mod3)
PseudoR2(mod.nb)


ajuste = c('mod2','mod3','mod.nb')
R2    = c(PseudoR2(mod2),
PseudoR2(mod3),
PseudoR2(mod.nb))

#verossimilhanca = c(logLik(glmod1),logLik(glmod2),logLik(glmod3))
df<-data.frame(ajuste, round(R2,4)) 
colnames(df)=c("modelos","R2")

fazendo a comaparação dos R2 percebemos que o único modelo capaz de explicar ao menos 70% da variabilidade presente em nossa variável de interesse é o modelo binomial negtaivo.

Conclusão

Como podemos ver mesmo o modelo mod sendo o melhor qualificado nas metricas do critério de informação de Akaique e na Deviance o grafico de envelope dos resíduos apresentou varios pontos fora do envelope, logo, não é o modelo mais indicado para explicar nossa variável de interesse,sendo assim escolhemos o modelo binomial negativo embora não tenha apresentado o melhor AIC nem a melhor deviance pussui o melhor R2, sendo assim o melhor modelo para explicar o preço dos carros.

Exibindo os Coeficientes do Nosso Modelo

mod.nb$coefficients
##       (Intercept) normalized.losses compression.ratio        horsepower 
##      8.4750819321      0.0008698883      0.0280664268      0.0115278524 
##          peak.rpm 
##     -0.0001330111

Referências

O primeiro automóvel-INFO.Escola , acessado em 10/10/2022 https://www.infoescola.com/curiosidades/historia-do-automovel/

1985 Modelo de Importação de Carro e Especificações de Caminhão, Anuário Automotivo de 1985 Ward.

Personal Auto Manuals, Insurance Services Office, 160 Water Street, Nova York, NY 10038

Insurance Collision Report, Insurance Institute for Highway Safety, Watergate 600, Washington, DC 20037

dados veículos- IBGE ,acessado em 09/10/2022 https://cidades.ibge.gov.br/brasil/pesquisa/22/28120