Crítica de “Projeção da dinâmica de transmissão do SARS-CoV-2 durante o período pós-pandemia” – Parte 1: Reproduzindo os resultados

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


[Esteartigofoipublicadopelaprimeiravezem[Thisarticlewasfirstpublishedon Programação R – Blog de Radford Neal, e gentilmente contribuiu para os R-blogueiros]. (Você pode relatar um problema sobre o conteúdo desta página aqui)


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

Eu estive analisando o artigo a seguir, de pesquisadores da escola de saúde pública de Harvard, publicada recentemente em Ciência:

Kissler, Tedijanto, Goldstein, Grad e Lipsitch (2020) Projetando a dinâmica de transmissão do SARS-CoV-2 durante o período pós-pandemia (também disponível aqui, com materiais suplementares aqui).

Este é um dos artigos mencionados no meu recente post sobre sazonalidade do COVID-19. O artigo faz várias coisas que parecem interessantes:

  • Ele analisa a incidência passada de coronavírus “resfriado comum”, estimando os números de reprodução (R) dos vírus ao longo do tempo e, a partir disso, seus graus de imunidade cruzada e o efeito sazonal na transmissão.
  • Ele se encaixa em um modelo ODE para os dois betacoronavírus comuns do resfriado, que estão relacionados ao SARS-CoV-2 (o vírus do COVID-19), usando os mesmos dados.
  • Em seguida, adiciona SARS-CoV-2 a este modelo de ODE e analisa vários cenários para o futuro, variando a duração da imunidade para SARS-CoV-2, o grau de imunidade cruzada de SARS-CoV-2 e betacoronavírus de resfriado comum e o efeito da estação na transmissão de SARS-CoV-2.

Nas próximas postagens, discutirei o conteúdo dessas contribuições. Neste post, falarei sobre meus esforços em reproduzir os resultados no documento a partir do código e dos dados disponíveis, que é um pré-requisito para examinar por que os resultados são como são e para analisar como os métodos usados ​​podem ser aprimorados .

Também falarei sobre um aspecto divertido / horrível do código R usado, que encontrei ao longo do caminho, sobre a política de compartilhamento de dados do CDC e sobre as escolhas dos autores em relação a algumas apresentações gráficas.

Os autores divulgaram alguns dos códigos e dados usados ​​no artigo em três repositórios do github:

  • https://github.com/c2-d2/CoV-seasonality
  • https://github.com/skissler/nCoV_introduction
  • https://github.com/c2-d2/postpandemic (no entanto, isso pode pertencer apenas a uma versão anterior do artigo)

Esses repositórios correspondem aproximadamente (mas incompletamente) às três partes do documento listadas acima. Vou falar sobre a reprodução dos resultados de cada parte, por sua vez.

Estimando e modelando a alteração de R ao longo do tempo para coronavírus comuns do resfriado

O documento usa dados do CDC para estimar como o número de reprodução, R, dos quatro coronavírus “resfriados comuns” mudou ao longo do tempo (do outono de 2014 à primavera de 2019, nos EUA) e usa essas estimativas para R para estimar o impacto da sazonalidade na transmissão desses coronavírus, o grau de imunidade que se desenvolve para eles e o grau de imunidade cruzada entre os dois betacoronavírus (HKU1 e OC43). Como o SARS-CoV-2 também é um betacoronavírus, pode-se esperar que ele se comporte pelo menos de maneira semelhante aos dois betacoronavírus comuns, e que talvez exista alguma imunidade cruzada entre o SARS-CoV-2 e os outros betacoronavírus.

Reproduzindo as estimativas para R

O procedimento usado no artigo para estimar R ao longo do tempo para cada um desses vírus tem várias etapas:

  • A incidência de infecção a cada semana com o coronavírus resfriado comum foi estimada (até um fator de escala desconhecido, relacionado à probabilidade de pessoas doentes visitarem um médico) multiplicando os relatórios semanais de consultas médicas por Influenza-Like Illness (ILI) por a porcentagem semanal de exames laboratoriais para os quatro coronavírus positivos para cada um deles.
  • Um procedimento baseado em spline foi usado para interpolar a incidência diária de cada vírus a partir desses números semanais.
  • A partir desses números de incidência diária, as estimativas de R para cada dia foram obtidas usando uma fórmula que analisa a incidência naquele dia e nos 19 dias anteriores.
  • As estimativas diárias de R foram usadas para produzir estimativas semanais de R, calculando a média geométrica dos 21 valores diários da semana em questão e das semanas anteriores e seguintes.

O primeiro problema com a reprodução dessas estimativas é que, embora os dados das consultas médicas para ILI estejam disponíveis a partir daqui (e incluídos no primeiro repositório acima), o CDC permite acesso apenas ao últimos dois anos de dados sobre testes positivos para coronavírus resfriado comum (daqui em diante). De acordo com o README do repositório, “Dados completos usados ​​em papel estão disponíveis por meio de um contrato de uso de dados com o CDC”.

Esse tipo de besteira faz pensar na mentalidade das pessoas que dirigem o CDC. Obviamente sem motivo por manter esses dados em sigilo. A confidencialidade do paciente não pode ser um problema, devido à natureza dos dados e ao fato de torná-los públicos nos últimos dois anos. Nem pode ser uma questão de minimizar o trabalho por parte do CDC – é necessário um esforço extra para continuar removendo os dados mais antigos, de modo que apenas dois anos estejam disponíveis, sem mencionar o esforço de processar contratos de uso de dados.

Essa política do CDC certamente resultou em trabalho extra para os autores deste artigo. Eles incluíram os últimos dois anos de dados publicamente disponíveis no primeiro repositório acima, juntamente com o código R que foi modificado para funcionar com apenas dois anos de dados em vez de cinco anos. Obviamente, os resultados produzidos não são os mesmos do artigo.

Leia Também  Visualizando componentes principais para imagens

Felizmente, o segundo repositório acima possui um arquivo de dados que de fato inclui os dados completos que foram omitidos no primeiro repositório. Os dados podem ser reformatados para o formulário necessário da seguinte maneira:

dfi <- read.csv("../nCoV_introduction-master/nrevssCDC_ILI.csv",head=TRUE)

dfo < as.data.frame(
         list(RepWeekDate=as.character(as.Date(dfi$WEEKEND),"%m/%d/%y"),
              CoVHKU1=round(100*dfi$HKU1,7), 
              CoVNL63=round(100*dfi$NL63,7),
              CoVOC43=round(100*dfi$OC43,7),
              CoV229E=round(100*dfi$E229,7)))

write.table (dfo, "full-Corona4PP_Nat.csv", sep=",",
             row.names=FALSE, col.names=TRUE, quote=FALSE)

Agora, a tarefa restante é modificar o código R fornecido no primeiro repositório para que funcione com os cinco anos completos de dados. Aqui estão as diferenças cruciais necessárias para fazer isso:

-# Data below shared by NREVSS team
-df.us_cov_national <- read.csv("Corona4PP_Nat.csv") #2018-03-10 through 2020-02-29
+# Reconstruction of full dataset used in paper.
+# Data is for 2014-07-05 through 2019-06-29.
+df.us_cov_national <- read.csv("full-Corona4PP_Nat.csv")

-    Week_start = "2018-07-01") & (Week_start = "2019-07-01") & (Week_start < "2020-07-01") ~ 2))  # 2018-2019 is the last season in our data
+    Week_start < "2014-07-06" ~ 0, # Before first season
+    Week_start < "2015-07-05" ~ 1,
+    Week_start < "2016-07-03" ~ 2,
+    Week_start < "2017-07-02" ~ 3,
+    Week_start < "2018-07-01" ~ 4,
+    Week_start < "2019-06-30" ~ 5, # 2018-2019 is the last season, last data is for 2019-06-29
+    TRUE ~ 0)) # after last season

