I agree with you about the chain wandering into some unrealistic regions. I think the prior might be something I want to look into. My code is very specific but I am sharing it nevertheless.
#function to generate epidemic using weibull as the hazard of contact
epi_weibull<-function(infector,n.pop,shape,rate.contact, rate.recovery)
{
#set.seed(10007)
#infector <- 1; n.pop <- 500; rate.contact<-0.0005; rate.recovery<-0.2
i <- 1
timemat <- result <- NULL
k <- 1:n.pop
susceptible <- k[k!=infector]
time_inf <- 0 # infection time of index case
inc_pd<-rgamma(1,shape=9,scale = 1/3) #incubation period of index case
sym_onset_time_own<- time_inf+inc_pd
IP <- 1/rate.recovery#rexp(1,rate.recovery) #major change here
CI <- rweibull(length(susceptible),shape=shape,scale=1/rate.contact)
contact_idx <- which(CI<=IP)
contact <- susceptible[contact_idx]
contact_time <- sym_onset_time_own + CI[contact_idx] #incubation period is present in contact time
n_contact <- length(contact_idx)
# add the index case to the output data frame.
result <- rbind(result, c(infector, 0, time_inf,sym_onset_time_own, IP, -1,-1, -1, -1))
if(n_contact > 0)
{
# timemat is a matrix containing susceptible individuals who had infectious contact with
# potential infectors.
timemat<-data.frame(infector=rep(infector, n_contact), contact=contact,
recovery_time_infector=rep(sym_onset_time_own + IP, n_contact),
contact_time=contact_time,time_inf_infector=rep(time_inf,n_contact),
sym_onset_infector=rep(sym_onset_time_own, n_contact), IP_infector=rep(IP, n_contact),inc_infector=rep(inc_pd,n_contact))
#contact_time<-round(timemat[,"contact_time"],digits=3)
}
while(!is.null(timemat) && nrow(timemat)>0)
{
i<-i+1
# a new infector is generated
index<-which.min(timemat$contact_time)
infector<-timemat$infector[index]
sym_onset_infector<-timemat$sym_onset_infector[index]
infectee<-timemat$contact[index]
time<-timemat$contact_time[index] #Time of contact
IP_infector<-timemat$IP_infector[index]
inc_infector<-timemat$inc_infector[index]
time_inf_infector=timemat$time_inf_infector[index]
IP <- 1/rate.recovery #rexp(1,rate.recovery) #generate ip for newly infected & merge the ip for previous person
inc_pd<-rgamma(1,shape=9,scale=1/3) #generate incubation period for newly infected
sym_onset_time= time+inc_pd
result<-rbind(result, c(infectee, i, time, sym_onset_time, IP,infector, time_inf_infector,sym_onset_infector, IP_infector))
#update timemat
timemat<-timemat[!(timemat$contact==infectee),] #remove newly infected from matrix
# update the population of susceptibles
susceptible<-susceptible[!(susceptible==infectee)]
if(length(susceptible) >0 )
{
CI <- rweibull(length(susceptible), shape=shape,scale=1/rate.contact) # generate ci with newly infected with susceptible list
contact_idx <- which(CI<=IP)
contact <- susceptible[contact_idx]
contact_time <-sym_onset_time + CI[contact_idx]
n_contact<-length(contact_idx)
if(n_contact > 0)
{
k<- which(!(contact %in% timemat$contact))
m <- length(k)
if(m > 0) timemat.new<-data.frame(infector=rep(infectee, m), contact=contact[k],
recovery_time_infector=rep(sym_onset_time+IP, m),
contact_time=contact_time[k], time_inf_infector= time,
sym_onset_infector=rep(sym_onset_time, m), IP_infector=rep(IP, m),inc_infector=rep(inc_pd,m))
k<- which(contact %in% timemat$contact)
n <- length(k)
if(n > 0)
{
for(j in k)
{
l <- which(timemat$contact == contact[j])
if(contact_time[j] < timemat$contact_time[l])
{
timemat$contact_time[l] <- contact_time[j]
timemat$infector[l] <- infectee
timemat$sym_onset_infector[l] <- sym_onset_time
timemat$recovery_time_infector[l] <- sym_onset_time+IP
timemat$IP_infector[l] <- IP
timemat$inc_infector[l]<-inc_pd
timemat$time_inf_infector<-time
}
}
}
if(!is.null(timemat))
{
if(m > 0) timemat<-rbind(timemat,timemat.new)
} else if(m > 0) timemat<-timemat.new
}
}
}
if(!is.null(result))
{
result <- as.data.frame(result)
colnames(result) <- c('infectee', 'generation', 'timeofinfection','sym_onset_time', 'IP','infector','time_inf_infector','sym_onset_time_infector', 'IP_infector')
}
return(result)
}
#log_likelihood function
log_likelihood_weibull_ind_an=function(shape,contact_rate,time_vec,npop,infector_vec,obs_data,indexcase){
#time is a vector containing times of infection for all infectees
#infector is a vector for infectors for all infectees
l=length(indexcase)
position=match(indexcase,obs_data$infectee)
total_inf= nrow(obs_data)-l
infected_seq= obs_data$infectee[-position]
seq=1:npop
escaped_inf_index= seq[!seq %in% obs_data$infectee]
escaped_inf_length= length(escaped_inf_index)
each_inf= numeric(total_inf)
for (i in 1:length(infected_seq)){
index=which(infected_seq==infected_seq[i]) #an index that runs from 1,2,
time_inf= time_vec[index]
infector= infector_vec[index]
inf_index=which(obs_data$infectee==infector)
sym_onsetinf= obs_data$sym_onset_time[inf_index]
exposed_time=time_inf-sym_onsetinf
escapeindex=which(!obs_data$infectee %in% infector & obs_data$sym_onset_time<time_inf)
escapetime= pmin(time_inf-obs_data$sym_onset_time[escapeindex],obs_data$IP[escapeindex])
inc_per= obs_data$sym_onset_time[index+1]-time_inf
each_inf[i]= -(contact_rate^shape)*exposed_time^shape+shape*log(contact_rate)+log(shape)+
shape*log(exposed_time)-(contact_rate^shape)*(sum(escapetime^shape))+dgamma(inc_per,shape=9,scale=1/3,log = T)
}
return(sum(each_inf)-escaped_inf_length*contact_rate^shape*sum(obs_data$IP))
}
#Running the MCMC
obs_data=epi_weibull(1,1000,5,0.055,0.2)
while(nrow(obs_data)<500) obs_data=epi_weibull(1,1000,5,0.055,0.2)
sd_log_scale=2
accept_shape=0
iter=1000
log_shape=numeric(iter)
log_shape[1]=log(5) #true shape paraneter is 5
for(i in 2:iter){
old_log_shape=log_shape[i-1]
new_log_shape=rnorm(1,mean=old_log_shape,sd=sd_log_scale)
log_lik_top=log_likelihood_weibull_ind_an(contact_rate=0.05,shape=exp(new_log_shape),npop=1000,infector_vec=obs_data$infector[-1],time_vec=obs_data$timeofinfection[-1],indexcase=1,obs_data=obs_data)
log_lik_bottom=log_likelihood_weibull_ind_an(contact_rate = 0.05,shape=exp(old_log_shape),npop=1000,infector_vec=obs_data$infector[-1],time_vec=obs_data$timeofinfection[-1],indexcase=1,obs_data=obs_data)
print(c(old_log_shape,new_log_shape,log_lik_top,log_lik_bottom))
ratio_lik=exp(log_lik_top-log_lik_bottom)
ratio_prior= dnorm(new_log_shape)/dnorm(old_log_shape)
ratio_prop= dnorm(old_log_shape,mean=new_log_shape,sd=sd_log_scale)/dnorm(new_log_shape,mean=old_log_shape,sd=sd_log_scale)
accept_prob= ratio_lik*ratio_prior*ratio_prop
u=runif(1)
if(u<accept_prob){
log_shape[i]=new_log_shape
accept_shape=accept_shape+1
} else log_shape[i]=old_log_shape
}
}