Calculating epidemic threshold using MEM package

Lendo dados e carregando pacotes

As séries temporais para semana epidêmica de dengue (periodo de 01/01/2010 até 31/12/2018) e de chikungunya (periodo de 01/01/2017 até 31/12/2018) para o municipio do Rio de Janeiro utilizados nas análises presentes neste documento foram obtidos atráves do InfoDengue. O número de casos foi contabilizado somando as 10 aps do Rio de Janeiro. Para o calculo da incidência foi utilizado a população estimada pelo IBGE para o ano de 2018, 6.688.927 habitantes.

Objetivo

O objetivo desse documento é modelar a incidência de dengue para detectar o início e prever a sua evolução, no que tange a intensidade e a duração das epidemias.

Visualizando a série histórica de dengue

Moving epidemic method (MEM)

Para a definição da limiar epidêmica foi utilizado o moving epidemic method (MEM). Esta metodologia define a atividade inicial da dengue (baseline) baseando-se nos dados históricos e estabelece um limiar epidêmico, sendo que acima desse limiar as taxas semanais são consideradas como um periodo epidêmico. Dessa forma, a intensidade da incidência de dengue pode ser descrita de acordo com as seguintes categorias:

1. Linha base (baseline): taxa semanal ≤ limiar epidemica;

2. Baixa: limiar epidemica < taxa semanal ≤ limiar de intensidade média;

3. Média: limiar de intensidade média < taxa semanal ≤ limiar de intensidade alta;

4. Alta: limiar de intensidade alta < taxa semanal ≤ muito limiar de intensidade alta;

5. Muito alta: taxa semanal > muito limiar de intensidade alta.

Evolução dos estimadores

# evolution of estimators
evolution <- memevolution(dados2,i.evolution.seasons = 10,i.evolution.method ="cross") #Default is 10.
evolution$evolution.seasons
##       2010  2011  2012  2013  2014  2015  2016  2017  2018
## 2010 FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## 2011  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## 2012  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## 2013  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
## 2014  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
## 2015  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
## 2016  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE
## 2017  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE
## 2018  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
## next  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

A evolution.seasons mostra em cada linha um ano e em cada coluna é mostrado se aquele ano especifico foi utilizada no calculo dos indicadores. Se TRUE, então foi usada, se FALSE, então nao foi usada. A ultima linha, chamada de NEXT, os indicadores são calculdados para o ano seguinte.

evolution$evolution.data
##      number durationll duration durationul startll start startul
## 2010      8   8.871287       12   13.34211       1    10      14
## 2011      8   6.871287       12   13.34211       1    10      16
## 2012      8   6.871287       12   12.00000       1    10      16
## 2013      8   6.871287       12   12.00000       1    10      16
## 2014      8   8.742574       12   13.34211       7    10      16
## 2015      8   6.871287       12   13.34211       1    10      14
## 2016      8   6.871287       12   13.34211       1    10      16
## 2017      8   6.871287       12   13.34211       7    10      16
## 2018      8   6.871287       12   13.34211       1    10      17
## next      9   7.455814       12   12.77209       3    10      15
##      percentagell percentage percentageul epidemic postepidemic    medium
## 2010     33.00929   65.46332     79.52526 44.68594     36.33922 14.720649
## 2011     25.10624   55.85102     79.82054 41.97071     36.87755  9.330666
## 2012     25.10624   55.85102     70.29221 22.87881     18.83228  8.777540
## 2013     25.10624   55.85102     70.29221 43.91114     37.47048  9.343711
## 2014     32.71446   65.46332     79.52526 42.61366     38.55813 14.934871
## 2015     25.10624   59.52151     79.52526 44.83965     38.59193 10.823445
## 2016     25.10624   55.85102     79.52526 44.71607     38.27734 10.027859
## 2017     25.10624   65.46332     79.52526 41.97317     38.56566 13.787831
## 2018     25.10624   65.46332     79.52526 44.76817     38.61988 13.271450
## next     27.25430   61.79282     75.47844 45.09744     40.61280 11.788840
##           high veryhigh
## 2010 167.77982 491.8607
## 2011 114.52909 346.9334
## 2012  85.80758 235.0526
## 2013 115.08167 349.1341
## 2014 166.28867 482.4787
## 2015 157.68985 515.2687
## 2016 139.30676 445.7223
## 2017 172.81217 528.3428
## 2018 174.19385 543.5382
## next 147.78916 451.8822