-for(s in 1:2){
-  temp.df % filter(season==s, epi_week>=season_start | epi_week<=season_end)
+for(s in 1:5){
+  temp.df % filter(season==s, epi_week>=season_start | epi_week%
-  mutate(season=factor(season, levels=c("1", "2"))) #Set season 1 as reference group in regression
-# Note: with this limited dataset, season 2 is incomplete. Full dataset has 5 complete seasons.
+    season==1 ~ "2014-15",
+    season==2 ~ "2015-16",
+    season==3 ~ "2016-17",
+    season==4 ~ "2017-18",
+    season==5 ~ "2018-19")) %>%
+  mutate(season=factor(season, levels=c("1", "2", "3", "4", "5"))) #Set season 1 as reference group in regression

Também adicionei código para produzir várias plotagens e outras saídas, algumas correspondendo a plotagens no papel ou em informações suplementares, e outras para meu uso em descobrir o que o código faz. O código original não vem com uma licença de código-fonte aberto, por isso não publicarei meu arquivo de origem modificado completo, mas parte do código que adicionei no final está aqui e alguns dos gráficos que ele produziu estão aqui e aqui.

Uma digressão sobre o código R

No entanto, vou falar sobre um pequeno trecho do programa original, cujo comportamento é ... interessante:

    RDaily <- numeric()
    for(u in 1:(length(week_list)*7)){ #Estimate for each day
      sumt <- 0
      for(t in u:(u+stop)){ #Look ahead starting at day u through (u+max SI)
        suma <- 0
        for(a in 0:(stop)){ #Calc denominator, from day t back through (t-max SI)
          suma = daily_inc[t-a,v]*func.SI_pull(a, serial_int) + suma
        }
        sumt = (daily_inc[t,v]*func.SI_pull(t-u, serial_int))/suma + sumt
      }
      RDaily[u] = sumt
    }

Esse código calcula estimativas diárias para R (colocando-as em RDaily), usando a seguinte fórmula a partir das informações suplementares:

Observe que o loop para você começa em 1, o loop para t dentro que começa às vocêe o loop para uma dentro que começa em 0 e sobe para Pare (Eumax na fórmula), cujo valor é 19. Para o primeiro acesso a daily_inc, o subscrito t-a será 1, da próxima vez, será 0, depois -1, -2,…, -18. Todos, exceto o primeiro desses valores de índice, parecem estar fora dos limites. Mas o programa é executado sem produzir um erro e produz resultados razoáveis. Como isso pode ser?

Bem, programadores R saberão que índices negativos são permitidos e extrairão todos os itens exceto aqueles identificados pelo subscrito negativo. assim daily_inc[-1,v] criará um vetor longo (números 1819) sem erro. Parece que um erro deve surgir mais tarde, no entanto, quando isso resulta em uma tentativa de armazenar 1819 números em RDaily[u], que tem espaço para apenas um.

Mas, crucialmente, antes que um índice negativo seja usado, há uma tentativa de acessar daily_inc[0,v]. Os programadores de R também podem saber que usar um índice zero não é um erro em R, mesmo que os vetores R sejam indexados a partir de 1 - índices zero são apenas ignorados. (Escrevi anteriormente sobre por que essa é uma má ideia.) Quando o índice subscrito é um único índice zero, ignorá-lo resulta na extração de um vetor de comprimento zero.

Agora, vetores de comprimento zero também parecem o tipo de coisa que levaria a algum tipo de erro mais tarde. Mas R está feliz (por um bom motivo) em multiplicar um vetor de comprimento zero por um escalar, com o resultado sendo outro vetor de comprimento zero. O mesmo vale para adição, portanto, quando t-a é 0, o efeito é que suma no loop mais interno é definido como um vetor de comprimento zero. (Não é o mesmo que 0, que foi inicializado!)

Somente depois de suma foi definido como um vetor de comprimento zero, ele é multiplicado por um vetor de comprimento 1891, acessando daily_inc[-1,v]. R também está feliz em multiplicar um vetor de comprimento zero por um vetor de comprimento maior que um (embora isso seja bastante dúbio), com o resultado sendo um vetor de comprimento zero. assim suma permanece um vetor de comprimento zero para o restante do loop interno, como daily_inc é acessado com índices de -1, -2,…, -18. Após a conclusão desse loop, suma é usado para calcular um termo para adicionar a sumt, com o tratamento da aritmética por R em vetores de comprimento zero, resultando em sumt sendo definido como um vetor de comprimento zero e permanecendo um vetor de comprimento zero, mesmo quando t torna-se grande o suficiente para que os acessos com índices menores que um não sejam mais realizados.

Mas ainda parece que devemos receber um erro! Após o loop t que calcula uma estimativa para R no momento você, essa estimativa é armazenada com a atribuição RDaily[u]= sumt. Desde a sumt é um vetor de comprimento zero, esperamos um erro - obtemos um com código como x = c (10,20); x[2]= numérico () por exemplo (observe que numérico()) cria um vetor numérico de comprimento zero). Agora, o código está realmente estendendo RDaily, em vez de substituir um elemento existente, mas isso não explica a falta de um erro, pois códigos como x = c (10,20); x[3]= numérico () também dá um erro.

