Tópicos 2 — Modelagem com apoio computacional | Prof. Dr. Helton Saulo Bezerra dos Santos
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, readxl, janitor,xgboost,tidymodels,vip,rpart.plot,
conflicted,tictoc,finetune,doParallel,DataExplorer,future,glmnet,
data.table,summarytools,knitr,compareGroups,bonsai,discrim,
baguette,naivebayes)
registerDoParallel(cores = parallel::detectCores())
tidymodels_prefer()
df = read_csv("data.csv")Modelagem estatística é uma ferramenta fundamental para compreensão da relação de variáveis, principalmente quando tratamos de conjuntos de dados grandes, com relações complexas e padrões ocultos quando analisados de forma univariado ou bivariados.
Ao lidar com dados de renda, é fundamental entender a relação de covariáveis que são pertinentes que expliquem esta variável, visto que sabemos que a renda média é influenciada por fatores socioculturais, educacionais e demográficos, por exemplo. Entender e traçar todas essas relações uma a uma pode ser uma tarefa complicada, ou talvez até impossível, visto que é possível que existam relações que só são explicadas se combinados dois ou mais fatores afim do clarividência do padrão.
O entendimento do funcionamento do modelo é fundamental para obter uma melhor qualidade de ajuste. Existem diversos trabalhos na área de estatística que estudam o funcionamento e melhor ajuste e calibração destes modelos, dos quais destacarei o trabalho de Pinheiro (2023) que foi uma das fontes de inspiração deste trabalho.
Neste trabalho, buscou-se estudar um conjunto censitário sintético com diversas características demográficas, com objetivo fim de entender a relação entre covariáveis o salário — variável esta que se encontra binarizada no conjunto de dados, sendo:
Os dados são públicos e podem ser acessados em Kaggle.
Apesar da natureza sintética dos dados, estes traduzem uma rotina cada vez mais comum do estatístico: lidar com conjuntos de dados massivos, complicados e de difícil separação.
A grande dificuldade para este conjunto de dados é justamente encontrar um modelo que performe bem ante a um conjunto de dados complexo. Temos categorias desbalanceadas, um alto número de observações e covariáveis relacionadas, fatores estes que podem levar a problemas no ajuste do modelo.
Existem inúmeras implementações de algoritmos que trabalham baseados na classificação binária de dados, sendo cada uma destas mais pertinente ou precisa para diferentes situações, em que se deve levar em consideração a natureza do dado com a qual se está trabalhando. No conjunto de dados utilizado no escopo deste trabalho, temos covariáveis numéricas e categóricas — sendo a última a responsável pela maioria das características existentes nos dados — o que possivelmente já descarta algumas possibilidades de algoritmos melhor preparados para lidar exclusivamente ou principalmente com dados numéricos, como k-means, análise de discriminante linear e quadrático, KNN, etc.
Além disso, está contido no escopo deste trabalho uma preocupação com a legibilidade e interpretação da aplicação dos modelos aos códigos. É fundamental no mercado de trabalho que o estatístico, ao trabalhar em grupo, produza análises e códigos que possam ser compreendidos e utilizados por outros técnicos da área, bem como garantir a fluidez e infiabilidade do código ante aos resultados publicados. Portanto, é importante utilizar um framework que seja simples, verborrágico e altamente interpretativo, trazendo ao técnico a usabilidade do produto, bem como a facilidade de interpretação para um terceiro que porventura necessite prestar suporte ao produto. Desta forma, optou-se por utilizar o tidymodels para a realização do trabalho, unificando a preparação do modelo afim de torná-los fidedignamente comparáveis ante aos seus resultados, bem como permitindo uma alta carga interpretativa sobre sua implementação.
Para este estudo, utilizou-se dados sintéticos disponíveis na plataforma Kaggle. Nesta plataforma, é comum os usuários postarem conjuntos de dados afim da comunidade trabalhar sobre eles possibilidades de modelagem, levando a competição entre usuários afim de encontrar o melhor ou mais eficiente modelo que descreva aqueles dados. Isso é interessante por um lado, pois permite a comparação dos resultados obtidos neste trabalho com o resultado do trabalho de outros usuários que tenham se debruçado sobre o mesmo conjunto e compartilhado seus resultados. Entretando, particularmente minha escolha sobre estes dados está mais relacionada ao fato da complexidade do conjunto ser proporcional ao que encontro no meu dia a dia. É muito comum, sobre um objetivo didático, deparar-se com conjuntos de dados pequenos, de fácil separação, ajuste e convergência quando em situação de aula na univerdade. Entretando, quando nos deparamos com dados reais, esta facilidade se esvai. Neste universo, nos vemos obrigados a aprender sobre pacotes e linguagens que trazem eficiencia e agilidade, bem como a modularização de nossas análises. Desta forma, tento trazer neste trabalho algumas preocupações relacionadas a esta rotina, como parcimônia na procura de hiperparâmetros de modelos, utilização de validação cruzada e paralelização do código sempre que possível.
Para realizar o ajuste deste modelo, selecionou-se três algoritmos dente os vários disponíveis para esta rotina. O algoritmo canônico, que quase sempre deve ser utilizado quando frente a um problema de classificação binária, é a regressão logística, visto que esta têm uma alta carga interpretativa sobre os parâmetros do modelo, sendo bastante flexível, além de computacionalmente eficiente.
Um possível problema da regressão logística é uma possível dificuldade na seleção de covariáveis pertinentes para a análise. É possível realizar procedimentos automáticos como o algoritmo stepwise para seleção destes parâmetros, mas podemos pensar também em técnicas que nativamente já são delineadas afim de considerar esta problemática. Uma dessas técnicas é a regressão Lasso, que realiza uma penalização em covariáveis pouco significativas afim de zerar seus coeficientes, tornando assim o modelo mais parcimonioso. Este procedimento aumenta o viés do modelo, porém diminui sua variância, estando portanto contido no paradigma do bias-variance tradeoff. Como neste exemplo estamos mais preocupados com a previsão do que com a interpretação dos parâmetros — ainda que seja interessante recuperar alguma carga interpretativa destes — podemos avaliar a eficácia desta penalização ante a comparação com a modelagem logística diretamente.
Será conduzido ainda um estudo de utilização do modelo XGBoosting, que é uma implementação eficiente de um modelo de gradient boosting afim de comparar com os dois modelos anteriores mais clássicos, visando uma melhora na qualidade da previsão do modelo. Por ser uma metodologia que utiliza da computação intensiva e menor carga interpretativa sobre os parâmetros do modelo, é esperado deste modelo que seja significativamente superior aos anteriores, visando a justificativa da sua utilização.
Do conjunto de dados recuperados do kaggle, realizou-se algumas rotinas de engenharia de features para melhor acomodação das covariáveis ao formato que se intendia realizar a modelagem, melhorando seu poder preditivo e também adaptando as necessidades as quais este trabalho visa responder. Existiam nos dados algumas covariáveis categóricas com diversos fatores possíveis, mas que para o escopo deste trabalho não seriam interessantes, portanto foram agrupadas segundo categorias mais convenientes. Os agrupamentos realizados podem ser observados abaixo.
# Transformando em fator a variável resposta, afim de evitar erros na modelagem
df$salary = factor(df$salary)
# tornando o nome das colunas mais trabalháveis
df = clean_names(df)
# Ajustando a covariável de peso
df$fnlwgt = df$fnlwgt/sum(df$fnlwgt)
# Realizando agrupamentos
df = df %>%
mutate(marital_status = factor(marital_status),
relationship = factor(relationship),
race = factor(race),
sex = factor(sex),
country = factor(ifelse(country == 'United-States',country,'Others')),
education = case_when(
education %in% c("Bachelors","Masters","Assoc-acdm",
"Assoc-voc","Doctorate","Prof-school") ~ "high_education",
.default = "low_education"),
education = factor(education),
occupation = case_when(
occupation %in% c("Adm-clerical", "Exec-managerial", "Prof-specialty",
"Tech-support", "Sales") ~ "white_collar",
# Ocupações que geralmente envolvem trabalho em escritório ou administrativo.
occupation %in% c("Craft-repair", "Farming-fishing", "Handlers-cleaners",
"Machine-op-inspct", "Transport-moving") ~ "blue_collar",
# Ocupações que envolvem trabalho manual ou técnico.
occupation %in% c("Other-service", "Priv-house-serv",
"Protective-serv") ~ "services",
# Ocupações no setor de serviços.
occupation %in% c("Armed-Forces","?") ~ "Others"),
occupation = factor(occupation),
civil_status = ifelse(marital_status %in% c('Never-married','Divorced',
'Separated', 'Widowed') |
relationship %in% c('Not-in-family','Unmarried'),
"Single","Non-single"),
civil_status = factor(civil_status),
workclass = case_when(
workclass %in% c("Federal-gov", "Local-gov", "State-gov") ~ "government",
workclass == "Private" ~ "private",
workclass %in% c("Self-emp-inc", "Self-emp-not-inc") ~ "self_employed",
workclass %in% c("Never-worked", "Without-pay") ~ "not_working",
workclass == "?" ~ "unknown"),
workclass = factor(workclass)
) %>%
select(-marital_status,-relationship)Primeiramente, utilizou-se do pacote compareGroups para elaborar uma tabela simples porém eficiente e didática para visualização das diferenças entre as faixas salariais discriminadas por covariáveis, afim de uma primeira análise exploratória sobre os dados.
data = df %>% mutate(salary = factor(ifelse(salary == 0,"<50k",">50k")))
comp = compareGroups(salary ~ . -fnlwgt, data=data, method = 4)
table = createTable(comp)export2md(table,
strip=TRUE,
first.strip=TRUE,
format='html',
size = 10,
caption = "")| <50k | >50k | p.overall | |
|---|---|---|---|
| N=24720 | N=7841 | ||
| age | 34.0 [25.0;46.0] | 44.0 [36.0;51.0] | 0.000 |
| workclass: | <0.001 | ||
| government | 3010 (12.2%) | 1341 (17.1%) | |
| not_working | 21 (0.08%) | 0 (0.00%) | |
| private | 17733 (71.7%) | 4963 (63.3%) | |
| self_employed | 2311 (9.35%) | 1346 (17.2%) | |
| unknown | 1645 (6.65%) | 191 (2.44%) | |
| education: | 0.000 | ||
| high_education | 5981 (24.2%) | 4535 (57.8%) | |
| low_education | 18739 (75.8%) | 3306 (42.2%) | |
| education_num | 9.00 [9.00;10.0] | 12.0 [10.0;13.0] | 0.000 |
| occupation: | 0.000 | ||
| blue_collar | 8362 (33.8%) | 1700 (21.7%) | |
| Others | 1660 (6.72%) | 192 (2.45%) | |
| services | 3744 (15.1%) | 349 (4.45%) | |
| white_collar | 10954 (44.3%) | 5600 (71.4%) | |
| race: | <0.001 | ||
| Amer-Indian-Eskimo | 275 (1.11%) | 36 (0.46%) | |
| Asian-Pac-Islander | 763 (3.09%) | 276 (3.52%) | |
| Black | 2737 (11.1%) | 387 (4.94%) | |
| Other | 246 (1.00%) | 25 (0.32%) | |
| White | 20699 (83.7%) | 7117 (90.8%) | |
| sex: | 0.000 | ||
| Female | 9592 (38.8%) | 1179 (15.0%) | |
| Male | 15128 (61.2%) | 6662 (85.0%) | |
| capital_gain | 0.00 [0.00;0.00] | 0.00 [0.00;0.00] | 0.000 |
| capital_loss | 0.00 [0.00;0.00] | 0.00 [0.00;0.00] | <0.001 |
| hours_per_week | 40.0 [35.0;40.0] | 40.0 [40.0;50.0] | 0.000 |
| country: | <0.001 | ||
| Others | 2721 (11.0%) | 670 (8.54%) | |
| United-States | 21999 (89.0%) | 7171 (91.5%) | |
| civil_status: | 0.000 | ||
| Non-single | 8357 (33.8%) | 6702 (85.5%) | |
| Single | 16363 (66.2%) | 1139 (14.5%) |
rm(data)Podemos observar já nesta tabela algumas diferenças entre os grupos de cada faixa salarial, que serão confirmadas ou não na etapa de modelagem. Notamos que para todas as covariáveis, os testes não paramétricos de diferença entre os grupos acusaram significância, ainda que esta possivelmente é inflacinada pelo tamanho do conjunto de dados, e o desbalanceamento de observações nestes.
df %>%
select(age, education_num, capital_gain, capital_loss, hours_per_week) %>%
gather() %>%
ggplot(aes(value)) +
geom_histogram(bins = 30) +
facet_wrap(~key, scales = 'free_x') +
labs(x = '', y = 'Frequência', title = 'Histogramas das Covariáveis Numéricas') +
theme_minimal()
kable(summary(select(df, age, education_num, capital_gain, capital_loss, hours_per_week)))| age | education_num | capital_gain | capital_loss | hours_per_week | |
|---|---|---|---|---|---|
| Min. :17.00 | Min. : 1.00 | Min. : 0 | Min. : 0.0 | Min. : 1.00 | |
| 1st Qu.:28.00 | 1st Qu.: 9.00 | 1st Qu.: 0 | 1st Qu.: 0.0 | 1st Qu.:40.00 | |
| Median :37.00 | Median :10.00 | Median : 0 | Median : 0.0 | Median :40.00 | |
| Mean :38.58 | Mean :10.08 | Mean : 1078 | Mean : 87.3 | Mean :40.44 | |
| 3rd Qu.:48.00 | 3rd Qu.:12.00 | 3rd Qu.: 0 | 3rd Qu.: 0.0 | 3rd Qu.:45.00 | |
| Max. :90.00 | Max. :16.00 | Max. :99999 | Max. :4356.0 | Max. :99.00 |
df %>%
select(salary, age, education_num, capital_gain, capital_loss, hours_per_week) %>%
gather(key = "variable", value = "value", -salary) %>%
ggplot(aes(x = factor(salary), y = value)) +
geom_boxplot() +
facet_wrap(~variable, scales = 'free_y') +
labs(x = 'Salary (0: <50k, 1: >50k)', y = 'Valor',
title = 'Boxplot das Variáveis Numéricas por Faixa de Salário') +
theme_minimal()
df %>%
select(salary, workclass, education, civil_status,
occupation, sex, race, country) %>%
gather(key = "variable", value = "value", -salary) %>%
mutate(variable = factor(variable,
levels = c("workclass", "education","civil_status",
"occupation", "sex", "race", "country"))) %>%
ggplot(aes(x = value, fill = factor(salary))) +
geom_bar(position = "dodge") +
facet_wrap(~variable, scales = 'free_x') +
labs(x = '', y = '', fill = 'Salário (0: <50k, 1: >50k)', title = '') +
scale_x_discrete(guide = guide_axis(n.dodge=3)) +
theme_minimal()
Podemos observar que para o grupo com renda anual acima de 50 mil dólares, categorias como a ocupação white collar são significativamente maiores que para o grupo de renda inferior a 50 mil dólares anuais, bem como negros representando mais que o dobro proporcionalmente do grupo de renda inferior a 50 mil dólares anuais ante ao grupo de renda superior a isso. Podemos observar ainda uma quantidade maior de anos de estudo, e portanto um nível de escolaridade significativamente superior no grupo de renda superior a 50 mil dólares. Destarte, é interessante observar ainda uma diferença possivelmente significativa de composição do grupo que trabalha na área privada ante as faixas salariais, sendo os trabalhadores desta categoria pior remunerados num geral. Aparenta haver ainda uma diferença significativa na composição de gênero entre as faixas salariais, com mulheres representando quase metade do grupo com renda inferior a 50 mil dólares anuais, e sendo menos da metade deste percentual no grupo com renda superior a 50 mil dólares anuais. Nesta tabela, a maior diferença significativa aparenta ser o estado civil, com não solteiros compondo a esmagadora maioria do grupo com renda superior a 50 mil dólares anuais, enquanto para o grupo com renda inferior a isto, os não solteiros também são maioria, porém uma maioria muito menos significativa, de apenas 66% ante a 85%.
Num contexto em que temos tantas covariáveis, tantas observações e, apesar de alguns indicativos observados na análise exploratória, não é possível observar um padrão óbvio que indique o salário do indivíduo. Para isso, utilizei do recurso da modelagem com apoio computacional afim de tornar possível esta análise. Um diferencial deste trabalho é a utilização do framework tidymodels, que é bastante verborrágico e permite uma compreensão das etapas do modelo pela leitura do código, além de eficiência e praticidade.
Irei testar diversos modelos e fazer comparação de resultados destes, mas utilizarei a mesma receita para todos.
set.seed(150167636)
split = initial_split(df, prop = .8, strata = salary)
train = training(split)
test = testing(split)
cv_folds <- vfold_cv(train,
v = 5,
strata = salary)
recipe <- recipe(salary ~ .,
data = train) %>%
update_role(fnlwgt, new_role = "case_weight") %>%
step_normalize(all_numeric_predictors()) %>%
step_dummy(all_nominal_predictors()) %>%
step_interact(terms = ~ starts_with("occupation"):starts_with("race") +
starts_with("occupation"):starts_with("sex") +
starts_with("hours_per_week"):starts_with("sex"))
dados_preparados <- recipe %>%
prep() %>%
juice()
kable(head(dados_preparados))| age | fnlwgt | education_num | capital_gain | capital_loss | hours_per_week | salary | workclass_not_working | workclass_private | workclass_self_employed | workclass_unknown | education_low_education | occupation_Others | occupation_services | occupation_white_collar | race_Asian.Pac.Islander | race_Black | race_Other | race_White | sex_Male | country_United.States | civil_status_Single | occupation_Others_x_race_Asian.Pac.Islander | occupation_Others_x_race_Black | occupation_Others_x_race_Other | occupation_Others_x_race_White | occupation_services_x_race_Asian.Pac.Islander | occupation_services_x_race_Black | occupation_services_x_race_Other | occupation_services_x_race_White | occupation_white_collar_x_race_Asian.Pac.Islander | occupation_white_collar_x_race_Black | occupation_white_collar_x_race_Other | occupation_white_collar_x_race_White | occupation_Others_x_sex_Male | occupation_services_x_sex_Male | occupation_white_collar_x_sex_Male | sex_Male_x_hours_per_week |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.0314592 | 1.25e-05 | 1.1354862 | 0.1510791 | -0.2161513 | -0.0353722 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | -0.0353722 |
| 0.8376443 | 1.35e-05 | 1.1354862 | -0.1461666 | -0.2161513 | -2.2198774 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | -2.2198774 |
| -0.0418303 | 3.49e-05 | -0.4209394 | -0.1461666 | -0.2161513 | -0.0353722 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -0.0353722 |
| -0.7747258 | 5.48e-05 | 1.1354862 | -0.1461666 | -0.2161513 | -0.0353722 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0.0000000 |
| 0.7643547 | 2.59e-05 | -1.9773649 | -0.1461666 | -0.2161513 | -1.9771546 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.0000000 |
| -1.1411735 | 1.98e-05 | 1.1354862 | -0.1461666 | -0.2161513 | -0.8444482 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0.0000000 |
Podemos ver um pedaço da forma geral dos dados preparados para serem inseridos no modelo.
Como a variável resposta é binária, o primeiro modelo que podemos tentar seria o logístico
glm_spec <- logistic_reg() %>%
set_engine("glm")
glm_wf <- workflow() %>%
add_recipe(recipe) %>%
add_model(glm_spec)
glm_fit <- glm_wf %>%
fit(data = train)
glm_fit %>%
extract_fit_parsnip() %>%
tidy() %>%
kable()| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -1.7372041 | 0.4866916 | -3.5694143 | 0.0003578 |
| age | 0.4055276 | 0.0227845 | 17.7984123 | 0.0000000 |
| education_num | 0.8441476 | 0.0438024 | 19.2717200 | 0.0000000 |
| capital_gain | 2.2344546 | 0.0800493 | 27.9134654 | 0.0000000 |
| capital_loss | 0.2691885 | 0.0165243 | 16.2904180 | 0.0000000 |
| hours_per_week | 0.3348932 | 0.0455127 | 7.3582312 | 0.0000000 |
| workclass_not_working | -11.7789355 | 117.4336842 | -0.1003029 | 0.9201039 |
| workclass_private | -0.0699743 | 0.0559738 | -1.2501244 | 0.2112541 |
| workclass_self_employed | -0.3486299 | 0.0730037 | -4.7755074 | 0.0000018 |
| workclass_unknown | -0.5554042 | 1.8010975 | -0.3083699 | 0.7578009 |
| education_low_education | 0.1502829 | 0.0759029 | 1.9799356 | 0.0477108 |
| occupation_Others | 0.1946984 | 2.0025527 | 0.0972251 | 0.9225476 |
| occupation_services | 0.8077583 | 0.8054694 | 1.0028417 | 0.3159372 |
| occupation_white_collar | 0.4311601 | 0.5736424 | 0.7516183 | 0.4522806 |
| race_Asian.Pac.Islander | 0.7500891 | 0.5035234 | 1.4896806 | 0.1363082 |
| race_Black | 0.4686787 | 0.4650463 | 1.0078109 | 0.3135452 |
| race_Other | 0.1107415 | 0.6625542 | 0.1671433 | 0.8672573 |
| race_White | 0.5799804 | 0.4454281 | 1.3020741 | 0.1928910 |
| sex_Male | 0.1722064 | 0.1804855 | 0.9541286 | 0.3400185 |
| country_United.States | 0.2594948 | 0.0762643 | 3.4025743 | 0.0006675 |
| civil_status_Single | -2.4131329 | 0.0517315 | -46.6472524 | 0.0000000 |
| occupation_Others_x_race_Asian.Pac.Islander | -1.1841893 | 1.4582220 | -0.8120775 | 0.4167472 |
| occupation_Others_x_race_Black | -0.4225307 | 1.3041276 | -0.3239949 | 0.7459419 |
| occupation_Others_x_race_Other | 1.4529448 | 1.6820682 | 0.8637847 | 0.3877062 |
| occupation_Others_x_race_White | 0.0266894 | 1.2087874 | 0.0220795 | 0.9823846 |
| occupation_services_x_race_Asian.Pac.Islander | -1.2872096 | 0.8711872 | -1.4775350 | 0.1395323 |
| occupation_services_x_race_Black | -1.1793043 | 0.8168805 | -1.4436680 | 0.1488324 |
| occupation_services_x_race_Other | -1.0749159 | 1.4074752 | -0.7637192 | 0.4450346 |
| occupation_services_x_race_White | -1.1017327 | 0.7831182 | -1.4068537 | 0.1594708 |
| occupation_white_collar_x_race_Asian.Pac.Islander | -0.3004502 | 0.6078681 | -0.4942687 | 0.6211164 |
| occupation_white_collar_x_race_Black | 0.0860022 | 0.5730514 | 0.1500777 | 0.8807033 |
| occupation_white_collar_x_race_Other | -0.3731938 | 0.8639278 | -0.4319734 | 0.6657608 |
| occupation_white_collar_x_race_White | 0.1513841 | 0.5481445 | 0.2761755 | 0.7824132 |
| occupation_Others_x_sex_Male | -0.3821967 | 0.3098797 | -1.2333712 | 0.2174373 |
| occupation_services_x_sex_Male | 0.2564560 | 0.2581336 | 0.9935012 | 0.3204658 |
| occupation_white_collar_x_sex_Male | 0.1581810 | 0.1879843 | 0.8414584 | 0.4000912 |
| sex_Male_x_hours_per_week | 0.0586262 | 0.0513488 | 1.1417232 | 0.2535691 |
Podemos observar as covariáveis de maior importância para este modelo, assim como sua matriz de confusão
glm_fit %>%
extract_fit_parsnip() %>%
vip(geom = "col") +
theme_minimal() +
theme(axis.text.x = element_text(size = 34),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
plot.title = element_text(size = 18, hjust = 0.5))
predictions <- glm_fit %>%
predict(new_data = test)
results <- bind_cols(test, predictions)
results %>%
conf_mat(salary,.pred_class) Truth
Prediction 0 1
0 4556 649
1 388 920metrics <- metric_set(accuracy, specificity, sensitivity)
resultados_logit = augment(glm_fit, new_data = test) %>%
metrics(truth = salary, estimate = .pred_class) %>%
select(.metric,.estimate)Como observado na receita do modelo, existem 38 covariáveis nesta modelagem. Diversas abordagens podem ser utilizadas para selecionar as covariáveis de maior importância, sendo uma dessas a regressão lasso, que penaliza coeficientes e torna-os 0 em caso de insignificância.
Este é um modelo que contém um hiperparâmetro, portanto iremos ajustar um grid para escolher o melhor possível.
lasso_spec <- logistic_reg(penalty = tune(),
mixture = 1) %>%
set_engine("glmnet")
lasso_wf <- workflow() %>%
add_recipe(recipe) %>%
add_model(lasso_spec)
grid <- grid_regular(penalty(),
levels = 100)
plan(multisession)
set.seed(150167636)
lasso_res <- lasso_wf %>%
tune_grid(resamples = cv_folds,
grid = grid,
metrics = metric_set(roc_auc))lasso_res %>%
collect_metrics() %>%
filter(.metric == "roc_auc") %>%
select(mean, penalty) %>%
pivot_longer(penalty,
names_to = "hiperparâmetro",
values_to = "valor") %>%
ggplot(aes(valor, mean)) +
geom_point()+
theme_minimal() +
theme(axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10),
axis.title.x = element_text(size = 10),
axis.title.y = element_text(size = 10),
plot.title = element_text(size = 10, hjust = 0.5))+
labs(x="",y="")
Podemos ver que um menor valor de penalização é benéfico ao modelo, visto a importância relativa das covariáveis serem altas neste caso
best_params = lasso_res %>%
select_best(metric = "roc_auc")
best_wf = finalize_workflow(lasso_wf, best_params)
final_fit <- best_wf %>%
fit(data = train)
final_fit %>%
extract_fit_parsnip() %>%
tidy() %>%
kable()| term | estimate | penalty |
|---|---|---|
| (Intercept) | -1.1948793 | 0.0007391 |
| age | 0.3943655 | 0.0007391 |
| education_num | 0.7870445 | 0.0007391 |
| capital_gain | 2.1066124 | 0.0007391 |
| capital_loss | 0.2622005 | 0.0007391 |
| hours_per_week | 0.3315385 | 0.0007391 |
| workclass_not_working | -1.9072610 | 0.0007391 |
| workclass_private | -0.0252262 | 0.0007391 |
| workclass_self_employed | -0.2793742 | 0.0007391 |
| workclass_unknown | -0.1482158 | 0.0007391 |
| education_low_education | 0.0470209 | 0.0007391 |
| occupation_Others | -0.1167739 | 0.0007391 |
| occupation_services | -0.0488290 | 0.0007391 |
| occupation_white_collar | 0.3271691 | 0.0007391 |
| race_Asian.Pac.Islander | 0.0246369 | 0.0007391 |
| race_Black | 0.0000000 | 0.0007391 |
| race_Other | -0.2075199 | 0.0007391 |
| race_White | 0.0681928 | 0.0007391 |
| sex_Male | 0.1888269 | 0.0007391 |
| country_United.States | 0.2273567 | 0.0007391 |
| civil_status_Single | -2.3832942 | 0.0007391 |
| occupation_Others_x_race_Asian.Pac.Islander | -0.5597923 | 0.0007391 |
| occupation_Others_x_race_Black | -0.1860392 | 0.0007391 |
| occupation_Others_x_race_Other | 0.5816320 | 0.0007391 |
| occupation_Others_x_race_White | 0.0000000 | 0.0007391 |
| occupation_services_x_race_Asian.Pac.Islander | 0.0000000 | 0.0007391 |
| occupation_services_x_race_Black | -0.0567795 | 0.0007391 |
| occupation_services_x_race_Other | 0.0000000 | 0.0007391 |
| occupation_services_x_race_White | 0.0000000 | 0.0007391 |
| occupation_white_collar_x_race_Asian.Pac.Islander | 0.0000000 | 0.0007391 |
| occupation_white_collar_x_race_Black | 0.1298919 | 0.0007391 |
| occupation_white_collar_x_race_Other | -0.2649066 | 0.0007391 |
| occupation_white_collar_x_race_White | 0.2695555 | 0.0007391 |
| occupation_Others_x_sex_Male | -0.3617209 | 0.0007391 |
| occupation_services_x_sex_Male | 0.0000000 | 0.0007391 |
| occupation_white_collar_x_sex_Male | 0.1420720 | 0.0007391 |
| sex_Male_x_hours_per_week | 0.0511361 | 0.0007391 |
Podemos observar as covariáveis de maior importância para este modelo, assim como sua matriz de confusão
final_fit %>%
extract_fit_parsnip() %>%
vip(geom = "col") +
theme_minimal() +
theme(axis.text.x = element_text(size = 34),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
plot.title = element_text(size = 18, hjust = 0.5))
predictions <- final_fit %>%
predict(new_data = test)
results <- bind_cols(test, predictions)
results %>%
conf_mat(salary,.pred_class) Truth
Prediction 0 1
0 4571 652
1 373 917metrics <- metric_set(accuracy, specificity, sensitivity)
resultados_lasso = augment(final_fit, new_data = test) %>%
metrics(truth = salary, estimate = .pred_class) %>%
select(.metric,.estimate)O XGBoost é um modelo de Gradient boosting baseado em árvores, que costuma performar bem em tarefas como esta, de classificação com diversas covariáveis
Este é um modelo de Boosting, ou seja, o “encaixe” de diversos modelos fracos afim da construção de um modelo robusto a partir da combinação destes resultados.
graph LR
A[Dados] --> B[Modelo 1];
B --> C[Modelo 2];
C --> D[...];
D --> E[Modelo m];
E --> F[Ensembling dos modelos];
F --> G[Modelo final];
Também iremos realizar o fine tuning de alguns hiperparâmetros deste modelo, no caso o número de árvores e a profundidade destas, afim de obter o melhor modelo.
model = boost_tree(mode = "classification",
trees = tune(),
tree_depth = tune()
) %>%
set_engine("xgboost")
wf = workflow() %>%
add_recipe(recipe) %>%
add_model(model)
grid = wf %>%
extract_parameter_set_dials() %>%
grid_regular(levels = 3)
plan(multisession)
set.seed(150167636)
tune_res = tune_grid(
wf,
resamples = cv_folds,
grid = grid,
metrics = metric_set(accuracy, roc_auc, sens,spec)
)tune_res %>%
collect_metrics() %>%
filter(.metric == "roc_auc") %>%
select(mean, trees:tree_depth) %>%
pivot_longer(trees:tree_depth,
names_to = "hiperparâmetro",
values_to = "valor") %>%
ggplot(aes(valor, mean, color = hiperparâmetro)) +
geom_point(show.legend = FALSE) +
facet_wrap(~hiperparâmetro, scales = "free_x") +
theme_minimal() +
theme(axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10),
axis.title.x = element_text(size = 10),
axis.title.y = element_text(size = 10),
plot.title = element_text(size = 10, hjust = 0.5))+
labs(x="",y="")
best_params = tune_res %>%
select_best(metric = "roc_auc")
best_params# A tibble: 1 × 3
trees tree_depth .config
<int> <int> <chr>
1 2000 1 Preprocessor1_Model3Vemos que a melhor combinação de hiperparâmetros encontrada é utilizando 2000 árvores de tamanho 1
Podemos observar as covariáveis de maior importância para este modelo, assim como sua matriz de confusão
best_wf = finalize_workflow(wf, best_params)
final_fit <- best_wf %>%
fit(data = train)final_fit %>%
extract_fit_parsnip() %>%
vip(geom = "col") +
theme_minimal() +
theme(axis.text.x = element_text(size = 34),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
plot.title = element_text(size = 18, hjust = 0.5))
predictions <- final_fit %>%
predict(new_data = test)
results <- bind_cols(test, predictions)
results %>%
conf_mat(salary,.pred_class) Truth
Prediction 0 1
0 4636 584
1 308 985metrics <- metric_set(accuracy, specificity, sensitivity)
resultados_xgboost = augment(final_fit, new_data = test) %>%
metrics(truth = salary, estimate = .pred_class) %>%
select(.metric,.estimate)Podemos verificar lado a lado a performance de cada modelo sobre este conjunto de dados.
Modelo logístico
kable(resultados_logit)| .metric | .estimate |
|---|---|
| accuracy | 0.8407800 |
| specificity | 0.5863607 |
| sensitivity | 0.9215210 |
Lasso
kable(resultados_lasso)| .metric | .estimate |
|---|---|
| accuracy | 0.8426224 |
| specificity | 0.5844487 |
| sensitivity | 0.9245550 |
XGBoost
kable(resultados_xgboost)| .metric | .estimate |
|---|---|
| accuracy | 0.8630431 |
| specificity | 0.6277884 |
| sensitivity | 0.9377023 |
O maior desafio para estes dados era capturar a especificidade (salário >50k). Vemos que todos os modelos tiveram dificuldade com esta métrica, porém houve um ganho sensível do modelo de árvore em relação aos modelos baseados em regressão. Entretando, a complexidade deste modelo possivelmente não justifica a sua utilização, visto um ganho tão baixo nas métricas ante aos modelos mais parcimoniosos.