Validação Deconvolução

Para validação da deconvolução, esta foi aplicada nos dados da amostra de notificados que forneceram informações da dada de início dos sintomas. A análise foi realizada para a amostra de Belo Horizonte, dada a esparcidade das amostras das cidades menores. Inicialmente foi realizada a análise da distribuição do atraso nas notificações.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(fitdistrplus)
## Carregando pacotes exigidos: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Carregando pacotes exigidos: survival
library(utils)
library(httr)
library(ggplot2)
library(EpiEstim)
library(incidental)
library(splines)
library(caret)
## Carregando pacotes exigidos: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:httr':
## 
##     progress
## The following object is masked from 'package:survival':
## 
##     cluster
data <- read.csv("INFLUD-15-02-2021.csv", sep = ";")

##### BH
# Filtrando dados apenas de BH
BH <- filter(data, ID_MUNICIP == "BELO HORIZONTE")

# Armazenando data de notificação
BH_dados <- as.data.frame(BH$ID_REGIONA)
BH_dados[, "DataNot"] <- as.data.frame(BH$DT_NOTIFIC)

# Armazenando data de início dos sintomas
BH_dados[, "DataInicio"] <- as.data.frame(BH$DT_SIN_PRI)

# Passando dados para o formato de data
BH_dados$DataNot <- as.Date(BH_dados$DataNot, "%d/%m/%y")
BH_dados$DataInicio <- as.Date(BH_dados$DataInicio, "%d/%m/%y")

# Criando coluna com atraso nos dias de notificação
BH_dados[, "BH_Atraso"] <- difftime(BH_dados$DataNot, BH_dados$DataInicio, units = "days")

# Armazenado dados (em valor absoluto) do atraso em um vetor
BH_Atraso <- abs(as.vector(BH_dados$BH_Atraso))

# Retirando valores negativos devido ao erro de ano na notificação dos dados
BH_Atraso2 <- BH_Atraso>0
BH_Atraso <- BH_Atraso[BH_Atraso2]

# Retirando valores maiores que 100 (prov?veis de erros de digitação)
BH_Atraso3 <- BH_Atraso<100
BH_Atraso <- BH_Atraso[BH_Atraso3]

Posteriormente foi realizada a deconvolução da amostra de Belo Horizonte.

### deconvolucao
BH_amostra_not <- tabulate(as.factor(BH_dados$DataNot))
BH_amostral <- fit_incidence(reported = BH_amostra_not, 
                             delay_dist = hist(BH_Atraso, 
                                               breaks = 100, 
                                               plot = F)$density,
                             regularization_order = 1,
                             dof_grid = 150)

### plotando dados amostra BH
plot(tabulate(as.factor(BH_dados$DataNot)), pch = 3, col = "red", xlab = "Tempo", 
     ylab = "Indivíduos", ylim = c(0,250))
points(tabulate(as.factor(BH_dados$DataInicio)), pch = 20, col = "black")
legend("topleft", legend = c("Notificados", "Início Sintomas"),
       col = c("red", "black"), pch = c(3,20), cex = 1, bty="n", y.intersp=2)
grid()

Finalmente foi realizada uma regressão via splines para aproximar os dados amostrais de início dos sintomas com o objetivo de averiguar o ajuste da deconvolução. Foi realizada uma validação cruzada para evitar overfitting da spline.

### plotando resultado da deconvolucao
plot(BH_amostral$Ihat, type = "l", col = "black", lty = 1, lwd = 2, xlab = "Tempo", 
     ylab = "Indiv?duos", ylim = c(0,250)) ### deconvolucao (sintomas)
lines(BH_amostral$Chat, col = "red", lty = 2, lwd = 2) ### convolucao (reportados)

data_BH_amostra <- data.frame(tempo = 1:length(tabulate(as.factor(BH_dados$DataInicio))),
                              doentes = tabulate(as.factor(BH_dados$DataInicio)))

set.seed(45)

# dados de treino (80%)
tr <- round(0.8*nrow(data_BH_amostra))
treino <- sample(nrow(data_BH_amostra), tr, replace = F)

# separando dados de treino e teste
dados.treino <- data_BH_amostra[treino,]
dados.teste <- data_BH_amostra[-treino,]

dados.treino <- dados.treino[order(dados.treino$tempo),]
dados.teste <- dados.teste[order(dados.teste$tempo),]

# treino spline
sp11 <- lm(doentes ~ ns(tempo, df = 25), dados.treino)
yhat.sp11 <- sp11$fitted.values
lines(dados.treino$tempo, yhat.sp11, col = "gray50", lwd = 2, lty = 3)
grid()
legend("topleft", legend = c("Deconvolução", "Spline (Notificados)", "Spline (Sintomas)"),
       col = c("black", "red", "gray50"), lty = c(1,2,3), cex = 1, bty="n", y.intersp=2)

# teste spline
yhat.sp11.teste <- predict(sp11, newdata = data.frame(tempo = dados.teste$tempo))

## resultado deconvolucao
res_deconv <- data.frame(obs = BH_amostral$Ihat,
                         pred = tabulate(as.factor(BH_dados$DataInicio))[1:length(BH_amostral$Ihat)])

## resultado spline treino
res_spline <- data.frame(obs = dados.treino$doentes,
                         pred = yhat.sp11)

## resultado spline teste
res_spline.te <- data.frame(obs = dados.teste$doentes,
                            pred = yhat.sp11.teste)

defaultSummary(res_deconv) # deconvolucao
##      RMSE  Rsquared       MAE 
## 21.603078  0.856698 14.475380
defaultSummary(res_spline) # treino spline
##       RMSE   Rsquared        MAE 
## 18.6776685  0.8972583 13.1093169
defaultSummary(res_spline.te) # teste spline
##       RMSE   Rsquared        MAE 
## 18.0856630  0.8802647 12.2104731