bayessir: PMMH algorithm for time-varying SIRS model

In vnminin/bayessir: Particle marginal Metropolis-Hastings (PMMH) algorithm for time-varying Susceptible-Infectious-Recovered-Susceptible (SIRS) model

Description

PMMH algorithm for time-varying SIRS model

Usage

bayessir(obscholLIST, obsdaysLIST, COVMATLIST, th, trans, invtrans, pmean, psd,
  betaIIndex, betaWIndex, gammaIndex, kappaIndex, etaIndex, muIndex,
  alphasIndex, rhoIndex, startmeansIndex, nu1Index, nu2Index, nu3Index,
  nu4Index, burn, prelim, iters, thin, tune, ll, psigma, deltavalue, critical,
  PopSize, theMU, numParticles, resultspath, UseGill, UseSIWR, setBetaW,
  setKappa, setEta, setRatio, setval, setAlpha0, setAlpha0val, setstartmeansval,
  usetprior, alphadf, uselaplaceprior, maxWval)

Arguments

obscholLIST	List containing the cholera counts for each phase of data
obsdaysLIST	List containing the observation days for each phase of data
COVMATLIST	List containing the matricies of daily covariates for each phase of data
th	Starting value for parameter vales
trans	function to transform the parameter values
invtrans	inverse transformation function
pmean	Prior means for Normal prior distributions
psd	Prior standard deviations for Normal prior distributions
betaIIndex	Index of which th values correspond to β_I
betaWIndex	Index of which th values correspond to β_W
gammaIndex	Index of which th values correspond to γ
kappaIndex	Index of which th values correspond to κ
etaIndex	Index of which th values correspond to η
muIndex	Index of which th values correspond to μ
alphasIndex	Index of which th values correspond to α parameters
rhoIndex	Vector of length equal to number of phases of data, Index of which th values correspond to ρ, the probability of infected individuals seeking treatment
startmeansIndex	Index of which th values correspond to means of initial distributions for the numbers of susceptible and infected individuals respectively
nu1Index	Index of which th value corresponds to ν_1 power
nu2Index	Index of which th value corresponds to ν_2 power
nu3Index	Index of which th value corresponds to ν_3 power
nu4Index	Index of which th value corresponds to ν_4 power
burn	Number of iterations for burn-in run
prelim	Number of total iterations for preliminary run, preliminary run = burn-in run + secondary run
iters	Number of iterations for final run
thin	Amount to thin the chain; only every thinth iteration of the PMMH algorithm is saved
tune	Tuning parameter for the covariance of the multivariate normal proposal distribution in the final run of the PMMH algorithm
ll	Starting value for log-likelihood
psigma	Standard deviation for independent normal proposal distribution in preliminary run
deltavalue	Initial value to use for tau in tau-leaping algorithm
critical	Critical size for modified tau-leaping algorithm; if the population of a compartment is lower than this number a single step algorithm is used until the population gets above the critical size.
PopSize	Population size
theMU	The rate at which immunity is lost, if setting this value
numParticles	Number of particles
resultspath	File path for results
UseGill	boolian; if 1, uses the gillespie algorithm. If 0, uses tau-leaping algorithm
UseSIWR	boolian; if 1, uses the SIWR model. If 0, uses SIRS model
setBetaW	boolian; if 1, sets BetaW parameter. If 0, estimates BetaW parameter.
setKappa	boolian; if 1, sets Kappa parameter. If 0, estimates Kappa parameter.
setEta	boolian; if 1, sets Eta parameter. If 0, estimates Eta parameter.
setRatio	boolian; if 1, sets ratio of kappa and eta.
setval
setAlpha0	boolian; if 1, sets alpha0 parameter. If 0, estimates alpha0 parameter.
setAlpha0val
setstartmeansval
usetprior
alphadf
uselaplaceprior
maxWval	Upper bound for W compartment. If 0, no bounding.

Value

Posterior samples from the final run of the PMMH algorithm

Also, writes 4 files which are updated every 100th iteration:

1. prelimpmcmctimes.csv: times and acceptance ratios for preliminary PMMH run

2. prelimthmat.csv: preliminary PMMH output

3. FINALpmcmctimes: times and acceptance ratios for final PMMH run

4. FINALthmat.csv: final PMMH output

Examples

## Not run: 

library(bayessir)
############################
## simulate data
############################

SimTimes=seq(0,365*4.5, by=14)

############################
# environmental force of infection
############################

int<- -6
A<-2
sincovAmp<- c(2.1,1.8,2,2.2,2)
wave<-pi/(365/2)
t<-0:(max(SimTimes))

sincov<-sin(wave*t)
allsincov=matrix(NA,nrow=length(t),ncol=length(sincovAmp))
for(i in 1:length(sincovAmp)){
  allsincov[,i]<-sincovAmp[i]*sincov
}

sincov[1:365]=allsincov[1:365,1]
sincov[366:(365*2)]=allsincov[366:(365*2),2]
sincov[(365*2+1):(365*3)]=allsincov[(365*2+1):(365*3),3]
sincov[(365*3+1):(365*4)]=allsincov[(365*3+1):(365*4),4]
sincov[(365*4+1):length(sincov)]=allsincov[(365*4+1):length(sincov),5]

alpha<-exp(int+A*sincov)

