Curso Visualização de dados de interesse para a vigilância em saúde

Módulo 7 - Modelando séries temporais com GAM

Agora vamos ver outra classe de modelos que também é capaz de capturar algumas estruturas de uma série temporal: os Modelos Aditivos Generalizados (GAMs).

Em um modelo linear generalizado (GLM), como uma regressão linear ou uma regressão logística, estamos sujeitos a explicar relações lineares entre os dados. O modelo aditivo generalizado (GAM) é uma extensão dessa classe de modelos que permite a incorporação de termos não lineares, chamados de funções suaves.

Por exemplo, em nossa série de SRAG no Paraná, vamos observar o efeito da variável mês em relação ao número de casos na Figura 42.

Figura 42: Distribuição do número de casos de SRAG no PR por mês de primeiros sintomas, 2013-2017.

Essa é, como vimos nos módulos anteriores, uma boa técnica para observar a sazonalidade dos dados. Vemos um aumento de casos nos meses intermediários do ano, principalmente nos meses de outono/inverno. Contudo, percebemos que a relação entre casos e meses não é linear, pois apresenta um aumento e depois uma diminuição, conforme avançamos os meses.

Mesmo assim, vamos tentar modelar essa relação por meio de um modelo linear generalizado (GLM), o modelo de Poisson:

mod <- glm(n ~ mes, data=srag_pr_mes, family = "poisson")



summary(mod)

#> 

#> Call:

#> glm(formula = n ~ mes, family = "poisson", data = srag_pr_mes)

#> 

#> Coefficients:

#>              Estimate Std. Error z value Pr(>|z|)    

#> (Intercept)  5.951726   0.014318 415.686  < 2e-16 ***

#> mes         -0.010897   0.001976  -5.514 3.51e-08 ***

#> ---

#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#> 

#> (Dispersion parameter for poisson family taken to be 1)

#> 

#>     Null deviance: 14492  on 59  degrees of freedom

#> Residual deviance: 14461  on 58  degrees of freedom

#> AIC: 14907

#> 

#> Number of Fisher Scoring iterations: 5

A função glm() ajusta um modelo de regressão linear generalizado, neste caso em particular será levado em consideração os dados com distribuição de Poisson.
O modelo tem a variável resposta n (o número de casos de SRAG) sendo explicada pela variável preditora mes (mês de primeiros sintomas) no conjunto de dados srag_pr_mes. O argumento family = "poisson" indica que o modelo deve usar a distribuição de Poisson, adequada para contagens.
A função summary(mod) exibe um resumo do modelo ajustado, incluindo coeficientes estimados, significância estatística, desvio residual e outras métricas relevantes para avaliar o ajuste do modelo.

Vemos que o mês de início de sintomas teve uma associação inversa estatisticamente significativa com o número de casos. Porém, vamos verificar os valores ajustados pelo modelo na Figura 43:

srag_pr_mes <- srag_pr_mes %>% 

  mutate(dt_mes = ymd(paste0(ano, "-", mes, "-01")))



srag_pr_mes$est_poisson <- fitted.values(mod)



ggplot(srag_pr_mes, aes(x=dt_mes)) + 

  geom_line(aes(y=n, color="Observado")) + 

  geom_point(aes(y=est_poisson, color = "Ajustado")) +

  scale_x_date(date_breaks = "1 year",

               date_labels = "%Y") + 

  labs(x="Mês/ano", y = "Casos de SRAG", color = "")

Figura 43: Casos de SRAG no Paraná observados versus ajustados por um modelo linear generalizado.

O script cria um gráfico utilizando a biblioteca ggplot2 para comparar os valores observados de casos de Síndrome Respiratória Aguda Grave (SRAG) ao longo do tempo com os valores ajustados por um modelo linear generalizado:

srag_pr_mes$est_poisson <- fitted.values(mod): Adiciona uma nova coluna ao dataset srag_pr_mes chamada est_poisson, que armazena os valores ajustados pelo modelo linear generalizado armazenado em mod.
ggplot(srag_pr_mes, aes(x=dt_mes)): Inicia a criação de um gráfico com ggplot2, definindo o eixo x como dt_mes, que representa as datas (mês/ano).
geom_line(aes(y=n, color="Observado")): Adiciona uma linha ao gráfico para os valores observados (variável n), atribuindo a cor “Observado”.
geom_point(aes(y=est_poisson, color="Ajustado")): Adiciona pontos representando os valores ajustados (variável est_poisson), atribuindo a cor “Ajustado”.
scale_x_date(): Define que o eixo x será formatado com intervalos anuais (date_breaks = "1 year") e rótulos no formato de ano (date_labels = "%Y").
labs(x="Mês/ano", y="Casos de SRAG", color=""): Adiciona rótulos aos eixos e remove o título da legenda de cores.

