HLM
Marcos H Silva
24/04/2022
Modelo de regressão multinível
Uma rápida introdução
O texto de hoje é sobre uma outra forma de criar um modelo capaz de “prever” um resultado baseado em dados. Meus últimos dois textos abordaram formas de criar redes neurais e como usá-las para prever resultados: Rede neural recorrente para identificação de dígitos e Rede LSTM para predição da cotação do dólar.
Hoje trago uma outra abordagem, tão importante quanto as redes neurais e capaz de criar resultados mais claros!
Imagine que você gostaria de avaliar qual seria o resultado de sua empresa usando por exemplo a nota do NPS que suas lojas (ou loja) recebeu. Para lhe ajudar com este cálculo você possui uma base de dados com registros históricos, acumulando as vendas dos últimos 12 meses e das notas que as lojas registraram.
Se colocarmos estes dados em um gráfico do excel com receita no eixo X e NPS no eixo Y, teremos o seguinte resultado:
Cada ponto preto representa o cruzamento de resultado e NPS de uma loja e a linha azul é a reta de regressão, uma linha que tenta passar o mais próximo possível de todos os pontos, numa tentativa de resumir o gráfico a uma só informação.
Um modelo de regressão tenta estimar qual seria a reta para dados que não estão em sua base, usando apenas as informações existentes. Por exemplo, qual será a receita de uma loja se ela tiver uma nota NPS 87,8?
A fórmula base para explicar um modelo de regressão linear com base neste exemplo seria:
\[ Y = \gamma_0 + NPS_1 * \gamma_1 + \epsilon \]
Onde :
\(\gamma_0\)(gama 0) representa uma a média da variável resposta
\(\gamma_1\) (gama 1) representa a participação do NPS na inclinação da reta de regressão
\(\epsilon\) (epsilon) representa a distância entre o valor real da receita e o valor que foi previsto pelos termos \(\gamma_1\), este valor também é chamado de erro ou resíduo
O desafio neste tipo de modelo é encontrar os valores de \(\gamma_0\) e \(\gamma_1\)para que possamos substituí-los em nossa fórmula e encontrar o valor da receita.
E quando o mundo não é linear? E quase sempre não é!
A modelagem linear cumpre bem o seu papel em determinados cenários, mas em muitos casos ela sozinha não é suficiente para explicar um mundo tão diverso. E é ai que entra a modelagem multinível! Um modelo que prevê que o mundo é feito de grupos e que estes grupos exercem influência sobre os dados. Este modelo considera que os dados podem ser agrupados em até três níveis:
Onde:
Grupo 1 (t), podemos trabalhar com uma escala de tempo, por exemplo, se houver dados de minha loja segmentados por mês;
Grupo 2 (i), temos o indivíduo, no nosso exemplo a loja;
Grupo 3 (j) temos o agrupamento, neste caso a cidade.
Neste cenário podemos avaliar que as lojas podem ter comportamento diferente dependendo da cidade onde estão localizadas e que variáveis específicas destas cidades podem influenciar em seu resultado, além do NPS.
Neste caso trabalharemos com os grupos i e j apenas, respectivamente os níveis 1 (loja) e 2 (cidade).
Preparação dos pacotes
pacotes = c("plotly","tidyverse","reshape2","knitr","kableExtra","rgl","car", "nlme","lmtest","fastDummies","msm","lmeInfo","jtools", "data.table", "olsrr")
if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
instalador <- pacotes[!pacotes %in% installed.packages()]
for(i in 1:length(instalador)) {
install.packages(instalador, dependencies = T)
break()}
sapply(pacotes, require, character = T)
} else {
sapply(pacotes, require, character = T)
}## plotly tidyverse reshape2 knitr kableExtra rgl
## TRUE TRUE TRUE TRUE TRUE TRUE
## car nlme lmtest fastDummies msm lmeInfo
## TRUE TRUE TRUE TRUE TRUE TRUE
## jtools data.table olsrr
## TRUE TRUE TRUEA função stderr_nlme, de autoria dos professores Fávero e Rafael Sousa ajuda a identificar a significância estatística dos dados de cada nível e de certa forma a calcular o quanto do comportamento dos dados é explicado por cada nível.
stderr_nlme <- function(model){
if(base::class(model) != "lme"){
base::message("Use a lme object model from nlme package")
stop()}
resume <- base::summary(model)
if(base::length(base::names(model$groups))==1){
m.type <- "HLM2"
} else if(base::length(base::names(model$groups))==2){
m.type <- "HLM3"
}
if(m.type == "HLM2"){
vcov_matrix <- model$apVar
logs_sd_re <- base::attr(vcov_matrix,"Pars")
if(base::length(logs_sd_re)==2){
stderr_tau00 <- msm::deltamethod(~exp(x1)^2,logs_sd_re,vcov_matrix)
stderr_sigma <- msm::deltamethod(~exp(x2)^2,logs_sd_re,vcov_matrix)
results <- base::data.frame(`RE Components`=base::c("Var(v0j)","Var(e)"),
`Variance Estimatives`= base::c(base::exp(logs_sd_re)[[1]]^2,
base::exp(logs_sd_re[[2]])^2),
`Std Err.`=base::c(stderr_tau00,
stderr_sigma),
z=base::c(base::exp(logs_sd_re)[[1]]^2/stderr_tau00,
base::exp(logs_sd_re[[2]])^2/stderr_sigma),
`p-value`=base::round(stats::pnorm(q=base::c(base::exp(logs_sd_re)[[1]]^2/stderr_tau00,
base::exp(logs_sd_re[[2]])^2/stderr_sigma),
lower.tail=F)*2,3))
return(results)
}
else{
stderr_tau00 <- msm::deltamethod(~exp(x1)^2,logs_sd_re,vcov_matrix)
stderr_tau01 <- msm::deltamethod(~exp(x2)^2,logs_sd_re,vcov_matrix)
stderr_sigma <- msm::deltamethod(~exp(x4)^2,logs_sd_re,vcov_matrix)
results <- base::data.frame(Components=base::c("Var(v0j)","Var(v1j)","Var(e)"),
`Variance Estimatives`= base::c(base::exp(logs_sd_re)[[1]]^2,
base::exp(logs_sd_re[[2]])^2,
base::exp(logs_sd_re[[4]])^2),
`Std Err.`=base::c(stderr_tau00,
stderr_tau01,
stderr_sigma),
z=base::c(base::exp(logs_sd_re)[[1]]^2/stderr_tau00,
base::exp(logs_sd_re[[2]])^2/stderr_tau01,
base::exp(logs_sd_re[[4]])^2/stderr_sigma),
`p-value`=base::round(stats::pnorm(q=base::c(base::exp(logs_sd_re)[[1]]^2/stderr_tau00,
base::exp(logs_sd_re[[2]])^2/stderr_tau01,
base::exp(logs_sd_re[[4]])^2/stderr_sigma),
lower.tail=F)*2,3))
return(results)
}
}
if(m.type == "HLM3"){
vcov_matrix <- model$apVar
logs_sd_re <- base::attr(vcov_matrix,"Pars")
if(base::length(logs_sd_re) == 3){
stderr_tau_r000 <- msm::deltamethod(~exp(x1)^2,logs_sd_re,vcov_matrix)
stderr_tau_u000 <- msm::deltamethod(~exp(x2)^2,logs_sd_re,vcov_matrix)
stderr_sigma <- msm::deltamethod(~exp(x3)^2,logs_sd_re,vcov_matrix)
results <- base::data.frame(Components=base::c("Var(t00k)","Var(v0jk)","Var(e)"),
Estimatives=base::c(base::exp(logs_sd_re)[[2]]^2,
base::exp(logs_sd_re)[[1]]^2,
base::exp(logs_sd_re)[[3]]^2),
Std_Err=base::c(stderr_tau_u000,
stderr_tau_r000,
stderr_sigma),
z=base::c(base::exp(logs_sd_re)[[2]]^2/stderr_tau_u000,
base::exp(logs_sd_re)[[1]]^2/stderr_tau_r000,
base::exp(logs_sd_re)[[3]]^2/stderr_sigma),
`p-value`=base::round(stats::pnorm(q=base::c(base::exp(logs_sd_re)[[2]]^2/stderr_tau_u000,
base::exp(logs_sd_re)[[1]]^2/stderr_tau_r000,
base::exp(logs_sd_re)[[3]]^2/stderr_sigma),
lower.tail=F)*2,3))
return(results)
}
else{
stderr_tau_r000 <- msm::deltamethod(~exp(x1)^2,logs_sd_re,vcov_matrix)
stderr_tau_r100 <- msm::deltamethod(~exp(x2)^2,logs_sd_re,vcov_matrix)
stderr_tau_u000 <- msm::deltamethod(~exp(x4)^2,logs_sd_re,vcov_matrix)
stderr_tau_u100 <- msm::deltamethod(~exp(x5)^2,logs_sd_re,vcov_matrix)
stderr_sigma <- msm::deltamethod(~exp(x7)^2,logs_sd_re,vcov_matrix)
results <- base::data.frame(`RE_Components`=base::c("Var(t00k)","Var(t10k)",
"Var(v0jk)","Var(v1jk)",
"Var(e)"),
`Variance Estimatives`=base::c(base::exp(logs_sd_re)[[4]]^2,
base::exp(logs_sd_re)[[5]]^2,
base::exp(logs_sd_re)[[1]]^2,
base::exp(logs_sd_re)[[2]]^2,
base::exp(logs_sd_re)[[7]]^2),
`Std Err.`=base::c(stderr_tau_u000,
stderr_tau_u100,
stderr_tau_r000,
stderr_tau_r100,
stderr_sigma),
z=base::c(base::exp(logs_sd_re)[[4]]^2/stderr_tau_u000,
base::exp(logs_sd_re)[[5]]^2/stderr_tau_u100,
base::exp(logs_sd_re)[[1]]^2/stderr_tau_r000,
base::exp(logs_sd_re)[[2]]^2/stderr_tau_r100,
base::exp(logs_sd_re)[[7]]^2/stderr_sigma),
`p-value`=base::round(stats::pnorm(q=base::c(base::exp(logs_sd_re)[[4]]^2/stderr_tau_u000,
base::exp(logs_sd_re)[[5]]^2/stderr_tau_u100,
base::exp(logs_sd_re)[[1]]^2/stderr_tau_r000,
base::exp(logs_sd_re)[[2]]^2/stderr_tau_r100,
base::exp(logs_sd_re)[[7]]^2/stderr_sigma),
lower.tail=F)*2,3))
return(results)
}
}
}Importação dos dados e análise exploratória
vendas =
fread("dados_vendas_totais_2021.csv", dec =",", sep=";") %>%
select(cidade, loja, receita, nps, inflacao) %>%
mutate(loja = as.factor(loja), cidade = as.factor(cidade)) Para testar o modelo vamos utilizar uma segunda base de dados semelhante à nossa base de treino (vendas) com 195 registros
teste =
fread("dados_para_teste.csv", dec =",", sep=";") %>%
select(cidade, loja, receita, nps, inflacao) %>%
mutate(loja = as.factor(loja), cidade = as.factor(cidade)) A visualização das primeiras 5 linhas das cidades 1 e 2, foi aplicado a função scalepara evitar que os dados sejam expostos
rbind(
vendas %>%
mutate_if(is.numeric, scale) %>%
filter(cidade == 1) %>%
filter(row_number() <= 5),
vendas %>%
mutate_if(is.numeric, scale) %>%
filter(cidade == 2) %>%
filter(row_number() <= 5)) %>%
kable() %>%
kable_styling(bootstrap_options = "striped",
position = "center",
full_width = F,
font_size = 14)| cidade | loja | receita | nps | inflacao |
|---|---|---|---|---|
| 1 | 1 | -0.6807167 | -1.0082727 | -0.0143209 |
| 1 | 2 | 0.4985884 | 0.7390536 | -0.0143209 |
| 1 | 3 | 0.4941708 | 1.2199451 | -0.0143209 |
| 1 | 4 | -0.4715108 | -0.4736577 | -0.0143209 |
| 1 | 5 | 0.4693952 | 1.0168438 | -0.0143209 |
| 2 | 168 | -0.5215818 | 0.4140915 | -0.6226181 |
| 2 | 169 | -1.3179497 | -1.4740954 | -0.6226181 |
| 2 | 170 | -1.3552432 | -0.7874820 | -0.6226181 |
| 2 | 171 | -0.7161907 | 1.4381152 | -0.6226181 |
| 2 | 172 | -0.9799720 | -0.3603141 | -0.6226181 |
Ao visualizarmos a função de densidade de probabilidade da função dependente (receita) observa-se distribuições diferentes dependendo da cidade, mostrando que temos uma dispersão grande dependendo da cidade.
vendas %>%
group_by(cidade) %>%
mutate(linhas = 1:n()) %>%
mutate(x = unlist(density(receita, n = max(linhas))["x"]),
y = unlist(density(receita, n = max(linhas))["y"])) %>%
mutate_if(is.numeric, scale) %>%
ggplot() +
geom_area(aes(x = x, y = y, group = cidade, fill = cidade), color = "black", alpha = 0.4) +
geom_histogram(aes(x = receita,
y = ..density..,
fill = cidade),
color = "black",
position = 'identity',
alpha = 0.3) +
facet_wrap(~ cidade) +
scale_fill_viridis_d() +
scale_color_viridis_d() +
theme_classic()## `mutate_if()` ignored the following grouping variables:
## Column `cidade`
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Verificando a inflação por cidade vemos que há uma dispersão significativa
vendas %>%
ggplot() +
geom_bar(aes(x = cidade, y=inflacao, fill = cidade), stat = "identity") +
scale_y_continuous() +
scale_fill_viridis_d() +
theme_minimal()
Temos aqui os indícios necessário para iniciar a exploração de modelos de regressão!
Criação dos modelos
Para a modelagem multinível utiliza-se a técnica step-up onde o modelo é otimizado a cada execução até obter-se sua versão final. Para verificar o modelo vamos utilizar o BIC (bayesian information criterion) critério que avalia o quão perto o modelo está dados reais considerando também as variáveis explicativas.
O primeiro modelo a ser criado é o modelo nulo, que não considera nenhuma variável preditora, apenas a receita e o agrupamento por cidade. Sua fórmula pode é:
\[ receita_{ij} = \gamma_{00} + \nu_{0j} + \epsilon_{ij} \]
Onde:
\(\gamma_{00}\) representa a média geral da variável receita
\(\nu_{0j}\) representa a influência da cidade sobre a receita
\(\nu_{ij}\) representa a distância entre a receita calculada para aquela observação e o valor real
Observe que todos os dados são subscritos de com i e j onde i representa o dado no nível loja e j no nível cidade.
No R a execução desta fórmula é feita pelo comando lme do pacote nlme
#fixed é a fórmula base, de regressão linear
#random é a fórmula de agrupamento dos dados, neste caso cidade
hlm_nulo =
lme(fixed = receita ~1,
random = ~1 | cidade,
data = vendas,
method = "REML")
BIC(hlm_nulo)## [1] 14212.24O BIC de nosso modelo base é 14212.24, ele será utilizado mais adiante para compararmos o ganho dos próximos modelos. Quando menor ele for, melhor!
stderr_nlme(hlm_nulo)## RE.Components Variance.Estimatives Std.Err. z p.value
## 1 Var(v0j) 626248.2 223061.02 2.80752 0.005
## 2 Var(e) 254473.4 11949.33 21.29604 0.000Outra informação importante vem através do método stderr_nlme, nele podemos ver que o p-value de nossa variável receita (v0j) é estatisticamente significativa (menor que 0,05), indicando que o modelo multinível é aplicável à este caso. Se ela fosse superior a 0,05 o modelo de regressão linear simples seria o mais indicado.
Utilizaremos agora a base de testes para executar nosso modelo e comparar se os resultados previstos estão próximos do resultado real.
teste$hlm_nulo_fitted = predict(hlm_nulo, teste)
teste %>%
mutate_if(is.numeric, scale) %>%
ggplot() +
geom_smooth(aes(x = receita, y = hlm_nulo_fitted),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y = receita), method = "lm",
color = "gray44", size = 1.05,
linetype = "longdash") +
labs(x = "Receita real", y = "Valores Preditos") +
theme_bw()## `geom_smooth()` using formula 'y ~ x'
Observe que a linha azul (modelo nulo) chega a tocar a linha pontilhada (resultado perfeito), mas apresenta considerável diferença no início e no final da reta.
Modelo com interceptos aleatórios
O próximo passo, é a modelo com interceptos aleatórios que considera apenas efeito que a cidade (intercepto) exerce sobre a receita :
\[ receita_{ij} = \gamma_{00} + \gamma_{10} * nps_{ij} + \nu_{0j} + \epsilon_{ij} \]
Onde:
\(\gamma_{00}\) representa a média geral da variável receita
\(\gamma_{10}\) representa a inclinação da reta de receita de acordo com o avanço do NPS
\(\nu_{0j}\) representa a influência da cidade sobre a receita
\(\nu_{ij}\) representa a distância entre a receita calculada para aquela observação e o valor real
No R a fórmula pode ser representada como no exemplo abaixo:
hlm_intercept_aleatorios =
lme(fixed = receita ~ nps + inflacao,
random = ~1 | cidade,
data = vendas,
method = "REML")
#observação do BIC do modelo
BIC(hlm_intercept_aleatorios)## [1] 13421.34| Modelo | BIC |
|---|---|
| Modelo Nulo | 14212.24 |
| Modelo com interceptos aleatórios | 13421,35 |
O indicador BIC, que no modelo nulo, é de 14212,24 foi reduzido para 13421,35 indicando que a variável NPS possui grande efeito sobre a receita.
stderr_nlme(hlm_intercept_aleatorios)## RE.Components Variance.Estimatives Std.Err. z p.value
## 1 Var(v0j) 142365.8 52902.868 2.691079 0.007
## 2 Var(e) 110375.8 5185.945 21.283637 0.000O p-value do nível v0j segue abaixo de 0,05 indicando que o modelo é estatisticamente significativo.
Vamos verificar agora como os modelos se aproximam do valor ideal:
teste$hlm_intercept_aleatorios_fitted = predict(hlm_intercept_aleatorios, teste)
teste %>%
mutate_if(is.numeric, scale) %>%
ggplot() +
geom_smooth(aes(x = receita, y = hlm_nulo_fitted, color = "Modelo Nulo"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y= hlm_intercept_aleatorios_fitted, color = "Interceptos Aletaórios"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y = receita), method = "lm",
color = "gray44", size = 1.05,
linetype = "longdash") +
scale_color_manual("Modelos:",
values = c("deepskyblue1", "maroon1")) +
labs(x = "Receita Real", y = "Receita Prevista") +
theme_bw()## `geom_smooth()` using formula 'y ~ x'
Observe que a linha do modelo com interceptos aleatórios chegou muito perto da reta ideal
Modelo com interceptos e inclinações aleatórias
O próximo passo é verificarmos o efeito que a cidade exerce sobre o NPS e a receita e pode ser descrita pela função:
\[ receita_{ij} = \gamma_{00} + \gamma_{10} * nps_{ij} + \nu_{0j} + \nu_{1j} * nps_{ij} + \epsilon_{ij} \]
Onde:
\(\gamma_{00}\) representa a média geral da variável receita
\(\gamma_{10}\) representa a inclinação da reta de receita de acordo com o avanço do NPS
\(\nu_{0j}\) representa a influência da cidade sobre a receita
\(\nu_{ij}\) representa a distância entre a receita calculada para aquela observação e o valor real
Observe que agora temos \(\nu_{ij}\) influenciando o valor de NPS: \(\nu_{ij} * nps_{ij}\)
No R teremos:
hlm_intercept_inclin =
lme(fixed = receita ~ nps + inflacao,
random = ~ nps | cidade,
data = vendas,
method = "REML")
BIC(hlm_intercept_inclin)## [1] 13299.13Nesta terceira execução observa-se uma nova queda no indicador BIC, mostrando que o volume de erros diminuiu em comparação com os modelos anteriores
| Modelo | BIC |
|---|---|
| Modelo Nulo | 14212 |
| Interceptos aleatórios | 13421 |
| Interceptos e inclinações aleatórias | 13299 |
stderr_nlme(hlm_intercept_inclin)## Components Variance.Estimatives Std.Err. z p.value
## 1 Var(v0j) 658818.6994 302905.1560 2.17500 0.030
## 2 Var(v1j) 123.7548 47.6011 2.59983 0.009
## 3 Var(e) 90701.3200 4301.1633 21.08763 0.000O p-value do nível v0j (nível loja) e v1j seguem abaixo de 0,05 indicando que o modelo continua estatisticamente significativo.
Por fim vamos comparar a performance dos modelos
teste$hlm_intercept_inclin_fitted = predict(hlm_intercept_inclin, teste)
teste %>%
mutate_if(is.numeric, scale) %>%
ggplot() +
geom_smooth(aes(x = receita, y = hlm_nulo_fitted, color = "Modelo Nulo"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y= hlm_intercept_aleatorios_fitted, color = "Interceptos Aletaórios"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y= hlm_intercept_inclin_fitted, color = "Interceptos e Inclinações aleatórias"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y = receita), method = "lm",
color = "gray44", size = 1.05,
linetype = "longdash") +
scale_color_manual("Modelos:",
values = c("deepskyblue1", "maroon1", "darkorchid2")) +
labs(x = "Receita Real", y = "Receita Prevista") +
theme_bw()## `geom_smooth()` using formula 'y ~ x'
Modelo Final
Agora que verificamos que o NPS tem efeito sobre a receita e que o agrupamento por cidade também exerce efeito sobre as demais informações, vamos verificar se as variáveis do nível de cidade, neste caso a inflação local também impacta no resultado. Teremos finalmente a fórmula completa:
\[ receita_{ij} = \gamma_{00} + \gamma_{10} * nps_{ij} + \gamma_{01} * inflacao_{j} + \gamma_{11} * nps_{ij} * inflacao_j + \nu_{0j} + \nu_{1j} * nps_{ij} + \epsilon_{ij} \]
No R teremos:
hlm_final =
lme(fixed = receita ~ nps + inflacao + nps:inflacao,
random = ~ nps | cidade,
data = vendas,
method = "REML")
BIC(hlm_final)## [1] 13281.57Observamos uma nova redução no BIC, indicando que o houve melhora no modelo ao incluir o efeito da inflação na receita e no NPS.
| Modelo | BIC |
|---|---|
| Modelo Nulo | 14212 |
| Interceptos aletórios | 13421 |
| Interceptos e inclinações aleatórias | 13299 |
| Modelo Final | 13281 |
stderr_nlme(hlm_final)## Components Variance.Estimatives Std.Err. z p.value
## 1 Var(v0j) 392327.54923 163961.75910 2.392799 0.017
## 2 Var(v1j) 56.82155 24.49188 2.320016 0.020
## 3 Var(e) 90634.81776 4294.83656 21.103205 0.000O p-value do nível v0j (nível loja) e v1j seguem abaixo de 0,05 ou seja, modelo significativo estatisticamente. Vamos verificar visualmente como essa nota reta se porta com a base de testes
teste$hlm_final_fitted = predict(hlm_final, teste)
teste %>%
mutate_if(is.numeric, scale) %>%
ggplot() +
geom_smooth(aes(x = receita, y = hlm_nulo_fitted, color = "1 - Modelo Nulo"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y= hlm_intercept_aleatorios_fitted, color = "2 - Interceptos Aletaórios"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y= hlm_intercept_inclin_fitted, color = "3 - Interceptos e Inclinações aleatórias"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y= hlm_final_fitted, color = "4 - Modelo Final"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y = receita), method = "lm",
color = "gray44", size = 1.05,
linetype = "longdash") +
scale_color_manual("Modelos:",
values = c("deepskyblue1", "maroon1", "darkorchid2", "SpringGreen3")) +
labs(x = "Receita Real", y = "Receita Prevista") +
theme_bw()## `geom_smooth()` using formula 'y ~ x'
Ao compararmos as retas de regressão observamos que o modelo final (verde) performa melhor em comparação com todos os demais, ao considerar o efeito de todas as variáveis disponíveis.
Bônus
Vamos comparar o modelo HLM com um modelo re regressão linar, utilizando todas as técnicas necessárias para aumentar sua acuracidade:
Sabendo que temos uma variável qualitativa (a cidade) utilizaremos o método dummy_cols para transformar as colunas de forma que possam ser utilizadas em uma regressão linar. A “dumização” transformará a coluna cidade em uma série de colunas (uma para cada cidade) preenchendo seu valor com 0 ou 1 dependendo do valor da coluna original
vendas_dummies =
dummy_cols(.data = vendas,
select_columns = "cidade",
remove_first_dummy = TRUE,
remove_selected_columns = TRUE)
#16 colunas cidade foram criadas, o método dummy_cols não inclui a primeira coluna "dumizada" já que ela é a auzência de dados em todas as demais 16Executamos o método lmpara que seja criado o modelo de regressão simples, observe que agora cada coluna cidade deve ser informado na fórmula que seria:
\[ Y = \gamma_0 + NPS_1 * \gamma_1 + inflacao * \gamma_2 + cidade_1 * \gamma_3 * cidade_2 * \gamma_4 + ... + cidade_n * \gamma_n + \epsilon \]
No R teremos a fórmula:
ols_final_dummies =
lm(formula = receita ~ nps + inflacao + cidade_2 + cidade_3 + cidade_4 + cidade_5 + cidade_6 + cidade_7 + cidade_8 + cidade_9 + cidade_10 + cidade_11 + cidade_12 + cidade_13+ cidade_14+ cidade_15+ cidade_16+ cidade_17,
data = vendas_dummies)Ao criar variáveis dummies é necessário verificar se uma ou mais desta variáveis possuem significância estatística, para isso usamos o comando summarypara verificar o p-value de cada um, toda variável com p-value maior que 0,05 deve ser removida do modelo
summary(ols_final_dummies)##
## Call:
## lm(formula = receita ~ nps + inflacao + cidade_2 + cidade_3 +
## cidade_4 + cidade_5 + cidade_6 + cidade_7 + cidade_8 + cidade_9 +
## cidade_10 + cidade_11 + cidade_12 + cidade_13 + cidade_14 +
## cidade_15 + cidade_16 + cidade_17, data = vendas_dummies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2514.31 -140.56 8.63 151.16 1333.93
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.344e+04 1.083e+03 -12.412 < 2e-16 ***
## nps 2.852e+01 8.283e-01 34.426 < 2e-16 ***
## inflacao 1.270e+05 1.026e+04 12.377 < 2e-16 ***
## cidade_2 1.896e+01 1.102e+02 0.172 0.8634
## cidade_3 -9.563e+01 1.103e+02 -0.867 0.3862
## cidade_4 -4.043e+02 8.725e+01 -4.634 4.11e-06 ***
## cidade_5 3.096e+02 1.257e+02 2.464 0.0139 *
## cidade_6 6.555e+01 6.813e+01 0.962 0.3362
## cidade_7 2.746e+01 7.334e+01 0.374 0.7082
## cidade_8 -1.934e+03 2.371e+02 -8.157 1.14e-15 ***
## cidade_9 9.821e+02 2.005e+02 4.899 1.14e-06 ***
## cidade_10 -7.141e+02 5.836e+01 -12.236 < 2e-16 ***
## cidade_11 -1.409e+02 6.263e+01 -2.249 0.0247 *
## cidade_12 5.940e+01 1.118e+02 0.531 0.5954
## cidade_13 -1.421e+02 1.131e+02 -1.256 0.2093
## cidade_14 -9.191e+02 1.178e+02 -7.799 1.71e-14 ***
## cidade_15 1.874e+03 2.508e+02 7.470 1.89e-13 ***
## cidade_16 3.923e+02 8.171e+01 4.801 1.84e-06 ***
## cidade_17 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 332.2 on 906 degrees of freedom
## Multiple R-squared: 0.8843, Adjusted R-squared: 0.8821
## F-statistic: 407.3 on 17 and 906 DF, p-value: < 2.2e-16Para fazer a remoção destas variáveis podemos utilizar o método stepwise, que fará a remoção e recalculará o modelo apenas com aquelas variáveis que possuem a capacidade de explicar o modelo
ols_final_dummies_step =
step(object = ols_final_dummies,
step = qchisq(p = 0.05, df = 1,lower.tail = FALSE))## Start: AIC=10746.98
## receita ~ nps + inflacao + cidade_2 + cidade_3 + cidade_4 + cidade_5 +
## cidade_6 + cidade_7 + cidade_8 + cidade_9 + cidade_10 + cidade_11 +
## cidade_12 + cidade_13 + cidade_14 + cidade_15 + cidade_16 +
## cidade_17
##
##
## Step: AIC=10746.98
## receita ~ nps + inflacao + cidade_2 + cidade_3 + cidade_4 + cidade_5 +
## cidade_6 + cidade_7 + cidade_8 + cidade_9 + cidade_10 + cidade_11 +
## cidade_12 + cidade_13 + cidade_14 + cidade_15 + cidade_16
##
## Df Sum of Sq RSS AIC
## - cidade_2 1 3270 100003054 10745
## - cidade_7 1 15473 100015257 10745
## - cidade_12 1 31152 100030936 10745
## - cidade_3 1 82951 100082735 10746
## - cidade_6 1 102172 100101956 10746
## - cidade_13 1 174243 100174027 10747
## <none> 99999784 10747
## - cidade_11 1 558331 100558114 10750
## - cidade_5 1 669874 100669658 10751
## - cidade_4 1 2370248 102370031 10767
## - cidade_16 1 2544440 102544224 10768
## - cidade_9 1 2648615 102648399 10769
## - cidade_15 1 6158538 106158321 10800
## - cidade_14 1 6714231 106714015 10805
## - cidade_8 1 7343946 107343730 10810
## - cidade_10 1 16525380 116525163 10886
## - inflacao 1 16909397 116909180 10889
## - nps 1 130813027 230812811 11518
##
## Step: AIC=10745.01
## receita ~ nps + inflacao + cidade_3 + cidade_4 + cidade_5 + cidade_6 +
## cidade_7 + cidade_8 + cidade_9 + cidade_10 + cidade_11 +
## cidade_12 + cidade_13 + cidade_14 + cidade_15 + cidade_16
##
## Df Sum of Sq RSS AIC
## - cidade_7 1 20346 100023400 10743
## - cidade_12 1 41780 100044834 10743
## - cidade_3 1 113279 100116333 10744
## - cidade_6 1 118867 100121921 10744
## <none> 100003054 10745
## - cidade_13 1 284689 100287743 10746
## - cidade_11 1 603296 100606350 10749
## - cidade_5 1 2096269 102099323 10762
## - cidade_4 1 2767891 102770945 10768
## - cidade_16 1 4742302 104745356 10786
## - cidade_9 1 9452588 109455642 10826
## - cidade_14 1 10763230 110766284 10838
## - cidade_10 1 16601943 116604997 10885
## - cidade_8 1 18475157 118478211 10900
## - cidade_15 1 21443025 121446079 10922
## - inflacao 1 56343070 156346124 11156
## - nps 1 130824194 230827247 11516
##
## Step: AIC=10743.2
## receita ~ nps + inflacao + cidade_3 + cidade_4 + cidade_5 + cidade_6 +
## cidade_8 + cidade_9 + cidade_10 + cidade_11 + cidade_12 +
## cidade_13 + cidade_14 + cidade_15 + cidade_16
##
## Df Sum of Sq RSS AIC
## - cidade_12 1 46998 100070398 10742
## - cidade_6 1 110240 100133639 10742
## - cidade_3 1 186673 100210073 10743
## <none> 100023400 10743
## - cidade_13 1 438006 100461406 10745
## - cidade_11 1 671247 100694647 10747
## - cidade_5 1 2230779 102254178 10762
## - cidade_4 1 3210063 103233463 10770
## - cidade_16 1 4734930 104758330 10784
## - cidade_9 1 10430538 110453938 10833
## - cidade_14 1 12824578 112847978 10853
## - cidade_10 1 18406072 118429472 10897
## - cidade_8 1 23029683 123053083 10933
## - cidade_15 1 23793717 123817117 10938
## - inflacao 1 67164122 167187522 11216
## - nps 1 130943242 230966641 11514
##
## Step: AIC=10741.63
## receita ~ nps + inflacao + cidade_3 + cidade_4 + cidade_5 + cidade_6 +
## cidade_8 + cidade_9 + cidade_10 + cidade_11 + cidade_13 +
## cidade_14 + cidade_15 + cidade_16BIC(ols_final_dummies_step)## [1] 13443.09| Modelo | BIC |
|---|---|
| Modelo Nulo | 14212 |
| Interceptos aleatórios | 13421 |
| Interceptos e inclinações aleatórias | 13299 |
| Modelo Final | 13281 |
| OLS com apenas c/ variáveis significativas | 13443 |
Para finalizar vamos comparar o modelo final com o modelo OLS com variáveis dummies e stepwise:
teste_dummies =
dummy_cols(.data = teste,
select_columns = "cidade",
remove_first_dummy = TRUE,
remove_selected_columns = TRUE)
teste_dummies$cidade_2 = NULL
teste_dummies$cidade_7 = NULL
teste_dummies$cidade_12 = NULL
teste_dummies$cidade_17 = NULL
teste$ols_completo_fitted = predict(ols_final_dummies_step, teste_dummies)
teste %>%
mutate_if(is.numeric, scale) %>%
ggplot() +
geom_smooth(aes(x = receita, y = hlm_final_fitted, color = "1 - HLM Modelo Final"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y= ols_completo_fitted, color = "2 - OLS c/ Dummies e Stepwise"),
method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
size = 1.5) +
geom_smooth(aes(x = receita, y = receita), method = "lm",
color = "gray44", size = 1.05,
linetype = "longdash") +
scale_color_manual("Modelos:",
values = c("deepskyblue1", "maroon1")) +
labs(x = "Receita Real", y = "Receita Prevista") +
theme_bw()## `geom_smooth()` using formula 'y ~ x'
Conclusão
A regressão multinível é capaz de considerar melhor todas os aspectos dos dados de treino do modelo, criando um modelo capaz de prever com maior exatidão valores futuros.
Uma outra forma de verificar a qualidade do modelo é através da verificação da raiz quadrada do erro médio (RMSE em inglês), uma métrica utilizada para comparar modelos indicando sua acuracidade (quando menor seu valor melhor o resultado) e reforça a capacidade do modelo multinível em prever resultados com menor erro, neste caso ele é 5x mais eficiente que o modelo tradicional.
teste %>%
mutate_if(is.numeric, scale) %>%
summarise(hlm = sqrt(mean((receita - hlm_final_fitted) ^ 2,na.rm = T)),
ols = sqrt(mean((receita - ols_completo_fitted) ^ 2,na.rm = T))) %>%
melt() %>%
ggplot(aes(y = value, x=variable, fill = factor(variable))) +
geom_bar(stat = "identity") +
geom_label(aes(label = (round(value,3))), hjust = 1.2, color = "white", size = 5) +
coord_flip() +
labs(title="RMSE",
x = "modelos",
y = "") +
theme_minimal()## Warning in melt(.): The melt generic in data.table has been passed a data.frame
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(.). In the
## next version, this warning will become an error.## No id variables; using all as measure variables
Infelizmente desta vez eu não posso divulgar os dados originais, mas com este mesmo código, apenas substituindo as variáveis receita, NPS, inflação e cidade é possível aplicá-la a diversos cenários de dados.
Caso encontre alguma falha nos algoritmos ou um erro no processo, por favor me avise no meu linkedin (https://www.linkedin.com/in/marcosmhs/) e ficarei feliz realizar as correções necessárias e incluí-lo neste texto.
Até a próxima!