O ponto crucial final é que todos esses acessos "fora dos limites" ocorrem no começando do procedimento, quando RDaily é um vetor de comprimento zero. Por nenhuma razão clara, R faz não sinalizar um erro para código como x = numérico (); x[3]= numérico (), mas simplesmente sai x como um vetor de comprimento zero. E assim é neste código, com o resultado que RDaily ainda é um vetor de comprimento zero depois que todas as operações com acesso fora dos limites zero e negativo foram concluídas. Nesse ponto, quando você é 20, um valor sensível para R será calculado e armazenado em RDaily[20]. R estenderá automaticamente RDaily do comprimento zero ao comprimento 20, com os 19 primeiros valores configurados para N / D, e os cálculos posteriores prosseguirão conforme o esperado.

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br

Portanto, no final, o resultado calculado é sensato, com as estimativas para R nos dias para os quais os dados dos 19 dias anteriores não estão disponíveis sendo definidos como N / D, embora por um mecanismo que tenho certeza de que não foi previsto pelo programador. Posteriormente, também existem acessos fora dos limites após o final do vetor, o que também resulta em N / D valores em vez de erros. Todas essas referências fora dos limites podem ser evitadas alterando o loop você do seguinte modo:

for (u in (stop+1):(length(week_list)*7-stop)) { #Estimate for each day

Modelando os efeitos em R de imunidade e sazonalidade

O código produz estimativas de R para cada semana da estação fria para todos os quatro coronavírus, mas a atenção se concentra principalmente nos dois betacoronavírus, HKU1 e OC43. Um modelo de regressão é construído para os valores R desses vírus em termos de um efeito sazonal (modelado como um spline, comum a ambos os vírus) e os efeitos da imunidade da exposição ao mesmo vírus e da imunidade cruzada da exposição ao outro dos dois betacoronavírus (quatro coeficientes). Os efeitos da imunidade podem ser estimados apenas até algum fator de escala desconhecido, com a suposição de que a soma dos números semanais de incidência até algum ponto da temporada é proporcional à fração da população que foi exposta a esse vírus.

Os resultados obtidos correspondem aos coeficientes de regressão na Tabela S1 das informações suplementares do artigo, e alguns gráficos adicionais também são esperados, dados os resultados do artigo.

Os efeitos sazonais e de imunidade ao longo do tempo estão resumidos na Figura 1 do artigo. Aqui está a parte desse valor referente ao HKU1 e à estação fria de 2015-2016:

A curva laranja mostra o efeito sazonal multiplicativo estimado em R (os pontos horizontais são um), a curva vermelha é o efeito estimado em R da imunidade a HKU1 e a curva azul é o efeito estimado da imunidade cruzada em OC43.

Aqui está a minha reprodução desta figura (sem tentar reproduzir as faixas de erro):

Isso parece corresponder perfeitamente ao gráfico no artigo, exceto que o gráfico no artigo mostra apenas 30 semanas, enquanto o modelo é adequado aos dados de 33 semanas, que também é o período de tempo do spline usado para modelar o efeito sazonal. Como se pode ver no meu gráfico, após a semana 30 (na barra vertical em negrito), o efeito sazonal modelado em R aumenta substancialmente. Mas esse recurso do ajuste do modelo não é visível nas figuras do artigo.

Os pesquisadores de Harvard realmente deveriam saber que não deveriam fazer isso. O aumento após a semana 30 que não é mostrado em suas parcelas é contrário à expectativa de que R diminuirá no verão e é uma indicação de que o procedimento de modelagem pode não ser bom. Em particular, depois de ver esse aumento no final da temporada, pode-se perguntar se o aumento acentuado no efeito sazonal no R visto no início da temporada é realmente real ou se é apenas um artefato do seu modelo de spline.

Um modelo ODE para incidência de betacoranavírus

O segundo tópico principal do artigo é a adaptação de um modelo ODE (Equação Diferencial Ordinária) para a incidência dos dois betacoronavirues comuns do resfriado. Os dados utilizados são os mesmos da primeira parte do artigo, mas, em vez de estimar diretamente R a cada momento, é utilizado um modelo subjacente do tipo SEIRS (suscetível a expostos a infectados e recuperados). números podem ser derivados e comparados com os dados.

De acordo com o artigo e as informações complementares, os parâmetros do modelo SEIRS (por exemplo, o grau de variação sazonal e a taxa na qual a imunidade diminui) foram ajustados por um procedimento que combinava a amostragem de hipercubo latino (implementada em R) e a otimização de Nelder-Mead ( implementado no Mathematica). O código para esses procedimentos não foi divulgado, no entanto. Portanto, não é possível reproduzir esta parte do papel.

O segundo repositório acima contém código para executar o modelo SEIRS, com parâmetros configurados para valores fixados no código (presumivelmente para os valores encontrados pelo procedimento de otimização que eles executaram).

Esse modelo SEIRS produz valores para R para cada vírus e ponto no tempo, que podem ser comparados com as estimativas da primeira parte do artigo. Para fazer isso, o código R para esta parte do artigo precisa ler as estimativas para R produzidas pelo código R para a primeira parte. Essas estimativas podem ser escritas da seguinte maneira:

rmv <- c(1:3,nrow(Reff.CoV_ili_x_pos_pct_SARS):(nrow(Reff.CoV_ili_x_pos_pct_SARS)-2))
write.table (Reff.CoV_ili_x_pos_pct_SARS[-rmv,],
             "R_ili_x_pos_pct_SARS.csv", row.names=FALSE, col.names=TRUE, quote=FALSE, sep=",")

Minha versão modificada do figuremaker.R O arquivo de origem R fornecido no segundo repositório acima está aqui. Possui pequenas modificações para ler os dados conforme descrito acima e para permitir a produção de plotagens.

Uma das parcelas produzidas pela execução desse código é uma reprodução exata da Figura 2A no artigo:

Este gráfico mostra a incidência real e simulada dos dois betacoronavírus comuns ao longo de cinco estações frias.

A execução do código também produz o que devem ser reproduções das Figuras 2B e 2C, nas quais os valores de R produzidos pelo modelo SEIRS de melhor ajuste (a curva) são comparados com as estimativas semanais de R da primeira parte do artigo. Mas essas reconstruções não correspondem ao papel. Aqui está a Figura 2B do artigo:

E aqui está o que o código produz:

As curvas são as mesmas (além da escala vertical), mas a figura no artigo está faltando as 12 primeiras estimativas para R, e as primeiras estimativas a seguir são visivelmente diferentes.

Eu descobri que uma reprodução exata das Figuras 2B e 2C no artigo lata ser obtido executando novamente o código para estimar R usando um arquivo de dados no qual as onze primeiras estimativas para R foram excluídas e no qual o código foi alterado para dizer que a primeira temporada começa em 28/09/2014 (em vez disso de 06/07/2014). Aqui está o resultado que combina perfeitamente:

Talvez as Figuras 2B e 2C do artigo tenham sido produzidas inadvertidamente usando um arquivo de estimativas R criadas usando um código preliminar que, por algum motivo, tratou o início da temporada de maneira diferente da versão final do código. É lamentável que o número publicado seja um pouco enganador em relação à correspondência entre as estimativas R da primeira parte do artigo e as estimativas R do modelo SEIRS, pois essa correspondência é significativamente pior para os pontos de dados ausentes do que para os outros.

Projetando o futuro curso da infecção por SARS-CoV-2

A parte final do trabalho estende o modelo ODE para os dois betacoronavírus comuns do resfriado a incluir o SARS-CoV-2, considerando várias possibilidades para as características do SARS-CoV-2, como grau de sazonalidade e duração da imunidade, bem como várias intervenções, como distanciamento social.

Embora a estrutura geral desse modelo estendido esteja documentada nas informações suplementares, o único código relacionado a essas simulações está no terceiro repositório acima, que parece ser uma versão preliminar do artigo. Este código está na forma de um notebook Mathematica (que pode ser visualizado, embora não seja executado, com o programa gratuito aqui). As figuras deste caderno se assemelham às da Figura 3 do artigo, mas não coincidem em detalhes.

Um modelo adicional é usado para modelar cenários relacionados à utilização dos cuidados de saúde e descrito nas informações suplementares. Nenhum código está disponível para este modelo.

Postagens futuras

Este post foi amplamente limitado a descobrir se os resultados no artigo podem ser reproduzidos e, em caso afirmativo, como.

Para a primeira parte do artigo, na qual foram feitas estimativas de R ao longo do tempo e usadas para examinar os efeitos sazonais e de imunidade, consegui reproduzir totalmente os resultados. Para a segunda parte, o modelo SEIRS para infecção por coronavirais comuns, podem ser reproduzidos os resultados para valores de parâmetros especificados, mas o método de otimização usado para encontrar parâmetros de melhor ajuste não é de todo reproduzível a partir das informações fornecidas. A terceira parte, na qual cenários futuros são simulados, também não é reproduzível.

Nas próximas postagens, discutirei os resultados substantivos do artigo, informados sempre que possível por experimentos que posso fazer agora que tenho um código que reproduz alguns dos resultados do artigo. Também considerarei possíveis melhorias nos métodos usados.

var vglnk = {key: '949efb41171ac6ec1bf7f206d57e90b8'};

(função (d, t) {
var s = d.createElement
s.src = '//cdn.viglink.com/api/vglnk.js';
var r = d.getElementsByTagName
} (documento, 'script'));

Para Deixe um comentário para o autor, siga o link e comente no blog: Programação R - Blog de Radford Neal.

R-bloggers.com oferece atualizações diárias por email sobre notícias e tutoriais do R sobre o aprendizado do R e muitos outros tópicos. Clique aqui se você deseja publicar ou encontrar um emprego em ciência da dados / R.


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

*As fotos exibidas neste post pertencem ao post www.r-bloggers.com

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br