Veja que o modelo linear generalizado não foi capaz de capturar o efeito dos meses sobre os dados, justamente porque esse efeito é não linear, como vimos no gráfico de boxplot anterior.

É aí que entra o modelo aditivo generalizado (GAM): ele permite que incluamos termos não lineares por meio de funções suaves cuja complexidade conseguimos controlar. Isso é útil para modelar relações como essa (entre mês e número de casos - sazonalidade), para modelar a tendência nos dados (seja ela linear ou não), e modelar o efeito de outras variáveis sobre nosso desfecho de interesse.

Vamos realizar o mesmo exemplo anterior com os dados de SRAG no Paraná, mas agora com um modelo GAM. Para isso, vamos utilizar a função gam() da biblioteca mgcv. Instale-a primeiro, se for preciso, utilizando o código install.packages("mgcv"). Agora, em vez de especificar o modelo como n ~ mes, utilizamos n ~ s(mes) para indicar que uma função suave não linear será aplicada à variável.

library(mgcv) 



mod_gam <-  gam(n ~ s(mes), data=srag_pr_mes, family = "poisson")

A biblioteca mgcv é carregada para possibilitar o uso de Modelos Aditivos Generalizados (GAM), que permitem modelar relações não lineares entre variáveis.
A função gam() ajusta um modelo GAM, onde:
- A variável resposta n (número de casos de SRAG) é modelada como uma função suave do termo mes (mês), representado por s(mes). Isso permite capturar relações não lineares entre o mês e os casos.
- O parâmetro family = "poisson" indica que os dados seguem uma distribuição de Poisson, comumente utilizada para modelar contagens.

Pronto! Agora vamos plotar o resultado da relação entre mês e casos de SRAG na Figura 44.

plot(mod_gam, select = 1)

Figura 44: Relação entre mês de primeiros sintomas e casos de SRAG no Paraná estimada com um modelo GAM.

A função plot() gera um gráfico para visualizar a relação não linear estimada entre o mês (mes) e os casos de SRAG (n). O argumento select = 1 garante que a suavização do primeiro termo seja exibida no gráfico.

Veja só! Agora o modelo, ainda que simples, conseguiu capturar melhor a relação entre o mês de primeiros sintomas e o número de casos de SRAG. Percebe-se que o aumento de casos, com ápice no mês 6, assim como o menor número no começo e final do ano, foram capturados pelo modelo.

Vamos na Figura 45 comparar os valores ajustados pelos modelos GLM e GAM com a série observada:

srag_pr_mes$est_gam <- fitted.values(mod_gam)



ggplot(srag_pr_mes, aes(x=dt_mes)) + 

  geom_line(aes(y=n, color="Observado")) + 

  geom_point(aes(y=est_poisson, color = "Ajustado (GLM)")) +

  geom_point(aes(y=est_gam, color = "Ajustado (GAM)")) +

  scale_x_date(date_breaks = "1 year",

               date_labels = "%Y") + 

  labs(x="Mês/ano", y = "Casos de SRAG", color = "")

Figura 45. Casos de SRAG no Paraná observados comparado com valores ajustados por modelo GAM e por modelo GLM.

Este trecho de código está criando um gráfico de linhas e pontos utilizando o pacote ggplot2 no R para comparar dados observados e ajustados em um conjunto de dados sobre casos de Síndrome Respiratória Aguda Grave (SRAG), contidos no objeto srag_pr_mes. A seguir, uma explicação dos elementos:

srag_pr_mes$est_gam <- fitted.values(mod_gam): Adiciona uma coluna chamada est_gam ao conjunto de dados srag_pr_mes, contendo os valores ajustados do modelo GAM armazenado em mod_gam.
ggplot(srag_pr_mes, aes(x=dt_mes)): Inicializa um gráfico com o eixo x representando os meses (dt_mes).
geom_line(aes(y=n, color="Observado")): Adiciona uma linha ao gráfico representando os valores observados de SRAG (n), com a cor definida como “Observado”.
geom_point(aes(y=est_poisson, color = "Ajustado (GLM)")): Adiciona pontos ao gráfico para os valores ajustados pelo modelo GLM, presentes na coluna est_poisson. A cor é definida como “Ajustado (GLM)”.
geom_point(aes(y=est_gam, color = "Ajustado (GAM)")): Adiciona pontos ao gráfico para os valores ajustados pelo modelo GAM, presentes na coluna est_gam. A cor é definida como “Ajustado (GAM)”.
scale_x_date(date_breaks = "1 year", date_labels = "%Y"): Configura o eixo x para exibir marcas de data com intervalos de 1 ano, formatadas para mostrar apenas o ano.
labs(x="Mês/ano", y = "Casos de SRAG", color = ""): Personaliza os rótulos dos eixos x e y, além da legenda de cores (color).

Resumindo, o gráfico resultante mostra a evolução dos casos de SRAG ao longo do tempo, comparando os dados observados e os valores ajustados pelos modelos GLM e GAM.

