#### Análise exploratória, manipulação e visualização de dados ####
# Instalação de pacotes
# install.packages("tidyverse")
# install.packages("plotly")
# Carregamento dos pacotes
library(tidyverse)
library(plotly)
#### Modelagem geral ####
# Instalação de pacotes
# install.packages("easystats")
# Carregamento dos pacotes
library(easystats)
#### Modelagem Frequentista ####
# Instalação de pacotes
# install.packages("nlme")
# install.packages("lme4")
# Carregamento dos pacotes
library(nlme)
library(lme4)
library(lmerTest)
#### Modelagem Bayesiana ####
# Instalação de pacotes
# install.packages("rstanarm")
# install.packages("brms")
# Carregamento dos pacotes
library(rstanarm)
library(brms)Análise Comparativa entre Abordagens Frequentistas e Bayesianas em Modelos Lineares Generalizados Hierárquicos: Um tutorial com R
1 Jornada do Cientista de Dados
Este tutorial pretende seguir a técnica do modelo narrativo da Jornada do Herói de Campbell (1949). O tutorial está dividido em etapas que, gradativamente, pretendem orientar o Cientista de Dados, que foi chamado para essa aventura. Ao longo da jornada, o cientista de dados deve vencer as etapas e definir o modelo campeão, para também ser o vencedor da jornada e estar preparado pra o próximo “chamado”.
O Cientista de Dados será um vencedor, quando eleger o modelo campeão com justificativas de valor, ou seja, quando agregar às análises, um (ou mais) insight (Dykes, 2023). Insight, do inglês medieval para inner sight, ilustra uma “visão interior”), ou, “visão com os olhos da mente” (Online Etymology Dictionary, 2019). Este tutorial, associa insight a “mudanças inesperadas” nas antigas crenças do Cientista de Dados, do Consumidor do Modelo Campeão, do Investigador Principal, e, quando existir, do Patrocinador (financeiro) do negócio. Crenças no processo de conhecimento. Há, portanto, de se considerar, insights valiosos (que agregam valores para todos os envilvidos na jornada do Cientista de Dados). E, descreve o livro como pretendemos desenvolver o curso, ou seja, este tutorial, pretende, focar em geração de insights significativos que ofereçam alguma promessa tangível de valor ao leitor e a todos os usuários do mesmo.
DYKES, Brent. Dados e Storytelling de Impacto (Effective Data Storytelling). 1a ed. 2023. São Paulo: Benvirá, 2023. Edição do Kindle.
1.1 Introdução – O Chamado à aventura
Este tutorial tem como objetivo apresentar uma análise comparativa entre abordagens frequentistas e bayesianas em Modelos Lineares Generalizados Hierárquicos (MLGHs) em um contexto aplicado, com intuito de ser utilizado como um guia.
Este estudo busca resolver um problema de negócio fictício:
Cenário:
“Um produtor agrícola, proprietário de várias terras no Brasil, e, necessita expandir sua produção para a indústria cafeeira. É importante conhecer a área total (conjunta) das propriedades que devem ser disponibilizadas para o plantio de café. A área destinada terá que ser calculada em função do retorno médio de produção do café. A pergunta de negócio então versa sobre a estimativa do valor da produção de café com base na área plantada. Este cenário, embora fictício, conta com não identificação dos Estados e no valor total de área das propriedades”.
Pergunta de Negócio:
“É possivel prever o valor da produção de café com base na área plantada, considerendo uma estrutura hierárquica de Estado e Município?”.
A base de dados de Produção Agrícola Municipal (PAM) Instituto Brasileiro de Geografia e Estatística (IBGE) foi utilizada para compor a Analytical Base Table (ABT) deste cenário.
Exploratory Data Analysis (EDA) e Data Viz (Tabelas e Gráficos):
Análise Exploratória dos Dados da Produção Agrícola Municipal (PAM).
Modelagem Hierárquica:
Ajuste do Modelo Linear Generalizado Hierárquico (MLGH) com efeitos fixos e aleatórios considerando “campeão”.
Para a escolha do modelo ótimo (campeão) serão consideradas as análises com as duas abordagens de estimação, quais sejam:
Modelagem Frequentista: Modelo Linear Generalizado Hierárquico (MLGH) ajustado com métodos e técnicas frequentistas, utilizando a função lme() do pacote nlme e a função glmer() do pacote lme4.
Modelagem Bayesiana: Modelo Linear Generalizado Hierárquico (MLGH) ajustado com métodos e técnicas bayesianas, utilizando a função stan_lmer() do pacote rstanarm e a função brm() do pacote brms.
2 A ajuda necessária aos dilemas técnicos da Jornada
2.1 Carregamento dos Pacotes
2.2 Carregamento dos Dados
# Carregamento dos dados
df <- read_csv("../data/pesquisa_agricola_municipal.zip")2.3 Data Wrangling – A construção da ABT
Vamos filtrar apenas os dados referentes ao produto “Café” e remover o objeto df para liberar memória.
# Filtrar apenas os dados de Café
df_cafe <- df %>%
filter(str_detect(produto, "Café"))
# Remover o objeto df para liberar memória
rm(df)Filtrar apenas dados de 2019 (último ano que consta no dataset).
# Filtrar apenas dados de 2019
df_cafe <- df_cafe %>%
filter(ano == 2019)Selecionar as variáveis de interesse para a análise.
# Selecionar variáveis de interesse
df_cafe <- df_cafe %>%
select(sigla_uf, id_municipio, area_plantada, valor_producao)Remover dados faltantes.
# Remover dados faltantes
df_cafe <- df_cafe %>%
drop_na()Filtrar os 5 estados com maior produção de café.
# Filtrar os 5 estados com maior produção de café
df_cafe %>%
group_by(sigla_uf) %>%
summarise(total_producao = sum(valor_producao)) %>%
top_n(5, total_producao) %>%
select(sigla_uf) %>%
pull() %>%
unique() -> estados
df_cafe <- df_cafe %>%
filter(sigla_uf %in% estados)Filtrar áreas plantadas menores que 30.000 hectares.
# Filtrar áreas plantadas menores que 30.000 hectares
df_cafe <- df_cafe %>%
filter(area_plantada <= 30000)3 Batalhas
3.1 Análise Descritiva – EDA e Data Viz
Variáveis do dataset:
# Estatísticas descritivas
df_cafe %>%
report() %>%
summaryThe data contains 2276 observations of the following 4 variables:
- sigla_uf: 5 entries, such as MG (47.36%); SP (27.86%); BA (11.99%) and 2
others
- id_municipio: Mean = 3.13e+06, SD = 4.85e+05, range: [1100015, 3557303]
- area_plantada: Mean = 1500.37, SD = 2844.43, range: [1, 20920]
- valor_producao: Mean = 14573.91, SD = 29316.26, range: [4, 209058]A variável que iremos prever é valor_producao, que representa o valor da produção de café por Município.
Vamos visualizar a distribuição da variável valor_producao:
# Visualização da distribuição da variável valor_producao
plot_1 <- df_cafe %>%
ggplot(aes(x = valor_producao/1000)) +
geom_histogram(bins = 50) +
labs(x = "Valor de Produção (R$)",
y = "Frequência",
title = "Distribuição do Valor de Produção",
subtitle = "Dados de Produção Agrícola Municipal (PAM)",
caption = "Fonte: IBGE") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5)) +
scale_x_continuous(labels = scales::number_format(prefix = "R$",
suffix = " Mil"))
#ggplotly(plot_1)
plot_1
Hipóteses:
- A maioria dos municípios tem um valor de produção na faixa mais baixa. Isso sugere que uma grande quantidade de municípios produz um valor relativamente pequeno de café.
- O gráfico apresenta uma cauda longa à direita, indicando que há poucos municípios com valores de produção muito altos. Isso pode sugerir que a indústria do café nos estados avaliados é composta por muitos pequenos produtores e alguns poucos produtores muito grandes.
A principal variável preditora que iremos utilizar é area_plantada, que representa a área plantada com café em cada município.
Vamos visualizar a distribuição da variável area_plantada:
# Visualização da distribuição da variável area_plantada
plot_2 <- df_cafe %>%
ggplot(aes(x = area_plantada)) +
geom_histogram(bins = 50) +
labs(x = "Área Plantada (ha)",
y = "Frequência",
title = "Distribuição da Área Plantada",
subtitle = "Dados de Produção Agrícola Municipal (PAM)",
caption = "Fonte: IBGE") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5))
#ggplotly(plot_2)
plot_2
Hipóteses:
- Assim como a distribuição do valor da produção, há uma desigualdade significativa nas áreas de plantio entre os municípios. A maioria tem uma área plantada pequena, enquanto poucos têm grandes extensões de terra dedicadas à agricultura.
Vamos avaliar a variável principal variável preditora, area_plantada, em relação à variável resposta valor_producao:
# Visualização da relação entre area_plantada e valor_producao
plot_3 <- df_cafe %>%
ggplot(aes(x = (area_plantada), y = (valor_producao)/1000)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(x = "Área Plantada (ha)",
y = "Valor de Produção (R$)",
title = "Relação entre Área Plantada e Valor de Produção",
subtitle = "Dados de Produção Agrícola Municipal (PAM)",
caption = "Fonte: IBGE") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5)) +
scale_x_continuous(labels = scales::number_format(suffix = " ha")) +
scale_y_continuous(labels = scales::number_format(prefix = "R$",
suffix = " Mil"))
#ggplotly(plot_3)
plot_3
Hipóteses:
- Há uma relação positiva entre a área plantada e o valor da produção. Isso sugere que, em geral, municípios com áreas maiores de plantio tendem a ter valores de produção mais altos, o que faz total sentido.
- No entanto, a dispersão dos pontos indica que a relação não é perfeita. Existem municípios com áreas menores que produzem mais do que municípios com áreas maiores. Isso pode ser devido a diferenças na produtividade, qualidade do solo, clima ou outros fatores.
- A linha vermelha representa a linha de regressão linear ajustada aos dados. Ela mostra a tendência geral da relação entre as variáveis, mas não captura toda a variação nos dados.
Insights:
- Ao que parece, podemos usar a área plantada como possível variável preditora para o valor da produção de café.
- No entanto, a relação entre essas variáveis não é perfeita, o que sugere que outros fatores também podem influenciar o valor da produção.
- É possível observar um padrão de dispersão dos pontos, indicando que a força de associação (inclinação - slope) entre alguma categoria pode variar. Uma observação interessante é que o intercepto é o mesmo para todas as categorias, e isso é uma característica da natureza dos dados: se a área plantada for zero, o valor da produção também será zero.
Vamos visualizar a quantidade de municípios em cada estado:
# Visualização da quantidade de municípios em cada estado
plot_4 <- df_cafe %>%
group_by(sigla_uf) %>%
summarise(n_municipios = n_distinct(id_municipio)) %>%
arrange(desc(n_municipios)) %>%
ggplot(aes(x = reorder(sigla_uf, n_municipios, decreasing = TRUE),
y = n_municipios)) +
geom_col() +
labs(x = "Estado",
y = "Quantidade de Municípios",
title = "Quantidade de Municípios por Estado",
subtitle = "Dados de Produção Agrícola Municipal (PAM)",
caption = "Fonte: IBGE") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5)) +
# adicionar labels com a quantidade de municipios ####
geom_text(aes(label = n_municipios),
vjust = -0.5,
size = 3)
#ggplotly(plot_4)
plot_4
Vamos avaliar a relação entre a variável resposta valor_producao e a variável preditora area_plantada por município:
# Visualização da relação entre area_plantada e valor_producao por município
plot_7 <- df_cafe %>%
ggplot(aes(x = (area_plantada),
y = (valor_producao)/1000,
color = factor(id_municipio))) +
geom_point(alpha = 0.1) +
geom_smooth(method = "lm", se = FALSE) +
labs(x = "Área Plantada (ha)",
y = "Valor de Produção (R$)",
title = "Relação entre Área Plantada e Valor de Produção por município",
subtitle = "Dados de Produção Agrícola Municipal (PAM)",
caption = "Fonte: IBGE") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5)) +
scale_x_continuous(labels = scales::number_format(suffix = " ha")) +
scale_y_continuous(labels = scales::number_format(prefix = "R$",
suffix = " Mil")) +
# excluir a legenda
theme(legend.position = "none")
#ggplotly(plot_7)
plot_7
Vamos avaliar a relação entre a variável resposta valor_producao e a variável preditora area_plantada em cada estado:
# Visualização da relação entre area_plantada e valor_producao por estado
plot_6 <- df_cafe %>%
ggplot(aes(x = (area_plantada),
y = (valor_producao)/1000,
color = sigla_uf)) +
geom_point(alpha = 0.1) +
geom_smooth(method = "lm", se = FALSE) +
labs(x = "Área Plantada (ha)",
y = "Valor de Produção (R$)",
title="Relação entre Área Plantada e Valor de Produção por Estado",
subtitle = "Dados de Produção Agrícola Municipal (PAM)",
caption = "Fonte: IBGE") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5)) +
scale_x_continuous(labels = scales::number_format(suffix = " ha")) +
scale_y_continuous(labels = scales::number_format(prefix = "R$",
suffix = " Mil"))
#ggplotly(plot_6)
plot_6
Hipóteses:
A relação entre a área plantada e o valor da produção parece ser consistente em todos os estados e municípios. A maioria dos estados e municípios mostram uma tendência positiva entre essas variáveis, o que sugere que áreas maiores de plantio tendem a ter valores de produção mais altos.
É possível observar que as inclinações das linhas de regressão variam entre os estados e municípios. Isso sugere que a força da associação entre a área plantada e o valor da produção pode ser diferente em cada estado e município.
De uma maneira qualitativa, as diferenças entre entre os municipios parecem ser maiores do que entre os estados. Isso sugere que a variabilidade entre os municípios é maior do que entre os estados. Iremos confirmar essa hipótese com a modelagem dos dados.
Hipóteses para a modelagem:
Com base nas análises exploratórias realizadas, podemos destacar alguns questionamentos importantes para o racional da modelagem dos dados de Produção Agrícola Municipal (PAM) de café:
- Avaliar se os níveis hierárquicos (estados e municípios) têm efeitos significativos na relação entre a área plantada e o valor da produção para intercepção e ou inclinação.
Com a análise exploratória e a limpeza dos dados realizadas, estamos prontos para avançar para a modelagem.
3.2 Modelagem
Para fins didáticos, vamos considerar apenas modelos de regressão linear com efeitos aleatórios para estados e municípios, a fim de avaliar a influência desses níveis hierárquicos na relação entre a área plantada e o valor da produção. Com isso, partimos do princípio que modelos mais simples como lineares de efeitos fixos e modelos nulos (sem variáveis independentes) já foram rodados. Vamos apenas considerar os modelos com efeitos aleatórios para entender a sintaxe e a interpretação dos resultados.
Planejamento
Vamos ajustar os modelos descritos a seguir:
A. Nível municípios - um nível
1. Modelo com interceptos aleatórios para municípios
2. Modelo com inclinações aleatórias para municípios
3. Modelo com interceptos e inclinações aleatórias para municípios
B. Nível estados - um nível
4. Modelo com interceptos aleatórios para estados
5. Modelo com inclinações aleatórias para estados
6. Modelo com interceptos e inclinações aleatórias para estados
C. Níveis municípios e estados - dois níveis aninhados*
7. Modelo com interceptos aleatórios para estados e municípios
8. Modelo com inclinações aleatórias para estados e municípios
9. Modelo com interceptos e inclinações aleatórias para estados e municípios
* O aninhamento se dá pelo fato de que os municípios estão contidos nos estados.
Vamos ajustar os modelos seguindo esses passos:
- Ajustar o modelo hierárquico
- Sumarizar os resultados
- Avaliar a qualidadea do ajuste do modelo
- Log-Likelihood para modelos frequentistas
- ELPD para modelos bayesianos
- Log-Likelihood para modelos frequentistas
Log-Likelihood (LL): é uma medida de ajuste do modelo frequentista. Ao ajustar um modelo a um conjunto de dados, o log-likelihood quantifica quão prováveis são os dados observados sob as suposições do modelo. Em termos simples, ele mede o quão bem o modelo explica os dados observados.
ELPD ELPD significa Expected Log Pointwise Predictive Density. Essencialmente, o ELPD avalia quão bem um modelo pode prever novos dados, levando em conta a incerteza do modelo. A métrica é baseada na densidade preditiva, que mede a probabilidade dos dados observados sob o modelo bayesiano.
Diferenças
Foco no Conjunto de Dados Atual vs. Previsão para Novos Dados: O log-likelihood foca na avaliação de quão bem o modelo se ajusta ao conjunto de dados atual. Ele olha para os dados já observados e avalia a probabilidade desses dados sob o modelo. Por outro lado, o ELPD está preocupado com a capacidade do modelo de prever novos dados, não apenas se ajustar bem aos dados existentes. Isso o torna especialmente útil em cenários de modelagem preditiva.
Incorporação da Incerteza do Modelo: O ELPD, sendo uma métrica comum em abordagens Bayesianas, naturalmente incorpora a incerteza nos parâmetros do modelo ao fazer previsões. Ele fornece uma medida de desempenho que leva em conta tanto o ajuste dos dados quanto a complexidade do modelo, penalizando modelos excessivamente complexos. O log-likelihood, em sua forma básica, não incorpora diretamente a incerteza dos parâmetros ou a penalidade por complexidade, embora variantes como o AIC (Akaike Information Criterion) e o BIC (Bayesian Information Criterion) derivem do log-likelihood e incluam termos de penalidade para a complexidade do modelo.
Vamos ajustar os modelos utilizando as seguintes bibliotecas:
- Paradigma frequentista:
nlme
lme4
- Paradigma bayesiano:
rstan_arm
brms
3.2.1 MUNICÍPIOS
3.2.1.1 Modelo com interceptos aleatórios
Ajuste do modelo
nlme_int_mun <- nlme::lme(
valor_producao ~ area_plantada, # efeitos fixos
random = list(~ 1 | id_municipio), # efeitos aleatórios
data = df_cafe
)Sumário do modelo
summary(nlme_int_mun)Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
45263.81 45286.73 -22627.91
Random effects:
Formula: ~1 | id_municipio
(Intercept) Residual
StdDev: 8045.244 2300.991
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) 823.0116 266.79904 1183 3.08476 0.0021
area_plantada 9.2266 0.06527 1183 141.35993 0.0000
Correlation:
(Intr)
area_plantada -0.366
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-8.566360477 -0.020662430 -0.014276444 -0.003001754 11.680325752
Number of Observations: 2276
Number of Groups: 1092 LogLik do modelo
df_ll <- tibble(
modelo = "municipios - intercepto",
loglik = logLik(nlme_int_mun)[1],
biblioteca = "nlme"
)
df_ll %>% filter(biblioteca == "nlme")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | nlme |
Ajuste do modelo
lme4_int_mun <- lme4::lmer(
valor_producao ~ area_plantada + (1 | id_municipio),
data = df_cafe
)Sumário do modelo
summary(lme4_int_mun)Linear mixed model fit by REML ['lmerMod']
Formula: valor_producao ~ area_plantada + (1 | id_municipio)
Data: df_cafe
REML criterion at convergence: 45255.8
Scaled residuals:
Min 1Q Median 3Q Max
-8.5664 -0.0207 -0.0143 -0.0030 11.6803
Random effects:
Groups Name Variance Std.Dev.
id_municipio (Intercept) 64726105 8045
Residual 5294549 2301
Number of obs: 2276, groups: id_municipio, 1092
Fixed effects:
Estimate Std. Error t value
(Intercept) 823.01439 266.79931 3.085
area_plantada 9.22656 0.06527 141.360
Correlation of Fixed Effects:
(Intr)
area_plantd -0.366
fit warnings:
Some predictor variables are on very different scales: consider rescalingLogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - intercepto",
loglik = logLik(lme4_int_mun)[1],
biblioteca = "lme4"
)
)
df_ll %>% filter(biblioteca == "lme4")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | lme4 |
Ajuste do modelo
rstanarm_int_mun <- rstanarm::stan_lmer(
valor_producao ~ area_plantada + (1 | id_municipio),
data = df_cafe,
# esta linha é apenas para manter o exemplo pequeno!
chains = 2, cores = 4, seed = 12345, iter = 1000
)Sumário do modelo
rstanarm_int_munstan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (1 | id_municipio)
observations: 2276
------
Median MAD_SD
(Intercept) 825.9 290.5
area_plantada 9.2 0.1
Auxiliary parameter(s):
Median MAD_SD
sigma 2303.1 45.0
Error terms:
Groups Name Std.Dev.
id_municipio (Intercept) 8052
Residual 2305
Num. levels: id_municipio 1092
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregLogLik do modelo
loo_rstanarm_int_mun <- rstanarm::loo(rstanarm_int_mun)
df_ll <- rbind(
df_ll,tibble(
modelo = "municipios - intercepto",
loglik = loo_rstanarm_int_mun$estimates[1,1],
biblioteca = "rstanarm"
)
)
df_ll %>% filter(biblioteca == "rstanarm")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21255.82 | rstanarm |
Ajuste do modelo
Sumário do modelo
brms_int_mun Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (1 | id_municipio)
Data: df_cafe (Number of observations: 2276)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Group-Level Effects:
~id_municipio (Number of levels: 1092)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 8034.53 204.10 7618.62 8436.88 1.01 183 439
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 813.33 280.20 251.30 1318.76 1.04 77 176
area_plantada 9.22 0.07 9.07 9.36 1.01 109 241
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 2296.61 50.27 2201.70 2404.16 1.00 320 352
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).LogLik do modelo
loo_brms_int_mun <- brms::loo(brms_int_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - intercepto",
loglik = loo_brms_int_mun$estimates[1,1],
biblioteca = "brms"
)
)
df_ll %>% filter(biblioteca == "brms")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21250.25 | brms |
3.2.1.2 Modelo com inclinações aleatórias
Ajuste do modelo
nlme_slope_mun <- nlme::lme(
valor_producao ~ area_plantada, # efeitos fixos
random = list(~ 0 + area_plantada | id_municipio), # efeitos aleatórios
data = df_cafe
)Sumário do modelo
summary(nlme_slope_mun)Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
39790.22 39813.14 -19891.11
Random effects:
Formula: ~0 + area_plantada | id_municipio
area_plantada Residual
StdDev: 2.940357 922.6896
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) -67.21994 28.455246 1183 -2.36230 0.0183
area_plantada 9.50261 0.136933 1183 69.39595 0.0000
Correlation:
(Intr)
area_plantada -0.355
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-16.71184579 -0.05711394 0.01791640 0.06819056 13.67393999
Number of Observations: 2276
Number of Groups: 1092 LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - inclinação",
loglik = logLik(nlme_slope_mun)[1],
biblioteca = "nlme"
)
)
df_ll %>% filter(biblioteca == "nlme")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | nlme |
| municipios - inclinação | -19891.11 | nlme |
Ajuste do modelo
lme4_slope_mun <- lme4::lmer(
valor_producao ~ area_plantada + (0 + area_plantada | id_municipio),
data = df_cafe
)Sumário do modelo
summary(lme4_slope_mun)Linear mixed model fit by REML ['lmerMod']
Formula: valor_producao ~ area_plantada + (0 + area_plantada | id_municipio)
Data: df_cafe
REML criterion at convergence: 39782.2
Scaled residuals:
Min 1Q Median 3Q Max
-16.7132 -0.0571 0.0179 0.0682 13.6750
Random effects:
Groups Name Variance Std.Dev.
id_municipio area_plantada 8.65 2.941
Residual 851218.39 922.615
Number of obs: 2276, groups: id_municipio, 1092
Fixed effects:
Estimate Std. Error t value
(Intercept) -67.213 28.454 -2.362
area_plantada 9.503 0.137 69.380
Correlation of Fixed Effects:
(Intr)
area_plantd -0.355
fit warnings:
Some predictor variables are on very different scales: consider rescaling
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00775324 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - inclinação",
loglik = logLik(lme4_slope_mun)[1],
biblioteca = "lme4"
)
)
df_ll %>% filter(biblioteca == "lme4")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | lme4 |
| municipios - inclinação | -19891.11 | lme4 |
Ajuste do modelo
rstanarm_slope_mun <- rstanarm::stan_lmer(
valor_producao ~ area_plantada + (0 + area_plantada | id_municipio),
data = df_cafe,
# esta linha é apenas para manter o exemplo pequeno!
chains = 2, cores = 4, seed = 12345, iter = 1000
)Sumário do modelo
rstanarm_slope_munstan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (0 + area_plantada | id_municipio)
observations: 2276
------
Median MAD_SD
(Intercept) -68.1 29.2
area_plantada 9.5 0.1
Auxiliary parameter(s):
Median MAD_SD
sigma 922.7 16.1
Error terms:
Groups Name Std.Dev.
id_municipio area_plantada 2.9
Residual 923.9
Num. levels: id_municipio 1092
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregLogLik do modelo
loo_rstanarm_slope_mun <- rstanarm::loo(rstanarm_slope_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - inclinação",
loglik = loo_rstanarm_slope_mun$estimates[1,1],
biblioteca = "rstanarm"
)
)
df_ll %>% filter(biblioteca == "rstanarm")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21255.82 | rstanarm |
| municipios - inclinação | -19013.34 | rstanarm |
Ajuste do modelo
Sumário do modelo
brms_slope_mun Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (0 + area_plantada | id_municipio)
Data: df_cafe (Number of observations: 2276)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Group-Level Effects:
~id_municipio (Number of levels: 1092)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(area_plantada) 2.96 0.10 2.78 3.16 1.03 165 381
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -68.13 27.85 -121.91 -13.91 1.00 762 746
area_plantada 9.51 0.14 9.26 9.82 1.01 86 173
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 922.23 15.96 892.31 954.34 1.01 783 665
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).LogLik do modelo
loo_brms_slope_mun <- brms::loo(brms_slope_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - inclinação",
loglik = loo_brms_slope_mun$estimates[1,1],
biblioteca = "brms"
)
)
df_ll %>% filter(biblioteca == "brms")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21250.25 | brms |
| municipios - inclinação | -19011.98 | brms |
3.2.1.3 Modelo com intercepto e inclinações aleatórias
Ajuste do modelo
nlme_int_slope_mun <- nlme::lme(
valor_producao ~ area_plantada, # efeitos fixos
random = list(~ area_plantada | id_municipio), # efeitos aleatórios
data = df_cafe
)Sumário do modelo
summary(nlme_int_slope_mun)Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
39794.22 39828.59 -19891.11
Random effects:
Formula: ~area_plantada | id_municipio
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 4.313362e-04 (Intr)
area_plantada 2.940357e+00 0.082
Residual 9.226897e+02
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) -67.21994 28.455247 1183 -2.36230 0.0183
area_plantada 9.50261 0.136933 1183 69.39597 0.0000
Correlation:
(Intr)
area_plantada -0.355
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-16.71184463 -0.05711397 0.01791640 0.06819056 13.67393906
Number of Observations: 2276
Number of Groups: 1092 LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - intercepto e inclinação",
loglik = logLik(nlme_int_slope_mun)[1],
biblioteca = "nlme"
)
)
df_ll %>% filter(biblioteca == "nlme")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | nlme |
| municipios - inclinação | -19891.11 | nlme |
| municipios - intercepto e inclinação | -19891.11 | nlme |
Ajuste do modelo
lme4_int_slope_mun <- lme4::lmer(
valor_producao ~ area_plantada + (area_plantada | id_municipio),
data = df_cafe
)Sumário do modelo
summary(lme4_int_slope_mun)Linear mixed model fit by REML ['lmerMod']
Formula: valor_producao ~ area_plantada + (area_plantada | id_municipio)
Data: df_cafe
REML criterion at convergence: 40264.3
Scaled residuals:
Min 1Q Median 3Q Max
-15.4175 -0.0254 0.0146 0.0355 15.0357
Random effects:
Groups Name Variance Std.Dev. Corr
id_municipio (Intercept) 873621.82 934.677
area_plantada 11.08 3.328 -0.51
Residual 764953.45 874.616
Number of obs: 2276, groups: id_municipio, 1092
Fixed effects:
Estimate Std. Error t value
(Intercept) -91.7078 44.1549 -2.077
area_plantada 9.4959 0.1507 62.992
Correlation of Fixed Effects:
(Intr)
area_plantd -0.427
fit warnings:
Some predictor variables are on very different scales: consider rescaling
optimizer (nloptwrap) convergence code: 0 (OK)
unable to evaluate scaled gradient
Model failed to converge: degenerate Hessian with 1 negative eigenvaluesLogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - intercepto e inclinação",
loglik = logLik(lme4_int_slope_mun)[1],
biblioteca = "lme4"
)
)
df_ll %>% filter(biblioteca == "lme4")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | lme4 |
| municipios - inclinação | -19891.11 | lme4 |
| municipios - intercepto e inclinação | -20132.16 | lme4 |
Ajuste do modelo
rstanarm_int_slope_mun <- rstanarm::stan_lmer(
valor_producao ~ area_plantada + (area_plantada | id_municipio),
data = df_cafe,
# esta linha é apenas para manter o exemplo pequeno!
chains = 2, cores = 4, seed = 12345, iter = 1000
)Sumário do modelo
rstanarm_int_slope_munstan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (area_plantada | id_municipio)
observations: 2276
------
Median MAD_SD
(Intercept) -67.4 28.7
area_plantada 9.5 0.1
Auxiliary parameter(s):
Median MAD_SD
sigma 922.6 15.4
Error terms:
Groups Name Std.Dev. Corr
id_municipio (Intercept) 32
area_plantada 3 -0.44
Residual 923
Num. levels: id_municipio 1092
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregLogLik do modelo
loo_rstanarm_int_slope_mun <- rstanarm::loo(rstanarm_int_slope_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - intercepto e inclinação",
loglik = loo_rstanarm_int_slope_mun$estimates[1,1],
biblioteca = "rstanarm"
)
)
df_ll %>% filter(biblioteca == "rstanarm")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21255.82 | rstanarm |
| municipios - inclinação | -19013.34 | rstanarm |
| municipios - intercepto e inclinação | -19014.47 | rstanarm |
Ajuste do modelo
Sumário do modelo
brms_int_slope_mun Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (area_plantada | id_municipio)
Data: df_cafe (Number of observations: 2276)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Group-Level Effects:
~id_municipio (Number of levels: 1092)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept) 8.89 8.45 0.28 31.66 1.00 1085
sd(area_plantada) 2.93 0.09 2.76 3.10 1.01 251
cor(Intercept,area_plantada) 0.94 0.06 0.77 1.00 1.51 4
Tail_ESS
sd(Intercept) 579
sd(area_plantada) 310
cor(Intercept,area_plantada) 14
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -67.10 28.95 -124.72 -11.63 1.00 827 876
area_plantada 9.50 0.13 9.25 9.76 1.01 124 301
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 923.81 15.93 893.72 954.71 1.00 854 759
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).LogLik do modelo
loo_brms_int_slope_mun <- brms::loo(brms_int_slope_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "municipios - intercepto e inclinação",
loglik = loo_brms_int_slope_mun$estimates[1,1],
biblioteca = "brms"
)
)
df_ll %>% filter(biblioteca == "brms")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21250.25 | brms |
| municipios - inclinação | -19011.98 | brms |
| municipios - intercepto e inclinação | -19034.87 | brms |
3.2.1.4 Comparação dos modelos para o nível de municípios
performance::compare_performance(nlme_int_mun,
nlme_slope_mun,
nlme_int_slope_mun, rank = TRUE)Random effect variances not available. Returned R2 does not account for random effects.| Name | Model | R2_conditional | R2_marginal | RMSE | Sigma | ICC | AIC_wt | AICc_wt | BIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|
| nlme_slope_mun | lme | NA | 0.9988361 | 808.4987 | 922.6896 | NA | 0.8807974 | 0.8818125 | 0.9995608 | 1.0000000 |
| nlme_int_slope_mun | lme | 0.9989628 | 0.8900703 | 808.4987 | 922.6897 | 0.9905650 | 0.1192026 | 0.1181875 | 0.0004392 | 0.3783003 |
| nlme_int_mun | lme | 0.9930223 | 0.9077204 | 1688.2508 | 2300.9911 | 0.9243856 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0270460 |
# intercept vs slope
lmtest::lrtest(nlme_int_mun, nlme_slope_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -22627.91 | NA | NA | NA |
| 4 | -19891.11 | 0 | 5473.592 | 0 |
# intercept vs intercept + slope
lmtest::lrtest(nlme_int_mun, nlme_int_slope_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -22627.91 | NA | NA | NA |
| 6 | -19891.11 | 2 | 5473.592 | 0 |
# slope vs intercept + slope
lmtest::lrtest(nlme_slope_mun, nlme_int_slope_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -19891.11 | NA | NA | NA |
| 6 | -19891.11 | 2 | 7.2e-06 | 0.9999964 |
performance::compare_performance(lme4_int_mun,
lme4_slope_mun,
lme4_int_slope_mun, rank = TRUE)| Name | Model | R2_conditional | R2_marginal | ICC | RMSE | Sigma | AIC_wt | AICc_wt | BIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|
| lme4_slope_mun | lmerMod | 0.9989630 | 0.8900178 | 0.9905715 | 808.4807 | 922.6150 | 1 | 1 | 1 | 0.92651 |
| lme4_int_slope_mun | lmerMod | 0.9990904 | 0.8675229 | 0.9931339 | 756.6665 | 874.6162 | 0 | 0 | 0 | 0.50000 |
| lme4_int_mun | lmerMod | 0.9930223 | 0.9077202 | 0.9243859 | 1688.2503 | 2300.9887 | 0 | 0 | 0 | 0.12500 |
# intercept vs slope
lmtest::lrtest(lme4_int_mun, lme4_slope_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -22627.91 | NA | NA | NA |
| 4 | -19891.11 | 0 | 5473.592 | 0 |
# intercept vs intercept + slope
lmtest::lrtest(lme4_int_mun, lme4_int_slope_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -22627.91 | NA | NA | NA |
| 6 | -20132.16 | 2 | 4991.486 | 0 |
# slope vs intercept + slope
lmtest::lrtest(lme4_slope_mun, lme4_int_slope_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -19891.11 | NA | NA | NA |
| 6 | -20132.16 | 2 | 482.1063 | 0 |
performance::compare_performance(rstanarm_int_mun,
rstanarm_slope_mun,
rstanarm_int_slope_mun, rank = TRUE)| Name | Model | R2 | R2_marginal | R2_adjusted | R2_adjusted_marginal | ICC | RMSE | Sigma | WAIC_wt | LOOIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|---|
| rstanarm_int_slope_mun | stanreg | 0.9990100 | 0.9988376 | 0.9990392 | 0.9990392 | 0.9907607 | 807.8838 | 922.6277 | 0.9999938 | 0.8866759 | 0.9999730 |
| rstanarm_slope_mun | stanreg | 0.9990092 | 0.9988391 | 0.9990364 | 0.9990364 | 0.9905499 | 808.9851 | 922.6888 | 0.0000062 | 0.1020848 | 0.7886586 |
| rstanarm_int_mun | stanreg | 0.9938302 | 0.9923635 | 0.9943984 | 0.9943984 | 0.9242599 | 1688.8128 | 2303.1155 | 0.0000000 | 0.0112393 | 0.0000000 |
rstanarm::loo_compare(loo_rstanarm_int_mun,
loo_rstanarm_slope_mun,
loo_rstanarm_int_slope_mun) elpd_diff se_diff
rstanarm_slope_mun 0.0 0.0
rstanarm_int_slope_mun -1.1 4.0
rstanarm_int_mun -2242.5 248.7performance::compare_performance(brms_int_mun,
brms_slope_mun,
brms_int_slope_mun, rank = TRUE)| Name | Model | R2 | R2_marginal | R2_adjusted | R2_adjusted_marginal | ICC | RMSE | Sigma | WAIC_wt | LOOIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|---|
| brms_slope_mun | brmsfit | 0.9990124 | 0.9163338 | 0.9990238 | 0.9013983 | 0.9907019 | 809.3677 | 927.7598 | 1 | 0.9889069 | 0.8917249 |
| brms_int_slope_mun | brmsfit | 0.9990070 | 0.9161204 | 0.9990054 | 0.9010739 | 0.9904779 | 809.8766 | 942.1286 | 0 | 0.0000036 | 0.6613630 |
| brms_int_mun | brmsfit | 0.9938522 | 0.9087231 | 0.9944675 | 0.9137811 | 0.9244654 | 1688.9258 | 2278.2195 | 0 | 0.0110894 | 0.1123567 |
brms::loo_compare(
brms::add_criterion(brms_int_mun,"loo"),
brms::add_criterion(brms_slope_mun,"loo"),
brms::add_criterion(brms_int_slope_mun,"loo")
) elpd_diff se_diff
brms::add_criterion(brms_slope_mun, "loo") 0.0 0.0
brms::add_criterion(brms_int_slope_mun, "loo") -22.9 9.2
brms::add_criterion(brms_int_mun, "loo") -2238.3 250.03.2.1.5 Sumarização dos resultados para nível de municípios
Os modelos com inclinação e intercepto+inclinação aleatórios se ajustaram melhor aos dados do que o modelo com apenas intercepto aleatório. Isso demonstra que existe varição importante entre os municípios para a força de associção (inclinação - slope) entre a área plantada e o valor da produção. O que corrobora com a nossa hipótese inicial para o nível de municípios.
Agora vamos avaliar o nível de estados.
3.2.2 ESTADOS
3.2.2.1 Modelo com interceptos aleatórios
Ajuste do modelo
nlme_int_est <- nlme::lme(
valor_producao ~ area_plantada,
random = list(~ 1 | sigla_uf),
data = df_cafe
)Sumário do modelo
summary(nlme_int_est)Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
47410.88 47433.8 -23701.44
Random effects:
Formula: ~1 | sigla_uf
(Intercept) Residual
StdDev: 529.9121 8072.302
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) -520.6662 328.7426 2270 -1.58381 0.1134
area_plantada 9.9110 0.0612 2270 161.90583 0.0000
Correlation:
(Intr)
area_plantada -0.297
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-8.495056130 -0.045866231 0.009214992 0.082672788 8.617890709
Number of Observations: 2276
Number of Groups: 5 LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - intercepto",
loglik = logLik(nlme_int_est)[1],
biblioteca = "nlme"
)
)
df_ll %>% filter(biblioteca == "nlme")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | nlme |
| municipios - inclinação | -19891.11 | nlme |
| municipios - intercepto e inclinação | -19891.11 | nlme |
| estados - intercepto | -23701.44 | nlme |
Ajuste do modelo
lme4_int_est <- lme4::lmer(
valor_producao ~ area_plantada + (1 | sigla_uf),
data = df_cafe
)Sumário do modelo
summary(lme4_int_est)Linear mixed model fit by REML ['lmerMod']
Formula: valor_producao ~ area_plantada + (1 | sigla_uf)
Data: df_cafe
REML criterion at convergence: 47402.9
Scaled residuals:
Min 1Q Median 3Q Max
-8.4951 -0.0459 0.0092 0.0827 8.6179
Random effects:
Groups Name Variance Std.Dev.
sigla_uf (Intercept) 280767 529.9
Residual 65162062 8072.3
Number of obs: 2276, groups: sigla_uf, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) -520.65426 328.72878 -1.584
area_plantada 9.91103 0.06121 161.906
Correlation of Fixed Effects:
(Intr)
area_plantd -0.297
fit warnings:
Some predictor variables are on very different scales: consider rescalingLogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - intercepto",
loglik = logLik(lme4_int_est)[1],
biblioteca = "lme4"
)
)
df_ll %>% filter(biblioteca == "lme4")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | lme4 |
| municipios - inclinação | -19891.11 | lme4 |
| municipios - intercepto e inclinação | -20132.16 | lme4 |
| estados - intercepto | -23701.44 | lme4 |
Ajuste do modelo
rstanarm_int_est <- rstanarm::stan_lmer(
valor_producao ~ area_plantada + (1 | sigla_uf),
data = df_cafe,
# esta linha é apenas para manter o exemplo pequeno!
chains = 2, cores = 4, seed = 12345, iter = 1000
)Sumário do modelo
rstanarm_int_eststan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (1 | sigla_uf)
observations: 2276
------
Median MAD_SD
(Intercept) -494.8 382.9
area_plantada 9.9 0.1
Auxiliary parameter(s):
Median MAD_SD
sigma 8072.8 127.3
Error terms:
Groups Name Std.Dev.
sigla_uf (Intercept) 1044
Residual 8074
Num. levels: sigla_uf 5
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregLogLik do modelo
loo_rstanarm_int_est <- rstanarm::loo(rstanarm_int_est)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - intercepto",
loglik = loo_rstanarm_int_est$estimates[1,1],
biblioteca = "rstanarm"
)
)
df_ll %>% filter(biblioteca == "rstanarm")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21255.82 | rstanarm |
| municipios - inclinação | -19013.34 | rstanarm |
| municipios - intercepto e inclinação | -19014.47 | rstanarm |
| estados - intercepto | -23720.41 | rstanarm |
Ajuste do modelo
brms_int_est <- brms::brm(
valor_producao ~ area_plantada + (1 | sigla_uf),
data = df_cafe,
# esta linha é apenas para manter o exemplo pequeno!
chains = 2, cores = 4, seed = 12345, iter = 1000
)Sumário do modelo
brms_int_est Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (1 | sigla_uf)
Data: df_cafe (Number of observations: 2276)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Group-Level Effects:
~sigla_uf (Number of levels: 5)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 719.54 427.48 94.12 1813.56 1.02 230 248
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -528.67 405.01 -1420.55 241.75 1.00 413 360
area_plantada 9.91 0.06 9.79 10.03 1.01 1242 764
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 8072.06 123.07 7839.32 8320.62 1.01 1031 806
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).LogLik do modelo
loo_brms_int_est <- brms::loo(brms_int_est)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - intercepto",
loglik = loo_brms_int_est$estimates[1,1],
biblioteca = "brms"
)
)
df_ll %>% filter(biblioteca == "brms")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21250.25 | brms |
| municipios - inclinação | -19011.98 | brms |
| municipios - intercepto e inclinação | -19034.87 | brms |
| estados - intercepto | -23719.77 | brms |
3.2.2.2 Modelo com inclinações aleatórias
Ajuste do modelo
nlme_slope_est <- nlme::lme(
valor_producao ~ area_plantada,
random = list(~ 0 + area_plantada | sigla_uf),
data = df_cafe
)Sumário do modelo
summary(nlme_slope_est)Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
47353.17 47376.09 -23672.59
Random effects:
Formula: ~0 + area_plantada | sigla_uf
area_plantada Residual
StdDev: 0.9265133 7952.216
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) -233.93797 189.0626 2270 -1.237357 0.2161
area_plantada 9.69433 0.4246 2270 22.831881 0.0000
Correlation:
(Intr)
area_plantada -0.075
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-8.69908387 -0.04987838 0.02418236 0.04234000 8.16877224
Number of Observations: 2276
Number of Groups: 5 LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - inclinação",
loglik = logLik(nlme_slope_est)[1],
biblioteca = "nlme"
)
)
df_ll %>% filter(biblioteca == "nlme")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | nlme |
| municipios - inclinação | -19891.11 | nlme |
| municipios - intercepto e inclinação | -19891.11 | nlme |
| estados - intercepto | -23701.44 | nlme |
| estados - inclinação | -23672.59 | nlme |
Ajuste do modelo
lme4_slope_est <- lme4::lmer(
valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf),
data = df_cafe
)Sumário do modelo
summary(lme4_slope_est)Linear mixed model fit by REML ['lmerMod']
Formula: valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf)
Data: df_cafe
REML criterion at convergence: 47345.2
Scaled residuals:
Min 1Q Median 3Q Max
-8.7011 -0.0500 0.0242 0.0423 8.1544
Random effects:
Groups Name Variance Std.Dev.
sigla_uf area_plantada 1.060e+00 1.029
Residual 6.322e+07 7950.927
Number of obs: 2276, groups: sigla_uf, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) -233.7417 189.0351 -1.236
area_plantada 9.6944 0.4697 20.640
Correlation of Fixed Effects:
(Intr)
area_plantd -0.067
fit warnings:
Some predictor variables are on very different scales: consider rescaling
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 1.9768 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - inclinação",
loglik = logLik(lme4_slope_est)[1],
biblioteca = "lme4"
)
)
df_ll %>% filter(biblioteca == "lme4")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | lme4 |
| municipios - inclinação | -19891.11 | lme4 |
| municipios - intercepto e inclinação | -20132.16 | lme4 |
| estados - intercepto | -23701.44 | lme4 |
| estados - inclinação | -23672.62 | lme4 |
Ajuste do modelo
rstanarm_slope_est <- rstanarm::stan_lmer(
valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf),
data = df_cafe,
# esta linha é apenas para manter o exemplo pequeno!
chains = 2, cores = 4, seed = 12345, iter = 1000
)Sumário do modelo
rstanarm_slope_eststan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf)
observations: 2276
------
Median MAD_SD
(Intercept) -229.9 178.4
area_plantada 9.7 0.5
Auxiliary parameter(s):
Median MAD_SD
sigma 7953.1 114.9
Error terms:
Groups Name Std.Dev.
sigla_uf area_plantada 1.5
Residual 7956.0
Num. levels: sigla_uf 5
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregLogLik do modelo
loo_rstanarm_slope_est <- rstanarm::loo(rstanarm_slope_est)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - inclinação",
loglik = loo_rstanarm_slope_est$estimates[1,1],
biblioteca = "rstanarm"
)
)
df_ll %>% filter(biblioteca == "rstanarm")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21255.82 | rstanarm |
| municipios - inclinação | -19013.34 | rstanarm |
| municipios - intercepto e inclinação | -19014.47 | rstanarm |
| estados - intercepto | -23720.41 | rstanarm |
| estados - inclinação | -23707.81 | rstanarm |
Ajuste do modelo
brms_slope_est <- brms::brm(
valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf),
data = df_cafe,
# esta linha é apenas para manter o exemplo pequeno!
chains = 2, cores = 4, seed = 12345, iter = 1000
)Sumário do modelo
brms_slope_est Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf)
Data: df_cafe (Number of observations: 2276)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Group-Level Effects:
~sigla_uf (Number of levels: 5)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(area_plantada) 1.95 1.67 0.57 6.99 1.04 44 28
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -253.51 194.62 -654.42 138.69 1.00 596 380
area_plantada 9.06 1.46 4.47 10.56 1.04 39 21
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 7954.09 113.67 7745.61 8191.49 1.01 406 286
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).LogLik do modelo
loo_brms_slope_est <- brms::loo(brms_slope_est)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - inclinação",
loglik = loo_brms_slope_est$estimates[1,1],
biblioteca = "brms"
)
)
df_ll %>% filter(biblioteca == "brms")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21250.25 | brms |
| municipios - inclinação | -19011.98 | brms |
| municipios - intercepto e inclinação | -19034.87 | brms |
| estados - intercepto | -23719.77 | brms |
| estados - inclinação | -23708.56 | brms |
3.2.2.3 Modelo com intercepto e inclinações aleatórias
Ajuste do modelo
nlme_int_slope_est <- nlme::lme(
valor_producao ~ area_plantada,
random = list(~ area_plantada | sigla_uf),
data = df_cafe
)Sumário do modelo
summary(nlme_int_slope_est)Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
47347.88 47382.26 -23667.94
Random effects:
Formula: ~area_plantada | sigla_uf
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 1384.940415 (Intr)
area_plantada 1.064845 -0.663
Residual 7923.204177
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) -845.2986 672.2123 2270 -1.257487 0.2087
area_plantada 9.7821 0.4862 2270 20.118614 0.0000
Correlation:
(Intr)
area_plantada -0.635
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-8.673989041 -0.065063553 0.001570653 0.045015337 8.070185710
Number of Observations: 2276
Number of Groups: 5 LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - intercepto e inclinação",
loglik = logLik(nlme_int_slope_est)[1],
biblioteca = "nlme"
)
)
df_ll %>% filter(biblioteca == "nlme")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | nlme |
| municipios - inclinação | -19891.11 | nlme |
| municipios - intercepto e inclinação | -19891.11 | nlme |
| estados - intercepto | -23701.44 | nlme |
| estados - inclinação | -23672.59 | nlme |
| estados - intercepto e inclinação | -23667.94 | nlme |
Ajuste do modelo
lme4_int_slope_est <- lme4::lmer(
valor_producao ~ area_plantada + (area_plantada | sigla_uf),
data = df_cafe
)Sumário do modelo
summary(lme4_int_slope_est)Linear mixed model fit by REML ['lmerMod']
Formula: valor_producao ~ area_plantada + (area_plantada | sigla_uf)
Data: df_cafe
REML criterion at convergence: 47346.4
Scaled residuals:
Min 1Q Median 3Q Max
-8.6841 -0.0666 -0.0017 0.0432 7.9872
Random effects:
Groups Name Variance Std.Dev. Corr
sigla_uf (Intercept) 3.426e+07 5852.782
area_plantada 3.845e+01 6.201 -0.97
Residual 6.262e+07 7913.033
Number of obs: 2276, groups: sigla_uf, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) -928.658 2631.252 -0.353
area_plantada 9.795 2.775 3.530
Correlation of Fixed Effects:
(Intr)
area_plantd -0.963
fit warnings:
Some predictor variables are on very different scales: consider rescaling
optimizer (nloptwrap) convergence code: 0 (OK)
unable to evaluate scaled gradient
Model failed to converge: degenerate Hessian with 1 negative eigenvaluesLogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - intercepto e inclinação",
loglik = logLik(lme4_int_slope_est)[1],
biblioteca = "lme4"
)
)
df_ll %>% filter(biblioteca == "lme4")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | lme4 |
| municipios - inclinação | -19891.11 | lme4 |
| municipios - intercepto e inclinação | -20132.16 | lme4 |
| estados - intercepto | -23701.44 | lme4 |
| estados - inclinação | -23672.62 | lme4 |
| estados - intercepto e inclinação | -23673.21 | lme4 |
Ajuste do modelo
rstanarm_int_slope_est <- rstanarm::stan_lmer(
valor_producao ~ area_plantada + (area_plantada | sigla_uf),
data = df_cafe,
# esta linha é apenas para manter o exemplo pequeno!
chains = 2, cores = 4, seed = 12345, iter = 1000
)Sumário do modelo
rstanarm_int_slope_eststan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (area_plantada | sigla_uf)
observations: 2276
------
Median MAD_SD
(Intercept) -248.8 206.3
area_plantada 9.7 0.5
Auxiliary parameter(s):
Median MAD_SD
sigma 7948.5 114.7
Error terms:
Groups Name Std.Dev. Corr
sigla_uf (Intercept) 286.2
area_plantada 1.6 -0.08
Residual 7948.7
Num. levels: sigla_uf 5
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregLogLik do modelo
loo_rstanarm_int_slope_est <- rstanarm::loo(rstanarm_int_slope_est)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - intercepto e inclinação",
loglik = loo_rstanarm_int_slope_est$estimates[1,1],
biblioteca = "rstanarm"
)
)
df_ll %>% filter(biblioteca == "rstanarm")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21255.82 | rstanarm |
| municipios - inclinação | -19013.34 | rstanarm |
| municipios - intercepto e inclinação | -19014.47 | rstanarm |
| estados - intercepto | -23720.41 | rstanarm |
| estados - inclinação | -23707.81 | rstanarm |
| estados - intercepto e inclinação | -23708.13 | rstanarm |
Ajuste do modelo
brms_int_slope_est <- brms::brm(
valor_producao ~ area_plantada + (area_plantada | sigla_uf),
data = df_cafe,
# esta linha é apenas para manter o exemplo pequeno!
chains = 2, cores = 4, seed = 12345, iter = 1000
)Sumário do modelo
brms_int_slope_est Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (area_plantada | sigla_uf)
Data: df_cafe (Number of observations: 2276)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Group-Level Effects:
~sigla_uf (Number of levels: 5)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept) 1776.34 986.23 484.10 4165.68 1.00 264
sd(area_plantada) 1.78 1.11 0.63 4.77 1.02 186
cor(Intercept,area_plantada) -0.31 0.43 -0.91 0.66 1.02 181
Tail_ESS
sd(Intercept) 350
sd(area_plantada) 112
cor(Intercept,area_plantada) 142
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1070.59 1160.22 -3582.53 802.27 1.01 315 327
area_plantada 9.41 1.15 6.28 11.34 1.02 118 70
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 7919.73 112.63 7705.72 8149.84 1.00 1056 667
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).LogLik do modelo
loo_brms_int_slope_est <- brms::loo(brms_int_slope_est)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados - intercepto e inclinação",
loglik = loo_brms_int_slope_est$estimates[1,1],
biblioteca = "brms"
)
)
df_ll %>% filter(biblioteca == "brms")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21250.25 | brms |
| municipios - inclinação | -19011.98 | brms |
| municipios - intercepto e inclinação | -19034.87 | brms |
| estados - intercepto | -23719.77 | brms |
| estados - inclinação | -23708.56 | brms |
| estados - intercepto e inclinação | -23702.82 | brms |
3.2.2.4 Comparação dos modelos para o nível de estados
performance::compare_performance(nlme_int_est,
nlme_slope_est,
nlme_int_slope_est, rank = TRUE)| Name | Model | R2_conditional | R2_marginal | ICC | RMSE | Sigma | AIC_wt | AICc_wt | BIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|
| nlme_int_slope_est | lme | 0.9259429 | 0.9133140 | 0.1456856 | 7907.524 | 7923.204 | 0.8171544 | 0.8157001 | 0.0143015 | 0.7518136 |
| nlme_slope_est | lme | 0.9240375 | 0.9133770 | 0.1230678 | 7942.026 | 7952.216 | 0.1828456 | 0.1842999 | 0.9856985 | 0.4852258 |
| nlme_int_est | lme | 0.9242470 | 0.9239205 | 0.0042909 | 8064.863 | 8072.302 | 0.0000000 | 0.0000000 | 0.0000000 | 0.1387412 |
# intercept vs slope
lmtest::lrtest(nlme_int_est, nlme_slope_est)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -23701.44 | NA | NA | NA |
| 4 | -23672.59 | 0 | 57.70796 | 0 |
# intercept vs intercept + slope
lmtest::lrtest(nlme_int_est, nlme_int_slope_est)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -23701.44 | NA | NA | NA |
| 6 | -23667.94 | 2 | 67.00051 | 0 |
# slope vs intercept + slope
lmtest::lrtest(nlme_slope_est, nlme_int_slope_est)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -23672.59 | NA | NA | NA |
| 6 | -23667.94 | 2 | 9.292549 | 0.0095973 |
performance::compare_performance(lme4_int_est,
lme4_slope_est,
lme4_int_slope_est, rank = TRUE)| Name | Model | R2_conditional | R2_marginal | ICC | RMSE | Sigma | AIC_wt | AICc_wt | BIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|
| lme4_slope_est | lmerMod | 0.9242501 | 0.9111208 | 0.1477210 | 7941.888 | 7950.927 | 0.996634 | 0.9966664 | 0.999989 | 0.7072082 |
| lme4_int_slope_est | lmerMod | 0.9462658 | 0.6661123 | 0.8390650 | 7905.832 | 7913.033 | 0.003366 | 0.0033336 | 0.000011 | 0.5008416 |
| lme4_int_est | lmerMod | 0.9242469 | 0.9239205 | 0.0042903 | 8064.863 | 8072.302 | 0.000000 | 0.0000000 | 0.000000 | 0.1250000 |
# intercept vs slope
lmtest::lrtest(lme4_int_est, lme4_slope_est)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -23701.44 | NA | NA | NA |
| 4 | -23672.62 | 0 | 57.63337 | 0 |
# intercept vs intercept + slope
lmtest::lrtest(lme4_int_est, lme4_int_slope_est)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -23701.44 | NA | NA | NA |
| 6 | -23673.21 | 2 | 56.45285 | 0 |
# slope vs intercept + slope
lmtest::lrtest(lme4_slope_est, lme4_int_slope_est)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -23672.62 | NA | NA | NA |
| 6 | -23673.21 | 2 | 1.180517 | 0.554184 |
performance::compare_performance(rstanarm_int_est,
rstanarm_slope_est,
rstanarm_int_slope_est, rank = TRUE)| Name | Model | R2 | R2_marginal | R2_adjusted | R2_adjusted_marginal | ICC | RMSE | Sigma | WAIC_wt | LOOIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|---|
| rstanarm_int_slope_est | stanreg | 0.9266089 | 0.9232193 | 0.9228432 | 0.9228432 | 0.2883107 | 7942.202 | 7948.457 | 0.7877877 | 0.4549389 | 0.6596446 |
| rstanarm_slope_est | stanreg | 0.9264159 | 0.9237139 | 0.9228559 | 0.9228559 | 0.2798354 | 7942.189 | 7953.085 | 0.2122045 | 0.0637836 | 0.5157594 |
| rstanarm_int_est | stanreg | 0.9241979 | 0.9242692 | 0.9233502 | 0.9233502 | 0.0164348 | 8065.017 | 8072.831 | 0.0000078 | 0.4812774 | 0.4444444 |
rstanarm::loo_compare(loo_rstanarm_int_est,
loo_rstanarm_slope_est,
loo_rstanarm_int_slope_est) elpd_diff se_diff
rstanarm_slope_est 0.0 0.0
rstanarm_int_slope_est -0.3 1.3
rstanarm_int_est -12.6 30.0 performance::compare_performance(brms_int_est,
brms_slope_est,
brms_int_slope_est, rank = TRUE)| Name | Model | R2 | R2_marginal | R2_adjusted | R2_adjusted_marginal | ICC | RMSE | Sigma | WAIC_wt | LOOIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|---|
| brms_int_slope_est | brmsfit | 0.9269569 | 0.9188242 | 0.9233108 | 0.5775771 | 0.4304307 | 7912.624 | 7705.927 | 0.8946403 | 0.5572164 | 0.8268800 |
| brms_int_est | brmsfit | 0.9241783 | 0.9240278 | 0.9231298 | 0.9236083 | 0.0078832 | 8059.369 | 7930.657 | 0.0000012 | 0.4427759 | 0.3914289 |
| brms_slope_est | brmsfit | 0.9264615 | 0.9158918 | 0.9226446 | 0.5468738 | 0.3842235 | 7960.086 | 7774.344 | 0.1053584 | 0.0000077 | 0.3558065 |
brms::loo_compare(
brms::add_criterion(brms_int_est,"loo"),
brms::add_criterion(brms_slope_est,"loo"),
brms::add_criterion(brms_int_slope_est,"loo")
) elpd_diff se_diff
brms::add_criterion(brms_int_slope_est, "loo") 0.0 0.0
brms::add_criterion(brms_slope_est, "loo") -5.7 5.4
brms::add_criterion(brms_int_est, "loo") -16.9 32.1 3.2.2.5 Sumarização dos resultados para nível de estados
Diferente do que ocorreu para o nível de municípios, para o nível de estados os modelos com intercepto, inclinação e intercepto+inclinação apresentaram desempenho semelhante. O que corrobora também, com a nossa hipótese de que, o efeito aleatório de estados é menos relevante do que o efeito aleatório de municípios, visto na interpretação qualitativa da análise exploratória.
Agora vamos avaliar se modelos com os dois níveis de efeito aleatório (estados e municípios) se ajustam melhor aos dados.
3.2.3 MUNICÍPIOS POR ESTADOS (aninhados)
3.2.3.1 Modelo com interceptos aleatórios
Ajuste do modelo
nlme_int_est_mun <- nlme::lme(
valor_producao ~ area_plantada,
random = list(~ 1 | sigla_uf, ~ 1 | id_municipio),
data = df_cafe
)Sumário do modelo
summary(nlme_int_est_mun)Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
45261.36 45290.01 -22625.68
Random effects:
Formula: ~1 | sigla_uf
(Intercept)
StdDev: 1104.33
Formula: ~1 | id_municipio %in% sigla_uf
(Intercept) Residual
StdDev: 8036.572 2294.035
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) 789.5451 598.4477 1183 1.31932 0.1873
area_plantada 9.1904 0.0665 1183 138.26113 0.0000
Correlation:
(Intr)
area_plantada -0.184
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-8.658080695 -0.024168104 -0.009842888 0.005329907 11.607880218
Number of Observations: 2276
Number of Groups:
sigla_uf id_municipio %in% sigla_uf
5 1092 LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - intercepto",
loglik = logLik(nlme_int_est_mun)[1],
biblioteca = "nlme"
)
)
df_ll %>% filter(biblioteca == "nlme")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | nlme |
| municipios - inclinação | -19891.11 | nlme |
| municipios - intercepto e inclinação | -19891.11 | nlme |
| estados - intercepto | -23701.44 | nlme |
| estados - inclinação | -23672.59 | nlme |
| estados - intercepto e inclinação | -23667.94 | nlme |
| estados e municípios - intercepto | -22625.68 | nlme |
Ajuste do modelo
lme4_int_est_mun <- lme4::lmer(
valor_producao ~ area_plantada + (1 | sigla_uf/id_municipio),
data = df_cafe
)Sumário do modelo
summary(lme4_int_est_mun)Linear mixed model fit by REML ['lmerMod']
Formula: valor_producao ~ area_plantada + (1 | sigla_uf/id_municipio)
Data: df_cafe
REML criterion at convergence: 45251.4
Scaled residuals:
Min 1Q Median 3Q Max
-8.6581 -0.0242 -0.0098 0.0053 11.6079
Random effects:
Groups Name Variance Std.Dev.
id_municipio:sigla_uf (Intercept) 64586494 8037
sigla_uf (Intercept) 1219545 1104
Residual 5262596 2294
Number of obs: 2276, groups: id_municipio:sigla_uf, 1092; sigla_uf, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 789.54515 598.44779 1.319
area_plantada 9.19039 0.06647 138.261
Correlation of Fixed Effects:
(Intr)
area_plantd -0.184
fit warnings:
Some predictor variables are on very different scales: consider rescalingLogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - intercepto",
loglik = logLik(lme4_int_est_mun)[1],
biblioteca = "lme4"
)
)
df_ll %>% filter(biblioteca == "lme4")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | lme4 |
| municipios - inclinação | -19891.11 | lme4 |
| municipios - intercepto e inclinação | -20132.16 | lme4 |
| estados - intercepto | -23701.44 | lme4 |
| estados - inclinação | -23672.62 | lme4 |
| estados - intercepto e inclinação | -23673.21 | lme4 |
| estados e municípios - intercepto | -22625.68 | lme4 |
Ajuste do modelo
rstanarm_int_est_mun <- rstanarm::stan_lmer(
valor_producao ~ area_plantada + (1 | sigla_uf/id_municipio),
data = df_cafe
)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000141 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.41 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 32.082 seconds (Warm-up)
Chain 1: 38.777 seconds (Sampling)
Chain 1: 70.859 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 7.8e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.78 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 28.903 seconds (Warm-up)
Chain 2: 38.111 seconds (Sampling)
Chain 2: 67.014 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 7.9e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.79 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 228.813 seconds (Warm-up)
Chain 3: 19.558 seconds (Sampling)
Chain 3: 248.371 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 8.2e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.82 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 30.499 seconds (Warm-up)
Chain 4: 19.625 seconds (Sampling)
Chain 4: 50.124 seconds (Total)
Chain 4: Sumário do modelo
rstanarm_int_est_munstan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (1 | sigla_uf/id_municipio)
observations: 2276
------
Median MAD_SD
(Intercept) 772.3 620.9
area_plantada 9.2 0.1
Auxiliary parameter(s):
Median MAD_SD
sigma 2299.3 49.2
Error terms:
Groups Name Std.Dev.
id_municipio:sigla_uf (Intercept) 8027
sigla_uf (Intercept) 1563
Residual 2300
Num. levels: id_municipio:sigla_uf 1092, sigla_uf 5
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregLogLik do modelo
loo_rstanarm_int_est_mun <- rstanarm::loo(rstanarm_int_est_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - intercepto",
loglik = loo_rstanarm_int_est_mun$estimates[1,1],
biblioteca = "rstanarm"
)
)
df_ll %>% filter(biblioteca == "rstanarm")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21255.82 | rstanarm |
| municipios - inclinação | -19013.34 | rstanarm |
| municipios - intercepto e inclinação | -19014.47 | rstanarm |
| estados - intercepto | -23720.41 | rstanarm |
| estados - inclinação | -23707.81 | rstanarm |
| estados - intercepto e inclinação | -23708.13 | rstanarm |
| estados e municípios - intercepto | -21286.08 | rstanarm |
Ajuste do modelo
brms_int_est_mun <- brms::brm(
valor_producao ~ area_plantada + (1 | sigla_uf/id_municipio),
data = df_cafe
)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000131 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.31 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 16.663 seconds (Warm-up)
Chain 1: 18.724 seconds (Sampling)
Chain 1: 35.387 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 7.9e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.79 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 24.786 seconds (Warm-up)
Chain 2: 18.552 seconds (Sampling)
Chain 2: 43.338 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 7.1e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.71 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 14.923 seconds (Warm-up)
Chain 3: 9.345 seconds (Sampling)
Chain 3: 24.268 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 7e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.7 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 19.036 seconds (Warm-up)
Chain 4: 9.32 seconds (Sampling)
Chain 4: 28.356 seconds (Total)
Chain 4: Sumário do modelo
brms_int_est_mun Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (1 | sigla_uf/id_municipio)
Data: df_cafe (Number of observations: 2276)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~sigla_uf (Number of levels: 5)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1439.28 826.62 350.33 3355.73 1.02 216 354
~sigla_uf:id_municipio (Number of levels: 1092)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 8035.61 184.37 7680.18 8401.81 1.01 473 874
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 613.80 887.89 -1251.45 2172.10 1.00 812 1027
area_plantada 9.19 0.07 9.05 9.33 1.01 350 885
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 2293.81 48.95 2202.42 2395.27 1.00 1392 2515
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).LogLik do modelo
loo_brms_int_est_mun <- brms::loo(brms_int_est_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - intercepto",
loglik = loo_brms_int_est_mun$estimates[1,1],
biblioteca = "brms"
)
)
df_ll %>% filter(biblioteca == "brms")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21250.25 | brms |
| municipios - inclinação | -19011.98 | brms |
| municipios - intercepto e inclinação | -19034.87 | brms |
| estados - intercepto | -23719.77 | brms |
| estados - inclinação | -23708.56 | brms |
| estados - intercepto e inclinação | -23702.82 | brms |
| estados e municípios - intercepto | -21276.93 | brms |
3.2.3.2 Modelo com inclinações aleatórias
Ajuste do modelo
nlme_slope_est_mun <- nlme::lme(
valor_producao ~ area_plantada,
random = list(~ 0 + area_plantada | sigla_uf, ~ 0 + area_plantada | id_municipio),
data = df_cafe
)Sumário do modelo
summary(nlme_slope_est_mun)Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
39764.01 39792.65 -19877
Random effects:
Formula: ~0 + area_plantada | sigla_uf
area_plantada
StdDev: 1.072088
Formula: ~0 + area_plantada | id_municipio %in% sigla_uf
area_plantada Residual
StdDev: 2.866569 919.7565
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) -61.00413 28.603700 1183 -2.132736 0.0332
area_plantada 8.91062 0.508967 1183 17.507280 0.0000
Correlation:
(Intr)
area_plantada -0.103
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-16.76832956 -0.05747788 0.01229872 0.06098704 13.71747550
Number of Observations: 2276
Number of Groups:
sigla_uf id_municipio %in% sigla_uf
5 1092 LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - inclinação",
loglik = logLik(nlme_slope_est_mun)[1],
biblioteca = "nlme"
)
)
df_ll %>% filter(biblioteca == "nlme")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | nlme |
| municipios - inclinação | -19891.11 | nlme |
| municipios - intercepto e inclinação | -19891.11 | nlme |
| estados - intercepto | -23701.44 | nlme |
| estados - inclinação | -23672.59 | nlme |
| estados - intercepto e inclinação | -23667.94 | nlme |
| estados e municípios - intercepto | -22625.68 | nlme |
| estados e municípios - inclinação | -19877.00 | nlme |
Ajuste do modelo
lme4_slope_est_mun <- lme4::lmer(
valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf/id_municipio),
data = df_cafe
)Sumário do modelo
summary(lme4_slope_est_mun)Linear mixed model fit by REML ['lmerMod']
Formula:
valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf/id_municipio)
Data: df_cafe
REML criterion at convergence: 39801.8
Scaled residuals:
Min 1Q Median 3Q Max
-16.7956 -0.0596 0.0113 0.0602 13.7391
Random effects:
Groups Name Variance Std.Dev.
id_municipio:sigla_uf area_plantada 8.241e+00 2.871
sigla_uf area_plantada 5.465e+05 739.262
Residual 8.433e+05 918.295
Number of obs: 2276, groups: id_municipio:sigla_uf, 1092; sigla_uf, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) -60.192 28.595 -2.105
area_plantada 8.845 330.609 0.027
Correlation of Fixed Effects:
(Intr)
area_plantd 0.000
fit warnings:
Some predictor variables are on very different scales: consider rescaling
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.884121 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?
Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - inclinação",
loglik = logLik(lme4_slope_est_mun)[1],
biblioteca = "lme4"
)
)
df_ll %>% filter(biblioteca == "lme4")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | lme4 |
| municipios - inclinação | -19891.11 | lme4 |
| municipios - intercepto e inclinação | -20132.16 | lme4 |
| estados - intercepto | -23701.44 | lme4 |
| estados - inclinação | -23672.62 | lme4 |
| estados - intercepto e inclinação | -23673.21 | lme4 |
| estados e municípios - intercepto | -22625.68 | lme4 |
| estados e municípios - inclinação | -19900.90 | lme4 |
Ajuste do modelo
rstanarm_slope_est_mun <- rstanarm::stan_lmer(
valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf/id_municipio),
data = df_cafe
)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000233 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 655.222 seconds (Warm-up)
Chain 1: 211.186 seconds (Sampling)
Chain 1: 866.408 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000226 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 2.26 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 1593.2 seconds (Warm-up)
Chain 2: 210.566 seconds (Sampling)
Chain 2: 1803.77 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000212 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.12 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 747.054 seconds (Warm-up)
Chain 3: 212.82 seconds (Sampling)
Chain 3: 959.874 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000227 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.27 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 713.266 seconds (Warm-up)
Chain 4: 412.875 seconds (Sampling)
Chain 4: 1126.14 seconds (Total)
Chain 4: Sumário do modelo
rstanarm_slope_est_munstan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf/id_municipio)
observations: 2276
------
Median MAD_SD
(Intercept) -60.7 28.4
area_plantada 8.9 0.7
Auxiliary parameter(s):
Median MAD_SD
sigma 919.4 15.6
Error terms:
Groups Name Std.Dev.
id_municipio:sigla_uf area_plantada 2.9
sigla_uf area_plantada 2.0
Residual 919.7
Num. levels: id_municipio:sigla_uf 1092, sigla_uf 5
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregLogLik do modelo
loo_rstanarm_slope_est_mun <- rstanarm::loo(rstanarm_slope_est_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - inclinação",
loglik = loo_rstanarm_slope_est_mun$estimates[1,1],
biblioteca = "rstanarm"
)
)
df_ll %>% filter(biblioteca == "rstanarm")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21255.82 | rstanarm |
| municipios - inclinação | -19013.34 | rstanarm |
| municipios - intercepto e inclinação | -19014.47 | rstanarm |
| estados - intercepto | -23720.41 | rstanarm |
| estados - inclinação | -23707.81 | rstanarm |
| estados - intercepto e inclinação | -23708.13 | rstanarm |
| estados e municípios - intercepto | -21286.08 | rstanarm |
| estados e municípios - inclinação | -19028.02 | rstanarm |
Ajuste do modelo
brms_slope_est_mun <- brms::brm(
valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf/id_municipio),
data = df_cafe
)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000214 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.14 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 79.03 seconds (Warm-up)
Chain 1: 90.026 seconds (Sampling)
Chain 1: 169.056 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 8.9e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.89 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 77.932 seconds (Warm-up)
Chain 2: 90.384 seconds (Sampling)
Chain 2: 168.316 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 9.1e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.91 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 78.68 seconds (Warm-up)
Chain 3: 89.998 seconds (Sampling)
Chain 3: 168.678 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 9.5e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.95 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 78.272 seconds (Warm-up)
Chain 4: 90.129 seconds (Sampling)
Chain 4: 168.401 seconds (Total)
Chain 4: Sumário do modelo
brms_slope_est_mun Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (0 + area_plantada | sigla_uf/id_municipio)
Data: df_cafe (Number of observations: 2276)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~sigla_uf (Number of levels: 5)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(area_plantada) 4.12 4.32 0.67 15.67 1.35 10 75
~sigla_uf:id_municipio (Number of levels: 1092)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(area_plantada) 2.88 0.09 2.72 3.07 1.08 65 34
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -57.98 28.75 -114.52 -1.22 1.02 279 2445
area_plantada 6.41 3.07 -0.11 9.56 1.42 9 32
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 919.70 16.18 889.51 952.34 1.01 915 1852
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).LogLik do modelo
loo_brms_slope_est_mun <- brms::loo(brms_slope_est_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - inclinação",
loglik = loo_brms_slope_est_mun$estimates[1,1],
biblioteca = "brms"
)
)
df_ll %>% filter(biblioteca == "brms")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21250.25 | brms |
| municipios - inclinação | -19011.98 | brms |
| municipios - intercepto e inclinação | -19034.87 | brms |
| estados - intercepto | -23719.77 | brms |
| estados - inclinação | -23708.56 | brms |
| estados - intercepto e inclinação | -23702.82 | brms |
| estados e municípios - intercepto | -21276.93 | brms |
| estados e municípios - inclinação | -19028.17 | brms |
3.2.3.3 Modelo com intercepto e inclinações aleatórias
Ajuste do modelo
nlme_int_slope_est_mun <- nlme::lme(
valor_producao ~ area_plantada,
random = list(~ area_plantada | sigla_uf, ~ area_plantada | id_municipio),
data = df_cafe
)Sumário do modelo
summary(nlme_int_slope_est_mun)Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
39769.95 39821.52 -19875.98
Random effects:
Formula: ~area_plantada | sigla_uf
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 57.8087428 (Intr)
area_plantada 0.9277332 0.998
Formula: ~area_plantada | id_municipio %in% sigla_uf
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 7.309286e-04 (Intr)
area_plantada 2.859809e+00 0.036
Residual 9.198860e+02
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) -94.53910 39.30005 1183 -2.405572 0.0163
area_plantada 8.97797 0.44543 1183 20.155718 0.0000
Correlation:
(Intr)
area_plantada 0.555
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-16.751955513 -0.065936387 0.007153838 0.062897106 13.719801256
Number of Observations: 2276
Number of Groups:
sigla_uf id_municipio %in% sigla_uf
5 1092 LogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - intercepto e inclinação",
loglik = logLik(nlme_int_slope_est_mun)[1],
biblioteca = "nlme"
)
)
df_ll %>% filter(biblioteca == "nlme")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | nlme |
| municipios - inclinação | -19891.11 | nlme |
| municipios - intercepto e inclinação | -19891.11 | nlme |
| estados - intercepto | -23701.44 | nlme |
| estados - inclinação | -23672.59 | nlme |
| estados - intercepto e inclinação | -23667.94 | nlme |
| estados e municípios - intercepto | -22625.68 | nlme |
| estados e municípios - inclinação | -19877.00 | nlme |
| estados e municípios - intercepto e inclinação | -19875.98 | nlme |
Ajuste do modelo
lme4_int_slope_est_mun <- lme4::lmer(
valor_producao ~ area_plantada + (area_plantada | sigla_uf/id_municipio),
data = df_cafe
)Sumário do modelo
summary(lme4_int_slope_est_mun)Linear mixed model fit by REML ['lmerMod']
Formula:
valor_producao ~ area_plantada + (area_plantada | sigla_uf/id_municipio)
Data: df_cafe
REML criterion at convergence: 40302.4
Scaled residuals:
Min 1Q Median 3Q Max
-15.3636 -0.0346 0.0012 0.0391 15.1178
Random effects:
Groups Name Variance Std.Dev. Corr
id_municipio:sigla_uf (Intercept) 894097.65 945.567
area_plantada 10.91 3.303 -0.57
sigla_uf (Intercept) 463261.63 680.633
area_plantada 751208.97 866.723 0.11
Residual 759530.29 871.510
Number of obs: 2276, groups: id_municipio:sigla_uf, 1092; sigla_uf, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) -200.793 311.704 -0.644
area_plantada 8.926 387.609 0.023
Correlation of Fixed Effects:
(Intr)
area_plantd 0.104
fit warnings:
Some predictor variables are on very different scales: consider rescaling
optimizer (nloptwrap) convergence code: 0 (OK)
unable to evaluate scaled gradient
Model failed to converge: degenerate Hessian with 2 negative eigenvaluesLogLik do modelo
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - intercepto e inclinação",
loglik = logLik(lme4_int_slope_est_mun)[1],
biblioteca = "lme4"
)
)
df_ll %>% filter(biblioteca == "lme4")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -22627.91 | lme4 |
| municipios - inclinação | -19891.11 | lme4 |
| municipios - intercepto e inclinação | -20132.16 | lme4 |
| estados - intercepto | -23701.44 | lme4 |
| estados - inclinação | -23672.62 | lme4 |
| estados - intercepto e inclinação | -23673.21 | lme4 |
| estados e municípios - intercepto | -22625.68 | lme4 |
| estados e municípios - inclinação | -19900.90 | lme4 |
| estados e municípios - intercepto e inclinação | -20151.21 | lme4 |
Ajuste do modelo
rstanarm_int_slope_est_mun <- rstanarm::stan_lmer(
valor_producao ~ area_plantada + (area_plantada | sigla_uf/id_municipio),
data = df_cafe
)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.00057 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.7 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 4612.36 seconds (Warm-up)
Chain 1: 1779.01 seconds (Sampling)
Chain 1: 6391.37 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000505 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 5.05 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 4705.05 seconds (Warm-up)
Chain 2: 1946.43 seconds (Sampling)
Chain 2: 6651.49 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.0005 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 5 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 4890.56 seconds (Warm-up)
Chain 3: 2295.11 seconds (Sampling)
Chain 3: 7185.68 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000484 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 4.84 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 3732.84 seconds (Warm-up)
Chain 4: 1904.58 seconds (Sampling)
Chain 4: 5637.41 seconds (Total)
Chain 4: Sumário do modelo
rstanarm_int_slope_est_munstan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (area_plantada | sigla_uf/id_municipio)
observations: 2276
------
Median MAD_SD
(Intercept) -61.9 28.9
area_plantada 8.9 0.7
Auxiliary parameter(s):
Median MAD_SD
sigma 919.4 16.0
Error terms:
Groups Name Std.Dev. Corr
id_municipio:sigla_uf (Intercept) 37.4
area_plantada 2.9 -0.42
sigla_uf (Intercept) 18.3
area_plantada 2.0 0.04
Residual 919.7
Num. levels: id_municipio:sigla_uf 1092, sigla_uf 5
------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanregLogLik do modelo
loo_rstanarm_int_slope_est_mun <- rstanarm::loo(rstanarm_int_slope_est_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - intercepto e inclinação",
loglik = loo_rstanarm_int_slope_est_mun$estimates[1,1],
biblioteca = "rstanarm"
)
)
df_ll %>% filter(biblioteca == "rstanarm")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21255.82 | rstanarm |
| municipios - inclinação | -19013.34 | rstanarm |
| municipios - intercepto e inclinação | -19014.47 | rstanarm |
| estados - intercepto | -23720.41 | rstanarm |
| estados - inclinação | -23707.81 | rstanarm |
| estados - intercepto e inclinação | -23708.13 | rstanarm |
| estados e municípios - intercepto | -21286.08 | rstanarm |
| estados e municípios - inclinação | -19028.02 | rstanarm |
| estados e municípios - intercepto e inclinação | -19029.49 | rstanarm |
Ajuste do modelo
brms_int_slope_est_mun <- brms::brm(
valor_producao ~ area_plantada + (area_plantada | sigla_uf/id_municipio),
data = df_cafe
)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000229 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.29 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 131.914 seconds (Warm-up)
Chain 1: 146.801 seconds (Sampling)
Chain 1: 278.715 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000161 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.61 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 130.646 seconds (Warm-up)
Chain 2: 144.426 seconds (Sampling)
Chain 2: 275.072 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000152 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.52 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 129.642 seconds (Warm-up)
Chain 3: 148.246 seconds (Sampling)
Chain 3: 277.888 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000153 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.53 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 128.764 seconds (Warm-up)
Chain 4: 146.174 seconds (Sampling)
Chain 4: 274.938 seconds (Total)
Chain 4: Sumário do modelo
brms_int_slope_est_mun Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (area_plantada | sigla_uf/id_municipio)
Data: df_cafe (Number of observations: 2276)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~sigla_uf (Number of levels: 5)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept) 177.03 173.10 2.19 632.02 1.04 59
sd(area_plantada) 12.15 5.39 5.86 26.63 1.08 44
cor(Intercept,area_plantada) 0.22 0.60 -0.93 0.99 1.16 19
Tail_ESS
sd(Intercept) 39
sd(area_plantada) 89
cor(Intercept,area_plantada) 66
~sigla_uf:id_municipio (Number of levels: 1092)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept) 49.56 43.38 0.36 130.43 1.80 6
sd(area_plantada) 2.91 0.11 2.72 3.14 1.17 20
cor(Intercept,area_plantada) -0.13 0.69 -0.87 0.82 2.96 5
Tail_ESS
sd(Intercept) 42
sd(area_plantada) 83
cor(Intercept,area_plantada) 15
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -128.55 124.83 -439.25 87.29 1.08 38 79
area_plantada 0.30 0.29 0.00 1.28 1.61 7 13
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 922.28 15.41 892.74 952.70 1.03 73 101
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).LogLik do modelo
loo_brms_int_slope_est_mun <- brms::loo(brms_int_slope_est_mun)
df_ll <- rbind(
df_ll,
tibble(
modelo = "estados e municípios - intercepto e inclinação",
loglik = loo_brms_int_slope_est_mun$estimates[1,1],
biblioteca = "brms"
)
)
df_ll %>% filter(biblioteca == "brms")| modelo | loglik | biblioteca |
|---|---|---|
| municipios - intercepto | -21250.25 | brms |
| municipios - inclinação | -19011.98 | brms |
| municipios - intercepto e inclinação | -19034.87 | brms |
| estados - intercepto | -23719.77 | brms |
| estados - inclinação | -23708.56 | brms |
| estados - intercepto e inclinação | -23702.82 | brms |
| estados e municípios - intercepto | -21276.93 | brms |
| estados e municípios - inclinação | -19028.17 | brms |
| estados e municípios - intercepto e inclinação | -19035.04 | brms |
3.2.3.4 Comparação dos modelos para o nível de estados e municípios
performance::compare_performance(nlme_int_est_mun,
nlme_slope_est_mun,
nlme_int_slope_est_mun, rank = TRUE)| Name | Model | RMSE | Sigma | AIC_wt | AICc_wt | BIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|
| nlme_slope_est_mun | lme | 806.4771 | 919.7565 | 0.9516232 | 0.9528287 | 0.9999995 | 1.0000000 |
| nlme_int_slope_est_mun | lme | 806.5823 | 919.8860 | 0.0483768 | 0.0471713 | 0.0000005 | 0.4200258 |
| nlme_int_est_mun | lme | 1682.9572 | 2294.0350 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 |
# intercept vs slope
lmtest::lrtest(nlme_int_est_mun, nlme_slope_est_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 5 | -22625.68 | NA | NA | NA |
| 5 | -19877.00 | 0 | 5497.356 | 0 |
# intercept vs intercept + slope
lmtest::lrtest(nlme_int_est_mun, nlme_int_slope_est_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 5 | -22625.68 | NA | NA | NA |
| 9 | -19875.98 | 4 | 5499.409 | 0 |
# slope vs intercept + slope
lmtest::lrtest(nlme_slope_est_mun, nlme_int_slope_est_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 5 | -19877.00 | NA | NA | NA |
| 9 | -19875.98 | 4 | 2.052831 | 0.7260423 |
performance::compare_performance(lme4_int_est_mun,
lme4_slope_est_mun,
lme4_int_slope_est_mun, rank = TRUE)| Name | Model | R2_conditional | R2_marginal | ICC | RMSE | Sigma | AIC_wt | AICc_wt | BIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|
| lme4_slope_est_mun | lmerMod | 0.9999999 | 0.0001120 | 0.9999999 | 806.1607 | 918.2947 | 1 | 1 | 1 | 0.8641952 |
| lme4_int_slope_est_mun | lmerMod | 0.9999999 | 0.0000830 | 0.9999999 | 756.5281 | 871.5103 | 0 | 0 | 0 | 0.5000000 |
| lme4_int_est_mun | lmerMod | 0.9930245 | 0.9058001 | 0.9259505 | 1682.9571 | 2294.0349 | 0 | 0 | 0 | 0.1250000 |
# intercept vs slope
lmtest::lrtest(lme4_int_est_mun, lme4_slope_est_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 5 | -22625.68 | NA | NA | NA |
| 5 | -19900.90 | 0 | 5449.573 | 0 |
# intercept vs intercept + slope
lmtest::lrtest(lme4_int_est_mun, lme4_int_slope_est_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 5 | -22625.68 | NA | NA | NA |
| 9 | -20151.21 | 4 | 4948.944 | 0 |
# slope vs intercept + slope
lmtest::lrtest(lme4_slope_est_mun, lme4_int_slope_est_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 5 | -19900.90 | NA | NA | NA |
| 9 | -20151.21 | 4 | 500.6291 | 0 |
performance::compare_performance(rstanarm_int_est_mun,
rstanarm_slope_est_mun,
rstanarm_int_slope_est_mun, rank = TRUE)| Name | Model | R2 | R2_marginal | R2_adjusted | R2_adjusted_marginal | ICC | RMSE | Sigma | WAIC_wt | LOOIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|---|
| rstanarm_slope_est_mun | stanreg | 0.9990165 | 0.9986835 | 0.9990052 | 0.9990052 | 0.9933007 | 808.1812 | 919.3872 | 0.8208663 | 0.9884992 | 0.9995313 |
| rstanarm_int_slope_est_mun | stanreg | 0.9990166 | 0.9986830 | 0.9990006 | 0.9990006 | 0.9934342 | 806.2692 | 919.3721 | 0.1791337 | 0.0001125 | 0.8018060 |
| rstanarm_int_est_mun | stanreg | 0.9938468 | 0.9923321 | 0.9941370 | 0.9941370 | 0.9266862 | 1685.8189 | 2299.3240 | 0.0000000 | 0.0113884 | 0.0012676 |
rstanarm::loo_compare(loo_rstanarm_int_est_mun,
loo_rstanarm_slope_est_mun,
loo_rstanarm_int_slope_est_mun) elpd_diff se_diff
rstanarm_slope_est_mun 0.0 0.0
rstanarm_int_slope_est_mun -1.5 3.1
rstanarm_int_est_mun -2258.1 254.9performance::compare_performance(brms_int_est_mun,
brms_slope_est_mun,
brms_int_slope_est_mun, rank = TRUE)| Name | Model | R2 | R2_marginal | R2_adjusted | R2_adjusted_marginal | ICC | RMSE | Sigma | WAIC_wt | LOOIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|---|
| brms_slope_est_mun | brmsfit | 0.9990168 | 0.8485788 | 0.9990204 | -0.8245032 | 0.9967655 | 806.9336 | 928.1870 | 0.9999991 | 0.9886731 | 0.8768554 |
| brms_int_slope_est_mun | brmsfit | 0.9990156 | 0.0003732 | 0.9989906 | -0.2567290 | 0.9995563 | 806.0623 | 922.9941 | 0.0000009 | 0.0000006 | 0.5924050 |
| brms_int_est_mun | brmsfit | 0.9938860 | 0.9076653 | 0.9941861 | 0.8550837 | 0.9268255 | 1685.7517 | 2284.0660 | 0.0000000 | 0.0113263 | 0.2234950 |
brms::loo_compare(
brms::add_criterion(brms_int_est_mun,"loo"),
brms::add_criterion(brms_slope_est_mun,"loo"),
brms::add_criterion(brms_int_slope_est_mun,"loo")
) elpd_diff se_diff
brms::add_criterion(brms_slope_est_mun, "loo") 0.0 0.0
brms::add_criterion(brms_int_slope_est_mun, "loo") -6.9 6.2
brms::add_criterion(brms_int_est_mun, "loo") -2248.8 254.73.2.3.5 Sumarização dos resultados para nível de estados e municípios
Os modelos que consideraram dois níveis hierárquicos (estados e municípios) tiveram o mesmo comportamento dos modelos que consideraram apenas o nível de municípios. O que sugere que o efeito aleatório do nível de municípios esteja explicando a maior parte da variabilidade dos dados dentro de uma estrutura hierárquica.
4 Batalha Final – Seleção dos melhores modelos estudados
Ao abordar a seleção de modelos hierárquicos, a escolha da menor log-verossimilhança (loglik) é pensada como uma abordagem fundamental e robusta. Os modelos hierárquicos, que incorporam a estrutura de agrupamento e aninhamento de dados, têm grande destaque devido à sua capacidade de lidar com a complexidade inerente aos conjuntos de dados observados em diversas áreas de estudo. Nesse contexto, a seleção do modelo mais adequado torna-se crucial para garantir uma inferência precisa e eficaz.
A menor log-verossimilhança, como critério de seleção, representa uma estratégia estatística que busca maximizar a adequação do modelo aos dados observados, ao mesmo tempo em que penaliza a complexidade do modelo. Essa abordagem é fundamentada no princípio da parcimônia, ou seja, a preferência por modelos mais simples quando estes são igualmente capazes de explicar os dados.
Ao escolher o modelo hierárquico com base na menor log-verossimilhança, estamos buscando um equilíbrio entre o ajuste do modelo aos dados e a sua complexidade. Modelos excessivamente complexos podem levar a uma superestimação dos parâmetros e, consequentemente, a uma menor capacidade de generalização para novos conjuntos de dados. Por outro lado, modelos muito simples podem não capturar adequadamente a estrutura subjacente nos dados, resultando em uma perda de informações importantes.
A utilização da menor log-verossimilhança como critério de seleção também oferece vantagens práticas. Ela é amplamente aplicável e fácil de interpretar, proporcionando uma base sólida para a comparação entre diferentes modelos hierárquicos. Além disso, sua formulação matemática permite a utilização de métodos de inferência estatística bem estabelecidos, facilitando a análise e interpretação dos resultados.
Por fim, a escolha da menor log-verossimilhança para selecionar os melhores modelos hierárquicos é justificada não apenas por sua fundamentação teórica sólida, mas também por sua aplicabilidade prática e capacidade de equilibrar o ajuste do modelo aos dados com sua complexidade. Essa abordagem representa um passo significativo em direção a uma inferência mais precisa e confiável em contextos de modelagem hierárquica.
Vamos ver os top 3 de cada biblioteca, considerando a menor log-verossimilhança (loglik).
df_ll %>%
filter(biblioteca == "nlme") %>%
group_by(biblioteca) %>%
slice_max(order_by = loglik, n=3) %>%
mutate(rank = row_number()) %>%
ungroup() | modelo | loglik | biblioteca | rank |
|---|---|---|---|
| estados e municípios - intercepto e inclinação | -19875.98 | nlme | 1 |
| estados e municípios - inclinação | -19877.00 | nlme | 2 |
| municipios - inclinação | -19891.11 | nlme | 3 |
performance::compare_performance(nlme_int_slope_est_mun,
nlme_slope_est_mun,
nlme_slope_mun, rank = TRUE)Random effect variances not available. Returned R2 does not account for random effects.| Name | Model | RMSE | Sigma | R2_conditional | R2_marginal | AIC_wt | AICc_wt | BIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|
| nlme_slope_est_mun | lme | 806.4771 | 919.7565 | NA | NA | 0.9516159 | 0.9528213 | 0.9998575 | 1.0000000 |
| nlme_int_slope_est_mun | lme | 806.5823 | 919.8860 | NA | NA | 0.0483764 | 0.0471709 | 0.0000005 | 0.4008277 |
| nlme_slope_mun | lme | 808.4987 | 922.6896 | NA | 0.9988361 | 0.0000077 | 0.0000077 | 0.0001419 | 0.0000283 |
lmtest::lrtest(nlme_int_slope_est_mun,
nlme_slope_est_mun,
nlme_slope_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 9 | -19875.98 | NA | NA | NA |
| 5 | -19877.00 | -4 | 2.052831 | 0.7260423 |
| 4 | -19891.11 | -1 | 28.210334 | 0.0000001 |
df_ll %>%
filter(biblioteca == "lme4") %>%
group_by(biblioteca) %>%
slice_max(order_by = loglik, n=3) %>%
mutate(rank = row_number()) %>%
ungroup() | modelo | loglik | biblioteca | rank |
|---|---|---|---|
| municipios - inclinação | -19891.11 | lme4 | 1 |
| estados e municípios - inclinação | -19900.90 | lme4 | 2 |
| municipios - intercepto e inclinação | -20132.16 | lme4 | 3 |
performance::compare_performance(lme4_slope_mun,
lme4_slope_est_mun,
lme4_int_slope_mun, rank = TRUE)| Name | Model | R2_conditional | R2_marginal | ICC | RMSE | Sigma | AIC_wt | AICc_wt | BIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|
| lme4_slope_mun | lmerMod | 0.9989630 | 0.8900178 | 0.9905715 | 808.4807 | 922.6150 | 1 | 1 | 1 | 0.5000000 |
| lme4_int_slope_mun | lmerMod | 0.9990904 | 0.8675229 | 0.9931339 | 756.6665 | 874.6162 | 0 | 0 | 0 | 0.4211671 |
| lme4_slope_est_mun | lmerMod | 0.9999999 | 0.0001120 | 0.9999999 | 806.1607 | 918.2947 | 0 | 0 | 0 | 0.2668478 |
lmtest::lrtest(lme4_slope_mun,
lme4_slope_est_mun,
lme4_int_slope_mun)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 4 | -19891.11 | NA | NA | NA |
| 5 | -19900.90 | 1 | 19.57232 | 9.7e-06 |
| 6 | -20132.16 | 1 | 462.53399 | 0.0e+00 |
df_ll %>%
filter(biblioteca == "rstanarm") %>%
group_by(biblioteca) %>%
slice_max(order_by = loglik, n=3) %>%
mutate(rank = row_number()) %>%
ungroup() | modelo | loglik | biblioteca | rank |
|---|---|---|---|
| municipios - inclinação | -19013.34 | rstanarm | 1 |
| municipios - intercepto e inclinação | -19014.47 | rstanarm | 2 |
| estados e municípios - inclinação | -19028.02 | rstanarm | 3 |
performance::compare_performance(rstanarm_slope_mun,
rstanarm_int_slope_mun,
rstanarm_int_slope_est_mun, rank = TRUE)| Name | Model | R2 | R2_marginal | R2_adjusted | R2_adjusted_marginal | ICC | RMSE | Sigma | WAIC_wt | LOOIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|---|
| rstanarm_int_slope_est_mun | stanreg | 0.9990166 | 0.9986830 | 0.9990006 | 0.9990006 | 0.9934342 | 806.2692 | 919.3721 | 0.9999769 | 0.0000050 | 0.5555556 |
| rstanarm_int_slope_mun | stanreg | 0.9990100 | 0.9988376 | 0.9990392 | 0.9990392 | 0.9907607 | 807.8838 | 922.6277 | 0.0000231 | 0.4194687 | 0.4803074 |
| rstanarm_slope_mun | stanreg | 0.9990092 | 0.9988391 | 0.9990364 | 0.9990364 | 0.9905499 | 808.9851 | 922.6888 | 0.0000000 | 0.5805263 | 0.4280889 |
rstanarm::loo_compare(loo_rstanarm_slope_mun,
loo_rstanarm_int_slope_mun,
loo_rstanarm_int_slope_est_mun) elpd_diff se_diff
rstanarm_slope_mun 0.0 0.0
rstanarm_int_slope_mun -1.1 4.0
rstanarm_int_slope_est_mun -16.1 8.5 df_ll %>%
filter(biblioteca == "brms") %>%
group_by(biblioteca) %>%
slice_max(order_by = loglik, n=3) %>%
mutate(rank = row_number()) %>%
ungroup() | modelo | loglik | biblioteca | rank |
|---|---|---|---|
| municipios - inclinação | -19011.98 | brms | 1 |
| estados e municípios - inclinação | -19028.17 | brms | 2 |
| municipios - intercepto e inclinação | -19034.87 | brms | 3 |
performance::compare_performance(brms_slope_mun,
brms_int_slope_mun,
brms_int_slope_est_mun, rank = TRUE)| Name | Model | R2 | R2_marginal | R2_adjusted | R2_adjusted_marginal | ICC | RMSE | Sigma | WAIC_wt | LOOIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|---|---|
| brms_slope_mun | brmsfit | 0.9990124 | 0.9163338 | 0.9990269 | 0.9009658 | 0.9907019 | 808.6345 | 927.7598 | 0.9999938 | 0.9164499 | 0.7730790 |
| brms_int_slope_est_mun | brmsfit | 0.9990156 | 0.0003732 | 0.9989981 | -0.2567222 | 0.9995563 | 807.3770 | 922.9941 | 0.0000062 | 0.0000067 | 0.4444451 |
| brms_int_slope_mun | brmsfit | 0.9990070 | 0.9161204 | 0.9990161 | 0.9013365 | 0.9904779 | 810.1788 | 942.1286 | 0.0000000 | 0.0835433 | 0.3019248 |
brms::loo_compare(
brms::add_criterion(brms_slope_mun,"loo"),
brms::add_criterion(brms_int_slope_mun,"loo"),
brms::add_criterion(brms_int_slope_est_mun,"loo")
) elpd_diff se_diff
brms::add_criterion(brms_slope_mun, "loo") 0.0 0.0
brms::add_criterion(brms_int_slope_mun, "loo") -22.9 9.2
brms::add_criterion(brms_int_slope_est_mun, "loo") -23.1 6.2 Entre as 4 bibliotecas, os modelos que consideraram apenas inclinação aleatória para o nível de municípios foram o que estavam entre os menores valores de log-verossimilhança. Pelo princípio da parcimônia, então, o modelo que considera apenas inclinação aleatória para o nível de municípios é o mais indicado, pois não apresenta grandes diferenças na qualidade de ajuste em comparação com modoelos com mais complexidade. Portanto, este modelo será considerado o mais adequado para a análise dos dados de produção de café no Brasil.
5 O prêmio: O modelo Campeão – compreensão estatística e interpretação dentro do cenário da pergunta de negócio
Vamos interpretar os resultados do modelo escolhido entre as bibliotecas.
5.1 Efeitos aleatórios
Em um modelo hierárquico, os efeitos aleatórios referem-se às variações entre os grupos ou unidades de observação que são modeladas como amostras aleatórias de uma população maior. Essas variações são representadas por parâmetros aleatórios que capturam a incerteza inerente à natureza heterogênea dos dados observados. Em outras palavras, os efeitos aleatórios são uma maneira de modelar a incerteza e a variabilidade inerentes aos dados, reconhecendo que as unidades de observação não são todas iguais e que há uma variação natural entre elas que não pode ser explicada apenas por variáveis fixas ou determinísticas.
Para avaliar o quanto a variabilidade dos dados é explicada pelo efeito aleatório do nível de municípios. Vamos calcular o ICC (Intra-class Correlation) para o modelo escolhido. O ICC, ou Coeficiente de Correlação Intraclasse, é uma medida estatística fundamental no contexto da modelagem hierárquica. Ele quantifica a proporção da variabilidade total da variável de interesse que é atribuída às diferenças entre os grupos ou níveis de aninhamento.
O ICC é calculado como a razão da variância entre os grupos pela variância total:
\[ ICC = \frac{\text{Variância entre grupos}}{\text{Variância total}} \] ou seja:
\[ ICC = \frac{\sigma^2_{\text{município}}}{\sigma^2_{\text{município}} + \sigma^2_{\text{resíduo}}} \]
Para interpretar o ICC, é importante entender que seu valor pode variar entre 0 e 1:
Um ICC próximo de 0 indica que a variação entre os grupos é pequena em relação à variação total, sugerindo que a variável de interesse é pouco influenciada pela estrutura hierárquica dos dados.
Um ICC próximo de 1 indica que a variação entre os grupos é grande em relação à variação total, o que significa que a estrutura hierárquica dos dados tem uma forte influência na variável de interesse.
stderr_nlme(nlme_slope_mun)| RE.Components | Variance.Estimatives | Std.Err. | z | p.value |
|---|---|---|---|---|
| Var(v0j) | 8.645702e+00 | 5.587057e-01 | 15.47452 | 0 |
| Var(e) | 8.513561e+05 | 2.923040e+04 | 29.12571 | 0 |
8.645702 / (8.645702 + 8.513561e+05)[1] 1.015511e-05performance::icc(lme4_slope_mun)| ICC_adjusted | ICC_unadjusted | optional |
|---|---|---|
| 0.9905715 | 0.1089453 | FALSE |
performance::icc(rstanarm_slope_mun)| ICC_adjusted | ICC_unadjusted | optional |
|---|---|---|
| 0.9905499 | 0.1088383 | FALSE |
performance::icc(brms_slope_mun)| ICC_adjusted | ICC_unadjusted | optional |
|---|---|---|
| 0.9907019 | 0.1100612 | FALSE |
5.1.1 Relatório do resultado estatístico do ICC
Enquanto o ICC ajustado se refere apenas aos efeitos aleatórios, o ICC não ajustado também leva em consideração as variâncias dos efeitos fixos, mais precisamente, a variância dos efeitos fixos é adicionada ao denominador da fórmula para calcular o ICC.
Para as bibliotecas lme4, rstanarm e brms, o ICC não ajustado foi de aproximadamente 0.11, o que indica que a variabilidade dos dados é explicada por volta de 11% pelo efeito aleatório do nível de municípios. Para a biblioteca nlme, o ICC ajustado foi de aproximadamente 0.01, o que indica que a variabilidade dos dados é explicada quase que inteiramente pelo efeito fixo do modelo.
5.2 Efeitos fixos
Os efeitos fixos em um modelo hierárquico representam os parâmetros que são considerados constantes e determinísticos em todos os níveis de aninhamento. Eles capturam as relações médias entre as variáveis independentes e dependentes, ou seja, como as mudanças nas variáveis independentes (área plantada) afetam a média da variável dependente (valor de produção), independentemente do grupo ao qual pertencem as observações (Municípios).
Uma forma de avaliar a significância dos efeitos fixos é por meio dos intervalos de confiança (frequentista) ou credibilidade (bayesiano) dos coeficientes estimados. Eles fornecem uma estimativa da incerteza associada aos parâmetros do modelo, indicando a faixa de valores plausíveis para os coeficientes com um certo nível (geralmente 95%).
Avaliar os intervalos de confiança e de credibilidade de 95% para os efeitos fixos é fundamental por diversas razões:
Incerteza nos parâmetros: Os intervalos de confiança e de credibilidade fornecem uma medida da incerteza associada aos parâmetros estimados. Eles indicam a faixa de valores plausíveis para os efeitos fixos, levando em consideração a variabilidade nos dados e o tamanho da amostra.
Inferência estatística: Os intervalos de confiança e de credibilidade ajudam na inferência estatística, permitindo determinar se os efeitos fixos são estatisticamente significativos. Se um intervalo de confiança ou de credibilidade não incluir zero, isso sugere que o efeito é significativamente diferente de zero e tem uma influência mensurável na variável de interesse.
Validação do modelo: Intervalos de confiança e de credibilidade estreitos e bem definidos para os efeitos fixos indicam uma boa precisão e ajuste do modelo aos dados observados. Por outro lado, intervalos amplos podem sugerir que o modelo pode não estar capturando adequadamente a variabilidade nos dados ou que são necessários mais dados para estimar os parâmetros com maior precisão.
nlme::fixed.effects(nlme_slope_mun, which = "fixed") (Intercept) area_plantada
-67.219936 9.502611 nlme::intervals(nlme_slope_mun, which = "fixed", level = 0.95)Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) -123.048312 -67.219936 -11.39156
area_plantada 9.233952 9.502611 9.77127Os resultados dos efeitos fixos mostram que a área plantada tem um efeito positivo sobre a produção de café de R$ 9,50 por hectare. Isso significa que, para cada hectare adicional de área plantada, a produção de café aumenta em R$ 9,50 e que pode variar entre R$ 9,25 e R$ 9,77.
lme4::fixef(lme4_slope_mun) (Intercept) area_plantada
-67.213372 9.502583 stats::confint(lme4_slope_mun, level = 0.95) 2.5 % 97.5 %
.sig01 2.759474 3.131920
.sigma 892.249236 954.348453
(Intercept) -123.040748 -11.451530
area_plantada 9.233775 9.771008Os resultados dos efeitos fixos mostram que a área plantada tem um efeito positivo sobre a produção de café de R$ 9,50 por hectare. Isso significa que, para cada hectare adicional de área plantada, a produção de café aumenta em R$ 9,50 e que pode variar entre R$ 9,23 e R$ 9,77.
rstanarm::fixef(rstanarm_slope_mun) (Intercept) area_plantada
-68.081834 9.509588 rstanarm::posterior_interval(rstanarm_slope_mun, pars = "(Intercept)", prob = 0.95) 2.5% 97.5%
(Intercept) -120.7219 -10.99258rstanarm::posterior_interval(rstanarm_slope_mun, pars = "area_plantada", prob = 0.95) 2.5% 97.5%
area_plantada 9.273183 9.797224Os resultados dos efeitos fixos mostram que a área plantada tem um efeito positivo sobre a produção de café de R$ 9,50 por hectare. Isso significa que, para cada hectare adicional de área plantada, a produção de café aumenta em R$ 9,50 e que pode variar entre R$ 9,27 e R$ 9,79.
brms::fixef(brms_slope_mun) Estimate Est.Error Q2.5 Q97.5
Intercept -68.133665 27.8531919 -121.912316 -13.908125
area_plantada 9.511073 0.1366858 9.256138 9.815771Os resultados dos efeitos fixos mostram que a área plantada tem um efeito positivo sobre a produção de café de R$ 9,50 por hectare. Isso significa que, para cada hectare adicional de área plantada, a produção de café aumenta em R$ 9,51 e que pode variar entre R$ 9,25 e R$ 9,81.
6 A volta da jornada – Discussão
6.1 Diferenças na implementação das bibliotecas
As bibliotecas nlme, lme4, rstanarm e brms são amplamente utilizadas para ajustar modelos de regressão linear misto. Cada uma delas tem suas particularidades e diferenças na implementação dos modelos.
6.1.1 Diferenças na sintaxe dos efeitos aleatórios
A sintaxe para especificar os efeitos aleatórios é diferente entre as bibliotecas. A função lme requer a especificação dos efeitos aleatórios em um agrumento separado na função, enquanto as outras bibliotecas requerem a especificação dos efeitos aleatórios na própria formula.
nlme
função(
formula = variável_dependente ~ variáveis_independentes,
random = ~ intercepto | inclinação,
data = dataset
)lme4, rstanarm e brms
função(
formula = variável_dependente ~ variáveis_independentes + (intercepto | inclinação),
data = dataset
)6.1.2 Diferenças na interpretação dos sumários
Ouput
Linear mixed-effects model fit by REML
Data: df_cafe
AIC BIC logLik
39790.22 39813.14 -19891.11
Random effects:
Formula: ~0 + area_plantada | id_municipio
area_plantada Residual
StdDev: 2.940357 922.6896
Fixed effects: valor_producao ~ area_plantada
Value Std.Error DF t-value p-value
(Intercept) -67.21994 28.455246 1183 -2.36230 0.0183
area_plantada 9.50261 0.136933 1183 69.39595 0.0000
Correlation:
(Intr)
area_plantada -0.355
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-16.71184579 -0.05711394 0.01791640 0.06819056 13.67393999
Number of Observations: 2276
Number of Groups: 1092 Interpretação
Critérios de informação: O AIC (Critério de Informação de Akaike) e o BIC (Critério de Informação Bayesiano) são medidas que avaliam a qualidade do modelo, levando em conta a complexidade do modelo e o ajuste aos dados. Valores menores indicam um melhor modelo. Neste caso, o AIC é 39790.22 e o BIC é 39813.14.
Efeitos aleatórios: O modelo inclui efeitos aleatórios para
area_plantadadentro deid_municipio. Isso significa que o modelo considera variações na relação entre área plantada e valor da produção que não são capturadas pelos efeitos fixos e que podem variar de um município para outro. A variância padrão para area_plantada é de 2.940357 e para o residual é de 922.6896, indicando a magnitude da variação nestes efeitos.Efeitos fixos: A relação principal de interesse é entre
valor_producao(a variável dependente) earea_plantada(a variável independente), ajustando-se por efeitos aleatórios de municípios. O intercepto (valor base devalor_producaoquandoarea_plantadaé 0) é -67.21994, com um erro padrão de 28.455246, significativamente diferente de zero (p-valor = 0.0183). O coeficiente paraarea_plantadaé 9.50261, com um erro padrão de 0.136933, mostrando uma relação positiva forte e significativa comvalor_producao(p-valor praticamente 0).Correlação: Existe uma correlação negativa entre o intercepto e
area_plantada(-0.355), indicando que ajustes na estimativa de um podem influenciar o outro.Resíduos padronizados: Os resíduos, ou diferenças entre os valores observados e os valores preditos pelo modelo, variam de -16.71 a 13.67, com a maioria concentrada perto de zero, o que é esperado em um bom modelo.
Observações e Grupos: O modelo foi ajustado usando 2276 observações, agrupadas em 1092
id_municipio, indicando uma estrutura de dados hierárquica ou agrupada.
Ouput
Linear mixed model fit by REML ['lmerMod']
Formula: valor_producao ~ area_plantada + (0 + area_plantada | id_municipio)
Data: df_cafe
REML criterion at convergence: 39782.2
Scaled residuals:
Min 1Q Median 3Q Max
-16.7132 -0.0571 0.0179 0.0682 13.6750
Random effects:
Groups Name Variance Std.Dev.
id_municipio area_plantada 8.65 2.941
Residual 851218.39 922.615
Number of obs: 2276, groups: id_municipio, 1092
Fixed effects:
Estimate Std. Error t value
(Intercept) -67.213 28.454 -2.362
area_plantada 9.503 0.137 69.380
Correlation of Fixed Effects:
(Intr)
area_plantd -0.355Interpretação
Critério REML na convergência: O valor do critério REML na convergência é 39782.2, que é uma medida de quão bem o modelo se ajusta aos dados, considerando a complexidade do modelo.
Resíduos padronizados: Os resíduos, que são as diferenças entre os valores observados e os valores preditos pelo modelo, variam de -16.7132 a 13.6750, com a maioria dos resíduos concentrados perto do zero (mediana = 0.0179). Isso sugere que o modelo é adequado para explicar a variabilidade nos dados de produção do café.
Efeitos aleatórios: O modelo inclui um efeito aleatório para a
area_plantadadentro de cadaid_municipio, permitindo que a relação entre área plantada e valor da produção varie entre municípios. A variância para esse efeito aleatório é 8.65, com um desvio padrão de 2.941. O grande desvio padrão residual (922.615) indica variação significativa no valor da produção que não é explicada pelo modelo.Observações e Grupos: O modelo foi ajustado usando 2276 observações, agrupadas em 1092 id_municipio. Isso reflete a estrutura de dados hierárquica, onde as medições estão agrupadas por município.
Efeitos fixos: O coeficiente para o intercepto é -67.213, com um erro padrão de 28.454, indicando que, quando a
area_plantadaé 0, o valor esperado da produção é negativo, o que pode ser interpretado como um ajuste do modelo. O coeficiente paraarea_plantadaé 9.503, com um erro padrão de 0.137, mostrando uma relação positiva forte e significativa com o valor da produção (t valor = 69.380). A presença de um coeficiente positivo paraarea_plantadasugere que aumentos na área plantada estão associados a aumentos no valor da produção.Correlação: A correlação entre o intercepto e a
area_plantadaé -0.355, indicando uma relação inversa moderada entre os ajustes desses parâmetros no modelo.
Output
stan_lmer
family: gaussian [identity]
formula: valor_producao ~ area_plantada + (0 + area_plantada | id_municipio)
observations: 2276
------
Median MAD_SD
(Intercept) -68.1 29.2
area_plantada 9.5 0.1
Auxiliary parameter(s):
Median MAD_SD
sigma 922.7 16.1
Error terms:
Groups Name Std.Dev.
id_municipio area_plantada 2.9
Residual 923.9
Num. levels: id_municipio 1092 Interpretação
Família e e Identidade: O modelo utiliza uma distribuição Gaussiana com função de ligação identidade, adequado para dados contínuos e respostas que variam linearmente com os preditores.
Fórmula: A relação entre
valor_producaoearea_plantadaé modelada, permitindo que a inclinação varie porid_municipiosem intercepto para os efeitos aleatórios, o que significa que cada município tem sua própria relação específica entre área plantada e valor da produção, mas o intercepto é comum.Observações: O modelo foi ajustado com 2276 observações.
Efeitos fixos: O valor mediano do intercepto é -68.1, com uma Desvio Médio Absoluto (MAD_SD) de 29.2. Isso indica a estimativa mediana do valor da produção quando a área plantada é 0, com uma variação representada pelo MAD_SD. O valor mediano para a inclinação é 9.5, com MAD_SD de 0.1, indicando uma forte relação positiva entre área plantada e valor da produção, onde aumentos na área plantada estão associados a aumentos no valor da produção.
Parâmetro auxiliar: O parâmetro sigma, que representa o desvio padrão dos erros (ou resíduos) do modelo, tem um valor mediano de 922.7 com MAD_SD de 16.1. Este valor indica a variabilidade dos dados que não é explicada pelo modelo.
Efeitos aleatórios: O desvio padrão para o efeito aleatório
area_plantadaé 2.9, com um valor mediano de 0. O que indica que a inclinação da relação entre área plantada e valor da produção varia entre municípios. O desvio padrão residual é 923.9, semelhante ao parâmetro sigma, indicando a variabilidade dos dados em torno da linha ajustada pelo modelo.Número de níveis: Há 1092 níveis distintos para o
id_municipio, refletindo a variabilidade entre municípios no dataset.
Output
Family: gaussian
Links: mu = identity; sigma = identity
Formula: valor_producao ~ area_plantada + (0 + area_plantada | id_municipio)
Data: df_cafe (Number of observations: 2276)
Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 1000
Group-Level Effects:
~id_municipio (Number of levels: 1092)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(area_plantada) 2.96 0.10 2.78 3.16 1.03 165 381
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -68.13 27.85 -121.91 -13.91 1.00 762 746
area_plantada 9.51 0.14 9.26 9.82 1.01 86 173
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 922.23 15.96 892.31 954.34 1.01 783 665
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Interpretação
Família e e Identidade: O modelo utiliza uma distribuição Gaussiana com função de ligação identidade, adequado para dados contínuos e respostas que variam linearmente com os preditores.
Fórmula: A relação entre valor_producao e area_plantada é modelada, permitindo que a inclinação varie por id_municipio sem intercepto para os efeitos aleatórios, o que significa que cada município tem sua própria relação específica entre área plantada e valor da produção, mas o intercepto é comum.
Dados e Amostragem: O modelo foi aplicado a 2276 observações com 2000 iterações de amostragem distribuídas em duas cadeias, das quais 1000 são pós-aquecimento.
Efeitos aleatórios - Group-Level: O desvio padrão para o efeito aleatório area_plantada é 2.96, com um intervalo de credibilidade de 95% entre 2.78 e 3.16. Isso indica que a inclinação da relação entre área plantada e valor da produção varia entre municípios.
Efeitos fixos - Population-Level: O valor do intercepto é -68.13, com um erro padrão de 27.85, indicando a estimativa do valor da produção quando a área plantada é 0. O coeficiente para area_plantada é 9.51, com um erro padrão de 0.14 e um IC de 95% entre 9.26 e 9.82, mostrando uma relação positiva forte e significativa com o valor da produção. A presença de um coeficiente positivo para area_plantada sugere que aumentos na área plantada estão associados a aumentos no valor da produção.
Parâmetro auxiliar: O parâmetro sigma, que representa o desvio padrão dos erros (ou resíduos) do modelo, tem um valor mediano de 922.23 com um erro padrão de 15.96. Este valor indica a variabilidade dos dados que não é explicada pelo modelo.
Diagnóstico de Convergência: O Rhat é próximo de 1 (idealmente, igual a 1) para todos os parâmetros, indicando boa convergência das cadeias de amostragem. O Bulk_ESS e o Tail_ESS (tamanho efetivo da amostra para a massa da distribuição e para a cauda, respectivamente) são indicadores da qualidade da amostra das simulações. Valores mais altos sugerem maior confiabilidade das estimativas. Neste caso, os valores variam, indicando que a qualidade da amostra pode ser melhor para alguns parâmetros do que para outros, é possível observar que na area_plantada com Bulk_ESS de 86, sugerere uma potencial necessidade de mais iterações nas simulações.
7 Conclusão
7.1 Problema de negócio
Neste relatório, ajustamos um modelo de regressão linear misto para avaliar o efeito da área plantada na produção de café no Brasil. O modelo escolhido foi o que considera apenas inclinação aleatória para o nível de municípios, pois apresentou um bom ajuste e foi o mais parcimonioso entre os modelos ajustados. Os resultados mostram que a área plantada tem um efeito positivo sobre a produção de café de R$ 9,50 por hectare. Isso significa que, para cada hectare adicional de área plantada, a produção de café aumenta em R$ 9,50 e que pode variar entre R$ 9,25 e R$ 9,81. A variabilidade dos dados é explicada por volta de 11% pelo efeito aleatório do nível de municípios.
7.2 Bibliotecas
A comparação entre as bibliotecas nlme, lme4, rstanarm e brms em R para a modelagem de dados hierárquicos apresenta diferenças em termos de usabilidade e funcionalidade. As bibliotecas nlme e lme4 são utilizadas para abordagens frequentistas, enquanto rstanarm e brms são utilizadas para abordagens bayesianas.
A biblioteca nlme é uma das mais antigas para modelagem linear e não-linear mista. É bastante consolidada e possui funções robustas para a estimativa de máxima verossimilhança. Uma particularidade da nlme é a necessidade de especificar os efeitos aleatórios separadamente, o que pode ser menos intuitivo para usuários não avançados.
A lme4, por outro lado, permite a inclusão de efeitos aleatórios diretamente na fórmula do modelo, o que é mais direto e muitas vezes mais conveniente. Ela é mais recente que a nlme e é projetada para ser mais eficiente e flexível, e é possível cálcular a correlação intraclasse (ICC) de forma fácil usando outras bibliotecas como a performance, o que pode ser uma pequena vantagem em comparação a biblioteca nlme.
Quando se trata das abordagens bayesianas, rstanarm e brms são as duas principais bibliotecas usadas como interfaces para Stan, que é uma plataforma para modelagem estatística bayesiana usando a técnica MCMC. rstanarm é projetado para ser mais fácil de usar para aqueles que já estão familiarizados com a sintaxe de fórmulas do lme4, oferecendo uma transição suave para modelagem bayesiana. No entanto, embora seja poderoso, ele pode ser menos flexível do que o brms em termos de especificação do modelo.
brms, embora use uma sintaxe similar à do lme4 e do rstanarm, oferece uma flexibilidade ainda maior para especificar modelos bayesianos complexos. Isso inclui a modelagem de efeitos aleatórios e fixos, distribuições de resposta não normais e relações não lineares. Além disso, brms fornece acesso direto a uma ampla gama de diagnósticos de modelos bayesianos e medidas de ajuste de modelo, como o ICC.
Em resumo, a escolha entre essas bibliotecas depende da complexidade do modelo, da preferência do paradigma estatístico (frequentista versus bayesiano), da familiaridade com a sintaxe da fórmula e da necessidade de flexibilidade no modelamento estatístico. Para modelos mais complexos e quando uma abordagem bayesiana é desejada, brms oferece a maior flexibilidade. Para um usuário menos familiarizado com a programação bayesiana, rstanarm pode ser uma boa introdução. Para aqueles que trabalham dentro do paradigma frequentista, lme4 e nlme oferecem opções robustas, com lme4 proporcionando uma abordagem mais moderna e intuitiva.