R/estimate_Re.R

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

#'Calculate an approximation of the effective reproduction rate
#'
#'For a given Epidemic object, return an object of type R0.R from the `R0' package that describes the effective reproduction number of the epidemic simulated over the course of the epidemic spreading. In effect, estimate_Re is a wrapper for feeding the json file generated by GMF to the est.R0.TD function of package `R0'. \cr
#'
#'@usage estimate_Re(epidemic)
#'@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.
#'@return An object of type R0.R that describes the effective reproductive number. Note that the time interval used is determined according to infection events.
#'@examples
#'my_epi <- Epidemic('example_gmf.json')
#'
#'estimate_Re(my_epi)
#'
#'@export

estimate_Re <- function(epidemic, trez=0.1, use_symptomatic_only=FALSE)
	{
	incidences <- calculate_incidence(epidemic, trez, use_symptomatic_only)
	incidences <- c(1,incidences)
	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::est.R0.TD(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"))), nsim=1e6,n.t0=1)
	}