Agora vemos claramente que os valores ajustados do modelo GAM se aproximam muito mais dos valores observados do que os valores ajustados do modelo GLM. Note que esse modelo GAM ainda é um modelo muito simples, sem incorporação da tendência ou qualquer outro componente, apenas para exemplificar a incorporação de termos não-lineares a um modelo. Vamos prosseguir para uma definição mais completa de um modelo GAM.

Definição de um modelo GAM

Um modelo GAM é, portanto, uma extensão do modelo linear generalizado (GLM) que permite a inclusão destas funções suaves não-lineares para capturar as relações presentes nos dados. Geralmente essas funções são funções de alisamento, como kernel, loess e splines. Neste curso não é nosso objetivo nos aprofundarmos sobre esse tipo de modelo, mas sim apenas contextualizar seu uso dentro do campo de séries temporais. Formalizando, para uma série temporal simples, temos:

η = β_{0} + β_{1} x_{1} + β_{2} x_{2} + \dots + β_{j} x_{j} + f_{1} (x_{j + 1}) + f_{2} (x_{j + 2}) + \dots

$\eta = \beta_0 + \beta_1x_{1} + \beta_2x_{2} + \ldots + \beta_jx_{j} + f_1(x_{j+1}) + f_2(x_{j+2}) + \ldots$

Aqui escrevemos a equação em termos de $\eta$ , que chamamos de preditor linear do modelo. Fazemos isso pois, a depender do tipo de modelo (linear, logístico, Poisson) haverá uma função de ligação entre o parâmetro modelado e a equação. No caso de modelos Poisson, essa função é o log, portanto: $\eta = log(\lambda)$ , onde $\lambda$ é a média da distribuição que costumamente se modela.

Vê-se que agora temos $j$ termos (referentes às variáveis $x_1, \ldots, x_j$ ) modelados de forma linear, como costumamos fazer em um GLM; seguidos de outros termos (referentes a $x_{j+1}, x_{j+2}, \ldots$ ) modelados através de funções $f_1, f_2, \ldots$ que introduzem a não-linearidade à equação.

Essas funções são representadas no R através da função s(), que possui três parâmetros principais:

s(var): o primeiro parâmetro a se especificar é sobre qual(is) variável(is) será aplicada a função de suavização. No nosso exemplo anterior essa variável era o mês de primeiros sintomas: s(mes);
s(var, bs="tp"): bs informa o tipo de função suave a ser usada. Mais detalhes sobre essas funções serão dados abaixo;
s(var, bs="tp", k=6): k é o número de nós, ou a dimensão utilizada para se ajustar o termo. Quanto menos nós, mais simples e parciminiosa nossa função; quanto mais nós, mais rigidamente ela se ajustará aos dados.

Tipos de Splines

Há diversas implementações de funções de suavização, sendo as mais comuns:

thin plate (bs="tp"),
cubic regression (bs="cr"),
cyclic cubic regression (bs="cc"),
B-splines (bs="bs"), e
P-splines, ou splines penalizadas (bs="ps").

De uma forma simplificada:

Thin plate – Baseia-se na obtenção do menor erro quadrático, e penaliza a spline conforme sua curvatura, o que evita complexidade nas curvas e o sobreajuste. Não precisa de muitos parâmetros, mas pode ser mais computacionalmente custosa.

Cubic regression spline – são baseadas em polinômios de grau 3 ajustados aos intervalos dos dados definidos pelos nós distribuídos regularmente em toda a amplitude dos dados. Tende a gerar transições suaves e sem curvas bruscas. No entanto, não se baseia em um estimador ótimo com as thin plate splines.

Cyclic cubic regression spline – semelhantes à spline anterior, há agora a restrição de ter o mesmo valor e mesmas derivadas no início e final, garantindo a ciclicidade. É interessante para séries com termos sazonais.

B-splines: Essas splines baseiam-se em funções-base locais que agem cada uma sobre uma parcela dos dados. Em cada segmento entre os nós especificados são ajustadas funções simples (como as funções cúbicas, mas não necessariamente). Geralmente possuem ajustes eficientes, e menos sensíveis a outliers.

P-splines – Baseiam-se na B-spline, porém aplicam uma penalização diferencial entre os termos adjacentes para controlar possíveis variações bruscas e suavizar a continuidade dos termos. Ideais para dados com variação complexa.

Mais detalhes sobre tais funções e modelos aditivos generalizados podem ser encontrados em Ross (2019).

Voltando ao exemplo das splines para o mês de primeiros sintomas na série de SRAG, vamos comparar as diferentes funções (bs) e dimensões (k) na Figura 46:

Figura 46: Relação entre mês de primeiros sintomas e casos de SRAG no Paraná estimada com diferentes funções e dimensões de modelo GAM.

<
1
2
3
4
5
7
8
Conclusão
>