Statthon da XLII SEST

Análise de Mortes Violentas em Manaus

Autores

Afiliação

Ingrid Martins

Universidade Federal do Amazonas (UFAM) - Departamento de Estatística

João Lucas

Rayssa Santana

Data de Publicação

22 de outubro de 2025

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.

Modelo 1: 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 2: 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 3: 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)

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.

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)

--- title: "Statthon da XLII SEST" subtitle: "Análise de Mortes Violentas em Manaus" author: - name: "Ingrid Martins" - name: "João Lucas" - name: "Rayssa Santana" institute: "Universidade Federal do Amazonas (UFAM) - Departamento de Estatística" date: today lang: pt-BR format: html: theme: cosmo toc: true toc-location: left toc-title: "Sumário" number-sections: true code-fold: true code-tools: true smooth-scroll: true df-print: paged page-layout: full include-in-header: | <style> h1, h2, h3 { color: #004d40; } .quarto-title-block { background-color: #e0f2f1; padding: 1.5em; border-radius: 10px; box-shadow: 0px 0px 5px rgba(0,0,0,0.1); } footer { text-align: center; color: gray; margin-top: 50px; font-size: 0.9em; } </style> --- # Introdução Este relatório foi desenvolvido no contexto da **Statthon**, por discentes do curso de Estatística da UFAM. # Parte dos Dados e Tratamento # Sistema de Informação sobre Mortalidade (SIM) ```{r} #| echo: false #| warning: false # Pacotes ---------------------------------------------------------------------- require(MASS) # Para glm.nb (Binomial Negativa), se precisarmos require(forecast) # Para auto.arima e previsão de séries temporais require(Matrix) require(hnp) #Carregando e Filtrando os Dados ------- # Carregando pacotes essenciais para a análise # Se não os tiver, instale com: install.packages(c("tidyverse", "lubridate")) require(tidyverse) require(lubridate) require(lme4) require(ggplot2) # Carregar o banco de dados require(gsheet) # Modelos ---------------------------------------------------------------------- mldPol = function(ordem = 2, desconto = .99){ if (ordem == 1){ Ft = matrix(1, c(1,1)) Gt = matrix(1, c(1,1)) Vt = matrix(10,c(1,1)) # variancia do Ruido Wt = matrix(1, c(1,1)) # variancia do sinal desconto = desconto } if (ordem == 2){ Ft = matrix(c(1,0), ncol = 1) Gt = matrix(c(1,0,1,1), 2,2) Vt = matrix(1000, c(1,1)) # variancia do Ruido Wt = diag(1, 2) # variancia do sinal desconto = desconto } quad <- list(Ft = Ft, Gt = Gt, Vt = Vt, Wt = Wt, q = dim(Gt)[1], desconto = desconto) class(quad) <- "eu_amo_mld" return(quad) } mldSaz = function(p, desconto = .9){ Ft = matrix(c(1, rep(0, p-2)), ncol = 1) Gt = rbind(rep(-1,p-1),cbind(diag(1,p-2),0)) Vt = matrix(1,1,1) # variancia do Ruido Wt = diag(.05, p-1) # variancia do sinal desconto = desconto quad <- list(Ft = Ft, Gt = Gt, Vt = Vt, Wt = Wt, q = dim(Gt)[1], desconto = desconto) class(quad) <- "eu_amo_mld" # tem que ser o mesmo nome return(quad) } mldHarm1 <- function(p){ d = 2*pi/p J2 = matrix(c(cos(d), sin(d), -sin(d), cos(d)), 2,2) # primiero harmonico é um padrao anual (h=1) # segundo harmonico é um padrao semestral (h=2) # repete tres vezes antes de fechar o ano, é quadrimestral (h=3) # repete 4 vezes antes de fechar o ano, é trimestral (h=4) #### caso par if ((p/2) == ceiling(p/2)){ quad = list( Ft = matrix(c(1,0,1), ncol = 1), Gt = as.matrix(bdiag(J2,1)), Vt = 1, # variancia do Ruido Wt = diag(1,3) # variancia do sinal ) } class(quad) <- "eu_amo_mld" # tem que ser o mesmo nome return(quad) } "+.eu_amo_mld" <- function(mod1,mod2){ Ft1 <- mod1$Ft Gt1 <- mod1$Gt Vt1 <- mod1$Vt Wt1 <- mod1$Wt Ft2 <- mod2$Ft Gt2 <- mod2$Gt Vt2 <- mod2$Vt Wt2 <- mod2$Wt desconto = c(mod1$desconto,mod2$desconto) q = c(mod1$q,mod2$q) Ft = rbind(Ft1,Ft2) Gt = as.matrix(bdiag(Gt1,Gt2)) Vt = Vt1 + Vt2 Wt = as.matrix(bdiag(Wt1,Wt2)) quad = list(Ft = Ft, Gt = Gt, Vt = Vt, Wt = Wt, q = q, desconto = desconto) class(quad) <- "eu_amo_mld" return(quad) } # W desconto w_desconto <- function(Pt, mod){ # Rt <- Gt*Ct*Gt' + Wt e Gt*Ct*Gt' = Pt delta <- mod$desconto q = mod$q # Wt <- Pt*(1/delta - 1) Wt <- matrix(0, sum(q), sum(q)) k <- 0 for (i in 1:length(q)){ ind = k + 1:q[i] Wt[ind, ind] <- Pt[ind, ind]*(1/delta[i] - 1) k <- cumsum(q)[i] } Wt } # Filtros ---------------------------------------------------------------------- # source("~/Filtros.R", echo=TRUE) filtro <- function(y, mod, Vt = "desconhecido", nt = 0.01, st = 0.01){# m0 = NULL, C0 = NULL, n <- length(y) # Tempo Ft <- mod$Ft Gt <- mod$Gt # Vt <- mod$Vt Wt <- mod$Wt q <- dim(Gt)[1] # Dimensão # Coisas importantes ------ # Priori no tempo t at <- array(NA_real_, c(n,q)) Rt <- array(NA_real_, c(n,q,q)) # Previsão no tempo t ft <- array(NA_real_, n) Qt <- array(NA_real_, n) # Posteriori no tempo t mt <- array(0, c(n+1,q)) Ct <- array(NA_real_, c(n+1,q,q)) # Informação inicial mt[1,][1] <- y[1] Ct[1,,] <- diag(100,q) # Filtro de kalman com v conhecido if (Vt != "desconhecido"){ for (t in 1:n){ # att a priori at[t,] <- Gt %*% mt[t,] # o mt1 é o m0 no filtro de kalman Rt[t,,] <- Gt %*% Ct[t,,] %*% t(Gt) + Wt # Previsão para o tempo t ft[t] <- t(Ft) %*% at[t,] Qt[t] <- t(Ft) %*% Rt[t,,] %*% Ft + Vt # atualiza a posteriori At <- Rt[t,,] %*% Ft %*% solve(Qt[t]) mt[t+1,] <- at[t,] + At %*% (y[t] - ft[t]) Ct[t+1,,] <- Rt[t,,] - At %*% Qt[t] %*% t(At) } } # Filtro de kalman com v desconhecido if (Vt == "desconhecido"){ for (t in 1:n){ # att a priori at[t,] <- Gt %*% mt[t,] # o mt1 é o m0 no filtro de kalman Pt <- Gt %*% Ct[t,,] %*% t(Gt) Rt[t,,] <- Pt + w_desconto(Pt, mod) # Previsão para o tempo t ft[t] <- t(Ft) %*% at[t,] Qt[t] <- t(Ft) %*% Rt[t,,] %*% Ft + 1 # atualiza a posteriori At <- Rt[t,,] %*% Ft %*% solve(Qt[t]) mt[t+1,] <- at[t,] + At %*% (y[t] - ft[t]) Ct[t+1,,] <- Rt[t,,] - At %*% Qt[t] %*% t(At) nt[t+1] <- nt[t] + 1 st[t+1] <- (nt[t]/nt[t+1]) * st[t] + (y[t] - ft[t])^2/(nt[t+1]*Qt[t]) } for(t in 1:n){ Rt[t,,] <- st[t]*Rt[t,,] Qt[t] <- st[t]*Qt[t] Ct[t+1,,] <- st[t+1]*Ct[t+1,,] } # Rt <- sapply(1:n, function(t){st[t]*Rt[t,,]}, simplify = FALSE) # Qt <- sapply(1:n, function(t){st[t]*Qt[t]}, simplify = TRUE) # Ct <- sapply(0:n, function(t){st[t+1]*Ct[t+1,,]}, simplify = FALSE) } resp <- list(y = y, mt = mt, Ct = Ct, at = at, Rt = Rt, ft = ft, Qt = Qt, Ft = Ft, Gt = Gt) return(resp) } # Suavização --------------- suave <- function(filtro){ Gt <- filtro$Gt mt <- filtro$mt Ct <- filtro$Ct at <- filtro$at Rt <- filtro$Rt n <- length(filtro$y) q <- dim(Gt)[1] as <- array(NA_real_, c(n+1,q)) # Media Rs <- array(NA_real_, c(n+1,q,q)) # Variancia as[n+1,] <- mt[n+1,] Rs[n+1,,] <- Ct[n+1,,] for (t in n:1){ Bt <- Ct[t,,] %*% t(Gt) %*% solve (Rt[t,,]) as[t,] <- mt[t,] + Bt %*% (as[t+1,] - at[t,]) Rs[t,,] <- Ct[t,,] - Bt %*% (Rt[t,,] - Rs[t+1,,]) %*% t(Bt) } resp <- list(as = as, Rs = Rs) return(resp) } ``` 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. ```{r} #| echo: false #| warning: false url <- 'https://docs.google.com/spreadsheets/d/1OkY8M8EvO3RD1AVpzz-jimINhC4yvHaDkNdudrhW9Hs/edit?usp=sharing' dados_brutos <- gsheet2tbl(url) # Definindo o intervalo de códigos CID-10 para agressões (homicídios) codigos_agressao <- paste0("X", 85:99) codigos_agressao <- c(codigos_agressao, paste0("Y0", 0:9)) # Filtrando e limpando os dados mortes_manaus <- dados_brutos %>% # 1. Filtrar óbitos ocorridos em Manaus (código 130260) filter(CODMUNOCOR == 130260) %>% # 2. Filtrar apenas as causas básicas que são agressões filter(substr(CAUSABAS, 1, 3) %in% codigos_agressao) # Visualizando as primeiras linhas do dataset filtrado head(mortes_manaus) # Tratamento das Variáveis ------------ # Função para decodificar a variável IDADE decodificar_idade <- function(idade_cod) { # A idade em anos começa com o dígito 4. Ex: 425 -> 25 anos # Outros prefixos (1=min, 2=h, 3=d) serão considerados como idade < 1 ano (0). idade_str <- as.character(idade_cod) ifelse( substr(idade_str, 1, 1) == "4", as.numeric(substr(idade_str, 2, 3)), 0 ) } # Aplicando as transformações mortes_manaus_tratado <- mortes_manaus %>% mutate( # Converter data do óbito para o formato de data data_obito = dmy(DTOBITO), # Decodificar idade para anos idade_anos = decodificar_idade(IDADE), # Decodificar sexo para fator com rótulos sexo = factor(SEXO, levels = c(1, 2), labels = c("Masculino", "Feminino")), # Decodificar raça/cor para fator com rótulos raca_cor = factor(RACACOR, levels = c(1, 2, 3, 4, 5), labels = c("Branca", "Preta", "Amarela", "Parda", "Indígena")), # Extrair o ano do óbito para análises temporais ano_obito = year(data_obito) ) %>% # Selecionar as colunas de interesse select(data_obito, ano_obito, idade_anos, sexo, raca_cor, causa_bas = CAUSABAS) # Visualizar o resultado final # glimpse(mortes_manaus_tratado) # Exemplo de Visualização Inicial --------- mortes_manaus_tratado %>% filter(!is.na(ano_obito)) %>% # Remover datas que não puderam ser convertidas count(ano_obito) %>% ggplot(aes(x = ano_obito, y = n)) + geom_line(color = "#0072B2", size = 1) + geom_point(color = "#0072B2", size = 2) + labs( title = "Evolução das Mortes por Agressão em Manaus", subtitle = "Dados do Sistema de Informação sobre Mortalidade (SIM)", x = "Ano do Óbito", y = "Número de Óbitos Registrados" ) + theme_minimal() + scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) ``` # 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. - Modelo 1: 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 2: 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 3: 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. ```{r} 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) BIC(modelo_poisson, modelo_nb, modelo_nb2) 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') hnp(modelo_nb, xlab = 'Percentil da N(0,1)', ylab = 'Resíduos', main = 'Binomial Negativa') ``` 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. # 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. ```{r} # --- 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. ```{r} # --- 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) ```