Statthon da XLII SEST
Análise de Mortes Violentas em Manaus
1 Introdução
Este relatório foi desenvolvido no contexto da Statthon, por discentes do curso de Estatística da UFAM.
2 Parte dos Dados e Tratamento
3 Sistema de Informação sobre Mortalidade (SIM)
Os dados utilizados nesta competição são uma amostra real do Sistema de Informações sobre Mortalidade (SIM), gerenciado pelo Ministério da Saúde do Brasil através do SUS.
O SIM foi criado em 1975 para coletar e consolidar dados de mortalidade de todo o país a partir das Declarações de Óbito (DO). Ele é uma das principais fontes de informação sobre a saúde da população brasileira, permitindo:
- Monitorar o perfil de mortalidade do país, estados e municípios.
- Analisar tendências de doenças, acidentes e violência.
- Subsidiar o planejamento de políticas públicas de saúde e segurança.
- Apoiar a pesquisa científica na área da saúde.
Os dados do SIM são essenciais para entender as causas de morte e os fatores demográficos associados a elas, como idade, sexo e raça/cor.
4 Análise utlizando o Modelo Poisson e Binomial Negativa
A distribuição binomial negativa é uma generalização da distribuição de Poisson, usada para modelar dados de contagem quando há superdispersão, ou seja, quando a variância é maior que a média. Enquanto a Poisson assume que a média e a variância são iguais, a binomial negativa introduz um parâmetro de dispersão que permite capturar a variabilidade extra nos dados.
Em termos práticos, ela é útil para modelar o número de ocorrências de um evento (como acidentes, mortes, falhas ou casos de doença) em situações em que há heterogeneidade não explicada ou observações mais variáveis do que o esperado sob a Poisson. Como os nossos dados são contagens de mortes violentas em Manaus, esse é o modelo adequado para capturar a variabilidade dos dados.
Para analisar a contagem de óbitos (n_mortes), foi utilizada uma abordagem progressiva de Modelos Lineares Generalizados (GLM), comparando diferentes distribuições e estruturas de regressores em seguida foi feito os resultados foram comparados utilizando Modelos Lineares Generalizados com Efeitos Mistos (GLMM).
4.1 Modelos Lineares Generalizados (GLM)
- Modelo Poisson (Modelo Aditivo)
O modelo de Poisson é frequentemente o ponto de partida para a análise de dados de contagem. A variável resposta \(Y_i\) (representando n_mortes) é assumida seguir uma distribuição de Poisson com média \(\mu_i\).
\[ Y_i \sim \text{Poisson}(\mu_i) \]
A média é então modelada em relação aos preditores através de uma função de ligação logarítmica, resultando em um modelo aditivo (sem interação):
\[ \log(\mu_i) = \beta_0 + \beta_1 \, \text{Sexo}_i + \beta_2 \, \text{RaçaCor}_i \]
Este modelo estima o efeito principal de \(\texttt{Sexo}\) e \(\texttt{RaçaCor}\) na contagem de óbitos. Sua principal limitação é a suposição de equidispersão, onde a variância é igual à média. Se a variabilidade dos dados for maior que a média, este modelo pode subestimar os erros padrão, sendo assim propomos um modelo Binomial Negativo com o intuito de ser mais apropriado.
Modelo de Regressão Binomial Negativa com Efeito Aleatório
O modelo ajustado é:
\[ Y_{ij} \sim \text{NegBin}(\mu_{ij}, \theta) \]
\[ \log(\mu_{ij}) = \beta_0 + \beta_1 \, \text{Sexo}_i + \beta_2 \, \text{RaçaCor}_i + u_j \]
\[ u_j \sim \mathcal{N}(0, \sigma^2_{\text{ano}}) \]
onde:
\(Y_{ij}\) = número de mortes da observação \(i\) no ano \(j\)
\(\mu_{ij}\) = média esperada de mortes
\(\beta_0\) = intercepto
\(\beta_1, \beta_2\) = efeitos fixos de sexo e raça/cor
\(u_j\) = efeito aleatório do ano de óbito
\(\theta\) = parâmetro de dispersão da binomial negativa
Modelo GLM Binomial Negativa (Aditivo)
O modelo ajustado é:
\[ Y_i \sim \text{NegBin}(\mu_i, \theta) \]
\[ \log(\mu_i) = \beta_0 + \beta_1 \, \text{Sexo}_i + \beta_2 \, \text{RaçaCor}_i \]
onde:
\(Y_i\) = número de mortes da observação \(i\)
\(\mu_i\) = média esperada de mortes
\(\beta_0\) = intercepto
\(\beta_1, \beta_2\) = efeitos fixos de sexo e raça/cor
\(\theta\) = parâmetro de dispersão da binomial negativa
Modelo GLM Binomial Negativa (Com Interação)
O modelo ajustado é:
\[ Y_i \sim \text{NegBin}(\mu_i, \theta) \]
\[ \log(\mu_i) = \beta_0 + \beta_1 \, \text{Sexo}_i + \beta_2 \, \text{RaçaCor}_i + \beta_3 \, (\text{Sexo}_i \times \text{RaçaCor}_i) \]
- Aqui incluímos a interação \(\beta_3\), permitindo que o efeito de Sexo dependa da RaçaCor.
A seleção do modelo final baseia-se na comparação de métricas de ajuste, como o AIC (Akaike Information Criterion) e o Teste da Razão de Verossimilhança (LRT), para identificar o modelo mais parcimonioso que melhor descreve a variabilidade dos dados.
Código
dados_agrupados <- mortes_manaus_tratado %>%
group_by(sexo, idade_anos, raca_cor, ano_obito) %>%
summarise(n_mortes = n(), .groups = "drop")
# Ajuste GLM Poisson
modelo_poisson <- glm(
n_mortes ~ sexo + raca_cor,
data = dados_agrupados,
family = poisson
)
# Ajuste GLM Binomial Negativa
modelo_nb <- glm.nb(
n_mortes ~ sexo + raca_cor,
data = dados_agrupados
)
# Ajuste GLM Binomial Negativa
modelo_nb2 <- glm.nb(
n_mortes ~ sexo * raca_cor,
data = dados_agrupados
)
# Comparação lado a lado
AIC(modelo_poisson, modelo_nb, modelo_nb2)Código
BIC(modelo_poisson, modelo_nb, modelo_nb2)Código
newdata <- expand.grid(
sexo = levels(dados_agrupados$sexo),
raca_cor = levels(dados_agrupados$raca_cor)
)
newdata <- newdata %>%
mutate(
taxa_poisson = predict(modelo_poisson, newdata, type = "response"),
taxa_nb = predict(modelo_nb, newdata, type = "response")
)
newdata$taxa_mortes <- predict(modelo_nb, newdata, type = "response")
tabela_taxas <- newdata %>%
arrange(sexo, raca_cor) %>%
mutate(taxa_mortes = round(taxa_mortes, 3)) # arredonda para 3 casas
par(mfrow = c(1,2))
hnp(modelo_poisson, xlab = 'Percentil da N(0,1)',
ylab = 'Resíduos',
main = 'Poisson')Poisson model
Código
hnp(modelo_nb, xlab = 'Percentil da N(0,1)',
ylab = 'Resíduos',
main = 'Binomial Negativa')Negative binomial model (using MASS package)
Código
hnp(modelo_nb2, xlab = 'Percentil da N(0,1)',
ylab = 'Resíduos',
main = 'Binomial Negativa')Negative binomial model (using MASS package)
A análise do modelo permitiu quantificar o impacto de cada variável na mortalidade por agressão. A taxa de referência (Intercepto) para o grupo base composto por indivíduos do sexo masculino e da raça/cor branca foi estimada em 2,31 (IC95% 2,10; 2,53).
A variável sexo mostrou-se um fator protetor significativo. Indivíduos do sexo feminino apresentaram uma Razão de Taxas (RR) de 0,154 (IC95% 0,139; 0,170). Isso indica que o risco de morte por agressão para mulheres é aproximadamente 85% menor em comparação direta com os homens.
A raça/cor também revelou diferenças substanciais. Em comparação com o grupo de referência (brancos), a maioria dos outros grupos apresentou um risco menor. Indivíduos da raça/cor preta (RR = 0,545; IC95% 0,434; 0,684) e amarela (RR = 0,553; IC95% 0,401; 0,763) tiveram uma redução de risco semelhante, cerca de 45%. O grupo indígena apresentou a maior redução, com um risco 51% menor (RR = 0,494; IC95% 0,320; 0,762).
Em forte contraste, a raça/cor parda emergiu como o principal fator de risco. Com uma Razão de Taxas de 5,589 (IC95% 5,03; 6,21), os resultados indicam que o risco de morte por agressão para um indivíduo pardo é aproximadamente 5,6 vezes maior do que para um indivíduo branco.
4.2 Modelos Lineares Generalizados com Efeitos Mistos (GLMM)
Código
# ---- Ajuste dos modelos ----
modelo_glmer_pois <- glmer(
n_mortes ~ sexo + raca_cor + (1 | ano_obito),
data = dados_agrupados,
family = poisson(link = "log")
)
modelo_glmer_nb <- glmer.nb(
n_mortes ~ sexo + raca_cor + (1 | ano_obito),
data = dados_agrupados
)
modelo_glmer_nb4 <- glmer.nb(
n_mortes ~ sexo * raca_cor * (1 | ano_obito),
data = dados_agrupados
)
# ---- Usar o mesmo data frame que o modelo usou ----
dados_modelo <- model.frame(modelo_glmer_pois)
#---- Adicionar resíduos e ajustes ----
dados_residuos <- dados_modelo %>%
mutate(
ajuste_pois = fitted(modelo_glmer_pois),
ajuste_nb = fitted(modelo_glmer_nb),
resid_pois = residuals(modelo_glmer_pois, type = "pearson"),
resid_nb = residuals(modelo_glmer_nb, type = "pearson"),
resid_nb4 = residuals(modelo_glmer_nb4, type = "pearson")
)
# ---- QQ-plot dos resíduos ----
res_pois <- qqnorm(dados_residuos$resid_pois, plot.it = FALSE)
res_nb <- qqnorm(dados_residuos$resid_nb, plot.it = FALSE)
res_nb4 <- qqnorm(dados_residuos$resid_nb4, plot.it = FALSE)
par(mfrow = c(1, 2))
plot(res_pois, main = "GLMM Poisson", xlab = "Quantis N(0,1)", ylab = "Resíduos de Pearson")
abline(0, 1, col = "red", lty = 2)
plot(res_nb, main = "GLMM Binomial Negativa", xlab = "Quantis N(0,1)", ylab = "Resíduos de Pearson")
abline(0, 1, col = "red", lty = 2)Código
# Comparação lado a lado
AIC(modelo_glmer_pois, modelo_glmer_nb, modelo_glmer_nb4)Código
BIC(modelo_glmer_pois, modelo_glmer_nb, modelo_glmer_nb4)Código
# ---- Criar novo conjunto com todas as combinações das variáveis ----
newdata <- expand.grid(
sexo = levels(dados_agrupados$sexo),
raca_cor = levels(dados_agrupados$raca_cor)
)
# ---- Prever taxas (médias esperadas) sem incluir efeitos aleatórios ----
newdata$taxa_pois <- predict(modelo_glmer_pois, newdata, type = "response", re.form = NA)
newdata$taxa_nb <- predict(modelo_glmer_nb, newdata, type = "response", re.form = NA)
# ---- Tabela formatada ----
tabela_taxas <- newdata %>%
arrange(sexo, raca_cor) %>%
rename(
"Sexo" = sexo,
"Raça/Cor" = raca_cor,
"Taxa Poisson" = taxa_pois,
"Taxa Binomial Negativa" = taxa_nb
)
# Mostrar tabela
tabela_taxasPara levar em conta a estrutura de agrupamento dos dados onde observações dentro do mesmo ano_obito podem estar correlacionadasos modelos GLM são estendidos para Modelos Lineares Generalizados Mistos (GLMM). Isso é feito pela inclusão de um efeito aleatório (um intercepto) para cada ano.
- Modelo GLMM Poisson (Efeito Aleatório Aditivo)
Este é o primeiro modelo misto, que estende o modelo Poisson para incluir um efeito aleatório \(u_j\) para cada ano \(j\). A variável \(Y_{ij}\) representa a contagem de óbitos para a \(i\)-ésima combinação de preditores no ano \(j\).
\[ Y_{ij} \sim \text{Poisson}(\mu_{ij}) \]
A média \(\mu_{ij}\) é modelada na escala logarítmica com efeitos fixos aditivos e o efeito aleatório do ano:
\[ \log(\mu_{ij}) = \beta_0 + \beta_1 \, \text{Sexo}_i + \beta_2 \, \text{RaçaCor}_i + u_j \]
\[ u_j \sim \mathcal{N}(0, \sigma^2_{\text{ano}}) \]
O termo \(u_j\) quantifica a variabilidade entre os anos que não é explicada pelos preditores fixos. Este modelo ainda assume equidispersão (média = variância) condicional ao efeito aleatório.
Modelo GLMM Binomial Negativa (Aditivo)
Este modelo combina a flexibilidade da distribuição Binomial Negativa (para tratar a superdispersão) com a estrutura de efeitos mistos (para tratar a correlação entre os anos).
\[ Y_{ij} \sim \text{NegBin}(\mu_{ij}, \theta) \]
O preditor linear é idêntico ao do modelo GLMM Poisson, controlando para os efeitos fixos aditivos e o efeito aleatório do ano:
\[ \log(\mu_{ij}) = \beta_0 + \beta_1 \, \text{Sexo}_i + \beta_2 \, \text{RaçaCor}_i + u_j \]
\[ u_j \sim \mathcal{N}(0, \sigma^2_{\text{ano}}) \]
Este modelo é comporta bem os dados, pois lida simultaneamente com a superdispersão (através do parâmetro \(\theta\)) e com a variabilidade entre os anos (através do efeito aleatório \(u_j\)).
Modelo GLMM Binomial Negativa (Com Interação)
Finalmente, o modelo mais complexo é ajustado para testar se a relação entre \(\text{RaçaCor}\) e a mortalidade difere entre os níveis de \(\text{Sexo}\), mantendo a estrutura de efeitos mistos.
\[ Y_{ij} \sim \text{NegBin}(\mu_{ij}, \theta) \]
O preditor linear agora inclui o termo de interação \(\beta_3\) entre os efeitos fixos:
\[ \log(\mu_{ij}) = \beta_0 + \beta_1 \, \text{Sexo}_i + \beta_2 \, \text{RaçaCor}_i + \beta_3 \, (\text{Sexo}_i \times \text{RaçaCor}_i) + u_j \]
\[ u_j \sim \mathcal{N}(0, \sigma^2_{\text{ano}}) \]
- Este modelo permite uma análise mais detalhada de como os fatores de risco se combinam, enquanto ainda controla pela heterogeneidade não observada entre os anos.
A análise dos Modelos Lineares Generalizados Mistos (GLMM) foi focada na seleção do modelo mais adequado e na interpretação dos fatores de risco. A comparação dos modelos, utilizando os critérios de informação AIC e BIC, foi conclusiva. Primeiramente, o critério AIC do modelo GLMM Poisson (AIC = 22384.75) foi massivamente superior ao do GLMM Binomial Negativo (AIC = 12489.26). Esta diferença é uma evidência estatística esmagadora de que os dados possuem severa superdispersão, tornando o modelo Poisson totalmente inadequado e o uso da distribuição Binomial Negativa obrigatório.
Posteriormente, o modelo Binomial Negativo aditivo (AIC = 12489.26) foi comparado a um modelo com interação entre sexo e raça/cor, que apresentou um AIC de 12392.99. A redução de aproximadamente 96 pontos no AIC, confirmada também pelo critério BIC (12462.54 vs 12535.63), indica que o modelo com interação é significativamente melhor. Isso demonstra que o risco associado à raça/cor não é o mesmo para homens e mulheres, e que os efeitos dependem um do outro.
Embora o modelo de interação seja o vencedor, a análise das taxas esperadas do modelo aditivo (Taxa Binomial Negativa) oferece uma visão clara dos efeitos principais. O sexo feminino apresenta um risco de mortalidade drasticamente menor em todos os grupos de raça/cor. Por exemplo, a taxa esperada para um indivíduo do grupo “Masculino/Branca” foi de 2,23, enquanto para o grupo “Feminino/Branca” foi de apenas 0,34.
O grupo “Parda” emergiu como o de maior risco. A taxa esperada para “Masculino/Parda” (12,68) é aproximadamente 5,7 vezes maior que a do grupo “Masculino/Branca” (2,23). Tendência similar foi observada no grupo feminino, onde a taxa para “Feminino/Parda” (1,94) foi 5,7 vezes maior que para “Feminino/Branca” (0,34). Os demais grupos (“Preta”, “Amarela” e “Indígena”) apresentaram taxas ligeiramente inferiores às do grupo “Branca” no modelo aditivo.
5 Modelo Linear Dinâmico
Um Modelo Linear Dinâmico é uma forma flexível de modelar séries temporais que mudam ao longo do tempo ele permite que os parâmetros (como a tendência e a sazonalidade) evoluam dinamicamente, em vez de serem fixos como nos modelos ARIMA tradicionais.
Em vez de assumir que o processo é estável, o MLD reconhece que o comportamento da série pode mudar por exemplo, uma tendência que cresce e depois estabiliza, ou uma sazonalidade que enfraquece com o tempo.
No nosso caso (óbitos em Manaus), o MLD tem vários propósitos práticos:
Capturar a estrutura da série: separar tendência, sazonalidade e ruído, para entender melhor o comportamento real da mortalidade. Isso responde perguntas como: “os óbitos aumentam de forma constante?” ou “há meses com picos sazonais recorrentes?”
Analisar estabilidade temporal: o fator de desconto do MLD permite validar quão rápido o sistema muda. Se for baixo, quer dizer que a tendência muda rapidamente; se for alto, a tendência é estável.
Filtrar e suavizar ruído: em séries com muita variação (como a mensal), o MLD fornece uma estimativa limpa do sinal subjacente sendo essencial para uma melhor interpretação.
Fazer previsão: como o modelo é baseado no filtro de Kalman, ele permite projeções futuras considerando incerteza e dinâmica do processo. Isso é útil para por exemplo o planejamento de saúde pública, monitoramento de epidemias ou avaliação de políticas sanitárias.
Código
# --- Dados MENSAIS ---
# (O 'y' são as contagens por mês)
y_mensal <- mortes_manaus_tratado %>%
filter(!is.na(data_obito)) %>%
# Criar uma coluna de 'ano_mes' para agrupar
mutate(ano_mes = floor_date(data_obito, "month")) %>%
count(ano_mes) %>%
pull(n)
y_mensal <- ts(y_mensal)
# --- Modelo Mensal (Tendência + Sazonalidade p=12) ---
mod_tend_mensal <- mldPol(ordem = 2, desconto = 0.98)
mod_saz_mensal <- mldSaz(p = 12, desconto = 0.99)
mod_mensal <- mod_tend_mensal + mod_saz_mensal
# --- Rodar o Filtro ---
filtro_mensal <- filtro(y_mensal, mod_mensal, Vt = "desconhecido")
# --- Plotar Resultado ---
plot(y_mensal, col = 'black',#type = 'l', #'grey80',
main = "Modelo Mensal (Tendência + Sazonalidade p=12)",
xlab = "Tempo (Meses)", ylab = "Nº de Óbitos")
# Nível (estado 1) + Sazonal (estado 2)
# mt[-1, 1] é a tendência (do mldPol)
# mt[-1, 2] é o primeiro componente sazonal (do mldSaz)
nivel_mensal <- filtro_mensal$mt[-1, 1]
sazonal_mensal <- filtro_mensal$mt[-1, 2]
# Série ajustada (Nível + Sazonalidade)
serie_ajustada_mensal <- nivel_mensal + sazonal_mensal
lines(serie_ajustada_mensal, col = "blue", lwd = 2)
legend("topleft", legend = c("Observado (Mensal)", "Ajuste (Nível + Sazonal)"),
col = c("black", "blue"), lty = 1, lwd = 2, cex = 0.8)Esse gráfico é uma excelente visualização do que um Modelo Linear Dinâmico (MLD) faz. Ele mostra o processo de “filtragem”, onde o modelo tenta separar o sinal verdadeiro do ruído aleatório.
Como podemos ver, ela tem um ruido muito alto. Há picos e vales extremos (como o pico agudo perto do tempo 130) que podem ser eventos atípicos.
A linha azul é a estimativa do MLD. É a “versão limpa” da série temporal, calculada pelo filtro de Kalman.
O modelo foi instruído a encontrar duas coisas nos dados:
mldPol(ordem = 2): Uma Tendência: um comportamento de crescimento ou decrescimento observados por um período não muito curto de tempo. A linha azul captura isso muito bem: ela sobe do início até ~o mês 75, depois se estabiliza em um patamar mais alto (entre 80 e 100 óbitos/mês).
mldSaz(p = 12): Uma Sazonalidade: um comportamento que se repetem em intervalos regulares de tempo. Esse intervalo de repetição é denominado período. As pequenas “ondulações” regulares na linha azul são a estimativa do modelo para esse padrão.
Além disso esse modelo possui O Fator de Desconto que controla a memória e a adaptabilidade do modelo.
Código
# --- Dados ANUAIS ---
# (O 'y' são as contagens por ano)
y_anual <- mortes_manaus_tratado %>%
filter(!is.na(ano_obito)) %>%
count(ano_obito) %>%
pull(n)
# --- Modelo de Tendência (Polinômio Ordem 2) ---
# Ordem 2 (nível e inclinação) captura melhor a subida e descida
mod_anual <- mldPol(ordem = 2, desconto = 0.8)
# --- Rodar o Filtro ---
filtro_anual <- filtro(y_anual, mod_anual, Vt = "desconhecido")
# --- Plotar Resultado ---
plot(y_anual, type = 'l', col = 'grey80', lwd = 2,
main = "Modelo Anual: Tendência (Nível)",
xlab = "Tempo (Anos)", ylab = "Nº de Óbitos")
# O nível é a primeira coluna do estado 'mt'
nivel_anual <- filtro_anual$mt[-1, 1]
lines(nivel_anual, col = "red", lwd = 2)
legend("topleft", legend = c("Observado (Anual)", "Tendência (Nível)"),
col = c("grey80", "red"), lty = 1, lwd = 2)