Functional changes for graphical hypothesis testing and clustering range estimation

.Rbuildignore

@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,7 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+src/*.o
+src/*.so
+src/*.dll

DESCRIPTION

@@ -1,16 +1,25 @@
 Package: IDSpatialStats
-Version: 0.3.11
-Date: 2019-11-13
+Version: 1.0.0
+Date: 2021-07-21
 Title: Estimate Global Clustering in Infectious Disease
-Author: John Giles <[email protected]>, Henrik Salje <[email protected]>, Justin Lessler <[email protected]> 
-Maintainer: Justin Lessler <[email protected]>
+Authors@R: c(person(given = "Justin Lessler", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-9741-8109", "Co-creator & maintainer"),
+                     email = "[email protected]"),
+              person(given = "Henrik Salje", role = "aut", comment = c(ORCID = "0000-0003-3626-4254", "Co-creator"),
+                     email = "[email protected]"),
+              person(given = "John Giles", role = "aut", comment = c(ORCID = "0000-0002-0954-4093", "Author & maintainer"),
+                     email = "[email protected]"),
+              person(given = "Timothy M Pollington", role = "aut", comment = c(ORCID = "0000-0002-9688-5960", "Author"),
+                     email = "[email protected]"))
 License: GPL (>=2)
 Description: Implements various novel and standard
     clustering statistics and other analyses useful for understanding the
     spread of infectious disease.
 RoxygenNote: 7.1.1
 Encoding: UTF-8
 LazyData: true
-Imports: igraph, spatstat.core, spatstat.geom
+Imports: igraph, spatstat.core, spatstat.geom, coxed, GET, scales
 Depends: doParallel, foreach, parallel, R (>= 2.10)
-Suggests: knitr, testthat
+Suggests: knitr, testthat, rmarkdown
+VignetteBuilder: knitr
+URL: https://github.com/HopkinsIDD/IDSpatialStats
+BugReports: https://github.com/HopkinsIDD/IDSpatialStats/issues

IDSpatialStats.Rproj

@@ -0,0 +1,17 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 1
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source

NAMESPACE

@@ -1,17 +1,20 @@
+import("coxed")
 import("parallel")
 import("doParallel")
 import("foreach")
 import("stats")
 import("graphics")
-importFrom("igraph", "graph.data.frame") # for get.transdist.theta
-importFrom("igraph", "shortest.paths") # for get.transdist.theta
-importFrom("igraph", "E") # for get.transdist.theta
-importFrom("spatstat.geom", "bounding.box.xy") # for est.transdist
-importFrom("spatstat.geom", "as.ppp") # for est.transdist
-importFrom("spatstat.geom", "ppp") # for est.transdist
-importFrom("spatstat.geom", "crossdist") # for est.transdist
-importFrom("spatstat.core", "Kcross") # for get.cross.K
-importFrom("spatstat.core", "pcf") # for get.cross.PCF
+import("scales")
+import("GET")
+importFrom("igraph", "graph.data.frame") # for get.transdist.theta()
+importFrom("igraph", "shortest.paths") # for get.transdist.theta()
+importFrom("igraph", "E") # for get.transdist.theta()
+importFrom("spatstat.geom", "bounding.box.xy") # for est.transdist()
+importFrom("spatstat.geom", "as.ppp") # for est.transdist()
+importFrom("spatstat.geom", "ppp") # for est.transdist()
+importFrom("spatstat.geom", "crossdist") # for est.transdist()
+importFrom("spatstat.core", "Kcross") # for get.cross.K()
+importFrom("spatstat.core", "pcf") # for get.cross.PCF()
 export(get.pi)
 export(get.pi.ci)
 export(get.pi.bootstrap)
@@ -22,6 +25,8 @@ export(get.theta.bootstrap)
 export(get.theta.permute)
 export(get.tau)
 export(get.tau.bootstrap)
+export(get.tau.GET)
+export(get.tau.D.param.est)
 export(get.tau.permute)
 export(get.pi.typed)
 export(get.pi.typed.bootstrap)
@@ -33,6 +38,7 @@ export(get.tau.typed)
 export(get.tau.ci)
 export(get.tau.typed.bootstrap)
 export(get.tau.typed.permute)
+S3method(plot, tau)
 export(sim.epidemic)
 export(sim.plot)
 export(est.wt.matrix)

NEWS.md

@@ -0,0 +1,48 @@
+# News on IDSpatialStats 1.0.0 release
+
+## Changes (top of list are most important)
+These changes mostly concern the tau statistic functions and have resulted in a major change revision from 0.3.9 to 1.0.0: these are big changes and we may have introduced bugs so please send us a `reprex()` example. We also note where function outputs are likely to have changed. 
+
+# Specific changes
+Percentile confidence intervals (CIs) replaced with BCa (bias-corrected and accelerated) CIs
+* `get.pi.ci()`, `get.theta.ci()`, `get.tau.ci()`: `quantile` method replaced with `coxed::bca`. This will change CI results.
+
+* `quantile` method also updated to `coxed::bca` where possible in `inst/tests/`. At times the `coxed::bca()` method gives slightly different test results if it is applied to asymmetric distributions.
+
+* New function `get.tau.GET()` performs graphical hypothesis tests with the tau statistic using the `GET` package, while `get.tau.D.param.est()` estimates the range of spatiotemporal clustering.
+
+* `plot.tau()` is a new method that can produce three types of tau(y-axis)-distance(x-axis) plots [[1: see Graphical abstract](https://doi.org/10.1016/j.spasta.2020.100438 "Developments in statistical inference when assessing spatiotemporal disease clustering with the tau statistic")]:
+  * Diagnostic to indicate the structure or magnitude of spatiotemporal clustering. Requires `tau` object; `tauCI` object optional to draw pointwise CIs. In this version we use piecewise error bars rather than continuous envelope lines. This reminds the user that this graph should not be used as a graphical hypothesis test for the whole distance range observed. It is only suitable for the purpose of a graphical hypothesis test for a specific distance band if that band is decided prior to graph creation.
+  * Graphical hypothesis tests to assess the evidence against the Null hypothesis (no spatiotemporal clustering nor inhibition). Requires `tau` and `tauGET` objects.
+  * Estimation of the clustering range (the distribution of the places on the horizontal tau=1 line, where decreasing bootstrap simulations first intercept). Requires `tau` and `tauparamest` objects; prior to this `get.tau.D.param.est()` requires a `taubstrap` object.
+
+* New S3 classes have been added to the return objects from the following functions. The purpose of the new classes is to encourage the use of functions in an ordered and principled way, in keeping with good practices of statistical inference [[1](https://doi.org/10.1016/j.spasta.2020.100438 "Developments in statistical inference when assessing spatiotemporal disease clustering with the tau statistic")]. It also means that the objects inputted as function arguments are of a known format:
+  * `get.tau()` returns a `tau` class
+  * `get.tau.ci()` returns a `tauCI` class
+  * `get.tau.GET()` returns a `tauGET` class
+  * `get.tau.bootstrap()` returns a `taubstrap` class
+  * `get.tau.D.param.est()` returns a `tauparamest` class. Requires a `taubstrap` object. Also requires a `tauGET` class to ensure the user has performed a graphical hypothesis test first, before considering parameter estimation.
+
+* CITATION file added
+* README.md formatting updated
+* `get.tau$tau` renamed to `get.tau$tau.pt.est`
+* Previously deprecated tests (already commented out in `inst/tests/`) have now been removed.
+* NEWS.md file added, but as you're reading this you probably knew that already :wink:
+
+# Generic changes
+* distance units can be defined on `r` and `r.low` and will be automatically feature in the x-axis label of `plot.tau()`
+* help files added/updated for `get.tau()`, `get.tau.ci()`, `get.tau.bootstrap()`, `get.tau.GET()`, `get.tau.d.param.est()` & `plot.tau()`
+* example files added/updated for `get_tau.R`, `get_tau_bootstrap.R`, `get_tau_ci.R`, `get_tau_GET.R`, `get_D_param_est.R` and `plot_tau.R`.
+
+## Bug fixes (top of list are most important)
+`get_tau.R`: using `geno.tau.R02$tau.pt.est` now allows the object to be accessed and the example run.
+
+# Release contributors
+Timothy M Pollington would like to thank the co-authors of the paper that informed this update[[1](https://doi.org/10.1016/j.spasta.2020.100438 "Developments in statistical inference when assessing spatiotemporal disease clustering with the tau statistic")] and the *essential* contribution of Peter J. Diggle (Lancaster) who advised on this principled inferential approach.
+
+# Next changes
+* The *Modified Marked Point Spatial Bootstrap* [[1](https://doi.org/10.1016/j.spasta.2020.100438 "Developments in statistical inference when assessing spatiotemporal disease clustering with the tau statistic")] has not yet been
+applied. In [[2](https://doi.org/10.5281/zenodo.2552850 "t-pollington/tau-statistic-speedup: First release of tau statistic speedup")] it was applied to the tau odds estimator however for consistency we have decided to delay its implementation so that we can apply it also to the tau prevalence estimator. So please note that `get.tau.bootstrap()` and `get.tau.D.param.est()` values are still likely to change.
+* `tau.GET()` examples increased to 2,500 permutations once speedups introduced.
+* Enable `plot.tau()` to accept `tauCI` objects alone, without need for `tau()` object.
+* Changes to the un-typed tau functions also applied to the typed tau functions.

R/examples/get_pi_ci.R

@@ -1,26 +1,27 @@
-\donttest{
+\donttest{# Simulate cases (type = 2, Normally-distributed points) and 
+# simulated non-cases (type = 1, Uniformally-distributed)
+X <- cbind(1, runif(100,-100,100), runif(100,-100,100))
+X <- rbind(X, cbind(2,rnorm(100,0,20), rnorm(100,0,20)))
+colnames(X) <- c("type","x","y")
 
-#compare normally distributed with uniform points
-x<-cbind(1,runif(100,-100,100), runif(100,-100,100))
-x<-rbind(x, cbind(2,rnorm(100,0,20), rnorm(100,0,20)))
-colnames(x) <- c("type","x","y")
-
-fun<-function(a,b) {
-    if(a[1]!=2) return(3)
-    if (b[1]==2) return(1)
-    return(2)
+fun <- function(a,b) {
+    # possible 'ab' pair types {'11'; '12'; '21'; '22'}
+    if(a[1]!=2) return(3) # it's {'11' or '12'} so ignore
+    if(b[1]==2) return(1) # it's '22' so count as a case-case pair in numerator and denominator
+    # else it's a '21' ie a case-non-case pair
+    return(2) # so count in denominator
 }
 
-r.max<-seq(10,100,10)
-r.min<-seq(0,90,10)
-r.mid <- (r.max+r.min)/2
-
+# define distance band set
+r.max <- seq(10,100,10)
+r.min <- seq(0,90,10)
+r.mid <- (r.max + r.min)/2
 
-pi<-get.pi(x,fun,r=r.max,r.low=r.min)
-pi.ci<-get.pi.ci(x,fun,r=r.max,r.low=r.min,boot.iter=100)
+# compute the pi point estimate and its 95% BCa CI
+pi <- get.pi(X, fun, r=r.max, r.low = r.min)
+pi.ci <- get.pi.ci(X, fun, r = r.max, r.low = r.min, boot.iter = 100)
 
+# plot the pi point estimate with its CI, at the midpoints of each distance band
 plot(r.mid, pi$pi, type="l")
-lines(r.mid, pi.ci[,2] , lty=2)
-lines(r.mid, pi.ci[,3] , lty=2)
-
-}
+lines(r.mid, pi.ci$ci.low, lty=2)
+lines(r.mid, pi.ci$ci.high, lty=2)}

R/examples/get_tau.R

@@ -7,39 +7,24 @@ data(DengueSimRepresentative)
 r.max<-seq(20,1000,20)
 r.min<-seq(0,980,20)
 r.mid<-(r.max+r.min)/2
-
 sero.type.func<-function(a,b,tlimit=20){
-     if(a[5]==b[5]&(abs(a[3]-b[3])<=tlimit)){rc=1}
-     else{rc=2}
-     return(rc)
+        if(a[5]==b[5]&(abs(a[3]-b[3])<=tlimit)){rc=1}
+        else{rc=2}
+        return(rc)
 }
-
 geno.type.func<-function(a,b,tlimit=20){
-     if(a[4]==b[4]&(abs(a[3]-b[3])<=tlimit)){rc=1}
-     else{rc=2}
-     return(rc)
-}
-
-sero.type.rep.func<-function(a,b,tlimit=20){
-     if(a[5]==1&b[5]==1&(abs(a[3]-b[3])<=tlimit)){rc=1}
-     else{if(a[5]==1&b[5]==-999){rc=2}else{rc=3}}
-     return(rc)
+        if(a[4]==b[4]&(abs(a[3]-b[3])<=tlimit)){rc=1}
+        else{rc=2}
+        return(rc)
 }
 
-sero.tau.R01 <- get.tau(DengueSimR01, sero.type.func, r=r.max, r.low=r.min, 
-                      comparison.type="independent")
-geno.tau.R01 <- get.tau(DengueSimR01, geno.type.func, r=r.max, r.low=r.min, 
-                      comparison.type="independent")
-
-sero.tau.R02 <- get.tau(DengueSimR02, sero.type.func, r=r.max, r.low=r.min, 
-                      comparison.type="independent")
-geno.tau.R02 <- get.tau(DengueSimR02, geno.type.func, r=r.max, r.low=r.min, 
-                      comparison.type="independent")
-
-sero.tau.representative <- get.tau(DengueSimRepresentative, sero.type.rep.func, 
-                                   r=r.max, r.low=r.min, comparison.type="representative")
-
 ## R0 of 1
+data(DengueSimulationR01)
+sero.tau.R01 <- get.tau(DengueSimR01, sero.type.func, r=r.max, r.low=r.min,
+                        comparison.type="independent")
+geno.tau.R01 <- get.tau(DengueSimR01, geno.type.func, r=r.max, r.low=r.min,
+                        comparison.type="independent")
+
 plot(r.mid,sero.tau.R01$tau,ylim=c(0.3,max(geno.tau.R01$tau)),log="y",
      cex.axis=1.25,col=rgb(t(col2rgb("blue")/255),alpha=0.6),
      xlab="Distance (m)",ylab="Tau",cex.main=0.9,lwd=2,type="l",las=1,cex.axis=0.75)
@@ -57,7 +42,13 @@ legend("topright",
        lty=c(1,1,2,1),bty="n")
 
 ## R0 of 2
-plot(r.mid,sero.tau.R02$tau,ylim=c(0.3,max(geno.tau.R02)),log="y",
+data(DengueSimulationR02)
+sero.tau.R02 <- get.tau(DengueSimR02, sero.type.func, r=r.max, r.low=r.min,
+                        comparison.type="independent")
+geno.tau.R02 <- get.tau(DengueSimR02, geno.type.func, r=r.max, r.low=r.min,
+                        comparison.type="independent")
+
+plot(r.mid,sero.tau.R02$tau,ylim=c(0.3,max(geno.tau.R02$tau.pt.est)),log="y",
      cex.axis=1.25,col=rgb(t(col2rgb("blue")/255),alpha=0.6),
      xlab="Distance (m)",ylab="Tau",cex.main=0.9,lwd=2,type="l",las=1,cex.axis=0.75)
 abline(h=1,lty=2)
@@ -69,4 +60,22 @@ legend("topright",
                 "Maximum transmission distance"),
        lwd=1,col=c("dark green","blue","black"),lty=1,bty="n")
 
+## Obtaining a diagnostic plot using plot.tau() with pointwise CIs
+data(DengueSimRepresentative)
+sero.type.rep.func<-function(a,b,tlimit=20){
+     if(a[5]==1&b[5]==1&(abs(a[3]-b[3])<=tlimit)){rc=1}
+     else{if(a[5]==1&b[5]==-999){rc=2}else{rc=3}}
+     return(rc)
+}
+
+# get point estimate
+Dengue.tau = get.tau(DengueSimRepresentative, sero.type.rep.func, r.max, r.min,
+                     "representative", data.frame = TRUE)
+
+# get 95% BCa CI
+CIs = get.tau.ci(DengueSimRepresentative, sero.type.rep.func, r.max, r.min, 25,
+                 "representative", ci.level = 0.95, data.frame = TRUE)
+
+#plot point estimate with CI
+plot.tau(x = Dengue.tau, r.mid = TRUE, ptwise.CI = CIs)
 }

R/examples/get_tau_D_param_est.R

@@ -0,0 +1,38 @@
+\donttest{
+# Load data
+r.max<-seq(20,1000,20)
+r.min<-seq(0,980,20)
+r.mid<-(r.max+r.min)/2
+sero.type.func<-function(a,b,tlimit=20){
+        if(a[5]==b[5]&(abs(a[3]-b[3])<=tlimit)){rc=1}
+        else{rc=2}
+        return(rc)
+}
+
+data(DengueSimRepresentative)
+sero.type.rep.func<-function(a,b,tlimit=20){
+     if(a[5]==1&b[5]==1&(abs(a[3]-b[3])<=tlimit)){rc=1}
+     else{if(a[5]==1&b[5]==-999){rc=2}else{rc=3}}
+     return(rc)
+}
+
+# get point estimate
+Dengue.tau = get.tau(DengueSimRepresentative, sero.type.rep.func, r.max, r.min, "representative",
+                     data.frame = TRUE)
+
+# perform graphical hypothesis test using a global envelope test
+Dengue.GET = get.tau.GET(DengueSimRepresentative, sero.type.rep.func, r.max,r.min, 
+                         permutations = 50, "representative")
+
+# plot point estimate with global envelope and simulation of the null distribution
+plot.tau(x = Dengue.tau, r.mid = TRUE, GET.res = Dengue.GET)
+
+# if the graphical hypothesis test and p-value interval suggests evidence against H_0, 
+# and the graph suggests clustering, the range of this can be estimated
+tausim = get.tau.bootstrap(DengueSimRepresentative, sero.type.rep.func, r.max, r.min, 100, 
+                           "representative", data.frame = FALSE)
+Dengue.dparam = get.tau.D.param.est(r = r.max, tausim = tausim, Dengue.GET)
+median(Dengue.dparam$envelope) # median estimate for the clustering endpoint
+Dengue.dparam$envelope # 95% BCa CI
+plot.tau(Dengue.tau, tausim = tausim, d.param.est = Dengue.dparam)
+}

R/examples/get_tau_GET.R

@@ -0,0 +1,29 @@
+\donttest{
+# Load data
+r.max<-seq(20,1000,20)
+r.min<-seq(0,980,20)
+r.mid<-(r.max+r.min)/2
+sero.type.func<-function(a,b,tlimit=20){
+        if(a[5]==b[5]&(abs(a[3]-b[3])<=tlimit)){rc=1}
+        else{rc=2}
+        return(rc)
+}
+
+data(DengueSimRepresentative)
+sero.type.rep.func<-function(a,b,tlimit=20){
+     if(a[5]==1&b[5]==1&(abs(a[3]-b[3])<=tlimit)){rc=1}
+     else{if(a[5]==1&b[5]==-999){rc=2}else{rc=3}}
+     return(rc)
+}
+
+# get point estimate
+Dengue.tau = get.tau(DengueSimRepresentative, sero.type.rep.func, r.max, r.min, "representative", 
+                     data.frame = TRUE)
+
+# perform graphical hypothesis test using a global envelope test
+Dengue.GET = get.tau.GET(DengueSimRepresentative, sero.type.rep.func, r.max,r.min, 
+                         permutations = 50, "representative")
+
+#plot point estimate with global envelope and simulation of the null distribution
+plot.tau(x = Dengue.tau, r.mid = TRUE, GET.res = Dengue.GET)
+}