# evaluating how well the method reconstructs the true trajectory
# set simulation parameters
niter <- 50
popsize = 500
samp_prob <- 0.25
censusInterval <- 0.05
tmax = 10
b <- 5/popsize
m <- 1
initdist <- c(0.95, 0.05, 0)
trajecs_every <- max(1,floor(niter/50))
resample_prop <- 1
burnin <- niter/2
obstimes <- seq(0, tmax, by = censusInterval)
# set.seed(1834127)
SIRres <- SIRsim(popsize = popsize, initdist = initdist, b = b, mu = m, a=0, tmax = tmax, censusInterval = censusInterval, sampprob = samp_prob, returnX = TRUE)
initdist.shift <- runif(1,-0.01,0.01)
inits <- list(beta.init = b - runif(1,0,b/3),
mu.init = m - runif(1,0,m/3),
alpha.init = 0,
probs.init = samp_prob + runif(1,-0.05,0.05),
initdist.init = initdist + c(initdist.shift,-initdist.shift,0))
# priors
priors <- list(beta.prior = c(5/popsize + runif(1,-0.1/popsize,0.1/popsize), 0.001),
mu.prior = c(0.001, 0.001),
alpha.prior = NULL,
p.prior = c(1,1),
init.prior = c(1,1))
tmax <- max(SIRres$results[,1])
# vectors for parameters
Beta <- vector(length=niter); Beta[1] <- inits$beta.init
Mu <- vector(length = niter); Mu[1] <- inits$mu.init
Alpha <- vector(length = niter); Alpha[1] <- inits$alpha.init
probs <- vector(length = niter); probs[1] <- inits$probs.init
p_initinfec <- vector(length = niter)
R_0 <- vector(length = niter)
suff.stats <- matrix(nrow = niter, ncol = 4); colnames(suff.stats) <- c("num_infec", "num_recov", "beta_suff.stat", "mu_suff.stat")
accepts <- vector(length = niter)
# vectors for parameters of distributions for beta, mu, and p. beta and mu have gamma distributions, p has beta distribution.
beta.prior <- priors$beta.prior
mu.prior <- priors$mu.prior
# alpha.prior <- c(6, 12000)
p.prior <- priors$p.prior
init.prior <- priors$init.prior
# log-likelihood vector
loglik <- vector(length=niter)
# list to store trajectories
trajectories <- vector(mode = "list", length = 1 + niter/trajecs_every)
# observation matrix
W.cur <- as.matrix(data.frame(time = SIRres$results$time, sampled = SIRres$results$Observed, augmented = 0))
# set initial distribution
initdist <- c(1 - W.cur[1,2]/probs[1]/popsize, W.cur[1,2]/probs[1]/popsize, 0)
initdist <- c(0.95, 0.05, 0)
p_initinfec[1] <- initdist[2]
# individual trajectories
X.cur <- initializeX(W = W.cur, b = Beta[1], mu = Mu[1], a = Alpha[1], samp_prob = probs[1], initdist = initdist, censusInterval = censusInterval, tmax = tmax, popsize = popsize)
# count matrix
Xcount.cur <- build_countmat(X = X.cur, popsize = popsize)
trajectories[[1]] <- Xcount.cur
# update observation matrix
W.cur <- updateW(W = W.cur, Xcount = Xcount.cur)
# calculate likelihood for initial trajectory
loglik[1] <- calc_loglike(Xcount = Xcount.cur, tmax = tmax, W = W.cur, b = Beta[1], m = Mu[1], a = Alpha[1], p = probs[1], initdist = initdist, popsize = popsize)
# calculate sufficient statistics for initialization
suff.stats[1,] <- c(numinfec = sum(diff(Xcount.cur[,2], lag = 1)>0),
numrecov = sum(diff(Xcount.cur[,2], lag = 1)<=0),
beta_suffstat = sum(Xcount.cur[1:(nrow(Xcount.cur) - 1),2]*Xcount.cur[1:(nrow(Xcount.cur)-1),3]*diff(Xcount.cur[,1], lag = 1)),
mu_suffstat = sum(Xcount.cur[1:(nrow(Xcount.cur) - 1),2]*diff(Xcount.cur[,1], lag = 1)))
R_0[1] <- Beta[1] * popsize / Mu[1]
writeLines(c(""), paste(paste("augSIR_log",popsize,censusInterval,samp_prob,resample_prop,sep="_"),".txt", sep=""))
start.time <- Sys.time()
# start sampler
for(k in 2:niter){
print(k)
if(k %% 5 == 0){
sink(paste(paste("log",popsize,censusInterval,samp_prob,resample_prop,sep="_"),".txt", sep=""), append=TRUE)
cat(paste("Starting iteration",k,"\n"))
sink()
}
subjects <- sample(unique(X.cur[,2]), floor(popsize*resample_prop), replace=TRUE)
pathirm.cur <- build_irm(Xcount = Xcount.cur, b = Beta[k-1], m = Mu[k-1], a = Alpha[k-1], pop = FALSE)
patheigen.cur <- irm_decomp(pathirm.cur = pathirm.cur)
pop_prob.cur <- pop_prob(Xcount = Xcount.cur, tmax = tmax, b = Beta[k-1], m = Mu[k-1], a = Alpha[k-1], initdist = initdist, popsize = popsize)
accepts.k <- 0
# redraw individuals
for(j in 1:length(subjects)){
# get current path
path.cur <- getpath(X.cur, subjects[j])
# get .other matrices
Xother <- X.cur[X.cur[,2]!=subjects[j],]
Xcount.other <- build_countmat(Xother, popsize - 1)
W.other <- get_W_other(W.cur = W.cur, path = path.cur)
# draw new path
path.new <- draw_path(Xcount = Xcount.other, irm = pathirm.cur, irm.eig = patheigen.cur, W = W.other, p = samp_prob, initdist = initdist, tmax = tmax)
# generate .new objects
X.new <- updateX(X = X.cur, path = path.new, j = subjects[j])
Xcount.new <- update_Xcount(Xcount.other = Xcount.other, path = path.new)
W.new <- updateW(W = W.other, path = path.new)
# check if update is needed to path.irm
if(max(Xcount.new[,2]) == pathirm.cur[4,4,dim(pathirm.cur)[3]]){
new.numinf <- pathirm.cur[4,4,dim(pathirm.cur)[3]]+1
pathirm.cur <- update_irm(irm = pathirm.cur, new.numinf = new.numinf, b = Beta[k-1], m = Mu[k-1], a = Alpha[k-1], popsize = popsize)
patheigen.cur <- update_eigen(patheigen = patheigen.cur, pathirm = pathirm.cur)
}
# Metropolis-Hastings
# population trajectory likelihoods
pop_prob.new <- pop_prob(Xcount = Xcount.new, tmax = tmax, b = b, m = m, a = 0, initdist = initdist, popsize = popsize)
# path likelihoods
path_prob.new <- path_prob(path = path.new, Xcount.other = Xcount.other, pathirm = pathirm.cur, initdist = initdist, tmax = tmax)
path_prob.cur <- path_prob(path = path.cur, Xcount.other = Xcount.other, pathirm = pathirm.cur, initdist = initdist, tmax = tmax)
# compute log acceptance ratio
a.prob <- accept_prob(pop_prob.new = pop_prob.new, pop_prob.cur = pop_prob.cur, path_prob.cur = path_prob.cur, path_prob.new = path_prob.new)
# accept or reject
if(a.prob >=0 || exp(a.prob) > runif(1)) {
X.cur <- X.new
Xcount.cur <- Xcount.new
W.cur <- W.new
pop_prob.cur <- pop_prob.new
accepts.k <- accepts.k + 1
}
}
# save proportion accepted
accepts[k-1] <- mean(accepts.k)
# draw new parameters
probs[k] <- update_prob(W = W.cur, p.prior = p.prior)
# new rate parameters
params.update <- update_rates2(Xcount = Xcount.cur, beta.prior = beta.prior, mu.prior = mu.prior, init.prior = init.prior, popsize = popsize)
params.new <- params.update[[1]]
suff.stats[k, ] <- params.update[[2]]
Beta[k] <- params.new[1]
Mu[k] <- params.new[2]
Alpha[k] <- params.new[3]
p_initinfec[k] <- params.new[5]
initdist <- params.new[4:6]
loglik[k] <- calc_loglike(Xcount = Xcount.cur, tmax = tmax, W = W.cur, b = Beta[k], m = Mu[k], a = Alpha[k], p = probs[k], initdist = initdist, popsize = popsize)
R_0[k] <- Beta[k] * popsize / Mu[k]
if(k%%trajecs_every == 0) {
trajectories[[1 + k/trajecs_every]] <- Xcount.cur
}
}
end.time <- Sys.time()
total.time <- difftime(end.time, start.time)
results <- list(total.time, quantities = cbind(loglik, Beta, Mu, probs, p_initinfec, suff.stats, R_0), trajectories)
assign(paste("augSIR_",popsize,censusInterval,samp_prob,resample_prop,sep="_"),results)
save(list = paste("augSIR_",popsize,censusInterval,samp_prob,resample_prop,sep="_"), file = paste(paste("augSIR_",popsize,censusInterval,samp_prob,resample_prop,sep="_"),".Rdata", sep=""))
# plot stuff
pdf(file = paste(paste("plots_augSIR",popsize,censusInterval,samp_prob,resample_prop,sep="_"),".pdf", sep=""))
par(mfrow = c(5,2))
ts.plot(results[[2]][,"Beta"]); plot(density(results[[2]][,"Beta"]), main = "Beta"); abline(v = b, col = "red")
ts.plot(results[[2]][,"Mu"]); plot(density(results[[2]][,"Mu"]), main = "Mu"); abline(v = m, col = "red")
ts.plot(results[[2]][,"probs"]); plot(density(results[[2]][,"probs"]), main = "Binomial probability"); abline(v = samp_prob, col = "red")
ts.plot(results[[2]][,"p_initinfec"]); plot(density(results[[2]][,"p_initinfec"]), main = "Prob. of infection at time 0"); abline(v = 0.05, col = "red")
ts.plot(results[[2]][,"R_0"]); plot(density(results[[2]][,"R_0"]), main = "R0"); abline(v = 0.05, col = "red")
par(mfrow = c(1,1))
pairs(results[[2]][,c(2:5, ncol(results[[2]]))])
ts.plot(loglik, main = "log-likelihood", xlab = "iteration")
truetrajec <- data.frame(time = unique(SIRres$trajectory[,1]),
infected = c(sum(SIRres$trajectory[SIRres$trajectory[,1]==0,3]), sum(SIRres$trajectory[SIRres$trajectory[,1]==0,3]) + cumsum(SIRres$trajectory[SIRres$trajectory[,1]!=0,3])),
simnum = 0)
trajectories <- results[[3]]
trajectories2 <- list(); observations2 <- list()
for(k in 1:length(results[[3]])){
if ((k%%1)==0){
traj <- trajectories[[k]]
Xobs <- data.frame(time = traj[,1],
infected = traj[,2],
simnum = k)
trajectories2[[k]] <- Xobs
obs <- data.frame(time = seq(0,max(Xobs$time),by=censusInterval),
truth = 0,
sampled = 0,
simnum = k)
for(j in 1:dim(obs)[1]){
obs$truth[j] <- traj[traj[,1] <= obs[j,1],2][sum(traj[,1] <= obs[j,1])]
}
obs$sampled <- rbinom(dim(obs)[1], obs$truth, prob = samp_prob)
observations2[[k]] <- obs
}
}
trajecs <- data.frame(do.call(rbind,trajectories2))
print(ggplot(data = trajecs, aes(x = time, y = infected, group = simnum)) + geom_path(alpha = 0.3) + geom_path(data = truetrajec, aes(x = time, y = infected),colour="red", size = 1) + geom_point(data=data.frame(SIRres$results,simnum=0), aes(x=time,y=Observed),size=4,colour="blue") + theme_bw() + scale_x_continuous(limits = c(0, tmax)))
dev.off()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.