R/estimate_R0.R

In kewok/Rgmf: R Toolkit Visualizing and Analyzing "Give Me Fever!" Epidemics

#'Estimate the basic reproduction rate
#'
#'For a given Epidemic object, return an object of type R0.sR using the estimate.R function from the package `R0'. In effect, estimate_R0 is a wrapper for feeding the json file generated by GMF to the estimate.R function of package `R0'. \cr
#'
#'@usage This function requires the R0 package. Estimate_R0(epidemic,  time_varying=FALSE,methods=c("ML","SB"))
#'@param epidemic An Epidemic object
#'@param trez The time resolution into which infection events will be sorted, as a fraction of the maximum time elapsed; default is 0.1 (i.e., in 10% time steps).
#'@param use_symptomatic_only A Boolean variable indicating whether only symptomatic infections should be considered; default is FALSE. Setting this to true can affect the calculation of the generation time.
#'@param trez The time resolution into which infection events will be sorted, as a fraction of the maximum time elapsed; default is 0.1 (i.e., in 10% time steps).
#'@details The distribution of the generation times (the time between a primary infection and a secondary infection) is estimated using the function `est.GT()' from the R0 package, which is fit by assuming that the serial interval (technically, the time between symptoms onset) is given by the recorded generation time.
#'@return An object of type R0.sR representing estimates of the basic reproduction number of the epidemic from the `estimate.R' function in the R0 package.
#'@seealso `estimate.R' from the package `R0'
#'@examples
#'my_epi <- Epidemic('example_gmf.json')
#'estimate_R0(my_epi)
#'@export

estimate_R0 <- function(epidemic,  trez=0.1, methods=c("ML","SB"), use_symptomatic_only=FALSE)
	{
	incidences <- calculate_incidence(epidemic, trez, use_symptomatic_only)
	dynamic_network <- contagion_network_evolution(epidemic, use_symptomatic_only)
	timings <- dynamic_network$sorted_infection_timings
	timings <- c(min(epidemic$hosts[[epidemic$patient_zero]]$host_switch_time), dynamic_network$sorted_infection_timings)
	genTimes <- c()

	# For each infection event, figure out how long it has been since the host was infected and since the host infected another host:
	# First, identify all the primary hosts:
	full_network <- reconstruct_epidemic_network(epidemic, use_symptomatic_only)
	primary_hosts <- unique(names(igraph::V(full_network)[igraph::get.edges(full_network, igraph::E(full_network))[,1]]))

	for (source_host in primary_hosts)
		{
		time_source_infected <- c() 
		# Identify the time-steps when that host got infected:
		for (this_edge in igraph::E(full_network))
			{
			if (names(igraph::V(full_network)[igraph::get.edges(full_network, this_edge)[,2]])==source_host)
				{
				time_source_infected <- c(time_source_infected, igraph::E(full_network)[this_edge]$infection_time)
				}			
			}
		if (source_host == epidemic$patient_zero)
			{
			time_source_infected <- c(time_source_infected, min(epidemic$hosts[[epidemic$patient_zero]]$host_switch_time))
			}
		# Make sure the timings of infections are ordered
		time_source_infected <- sort(unlist(time_source_infected))

		time_source_infects <- c()
		# for each source host, determine the times it infected another host
		endemic_cases <- epidemic$epidemic_network
		endemic_cases[[1]] <- NULL # Exclude the first case
		for (infection_event in endemic_cases)
			{
			if (source_host==infection_event$source_host)
				{
				time_source_infects <- c(time_source_infects, infection_event$infection_date)
				}
			}

		# The generation time is therefore the time between the source's most recent own infection event and the minimum time in which the source caused a subsequent infection 
		for (i in 1:length(time_source_infected))
			{
			onward_infections <- c()
			if (i < length(time_source_infected))
				{
				onward_infections <- intersect(which(time_source_infects >= time_source_infected[i]), which(time_source_infects < time_source_infected[i+1]))
				}
			else 
				{
				onward_infections <- which(time_source_infects > time_source_infected[i])
				}
			if (length(onward_infections)) # if there were any onward infections
				{
				genTimes <- c(genTimes, as.POSIXct(min(time_source_infects[onward_infections]), format="%Y-%m-%dT%H:%M:%OS") - as.POSIXct(time_source_infected[i], format="%Y-%m-%dT%H:%M:%OS"))
				}
			}
		}		

	# Fit the generation time distribution
	gentime <- R0::est.GT(serial.interval=genTimes)
	R0::estimate.R(incidences, GT=gentime, end=as.numeric(difftime(as.POSIXct(max(timings),format="%Y-%m-%dT%H:%M:%OS"), as.POSIXct(min(timings),format="%Y-%m-%dT%H:%M:%OS"))), methods=methods, pop.size=epidemic$N_hosts, S0=1)
	}