HollTestFit.R

####################################################################################################
## Makes control files for each analysis within which each line giving one R CMD BATCH command line
## to run on a cluster.
####################################################################################################
## Create simulations that fit the Hollingsworth model via MCMC to
## data generated by the same Hollingsworth model.
####################################################################################################
## rm(list=ls())                                  # clear workspace
load("data files/ds.nm.all.Rdata") # country names
load('data files/pars.arr.ac.Rdata')    # load acute phase relative hazards used to fit (in.arr[,,2])
## load('data files/CFJobsToDo.Rdata') ## for finishing up jobs from last run that didn't get finished due to cluster problems.
hazs <- c('bmb','bfb','bme','bfe','bmp','bfp') #  transmission coefficient names, for convenience
nc <- 12                                       # core per simulation
## source('HollTestFit.R')

sb.sim <- FALSE ## don't use SB cohort model, use HOLL model to generate data
blocks <- expand.grid(acute.sc = c(1,2,5,7, seq(10,50,by=5)), late.sc = c(2,5,10),
                      bp = .007, 
                      dur.ac = seq(.5,5, by = .5),
                      dur.lt = c(5,10), dur.aids = 10,  aids.sc=0)
blocks.add <- expand.grid(acute.sc = c(1,2,5,7, seq(10,50,by=5)), late.sc = c(1), ## added late.sc=1, rbind() to keep jobnums the same
                      bp = .007, 
                      dur.ac = seq(.5,5, by = .5),
                      dur.lt = c(5,10), dur.aids = 10,  aids.sc=0)
blocks <- rbind(blocks, blocks.add)

blocks$job <- 1:nrow(blocks)
i.n <- 600
p.n <- 2500
l.n <- 120

####################################################################################################
####################################################################################################
##   MCMC Fitting of Holl Cohort Simulations
####################################################################################################
## Fitting simulated couples data with Wawer & Hollingsworth models to explore potential biases.
####################################################################################################
batchdirnm <- file.path('results','RakAcute','HollModTestFits')
if(!file.exists(batchdirnm))      dir.create(batchdirnm) # create directory if necessary
if(!file.exists(file.path(batchdirnm,'routs')))      dir.create(file.path(batchdirnm, 'routs')) # setup directory to store Rout files
max.vis <- 5
interv <- 10
ltf.prob <- -log(1-.25)/10 ## loss to follow-up probability of each couple, this yields 35% LOSS over 10 months as observed
rr.ltf.ff <- 1.5     ## yields loss to follow-up of 50% for SDCs
rr.ltf.mm <- 1.5     ## yields loss to follow-up of 50% for SDCs
rr.ltf.hh <- 1 ## if ++?
rr.ltf.d <- 1 ## if dead?
rr.inc.sdc <- 1.5 ## how much faster is the rate at which incident SDC (ie who were previously -- in a survey visit) are LTF than regular SDC?

excl.extram <- T ## exclude couples with 2nd partner infected extra-couply?
decont <- F      ## decontaminate? (remove prevalent couples with person-time exposed to acute/late partners etc)
seed.bump <- 0
aniter <- 5000
anburn <- 1000
niter <- 10000
nburn <- 1500
init.jit <- .45
nc <- 12
to.do <- 1:nrow(blocks)
to.do <- which(blocks$late.sc==1)

num.doing <- 0
totn <- 0
sink("HollFitTest.txt")         # create a control file to send to the cluster
## ####################################################################
for(ii in to.do) {
  jb <- ii                   # job num
  totn <- totn+1             # total jobs
  cmd <- paste("R CMD BATCH '--args jobnum=", blocks$job[ii], " simj=", ii, " batchdirnm=\"", batchdirnm, "\"", " nc=", nc,
               " sim.nm=NULL", " simul=T sb.sim=", sb.sim, " i.n=", i.n, " p.n=", p.n, " l.n=", l.n, " nc=12 seed.bump=", seed.bump,
               " interv=", interv, " max.vis=", max.vis, " ltf.prob=", ltf.prob,
               " rr.ltf.ff=", rr.ltf.ff, " rr.ltf.mm=", rr.ltf.mm, " rr.ltf.hh=", rr.ltf.hh, " rr.ltf.d=", rr.ltf.d, " rr.inc.sdc=", rr.inc.sdc,
               " aniter=", aniter, " anburn=", anburn, " niter=", niter, " nburn=", nburn,
               " init.jit=", init.jit, " excl.extram=", excl.extram, " decont=", decont,
               "' RakaiCohortSimulator.R ", file.path(batchdirnm, "routs", paste0('Fit', ii, ".Rout")), sep='')
  num.doing <- num.doing+1
  cat(cmd)               ## add command
  cat('\n')              ## add new line
}
sink()
totn

####################################################################################################
print(totn)
print(num.doing)
save(blocks, file = file.path(batchdirnm, 'blocksHollFitTest.Rdata'))