O objeto acima obtido através da função memevolution analisa a evolução dos estimadores por meio de 4 indicadores, são eles: duration:,start:, percentage: e as limiares. As saidas do evolution$evolution.data representam:

  • durationll: limite inferior do IC para a duração média da epidemia;
  • duration: duração média da epidemia;
  • durationul: limite superior do IC para a duração média da epidemia;
  • startll: limite inferior do IC para a duração média do inicio da epidemia;
  • start: duração média do inicio da epidemia;
  • startul: limite superior do IC para a duração média do inicio da epidemia;
  • percentagell: limite inferior do IC para o percentage;
  • percentage: Soma da incidência no periodo epidêmico dividido pela soma total da incidência de todo periodo observado;
  • percentageul: limite superior do IC para o percentage.

e os estimadores das 5 limiares:

  • epidemic: limiar epidemica;
  • postepidemic: limiar epidemica baixa;
  • medium: limiar epidemica média;
  • high: limiar epidemica alta;
  • veryhigh: limiar epidemica muito alta.

Qualidade do ajuste

A função *memgoodness** calcula diferentes indicadores relacionados à qualidade do MEM para detectar epidemias. Isto é feito usando os dados vindos do modelo e usando todos os dados do conjunto de dados originais.

# Goodness of fit
epi.good <- memgoodness(dados2[c(1,5:9)],
                        i.seasons = 6,
                        i.detection.values = seq(2.5, 2.8,0.1),
                        i.goodness.method = "threshold",
                        i.output="~/mem_report_files")
#i.goodness.method = sequential,cross ou threshold
epi.good$results
##                            Weeks                Non-missing weeks 
##                      312.0000000                      312.0000000 
##                   True positives                  False positives 
##                      100.0000000                       44.0000000 
##                   True negatives                  False negatives 
##                      954.0000000                      150.0000000 
##                      Sensitivity                      Specificity 
##                        0.4000000                        0.9559118 
##        Positive predictive value        Negative predictive value 
##                        0.6944444                        0.8641304 
##          Positive likehood ratio          Negative likehood ratio 
##                        9.0727273                        0.6276730 
##                Percent agreement Matthews correlation coefficient 
##                        0.8445513                        0.4458738 
##                    Youdens Index 
##                        0.3559118
#CCM = (TP * TN - FP * FN)/(((TP + FP) * (FN + TN) * (FP + TN) * (TP + FN)(])^(1/2))

Os indicadores são sensibilidade, especificidade, valor preditivo positivo, valor preditivo negativo, acurácia e coeficiente de correlação de Matthews.

Os calculos são feitos a partir de processos iterativos onde cada iteração:

  1. Para uma especifico ano o tempo é calculado visando determinar quais semanas estão inseridas nos periodos pré, pós e período epidêmico. Isto é usado como dado real: os reais positivos (semanas epidemicas) e os reais negativos (pré e pós semana epidemica).

  2. Para cada ano a limiar pré-epidemica é calculada. Essa limiar é comparada com os valores dos anos selecionados no primeiro passo e vê-se os valores acima e abaixo da limiar. Isso é utilizado como dado observado: um valor positivo observado (semana acima da limiar) e um valor negativo observado (semana abaixo da limiar).

  3. Cada semana tem um valor real e observado, assim pode-se classificar-las em:
  • Verdadeiro positivo (VP): positivo real, positivo observado: valores do periodo epidêmico acima da limiar.
  • Verdadeiro negativo (VN): negativo real, negativo observado: valores fora do periodo epidêmico abaixo da limiar.
  • Falso positivo (FP): negativo real, positivo observado: valores fora do periodo epidêmico acima da limiar.
  • Falso negativo (FN): positivo real, negativo observado: valores do periodo epidêmico abaixo da limiar.
  1. O processo é repetido para cada ano no conjunto de dados até todas os anos tenham sido processadas.

  2. Todos os VP, VN, FP and FN são contabilizados e a sensibilidade, especificidade, valor preditivo positivo, valor preditivo negativo, a acurácia e o coeficiente de correlação de Matthews são calculados.

  • Sensibilidade: O número de semanas epidêmicas acima do limiar pré-epidemia (antes do pico) e acima o limiar pós-epidemia (após o pico) dividido por o número de semanas epidêmicas (duração da epidemia).

  • Especificidade: O número de semanas não-epidêmicas abaixo limiar pré-epidemia (antes do pico) e abaixo o limiar pós-epidemia (após o pico) dividido por o número de semanas não-epidêmicas.

  • Valor preditivo positivo (VPP): o número de semanas epidêmicas acima do limiar dividido pelo número de semanas acima do limiar.

  • Valor preditivo negativo (NPV): o número de semanas não-epidêmicas abaixo do limiar dividido pelo número de semanas abaixo do limite.

  • Coeficiente de correlação de Matthew (CCM): Varia de -1 a 1, indicando quão bem o modelo está classificando semanas epidêmicas. Quanto mais próximo a 1, melhor esta sendo a classificação, quanto mais próximo a -1 pior esta sendo a classificação. Assim, o MCC leva em conta os verdadeiros e falsos positivos e negativos

