R/Rj.R

In wenlab501/EpiTrans: R package for estimation of epidemic transmissibility in time and space

#' Calculate reproductive numbers
#'
#' Calculate individual reproductive numbers given data (n-length vectors : t, x, y) and prespecified likelihood functions (lpdf_GI, lpdf_SP)
#'
#' @param t Vector of time of data points.
#' @param x Vector of x coordinates of data points.
#' @param y Vector of y coordinates of data points.
#' @param adjSP Logical; if TRUE, reproductive numbers will be calculated with spatial adjustment; if FALSE, static plotting non-interactive thematic map reproductive numbers will be calculated only considering temporal relation.
#' @param GI.pdf Probability distribution of generation interval; it can be generated by function 'lpdf_GI'.
#' @param SW.pdf Probability distribution (exponential distribution) of tranmission distance; it can be generated by function 'lpdf_SW'.
#' @param unit_coord The unit of the coordination; if 'degree' then calculate distance througth 'distLongLat'; if 'meter' then calculate distance directly.
#' @return
#'
#' @examples
#' data("EpiTrans")
#'
#' res_adj = Rj(t = dengue$date, x = dengue$long, y = dengue$lat,GI.pdf = lpdf_GI(),SW.pdf = lpdf_SW(), adjSP = TRUE)
#' res = Rj(t = dengue$date, GI.pdf = lpdf_GI(), adjSP = FALSE)
#' @export

Rj <- function(t, x, y, adjSP=TRUE, GI.pdf=lpdf_GI(), SW.pdf=lpdf_SW(), unit_coord=c("degree","meter")){
  N <- length(t)
  data <- data.frame(t = as.numeric(t))
  if(adjSP){
    if(missing(x)) stop("missing argument x")
    if(missing(x)) stop("missing argument y")
    if(N!=length(x)) stop("x have different length from t")
    if(N!=length(y)) stop("y have different length from t")
    data$x <- as.numeric(x)
    data$y <- as.numeric(y)
    unit_coord <- unit_coord[1]
  }
  data <- data[order(data$t),]

  #calculate Rt consider time
  dist.time <- as.vector(dist(t))
  lik <- GI.pdf(dist.time)

  #calculate Rt consider space
  if(adjSP){
    if(unit_coord=="degree") dist.space <- distLongLat(data$x,data$y)
    else if(unit_coord=="meter") dist.space <-as.vector(dist(cbind(data$x,data$y)))
    else stop("set space unit to degree or meter")
    dist.space[dist.space==0] <- 1
    likSP <- SW.pdf(dist.space)
    likSP[lik==0] <- 0
    lik <- lik + likSP
  }

  lik <- exp(lik)

  mat <- matrix(0,N,N)
  mat[lower.tri(mat)] <- lik
  mat <- t(mat)

  matColSum <- 1 / colSums(mat,na.rm = T)
  Pij <- t(apply(mat,1,function(x) x*matColSum))
  Rj <- rowSums(Pij,na.rm = T)
  return(list(Rj=Rj, Pij=Pij))
}