###########
pop=10000 #population size

phiS=2900
phiI=84

th1=.5/10000 #beta
th2=0.12 #gamma
th3=.0018 #mu
rho=90/10000 #reporting rate

nu1=1
nu2=1
nu3=0
nu4=0

set.seed(10)
sus0=rpois(1,phiS)
inf0=rpois(1,phiI)
shortstart<-as.matrix(c(sus0,inf0))

allcovs<-enviforce(as.matrix(c(sincov)),SimTimes,c(int,A))

simstates<-matrix(NA,nrow=(length(SimTimes)),ncol=2)
simstates[1,]<-shortstart

for (i in 2:length(SimTimes)){
  simstates[i,]<-inhomoSIRSGillespie(simstates[i-1,],pop,SimTimes[i-1],SimTimes[i]-SimTimes[i-1],
                                     c(th1,th2,th3,nu1,nu2,nu3,nu4),allcovs[[i-1]][,2],allcovs[[i-1]][,1])
}
set.seed(9)
SimData<-c()
for(i in 1:dim(simstates)[1]) SimData[i]<-rbinom(1,simstates[i,2],rho)

##################################################

#####
# data for inference
#####
COVMATLIST=list(as.matrix(sincov))
obscholLIST=list(SimData)
obsdaysLIST=list(SimTimes)

numofcovs=1
#################################################
trans=function(p){
  c(log(p[1]),  #beta
    log(p[2]),  #gamma
    log(p[3]),  #mu
    p[4],       #alpha0
    p[5],       #alpha1
    logit(p[6]))#rho
}

invtrans=function(p){
  c(exp(p[1]),  #beta
    exp(p[2]),  #gamma
    exp(p[3]),  #mu
    p[4],       #alpha0
    p[5],       #alpha1
    expit(p[6]))#rho
}

#prior means
pbetaI=log(1.25e-04)
pgamma=log(.1)
pmu=log(.0009)
palpha0=-8
palphas=rep(0,numofcovs)
prho=logit(.03)

pmean=c(pbetaI,pgamma,pmu,palpha0,palphas,prho)

#prior standard deviations
psd=c(5,  #beta
      .09,#gamma
      .3, #mu
      5,  #alpha0
      5,  #alpha1
      2)  #rho

betaIIndex=1 #need one for each phase of data collection, we only simulated one phase
gammaIndex=2
muIndex=3
alphasIndex=4:5
rhoIndex=6

startmeansIndex=nu1Index=nu2Index=nu3Index=nu4Index=NA

betaWIndex=kappaIndex=etaIndex=NA

# Iterations set small for example purposes; increase for applications
burn = 0
prelim = 10
iters =10
thin =1

tune=1

psigma<-diag(c(0.012, #beta
            0.012, #gamma
            0.180, #mu
            0.120, #alpha0
            0.120, #alpha1
            0.012)) #rho


#start values
#Names of th input are used for the column names in the matrix output
th=c(
  betaI=abs(rnorm(1,th1,th1/3)),
  gamma=abs(rnorm(1,th2,th2/10)),
  mu=abs(rnorm(1,th3,th3/10)),
  alpha0=rnorm(1,int,1),
  alpha1=rnorm(1,A,1),
  rho=abs(rnorm(1,rho,rho)))



resultspath<-getwd()

deltavalue=1
critical=10
numParticles=100
UseGill=0
UseSIWR=0

#Set the population size for inference
PopSize=10000
#Set the rate immunity is lost
theMU=NA #doesn't matter what this is since we are estimating mu in this example


setstartmeansval=list(c(10000*.21,10000*.0015))

setBetaW=setKappa=setEta=setRatio=setAlpha0=0 #don't want to set these right now
setval=setAlpha0val=NA #don't want to set these right now


uset=uselaplace=0 # not using shrinkage priors for alpha parameters
alphadf=5

maxW=50000

ll=-50000


bayessirOUT=bayessir(obscholLIST,obsdaysLIST,COVMATLIST,
                    th,trans,invtrans,pmean,psd,
                    betaIIndex,betaWIndex,gammaIndex,kappaIndex,etaIndex,muIndex,alphasIndex,rhoIndex,startmeansIndex,nu1Index,nu2Index,nu3Index,nu4Index,
                    burn,prelim,iters,thin,tune,ll,psigma,
                    deltavalue,critical,PopSize,theMU,numParticles,resultspath,UseGill,UseSIWR,setBetaW,setKappa,setEta,setRatio,setval,setAlpha0,
                    setAlpha0val,setstartmeansval,uset,alphadf,uselaplace,maxW)

#Output columns are posterior samples for parameters in th, in addition to the log-likelihood and accepted values of the hidden states susT and infT at the final observation time T

#Posterior histograms for parameter values
nvars=dim(bayessirOUT)[2]
par(mfrow=c(1,nvars-3))
for(i in 1:(nvars-3)) hist(bayessirOUT[,i],main="",xlab=colnames(bayessirOUT)[i])

#Trace plots for all output
par(mfrow=c(1,nvars))
for(i in 1:(nvars)) plot(ts(bayessirOUT[,i]),xlab=colnames(bayessirOUT)[i],ylab="")

## End(Not run)

vnminin/bayessir documentation built on May 3, 2019, 6:17 p.m.