epi.good$peaks
##   Level     Description Count Percentage
## 1     1        Baseline     4  0.6666667
## 2     2             Low     0  0.0000000
## 3     3          Medium     0  0.0000000
## 4     4            High     1  0.1666667
## 5     5       Very high     1  0.1666667
## 6     0 No data seasons     0  0.0000000
## 7    -1   Total seasons     6  1.0000000

O objeto epi.good$peaks mostra a frequencia absoluta e relativa do nível de intensidade dos picos, sendo a cada unidade de medida um ano. Desse modo, as séries de incidência de dengue foram subdivididas por ano e classificados de acordo com as limiares. Os anos 2010, 2014, 2015, 2016, 2017, 2018 foram classificados como baseline; os anos 2011 e 2013 com intensidade média; e apenas 2012 com intensidade alta para a incidência de dengue.

epi.good$peaks.data
##           Peak Peak week Epidemic threshold Medium threshold
## 2010  2.915266        50          10.054669        10.054670
## 2014  2.616264         3           9.248693         9.248694
## 2015 22.484922        20           9.798952         9.798953
## 2016 44.297090        14           4.601117         4.601117
## 2017  4.320573         4           9.110887         9.110888
## 2018  5.217578        20          10.154601        10.154602
##      High threshold Very high threshold Level Description
## 2010       30.98920            66.30387     1    Baseline
## 2014       31.99924            67.60030     1    Baseline
## 2015       19.86320            41.63736     4        High
## 2016       12.85631            22.89212     5   Very high
## 2017       32.60144            73.22303     1    Baseline
## 2018       31.94531            73.33375     1    Baseline

Para cada ano, o objeto epi.good$peaks.data mostra o valor do pico (incidência máxima para aquele ano), qual foi a semana que apresentou pico, e os valores das limiares epidemicas e de intensidade/nível de incidência cada ano de acordo com as limiares.

Modelo MEM

# best <- invisible(lapply(1:6, function(j) lapply(2:ncol(dados2), function(i) memmodel(dados2,i.seasons = i, i.method = j))))
# best2 <- lapply(1:(ncol(dados2)-1), function(i) best %>% lapply('[[', i) %>%
#                   lapply('[[', 3)) %>% data.frame() %>% t() %>% `rownames<-`(1:48) %>%
#   `colnames<-`(c("percentagell","percentage","percentageul")) %>% data.frame() %>%
#   mutate(i.method = rep(1:6,(ncol(dados2)-1)),
#          i.season = sort(rep(2:ncol(dados2),6))
#          )
# best2$weeks <- lapply(1:(ncol(dados2)-1), function(i) best %>% lapply('[[', i) %>%
#                        lapply('[[', 5)) %>% data.frame() %>% t() %>% `rownames<-`(1:48) %>%
#   `colnames<-`(c("weeks")) %>% data.frame()
# best2$indice <- best2$percentage/best2$weeks
# best2[which(best2$indice==max(best2$indice)),]
# rm(best)
# mem model
epi <- memmodel(dados2[c(2:7)],i.method = 2,i.seasons = 6)
# 1 Arithmetic mean and mean confidence interval.
# 2 Geometric mean and mean confidence interval.
# 3 Median and the Nyblom/normal aproximation confidence interval.
# 4 Median and bootstrap confidence interval.
# 5 Arithmetic mean and point confidence interval (standard deviations).
# 6 Geometric mean and point confidence interval (standard deviations)

