adding mean.equals.sd arg and updating docs

Author: gilesjohnr

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: IDSpatialStats
-Version: 0.3.4
-Date: 2018-06-02
+Version: 0.3.5
+Date: 2018-06-11
 Title: Estimate Global Clustering in Infectious Disease
 Author: Justin Lessler <[email protected]>, Henrik Salje <[email protected]>, John Giles <[email protected]>
 Maintainer: Justin Lessler <[email protected]>

R/examples/est_transdist_bootstrap_ci.R

@@ -1,6 +1,6 @@
 set.seed(1)
 
-# Normally distributed transmission kernel with mean and standard deviation = 100
+# Exponentially distributed transmission kernel with mean and standard deviation = 100
 dist.func <- alist(n=1, a=1/100, rexp(n, a)) 
 
 # Simulate epidemic
@@ -19,6 +19,7 @@ b <- est.transdist.bootstrap.ci(epi.data=a,
                                 max.sep=1e10,
                                 max.dist=1e10,
                                 n.transtree.reps=10,
+                                mean.equals.sd=TRUE,
                                 boot.iter=10,
                                 ci.low=0.025,
                                 ci.high=0.975,

R/examples/est_transdist_temporal.R

@@ -20,12 +20,15 @@ b <- est.transdist.temporal(epi.data=a,
                             t1=0,
                             max.sep=1e10,
                             max.dist=1e10,
-                            n.transtree.reps=3,
+                            n.transtree.reps=5,
+                            mean.equals.sd=TRUE,
                             parallel=FALSE)
 
-plot(b, pch=19, col='grey', ylim=c(min(b, na.rm=TRUE), max(b, na.rm=TRUE)), 
+plot(b[,1], pch=19, col='grey', ylim=c(min(b[,1], na.rm=TRUE), max(b[,1], na.rm=TRUE)), 
      xlab='Time step', ylab='Estimated mean of transmission kernel')
+abline(h=100, col='red', lty=2)
+axis(3, b[,2])
 
-low <- loess(b ~ as.vector(1:length(b)))
-low <- predict(low, newdata=data.frame(as.vector(1:length(b))))
+low <- loess(b[,1] ~ as.vector(1:length(b[,1])))
+low <- predict(low, newdata=data.frame(as.vector(1:length(b[,1]))))
 lines(low, lwd=3, col='blue')

R/examples/est_transdist_temporal_bootstrap_ci.R

@@ -1,5 +1,5 @@
 \donttest{
-  set.seed(1)
+  set.seed(123)
 
   # Exponentially distributed transmission kernel with mean and standard deviation = 100
   dist.func <- alist(n=1, a=1/100, rexp(n, a)) 
@@ -12,6 +12,8 @@
                     min.cases=30,
                     trans.kern.func=dist.func)
 
+  a <- a[sample(1:nrow(a), 100),] # subsample a to 100 observations
+  
   # Estimate change in mean transmission kernel over time with confidence intervals
   nc <- 4 # Run in parallel
 
@@ -21,15 +23,18 @@
                                            t1=0,
                                            max.sep=1e10,
                                            max.dist=1e10,
-                                           n.transtree.reps=3,
-                                           boot.iter=5,
+                                           n.transtree.reps=10,
+                                           mean.equals.sd=TRUE,
+                                           boot.iter=10,
                                            ci.low=0.025,
                                            ci.high=0.975,
                                            parallel=TRUE,
                                            n.cores=nc)
 
-  plot(b[,1], pch=19, col='grey', ylim=c(min(b, na.rm=TRUE), max(b, na.rm=TRUE)), 
+  plot(b[,1], pch=19, col='grey', ylim=c(min(b[,1:3], na.rm=TRUE), max(b[,1:3], na.rm=TRUE)), 
        xlab='Time step', ylab='Estimated mean of transmission kernel')
+  abline(h=100, col='red', lty=2)
+  axis(3, 1:nrow(b), b[,4])
 
   low <- loess(b[,1] ~ as.vector(1:nrow(b)), span=1)
   low <- predict(low, newdata=data.frame(as.vector(1:nrow(b))))

R/transdistfuncs.r

@@ -268,6 +268,7 @@ get.transdist.theta <-function(wal.teun.mat,
   # Clean matrix containing theta estimates
   a <- order(as.numeric(rownames(theta.mat))) # reorder by cases
   suppressWarnings(theta.mat <- matrix(as.integer(floor(theta.mat[a,a])), n, n))
+  #suppressWarnings(theta.mat <- matrix(as.integer(theta.mat[a,a])), n, n)
   diag(theta.mat) <- NA
 
   # Return theta matrix for unit testing
@@ -505,6 +506,7 @@ est.transdist <- function(
 ##' @param max.sep maximum number of time steps allowed between two cases (passed to the \code{get.transdist.theta} function)
 ##' @param max.dist maximum spatial distance between two cases considered in calculation
 ##' @param n.transtree.reps number of time to simulate transmission trees when estimating the weights of theta (passed to the \code{est.transdist.theta.weights} function, default = 10). Warning: higher values of this parameter cause significant increases in computation time.
+##' @param mean.equals.sd logical term indicating if the mean and standard deviation of the transmission kernel are expected to be equal (default = FALSE)
 ##' @param theta.weights use external matrix of theta weights. If NULL (default) the matrix of theta weights is automatically estimated by calling the \code{est.transdist.theta.weights} function
 ##' @param boot.iter the number of bootstrapped iterations to perform
 ##' @param ci.low low end of the confidence interval (default = 0.025)
@@ -532,6 +534,7 @@ est.transdist.bootstrap.ci <- function(
   max.sep, 		                         # Maximum time between cases considered in calculation
   max.dist, 	                         # Maximum distance considered in calculation
   n.transtree.reps=100, 			         # How many times to simulate transmission trees
+  mean.equals.sd=FALSE,                     # logical term indicating if the mean and sd of the transmission kernel are expected to be equal
   theta.weights=NULL,  			           # External matrix of theta weights
   boot.iter,
   ci.low=0.025,
@@ -551,6 +554,9 @@ est.transdist.bootstrap.ci <- function(
                      theta.weights=theta.weights,
                      silent=TRUE)
 
+  if(mean.equals.sd == TRUE) pt.est <- p$mu.est
+  if(mean.equals.sd == FALSE) pt.est <- p$bound.mu.est
+  
   n <- nrow(epi.data)
 
   if(parallel == FALSE) {
@@ -569,7 +575,10 @@ est.transdist.bootstrap.ci <- function(
                          n.transtree.reps=n.transtree.reps,
                          theta.weights=theta.weights,
                          silent=TRUE)
-      a$mu.est
+      
+      if(mean.equals.sd == TRUE) a <- a$mu.est
+      if(mean.equals.sd == FALSE) a <- a$bound.mu.est
+      a
     }
   }
 
@@ -596,14 +605,16 @@ est.transdist.bootstrap.ci <- function(
                          theta.weights=theta.weights,
                          silent=TRUE)
 
-      a$mu.est
+      if(mean.equals.sd == TRUE) a <- a$mu.est
+      if(mean.equals.sd == FALSE) a <- a$bound.mu.est
+      a
     }
     parallel::stopCluster(clust)
   }
 
   ci <- quantile(bs, probs=c(ci.low, ci.high))
 
-  return(list(mu.est=p$mu.est,
+  return(list(mu.est=pt.est,
               mu.ci.low=ci[1],
               mu.ci.high=ci[2]))
 }
@@ -620,11 +631,12 @@ est.transdist.bootstrap.ci <- function(
 ##' @param max.sep maximum number of time steps allowed between two cases (passed to the \code{get.transdist.theta} function)
 ##' @param max.dist maximum spatial distance between two cases considered in calculation
 ##' @param n.transtree.reps number of time to simulate transmission trees when estimating the weights of theta (passed to the \code{est.transdist.theta.weights} function, default = 10). Higher values of this parameter cause significant increases in computation time.
+##' @param mean.equals.sd logical term indicating if the mean and standard deviation of the transmission kernel are expected to be equal (default = FALSE)
 ##' @param theta.weights use external matrix of theta weights. If NULL (default) the matrix of theta weights is automatically estimated by calling the \code{est.transdist.theta.weights} function
 ##' @param parallel run time steps in parallel (default = FALSE)
 ##' @param n.cores number of cores to use when \code{parallel} = TRUE (default = NULL, which uses half the available cores)
 ##' 
-##' @return a vector containing the point estimate for mean transmission distance for each unique time step of the epidemic. 
+##' @return a numeric matrix containing the point estimate for mean transmission distance for each unique time step of the epidemic and the sample size $n$ used to make the estimate 
 ##' NAs are returned for time steps which contain fewer than three cases
 ##'
 ##' @author Justin Lessler, Henrik Salje, and John Giles
@@ -645,6 +657,7 @@ est.transdist.temporal <- function(
   max.sep, 		                         # Maximum time between cases considered in calculation
   max.dist, 	                         # Maximum distance considered in calculation
   n.transtree.reps=10, 			           # How many times to simulate transmission trees
+  mean.equals.sd=FALSE,
   theta.weights=NULL,  			           # External matrix of theta weights
   parallel=FALSE,
   n.cores=NULL
@@ -669,7 +682,6 @@ est.transdist.temporal <- function(
     out <- foreach::foreach(i=seq_along(unique.times), .combine='rbind') %do% {
 
       d <- epi.data[epi.data[,3] <= unique.times[i],]
-      n <- nrow(matrix(d, ncol=ncol(epi.data)))
 
       pt.est <- try(est.transdist(epi.data=d,
                                   gen.t.mean=gen.t.mean,
@@ -682,9 +694,10 @@ est.transdist.temporal <- function(
                     silent=TRUE)
       if(class(pt.est) == "try-error") {
         pt.est <- NA } else {
-          pt.est <- pt.est$mu.est
+          if(mean.equals.sd == TRUE) pt.est <- pt.est$mu.est
+          if(mean.equals.sd == FALSE) pt.est <- pt.est$bound.mu.est
         }
-      pt.est
+      x <- c(pt.est, nrow(d))
     }
   }
 
@@ -703,7 +716,6 @@ est.transdist.temporal <- function(
     out <- foreach::foreach(i=seq_along(unique.times), .combine='rbind') %dopar% {
 
       d <- epi.data[epi.data[,3] <= unique.times[i],]
-      n <- nrow(matrix(d, ncol=ncol(epi.data)))
 
       pt.est <- try(est.transdist(epi.data=d,
                                   gen.t.mean=gen.t.mean,
@@ -717,13 +729,16 @@ est.transdist.temporal <- function(
       if(class(pt.est) == "try-error") {
         pt.est <- NA 
       } else {
-        pt.est <- pt.est$mu.est
+           if(mean.equals.sd == TRUE) pt.est <- pt.est$mu.est
+           if(mean.equals.sd == FALSE) pt.est <- pt.est$bound.mu.est
       }
-      pt.est
+      x <- c(pt.est, nrow(d))
     }
     parallel::stopCluster(clust)
   }
-  return(as.vector(out))
+  rownames(out) <- NULL
+  colnames(out) <- c('pt.est', 'n')
+  return(out)
 }
 
 ##' Bootstrapped confidence intervals for the change in mean transmission distance over time
@@ -738,14 +753,15 @@ est.transdist.temporal <- function(
 ##' @param max.sep maximum number of time steps allowed between two cases (passed to the \code{get.transdist.theta} function)
 ##' @param max.dist maximum spatial distance between two cases considered in calculation
 ##' @param n.transtree.reps number of time to simulate transmission trees when estimating the weights of theta (passed to the \code{est.transdist.theta.weights} function, default = 10). Warning: higher values of this parameter cause significant increases in computation time.
+##' @param mean.equals.sd logical term indicating if the mean and standard deviation of the transmission kernel are expected to be equal (default = FALSE)
 ##' @param theta.weights use external matrix of theta weights. If NULL (default) the matrix of theta weights is automatically estimated by calling the \code{est.transdist.theta.weights} function
 ##' @param boot.iter the number of bootstrapped iterations to perform
 ##' @param ci.low low end of the confidence interval (default = 0.025)
 ##' @param ci.high high end of the confidence interval (default = 0.975)
 ##' @param parallel run bootstraps in parallel (default = FALSE)
 ##' @param n.cores number of cores to use when \code{parallel} = TRUE (default = NULL, which uses half the available cores)
 ##' 
-##' @return a three-column matrix containing the point estimate for mean transmission distance and low and high bootstrapped confidence intervals
+##' @return a four-column numeric matrix containing the point estimate for mean transmission distance, low and high bootstrapped confidence intervals, and the sample size up to each time step
 ##'
 ##' @author Justin Lessler, Henrik Salje, and John Giles
 ##' 
@@ -757,7 +773,6 @@ est.transdist.temporal <- function(
 ##' @example R/examples/est_transdist_temporal_bootstrap_ci.R
 ##'
 
-
 est.transdist.temporal.bootstrap.ci <- function(
   epi.data,				                     # Three column matrix: coordinates and time of case (must be in x, y, t, order)
   gen.t.mean,				                   # Mean of generation time
@@ -766,6 +781,7 @@ est.transdist.temporal.bootstrap.ci <- function(
   max.sep, 		                         # Maximum time between cases considered in calculation
   max.dist, 	                         # Maximum distance considered in calculation
   n.transtree.reps=100, 			         # How many times to simulate transmission trees
+  mean.equals.sd=FALSE,
   theta.weights=NULL,  			           # External matrix of theta weights
   boot.iter,
   ci.low=0.025,
@@ -806,7 +822,8 @@ est.transdist.temporal.bootstrap.ci <- function(
                     silent=TRUE)
       if(class(pt.est) == "try-error") {
         pt.est <- NA } else {
-          pt.est <- pt.est$mu.est
+             if(mean.equals.sd == TRUE) pt.est <- pt.est$mu.est
+             if(mean.equals.sd == FALSE) pt.est <- pt.est$bound.mu.est
         }
 
       bs <- rep(NA, boot.iter)
@@ -825,12 +842,13 @@ est.transdist.temporal.bootstrap.ci <- function(
                  silent=TRUE)
         if(class(a) == "try-error") {
           bs[j] <- NA} else {
-            bs[j] <- a$mu.est
+               if(mean.equals.sd == TRUE) bs[j] <- a$mu.est
+               if(mean.equals.sd == FALSE) bs[j] <- a$bound.mu.est
           }
       }
 
       suppressWarnings(bs.probs <- quantile(as.numeric(bs), probs=c(ci.low, ci.high), na.rm=TRUE))
-      x <- c(pt.est, bs.probs[1], bs.probs[2])
+      x <- c(pt.est, bs.probs[1], bs.probs[2], n)
     }
   }
 
@@ -863,7 +881,8 @@ est.transdist.temporal.bootstrap.ci <- function(
       if(class(pt.est) == "try-error") {
         pt.est <- NA 
         } else {
-          pt.est <- pt.est$mu.est
+             if(mean.equals.sd == TRUE) pt.est <- pt.est$mu.est
+             if(mean.equals.sd == FALSE) pt.est <- pt.est$bound.mu.est
         }
 
       bs <- rep(NA, boot.iter)
@@ -883,18 +902,19 @@ est.transdist.temporal.bootstrap.ci <- function(
         if(class(a) == "try-error") {
           bs[j] <- NA
           } else {
-            bs[j] <- a$mu.est
+               if(mean.equals.sd == TRUE) bs[j] <- a$mu.est
+               if(mean.equals.sd == FALSE) bs[j] <- a$bound.mu.est
           }
       }
 
       suppressWarnings(bs.probs <- quantile(as.numeric(bs), probs=c(ci.low, ci.high), na.rm=TRUE))
-      x <- c(pt.est, bs.probs[1], bs.probs[2])
+      x <- c(pt.est, bs.probs[1], bs.probs[2], n)
     }
     parallel::stopCluster(clust)
   }
 
   row.names(out) <- NULL
-  colnames(out) <- c('mu.est', 'mu.ci.low', 'mu.ci.high')
+  colnames(out) <- c('mu.est', 'mu.ci.low', 'mu.ci.high', 'n')
   return(out)
 }