diff --git a/NAMESPACE b/NAMESPACE index e6cb361..439552f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,7 +6,6 @@ export(Generate_Multiple_Datasets) export(Generate_Sero_Dataset) export(Generate_VIMC_Burden_Dataset) export(MCMC) -export(MCMC2) export(Model_Run) export(Model_Run_Many_Reps) export(Model_Run_Multi_Input) @@ -22,7 +21,6 @@ export(input_data_truncate) export(mcmc_FOI_R0_setup) export(mcmc_checks) export(mcmc_prelim_fit) -export(mcmc_prelim_fit2) export(param_calc_enviro) export(param_prop_setup) export(parameter_setup) @@ -31,7 +29,6 @@ export(sero_calculate) export(sero_calculate2) export(sero_data_compare) export(single_posterior_calc) -export(single_posterior_calc2) export(template_region_xref) import(dde) import(odin.dust) diff --git a/R/main.R b/R/main.R index bec64ef..6894b46 100644 --- a/R/main.R +++ b/R/main.R @@ -562,18 +562,18 @@ parameter_setup <- function(FOI_spillover=0.0,R0=1.0,vacc_data=list(),pop_data=l #' @param input_data List of population and vaccination data for multiple regions (created using data input creation #' code and usually loaded from RDS file) #' @param enviro_data Environmental data frame, containing only relevant environmental variables -#' @param extra_params Vector of strings listing parameters besides ones determining FOI/R0 (may include vaccine -#' efficacy and/or infection/death reporting probabilities) +#' @param extra_estimated_params Vector of strings listing variable parameters besides ones determining FOI/R0 (may include +#' vaccine efficacy and/or infection/death reporting probabilities and/or Brazil FOI adjustment factor) #' #' @export #' -create_param_labels <- function(type="FOI",input_data=list(),enviro_data=NULL,extra_params=c("vacc_eff")){ +create_param_labels <- function(type="FOI",input_data=list(),enviro_data=NULL,extra_estimated_params=c("vacc_eff")){ assert_that(type %in% c("FOI","FOI+R0","FOI enviro","FOI+R0 enviro")) assert_that(input_data_check(input_data),msg=paste("Input data must be in standard format", " (see https://mrc-ide.github.io/YEP/articles/CGuideAInputs.html)")) - n_extra=length(extra_params) + n_extra=length(extra_estimated_params) if(type %in% c("FOI","FOI+R0")){ regions=input_data$region_labels @@ -595,7 +595,7 @@ create_param_labels <- function(type="FOI",input_data=list(),enviro_data=NULL,ex if(type=="FOI+R0 enviro"){param_names[i+n_env_vars]=paste("R0_",env_vars[i],sep="")} } } - if(n_extra>0){param_names[(n_params-n_extra+1):n_params]=extra_params} + if(n_extra>0){param_names[(n_params-n_extra+1):n_params]=extra_estimated_params} return(param_names) } diff --git a/R/mcmc.R b/R/mcmc.R index fbbb4bf..3a2a2a2 100644 --- a/R/mcmc.R +++ b/R/mcmc.R @@ -35,32 +35,30 @@ #' @param Niter Total number of steps to run #' @param type Type of parameter set (FOI only, FOI+R0, FOI and/or R0 coefficients associated with environmental #' covariates); choose from "FOI","FOI+R0","FOI enviro","FOI+R0 enviro" -#' @param log_params_min Lower limits of varied parameter values if specified -#' @param log_params_max Upper limits of varied parameter values if specified +#' @param prior_settings List containing settings for priors: must contain text named "type": +#' If type = "zero", prior probability is always zero +#' If type = "flat", prior probability is zero if log parameter values in designated ranges log_params_min and log_params_max, +#' -Inf otherwise; log_params_min and log_params_max included in prior_settings as vectors of same length as log_params_ini +#' If type = "exp", prior probability is given by dexp calculation on FOI/R0 values +#' If type = "norm", prior probability is given by dnorm calculation on parameter values with settings based on vectors of values +#' in prior_settings; norm_params_mean and norm_params_sd (mean and standard deviation values applied to log FOI/R0 +#' parameters and to actual values of additional parameters) #' @param mode_start Flag indicating how to set initial population immunity level in addition to vaccination #' If mode_start=0, only vaccinated individuals #' If mode_start=1, shift some non-vaccinated individuals into recovered to give herd immunity (uniform by age, R0 based only) #' If mode_start=3, shift some non-vaccinated individuals into recovered to give herd immunity (stratified by age) -#' @param prior_type Text indicating which type of calculation to use for prior probability -#' If prior_type = "zero", prior probability is always zero -#' If prior_type = "flat", prior probability is zero if FOI/R0 in designated ranges, -Inf otherwise -#' If prior_type = "exp", prior probability is given by dexp calculation on FOI/R0 values -#' If prior_type = "norm", prior probability is given by dnorm calculation on parameter values using norm_prior_sd values -#' @param norm_prior_sd Vector of 3 values of standard deviations to use in dnorm calculations -#' Value 1 = SD to use for FOI values/coefficients if prior_type = "norm" -#' Value 2 = SD to use for R0 values/coefficients if prior_type = "norm" (if relevant) -#' Value 3 = SD to use for reporting probability prior calculations for all prior_type values (if relevant) #' @param dt time increment in days (must be 1 or 5) #' @param n_reps Number of times to repeat calculations to get average likelihood at each step #' @param enviro_data Data frame containing values of environmental covariates; set to NULL if not in use #' @param R0_fixed_values Values of R0 to use if only FOI is subject to fitting (i.e. type set to "FOI" or "FOI #' enviro"); set to NULL if not in use -#' @param vaccine_efficacy Vaccine efficacy (set to NULL if being varied as a parameter) #' @param p_severe_inf Probability of an infection being severe #' @param p_death_severe_inf Probability of a severe infection resulting in death -#' @param p_rep_severe Probability of observation of severe infection (set to NULL if being varied as a parameter) -#' @param p_rep_death Probability of observation of death (set to NULL if being varied as a parameter) -#' @param m_FOI_Brazil Multiplier of spillover FOI for Brazil regions (set to NULL if being varied as a parameter) +#' @param add_values List containing values of additional inputs or NULL to indicate they are being varied as parameters: +#' vaccine_efficacy Vaccine efficacy (set to NULL if being varied as a parameter) +#' p_rep_severe Probability of observation of severe infection (set to NULL if being varied as a parameter) +#' p_rep_death Probability of observation of death (set to NULL if being varied as a parameter) +#' m_FOI_Brazil Multiplier of spillover FOI for Brazil regions (set to NULL if being varied as a parameter) #' @param deterministic TRUE/FALSE - set model to run in deterministic mode if TRUE #' @param mode_parallel Set mode for parallelization, if any: #' If mode_parallel="none", no parallelization @@ -71,20 +69,19 @@ #' @export #' MCMC <- function(log_params_ini=c(),input_data=list(),obs_sero_data=NULL,obs_case_data=NULL,filename_prefix="Chain", - Niter=1,type=NULL,log_params_min=c(),log_params_max=c(),mode_start=0,prior_type="zero", - norm_prior_sd=c(30,10,30),dt=1.0,n_reps=1,enviro_data=NULL,R0_fixed_values=NULL, - vaccine_efficacy=NULL,p_severe_inf=0.12,p_death_severe_inf=0.39,p_rep_severe=NULL,p_rep_death=NULL, - m_FOI_Brazil=1.0,deterministic=FALSE,mode_parallel="none",cluster=NULL){ + Niter=1,type=NULL,mode_start=0,prior_settings=list(type="zero"),dt=1.0,n_reps=1,enviro_data=NULL, + R0_fixed_values=NULL,p_severe_inf=0.12,p_death_severe_inf=0.39, + add_values=list(vaccine_efficacy=1.0,p_rep_severe=1.0,p_rep_death=1.0,m_FOI_Brazil=1.0), + deterministic=FALSE,mode_parallel="none",cluster=NULL){ assert_that(is.logical(deterministic)) assert_that(mode_start %in% c(0,1,3),msg="mode_start must have value 0, 1 or 3") n_params=length(log_params_ini) - extra_params=c() - if(is.null(vaccine_efficacy)==TRUE){extra_params=append(extra_params,"vaccine_efficacy")} - if(is.null(p_rep_severe)==TRUE){extra_params=append(extra_params,"p_rep_severe")} - if(is.null(p_rep_death)==TRUE){extra_params=append(extra_params,"p_rep_death")} - if(is.null(m_FOI_Brazil)==TRUE){extra_params=append(extra_params,"m_FOI_Brazil")} + extra_estimated_params=c() + for(var_name in names(add_values)){ + if(is.null(add_values[[var_name]])==TRUE){extra_estimated_params=append(extra_estimated_params,var_name)} + } #Process input data to check that all regions with sero and/or case data supplied are present, remove #regions without any supplied data, and add cross-referencing tables for use when calculating likelihood. Take @@ -98,22 +95,19 @@ MCMC <- function(log_params_ini=c(),input_data=list(),obs_sero_data=NULL,obs_cas } #Label parameters according to order and fitting type - param_names=create_param_labels(type,input_data,enviro_data,extra_params) + param_names=create_param_labels(type,input_data,enviro_data,extra_estimated_params) names(log_params_ini)=param_names - #Run checks to ensure that number of parameters is correct for fitting type and number of regions/environmental - #covariates - checks<-mcmc_checks(log_params_ini,n_regions,type,log_params_min,log_params_max,prior_type,norm_prior_sd, - enviro_data,R0_fixed_values,vaccine_efficacy,p_rep_severe,p_rep_death,m_FOI_Brazil) - if(prior_type=="flat"){names(log_params_min)=names(log_params_max)=param_names} + #Run checks on inputs + checks<-mcmc_checks(log_params_ini,n_regions,type,prior_settings,enviro_data,R0_fixed_values,add_values, + extra_estimated_params) + if(prior_settings$type=="flat"){names(prior_settings$log_params_min)=names(prior_settings$log_params_max)=param_names} #Set up list of invariant parameter values to supply to other functions - consts=list(type=type,log_params_min=log_params_min,log_params_max=log_params_max,mode_start=mode_start, - prior_type=prior_type,norm_prior_sd=norm_prior_sd, - dt=dt,n_reps=n_reps,enviro_data=enviro_data,R0_fixed_values=R0_fixed_values, - vaccine_efficacy=vaccine_efficacy,p_severe_inf=p_severe_inf,p_death_severe_inf=p_death_severe_inf, - p_rep_severe=p_rep_severe,p_rep_death=p_rep_death,m_FOI_Brazil=m_FOI_Brazil,deterministic=deterministic, - mode_parallel=mode_parallel,cluster=cluster) + consts=list(type=type,mode_start=mode_start,prior_settings=prior_settings,dt=dt,n_reps=n_reps,enviro_data=enviro_data, + R0_fixed_values=R0_fixed_values,p_severe_inf=p_severe_inf,p_death_severe_inf=p_death_severe_inf, + add_values=add_values,extra_estimated_params=extra_estimated_params, + deterministic=deterministic,mode_parallel=mode_parallel,cluster=cluster) #MCMC setup chain=chain_prop=posterior_current=posterior_prop=flag_accept=chain_cov_all=NULL @@ -127,14 +121,11 @@ MCMC <- function(log_params_ini=c(),input_data=list(),obs_sero_data=NULL,obs_cas #Iterative estimation for (iter in 1:Niter){ - #cat("\n\tGenerating new parameter values") #Propose new parameter values log_params_prop=param_prop_setup(log_params,chain_cov,adapt) - #cat("\n\tCalculating likelihood") #Calculate likelihood using single_posterior_calc function posterior_value_prop=single_posterior_calc(log_params_prop,input_data,obs_sero_data,obs_case_data,consts) - #cat("\n\tLikelihood calculated") if(is.finite(posterior_value_prop)==FALSE) { p_accept = -Inf @@ -182,7 +173,7 @@ MCMC <- function(log_params_ini=c(),input_data=list(),obs_sero_data=NULL,obs_cas } if(file.exists(filename)==FALSE){file.create(filename)} lines=min((fileIndex * 10000+1),iter):iter - #cat("\nIteration ",iter,sep="") + data_out<-cbind(posterior_current,posterior_prop,exp(chain),flag_accept,exp(chain_prop),chain_cov_all)[lines,] if(fileAccess(filename,2)==0){write.csv(data_out,filename,row.names=FALSE)} } @@ -219,91 +210,67 @@ MCMC <- function(log_params_ini=c(),input_data=list(),obs_sero_data=NULL,obs_cas #' positives #' @param obs_case_data Annual reported case/death data for comparison, by region and year, in format no. cases/no. #' deaths -#' @param consts = List of constant parameters/flags/etc. loaded to mcmc() (type,log_params_min,log_params_max, -#' mode_start,prior_type,norm_prior_sd,dt,n_reps,enviro_data,R0_fixed_values,vaccine_efficacy,p_severe_inf, -#' p_death_severe_inf,p_rep_severe,p_rep_death,m_FOI_Brazil,deterministic, mode_parallel, cluster) +#' @param consts = List of constant parameters/flags/etc. loaded to mcmc() (type, +#' mode_start,prior_settings,dt,n_reps,enviro_data,R0_fixed_values,p_severe_inf, +#' p_death_severe_inf,add_values list,extra_estimated_params,deterministic, mode_parallel, cluster) #' #' @export #' single_posterior_calc <- function(log_params_prop=c(),input_data=list(),obs_sero_data=NULL,obs_case_data=NULL, - consts=list()) { + consts=list()){ regions=input_data$region_labels n_regions=length(regions) - #Get vaccine efficacy and calculate associated prior - if(is.numeric(consts$vaccine_efficacy)==FALSE){ - vaccine_efficacy=exp(log_params_prop[names(log_params_prop)=="vaccine_efficacy"]) - prior_vacc=log(dtrunc(vaccine_efficacy,"norm",a=0,b=1,mean=0.975,sd=0.05)) - } else { - vaccine_efficacy=consts$vaccine_efficacy - prior_vacc=0 - } - - #Get reporting probabilities and check they are within bounds - prior_report=0 - if(is.numeric(consts$p_rep_severe)==FALSE){ - p_rep_severe=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_severe"])) - if(consts$prior_type=="flat"){ - if(p_rep_severeexp(consts$log_params_max[names(consts$log_params_max)=="p_rep_severe"])){prior_report=-Inf} + #Get additional values and calculate associated priors + vaccine_efficacy=p_rep_severe=p_rep_death=m_FOI_Brazil=1.0 + prior_add=0 + for(var_name in names(consts$add_values)){ + if(is.numeric(consts$add_values[[var_name]])==FALSE){ + i=match(var_name,names(log_params_prop)) + value=exp(as.numeric(log_params_prop[i])) + assign(var_name,value) + if(consts$prior_settings$type=="norm"){ + prior_add=prior_add+log(dtrunc(value,"norm",a=0,b=1,mean=consts$prior_settings$norm_params_mean[i], + sd=consts$prior_settings$norm_params_sd[i])) + # if(var_name=="vaccine_efficacy"){ + # prior_add=prior_add+log(dtrunc(value,"norm",a=0,b=1,mean=consts$prior_settings$norm_params_mean[i], + # sd=consts$prior_settings$norm_params_sd[i])) + # } else {prior_add=prior_add+dnorm(log(value),mean = 0,sd = 30,log = TRUE)} + } else { + if(consts$prior_settings$type=="flat"){ + if(valueconsts$prior_settings$log_params_max[i]){prior_add=-Inf} + } + } } else { - if(p_rep_severe<=0.0 || p_rep_severe>1.0){prior_report=-Inf} + assign(var_name,consts$add_values[[var_name]]) } - } else { - p_rep_severe=consts$p_rep_severe } - if(is.numeric(consts$p_rep_death)==FALSE){ - p_rep_death=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_death"])) - if(consts$prior_type=="flat"){ - if(p_rep_deathexp(consts$log_params_max[names(consts$log_params_max)=="p_rep_death"])){prior_report=-Inf} - } else { - if(p_rep_death<=0.0 || p_rep_death>1.0){prior_report=-Inf} + # prior_add=log(dtrunc(vaccine_efficacy,"norm",a=0,b=1,mean=0.975,sd=0.05))+sum(dnorm(log(c(p_rep_severe,p_rep_death)), + # mean = 0,sd = 30, + # log = TRUE)) + + #If additional values give finite prior, get FOI and R0 values and calculate associated prior + if(is.finite(prior_add)){ + FOI_R0_data=mcmc_FOI_R0_setup(consts$type,consts$prior_settings,regions,log_params_prop, + consts$enviro_data,consts$R0_fixed_values) + FOI_values=FOI_R0_data$FOI_values + for(n_region in 1:n_regions){ + if(substr(regions[n_region],1,3)=="BRA"){FOI_values[n_region]=FOI_values[n_region]*m_FOI_Brazil} } - } else { - p_rep_death=consts$p_rep_death - } + R0_values=FOI_R0_data$R0_values + prior_prop=FOI_R0_data$prior+prior_add + }else{prior_prop=-Inf} - #Get value of Brazil multiplier and check bounds - if(is.numeric(consts$m_FOI_Brazil)==FALSE){ - m_FOI_Brazil=as.numeric(exp(log_params_prop[names(log_params_prop)=="m_FOI_Brazil"])) - if(consts$prior_type=="flat"){ - if(m_FOI_Brazilexp(consts$log_params_max[names(consts$log_params_max)=="m_FOI_Brazil"])){prior_report=-Inf} - } else { - if(m_FOI_Brazil<=0.0 || m_FOI_Brazil>1.0){prior_report=-Inf} - } - } else { - m_FOI_Brazil=consts$m_FOI_Brazil - } - - #Get FOI and R0 values and calculate associated prior - FOI_R0_data=mcmc_FOI_R0_setup(consts$type,consts$prior_type,consts$norm_prior_sd,regions,log_params_prop, - consts$enviro_data,consts$R0_fixed_values, - consts$log_params_min,consts$log_params_max) - FOI_values=FOI_R0_data$FOI_values - for(n_region in 1:n_regions){ - if(substr(regions[n_region],1,3)=="BRA"){FOI_values[n_region]=FOI_values[n_region]*m_FOI_Brazil} - } - R0_values=FOI_R0_data$R0_values - prior_prop=prior_vacc+prior_report+FOI_R0_data$prior+sum(dnorm(log(c(p_rep_severe,p_rep_death)), - mean = 0,sd = consts$norm_prior_sd[3],log = TRUE)) - - if(is.null(obs_sero_data)){sero_like_values=0} - if(is.null(obs_case_data)){cases_like_values=deaths_like_values=0} - #cat("\n",prior_prop,sep="") - - #cat("\n\t\tGenerating dataset for regions: ") ### If prior finite, evaluate likelihood ### if (is.finite(prior_prop)) { + if(is.null(obs_sero_data)){sero_like_values=0} + if(is.null(obs_case_data)){cases_like_values=deaths_like_values=0} - #cat("\n\tGenerating dataset") #Generate modelled data over all regions dataset <- Generate_Dataset(input_data,FOI_values,R0_values,obs_sero_data,obs_case_data,vaccine_efficacy,consts$p_severe_inf, consts$p_death_severe_inf,p_rep_severe,p_rep_death,consts$mode_start,start_SEIRV=NULL,consts$dt, consts$n_reps,consts$deterministic,consts$mode_parallel,consts$cluster) - #cat("\n\tDataset generation completed") #Likelihood of observing serological data if(is.null(obs_sero_data)==FALSE){ @@ -319,8 +286,6 @@ single_posterior_calc <- function(log_params_prop=c(),input_data=list(),obs_sero posterior=prior_prop+mean(c(sum(sero_like_values,na.rm=TRUE),sum(cases_like_values,na.rm=TRUE), sum(deaths_like_values,na.rm=TRUE)),na.rm=TRUE) - # dataset<-NULL - # gc() } else {posterior=-Inf} @@ -339,75 +304,46 @@ single_posterior_calc <- function(log_params_prop=c(),input_data=list(),obs_sero #' @param n_regions Number of regions #' @param type Type of parameter set (FOI only, FOI+R0, FOI and/or R0 coefficients associated with environmental #' covariates); choose from "FOI","FOI+R0","FOI enviro","FOI+R0 enviro" -#' @param log_params_min Lower limits of varied parameter values if specified -#' @param log_params_max Upper limits of varied parameter values if specified -#' @param prior_type Text indicating which type of calculation to use for prior probability -#' If prior_type = "zero", prior probability is always zero -#' If prior_type = "flat", prior probability is zero if FOI/R0 in designated ranges, -Inf otherwise -#' If prior_type = "exp", prior probability is given by dexp calculation on FOI/R0 values -#' If prior_type = "norm", prior probability is given by dnorm calculation on parameter values -#' @param norm_prior_sd Vector of 3 values of standard deviations to use in dnorm calculations -#' Value 1 = SD to use for FOI values/coefficients if prior_type = "norm" -#' Value 2 = SD to use for R0 values/coefficients if prior_type = "norm" (if relevant) -#' Value 3 = SD to use for reporting probability prior calculations for all prior_type values (if relevant) +#' @param prior_settings TBA #' @param enviro_data Values of environmental covariates (if in use) #' @param R0_fixed_values Values of R0 to use if not being varied (set to NULL if R0 is being varied) -#' @param vaccine_efficacy Vaccine efficacy (set to NULL if being varied as a parameter) -#' @param p_rep_severe Probability of observation of severe infection (set to NULL if being varied as a parameter) -#' @param p_rep_death Probability of observation of death (set to NULL if being varied as a parameter) -#' @param m_FOI_Brazil FOI multiplier for Brazil regions +#' @param add_values TBA +#' @param extra_estimated_params TBA #' #' @export #' -mcmc_checks <- function(log_params_ini=c(),n_regions=1,type=NULL,log_params_min=c(),log_params_max=c(), - prior_type=NULL,norm_prior_sd=c(30,10,30),enviro_data=NULL,R0_fixed_values=NULL, - vaccine_efficacy=NULL,p_rep_severe=NULL,p_rep_death=NULL,m_FOI_Brazil=1.0){ +mcmc_checks <- function(log_params_ini=c(),n_regions=1,type=NULL,prior_settings=list(type="zero"),enviro_data=NULL, + R0_fixed_values=NULL, + add_values=list(vaccine_efficacy=1.0,p_rep_severe=1.0,p_rep_death=1.0,m_FOI_Brazil=1.0), + extra_estimated_params=list()){ param_names=names(log_params_ini) n_params=length(log_params_ini) assert_that(is.null(param_names)==FALSE) #Check that parameters have been named (should always be done in MCMC()) assert_that(type %in% c("FOI+R0","FOI","FOI+R0 enviro","FOI enviro")) #Check that type is one of those allowed - assert_that(prior_type %in% c("zero","flat","exp","norm")) #Check that prior type is one of those allowed - if(prior_type=="flat"){ - assert_that(length(log_params_min)==n_params) - assert_that(length(log_params_max)==n_params) + assert_that(prior_settings$type %in% c("zero","flat","exp","norm")) #Check that prior type is one of those allowed + if(prior_settings$type=="flat"){ + assert_that(length(prior_settings$log_params_min)==n_params) + assert_that(length(prior_settings$log_params_max)==n_params) } - assert_that(length(norm_prior_sd)==3) - assert_that(all(norm_prior_sd>0)) - - #If vaccine efficacy, severe case reporting probability and/or fatal case reporting probability have been set to - #NULL, check that they appear among the parameters (should always have been set up in MCMC() if - #create_param_labels() ran correctly ) - if(is.numeric(vaccine_efficacy)==TRUE){ - flag_vacc_eff=0 - assert_that(0.0<=vaccine_efficacy && vaccine_efficacy<=1.0) - } else { - flag_vacc_eff=1 - assert_that("vaccine_efficacy" %in% param_names) + if(prior_settings$type=="norm"){ + assert_that(length(prior_settings$norm_params_mean)==n_params) + assert_that(length(prior_settings$norm_params_sd)==n_params) } - if(is.numeric(p_rep_severe)==TRUE){ - flag_severe=0 - assert_that(0.0<=p_rep_severe && p_rep_severe<=1.0) - } else { - flag_severe=1 - assert_that("p_rep_severe" %in% param_names) - } - if(is.numeric(p_rep_death)==TRUE){ - flag_death=0 - assert_that(0.0<=p_rep_death && p_rep_death<=1.0) - } else { - flag_death=1 - assert_that("p_rep_death" %in% param_names) - } - if(is.numeric(m_FOI_Brazil)==TRUE){ - flag_br=0 - assert_that(0.0<=m_FOI_Brazil && m_FOI_Brazil<=1.0) - } else { - flag_br=1 - assert_that("m_FOI_Brazil" %in% param_names) + + # Check additional values + add_value_names=names(add_values) + assert_that(add_value_names==c("vaccine_efficacy","p_rep_severe","p_rep_death","m_FOI_Brazil")) #TBV for alt version + assert_that(all(extra_estimated_params %in% add_value_names)) + for(var_name in add_value_names){ + if(var_name %in% extra_estimated_params){ + assert_that(is.null(add_values[[var_name]])) + }else{ + assert_that(add_values[[var_name]]<=1.0 && add_values[[var_name]]>=0.0) + } } - #If environmental data has been supplied, get names of variables + # If environmental data has been supplied, get names of variables if(is.null(enviro_data)==FALSE){ env_vars=names(enviro_data[c(2:ncol(enviro_data))]) n_env_vars=length(env_vars) @@ -415,51 +351,65 @@ mcmc_checks <- function(log_params_ini=c(),n_regions=1,type=NULL,log_params_min= assert_that(type %in% c("FOI+R0","FOI")) } - #Check that total number of parameters is correct based on "type" and on number of additional parameters (vaccine - #efficacy, reporting probabilities); check parameters named in correct order + # Check that total number of parameters is correct based on "type" and on number of additional parameters (vaccine + # efficacy, reporting probabilities); check parameters named in correct order (TBA); all should be correct if parameter + # names created using create_param_labels if(type=="FOI+R0"){ - assert_that(n_params==(2*n_regions)+flag_vacc_eff+flag_severe+flag_death+flag_br) - if(flag_vacc_eff==1){assert_that(param_names[1+(2*n_regions)]=="vaccine_efficacy")} - if(flag_severe==1){assert_that(param_names[1+flag_vacc_eff+(2*n_regions)]=="p_rep_severe")} - if(flag_death==1){assert_that(param_names[1+flag_vacc_eff+flag_severe+(2*n_regions)]=="p_rep_death")} - if(flag_br==1){assert_that(param_names[1+flag_vacc_eff+flag_severe+flag_death+(2*n_regions)]=="m_FOI_Brazil")} + assert_that(n_params==(2*n_regions)+length(extra_estimated_params)) } if(type=="FOI"){ + assert_that(n_params==n_regions+length(extra_estimated_params)) assert_that(is.null(R0_fixed_values)==FALSE) assert_that(length(R0_fixed_values)==n_regions) - assert_that(n_params==n_regions+flag_vacc_eff+flag_severe+flag_death+flag_br) - if(flag_vacc_eff==1){assert_that(param_names[1+n_regions]=="vaccine_efficacy")} - if(flag_severe==1){assert_that(param_names[1+flag_vacc_eff+n_regions]=="p_rep_severe")} - if(flag_death==1){assert_that(param_names[1+flag_vacc_eff+flag_severe+n_regions]=="p_rep_death")} - if(flag_br==1){assert_that(param_names[1+flag_vacc_eff+flag_severe+flag_death+n_regions]=="m_FOI_Brazil")} } if(type=="FOI+R0 enviro"){ assert_that(is.null(enviro_data)==FALSE) - assert_that(n_params==(2*n_env_vars)+flag_vacc_eff+flag_severe+flag_death+flag_br) + assert_that(n_params==(2*n_env_vars)+length(extra_estimated_params)) for(i in 1:n_env_vars){ assert_that(param_names[i]==paste("FOI_",env_vars[i],sep="")) assert_that(param_names[i+n_env_vars]==paste("R0_",env_vars[i],sep="")) } - if(flag_vacc_eff==1){assert_that(param_names[1+(2*n_env_vars)]=="vaccine_efficacy")} - if(flag_severe==1){assert_that(param_names[1+flag_vacc_eff+(2*n_env_vars)]=="p_rep_severe")} - if(flag_death==1){assert_that(param_names[1+flag_vacc_eff+flag_severe+(2*n_env_vars)]=="p_rep_death")} - if(flag_br==1){assert_that(param_names[1+flag_vacc_eff+flag_severe+flag_death+(2*n_env_vars)]=="m_FOI_Brazil")} } if(type=="FOI enviro"){ assert_that(is.null(enviro_data)==FALSE) + assert_that(n_params==n_env_vars+length(extra_estimated_params)) assert_that(is.null(R0_fixed_values)==FALSE) assert_that(length(R0_fixed_values)==n_regions) - assert_that(n_params==n_env_vars+flag_vacc_eff+flag_severe+flag_death+flag_br) - for(i in 1:n_env_vars){assert_that(param_names[i]==paste("FOI_",env_vars[i],sep=""))} - if(flag_vacc_eff==1){assert_that(param_names[1+n_env_vars]=="vaccine_efficacy")} - if(flag_severe==1){assert_that(param_names[1+flag_vacc_eff+n_env_vars]=="p_rep_severe")} - if(flag_death==1){assert_that(param_names[1+flag_vacc_eff+flag_severe+n_env_vars]=="p_rep_death")} - if(flag_br==1){assert_that(param_names[1+flag_vacc_eff+flag_severe+flag_death+n_env_vars]=="m_FOI_Brazil")} } return(NULL) } #------------------------------------------------------------------------------- +#' @title param_prop_setup +#' +#' @description Set up proposed new log parameter values for next step in chain +#' +#' @details Takes in current values of parameter set used for Markov Chain Monte Carlo fitting and proposes new values +#' from multivariate normal distribution where the existing values form the mean and the standard deviation is +#' based on the chain covariance or (if the flag "adapt" is set to 1) a flat value based on the number of parameters. +#' +#' @param log_params Previous log parameter values used as input +#' @param chain_cov Covariance calculated from previous steps in chain +#' @param adapt 0/1 flag indicating which type of calculation to use for proposition value +#' ' +#' @export +#' +param_prop_setup <- function(log_params=c(),chain_cov=1,adapt=0){ + + n_params = length(log_params) + if (adapt==1) { + sigma = (2.38 ^ 2) * chain_cov / n_params #'optimal' scaling of chain covariance + log_params_prop_a = rmvnorm(n = 1, mean = log_params, sigma = sigma) + } else { + sigma = ((1e-2) ^ 2) * diag(n_params) / n_params #this is an inital proposal covariance, see [Mckinley et al 2014] + log_params_prop_a = rmvnorm(n = 1, mean = log_params, sigma = sigma) + } + log_params_prop = log_params_prop_a[1,] + #names(log_params_prop)=names(log_params) + + return(log_params_prop) +} +#------------------------------------------------------------------------------- #' @title mcmc_FOI_R0_setup #' #' @description Set up FOI and R0 values and calculate some prior probability values for MCMC calculation @@ -470,23 +420,16 @@ mcmc_checks <- function(log_params_ini=c(),n_regions=1,type=NULL,log_params_min= #' #' @param type Type of parameter set (FOI only, FOI+R0, FOI and/or R0 coefficients associated with environmental #' covariates); choose from "FOI","FOI+R0","FOI enviro","FOI+R0 enviro" -#' @param prior_type Text indicating which type of calculation to use for prior probability -#' If prior_type = "zero", prior probability is always zero -#' If prior_type = "flat", prior probability is zero if FOI/R0 in designated ranges, -Inf otherwise -#' If prior_type = "exp", prior probability is given by dexp calculation on FOI/R0 values -#' If prior_type = "norm", prior probability is given by dnorm calculation on parameter values -#' @param norm_prior_sd TBA +#' @param prior_settings TBA #' @param regions Vector of region names #' @param log_params_prop Proposed values of varied parameters (natural logarithm of actual parameters) #' @param enviro_data Environmental data frame, containing only relevant environmental variables #' @param R0_fixed_values Values of R0 to use if not being varied -#' @param log_params_min Lower limits of varied parameter values if specified -#' @param log_params_max Upper limits of varied parameter values if specified #' ' #' @export #' -mcmc_FOI_R0_setup <- function(type="",prior_type="",norm_prior_sd=c(30,10,30),regions="",log_params_prop=c(),enviro_data=list(), - R0_fixed_values=c(),log_params_min=c(),log_params_max=c()){ +mcmc_FOI_R0_setup <- function(type="",prior_settings=list(type="zero"),regions="",log_params_prop=c(),enviro_data=list(), + R0_fixed_values=c()){ n_params=length(log_params_prop) n_regions=length(regions) @@ -494,55 +437,42 @@ mcmc_FOI_R0_setup <- function(type="",prior_type="",norm_prior_sd=c(30,10,30),re if(type %in% c("FOI+R0 enviro","FOI enviro")){ n_env_vars=ncol(enviro_data)-1 - if(type=="FOI+R0 enviro"){ - enviro_coeffs=exp(log_params_prop[c(1:(2*n_env_vars))]) - } else { - enviro_coeffs=exp(log_params_prop[c(1:n_env_vars)])} + if(type=="FOI+R0 enviro"){n_values=2*n_env_vars}else{n_values=n_env_vars} + enviro_coeffs=exp(log_params_prop[c(1:n_values)]) + for(i in 1:n_regions){ model_params=param_calc_enviro(enviro_coeffs, as.numeric(enviro_data[enviro_data$region==regions[i],1+c(1:n_env_vars)])) FOI_values[i]=model_params$FOI if(type=="FOI+R0 enviro"){R0_values[i]=model_params$R0} else {R0_values[i]=R0_fixed_values[i]} } - } - if(type %in% c("FOI+R0","FOI")){ + } else { FOI_values=exp(log_params_prop[c(1:n_regions)]) - if(type=="FOI+R0"){R0_values=exp(log_params_prop[c((n_regions+1):(2*n_regions))]) - } else {R0_values=R0_fixed_values} + if(type=="FOI+R0"){ + n_values=2*n_regions + R0_values=exp(log_params_prop[c((n_regions+1):n_values)]) + } else { + n_values=n_regions + R0_values=R0_fixed_values + } } prior=0 - if(prior_type=="exp"){ - prior_FOI=dexp(FOI_values,rate=1,log=TRUE) - if(type %in% c("FOI+R0","FOI+R0 enviro")){prior_R0=dexp(R0_values,rate=1,log=TRUE)} else {prior_R0=0} - prior = prior+sum(prior_FOI)+sum(prior_R0) - } - if(prior_type=="flat"){ - if(is.null(log_params_min)==FALSE){ - for(i in 1:n_params){ - if(log_params_prop[i]log_params_max[i]){prior=-Inf} + if(log_params_prop[i]prior_settings$log_params_max[i]){prior=-Inf} } - } - } - if(prior_type=="norm"){ - if(type=="FOI"){ - prior=sum(dnorm(log_params_prop[c(1:n_regions)],mean = 0,sd = norm_prior_sd[1],log = TRUE)) - } - if(type=="FOI+R0"){ - prior=sum(dnorm(log_params_prop[c(1:n_regions)],mean = 0,sd = norm_prior_sd[1],log = TRUE)) + - sum(dnorm(log_params_prop[c(1:n_regions)+n_regions],mean = 0,sd = norm_prior_sd[2],log = TRUE)) - } - if(type=="FOI enviro"){ - prior=sum(dnorm(log_params_prop[c(1:n_env_vars)],mean = 0,sd = norm_prior_sd[1],log = TRUE)) + } else { + if(prior_settings$type=="exp"){ + prior_FOI=dexp(FOI_values,rate=1,log=TRUE) + if(type %in% c("FOI+R0","FOI+R0 enviro")){prior_R0=dexp(R0_values,rate=1,log=TRUE)} else {prior_R0=0} + prior = prior+sum(prior_FOI)+sum(prior_R0) } - if(type=="FOI+R0 enviro"){ - prior=sum(dnorm(log_params_prop[c(1:n_env_vars)],mean = 0,sd = norm_prior_sd[1],log = TRUE)) + - sum(dnorm(log_params_prop[c(1:n_env_vars)+n_env_vars],mean = 0,sd = norm_prior_sd[2],log = TRUE)) } } @@ -552,36 +482,6 @@ mcmc_FOI_R0_setup <- function(type="",prior_type="",norm_prior_sd=c(30,10,30),re return(output) } #------------------------------------------------------------------------------- -#' @title param_prop_setup -#' -#' @description Set up proposed new log parameter values for next step in chain -#' -#' @details Takes in current values of parameter set used for Markov Chain Monte Carlo fitting and proposes new values -#' from multivariate normal distribution where the existing values form the mean and the standard deviation is -#' based on the chain covariance or (if the flag "adapt" is set to 1) a flat value based on the number of parameters. -#' -#' @param log_params Previous log parameter values used as input -#' @param chain_cov Covariance calculated from previous steps in chain -#' @param adapt 0/1 flag indicating which type of calculation to use for proposition value -#' ' -#' @export -#' -param_prop_setup <- function(log_params=c(),chain_cov=1,adapt=0){ - - n_params = length(log_params) - if (adapt==1) { - sigma = (2.38 ^ 2) * chain_cov / n_params #'optimal' scaling of chain covariance - log_params_prop_a = rmvnorm(n = 1, mean = log_params, sigma = sigma) - } else { - sigma = ((1e-2) ^ 2) * diag(n_params) / n_params #this is an inital proposal covariance, see [Mckinley et al 2014] - log_params_prop_a = rmvnorm(n = 1, mean = log_params, sigma = sigma) - } - log_params_prop = log_params_prop_a[1,] - #names(log_params_prop)=names(log_params) - - return(log_params_prop) -} -#------------------------------------------------------------------------------- #' @title mcmc_prelim_fit #' #' @description Test multiple sets of parameters randomly drawn from range between maximum and minimum @@ -609,21 +509,18 @@ param_prop_setup <- function(log_params=c(),chain_cov=1,adapt=0){ #' If mode_start=0, only vaccinated individuals #' If mode_start=1, shift some non-vaccinated individuals into recovered to give herd immunity (uniform by age, R0 based only) #' If mode_start=3, shift some non-vaccinated individuals into recovered to give herd immunity (stratified by age) -#' @param prior_type Text indicating which type of calculation to use for prior probability -#' If prior_type = "zero", prior probability is always zero -#' If prior_type = "exp", prior probability is given by dexp calculation on FOI/R0 values -#' If prior_type = "norm", prior probability is given by dnorm calculation on parameter values -#' @param norm_prior_sd TBA +#' @param prior_settings TBA #' @param dt time increment in days (must be 1 or 5) #' @param n_reps Number of repetitions #' @param enviro_data Values of environmental variables (if in use) #' @param R0_fixed_values Values of R0 to use if not being varied -#' @param vaccine_efficacy Vaccine efficacy (set to NULL if being varied as a parameter) #' @param p_severe_inf TBA #' @param p_death_severe_inf TBA -#' @param p_rep_severe Probability of observation of severe infection (set to NULL if being varied as a parameter) -#' @param p_rep_death Probability of observation of death (set to NULL if being varied as a parameter) -#' @param m_FOI_Brazil Multiplier of FOI in Brazil regions +#' @param add_values List containing values of additional inputs or NULL to indicate they are being varied as parameters: +#' vaccine_efficacy Vaccine efficacy (set to NULL if being varied as a parameter) +#' p_rep_severe Probability of observation of severe infection (set to NULL if being varied as a parameter) +#' p_rep_death Probability of observation of death (set to NULL if being varied as a parameter) +#' m_FOI_Brazil Multiplier of spillover FOI for Brazil regions (set to NULL if being varied as a parameter) #' @param deterministic TBA #' @param mode_parallel TBA #' @param cluster Cluster of threads to use if multithreading to be used; set to NULL otherwise @@ -632,25 +529,28 @@ param_prop_setup <- function(log_params=c(),chain_cov=1,adapt=0){ #' mcmc_prelim_fit <- function(n_iterations=1,n_param_sets=1,n_bounds=1,type=NULL,log_params_min=NULL, log_params_max=NULL,input_data=list(),obs_sero_data=list(),obs_case_data=list(), - mode_start=0,prior_type="zero",norm_prior_sd=c(30,10,30), - dt=1.0,n_reps=1,enviro_data=NULL,R0_fixed_values=c(),vaccine_efficacy=NULL, - p_severe_inf = 0.12, p_death_severe_inf = 0.39,p_rep_severe=NULL,p_rep_death=NULL,m_FOI_Brazil=1.0, + mode_start=0,prior_settings=list(type="zero"),dt=1.0,n_reps=1,enviro_data=NULL,R0_fixed_values=c(), + p_severe_inf = 0.12, p_death_severe_inf = 0.39, + add_values=list(vaccine_efficacy=1.0,p_rep_severe=1.0,p_rep_death=1.0,m_FOI_Brazil=1.0), deterministic=TRUE,mode_parallel="none",cluster=NULL){ #TODO - Add assertthat functions assert_that(mode_start %in% c(0,1,3),msg="mode_start must have value 0, 1 or 3") assert_that(length(log_params_min)==length(log_params_max),msg="Parameter limit vectors must have same lengths") assert_that(type %in% c("FOI+R0","FOI","FOI+R0 enviro","FOI enviro")) - assert_that(prior_type %in% c("zero","exp","norm")) + assert_that(prior_settings$type %in% c("zero","exp","norm")) best_fit_results=list() n_params=length(log_params_min) - extra_params=c() - if(is.null(vaccine_efficacy)==TRUE){extra_params=append(extra_params,"vaccine_efficacy")} - if(is.null(p_rep_severe)==TRUE){extra_params=append(extra_params,"p_rep_severe")} - if(is.null(p_rep_death)==TRUE){extra_params=append(extra_params,"p_rep_death")} - if(is.null(m_FOI_Brazil)==TRUE){extra_params=append(extra_params,"m_FOI_Brazil")} - param_names=create_param_labels(type,input_data,enviro_data,extra_params) + + #Get additional values + extra_estimated_params=c() + add_value_names=names(add_values) + assert_that(add_value_names==c("vaccine_efficacy","p_rep_severe","p_rep_death","m_FOI_Brazil")) #TBV in alt version + for(var_name in add_value_names){ + if(is.null(add_values[[var_name]])==TRUE){extra_estimated_params=append(extra_estimated_params,var_name)} + } + param_names=create_param_labels(type,input_data,enviro_data,extra_estimated_params) #TODO - Additional assert_that checks assert_that(length(param_names)==n_params) @@ -665,24 +565,17 @@ mcmc_prelim_fit <- function(n_iterations=1,n_param_sets=1,n_bounds=1,type=NULL,l cat("\nIteration: ",iteration,"\n",sep="") all_param_sets <- lhs(n=n_param_sets,rect=cbind(log_params_min,log_params_max)) results=data.frame() - consts=list(type=type,log_params_min=log_params_min,log_params_max=log_params_max, - mode_start=mode_start,prior_type=prior_type,norm_prior_sd=norm_prior_sd, - dt=dt,n_reps=n_reps,enviro_data=enviro_data, - R0_fixed_values=R0_fixed_values,vaccine_efficacy=vaccine_efficacy, - p_severe_inf = p_severe_inf, p_death_severe_inf = p_death_severe_inf, - p_rep_severe=p_rep_severe,p_rep_death=p_rep_death,m_FOI_Brazil=m_FOI_Brazil, + consts=list(type=type,mode_start=mode_start,prior_settings=prior_settings, + dt=dt,n_reps=n_reps,enviro_data=enviro_data,R0_fixed_values=R0_fixed_values, + p_severe_inf = p_severe_inf, p_death_severe_inf = p_death_severe_inf,add_values=add_values, deterministic=deterministic,mode_parallel=mode_parallel,cluster=cluster) for(set in 1:n_param_sets){ - #if(set %% 10 == 0){cat("\n")} - #cat(set,"\t",sep="") cat("\n\tSet: ",set,sep="") log_params_prop=all_param_sets[set,] cat("\n\tParams: ",signif(log_params_prop,3)) - #cat(log_params_prop,file=progress_file,sep="\t",append=TRUE) - names(log_params_prop)=param_names posterior_value=single_posterior_calc(log_params_prop,input_data,obs_sero_data,obs_case_data,consts) results<-rbind(results,c(set,exp(log_params_prop),posterior_value)) diff --git a/R/mcmc2.R b/R/mcmc2.R index 86350b6..1e4e116 100644 --- a/R/mcmc2.R +++ b/R/mcmc2.R @@ -1,615 +1,615 @@ -# R file for functions used for Markov Chain Monte Carlo fitting (and preliminary maximum-likelihood fitting) in -# YEP package - alternate versions -#------------------------------------------------------------------------------- -#------------------------------------------------------------------------------- -#' @title MCMC2 -#' -#' @description TBA -#' -#' @details TBA -#' -#' @param log_params_ini TBA -#' @param input_data TBA -#' @param obs_sero_data TBA -#' @param obs_case_data TBA -#' @param filename_prefix TBA -#' @param Niter TBA -#' @param type TBA -#' @param log_params_min TBA -#' @param log_params_max TBA -#' @param mode_start TBA -#' @param prior_type TBA -#' @param norm_prior_sd TBA -#' @param dt TBA -#' @param n_reps TBA -#' @param enviro_data TBA -#' @param R0_fixed_values TBA -#' @param vaccine_efficacy TBA -#' @param p_severe_inf TBA -#' @param p_death_severe_inf TBA -#' @param p_rep_severe_af TBA -#' @param p_rep_death_af TBA -#' @param p_rep_severe_sa TBA -#' @param p_rep_death_sa TBA -#' @param m_FOI_Brazil TBA -#' @param deterministic TBA -#' @param mode_parallel TBA -#' @param cluster TBA -#' ' -#' @export -#' -MCMC2 <- function(log_params_ini=c(),input_data=list(),obs_sero_data=NULL,obs_case_data=NULL,filename_prefix="Chain", - Niter=1,type=NULL,log_params_min=c(),log_params_max=c(),mode_start=0,prior_type="zero", - norm_prior_sd=c(30,10,30),dt=1.0,n_reps=1,enviro_data=NULL,R0_fixed_values=NULL, - vaccine_efficacy=NULL,p_severe_inf=0.12,p_death_severe_inf=0.39, - p_rep_severe_af=NULL,p_rep_death_af=NULL,p_rep_severe_sa=NULL,p_rep_death_sa=NULL, - m_FOI_Brazil=1.0,deterministic=FALSE,mode_parallel="none",cluster=NULL){ - - assert_that(is.logical(deterministic)) - assert_that(mode_start %in% c(0,1,3),msg="mode_start must have value 0, 1 or 3") - n_params=length(log_params_ini) - - extra_params=c() - if(is.null(vaccine_efficacy)==TRUE){extra_params=append(extra_params,"vaccine_efficacy")} - if(is.null(p_rep_severe_af)==TRUE){extra_params=append(extra_params,"p_rep_severe_af")} - if(is.null(p_rep_death_af)==TRUE){extra_params=append(extra_params,"p_rep_death_af")} - if(is.null(p_rep_severe_sa)==TRUE){extra_params=append(extra_params,"p_rep_severe_sa")} - if(is.null(p_rep_death_sa)==TRUE){extra_params=append(extra_params,"p_rep_death_sa")} - if(is.null(m_FOI_Brazil)==TRUE){extra_params=append(extra_params,"m_FOI_Brazil")} - - #Process input data to check that all regions with sero and/or case data supplied are present, remove - #regions without any supplied data, and add cross-referencing tables for use when calculating likelihood. Take - #subset of environmental data (if used) and check that environmental data available for all regions - input_data=input_data_process(input_data,obs_sero_data,obs_case_data) - regions=names(table(input_data$region_labels)) #Regions in new processed input data list - n_regions=length(regions) - if(is.null(enviro_data)==FALSE){ - for(region in regions){assert_that(region %in% enviro_data$region)} - enviro_data=subset(enviro_data,enviro_data$region %in% regions) - } - - #Label parameters according to order and fitting type - param_names=create_param_labels(type,input_data,enviro_data,extra_params) - names(log_params_ini)=param_names - - #Set up list of invariant parameter values to supply to other functions - consts=list(type=type,log_params_min=log_params_min,log_params_max=log_params_max,mode_start=mode_start, - prior_type=prior_type,norm_prior_sd=norm_prior_sd, - dt=dt,n_reps=n_reps,enviro_data=enviro_data,R0_fixed_values=R0_fixed_values, - vaccine_efficacy=vaccine_efficacy,p_severe_inf=p_severe_inf,p_death_severe_inf=p_death_severe_inf, - p_rep_severe_af=p_rep_severe_af,p_rep_death_af=p_rep_death_af, - p_rep_severe_sa=p_rep_severe_sa,p_rep_death_sa=p_rep_death_sa, - m_FOI_Brazil=m_FOI_Brazil,deterministic=deterministic, - mode_parallel=mode_parallel,cluster=cluster) - - #MCMC setup - chain=chain_prop=posterior_current=posterior_prop=flag_accept=chain_cov_all=NULL - burnin = min(2*n_params, Niter) - fileIndex = 0 - log_params=log_params_ini - chain_cov=1 - adapt=0 - posterior_value_current=-Inf - - #Iterative estimation - for (iter in 1:Niter){ - - #cat("\n\tGenerating new parameter values") - #Propose new parameter values - log_params_prop=param_prop_setup(log_params,chain_cov,adapt) - - #cat("\n\tCalculating likelihood") - #Calculate likelihood using single_posterior_calc2 function - posterior_value_prop=single_posterior_calc2(log_params_prop,input_data,obs_sero_data,obs_case_data,consts) - #cat("\n\tLikelihood calculated") - - if(is.finite(posterior_value_prop)==FALSE) { - p_accept = -Inf - } else { - p_accept = posterior_value_prop - posterior_value_current - if(is.na(p_accept) ){ p_accept = -Inf} - } - - ## accept/reject step: - tmp = runif(1) - if(tmp=10){ - filename=paste0(filename_prefix,fileIndex,".csv") - } else { - filename=paste0(filename_prefix,"0",fileIndex,".csv") - } - if(file.exists(filename)==FALSE){file.create(filename)} - lines=min((fileIndex * 10000+1),iter):iter - #cat("\nIteration ",iter,sep="") - data_out<-cbind(posterior_current,posterior_prop,exp(chain),flag_accept,exp(chain_prop),chain_cov_all)[lines,] - if(fileAccess(filename,2)==0){write.csv(data_out,filename,row.names=FALSE)} - } - - #Decide whether next iteration will be adaptive - if (iter>burnin & runif(1)<0.9){ #adapt - adapt = 1 - chain_cov = cov(chain[max(nrow(chain)-10000, 1):nrow(chain),]) - } else { - adapt = 0 - chain_cov = 1 - } - } - - #Get final parameter values - param_out=exp(log_params) - names(param_out)=names(log_params_ini) - - return(param_out) -} -#------------------------------------------------------------------------------- -#' @title single_posterior_calc2 -#' -#' @description TBA -#' -#' @details TBA -#' -#' @param log_params_prop TBA -#' @param input_data TBA -#' @param obs_sero_data TBA -#' @param obs_case_data TBA -#' @param consts TBA -#' -#' @export -#' -single_posterior_calc2 <- function(log_params_prop=c(),input_data=list(),obs_sero_data=NULL,obs_case_data=NULL, - consts=list()) { - - regions=input_data$region_labels - n_regions=length(regions) - - #Get vaccine efficacy and calculate associated prior - if(is.numeric(consts$vaccine_efficacy)==FALSE){ - vaccine_efficacy=exp(log_params_prop[names(log_params_prop)=="vaccine_efficacy"]) - prior_vacc=log(dtrunc(vaccine_efficacy,"norm",a=0,b=1,mean=0.975,sd=0.05)) - } else { - vaccine_efficacy=consts$vaccine_efficacy - prior_vacc=0 - } - - #Get reporting probabilities and check they are within specified bounds - prior_report=0 - if(is.numeric(consts$p_rep_severe_af)==FALSE){ - p_rep_severe_af=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_severe_af"])) - if(p_rep_severe_af<=0.0 || p_rep_severe_af>1.0){prior_report=-Inf} - } else { - p_rep_severe_af=consts$p_rep_severe_af - } - if(is.numeric(consts$p_rep_severe_sa)==FALSE){ - p_rep_severe_sa=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_severe_sa"])) - if(p_rep_severe_sa<=0.0 || p_rep_severe_sa>1.0){prior_report=-Inf} - } else { - p_rep_severe_sa=consts$p_rep_severe_sa - } - if(is.numeric(consts$p_rep_death_af)==FALSE){ - p_rep_death_af=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_death_af"])) - if(p_rep_death_af<=0.0 || p_rep_death_af>1.0){prior_report=-Inf} - } else { - p_rep_death_af=consts$p_rep_death_af - } - if(is.numeric(consts$p_rep_death_sa)==FALSE){ - p_rep_death_sa=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_death_sa"])) - if(p_rep_death_sa<=0.0 || p_rep_death_sa>1.0){prior_report=-Inf} - } else { - p_rep_death_sa=consts$p_rep_death_sa - } - - #Get value of Brazil multiplier and check bounds - if(is.numeric(consts$m_FOI_Brazil)==FALSE){ - m_FOI_Brazil=as.numeric(exp(log_params_prop[names(log_params_prop)=="m_FOI_Brazil"])) - if(m_FOI_Brazil<=0.0 || m_FOI_Brazil>1.0){prior_report=-Inf} - } else { - m_FOI_Brazil=consts$m_FOI_Brazil - } - - #Get FOI and R0 values and calculate associated prior - FOI_R0_data=mcmc_FOI_R0_setup(consts$type,consts$prior_type,consts$norm_prior_sd,regions,log_params_prop, - consts$enviro_data,consts$R0_fixed_values, - consts$log_params_min,consts$log_params_max) - FOI_values=FOI_R0_data$FOI_values - for(n_region in 1:n_regions){ - if(substr(regions[n_region],1,3)=="BRA"){FOI_values[n_region]=FOI_values[n_region]*m_FOI_Brazil} - } - R0_values=FOI_R0_data$R0_values - prior_prop=prior_vacc+prior_report+FOI_R0_data$prior + - sum(dnorm(log(c(p_rep_severe_af,p_rep_death_af,p_rep_severe_sa,p_rep_death_sa)), - mean = 0,sd = consts$norm_prior_sd[3],log = TRUE)) - - if(is.null(obs_sero_data)){sero_like_values=0} - if(is.null(obs_case_data)){cases_like_values=deaths_like_values=0} - #cat("\n",prior_prop,sep="") - - #cat("\n\t\tGenerating dataset for regions: ") - ### If prior finite, evaluate likelihood ### - if (is.finite(prior_prop)) { - - #cat("\n\tGenerating dataset") - #Generate modelled data over all regions - dataset <- Generate_Dataset2(input_data,FOI_values,R0_values,obs_sero_data,obs_case_data,vaccine_efficacy,consts$p_severe_inf, - consts$p_death_severe_inf,p_rep_severe_af,p_rep_death_af,p_rep_severe_sa,p_rep_death_sa,consts$mode_start,start_SEIRV=NULL,consts$dt, - consts$n_reps,consts$deterministic,consts$mode_parallel,consts$cluster) - #cat("\n\tDataset generation completed") - - #Likelihood of observing serological data - if(is.null(obs_sero_data)==FALSE){ - sero_like_values=sero_data_compare(dataset$model_sero_values,obs_sero_data) - } - #Likelihood of observing annual case/death data - if(is.null(obs_case_data)==FALSE){ - cases_like_values=case_data_compare(dataset$model_case_values,obs_case_data$cases) - if(is.null(obs_case_data$deaths)==FALSE){ - deaths_like_values=case_data_compare(dataset$model_death_values,obs_case_data$deaths) - } else {deaths_like_values=0} - } - - posterior=prior_prop+mean(c(sum(sero_like_values,na.rm=TRUE),sum(cases_like_values,na.rm=TRUE), - sum(deaths_like_values,na.rm=TRUE)),na.rm=TRUE) - # dataset<-NULL - # gc() - - } else {posterior=-Inf} - - return(posterior) -} -#------------------------------------------------------------------------------- -#' @title mcmc_prelim_fit2 -#' -#' @description TBA -#' -#' @details TBA -#' -#' @param n_iterations TBA -#' @param n_param_sets TBA -#' @param n_bounds TBA -#' @param type TBA -#' @param log_params_min TBA -#' @param log_params_max TBA -#' @param input_data TBA -#' @param obs_sero_data TBA -#' @param obs_case_data TBA -#' @param mode_start TBA -#' @param prior_type TBA -#' @param norm_prior_sd TBA -#' @param dt TBA -#' @param n_reps TBA -#' @param enviro_data TBA -#' @param R0_fixed_values TBA -#' @param vaccine_efficacy TBA -#' @param p_severe_inf TBA -#' @param p_death_severe_inf TBA -#' @param p_rep_severe_af TBA -#' @param p_rep_death_af TBA -#' @param p_rep_severe_sa TBA -#' @param p_rep_death_sa TBA -#' @param m_FOI_Brazil TBA -#' @param deterministic TBA -#' @param mode_parallel TBA -#' @param cluster TBA -#' ' -#' @export -#' -mcmc_prelim_fit2 <- function(n_iterations=1,n_param_sets=1,n_bounds=1,type=NULL,log_params_min=NULL, - log_params_max=NULL,input_data=list(),obs_sero_data=list(),obs_case_data=list(), - mode_start=0,prior_type="zero",norm_prior_sd=c(30,10,30), - dt=1.0,n_reps=1,enviro_data=NULL,R0_fixed_values=c(),vaccine_efficacy=NULL, - p_severe_inf = 0.12, p_death_severe_inf = 0.39, - p_rep_severe_af=NULL,p_rep_death_af=NULL,p_rep_severe_sa=NULL,p_rep_death_sa=NULL, - m_FOI_Brazil=1.0,deterministic=TRUE,mode_parallel="none",cluster=NULL){ - - #TODO - Add assertthat functions - assert_that(mode_start %in% c(0,1,3),msg="mode_start must have value 0, 1 or 3") - assert_that(length(log_params_min)==length(log_params_max),msg="Parameter limit vectors must have same lengths") - assert_that(type %in% c("FOI+R0","FOI","FOI+R0 enviro","FOI enviro")) - assert_that(prior_type %in% c("zero","exp","norm")) - - best_fit_results=list() - n_params=length(log_params_min) - extra_params=c() - if(is.null(vaccine_efficacy)==TRUE){extra_params=append(extra_params,"vaccine_efficacy")} - if(is.null(p_rep_severe_af)==TRUE){extra_params=append(extra_params,"p_rep_severe_af")} - if(is.null(p_rep_death_af)==TRUE){extra_params=append(extra_params,"p_rep_death_af")} - if(is.null(p_rep_severe_sa)==TRUE){extra_params=append(extra_params,"p_rep_severe_sa")} - if(is.null(p_rep_death_sa)==TRUE){extra_params=append(extra_params,"p_rep_death_sa")} - if(is.null(m_FOI_Brazil)==TRUE){extra_params=append(extra_params,"m_FOI_Brazil")} - param_names=create_param_labels(type,input_data,enviro_data,extra_params) - - #TODO - Additional assert_that checks - assert_that(length(param_names)==n_params) - names(log_params_min)=names(log_params_max)=param_names - xlabels=param_names - for(i in 1:n_params){xlabels[i]=substr(xlabels[i],1,15)} - ylabels=10^c(-8,-6,-4,-3,-2,-1,0,1) - par(mar=c(6,2,1,1)) - ylim=c(min(log_params_min),max(log_params_max)) - - for(iteration in 1:n_iterations){ - cat("\nIteration: ",iteration,"\n",sep="") - all_param_sets <- lhs(n=n_param_sets,rect=cbind(log_params_min,log_params_max)) - results=data.frame() - consts=list(type=type,log_params_min=log_params_min,log_params_max=log_params_max, - mode_start=mode_start,prior_type=prior_type,norm_prior_sd=norm_prior_sd, - dt=dt,n_reps=n_reps,enviro_data=enviro_data, - R0_fixed_values=R0_fixed_values,vaccine_efficacy=vaccine_efficacy, - p_severe_inf = p_severe_inf, p_death_severe_inf = p_death_severe_inf, - p_rep_severe_af=p_rep_severe_af,p_rep_death_af=p_rep_death_af, - p_rep_severe_sa=p_rep_severe_sa,p_rep_death_sa=p_rep_death_sa, - m_FOI_Brazil=m_FOI_Brazil, - deterministic=deterministic,mode_parallel=mode_parallel,cluster=cluster) - - for(set in 1:n_param_sets){ - #if(set %% 10 == 0){cat("\n")} - #cat(set,"\t",sep="") - cat("\n\tSet: ",set,sep="") - log_params_prop=all_param_sets[set,] - - cat("\n\tParams: ",signif(log_params_prop,3)) - - #cat(log_params_prop,file=progress_file,sep="\t",append=TRUE) - - names(log_params_prop)=param_names - posterior_value=single_posterior_calc2(log_params_prop,input_data,obs_sero_data,obs_case_data,consts) - results<-rbind(results,c(set,exp(log_params_prop),posterior_value)) - if(set==1){colnames(results)=c("set",param_names,"posterior")} - - cat("\n\tPosterior likelihood = ",posterior_value,sep="") - - } - results<-results[order(results$posterior,decreasing=TRUE), ] - best_fit_results[[iteration]]=results - - log_params_min_new=log_params_max_new=rep(0,n_params) - for(i in 1:n_params){ - log_params_min_new[i]=min(log(results[c(1:n_bounds),i+1])) - log_params_max_new[i]=max(log(results[c(1:n_bounds),i+1])) - } - names(log_params_min_new)=names(log_params_max_new)=param_names - - matplot(x=c(1:n_params),y=log(t(results[c(1:n_bounds),c(1:n_params)+1])),type="p",pch=16,col=1, - xaxt="n",yaxt="n",xlab="",ylab="",ylim=ylim) - axis(side=1,at=c(1:n_params),labels=xlabels,las=2,cex.axis=0.7) - axis(side=2,at=log(ylabels),labels=ylabels) - matplot(x=c(1:n_params),y=log_params_min,type="l",col=1,lty=2,add=TRUE) - matplot(x=c(1:n_params),y=log_params_max,type="l",col=1,lty=2,add=TRUE) - matplot(x=c(1:n_params),y=log_params_min_new,type="l",col=2,add=TRUE) - matplot(x=c(1:n_params),y=log_params_max_new,type="l",col=2,add=TRUE) - - log_params_min=log_params_min_new - log_params_max=log_params_max_new - } - - return(best_fit_results) -} -#------------------------------------------------------------------------------- -Generate_Dataset2 <- function(input_data = list(),FOI_values = c(),R0_values = c(),sero_template = NULL, - case_template = NULL, vaccine_efficacy = 1.0, - p_severe_inf = 0.12, p_death_severe_inf = 0.39, - p_rep_severe_af = 1.0,p_rep_death_af = 1.0,p_rep_severe_sa = 1.0,p_rep_death_sa = 1.0, - mode_start = 1,start_SEIRV = NULL, dt = 1.0,n_reps = 1, deterministic = FALSE, - mode_parallel = "none",cluster = NULL,output_frame=FALSE){ - - assert_that(input_data_check(input_data),msg=paste("Input data must be in standard format", - " (see https://mrc-ide.github.io/YEP/articles/CGuideAInputs.html)")) - assert_that(any(is.null(sero_template)==FALSE,is.null(case_template)==FALSE),msg="Need at least one template") - if(is.null(sero_template)==FALSE){ - assert_that(all(c("region","year","age_min","age_max","samples","vc_factor") %in% names(sero_template))) - } - if(is.null(case_template)==FALSE){ - assert_that(all(c("region","year") %in% names(case_template))) - assert_that(p_severe_inf>=0.0 && p_severe_inf<=1.0,msg="Severe infection rate must be between 0-1") - assert_that(p_death_severe_inf>=0.0 && p_death_severe_inf<=1.0, - msg="Fatality rate of severe infections must be between 0-1") - assert_that(all(c(p_rep_severe_af,p_rep_severe_sa,p_rep_death_af,p_rep_death_sa)>=0.0)) - assert_that(all(c(p_rep_severe_af,p_rep_severe_sa,p_rep_death_af,p_rep_death_sa)<=1.0)) - } - assert_that(mode_parallel %in% c("none","pars_multi","clusterMap")) - if(mode_parallel=="clusterMap"){assert_that(is.null(cluster)==FALSE)} - - #Prune input data based on regions - regions=regions_breakdown(c(sero_template$region,case_template$region)) - input_data=input_data_truncate(input_data,regions) - n_regions=length(input_data$region_labels) - - #Cross-reference templates with input regions - if(is.null(sero_template)==FALSE){ - xref_sero=template_region_xref(sero_template,input_data$region_labels) - sero_line_list=xref_sero$line_list - } else { - xref_sero=data.frame(year_data_begin=rep(Inf,n_regions),year_end=rep(-Inf,n_regions)) - sero_line_list=rep(NA,n_regions) - } - if(is.null(case_template)==FALSE){ - country_list_af=c("AGO","BDI","BEN","BFA","CAF","CIV","CMR","COD","COG","ERI","ETH","GAB","GHA","GIN","GMB","GNB","GNQ", - "KEN","LBR","MLI","MRT","NER","NGA","RWA","SDN","SEN","SLE","SOM","SSD","TCD","TGO","TZA","UGA","ZMB") - country_list_sa=c("ARG","BOL","BRA","COL","ECU","GUF","GUY","PER","PRY","SUR","VEN") - countries=substr(regions,1,3) - assert_that(all(countries %in% c(country_list_af,country_list_sa))) - xref_case=template_region_xref(case_template,input_data$region_labels) - case_line_list=xref_case$line_list - } else { - xref_case=data.frame(year_data_begin=rep(Inf,n_regions),year_end=rep(-Inf,n_regions)) - case_line_list=rep(NA,n_regions) - } - year_data_begin=year_end=rep(NA,length(input_data$region_labels)) - for(i in 1:length(year_data_begin)){ - year_data_begin[i]=min(xref_sero$year_data_begin[i],xref_case$year_data_begin[i]) - year_end[i]=max(xref_sero$year_end[i],xref_case$year_end[i]) - } - - inv_reps=1/n_reps - assert_that(length(FOI_values)==n_regions,msg="Length of FOI_values must match number of regions to be modelled") - assert_that(length(R0_values)==n_regions,msg="Length of R0_values must match number of regions to be modelled") - if(mode_start==2){assert_that(length(start_SEIRV)==n_regions, - msg="Number of start_SEIRV datasets must match number of regions")} - - #Set up data structures to take modelled data corresponding to observed data - if(is.null(sero_template)){model_sero_data=NULL} else { - blank=rep(0,nrow(sero_template)) - model_sero_data=data.frame(samples=blank,positives=blank,sero=blank) - } - if(is.null(case_template)){model_case_values=model_death_values=NA} else { - model_case_values=model_death_values=rep(0,nrow(case_template)) - } - - #Set up vector of output types to get from model if needed - if(mode_parallel %in% c("none","clusterMap")){ - output_types=rep(NA,n_regions) - for(n_region in 1:n_regions){ - if(is.na(case_line_list[[n_region]][1])==FALSE){ - if(is.na(sero_line_list[[n_region]][1])==FALSE){ - output_types[n_region]="case+sero" - } else{ - output_types[n_region]="case" - } - } else {output_types[n_region]="sero"} - } - } - - #Model all regions in parallel if parallel modes in use - if(mode_parallel=="pars_multi"){ - years_data_all=c(min(year_data_begin):max(year_end)) - if(is.null(sero_template)==FALSE){if(is.null(case_template)==FALSE){output_type="case+sero"} else {output_type="sero"} - } else {output_type="case"} - model_output_all=Model_Run_Multi_Input(FOI_spillover = FOI_values,R0 = R0_values, - vacc_data = input_data$vacc_data, pop_data = input_data$pop_data, - years_data = years_data_all, start_SEIRV=start_SEIRV,output_type = output_type, - year0 = input_data$years_labels[1],mode_start = mode_start, - vaccine_efficacy = vaccine_efficacy, dt = dt, n_particles = n_reps, - n_threads = n_reps*n_regions,deterministic = deterministic) - } - if(mode_parallel=="clusterMap"){ - vacc_data_subsets=pop_data_subsets=years_data_sets=list() #TODO - change input data? - for(n_region in 1:n_regions){ - vacc_data_subsets[[n_region]]=input_data$vacc_data[n_region,,] - pop_data_subsets[[n_region]]=input_data$pop_data[n_region,,] - years_data_sets[[n_region]]=c(year_data_begin[n_region]:year_end[n_region]) - } - if(is.null(start_SEIRV)){start_SEIRV=rep(NA,n_regions)} - model_output_all=clusterMap(cl = cluster,fun = Model_Run, FOI_spillover = FOI_values, R0 = R0_values, - vacc_data = vacc_data_subsets,pop_data = pop_data_subsets, - years_data = years_data_sets, start_SEIRV = start_SEIRV, output_type = output_types, - MoreArgs=list(year0 = input_data$years_labels[1],mode_start = mode_start, - vaccine_efficacy = vaccine_efficacy, dt = dt, n_particles = n_reps, - n_threads = 1 ,deterministic = deterministic)) - } - #if(mode_parallel=="hybrid") #Potential future option combining parallelization types - - #Save relevant output data from each region - for(n_region in 1:n_regions){ - - #Run model if not using parallelization - if(mode_parallel=="none"){ - #cat("\n\t\tBeginning modelling region ",input_data$region_labels[n_region]) - model_output = Model_Run(FOI_spillover = FOI_values[n_region],R0 = R0_values[n_region], - vacc_data = input_data$vacc_data[n_region,,],pop_data = input_data$pop_data[n_region,,], - years_data = c(year_data_begin[n_region]:year_end[n_region]), - start_SEIRV=start_SEIRV[[n_region]],output_type = output_types[n_region], - year0 = input_data$years_labels[1],mode_start = mode_start, - vaccine_efficacy = vaccine_efficacy, dt = dt, n_particles = n_reps,n_threads = n_reps, - deterministic = deterministic) - #cat("\n\t\tFinished modelling region ",n_region) - } else { - model_output = model_output_all[[n_region]] - } - t_pts=length(model_output$year) - - #Compile case data if needed - if(is.na(case_line_list[[n_region]][1])==FALSE){ - case_line_list_region=case_line_list[[n_region]] - years_case=case_template$year[case_line_list_region] - n_lines=length(case_line_list_region) - - #Get reporting probabilities based on region - if(countries[n_region] %in% country_list_af){ - p_rep_severe=p_rep_severe_af - p_rep_death=p_rep_death_af - } else { - p_rep_severe=p_rep_severe_sa - p_rep_death=p_rep_death_sa - } - - for(n_rep in 1:n_reps){ - rep_cases=rep_deaths=rep(0,n_lines) - for(n_line in 1:n_lines){ - #TODO - Set to adjust reporting probabilities based on region - pts=c(1:t_pts)[model_output$year==years_case[n_line]] - infs=sum(model_output$C[n_rep,pts]) - if(deterministic){ - severe_infs=floor(infs)*p_severe_inf - deaths=severe_infs*p_death_severe_inf - rep_deaths[n_line]=round(deaths*p_rep_death) - rep_cases[n_line]=rep_deaths[n_line]+round((severe_infs-deaths)*p_rep_severe) - - } else { - severe_infs=rbinom(1,floor(infs),p_severe_inf) - deaths=rbinom(1,severe_infs,p_death_severe_inf) - rep_deaths[n_line]=rbinom(1,deaths,p_rep_death) - rep_cases[n_line]=rep_deaths[n_line]+rbinom(1,floor(severe_infs-deaths),p_rep_severe) - } - } - - model_case_values[case_line_list_region]=model_case_values[case_line_list_region]+rep_cases - model_death_values[case_line_list_region]=model_death_values[case_line_list_region]+rep_deaths - } - } - - #Compile seroprevalence data if necessary - if(is.na(sero_line_list[[n_region]][1])==FALSE){ - sero_line_list_region=sero_line_list[[n_region]] - for(n_rep in 1:n_reps){ - sero_results=sero_calculate2(sero_template[sero_line_list_region,],model_output,n_rep) - model_sero_data$samples[sero_line_list_region]=model_sero_data$samples[sero_line_list_region]+sero_results$samples - model_sero_data$positives[sero_line_list_region]=model_sero_data$positives[sero_line_list_region]+sero_results$positives - } - } - } - - if(is.null(sero_template)==FALSE){model_sero_data$sero=model_sero_data$positives/model_sero_data$samples} - if(is.null(case_template)==FALSE){ - model_case_values=model_case_values*inv_reps - model_death_values=model_death_values*inv_reps - } - - if(output_frame) { #Output complete frames of data - return(list(model_sero_data=data.frame(region=sero_template$region,year=sero_template$year, - age_min=sero_template$age_min,age_max=sero_template$age_max, - samples=sero_template$samples,positives=sero_template$samples*model_sero_data$sero), - model_case_data=data.frame(region=case_template$region,year=case_template$year, - cases=model_case_values,deaths=model_death_values))) - } else { #Minimal output for MCMC - return(list(model_sero_values=model_sero_data$sero,model_case_values=model_case_values, - model_death_values=model_death_values)) - } -} +#' # R file for functions used for Markov Chain Monte Carlo fitting (and preliminary maximum-likelihood fitting) in +#' # YEP package - alternate versions +#' #------------------------------------------------------------------------------- +#' #------------------------------------------------------------------------------- +#' #' @title MCMC2 +#' #' +#' #' @description TBA +#' #' +#' #' @details TBA +#' #' +#' #' @param log_params_ini TBA +#' #' @param input_data TBA +#' #' @param obs_sero_data TBA +#' #' @param obs_case_data TBA +#' #' @param filename_prefix TBA +#' #' @param Niter TBA +#' #' @param type TBA +#' #' @param log_params_min TBA +#' #' @param log_params_max TBA +#' #' @param mode_start TBA +#' #' @param prior_type TBA +#' #' @param norm_prior_sd TBA +#' #' @param dt TBA +#' #' @param n_reps TBA +#' #' @param enviro_data TBA +#' #' @param R0_fixed_values TBA +#' #' @param vaccine_efficacy TBA +#' #' @param p_severe_inf TBA +#' #' @param p_death_severe_inf TBA +#' #' @param p_rep_severe_af TBA +#' #' @param p_rep_death_af TBA +#' #' @param p_rep_severe_sa TBA +#' #' @param p_rep_death_sa TBA +#' #' @param m_FOI_Brazil TBA +#' #' @param deterministic TBA +#' #' @param mode_parallel TBA +#' #' @param cluster TBA +#' #' ' +#' #' @export +#' #' +#' MCMC2 <- function(log_params_ini=c(),input_data=list(),obs_sero_data=NULL,obs_case_data=NULL,filename_prefix="Chain", +#' Niter=1,type=NULL,log_params_min=c(),log_params_max=c(),mode_start=0,prior_type="zero", +#' norm_prior_sd=c(30,10,30),dt=1.0,n_reps=1,enviro_data=NULL,R0_fixed_values=NULL, +#' vaccine_efficacy=NULL,p_severe_inf=0.12,p_death_severe_inf=0.39, +#' p_rep_severe_af=NULL,p_rep_death_af=NULL,p_rep_severe_sa=NULL,p_rep_death_sa=NULL, +#' m_FOI_Brazil=1.0,deterministic=FALSE,mode_parallel="none",cluster=NULL){ +#' +#' assert_that(is.logical(deterministic)) +#' assert_that(mode_start %in% c(0,1,3),msg="mode_start must have value 0, 1 or 3") +#' n_params=length(log_params_ini) +#' +#' extra_params=c() +#' if(is.null(vaccine_efficacy)==TRUE){extra_params=append(extra_params,"vaccine_efficacy")} +#' if(is.null(p_rep_severe_af)==TRUE){extra_params=append(extra_params,"p_rep_severe_af")} +#' if(is.null(p_rep_death_af)==TRUE){extra_params=append(extra_params,"p_rep_death_af")} +#' if(is.null(p_rep_severe_sa)==TRUE){extra_params=append(extra_params,"p_rep_severe_sa")} +#' if(is.null(p_rep_death_sa)==TRUE){extra_params=append(extra_params,"p_rep_death_sa")} +#' if(is.null(m_FOI_Brazil)==TRUE){extra_params=append(extra_params,"m_FOI_Brazil")} +#' +#' #Process input data to check that all regions with sero and/or case data supplied are present, remove +#' #regions without any supplied data, and add cross-referencing tables for use when calculating likelihood. Take +#' #subset of environmental data (if used) and check that environmental data available for all regions +#' input_data=input_data_process(input_data,obs_sero_data,obs_case_data) +#' regions=names(table(input_data$region_labels)) #Regions in new processed input data list +#' n_regions=length(regions) +#' if(is.null(enviro_data)==FALSE){ +#' for(region in regions){assert_that(region %in% enviro_data$region)} +#' enviro_data=subset(enviro_data,enviro_data$region %in% regions) +#' } +#' +#' #Label parameters according to order and fitting type +#' param_names=create_param_labels(type,input_data,enviro_data,extra_params) +#' names(log_params_ini)=param_names +#' +#' #Set up list of invariant parameter values to supply to other functions +#' consts=list(type=type,log_params_min=log_params_min,log_params_max=log_params_max,mode_start=mode_start, +#' prior_type=prior_type,norm_prior_sd=norm_prior_sd, +#' dt=dt,n_reps=n_reps,enviro_data=enviro_data,R0_fixed_values=R0_fixed_values, +#' vaccine_efficacy=vaccine_efficacy,p_severe_inf=p_severe_inf,p_death_severe_inf=p_death_severe_inf, +#' p_rep_severe_af=p_rep_severe_af,p_rep_death_af=p_rep_death_af, +#' p_rep_severe_sa=p_rep_severe_sa,p_rep_death_sa=p_rep_death_sa, +#' m_FOI_Brazil=m_FOI_Brazil,deterministic=deterministic, +#' mode_parallel=mode_parallel,cluster=cluster) +#' +#' #MCMC setup +#' chain=chain_prop=posterior_current=posterior_prop=flag_accept=chain_cov_all=NULL +#' burnin = min(2*n_params, Niter) +#' fileIndex = 0 +#' log_params=log_params_ini +#' chain_cov=1 +#' adapt=0 +#' posterior_value_current=-Inf +#' +#' #Iterative estimation +#' for (iter in 1:Niter){ +#' +#' #cat("\n\tGenerating new parameter values") +#' #Propose new parameter values +#' log_params_prop=param_prop_setup(log_params,chain_cov,adapt) +#' +#' #cat("\n\tCalculating likelihood") +#' #Calculate likelihood using single_posterior_calc2 function +#' posterior_value_prop=single_posterior_calc2(log_params_prop,input_data,obs_sero_data,obs_case_data,consts) +#' #cat("\n\tLikelihood calculated") +#' +#' if(is.finite(posterior_value_prop)==FALSE) { +#' p_accept = -Inf +#' } else { +#' p_accept = posterior_value_prop - posterior_value_current +#' if(is.na(p_accept) ){ p_accept = -Inf} +#' } +#' +#' ## accept/reject step: +#' tmp = runif(1) +#' if(tmp=10){ +#' filename=paste0(filename_prefix,fileIndex,".csv") +#' } else { +#' filename=paste0(filename_prefix,"0",fileIndex,".csv") +#' } +#' if(file.exists(filename)==FALSE){file.create(filename)} +#' lines=min((fileIndex * 10000+1),iter):iter +#' #cat("\nIteration ",iter,sep="") +#' data_out<-cbind(posterior_current,posterior_prop,exp(chain),flag_accept,exp(chain_prop),chain_cov_all)[lines,] +#' if(fileAccess(filename,2)==0){write.csv(data_out,filename,row.names=FALSE)} +#' } +#' +#' #Decide whether next iteration will be adaptive +#' if (iter>burnin & runif(1)<0.9){ #adapt +#' adapt = 1 +#' chain_cov = cov(chain[max(nrow(chain)-10000, 1):nrow(chain),]) +#' } else { +#' adapt = 0 +#' chain_cov = 1 +#' } +#' } +#' +#' #Get final parameter values +#' param_out=exp(log_params) +#' names(param_out)=names(log_params_ini) +#' +#' return(param_out) +#' } +#' #------------------------------------------------------------------------------- +#' #' @title single_posterior_calc2 +#' #' +#' #' @description TBA +#' #' +#' #' @details TBA +#' #' +#' #' @param log_params_prop TBA +#' #' @param input_data TBA +#' #' @param obs_sero_data TBA +#' #' @param obs_case_data TBA +#' #' @param consts TBA +#' #' +#' #' @export +#' #' +#' single_posterior_calc2 <- function(log_params_prop=c(),input_data=list(),obs_sero_data=NULL,obs_case_data=NULL, +#' consts=list()) { +#' +#' regions=input_data$region_labels +#' n_regions=length(regions) +#' +#' #Get vaccine efficacy and calculate associated prior +#' if(is.numeric(consts$vaccine_efficacy)==FALSE){ +#' vaccine_efficacy=exp(log_params_prop[names(log_params_prop)=="vaccine_efficacy"]) +#' prior_vacc=log(dtrunc(vaccine_efficacy,"norm",a=0,b=1,mean=0.975,sd=0.05)) +#' } else { +#' vaccine_efficacy=consts$vaccine_efficacy +#' prior_vacc=0 +#' } +#' +#' #Get reporting probabilities and check they are within specified bounds +#' prior_report=0 +#' if(is.numeric(consts$p_rep_severe_af)==FALSE){ +#' p_rep_severe_af=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_severe_af"])) +#' if(p_rep_severe_af<=0.0 || p_rep_severe_af>1.0){prior_report=-Inf} +#' } else { +#' p_rep_severe_af=consts$p_rep_severe_af +#' } +#' if(is.numeric(consts$p_rep_severe_sa)==FALSE){ +#' p_rep_severe_sa=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_severe_sa"])) +#' if(p_rep_severe_sa<=0.0 || p_rep_severe_sa>1.0){prior_report=-Inf} +#' } else { +#' p_rep_severe_sa=consts$p_rep_severe_sa +#' } +#' if(is.numeric(consts$p_rep_death_af)==FALSE){ +#' p_rep_death_af=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_death_af"])) +#' if(p_rep_death_af<=0.0 || p_rep_death_af>1.0){prior_report=-Inf} +#' } else { +#' p_rep_death_af=consts$p_rep_death_af +#' } +#' if(is.numeric(consts$p_rep_death_sa)==FALSE){ +#' p_rep_death_sa=as.numeric(exp(log_params_prop[names(log_params_prop)=="p_rep_death_sa"])) +#' if(p_rep_death_sa<=0.0 || p_rep_death_sa>1.0){prior_report=-Inf} +#' } else { +#' p_rep_death_sa=consts$p_rep_death_sa +#' } +#' +#' #Get value of Brazil multiplier and check bounds +#' if(is.numeric(consts$m_FOI_Brazil)==FALSE){ +#' m_FOI_Brazil=as.numeric(exp(log_params_prop[names(log_params_prop)=="m_FOI_Brazil"])) +#' if(m_FOI_Brazil<=0.0 || m_FOI_Brazil>1.0){prior_report=-Inf} +#' } else { +#' m_FOI_Brazil=consts$m_FOI_Brazil +#' } +#' +#' #Get FOI and R0 values and calculate associated prior +#' FOI_R0_data=mcmc_FOI_R0_setup(consts$type,consts$prior_type,consts$norm_prior_sd,regions,log_params_prop, +#' consts$enviro_data,consts$R0_fixed_values, +#' consts$log_params_min,consts$log_params_max) +#' FOI_values=FOI_R0_data$FOI_values +#' for(n_region in 1:n_regions){ +#' if(substr(regions[n_region],1,3)=="BRA"){FOI_values[n_region]=FOI_values[n_region]*m_FOI_Brazil} +#' } +#' R0_values=FOI_R0_data$R0_values +#' prior_prop=prior_vacc+prior_report+FOI_R0_data$prior + +#' sum(dnorm(log(c(p_rep_severe_af,p_rep_death_af,p_rep_severe_sa,p_rep_death_sa)), +#' mean = 0,sd = consts$norm_prior_sd[3],log = TRUE)) +#' +#' if(is.null(obs_sero_data)){sero_like_values=0} +#' if(is.null(obs_case_data)){cases_like_values=deaths_like_values=0} +#' #cat("\n",prior_prop,sep="") +#' +#' #cat("\n\t\tGenerating dataset for regions: ") +#' ### If prior finite, evaluate likelihood ### +#' if (is.finite(prior_prop)) { +#' +#' #cat("\n\tGenerating dataset") +#' #Generate modelled data over all regions +#' dataset <- Generate_Dataset2(input_data,FOI_values,R0_values,obs_sero_data,obs_case_data,vaccine_efficacy,consts$p_severe_inf, +#' consts$p_death_severe_inf,p_rep_severe_af,p_rep_death_af,p_rep_severe_sa,p_rep_death_sa,consts$mode_start,start_SEIRV=NULL,consts$dt, +#' consts$n_reps,consts$deterministic,consts$mode_parallel,consts$cluster) +#' #cat("\n\tDataset generation completed") +#' +#' #Likelihood of observing serological data +#' if(is.null(obs_sero_data)==FALSE){ +#' sero_like_values=sero_data_compare(dataset$model_sero_values,obs_sero_data) +#' } +#' #Likelihood of observing annual case/death data +#' if(is.null(obs_case_data)==FALSE){ +#' cases_like_values=case_data_compare(dataset$model_case_values,obs_case_data$cases) +#' if(is.null(obs_case_data$deaths)==FALSE){ +#' deaths_like_values=case_data_compare(dataset$model_death_values,obs_case_data$deaths) +#' } else {deaths_like_values=0} +#' } +#' +#' posterior=prior_prop+mean(c(sum(sero_like_values,na.rm=TRUE),sum(cases_like_values,na.rm=TRUE), +#' sum(deaths_like_values,na.rm=TRUE)),na.rm=TRUE) +#' # dataset<-NULL +#' # gc() +#' +#' } else {posterior=-Inf} +#' +#' return(posterior) +#' } +#' #------------------------------------------------------------------------------- +#' #' @title mcmc_prelim_fit2 +#' #' +#' #' @description TBA +#' #' +#' #' @details TBA +#' #' +#' #' @param n_iterations TBA +#' #' @param n_param_sets TBA +#' #' @param n_bounds TBA +#' #' @param type TBA +#' #' @param log_params_min TBA +#' #' @param log_params_max TBA +#' #' @param input_data TBA +#' #' @param obs_sero_data TBA +#' #' @param obs_case_data TBA +#' #' @param mode_start TBA +#' #' @param prior_type TBA +#' #' @param norm_prior_sd TBA +#' #' @param dt TBA +#' #' @param n_reps TBA +#' #' @param enviro_data TBA +#' #' @param R0_fixed_values TBA +#' #' @param vaccine_efficacy TBA +#' #' @param p_severe_inf TBA +#' #' @param p_death_severe_inf TBA +#' #' @param p_rep_severe_af TBA +#' #' @param p_rep_death_af TBA +#' #' @param p_rep_severe_sa TBA +#' #' @param p_rep_death_sa TBA +#' #' @param m_FOI_Brazil TBA +#' #' @param deterministic TBA +#' #' @param mode_parallel TBA +#' #' @param cluster TBA +#' #' ' +#' #' @export +#' #' +#' mcmc_prelim_fit2 <- function(n_iterations=1,n_param_sets=1,n_bounds=1,type=NULL,log_params_min=NULL, +#' log_params_max=NULL,input_data=list(),obs_sero_data=list(),obs_case_data=list(), +#' mode_start=0,prior_type="zero",norm_prior_sd=c(30,10,30), +#' dt=1.0,n_reps=1,enviro_data=NULL,R0_fixed_values=c(),vaccine_efficacy=NULL, +#' p_severe_inf = 0.12, p_death_severe_inf = 0.39, +#' p_rep_severe_af=NULL,p_rep_death_af=NULL,p_rep_severe_sa=NULL,p_rep_death_sa=NULL, +#' m_FOI_Brazil=1.0,deterministic=TRUE,mode_parallel="none",cluster=NULL){ +#' +#' #TODO - Add assertthat functions +#' assert_that(mode_start %in% c(0,1,3),msg="mode_start must have value 0, 1 or 3") +#' assert_that(length(log_params_min)==length(log_params_max),msg="Parameter limit vectors must have same lengths") +#' assert_that(type %in% c("FOI+R0","FOI","FOI+R0 enviro","FOI enviro")) +#' assert_that(prior_type %in% c("zero","exp","norm")) +#' +#' best_fit_results=list() +#' n_params=length(log_params_min) +#' extra_params=c() +#' if(is.null(vaccine_efficacy)==TRUE){extra_params=append(extra_params,"vaccine_efficacy")} +#' if(is.null(p_rep_severe_af)==TRUE){extra_params=append(extra_params,"p_rep_severe_af")} +#' if(is.null(p_rep_death_af)==TRUE){extra_params=append(extra_params,"p_rep_death_af")} +#' if(is.null(p_rep_severe_sa)==TRUE){extra_params=append(extra_params,"p_rep_severe_sa")} +#' if(is.null(p_rep_death_sa)==TRUE){extra_params=append(extra_params,"p_rep_death_sa")} +#' if(is.null(m_FOI_Brazil)==TRUE){extra_params=append(extra_params,"m_FOI_Brazil")} +#' param_names=create_param_labels(type,input_data,enviro_data,extra_params) +#' +#' #TODO - Additional assert_that checks +#' assert_that(length(param_names)==n_params) +#' names(log_params_min)=names(log_params_max)=param_names +#' xlabels=param_names +#' for(i in 1:n_params){xlabels[i]=substr(xlabels[i],1,15)} +#' ylabels=10^c(-8,-6,-4,-3,-2,-1,0,1) +#' par(mar=c(6,2,1,1)) +#' ylim=c(min(log_params_min),max(log_params_max)) +#' +#' for(iteration in 1:n_iterations){ +#' cat("\nIteration: ",iteration,"\n",sep="") +#' all_param_sets <- lhs(n=n_param_sets,rect=cbind(log_params_min,log_params_max)) +#' results=data.frame() +#' consts=list(type=type,log_params_min=log_params_min,log_params_max=log_params_max, +#' mode_start=mode_start,prior_type=prior_type,norm_prior_sd=norm_prior_sd, +#' dt=dt,n_reps=n_reps,enviro_data=enviro_data, +#' R0_fixed_values=R0_fixed_values,vaccine_efficacy=vaccine_efficacy, +#' p_severe_inf = p_severe_inf, p_death_severe_inf = p_death_severe_inf, +#' p_rep_severe_af=p_rep_severe_af,p_rep_death_af=p_rep_death_af, +#' p_rep_severe_sa=p_rep_severe_sa,p_rep_death_sa=p_rep_death_sa, +#' m_FOI_Brazil=m_FOI_Brazil, +#' deterministic=deterministic,mode_parallel=mode_parallel,cluster=cluster) +#' +#' for(set in 1:n_param_sets){ +#' #if(set %% 10 == 0){cat("\n")} +#' #cat(set,"\t",sep="") +#' cat("\n\tSet: ",set,sep="") +#' log_params_prop=all_param_sets[set,] +#' +#' cat("\n\tParams: ",signif(log_params_prop,3)) +#' +#' #cat(log_params_prop,file=progress_file,sep="\t",append=TRUE) +#' +#' names(log_params_prop)=param_names +#' posterior_value=single_posterior_calc2(log_params_prop,input_data,obs_sero_data,obs_case_data,consts) +#' results<-rbind(results,c(set,exp(log_params_prop),posterior_value)) +#' if(set==1){colnames(results)=c("set",param_names,"posterior")} +#' +#' cat("\n\tPosterior likelihood = ",posterior_value,sep="") +#' +#' } +#' results<-results[order(results$posterior,decreasing=TRUE), ] +#' best_fit_results[[iteration]]=results +#' +#' log_params_min_new=log_params_max_new=rep(0,n_params) +#' for(i in 1:n_params){ +#' log_params_min_new[i]=min(log(results[c(1:n_bounds),i+1])) +#' log_params_max_new[i]=max(log(results[c(1:n_bounds),i+1])) +#' } +#' names(log_params_min_new)=names(log_params_max_new)=param_names +#' +#' matplot(x=c(1:n_params),y=log(t(results[c(1:n_bounds),c(1:n_params)+1])),type="p",pch=16,col=1, +#' xaxt="n",yaxt="n",xlab="",ylab="",ylim=ylim) +#' axis(side=1,at=c(1:n_params),labels=xlabels,las=2,cex.axis=0.7) +#' axis(side=2,at=log(ylabels),labels=ylabels) +#' matplot(x=c(1:n_params),y=log_params_min,type="l",col=1,lty=2,add=TRUE) +#' matplot(x=c(1:n_params),y=log_params_max,type="l",col=1,lty=2,add=TRUE) +#' matplot(x=c(1:n_params),y=log_params_min_new,type="l",col=2,add=TRUE) +#' matplot(x=c(1:n_params),y=log_params_max_new,type="l",col=2,add=TRUE) +#' +#' log_params_min=log_params_min_new +#' log_params_max=log_params_max_new +#' } +#' +#' return(best_fit_results) +#' } +#' #------------------------------------------------------------------------------- +#' Generate_Dataset2 <- function(input_data = list(),FOI_values = c(),R0_values = c(),sero_template = NULL, +#' case_template = NULL, vaccine_efficacy = 1.0, +#' p_severe_inf = 0.12, p_death_severe_inf = 0.39, +#' p_rep_severe_af = 1.0,p_rep_death_af = 1.0,p_rep_severe_sa = 1.0,p_rep_death_sa = 1.0, +#' mode_start = 1,start_SEIRV = NULL, dt = 1.0,n_reps = 1, deterministic = FALSE, +#' mode_parallel = "none",cluster = NULL,output_frame=FALSE){ +#' +#' assert_that(input_data_check(input_data),msg=paste("Input data must be in standard format", +#' " (see https://mrc-ide.github.io/YEP/articles/CGuideAInputs.html)")) +#' assert_that(any(is.null(sero_template)==FALSE,is.null(case_template)==FALSE),msg="Need at least one template") +#' if(is.null(sero_template)==FALSE){ +#' assert_that(all(c("region","year","age_min","age_max","samples","vc_factor") %in% names(sero_template))) +#' } +#' if(is.null(case_template)==FALSE){ +#' assert_that(all(c("region","year") %in% names(case_template))) +#' assert_that(p_severe_inf>=0.0 && p_severe_inf<=1.0,msg="Severe infection rate must be between 0-1") +#' assert_that(p_death_severe_inf>=0.0 && p_death_severe_inf<=1.0, +#' msg="Fatality rate of severe infections must be between 0-1") +#' assert_that(all(c(p_rep_severe_af,p_rep_severe_sa,p_rep_death_af,p_rep_death_sa)>=0.0)) +#' assert_that(all(c(p_rep_severe_af,p_rep_severe_sa,p_rep_death_af,p_rep_death_sa)<=1.0)) +#' } +#' assert_that(mode_parallel %in% c("none","pars_multi","clusterMap")) +#' if(mode_parallel=="clusterMap"){assert_that(is.null(cluster)==FALSE)} +#' +#' #Prune input data based on regions +#' regions=regions_breakdown(c(sero_template$region,case_template$region)) +#' input_data=input_data_truncate(input_data,regions) +#' n_regions=length(input_data$region_labels) +#' +#' #Cross-reference templates with input regions +#' if(is.null(sero_template)==FALSE){ +#' xref_sero=template_region_xref(sero_template,input_data$region_labels) +#' sero_line_list=xref_sero$line_list +#' } else { +#' xref_sero=data.frame(year_data_begin=rep(Inf,n_regions),year_end=rep(-Inf,n_regions)) +#' sero_line_list=rep(NA,n_regions) +#' } +#' if(is.null(case_template)==FALSE){ +#' country_list_af=c("AGO","BDI","BEN","BFA","CAF","CIV","CMR","COD","COG","ERI","ETH","GAB","GHA","GIN","GMB","GNB","GNQ", +#' "KEN","LBR","MLI","MRT","NER","NGA","RWA","SDN","SEN","SLE","SOM","SSD","TCD","TGO","TZA","UGA","ZMB") +#' country_list_sa=c("ARG","BOL","BRA","COL","ECU","GUF","GUY","PER","PRY","SUR","VEN") +#' countries=substr(regions,1,3) +#' assert_that(all(countries %in% c(country_list_af,country_list_sa))) +#' xref_case=template_region_xref(case_template,input_data$region_labels) +#' case_line_list=xref_case$line_list +#' } else { +#' xref_case=data.frame(year_data_begin=rep(Inf,n_regions),year_end=rep(-Inf,n_regions)) +#' case_line_list=rep(NA,n_regions) +#' } +#' year_data_begin=year_end=rep(NA,length(input_data$region_labels)) +#' for(i in 1:length(year_data_begin)){ +#' year_data_begin[i]=min(xref_sero$year_data_begin[i],xref_case$year_data_begin[i]) +#' year_end[i]=max(xref_sero$year_end[i],xref_case$year_end[i]) +#' } +#' +#' inv_reps=1/n_reps +#' assert_that(length(FOI_values)==n_regions,msg="Length of FOI_values must match number of regions to be modelled") +#' assert_that(length(R0_values)==n_regions,msg="Length of R0_values must match number of regions to be modelled") +#' if(mode_start==2){assert_that(length(start_SEIRV)==n_regions, +#' msg="Number of start_SEIRV datasets must match number of regions")} +#' +#' #Set up data structures to take modelled data corresponding to observed data +#' if(is.null(sero_template)){model_sero_data=NULL} else { +#' blank=rep(0,nrow(sero_template)) +#' model_sero_data=data.frame(samples=blank,positives=blank,sero=blank) +#' } +#' if(is.null(case_template)){model_case_values=model_death_values=NA} else { +#' model_case_values=model_death_values=rep(0,nrow(case_template)) +#' } +#' +#' #Set up vector of output types to get from model if needed +#' if(mode_parallel %in% c("none","clusterMap")){ +#' output_types=rep(NA,n_regions) +#' for(n_region in 1:n_regions){ +#' if(is.na(case_line_list[[n_region]][1])==FALSE){ +#' if(is.na(sero_line_list[[n_region]][1])==FALSE){ +#' output_types[n_region]="case+sero" +#' } else{ +#' output_types[n_region]="case" +#' } +#' } else {output_types[n_region]="sero"} +#' } +#' } +#' +#' #Model all regions in parallel if parallel modes in use +#' if(mode_parallel=="pars_multi"){ +#' years_data_all=c(min(year_data_begin):max(year_end)) +#' if(is.null(sero_template)==FALSE){if(is.null(case_template)==FALSE){output_type="case+sero"} else {output_type="sero"} +#' } else {output_type="case"} +#' model_output_all=Model_Run_Multi_Input(FOI_spillover = FOI_values,R0 = R0_values, +#' vacc_data = input_data$vacc_data, pop_data = input_data$pop_data, +#' years_data = years_data_all, start_SEIRV=start_SEIRV,output_type = output_type, +#' year0 = input_data$years_labels[1],mode_start = mode_start, +#' vaccine_efficacy = vaccine_efficacy, dt = dt, n_particles = n_reps, +#' n_threads = n_reps*n_regions,deterministic = deterministic) +#' } +#' if(mode_parallel=="clusterMap"){ +#' vacc_data_subsets=pop_data_subsets=years_data_sets=list() #TODO - change input data? +#' for(n_region in 1:n_regions){ +#' vacc_data_subsets[[n_region]]=input_data$vacc_data[n_region,,] +#' pop_data_subsets[[n_region]]=input_data$pop_data[n_region,,] +#' years_data_sets[[n_region]]=c(year_data_begin[n_region]:year_end[n_region]) +#' } +#' if(is.null(start_SEIRV)){start_SEIRV=rep(NA,n_regions)} +#' model_output_all=clusterMap(cl = cluster,fun = Model_Run, FOI_spillover = FOI_values, R0 = R0_values, +#' vacc_data = vacc_data_subsets,pop_data = pop_data_subsets, +#' years_data = years_data_sets, start_SEIRV = start_SEIRV, output_type = output_types, +#' MoreArgs=list(year0 = input_data$years_labels[1],mode_start = mode_start, +#' vaccine_efficacy = vaccine_efficacy, dt = dt, n_particles = n_reps, +#' n_threads = 1 ,deterministic = deterministic)) +#' } +#' #if(mode_parallel=="hybrid") #Potential future option combining parallelization types +#' +#' #Save relevant output data from each region +#' for(n_region in 1:n_regions){ +#' +#' #Run model if not using parallelization +#' if(mode_parallel=="none"){ +#' #cat("\n\t\tBeginning modelling region ",input_data$region_labels[n_region]) +#' model_output = Model_Run(FOI_spillover = FOI_values[n_region],R0 = R0_values[n_region], +#' vacc_data = input_data$vacc_data[n_region,,],pop_data = input_data$pop_data[n_region,,], +#' years_data = c(year_data_begin[n_region]:year_end[n_region]), +#' start_SEIRV=start_SEIRV[[n_region]],output_type = output_types[n_region], +#' year0 = input_data$years_labels[1],mode_start = mode_start, +#' vaccine_efficacy = vaccine_efficacy, dt = dt, n_particles = n_reps,n_threads = n_reps, +#' deterministic = deterministic) +#' #cat("\n\t\tFinished modelling region ",n_region) +#' } else { +#' model_output = model_output_all[[n_region]] +#' } +#' t_pts=length(model_output$year) +#' +#' #Compile case data if needed +#' if(is.na(case_line_list[[n_region]][1])==FALSE){ +#' case_line_list_region=case_line_list[[n_region]] +#' years_case=case_template$year[case_line_list_region] +#' n_lines=length(case_line_list_region) +#' +#' #Get reporting probabilities based on region +#' if(countries[n_region] %in% country_list_af){ +#' p_rep_severe=p_rep_severe_af +#' p_rep_death=p_rep_death_af +#' } else { +#' p_rep_severe=p_rep_severe_sa +#' p_rep_death=p_rep_death_sa +#' } +#' +#' for(n_rep in 1:n_reps){ +#' rep_cases=rep_deaths=rep(0,n_lines) +#' for(n_line in 1:n_lines){ +#' #TODO - Set to adjust reporting probabilities based on region +#' pts=c(1:t_pts)[model_output$year==years_case[n_line]] +#' infs=sum(model_output$C[n_rep,pts]) +#' if(deterministic){ +#' severe_infs=floor(infs)*p_severe_inf +#' deaths=severe_infs*p_death_severe_inf +#' rep_deaths[n_line]=round(deaths*p_rep_death) +#' rep_cases[n_line]=rep_deaths[n_line]+round((severe_infs-deaths)*p_rep_severe) +#' +#' } else { +#' severe_infs=rbinom(1,floor(infs),p_severe_inf) +#' deaths=rbinom(1,severe_infs,p_death_severe_inf) +#' rep_deaths[n_line]=rbinom(1,deaths,p_rep_death) +#' rep_cases[n_line]=rep_deaths[n_line]+rbinom(1,floor(severe_infs-deaths),p_rep_severe) +#' } +#' } +#' +#' model_case_values[case_line_list_region]=model_case_values[case_line_list_region]+rep_cases +#' model_death_values[case_line_list_region]=model_death_values[case_line_list_region]+rep_deaths +#' } +#' } +#' +#' #Compile seroprevalence data if necessary +#' if(is.na(sero_line_list[[n_region]][1])==FALSE){ +#' sero_line_list_region=sero_line_list[[n_region]] +#' for(n_rep in 1:n_reps){ +#' sero_results=sero_calculate2(sero_template[sero_line_list_region,],model_output,n_rep) +#' model_sero_data$samples[sero_line_list_region]=model_sero_data$samples[sero_line_list_region]+sero_results$samples +#' model_sero_data$positives[sero_line_list_region]=model_sero_data$positives[sero_line_list_region]+sero_results$positives +#' } +#' } +#' } +#' +#' if(is.null(sero_template)==FALSE){model_sero_data$sero=model_sero_data$positives/model_sero_data$samples} +#' if(is.null(case_template)==FALSE){ +#' model_case_values=model_case_values*inv_reps +#' model_death_values=model_death_values*inv_reps +#' } +#' +#' if(output_frame) { #Output complete frames of data +#' return(list(model_sero_data=data.frame(region=sero_template$region,year=sero_template$year, +#' age_min=sero_template$age_min,age_max=sero_template$age_max, +#' samples=sero_template$samples,positives=sero_template$samples*model_sero_data$sero), +#' model_case_data=data.frame(region=case_template$region,year=case_template$year, +#' cases=model_case_values,deaths=model_death_values))) +#' } else { #Minimal output for MCMC +#' return(list(model_sero_values=model_sero_data$sero,model_case_values=model_case_values, +#' model_death_values=model_death_values)) +#' } +#' } diff --git a/man/MCMC.Rd b/man/MCMC.Rd index 53abb02..ee4260f 100644 --- a/man/MCMC.Rd +++ b/man/MCMC.Rd @@ -12,21 +12,16 @@ MCMC( filename_prefix = "Chain", Niter = 1, type = NULL, - log_params_min = c(), - log_params_max = c(), mode_start = 0, - prior_type = "zero", - norm_prior_sd = c(30, 10, 30), + prior_settings = list(type = "zero"), dt = 1, n_reps = 1, enviro_data = NULL, R0_fixed_values = NULL, - vaccine_efficacy = NULL, p_severe_inf = 0.12, p_death_severe_inf = 0.39, - p_rep_severe = NULL, - p_rep_death = NULL, - m_FOI_Brazil = 1, + add_values = list(vaccine_efficacy = 1, p_rep_severe = 1, p_rep_death = 1, m_FOI_Brazil + = 1), deterministic = FALSE, mode_parallel = "none", cluster = NULL @@ -67,25 +62,19 @@ cases/no. deaths} \item{type}{Type of parameter set (FOI only, FOI+R0, FOI and/or R0 coefficients associated with environmental covariates); choose from "FOI","FOI+R0","FOI enviro","FOI+R0 enviro"} -\item{log_params_min}{Lower limits of varied parameter values if specified} - -\item{log_params_max}{Upper limits of varied parameter values if specified} - \item{mode_start}{Flag indicating how to set initial population immunity level in addition to vaccination If mode_start=0, only vaccinated individuals If mode_start=1, shift some non-vaccinated individuals into recovered to give herd immunity (uniform by age, R0 based only) If mode_start=3, shift some non-vaccinated individuals into recovered to give herd immunity (stratified by age)} -\item{prior_type}{Text indicating which type of calculation to use for prior probability -If prior_type = "zero", prior probability is always zero -If prior_type = "flat", prior probability is zero if FOI/R0 in designated ranges, -Inf otherwise -If prior_type = "exp", prior probability is given by dexp calculation on FOI/R0 values -If prior_type = "norm", prior probability is given by dnorm calculation on parameter values using norm_prior_sd values} - -\item{norm_prior_sd}{Vector of 3 values of standard deviations to use in dnorm calculations -Value 1 = SD to use for FOI values/coefficients if prior_type = "norm" -Value 2 = SD to use for R0 values/coefficients if prior_type = "norm" (if relevant) -Value 3 = SD to use for reporting probability prior calculations for all prior_type values (if relevant)} +\item{prior_settings}{List containing settings for priors: must contain text named "type": +If type = "zero", prior probability is always zero +If type = "flat", prior probability is zero if log parameter values in designated ranges log_params_min and log_params_max, + -Inf otherwise; log_params_min and log_params_max included in prior_settings as vectors of same length as log_params_ini +If type = "exp", prior probability is given by dexp calculation on FOI/R0 values +If type = "norm", prior probability is given by dnorm calculation on parameter values with settings based on vectors of values + in prior_settings; norm_params_mean and norm_params_sd (mean and standard deviation values applied to log FOI/R0 + parameters and to actual values of additional parameters)} \item{dt}{time increment in days (must be 1 or 5)} @@ -96,17 +85,15 @@ Value 3 = SD to use for reporting probability prior calculations for all prior_t \item{R0_fixed_values}{Values of R0 to use if only FOI is subject to fitting (i.e. type set to "FOI" or "FOI enviro"); set to NULL if not in use} -\item{vaccine_efficacy}{Vaccine efficacy (set to NULL if being varied as a parameter)} - \item{p_severe_inf}{Probability of an infection being severe} \item{p_death_severe_inf}{Probability of a severe infection resulting in death} -\item{p_rep_severe}{Probability of observation of severe infection (set to NULL if being varied as a parameter)} - -\item{p_rep_death}{Probability of observation of death (set to NULL if being varied as a parameter)} - -\item{m_FOI_Brazil}{Multiplier of spillover FOI for Brazil regions (set to NULL if being varied as a parameter)} +\item{add_values}{List containing values of additional inputs or NULL to indicate they are being varied as parameters: +vaccine_efficacy Vaccine efficacy (set to NULL if being varied as a parameter) +p_rep_severe Probability of observation of severe infection (set to NULL if being varied as a parameter) +p_rep_death Probability of observation of death (set to NULL if being varied as a parameter) +m_FOI_Brazil Multiplier of spillover FOI for Brazil regions (set to NULL if being varied as a parameter)} \item{deterministic}{TRUE/FALSE - set model to run in deterministic mode if TRUE} diff --git a/man/MCMC2.Rd b/man/MCMC2.Rd deleted file mode 100644 index eb3f574..0000000 --- a/man/MCMC2.Rd +++ /dev/null @@ -1,98 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mcmc2.R -\name{MCMC2} -\alias{MCMC2} -\title{MCMC2} -\usage{ -MCMC2( - log_params_ini = c(), - input_data = list(), - obs_sero_data = NULL, - obs_case_data = NULL, - filename_prefix = "Chain", - Niter = 1, - type = NULL, - log_params_min = c(), - log_params_max = c(), - mode_start = 0, - prior_type = "zero", - norm_prior_sd = c(30, 10, 30), - dt = 1, - n_reps = 1, - enviro_data = NULL, - R0_fixed_values = NULL, - vaccine_efficacy = NULL, - p_severe_inf = 0.12, - p_death_severe_inf = 0.39, - p_rep_severe_af = NULL, - p_rep_death_af = NULL, - p_rep_severe_sa = NULL, - p_rep_death_sa = NULL, - m_FOI_Brazil = 1, - deterministic = FALSE, - mode_parallel = "none", - cluster = NULL -) -} -\arguments{ -\item{log_params_ini}{TBA} - -\item{input_data}{TBA} - -\item{obs_sero_data}{TBA} - -\item{obs_case_data}{TBA} - -\item{filename_prefix}{TBA} - -\item{Niter}{TBA} - -\item{type}{TBA} - -\item{log_params_min}{TBA} - -\item{log_params_max}{TBA} - -\item{mode_start}{TBA} - -\item{prior_type}{TBA} - -\item{norm_prior_sd}{TBA} - -\item{dt}{TBA} - -\item{n_reps}{TBA} - -\item{enviro_data}{TBA} - -\item{R0_fixed_values}{TBA} - -\item{vaccine_efficacy}{TBA} - -\item{p_severe_inf}{TBA} - -\item{p_death_severe_inf}{TBA} - -\item{p_rep_severe_af}{TBA} - -\item{p_rep_death_af}{TBA} - -\item{p_rep_severe_sa}{TBA} - -\item{p_rep_death_sa}{TBA} - -\item{m_FOI_Brazil}{TBA} - -\item{deterministic}{TBA} - -\item{mode_parallel}{TBA} - -\item{cluster}{TBA -'} -} -\description{ -TBA -} -\details{ -TBA -} diff --git a/man/create_param_labels.Rd b/man/create_param_labels.Rd index 2cf3a7e..1b9273f 100644 --- a/man/create_param_labels.Rd +++ b/man/create_param_labels.Rd @@ -8,7 +8,7 @@ create_param_labels( type = "FOI", input_data = list(), enviro_data = NULL, - extra_params = c("vacc_eff") + extra_estimated_params = c("vacc_eff") ) } \arguments{ @@ -20,8 +20,8 @@ code and usually loaded from RDS file)} \item{enviro_data}{Environmental data frame, containing only relevant environmental variables} -\item{extra_params}{Vector of strings listing parameters besides ones determining FOI/R0 (may include vaccine -efficacy and/or infection/death reporting probabilities)} +\item{extra_estimated_params}{Vector of strings listing variable parameters besides ones determining FOI/R0 (may include +vaccine efficacy and/or infection/death reporting probabilities and/or Brazil FOI adjustment factor)} } \description{ Apply names to the parameters in a set used for data matching and parameter fitting diff --git a/man/mcmc_FOI_R0_setup.Rd b/man/mcmc_FOI_R0_setup.Rd index d7fd6b2..9f5498d 100644 --- a/man/mcmc_FOI_R0_setup.Rd +++ b/man/mcmc_FOI_R0_setup.Rd @@ -6,27 +6,18 @@ \usage{ mcmc_FOI_R0_setup( type = "", - prior_type = "", - norm_prior_sd = c(30, 10, 30), + prior_settings = list(type = "zero"), regions = "", log_params_prop = c(), enviro_data = list(), - R0_fixed_values = c(), - log_params_min = c(), - log_params_max = c() + R0_fixed_values = c() ) } \arguments{ \item{type}{Type of parameter set (FOI only, FOI+R0, FOI and/or R0 coefficients associated with environmental covariates); choose from "FOI","FOI+R0","FOI enviro","FOI+R0 enviro"} -\item{prior_type}{Text indicating which type of calculation to use for prior probability -If prior_type = "zero", prior probability is always zero -If prior_type = "flat", prior probability is zero if FOI/R0 in designated ranges, -Inf otherwise -If prior_type = "exp", prior probability is given by dexp calculation on FOI/R0 values -If prior_type = "norm", prior probability is given by dnorm calculation on parameter values} - -\item{norm_prior_sd}{TBA} +\item{prior_settings}{TBA} \item{regions}{Vector of region names} @@ -34,11 +25,7 @@ If prior_type = "norm", prior probability is given by dnorm calculation on param \item{enviro_data}{Environmental data frame, containing only relevant environmental variables} -\item{R0_fixed_values}{Values of R0 to use if not being varied} - -\item{log_params_min}{Lower limits of varied parameter values if specified} - -\item{log_params_max}{Upper limits of varied parameter values if specified +\item{R0_fixed_values}{Values of R0 to use if not being varied '} } \description{ diff --git a/man/mcmc_checks.Rd b/man/mcmc_checks.Rd index 76fe83b..f2dfd49 100644 --- a/man/mcmc_checks.Rd +++ b/man/mcmc_checks.Rd @@ -8,16 +8,12 @@ mcmc_checks( log_params_ini = c(), n_regions = 1, type = NULL, - log_params_min = c(), - log_params_max = c(), - prior_type = NULL, - norm_prior_sd = c(30, 10, 30), + prior_settings = list(type = "zero"), enviro_data = NULL, R0_fixed_values = NULL, - vaccine_efficacy = NULL, - p_rep_severe = NULL, - p_rep_death = NULL, - m_FOI_Brazil = 1 + add_values = list(vaccine_efficacy = 1, p_rep_severe = 1, p_rep_death = 1, m_FOI_Brazil + = 1), + extra_estimated_params = list() ) } \arguments{ @@ -28,32 +24,15 @@ mcmc_checks( \item{type}{Type of parameter set (FOI only, FOI+R0, FOI and/or R0 coefficients associated with environmental covariates); choose from "FOI","FOI+R0","FOI enviro","FOI+R0 enviro"} -\item{log_params_min}{Lower limits of varied parameter values if specified} - -\item{log_params_max}{Upper limits of varied parameter values if specified} - -\item{prior_type}{Text indicating which type of calculation to use for prior probability -If prior_type = "zero", prior probability is always zero -If prior_type = "flat", prior probability is zero if FOI/R0 in designated ranges, -Inf otherwise -If prior_type = "exp", prior probability is given by dexp calculation on FOI/R0 values -If prior_type = "norm", prior probability is given by dnorm calculation on parameter values} - -\item{norm_prior_sd}{Vector of 3 values of standard deviations to use in dnorm calculations -Value 1 = SD to use for FOI values/coefficients if prior_type = "norm" -Value 2 = SD to use for R0 values/coefficients if prior_type = "norm" (if relevant) -Value 3 = SD to use for reporting probability prior calculations for all prior_type values (if relevant)} +\item{prior_settings}{TBA} \item{enviro_data}{Values of environmental covariates (if in use)} \item{R0_fixed_values}{Values of R0 to use if not being varied (set to NULL if R0 is being varied)} -\item{vaccine_efficacy}{Vaccine efficacy (set to NULL if being varied as a parameter)} - -\item{p_rep_severe}{Probability of observation of severe infection (set to NULL if being varied as a parameter)} - -\item{p_rep_death}{Probability of observation of death (set to NULL if being varied as a parameter)} +\item{add_values}{TBA} -\item{m_FOI_Brazil}{FOI multiplier for Brazil regions} +\item{extra_estimated_params}{TBA} } \description{ Perform checks on MCMC inputs diff --git a/man/mcmc_prelim_fit.Rd b/man/mcmc_prelim_fit.Rd index 6012db4..215fcd1 100644 --- a/man/mcmc_prelim_fit.Rd +++ b/man/mcmc_prelim_fit.Rd @@ -15,18 +15,15 @@ mcmc_prelim_fit( obs_sero_data = list(), obs_case_data = list(), mode_start = 0, - prior_type = "zero", - norm_prior_sd = c(30, 10, 30), + prior_settings = list(type = "zero"), dt = 1, n_reps = 1, enviro_data = NULL, R0_fixed_values = c(), - vaccine_efficacy = NULL, p_severe_inf = 0.12, p_death_severe_inf = 0.39, - p_rep_severe = NULL, - p_rep_death = NULL, - m_FOI_Brazil = 1, + add_values = list(vaccine_efficacy = 1, p_rep_severe = 1, p_rep_death = 1, m_FOI_Brazil + = 1), deterministic = TRUE, mode_parallel = "none", cluster = NULL @@ -61,12 +58,7 @@ If mode_start=0, only vaccinated individuals If mode_start=1, shift some non-vaccinated individuals into recovered to give herd immunity (uniform by age, R0 based only) If mode_start=3, shift some non-vaccinated individuals into recovered to give herd immunity (stratified by age)} -\item{prior_type}{Text indicating which type of calculation to use for prior probability -If prior_type = "zero", prior probability is always zero -If prior_type = "exp", prior probability is given by dexp calculation on FOI/R0 values -If prior_type = "norm", prior probability is given by dnorm calculation on parameter values} - -\item{norm_prior_sd}{TBA} +\item{prior_settings}{TBA} \item{dt}{time increment in days (must be 1 or 5)} @@ -76,17 +68,15 @@ If prior_type = "norm", prior probability is given by dnorm calculation on param \item{R0_fixed_values}{Values of R0 to use if not being varied} -\item{vaccine_efficacy}{Vaccine efficacy (set to NULL if being varied as a parameter)} - \item{p_severe_inf}{TBA} \item{p_death_severe_inf}{TBA} -\item{p_rep_severe}{Probability of observation of severe infection (set to NULL if being varied as a parameter)} - -\item{p_rep_death}{Probability of observation of death (set to NULL if being varied as a parameter)} - -\item{m_FOI_Brazil}{Multiplier of FOI in Brazil regions} +\item{add_values}{List containing values of additional inputs or NULL to indicate they are being varied as parameters: +vaccine_efficacy Vaccine efficacy (set to NULL if being varied as a parameter) +p_rep_severe Probability of observation of severe infection (set to NULL if being varied as a parameter) +p_rep_death Probability of observation of death (set to NULL if being varied as a parameter) +m_FOI_Brazil Multiplier of spillover FOI for Brazil regions (set to NULL if being varied as a parameter)} \item{deterministic}{TBA} diff --git a/man/mcmc_prelim_fit2.Rd b/man/mcmc_prelim_fit2.Rd deleted file mode 100644 index 97e9341..0000000 --- a/man/mcmc_prelim_fit2.Rd +++ /dev/null @@ -1,98 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mcmc2.R -\name{mcmc_prelim_fit2} -\alias{mcmc_prelim_fit2} -\title{mcmc_prelim_fit2} -\usage{ -mcmc_prelim_fit2( - n_iterations = 1, - n_param_sets = 1, - n_bounds = 1, - type = NULL, - log_params_min = NULL, - log_params_max = NULL, - input_data = list(), - obs_sero_data = list(), - obs_case_data = list(), - mode_start = 0, - prior_type = "zero", - norm_prior_sd = c(30, 10, 30), - dt = 1, - n_reps = 1, - enviro_data = NULL, - R0_fixed_values = c(), - vaccine_efficacy = NULL, - p_severe_inf = 0.12, - p_death_severe_inf = 0.39, - p_rep_severe_af = NULL, - p_rep_death_af = NULL, - p_rep_severe_sa = NULL, - p_rep_death_sa = NULL, - m_FOI_Brazil = 1, - deterministic = TRUE, - mode_parallel = "none", - cluster = NULL -) -} -\arguments{ -\item{n_iterations}{TBA} - -\item{n_param_sets}{TBA} - -\item{n_bounds}{TBA} - -\item{type}{TBA} - -\item{log_params_min}{TBA} - -\item{log_params_max}{TBA} - -\item{input_data}{TBA} - -\item{obs_sero_data}{TBA} - -\item{obs_case_data}{TBA} - -\item{mode_start}{TBA} - -\item{prior_type}{TBA} - -\item{norm_prior_sd}{TBA} - -\item{dt}{TBA} - -\item{n_reps}{TBA} - -\item{enviro_data}{TBA} - -\item{R0_fixed_values}{TBA} - -\item{vaccine_efficacy}{TBA} - -\item{p_severe_inf}{TBA} - -\item{p_death_severe_inf}{TBA} - -\item{p_rep_severe_af}{TBA} - -\item{p_rep_death_af}{TBA} - -\item{p_rep_severe_sa}{TBA} - -\item{p_rep_death_sa}{TBA} - -\item{m_FOI_Brazil}{TBA} - -\item{deterministic}{TBA} - -\item{mode_parallel}{TBA} - -\item{cluster}{TBA -'} -} -\description{ -TBA -} -\details{ -TBA -} diff --git a/man/single_posterior_calc.Rd b/man/single_posterior_calc.Rd index fc8c24c..c8fd914 100644 --- a/man/single_posterior_calc.Rd +++ b/man/single_posterior_calc.Rd @@ -25,9 +25,9 @@ positives} \item{obs_case_data}{Annual reported case/death data for comparison, by region and year, in format no. cases/no. deaths} -\item{consts}{= List of constant parameters/flags/etc. loaded to mcmc() (type,log_params_min,log_params_max, -mode_start,prior_type,norm_prior_sd,dt,n_reps,enviro_data,R0_fixed_values,vaccine_efficacy,p_severe_inf, -p_death_severe_inf,p_rep_severe,p_rep_death,m_FOI_Brazil,deterministic, mode_parallel, cluster)} +\item{consts}{= List of constant parameters/flags/etc. loaded to mcmc() (type, +mode_start,prior_settings,dt,n_reps,enviro_data,R0_fixed_values,p_severe_inf, +p_death_severe_inf,add_values list,extra_estimated_params,deterministic, mode_parallel, cluster)} } \description{ Function which calculates and outputs posterior likelihood of observing simulated data diff --git a/man/single_posterior_calc2.Rd b/man/single_posterior_calc2.Rd deleted file mode 100644 index 4ee3138..0000000 --- a/man/single_posterior_calc2.Rd +++ /dev/null @@ -1,31 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mcmc2.R -\name{single_posterior_calc2} -\alias{single_posterior_calc2} -\title{single_posterior_calc2} -\usage{ -single_posterior_calc2( - log_params_prop = c(), - input_data = list(), - obs_sero_data = NULL, - obs_case_data = NULL, - consts = list() -) -} -\arguments{ -\item{log_params_prop}{TBA} - -\item{input_data}{TBA} - -\item{obs_sero_data}{TBA} - -\item{obs_case_data}{TBA} - -\item{consts}{TBA} -} -\description{ -TBA -} -\details{ -TBA -}