# Calculates intensity thresholds
intensity <- memintensity(epi)
#print(epi)
summary(epi)
## Call:
## memmodel(i.data = dados2[c(2:7)], i.seasons = 6, i.method = 2)
## 
## Parameters:
##  - General:
##      + Number of seasons restriction: Restricted to 6
##      + Number of seasons used: 6
##      + Seasons used: 2011, 2012, 2013, 2014, 2015, 2016
##      + Number of weeks: 52
##  - Confidence intervals:
##      + Epidemic threshold: Arithmetic mean and its one sided 95% CI using 2*SD
##      + Intensity: Geometric mean and its one sided 40, 90, 97.5% CI using (log) 2*SD
##      + Curve: Geometric mean and its two sided 95% CI using the (log) normal approximation
##      + Others: Median and its two sided 95% CI using the KC Method
##  - Epidemic timing calculation:
##      + Method: 2
##      + Parameter: 2.8
##  - Epidemic threshold calculation:
##      + Pre-epidemic values: Optimized: 5
##      + Tails of CI: 1
##  - Intensity thresholds calculation:
##      + Number of values: Optimized: 5
##      + Tails of CI: 1
##      + Levels of CI:  Medium: 40%, High: 90%, Very high: 97.5% 
##  - Bootstrap (if used):
##      + Technique: norm
##      + Bootstrap samples: 10000
## 
## Epidemic description:
##  - Typical influenza epidemics starts at week 10. 95% CI: [4,14]
##  - Typical influenza season lasts 12 weeks. 95% CI: [8.79,13.64]
##  - This optimal 12 weeks influenza season includes the 69.57% of the total sum of rates
## 
## 
## Epidemic threshold:
##             Pre  Post
## Threshold 44.84 37.67
## 
## 
## Intensity thresholds:
##                   Threshold
## Medium (40%)          25.06
## High (90%)           236.71
## Very high (97.5%)    638.59

Pareando as semanas epidemiológicas de acordo com o periodo epidêmico de dengue obtêm-se o modelo epidemiológico e a limiar para determinar o inicio da epidemia. Uma vez ultrapassado esse limiar, os intervalos de confiança da curva estimada são usados para verificar a intensidade e predizer o fim do periodo epidêmico, monitorando e ajustando o nivel de confiança de acordo com a variância e o tamanho da amostra.

Moving epidemics

A duração do periodo epidêmico obtida pela modelagem tem duração estimada em 12 semanas [8.87,13.34], ao qual contabiliza cerca de 65.46% do total da incidência de dengue.

Average Curve

Performance da estabilidade dos estimadores

# Stability - Function memstability perform an analisys of stability of estimators
stability <- memstability(dados2)
stability$stability.data
##   durationll duration durationul startll start startul percentagell
## 2   9.000000     10.5   12.00000       1     7      13     33.53261
## 3   9.000000     12.0   12.00000       1     9      13     33.53261
## 4   9.000000     12.0   12.00000       1    11      15     33.53261
## 5   7.000000     12.0   12.00000       1     9      15     25.40107
## 6   7.714286     12.0   12.64286       1     8      15     28.30519
## 7   8.466667     12.0   13.26667       1     9      14     31.36420
## 8   8.871287     12.0   13.34211       1    10      14     33.00929
## 9   7.455814     12.0   12.77209       3    10      15     27.25430
##   percentage percentageul  epidemic postepidemic    medium       high
## 2   41.72091     49.90921  2.877821     3.051424  2.875433   4.946559
## 3   49.90921     69.13381 10.351150     8.271990  5.392399  28.360421
## 4   55.85102     69.13381  9.548386     8.993753  7.154377  33.738456
## 5   49.90921     69.13381 10.054669     9.121433  5.544298  30.989198
## 6   55.85102     76.05173 16.300397    14.026405  8.159174  66.336547
## 7   61.79282     77.84275 45.502671    36.877551 12.034562 140.311424
## 8   65.46332     79.52526 44.685935    36.339215 14.720649 167.779821
## 9   61.79282     75.47844 45.097443    40.612796 11.788840 147.789160
##     veryhigh
## 2   6.286898
## 3  59.068784
## 4  66.962274
## 5  66.303869
## 6 167.499083
## 7 415.478996
## 8 491.860743
## 9 451.882203

