################################################################################ ################################### Figure 2 ################################### ################################################################################ data <- structure(list(Date = c("2019-11-27", "2019-11-28", "2019-11-29", "2019-11-30", "2019-12-01", "2019-12-02", "2019-12-03", "2019-12-04", "2019-12-05", "2019-12-06", "2019-12-07", "2019-12-08", "2019-12-09", "2019-12-10", "2019-12-11", "2019-12-12", "2019-12-13", "2019-12-14", "2019-12-15", "2019-12-16", "2019-12-17", "2019-12-18", "2019-12-19", "2019-12-20", "2019-12-21", "2019-12-22", "2019-12-23", "2019-12-24", "2019-12-25", "2019-12-26", "2019-12-27", "2019-12-28", "2019-12-29", "2019-12-30", "2019-12-31", "2020-01-01", "2020-01-02", "2020-01-03", "2020-01-04", "2020-01-05", "2020-01-06", "2020-01-07", "2020-01-08", "2020-01-09", "2020-01-10", "2020-01-11", "2020-01-12", "2020-01-13", "2020-01-14", "2020-01-15", "2020-01-16", "2020-01-17", "2020-01-18", "2020-01-19", "2020-01-20", "2020-01-21"), Total = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 2L, 1L, 1L, 0L, 2L, 4L, 4L, 3L, 4L, 0L, 4L, 1L, 3L, 2L, 6L, 5L, 2L, 10L, 7L, 13L, 13L, 15L, 19L, 31L, 44L, 32L, 33L, 33L, 30L, 22L, 23L, 23L, 13L, 6L, 6L, 2L, 3L, 2L), Huanan = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 2L, 3L, 3L, 2L, 3L, 0L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 4L, 1L, 2L, 1L, 1L, 2L, 4L, 1L, 1L, 2L, 0L, 4L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), class = "data.frame", row.names = c(NA, -56L)) data$Date <- as.Date(data$Date) data$non.Huanan <- data$Total - data$Huanan data$Date.numeric <- as.numeric(data$Date - as.Date("2019-11-26")) fit <- glm(non.Huanan ~ Date, data, family = poisson, subset = (Date >= as.Date("2019-12-10") & Date <= as.Date("2020-01-04"))) fit2 <- glm(non.Huanan ~ 1, data, family = poisson, subset = (Date >= as.Date("2019-12-10") & Date <= as.Date("2020-01-04")), offset = I(log(2) / 7.4 * Date.numeric)) pdf("fig2.pdf", width = 8, height = 5) par(mfrow = c(1, 2)) plot(non.Huanan ~ Date, data) lines(data$Date, predict(fit, data, type = "response"), col = "red") lines(data$Date, predict(fit2, data, type = "response"), col = "blue") abline(v = as.Date("2019-12-10"), lty = "dashed") abline(v = as.Date("2020-01-04"), lty = "dashed") plot(log(cumsum(non.Huanan)) ~ Date, data) lines(data$Date, log(cumsum(predict(fit, data, type = "response"))), col = "red") lines(data$Date, log(cumsum(predict(fit2, data, type = "response"))), col = "blue") abline(v = as.Date("2019-12-10"), lty = "dashed") abline(v = as.Date("2020-01-04"), lty = "dashed") dev.off() ## library(R0) ## GT.covid <- generation.time("gamma", c(7.5, 3.4)) ## COVID.Wuhan <- data$non.Huanan ## names(COVID.Wuhan) <- data$Date ## estimate.R(COVID.Wuhan[14:39], GT.covid, methods = c("EG", "ML")) ## str(est.R0.EG(COVID.Wuhan[14:39], GT.covid)) beta <- 3.4^2 / 7.5 alpha <- 7.5 * beta r.to.R0 <- function(r) { 1 / (1 + r / beta)^(- alpha) } r.to.R0(log(2) / 7.4) r.to.R0(fit$coef[2]) ################################################################################ ################################### Figure 3 ################################### ################################################################################ onset <- c(5, 5, 19, 0, 21, 17, 20, 3, 14, 21, 26, 23, 14, 20, 21, 21, 23, 25, 24, 17, 10, 22, 21, 19, 16, 22, 23, 19, 23, 13, 20, 20, 7, 23, 23, 23, 26, 25, 11, 20, 21, 22, 25, 25, 25, 19, 20, 23, 23, 21, 22, 23, 25, 25, 26) data <- data.frame(Date = as.Date("2019-12-25") + (0:(as.Date("2020-01-26") - as.Date("2019-12-25"))), Exported = 0) rownames(data) <- data$Date tab <- table(as.Date("2020-01-01") - 1 + onset) data[names(tab), ]$Exported <- tab rownames(data) <- NULL fit <- glm(Exported ~ Date, data, family = poisson, subset = (Date >= as.Date("2019-12-25") & Date <= as.Date("2020-01-19"))) log(2) / fit$coef[2] fit2 <- glm(Exported ~ Date, data, family = poisson, subset = (Date >= as.Date("2019-12-25") & Date <= as.Date("2020-01-23"))) log(2) / fit2$coef[2] pdf("fig3.pdf", width = 8, height = 5) par(mfrow = c(1, 2)) plot(Exported ~ Date, data) lines(data$Date, predict(fit, data, type = "response"), col = "blue") lines(data$Date, predict(fit2, data, type = "response"), col = "red") abline(v = as.Date("2020-01-19"), lty = "dashed", col = "blue") abline(v = as.Date("2020-01-23"), lty = "dashed", col = "red") plot(log(cumsum(Exported)) ~ Date, data) lines(data$Date, log(cumsum(predict(fit, data, type = "response"))), col = "blue") lines(data$Date, log(cumsum(predict(fit2, data, type = "response"))), col = "red") abline(v = as.Date("2020-01-19"), lty = "dashed", col = "blue") abline(v = as.Date("2020-01-23"), lty = "dashed", col = "red") dev.off() beta <- 3.8^2 / 8.4 alpha <- 8.4 * beta r.to.R0 <- function(r) { 1 / (1 + r / beta)^(- alpha) } r.to.R0(fit$coef[2]) r.to.R0(confint(fit)[2, ]) ################################################################################ ################################### Figure 5 ################################### ################################################################################ library(bets.covid19) data(wuhan_exported) data <- wuhan_exported[1:46, ] T <- 54 # January 23 library(R0) incubation.covid <- generation.time("lognormal", c(5.2, 4.9), truncate = 100)$GT infection.dist <- function(entry) { dist <- rep(0, max(T, entry$S)) dist[1:entry$S] <- incubation.covid[entry$S:1] dist <- dist[1:T] if (entry$B > 1) { dist[1:(entry$B-1)] <- 0 } if (entry$E < T) { dist[(entry$E+1):T] <- 0 } dist / sum(dist) } dist <- do.call(rbind, lapply(1:46, function(i) infection.dist(data[i, ]))) data <- data.frame(Date = as.Date("2019-12-01") + 1:T - 1, Infection = apply(dist, 2, sum)) data <- subset(data, Date >= as.Date("2020-01-01")) fit <- glm(Infection ~ Date, data, family = quasipoisson, subset = (Date <= as.Date("2020-01-20"))) data$Days.before.Jan24 <- as.numeric(as.Date("2020-01-24") - data$Date) fit2 <- glm(I(Infection / Days.before.Jan24) ~ Date, data, family = quasipoisson, subset = (Date <= as.Date("2020-01-20"))) pdf("fig5.pdf", width = 8, height = 5) par(mfrow = c(1, 2)) plot(log(Infection) ~ Date, data) lines(data$Date, predict(fit, data, type = "link"), col = "blue") abline(v = as.Date("2020-01-20"), lty = "dashed", col = "blue") text(x = as.Date("2020-01-13"), y = 1.5, labels = paste("Slope =", signif(fit$coef[2], 2)), col = "blue") plot(log(Infection / Days.before.Jan24) ~ Date, data) lines(data$Date, predict(fit2, data, type = "link"), col = "blue") abline(v = as.Date("2020-01-20"), lty = "dashed", col = "blue") text(x = as.Date("2020-01-14"), y = 0, labels = paste("Slope =", signif(fit2$coef[2], 2)), col = "blue") dev.off() log(2) / fit$coef[2] log(2) / fit2$coef[2] sample.infected <- function() { data2 <- data rownames(data2) <- data2$Date simulated.infections <- sapply(1:nrow(dist), function(i) sample(1:T, 1, prob = dist[i, ])) + as.Date("2019-12-01") - 1 simulated.infections <- simulated.infections[simulated.infections >= as.Date("2020-01-01")] tab <- table(simulated.infections) data2$Infection <- 0 data2[names(tab), ]$Infection <- tab fit2 <- glm(I(Infection / Days.before.Jan24) ~ Date, data2, family = quasipoisson, subset = (Date <= as.Date("2020-01-20"))) fit2$coef[2] } quantile(replicate(1000, sample.infected()), c(0.025, 0.975)) ################################################################################ ################################### Figure 6 ################################### ################################################################################ ## This replicates the analysis in https://github.com/cecilekremer/COVID19/blob/8edca3930e0b68da801fc7b51dfb623232e19a75/MCMC_Singapore_update.R rm(list = ls()) Cluster <- c(rep(1,8), rep(5,23), rep(2,8), rep(3,3), rep(4,5), rep(6,4), rep(7,3)) Case.ID <- c(8,9,31,33,38,83,90,91,66,68,70,71,80,84,88,48,49,51,54,57,58,60,53,61,62,63,67,73,74,78,81,19,20,21,24,25,27,34,40,30,36,39,42,47,52,56,69,2,13,26,44,50,55,77) Case <- c(1:length(Case.ID)) Time <- c(4,4,3,10,14,8,20,3,9,10,14,12,15,15,27,12,14,15,21,22,21,19,21,17,20,21,20,20,23,20,27,9,5,13,10,4,12,7,10,1,4,9,12,17,18,23,22,1,8,8,11,18,10,21) vi <- c(0,0,rep(NA,6), rep(NA,18),24,18,rep(NA,3),0,0,32,0,0,32,rep(NA,8),44,NA,0,0,0,NA,NA,0,52) data <- data.frame(Cluster,Case,Case.ID,Time,vi) NCases <- length(Case) # Which cases are possible? vi.list <- list() vi.list[[3]] <- c(1,2,4:8) vi.list[[4]] <- c(1,2,3,5:8) vi.list[[5]] <- c(1:4,6,7,8) vi.list[[6]] <- c(1,2) vi.list[[7]] <- c(1:6,8) vi.list[[8]] <- c(1,2) vi.list[[9]] = vi.list[[10]] = vi.list[[11]] = vi.list[[12]] = vi.list[[13]] = vi.list[[14]] = vi.list[[15]] <- c(6,8) vi.list[[16]] <- c(9, 17:22) vi.list[[17]] <- c(9, 16,18:22) vi.list[[18]] <- c(9, 16,17,19:22) vi.list[[19]] <- c(9, 16:18,20:22) vi.list[[20]] <- c(9, 16:19,21,22) vi.list[[21]] <- c(9, 16:20,22) vi.list[[22]] <- c(9, 16:21) vi.list[[23]] <- c(9, 16:22, 24:31) vi.list[[24]] <- c(9, 27) vi.list[[25]] <- c(9, 16:22, 23,24,26:31) vi.list[[26]] <- c(9, 16:22, 23:25,27:31) vi.list[[29]] <- c(9, 16:22, 23:28,30,31) vi.list[[30]] <- c(9, 16:22,23:29,31) vi.list[[31]] <- c(9, 16:22, 23:30) vi.list[[38]] <- c(32,33) vi.list[[39]] <- c(32,33) #vi.list[[40]] <- c(41,42,0) vi.list[[40]] <- c(0) # for SI > -3 vi.list[[41]] <- c(40,42,0) vi.list[[42]] <- c(40,41,0) vi.list[[43]] <- c(44,45,0) vi.list[[44]] <- c(43,45,46,0) vi.list[[45]] <- c(43,44,0) vi.list[[47]] <- c(43:46,0) vi.list[[51]] <- c(49,50) vi.list[[52]] <- c(53,54) vi.list[[53]] <- c(0) vi.list[[54]] <- c(52) for(i in 1:NCases){ if(is.null(vi.list[[i]])) vi.list[[i]] <- vi[i] } PossibleInfector <- matrix(nrow=NCases,ncol=max(lengths(vi.list))) for(i in 1:NCases){ PossibleInfector[i,1:lengths(vi.list)[i]] <- vi.list[[i]] } NPossibleInfector <- rowSums(!is.na(PossibleInfector)) IsNotContributorToLikelorg <- c(which(vi==0)) # known index cases IsContributorToLikelorg <- Case[!Case%in%IsNotContributorToLikelorg] ## Constrain neg SI to be max 3 days x = -3 inf.mat = matrix(nrow=NCases, ncol=NCases) for(i in 1:NCases){ for(j in 1:NCases){ inf.mat[i,j] = ifelse(j%in%PossibleInfector[i,],1,0) } } # all possible SI time.mat = matrix(nrow=NCases, ncol=NCases) for(i in 1:NCases){ for(j in 1:NCases){ time.mat[i,j] = Time[i] - Time[j] } } SI.mat = inf.mat*time.mat # min(SI.mat) = -17 (lowest SI) # weigh by number of infectors w.SI <- NULL for(i in IsContributorToLikelorg){ w.SI[i] = sum(SI.mat[i,]) / NPossibleInfector[i] } # Select SI >= x for(i in 1:NCases){ for(j in 1:NCases){ inf.mat[i,j]=ifelse(time.mat[i,j] < x, 0, inf.mat[i,j]) } } SI.mat2 = inf.mat*time.mat # weigh by number of infectors w.SI2 <- NULL for(i in IsContributorToLikelorg){ w.SI2[i] = sum(SI.mat2[i,]) / sum(inf.mat[i,]) } # change vi.list vi.list2 <- rep(list(NULL),NCases) #k=1 for(i in Case){ for(j in Case){ if(inf.mat[i,j]==1) vi.list2[[i]] <- c(vi.list2[[i]],Case[j]) # k=k+1 } } for(i in 1:NCases){ if(is.null(vi.list2[[i]])) vi.list2[[i]]=0 } PossibleInfector <- matrix(nrow=NCases,ncol=max(lengths(vi.list2))) for(i in 1:NCases){ PossibleInfector[i,1:lengths(vi.list2)[i]] <- vi.list2[[i]] } NPossibleInfector <- rowSums(!is.na(PossibleInfector)) IsNotContributorToLikelorg <- c(which(vi==0 | Case==40 | Case==43)) # when allowing SI >= -3 IsContributorToLikelorg <- Case[!Case%in%IsNotContributorToLikelorg] library(igraph) # Function to find cycles in network FindCycles = function(g) { Cycles = NULL for(v1 in V(g)) { if(degree(g, v1, mode="in") == 0) { next } GoodNeighbors = neighbors(g, v1, mode="out") GoodNeighbors = GoodNeighbors[GoodNeighbors > v1] for(v2 in GoodNeighbors) { TempCyc = lapply(all_simple_paths(g, v2,v1, mode="out"), function(p) c(v1,p)) TempCyc = TempCyc[sapply(TempCyc, min) == sapply(TempCyc, `[`, 1)] Cycles = c(Cycles, TempCyc) } } Cycles } ##############-----MCMC----------------------- # likelihood likelihood <- function(){ alpha<-c(); Beta<-c() # alpha[1] & Beta[1] are shape & rate of the serial interval distriution alpha[1] <- theta[1]^2/theta[2] # shape = mean^2/var Beta[1] <- theta[1]/theta[2]; # rate = mean/var # alpha[2] & Beta[2] are shape & rate of the distribution of sum of 2 independent gamma variables alpha[2] <- theta[3]^2/theta[4] # shape = mean^2/var Beta[2] <- theta[3]/theta[4] # rate = mean/var # evaluate density f_Z <- function(Z){ monteCarloN <- 300 delta_i <- rgamma(monteCarloN, shape=alpha[2], rate=Beta[2]) delta_j <- rgamma(monteCarloN, shape=alpha[2], rate=Beta[2]) Y <- delta_i-delta_j Z_Y <- Z - Y return(mean(dgamma(Z_Y, shape=alpha[1], rate=Beta[1], log=F))) } SerialInterval <- Time[IsContributorToLikel] - Time[Network[IsContributorToLikel]] # serial interval return(sum(log(1e-50+sapply(SerialInterval, function(x) f_Z(x))))) } # prior prior <- function(){ alpha<-c(); Beta<-c() alpha[1] <- theta[1]^2/theta[2] # shape = mean^2/var Beta[1] <- theta[1]/theta[2]; # rate = mean/var alpha.prior <- dunif(alpha[1], 0, 30) Beta.prior <- dunif(Beta[1], 0, 20) return(log(1e-50+alpha.prior)+log(1e-50+Beta.prior)) } # posterior posterior <- function(){ return (likelihood()+prior()) } #-------------------------------------------------- #----MCMC--------------------- #-------------------------------------------------- Network <- numeric(NCases)+0 Update <- IsContributorToLikelorg # Resample network until there are no cylces n.cycles=1 while(n.cycles>0){ Draw <- round(runif(length(Update),min=0.5,max=NPossibleInfector[Update]+0.5)) for(i in 1:length(IsContributorToLikelorg)){ Network[Update[i]] <- PossibleInfector[Update[i],Draw[i]] } AcceptedNetwork <- Network Infector=Network[IsContributorToLikelorg] Infectee=Case[IsContributorToLikelorg] tree=cbind(Infector,Infectee); tree=tree[which(Infector!=0),] g = graph_from_edgelist(tree) n.cycles = length(FindCycles(g)) } IsContributorToLikel <- Case[Network!=0] IsNotContributorToLikel <- Case[!Case%in%IsContributorToLikel] AcceptedTheta=theta <- c(1, 1, 5.2, 2.8^2) # incubation: mean=5.2 days & sd=2.8 days SerialInterval <- Time[IsContributorToLikel] - Time[Network[IsContributorToLikel]] # serial interval P <- posterior() AcceptedP <- P NRuns <- 3000000 NUpdate <- length(IsContributorToLikelorg) Burnin <- 500000 Thinning<-200 SaveP <- numeric() SaveNetwork <- matrix(nrow=NCases,ncol=(NRuns-Burnin)/Thinning) Savetheta <- matrix(nrow=(NRuns-Burnin)/Thinning,ncol=(2+length(theta))) # theta + var SI + mean SI anetwork=asd<-0 tuning <- c(0.5, 0.5) a <- 0 sd = 0.5 for(b in 1:NRuns){ if(b%%100 == 0){ theta <- AcceptedTheta Update <- IsContributorToLikelorg n.cycles=1 while(n.cycles>0){ Draw <- round(runif(length(Update),min=0.5,max=NPossibleInfector[Update]+0.5)) for(i in 1:length(IsContributorToLikelorg)){ Network[Update[i]] <- PossibleInfector[Update[i],Draw[i]] } AcceptedNetwork <- Network Infector=Network[IsContributorToLikelorg] Infectee=Case[IsContributorToLikelorg] tree=cbind(Infector,Infectee); tree=tree[which(Infector!=0),] g = graph_from_edgelist(tree) n.cycles = length(FindCycles(g)) } IsContributorToLikel <- Case[Network!=0] IsNotContributorToLikel <- Case[!Case%in%IsContributorToLikel] meanSI = mean(Time[IsContributorToLikel] - Time[Network[IsContributorToLikel]]) } if(b%%100 != 0){ Network <- AcceptedNetwork theta[1] <- runif(1, (AcceptedTheta[1]-tuning[1]), (AcceptedTheta[1]+tuning[1])) if(theta[1]<0){theta[1] <- AcceptedTheta[1]} theta[2] <- runif(1, (AcceptedTheta[2]-tuning[2]), (AcceptedTheta[2]+tuning[2])) if(theta[2]<0){theta[2] <- AcceptedTheta[2]} } P <- posterior() AcceptYN <- runif(1,min=0,max=1) <= exp(P-AcceptedP) if(AcceptYN){ if(b%%100 == 0){ anetwork <- anetwork + 1 AcceptedNetwork <- Network } if(b%%100 != 0){ asd <- asd+1 AcceptedTheta <- theta } AcceptedP <- P } if(b%%Thinning == 0 & b>Burnin){ a <- a + 1 Savetheta[a,] <- c(AcceptedTheta, (AcceptedTheta[2]+2*AcceptedTheta[4]), meanSI) SaveNetwork[,a] <- AcceptedNetwork SaveP[a] <- AcceptedP } } save.image("ganyani20.rda") data <- data.frame(mean = Savetheta[, 1], sd = sqrt(Savetheta[, 2])) library(ggplot2) pl <- ggplot(data) + aes(x = mean, y = sd) + geom_point(data = data[sample(1:nrow(data), 100), ], alpha = 0.7) + geom_density_2d(size = 0.5, alpha = 0.7, color = "grey50") + theme_bw(base_size = 15) points <- data.frame(x = c(4, 6), y = c(2, 4), col = c("red", "blue")) pl <- pl + geom_point(aes(x = x, y = y, color = col), data = points, shape = 4, size = 5, stroke = 2) ## + geom_text(aes(x = x, y = y - 0.4, color = col, label = paste0("(", x, ",", y, ")")), data = points, size = 6, fontface = "bold") pl <- pl + xlab("Mean") + ylab("Standard deviation") + theme(legend.position = "none") ## ggsave("generation-time.pdf", pl) r.to.R0 <- function(r, mean, sd) { beta <- sd^2 / mean alpha <- mean * beta 1 / (1 + r / beta)^(- alpha) } r.to.R0(log(2) / 2.5, 4, 2) r.to.R0(log(2) / 2.5, 6, 4) r <- seq(0.05, 0.42, 0.01) data <- rbind(data.frame(r = r, mean = 4, sd = 2, R0 = r.to.R0(r, 4, 2)), data.frame(r = r, mean = 6, sd = 4, R0 = r.to.R0(r, 6, 4))) data$doubling.time <- log(2) / data$r data$parameter <- paste0("mean = ", data$mean, ", sd = ", data$sd) dt <- data.frame(doubling.time = c(1.8, 2, 2.5, 3, 4, 6)) pl2 <- ggplot(data) + aes(x = r, y = R0) + geom_line(aes(color = parameter)) + geom_vline(aes(xintercept = log(2) / doubling.time), data = dt, linetype = "dashed", color = "gray50") + geom_text(aes(x = log(2) / doubling.time + 0.01, y = 9.5, label = paste(doubling.time)), data = dt, size = 5) + geom_text(aes(x = 0.12, y = 10, label = "Doubling days = "), size = 5) + theme_bw(base_size = 15) + theme(panel.grid.major.x = element_blank(), legend.title = element_blank(), legend.position = "bottom", legend.text = element_text(size = 20)) + ylab(expression(R[0])) + xlab(expression(r)) + expand_limits(y = 1) legend2 <- get_legend(pl2) library(cowplot) pdf("fig6.pdf", width = 10, height = 6) pll <- plot_grid(pl, pl2 + theme(legend.position = "none"), nrow = 1, align = "h", axis = "b") plot_grid(pll, legend2, ncol = 1, rel_heights = c(1, .1)) dev.off()