Mean of the Maximum Percentage - MAP

Ainda vou escrever sobre Mean of the maximum percentage.

Incidência de dengue considerando as limiares

Ainda vou escrever sobre Mean of the maximum percentage. Mas em resumo, é o gráfico de incidência de dengue, mas considerando as thresholds de acordo com o MEM.

Plot Slope

#Plot Slope Farei dia 21/07
# temp <- data.frame(tim$slope.curve[1],tim$slope.curve[2])
# ggplot(data = temp, aes(weeks, slope)) +
#   geom_line(color = "darkgray",size = 1) +
#   geom_point(color="black",size = 3) + 
#   scale_x_continuous(breaks=seq(0,nrow(temp),5)) +
#   scale_y_continuous(breaks=seq(0,5,0.5)) +
#   theme_classic() +
#   scale_fill_manual(values=c("#ffffff")) +
#   theme(plot.title = element_text(size = 14,face = "bold"),
#         text = element_text(size = 12),
#         axis.title = element_text(face="bold"),
#         axis.text.x=element_text(size = 12)) +
#   ggtitle(paste0("Slope curve - ",names(dados2)[i])) + 
#   xlab("Semana epidemiológica") + ylab("Limiar")+
#   geom_hline(yintercept=tim$slope.threshold, linetype="dashed", 
#              color = "purple", size=1.5) +
#   geom_segment(aes(x = 1, y = tim$slope.threshold,
#                    xend = nrow(temp), yend = tim$slope.threshold),
#                color = "orange", size=1) +
#   geom_point(aes(x=tim$optimum.map[1], y=tim$slope.threshold), colour="red",size=4)

O slope method calcula a inclinação da curva MAP, embora o valor ótimo é o que corresponde ao valor global. Este é o método 3, o default do pacote é o método 2, chamado de * fixed criterium method* é uma atualização do MEM que considera a inclinação da curva para encontrar o valor ótimo.

# Methods for influenza trend calculation
# memtrend is used to calculate the two parameters for defining the current influenza trend.

# Calculates trend thresholds
trend <- memtrend(epi)
trend
## $trend.thresholds
##                  Ascending Threshold Descending Threshold
## Trend Thresholds            8.392185            -7.281236
## 
## $param.flu
## Call:
## memmodel(i.data = dados2[c(2:7)], i.seasons = 6, i.method = 2)
## 
## Epidemic threshold:
##             Pre  Post
## Threshold 44.84 37.67
## 
## Intensity thresholds:
##                   Threshold
## Medium (40%)          25.06
## High (90%)           236.71
## Very high (97.5%)    638.59
## 
## $param.type
## [1] 1
## 
## $param.level
## [1] 0.95
## 
## $param.type.boot
## [1] "norm"
## 
## $param.iter.boot
## [1] 10000
## 
## $call
## memtrend(i.flu = epi)
#Full process plots for mem (creates graphs of all mem process)
#processPlots(epi,i.output = "arquivos extras")

Curva ROC

#Analysis of different indicators to find the optimum value of the window parameter
#roc.analysis perform a ROC analysis
# epi.roc <- roc.analysis(dados2)
#epi.roc$roc.data
# epi.roc$optimum
# epi.roc$rankings

Series pelo treshold

ggplot(data = temp, aes(se, valor,color=limiar)) +
  geom_line(color = "darkgray",size = 1) +
  geom_point(size = 2) + 
  theme_classic() +
  scale_color_manual(name="",values = c("#00C000","#d1e802","#e8a702","#e82102"),
                     labels = c("Normal","Alerta","Emergência","Calamidade"),
                     drop=FALSE) +
  theme(plot.title = element_text(size = 14,face = "bold"),
        text = element_text(size = 12),
        axis.title = element_text(face="bold"),
        axis.text.x=element_text(size = 12),
        legend.position="bottom") +
  ggtitle("Incidência por ano e nível de intensidade") +
  xlab("Semana epidemiológica") + ylab("Incidência (por 100mil)") + 
  # geom_hline(yintercept=epi$intensity.thresholds[1], linetype="dashed",
  #            color = "#d1e802", size=0.5) +
  facet_wrap(~ ano)

Avatar
Lucas Bianchi
PhD Student

My research interests include statistics models applied to arbovirus epidemics.

comments powered by Disqus