diff --git a/.Rbuildignore b/.Rbuildignore
new file mode 100644
index 0000000..91114bf
--- /dev/null
+++ b/.Rbuildignore
@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..56843bc
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,7 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+src/*.o
+src/*.so
+src/*.dll
diff --git a/DESCRIPTION b/DESCRIPTION
index 5b83bf9..f5e881c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,9 +1,15 @@
 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 <jrgiles@uw.edu>, Henrik Salje <hs743@cam.ac.uk>, Justin Lessler <jlessler@unc.edu> 
-Maintainer: Justin Lessler <jlessler@unc.edu>
+Authors@R: c(person(given = "Justin Lessler", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-9741-8109", "Co-creator & maintainer"),
+                     email = "justin@jhu.edu"),
+              person(given = "Henrik Salje", role = "aut", comment = c(ORCID = "0000-0003-3626-4254", "Co-creator"),
+                     email = "hs743@cam.ac.uk"),
+              person(given = "John Giles", role = "aut", comment = c(ORCID = "0000-0002-0954-4093", "Author & maintainer"),
+                     email = "hs743@cam.ac.uk"),
+              person(given = "Timothy M Pollington", role = "aut", comment = c(ORCID = "0000-0002-9688-5960", "Author"),
+                     email = "timothy.pollington@gmail.com"))
 License: GPL (>=2)
 Description: Implements various novel and standard
     clustering statistics and other analyses useful for understanding the
@@ -11,6 +17,9 @@ Description: Implements various novel and standard
 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
diff --git a/IDSpatialStats.Rproj b/IDSpatialStats.Rproj
new file mode 100644
index 0000000..c52d81d
--- /dev/null
+++ b/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
diff --git a/NAMESPACE b/NAMESPACE
index c50fcf0..bc51bf8 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/NEWS.md b/NEWS.md
new file mode 100644
index 0000000..b84ab4a
--- /dev/null
+++ b/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.
\ No newline at end of file
diff --git a/R/examples/get_pi_ci.R b/R/examples/get_pi_ci.R
index 59c2e0d..cf1aea9 100644
--- a/R/examples/get_pi_ci.R
+++ b/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)
-
-}
\ No newline at end of file
+lines(r.mid, pi.ci$ci.low, lty=2)
+lines(r.mid, pi.ci$ci.high, lty=2)}
\ No newline at end of file
diff --git a/R/examples/get_tau.R b/R/examples/get_tau.R
index ce99e51..ed842d5 100644
--- a/R/examples/get_tau.R
+++ b/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)
 }
diff --git a/R/examples/get_tau_D_param_est.R b/R/examples/get_tau_D_param_est.R
new file mode 100644
index 0000000..504958d
--- /dev/null
+++ b/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)
+}
\ No newline at end of file
diff --git a/R/examples/get_tau_GET.R b/R/examples/get_tau_GET.R
new file mode 100644
index 0000000..5384683
--- /dev/null
+++ b/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)
+}
\ No newline at end of file
diff --git a/R/examples/get_tau_bootstrap.R b/R/examples/get_tau_bootstrap.R
index d0eaab2..631a295 100644
--- a/R/examples/get_tau_bootstrap.R
+++ b/R/examples/get_tau_bootstrap.R
@@ -1,6 +1,6 @@
 \donttest{
 
-#compare normally distributed with uniform points
+# 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")
@@ -15,14 +15,17 @@ r.max<-seq(10,100,10)
 r.min<-seq(0,90,10)
 r.mid <- (r.max+r.min)/2
 
-tau<-get.tau(x,fun,r=r.max,r.low=r.min)
-tau.boot<-get.tau.bootstrap(x,fun,r=r.max,r.low=r.min,boot.iter=50)
+tau<-get.tau(x,fun,r=r.max,r.low=r.min,"representative", data.frame = TRUE)
+tau.ci = get.tau.ci(x, fun, r.max, r.min, 50, "representative", 0.95, data.frame = TRUE)
 
-tau.ci<-apply(tau.boot[,-(1:2)],1,quantile,probs=c(0.25,0.75))
+## plot.tau() method
+plot.tau(tau, r.mid = TRUE, ptwise.CI = tau.ci)
 
-plot(r.mid, tau$tau ,ylim=c(min(tau.ci),max(tau.ci)), type="l", log="y")
+## previous plot() method using connected lines to join the top and bottoms of the pointwise CIs. 
+#This may lead the user to perform graphical hypothesis testing using this plot without considering 
+#the specific distance band of interest before plotting.
+plot(r.mid, tau$tau.pt.est ,ylim=c(min(tau.ci$ci.low),max(tau.ci$ci.high)), type="l", log="y")
 lines(c(0,100),c(1,1), lty=3, col="grey")
-lines(r.mid, tau.ci[1,] , lty=2)
-lines(r.mid, tau.ci[2,] , lty=2)
-
+lines(r.mid, tau.ci$ci.low, lty=2)
+lines(r.mid, tau.ci$ci.high, lty=2)
 }
\ No newline at end of file
diff --git a/R/examples/get_tau_ci.R b/R/examples/get_tau_ci.R
index 8400219..846ab52 100644
--- a/R/examples/get_tau_ci.R
+++ b/R/examples/get_tau_ci.R
@@ -15,12 +15,17 @@ r.max<-seq(10,100,10)
 r.min<-seq(0,90,10)
 r.mid <- (r.max+r.min)/2
 
-tau <- get.tau.ci(x,fun,r=r.max,r.low=r.min,boot.iter=50)
+tau.CI <- get.tau.ci(x,fun,r=r.max,r.low=r.min,boot.iter=50, comparison.type = "representative")
 
-plot(r.mid, tau$pt.est, ylim=c(1/max(tau[,3:5]), max(tau[,3:5])), type="l", log="y",
-     xlab="Distance", ylab="Tau")
-lines(r.mid, tau$ci.low , lty=2)
-lines(r.mid, tau$ci.high, lty=2)
-lines(c(0,100),c(1,1), lty=3, col="grey")
+## plot.tau() method
+tau = get.tau(x,fun,r=r.max,r.low=r.min, comparison.type = "representative")
+plot.tau(x = tau, ptwise.CI = tau.CI)
 
+## previous plot() method
+plot(r.mid, tau.CI$pt.est, ylim=c(min(tau.CI$pt.est,tau.CI$ci.low), 
+                                  max(tau.CI$pt.est,tau.CI$ci.high)), type="l", xlab="Distance", 
+     ylab="Tau")
+lines(r.mid, tau.CI$ci.low , lty=2)
+lines(r.mid, tau.CI$ci.high, lty=2)
+lines(c(0,100),c(1,1), lty=3, col="grey")
 }
\ No newline at end of file
diff --git a/R/examples/get_tau_permute.R b/R/examples/get_tau_permute.R
index a1bfa42..b2d1ccc 100644
--- a/R/examples/get_tau_permute.R
+++ b/R/examples/get_tau_permute.R
@@ -20,9 +20,11 @@ tau.null<-get.tau.permute(x,fun,r=r.max,r.low=r.min,permutations=50,comparison.t
 
 null.ci<-apply(tau.null[,-(1:2)],1,quantile,probs=c(0.25,0.75))
 
+# note these plots are only for illustrative purposes to show how get.tau.permute() can generate 
+# the null distribution. These should not be used for graphical hypothesis testing nor parameter
+# estimation of the clustering endpoint.
 plot(r.mid, tau$tau, ylim=c(1/max(tau$tau),max(tau$tau)), type="l", log="y")
 lines(c(0,100),c(1,1), lty=3, col="grey")
 lines(r.mid, null.ci[1,] , lty=2)
 lines(r.mid, null.ci[2,] , lty=2)
-
 }
\ No newline at end of file
diff --git a/R/spatialfuncs.r b/R/spatialfuncs.r
index 2e791b6..349d688 100644
--- a/R/spatialfuncs.r
+++ b/R/spatialfuncs.r
@@ -1,20 +1,28 @@
+applyBCa <- function(boots, ci.level){
+  boots = boots[!is.na(boots)]
+  CI = coxed::bca(boots, conf.level = ci.level)
+  return(CI)
+}
+
 ##' Generalized version of \code{get.pi}
 ##'
 ##' Generalized version of the \code{get.pi} function that takes in an arbitrary function and
-##' returns the probability that a point within a particular range of a point of interest shares the relationship
+##' returns the probability that a point within a particular range of a point of interest shares 
+##' the relationship
 ##' specified by the passed in function with that point.
 ##'
-##' @param posmat a matrix with columns x, y and any other named columns
+##' @param posmat a matrix with columns x, y and any other named
 ##'    columns needed by \code{fun}
 ##' @param fun a function that takes in two rows of \code{posmat} and returns:
 ##' \enumerate{
 ##'      \item for pairs included in the numerator and denominator
 ##'      \item for pairs that should only be included in the denominator
 ##'      \item for pairs that should be ignored all together}
-##' Note that names from \code{posmat} are not preserved in calls to \code{fun}, so the columns of the matrix should be
+##' Note that names from \code{posmat} are not preserved in calls to \code{fun}, so the columns of 
+##' the matrix should be
 ##' referenced numerically
 ##' so this is not available to the \code{fun}
-##' @param r the series of spatial distances (or there maximums) we are
+##' @param r the series of spatial distances (or their maximums) we are
 ##'          interested in
 ##' @param r.low the low end of each range, 0 by default
 ##' @param data.frame logical indicating whether to return results as a data frame (default = TRUE)
@@ -37,7 +45,7 @@ get.pi <- function(posmat,
                    r.low=rep(0,length(r)),
                    data.frame=TRUE) {
 
-  xcol <-  which(colnames(posmat) == "x")
+  xcol <- which(colnames(posmat) == "x")
   ycol <- which(colnames(posmat) == "y")
 
   #check that both columns exist
@@ -70,8 +78,8 @@ get.pi <- function(posmat,
 ##' returns the odds that a point within a particular range of a point of interest shares the relationship
 ##' specified by the passed in function with that point.
 ##'
-##' @param posmat a matrix with columns x, y and any other named columns
-##'    columns needed by fun
+##' @param posmat a matrix with columns x, y and any other named
+##'    columns needed by \code{fun}
 ##' @param fun a function that takes in two rows of posmat and returns:
 ##' \enumerate{
 ##'      \item  for pairs that are (potentially) related
@@ -80,7 +88,7 @@ get.pi <- function(posmat,
 ##' Note that names from \code{posmat} are not preserved in calls to \code{fun}, so the columns of the matrix should be
 ##' referenced numerically
 ##' so this is not available to the fun
-##' @param r the series of spatial distances (or there maximums) we are
+##' @param r the series of spatial distances (or their maximums) we are
 ##'          interested in
 ##' @param r.low the low end of each range, 0 by default
 ##' @param data.frame logical indicating whether to return results as a data frame (default = TRUE)
@@ -230,93 +238,90 @@ get.theta.typed <- function(posmat,
 }
 
 
-##' Calculate bootstrapped confidence intervals for \code{get.pi} values.
+##' Calculate bootstrapped BCa confidence intervals from \code{get.pi} values.
 ##'
-##' Wrapper to \code{get.pi.bootstrap} that takes care of calculating the
-##' confidence intervals based on the bootstrapped values..
+##' Wrapper using \pkg{coxed} package to calculate the
+##' BCa (bias-corrected and accelerated) confidence interval (CI) for \eqn{\pi}(\code{r.low}, \code{r}), based on bootstrapped values from \code{get.pi.bootstrap}.
 ##'
+##' @param posmat a matrix with named columns x and y for 2D individual location
+##' @param fun the function to decide transmission-related pairs
+##' @param r the upper end of each distance band
+##' @param r.low the low end of each distance band (default: a vector of zeroes)
+##' @param boot.iter the number of bootstrap iterations (default = 1000)
+##' @param ci.level the level of the desired BCa CI (default = 0.95)
+##' @param data.frame logical: indicating whether to return results as a data frame (default = TRUE)
 ##'
-##' @param posmat a matrix with columns type, x and y
-##' @param fun the function to decide relationships
-##' @param r the series of spatial distances wer are interested in
-##' @param r.low the low end of each range. 0 by default
-##' @param boot.iter the number of bootstrap iterations
-##' @param ci.low the low end of the ci...0.025 by default
-##' @param ci.high the high end of the ci...0.975 by default
-##' @param data.frame logical indicating whether to return results as a data frame (default = TRUE)
+##' @return If \code{data.frame = TRUE} then a data frame of 5 variables \code{r.low}, \code{r}, \code{pt.est} (the point estimate from \code{get.pi}), the confidence envelope as \code{ci.low} and \code{ci.high}, with the observations representing ascending distance bands. Else a matrix with first row \code{ci.low} and second row \code{ci.high} with columns representing ascending distance bands.
 ##'
-##' @return a matrix with a row for the high and low values and
-##'     a column per distance
+##' @author Justin Lessler and Timothy M Pollington
 ##'
-##' @author Justin Lessler
+##' @references \href{https://arxiv.org/pdf/1911.08022v4.pdf#page=12}{Rationale for BCa rather than percentile CIs} is described in Pollington et al. (2020)
+##' Developments in statistical inference when assessing 
+##' spatiotemporal disease clustering with the tau statistic.
+##' *arXiv/stat.ME: 1911.08022v4*.
 ##'
 ##' @family get.pi
+##' 
+##' @section Depends on:
+##' coxed::bca()
 ##'
 ##' @example R/examples/get_pi_ci.R
-##'
+##' @md
 
 get.pi.ci <- function(posmat,
                       fun,
-                      r=1,
-                      r.low=rep(0,length(r)),
+                      r = 1,
+                      r.low = rep(0,length(r)),
                       boot.iter = 1000,
-                      ci.low=0.025,
-                      ci.high=0.975,
-                      data.frame=TRUE) {
+                      ci.level = 0.95,
+                      data.frame = TRUE) {
      
   boots <- get.pi.bootstrap(posmat, fun, r, r.low, boot.iter)
-
-  rc <- apply(boots[,-(1:2)], 1, quantile, probs=c(ci.low, ci.high))
+  boots = boots[,-(1:2)]
+  
+  rc <- apply(boots, 1, applyBCa, ci.level = 0.95)
   
   if (data.frame == FALSE) {
        return(rc)
   } else if (data.frame == TRUE) {
-       return(data.frame(r.low=r.low, 
-                         r=r, 
-                         pt.est=get.pi(posmat, fun, r, r.low)$pi, 
-                         ci.low=rc[1,], 
-                         ci.high=rc[2,]))
+       return(data.frame(r.low = r.low, 
+                         r = r, 
+                         pt.est = get.pi(posmat, fun, r, r.low)$pi, 
+                         ci.low = rc[1,], 
+                         ci.high = rc[2,]))
   }
 }
- 
+
 
 ##' Calculate bootstrapped confidence intervals for \code{get.theta} values.
 ##'
 ##' Wrapper to \code{get.theta.bootstrap} that takes care of calculating the
-##' confience intervals based on the bootstrapped values.
-##'
+##' confidence intervals based on the bootstrapped values.
 ##'
 ##' @param posmat a matrix with columns type, x and y
 ##' @param fun the function to decide relationships
 ##' @param r the series of spatial distances we are interested in
 ##' @param r.low the low end of each range. 0 by default
 ##' @param boot.iter the number of bootstrap iterations
-##' @param ci.low the low end of the ci...0.025 by default
-##' @param ci.high the high end of the ci...0.975 by default
+##' @param ci.level significance level of the 95% BCa CI, default = 0.95
 ##' @param data.frame logical indicating whether to return results as a data frame (default = TRUE)
 ##'
-##' @return a matrix with a row for the high and low values and
-##'     a column per distance
-##'
+##' @return a matrix with a row for the high and low values and a column per distance
 ##' @author Justin Lessler
-##'
 ##' @family get.theta
-##'
 ##' @example R/examples/get_theta_ci.R
-##'
 
 get.theta.ci <- function(posmat,
                          fun,
                          r=1,
                          r.low=rep(0,length(r)),
                          boot.iter = 1000,
-                         ci.low=0.025,
-                         ci.high=0.975,
+                         ci.level=0.95,
                          data.frame=TRUE) {
      
   boots <- get.theta.bootstrap(posmat, fun, r, r.low, boot.iter)
-
-  rc <- apply(boots[,-(1:2)], 1, quantile, probs=c(ci.low, ci.high))
+  boots = boots[,-(1:2)]
+  rc <- apply(boots, 1, applyBCa, ci.level = 0.95)
   
   if (data.frame == FALSE) {
        return(rc)
@@ -343,7 +348,9 @@ get.theta.ci <- function(posmat,
 ##' @param boot.iter the number of bootstrap iterations
 ##' @param data.frame logical indicating whether to return results as a data frame (default = TRUE)
 ##'
-##' @return pi values for all the distances we looked at
+##' @return Values of pi for all distance bands. Return value dependent on data.frame argument.
+##' Asa matrix (rows = bootstrap samples, columns = increasing distance bands)
+##' or a data.frame (r.low, r and increasing distance bands)
 ##'
 ##' @note In each bootstrap iteration N observations are drawn from the existing data with replacement. To avoid errors in
 ##' inference resulting from the same observatin being compared with itself in the bootstrapped data set, original indices
@@ -363,8 +370,7 @@ get.pi.bootstrap <- function(posmat,
                              boot.iter=500,
                              data.frame=TRUE) {
 
-
-  xcol <-  which(colnames(posmat)=="x")
+  xcol <- which(colnames(posmat)=="x")
   ycol <- which(colnames(posmat)=="y")
 
   #check that both columns exist
@@ -825,8 +831,8 @@ get.theta.typed.permute <- function(posmat,
 ##' from an index point share some relationship with that point versus
 ##' the probability (or odds) any point shares that relationship with that point.
 ##'
-##' @param posmat a matrix with columns x, y and any other named columns
-##'    columns needed by fun
+##' @param posmat a matrix with columns x, y and any other named
+##'    columns needed by \code{fun}
 ##' @param fun a function that takes in two rows of posmat and returns:
 ##' \enumerate{
 ##'      \item for pairs included in the numerator (and the denominator for independent data)
@@ -835,7 +841,7 @@ get.theta.typed.permute <- function(posmat,
 ##' Note that names from \code{posmat} are not preserved in calls to
 ##' \code{fun}, so the columns of the matrix should be referenced numerically
 ##' so this is not available to fun
-##' @param r the series of spatial distances (or there maximums) we are
+##' @param r the series of spatial distances (or their maximums) we are
 ##'          interested in
 ##' @param r.low the low end of each range, 0 by default
 ##' @param comparison.type what type of points are included in the comparison set.
@@ -843,9 +849,9 @@ get.theta.typed.permute <- function(posmat,
 ##'   \item "representative" if comparison set is representative of the underlying population
 ##'   \item "independent" if comparison set is cases/events coming from an indepedent process
 ##' }
-##' @param data.frame logical indicating whether to return results as a data frame (default = TRUE)
+##' @param data.frame logical indicating whether to return results 'like' a data frame format (default = TRUE)
 ##'
-##' @return The tau value for each distance we look at. If \code{comparison.type} is "representative", this is:
+##' @return The tau value for each distance we look at as a tau class with a matrix or data frame style. If \code{comparison.type} is "representative", this is:
 ##'
 ##' \code{tau = get.pi(posmat, fun, r, r.low)/get.pi(posmat,fun,infinity,0)}
 ##'
@@ -853,7 +859,7 @@ get.theta.typed.permute <- function(posmat,
 ##'
 ##' \code{tau = get.theta(posmat, fun, r, r.low)/get.theta(posmat,fun,infinity,0)}
 ##'
-##' @author Justin Lessler and Henrik Salje
+##' @author Justin Lessler, Timothy M Pollington and Henrik Salje
 ##'
 ##' @family get.tau
 ##' @family spatialtau
@@ -868,7 +874,7 @@ get.tau <- function(posmat,
                     comparison.type = "representative",
                     data.frame=TRUE) {
 
-  xcol <-  which(colnames(posmat)=="x")
+  xcol <- which(colnames(posmat)=="x")
   ycol <- which(colnames(posmat)=="y")
 
   #check that both columns exist
@@ -881,8 +887,8 @@ get.tau <- function(posmat,
   } else if (comparison.type == "independent") {
     comp.type.int <- 1
   } else {
-    stop("unkown comparison type specified")
-  }
+    stop("unknown comparison.type specified")
+  } 
 
   rc <- .Call("get_tau",
               posmat,
@@ -895,12 +901,306 @@ get.tau <- function(posmat,
               ycol)
   
   if (data.frame == FALSE) {
+       class(rc) <- "tau"
+       attr(rc, "comparison.type") = comparison.type
        return(rc)
   } else if (data.frame == TRUE) {
-       return(data.frame(r.low=r.low, r=r, tau=rc))
+       rc = data.frame(r.low=r.low, r=r, tau.pt.est=rc)
+       class(rc) <- "tau"
+       attr(rc, "comparison.type") = comparison.type
+       return(rc)
   }
 }
 
+##' Global hypothesis testing for the tau statistic
+##'
+##' Performs a graphical hypothesis test to assess the evidence against the null hypothesis (H_0: tau = 1 i.e. no spatiotemporal clustering nor inhibition). A global envelope test from the \code{GET} package is used to see if any part of the point estimate connected line is outside the lower or upper bounds of the global envelope. The global envelope is formed on the tau estimator acting on time-permuted data to simulate H_0. The global envelope test is of 'extreme rank type' i.e. minimum of pointwise ranks with 95\% significance level. 
+##'
+##' @param posmat a matrix with columns x, y and any other named
+##'    columns needed by \code{fun}
+##' @param fun a function that takes in two rows of posmat and returns:
+##' \enumerate{
+##'      \item for pairs included in the numerator (and the denominator for independent data)
+##'      \item for pairs that should only be included in the denominator
+##'      \item for pairs that should be ignored all together}
+##' Note that names from \code{posmat} are not preserved in calls to
+##' \code{fun}, so the columns of the matrix should be referenced numerically
+##' so this is not available to fun
+##' @param r the series of spatial distances (or their maximums) we are
+##'          interested in
+##' @param r.low the low end of each range, 0 by default
+##' @param permutations number of simulations of H_0. 2,500 is an optimal number according to Myllymäki et al. (2017).
+##' @param comparison.type what type of points are included in the comparison set.
+##' \itemize{
+##'   \item "representative" if comparison set is representative of the underlying population
+##'   \item "independent" if comparison set is cases/events coming from an indepedent process
+##' }
+##' @return An object of class \code{tauGET} which can then be plotted using \code{plot.tau()} and an additional \code{tau} class object. The object consists of:
+##' \itemize{
+##'   \item r that inputted earlier
+##'   \item obs the tau point estimate computed internally using \code{get.tau()}
+##'   \item central the median estimate of all simulation curves that represent the null hypothesis. Comparing this to the tau=1 line indicates if it is reasonable to assume that H_0 was adequately simulated.
+##'   \item lo the lower bound of the global envelope
+##'   \item hi the upper bound of the global envelope
+##'   \item tau.permute the entire record of simulations of H_0, to plot with the global envelope using \code{plot.tau()}.
+##' }
+##' @section Attributes:
+##' \itemize{
+##'   \item p_interval represents a range rather than a single p-value to assess the evidence against H_0. Accessed using \code{attr(x,"p_interval")}.
+##' }
+##' @author Timothy M Pollington
+##'
+##' @family get.tau
+##' @family spatialtau
+##' @example R/examples/get_tau_GET.R
+##'
+
+get.tau.GET <- function(posmat, fun, r, r.low, permutations = 2500, comparison.type){
+  get.tau = IDSpatialStats::get.tau(posmat = posmat, fun = fun, r = r, r.low = r.low, comparison.type = comparison.type, data.frame = FALSE)
+  tau.permute = IDSpatialStats::get.tau.permute(posmat = posmat, fun = fun, r = r, r.low = r.low, permutations = permutations, comparison.type = comparison.type, data.frame = FALSE)
+  curveset = GET::create_curve_set(list(r = r, obs = as.numeric(get.tau), sim_m = t(tau.permute)))
+  GET.res = GET::global_envelope_test(curve_sets = curveset, type = "rank", alpha = 0.05,
+           alternative = c("two.sided"), ties = "erl", probs = c(0.025, 0.975), quantile.type = 7, 
+           central = "median")
+  GET.res = list(GET.res, tau.permute)
+  class(GET.res) <- "tauGET"
+  return(GET.res)
+}
+
+##' Cluster range estimation using \code{get.tau.D.param.est}
+##'
+##' Estimates the range of spatiotemporal clustering. It records the place on the horizontal tau=1 line where each spatially bootstrapped simulation touches. This distribution then represents an empirical distribution for the clustering range and a confidence interval can be computed.  
+##'
+##' @param r the series of spatial distances (or their maximums) we are
+##'          interested in
+##' @param tausim the set of spatially-bootstrapped simulations. Has to be \code{taubstrap} class; use \code{get.tau.bootstrap(..., data.frame = FALSE)} to obtain this. 
+##' @param GETres is a required object and is obtained from a previous global hypothesis test using \code{get.tau.GET}. It ensures that the user has performed a graphical hypothesis test first and has considered there is evidence against H_0, before deciding to estimate the clustering range.
+##' @return An object of class \code{tauparamest} which can then be plotted using \code{plot.tau()}. The object consists of:
+##' \itemize{
+##'   \item envelope the distribution of clustering range estimates
+##' }
+##' @section Attributes:
+##' \itemize{
+##'   \item BCaCI the BCa CI for the distribution of clustering range estimates
+##' }
+##' @author Timothy M Pollington
+##'
+##' @family get.tau
+##' @family spatialtau
+##' @example R/examples/get_tau_D_param_est.R
+
+get.tau.D.param.est <- function(r, tausim, GETres = NULL){
+  stopifnot(!is.null(GETres)) # makes sure the user has been principled and performed a global
+  # hypothesis test using get.tau() before estimating D
+  stopifnot(length(r)>1)
+  stopifnot(class(tausim)=="taubstrap")
+  if(!is.null(names(tausim))){ # ie if tausim is like a 'data.frame despite having a taubstrap class
+    tausim = t(tausim[,-c(1,2)])
+  }
+  boot.iter = dim(tausim)[1]
+  ciIntercept <- function(boot.iter, r, tausim) {
+    j.max = length(r)
+    # define d.envelope by finding for each bootstrap sample the (interpolated) d-intercept point
+    alwaysabove1 = 0
+    d.envelope = NULL
+    for (i in 1:boot.iter) {
+      j = 1 # first distance band
+      if(tausim[i,j] > 1){ # else ignore simulation as starting from below tau = 1
+        stillabove1 = TRUE
+          while (stillabove1 & (j < j.max)) {
+            j = j + 1
+            if(tausim[i,j] <= 1){ # else it stays above tau = 1 until the next j is tested
+              stillabove1 = FALSE
+              root.tau1 = ((1-tausim[i,(j-1)])*(r[j]-r[j-1])/(tausim[i,j]-tausim[i,(j-1)]))+r[j-1]
+              d.envelope = c(d.envelope, root.tau1)
+            }
+          }
+          if(stillabove1 & j==j.max){
+            alwaysabove1 = alwaysabove1 + 1
+          }
+      }
+    }
+    print(paste0(length(d.envelope)/boot.iter*100, "% of boostrap sims crossing tau = 1 from above"))
+    print(paste0(alwaysabove1/boot.iter*100, "% of bootstrap sims always above tau = 1"))
+    if(alwaysabove1>0){
+      warning("Note that there are some bootstrap sims that stay above tau = 1 for the entire distance band set. If more than a few percent of these are above tau = 1 then a reliable CI cannot be constructed as it will have not have come from a random sample.")
+    }
+    return(d.envelope)
+  }
+  envelope = ciIntercept(boot.iter,r,tausim)
+  d.envelope = as.data.frame(envelope)
+  attr(d.envelope,"BCaCI") = coxed::bca(d.envelope$envelope, conf.level = 0.95)
+  
+  class(d.envelope) <- "tauparamest"
+  return(d.envelope)
+}
+  
+##' Plotting the results from tau functions
+##'
+##' Three types of plots:
+##' \enumerate{
+##' \item Diagnostic plot to indicate the structure or magnitude of spatiotemporal clustering. Requires \code{tau} object; \code{tauCI} object optional to draw pointwise CIs. This plot is only suitable for the purpose of a graphical hypothesis test in the situation that a specific distance band is selected prior to graph creation.
+##' \item Graphical hypothesis test to assess the evidence against the null hypothesis (no spatiotemporal clustering nor inhibition). Requires \code{tau} and \code{tauGET} objects.
+##' \item Estimation of the clustering range (the distribution of the places on the horizontal tau=1 line, where decreasing bootstrap simulations first intercept). Requires \code{tau}, \code{tauparamest} and \code{taubstrap} objects.
+##' }
+##'
+##' @param x \code{tau} object; create using \code{get.tau(..., data.frame = TRUE)}. Required for all plots.
+##' @param r.mid If \code{TRUE}(default) then for each point the x-coordinate of the midpoint of a distance band is plotted and if \code{FALSE} the endpoint of the distance band is plotted.
+##' @param tausim the set of spatially-bootstrapped simulations of \code{taubstrap} class; use \code{get.tau.bootstrap()} to obtain this. Required for Estimation of the clustering range plot.
+##' @param ptwise.CI the set of pointwise CIs of \code{tauCI} class; create using \code{get.tau(..., data.frame = TRUE)}. Optional for the diagnostic plot but should not be supplied for the other plots.
+##' @param GET.res is a required object for the graphical hypothesis test plot but should not be supplied for the other plots. It is obtained from \code{get.tau.GET(..., data.frame = TRUE)}. It ensures that the user has performed a graphical hypothesis test first and has considered there is evidence against H_0, before deciding to estimate the clustering range.
+##' @param d.param.est a required object for Estimating the clustering range plot from \code{get.tau.D.param(..., data.frame = TRUE)}, but should not be supplied for the other plots. A \code{taubstrap} object will also be necessary.
+##' @param ... other arguments which are standard for \code{plot()} for plot customisation
+##' @author Timothy M Pollington
+##'
+##' @family get.tau
+##' @family spatialtau
+##'
+
+plot.tau <- function(x, r.mid = TRUE, tausim = NULL, ptwise.CI = NULL, GET.res = NULL, d.param.est = NULL, ...)
+{
+  stopifnot(class(x)=="tau")
+  if(!is.null(ptwise.CI)){
+    stopifnot(class(ptwise.CI)=="tauCI")
+  }
+  if(!is.null(GET.res)){
+    stopifnot(class(GET.res)=="tauGET")
+  }
+  if(!is.null(d.param.est)){
+    stopifnot(class(d.param.est)=="tauparamest")
+  }
+  if(!is.null(tausim)){
+    stopifnot(class(tausim)=="taubstrap")
+  }
+  if(!is.null(ptwise.CI) & !is.null(GET.res)){
+    stop("To avoid misinterpretation of visual results we do not allow pointwise CIs and global envelopes to be plotted on the same graph")
+  }
+  if(!is.null(ptwise.CI) & !is.null(d.param.est)){
+    stop("To avoid misinterpretation of visual results we do not allow pointwise CIs and clustering range estimates to be plotted on the same graph")
+  }
+  if(!is.null(GET.res) & !is.null(d.param.est)){
+    stop("To avoid misinterpretation of visual results we do not allow global envelopes and clustering range estimates to be plotted on the same graph")
+  }
+  if(is.null(tausim) & !is.null(d.param.est)){
+    stop("Need tausim and d.param.est class objects to plot clustering range estimates")
+  }
+  if(r.mid==TRUE){
+    r.end = 0.5*(x$r.low + x$r)
+    midorend = "at distance band midpoint"
+    xlim = c(0,(max(r.end)*1.01))
+  }
+  else{
+    r.end = x$r
+    midorend = "at distance band endpoint"
+    xlim = c(0,(max(x$r)*1.01))
+  }
+  
+  # identify if the lower bound of each distance band contains zero or not, 
+  # and label graph appropriately, with correct units if provided
+  if(!is.null(attr(x$r.low, "units")) & !is.null(attr(x$r, "units")) & 
+     identical(attr(x$r.low, "units"), attr(x$r, "units"))){
+    unitslabel = c("(", attr(x$r.low, "units"), ")")
+  }
+  else{
+    unitslabel = ""
+  }
+  
+  if(all(x$r.low==0)){
+    xlab = bquote("Distance [0," * d[m] * ") from an average case " * .(unitslabel))
+  }
+  else{
+    xlab = bquote("Distance [" * d[l] * "," * d[m] * ") from an average case " * .(unitslabel))
+  }
+  
+  if(is.null(GET.res) | is.null(d.param.est)){
+    if(is.null(ptwise.CI)){
+    plot(x = r.end, y = x$tau.pt.est, xlim=xlim,
+       ylim=range(x$tau.pt.est, na.rm = TRUE)+diff(range(x$tau.pt.est, na.rm = TRUE))*c(-0.05,0.05),
+       cex.axis=1,col="black", xlab=xlab, ylab="Tau", 
+       cex.main=1, lwd=2, type="p", las=1, cex.axis=1, xaxs = "i", yaxs = "i", pch = 16)
+       abline(h=1,lty=2)
+       legend("topright",legend=bquote("point estimate" ~ hat(tau) * "," ~ .(midorend)),
+         col="black", pch=16)
+    }   
+    if(!is.null(ptwise.CI)){
+      ylimrange = range(c(x$tau.pt.est,ptwise.CI$ci.low,ptwise.CI$ci.high), na.rm = TRUE)
+      plot(x = r.end, y = x$tau.pt.est, xlim=xlim, 
+           ylim=ylimrange+diff(ylimrange)*c(-0.05,0.05),
+           cex.axis=1,col="black", xlab=xlab, ylab="Tau", 
+           cex.main=1, lwd=2, type="p", las=1, cex.axis=1, xaxs = "i", yaxs = "i", pch = 16)
+      arrows(r.end, ptwise.CI$ci.low, r.end, ptwise.CI$ci.high, length = 0.04, angle = 90, code = 3)
+      abline(h=1,lty=2)
+      legend("topright",legend=bquote("point estimate" ~ hat(tau) * "," ~ .(midorend)),
+        col="black", pch=16)
+    }
+  }
+  
+  if(!is.null(GET.res)){
+  permutations = dim(GET.res[[2]])[1]
+  plot(NULL, xlim = c(0,max(x$r, na.rm = TRUE)), ylim = c(min(GET.res[[1]]$lo, GET.res[[1]]$obs, na.rm = TRUE),max(GET.res[[1]]$hi, GET.res[[1]]$obs, na.rm = TRUE)), xaxt = "n", yaxt = "n", xaxs = "i", yaxs = "i",
+       ylab = "Tau", xlab = xlab, lwd = 4, cex.lab = 1.5)
+  for (i in 1:permutations) {
+    lines(x$r, GET.res[[2]][i,], col = scales::alpha("grey", alpha = 0.3), lwd = 1)
+  }
+  yaxis.range = c(min(GET.res[[1]]$lo, GET.res[[1]]$obs, na.rm = TRUE),max(GET.res[[1]]$hi, GET.res[[1]]$obs, na.rm = TRUE))
+  yaxis.lab = c(seq(yaxis.range[1],yaxis.range[2],length.out = 5),1)
+  yaxis.lab = sort(yaxis.lab)
+  yaxis.lab = round(yaxis.lab,digits = 1)
+  yaxis.lab = unique(yaxis.lab) # prevents more than one 1.0 value
+  yaxis.lab[which(yaxis.lab==1)] = round(yaxis.lab[which(yaxis.lab==1)],digits = 0)
+  axis(2, las=1, at=yaxis.lab, labels = as.character(yaxis.lab), lwd = 1)
+  lines(GET.res[[1]]$r, GET.res[[1]]$lo, col = "slategrey", lwd = 3)
+  lines(GET.res[[1]]$r, GET.res[[1]]$hi, col = "slategrey", lwd = 3)
+  lines(GET.res[[1]]$r, GET.res[[1]]$central, col = "red", lwd = 3)
+  lines(GET.res[[1]]$r, GET.res[[1]]$obs, lwd = 3)
+  axis(1, lwd = 1)
+  abline(h=1, lty = 2, lwd = 3)
+  legend("topright", legend=c(as.expression(bquote(~ hat(tau) ~ "point estimate")),
+                              "95% global envelope",as.expression(bquote("simulations of " ~ H[0])),
+                              "median simulation",
+                              as.expression(bquote(~ tau == 1)) ),
+         col=c("black", "slategrey", "grey", "red", "black"),
+         lty=c(1,1,1,1,2), cex=1.05, yjust = 0.5, lwd = 6)
+  par(xpd = TRUE)
+  pint.lo = round(attr(GET.res[[1]],"p_interval"), digits = 3)[1]
+  pint.hi = round(attr(GET.res[[1]],"p_interval"), digits = 3)[2]
+  pint.x = 0.5 * max(x$r, na.rm = TRUE)
+  pint.y = c(min(GET.res[[1]]$lo, GET.res[[1]]$obs, na.rm = TRUE),max(GET.res[[1]]$hi, GET.res[[1]]$obs, na.rm = TRUE))[1] + 0.5*diff(c(min(GET.res[[1]]$lo, GET.res[[1]]$obs, na.rm = TRUE),max(GET.res[[1]]$hi, GET.res[[1]]$obs, na.rm = TRUE)))
+  text(bquote("p-value in [" ~ .(pint.lo) * "," * .(pint.hi) * "]"), x = pint.x, y = pint.y)
+  }
+
+  if(!is.null(d.param.est) & !is.null(tausim)){
+    yaxis.range = c(min(x$tau.pt.est, tausim, na.rm = TRUE),max(x$tau.pt.est, tausim, 
+    na.rm = TRUE))
+    yaxis.lab = c(seq(yaxis.range[1],yaxis.range[2],length.out = 5),1)
+    yaxis.lab = sort(yaxis.lab)
+    yaxis.lab = round(yaxis.lab,digits = 1)
+    yaxis.lab = unique(yaxis.lab) # prevents more than one 1.0 value
+    yaxis.lab[which(yaxis.lab==1)] = round(yaxis.lab[which(yaxis.lab==1)],digits = 0)
+    plot(NULL, xlim = xlim, ylim = yaxis.range, xaxt = "n", yaxt = "n", xaxs = "i", yaxs = "i",
+    ylab = "", xlab = "")
+    mtext("Tau", side=2, line=3, cex = 1.5)
+    mtext(xlab, side=1, line=3, cex = 1.5)
+    for (i in 1:dim(tausim)[1]) {
+      lines(x$r, tausim[i,], col = scales::alpha("grey", alpha = 0.2), lwd = 4)
+    }
+    axis(2, las=1, at=yaxis.lab, labels = as.character(yaxis.lab), lwd = 1)
+    axis(1, lwd = 1)
+    lines(x = c(0,max(x$r, na.rm = TRUE)), y = c(1,1), lty = 2, lwd = 1) # as abline seems to overlap
+    par(lend=1);
+    lines(x = attr(d.param.est,"BCaCI"), y=c(1.03,1.03),
+          type = "l", lwd = 20, col = "red")
+    dintercept.ptest = median(d.param.est$envelope)
+    lines(x=c(dintercept.ptest,dintercept.ptest), y = c(0.9,1.1), lwd = 4)
+    lines(x$r, x$tau.pt.est, lwd = 4, col = "black")
+    legend("topright",
+           legend=c(as.expression(bquote(hat(tau) ~ "point estimate & " ~ hat(D) ~ "estimate")),
+                    as.expression(bquote(underline(tau)^"*" ~ "bootstrap estimate (N=" ~ .(dim(tausim)[1]) * ")")),
+                    as.expression(bquote("95% BCa CI of " ~ underline(D))),"tau = 1"), 
+           col=c("black", "grey", "red", "black"),
+           lty=c(1,1,1,2), lwd = c(2,2,10,1), pch = c(124,NA,NA,NA), cex=1.05, xjust = 1, yjust = 0.5)
+  }
+}
 
 ##' Optimized version of \code{get.tau} for typed data
 ##'
@@ -942,7 +1242,7 @@ get.tau.typed <- function(posmat,
      } else if (comparison.type == "independent") {
           comp.type.int <- 1
      } else {
-          stop("unkown comparison type specified")
+          stop("unknown comparison.type specified")
      }
      
      rc <- .C("get_tau_typed",
@@ -970,7 +1270,7 @@ get.tau.typed <- function(posmat,
 ##' Bootstrap confidence interval for the \code{get.tau} values
 ##'
 ##' Wrapper to \code{get.tau.bootstrap} that takes care of calulating
-##' the confidence intervals based on the bootstrapped values
+##' the confidence intervals based on the bootstrapped values.
 ##'
 ##' @param posmat a matrix appropriate for input to \code{get.tau}
 ##' @param fun a function appropriate as input to \code{get.pi}
@@ -978,8 +1278,7 @@ get.tau.typed <- function(posmat,
 ##' @param r.low the low end of each range....0  by default
 ##' @param boot.iter the number of bootstrap iterations
 ##' @param comparison.type the comparison type to pass to get.tau
-##' @param ci.low the low end of the ci...0.025 by default
-##' @param ci.high the high end of the ci...0.975 by default
+##' @param ci.level significance level of the BCa CI, default = 0.95
 ##' @param data.frame logical indicating whether to return results as a data frame (default = TRUE)
 ##'
 ##' @return a data frame with the point estimate of tau and its low and high confidence interval at each distance
@@ -997,24 +1296,21 @@ get.tau.ci <- function(posmat,
                        r.low=rep(0,length(r)),
                        boot.iter = 1000,
                        comparison.type = "representative",
-                       ci.low=0.025,
-                       ci.high=0.975,
+                       ci.level = 0.95,
                        data.frame=TRUE) {
      
-     boots <- get.tau.bootstrap(posmat, fun,
-                                r, r.low, boot.iter,
-                                comparison.type)
-     
-     rc <- apply(boots[,-(1:2)], 1, quantile, probs=c(ci.low, ci.high))
+     boots <- get.tau.bootstrap(posmat, fun, r, r.low, boot.iter, comparison.type, 
+                                data.frame = FALSE)
+     rc <- apply(boots, 2, applyBCa, ci.level)
      
      if (data.frame == FALSE) {
+          class(rc) <- "tauCI"
           return(rc)
      } else if (data.frame == TRUE) {
-          return(data.frame(r.low=r.low, 
-                            r=r, 
-                            pt.est=get.tau(posmat, fun, r, r.low)$tau, 
-                            ci.low=rc[1,], 
-                            ci.high=rc[2,]))
+          rc = data.frame(r.low=r.low, r=r, pt.est=get.tau(posmat, fun, r, r.low)$tau, 
+                  ci.low=rc[1,], ci.high=rc[2,])
+          class(rc) <- "tauCI"
+          return(rc)
      }
 }
 
@@ -1051,8 +1347,7 @@ get.tau.bootstrap <- function(posmat,
                               comparison.type = "representative",
                               data.frame=TRUE) {
 
-
-  xcol <-  which(colnames(posmat)=="x")
+  xcol <- which(colnames(posmat)=="x")
   ycol <- which(colnames(posmat)=="y")
 
   #check that both columns exist
@@ -1065,7 +1360,7 @@ get.tau.bootstrap <- function(posmat,
   } else if (comparison.type == "independent") {
     comp.type.int <- 1
   } else {
-    stop("unkown comparison type specified")
+    stop("unknown comparison type specified")
   }
 
   rc <- matrix(nrow=boot.iter, ncol=length(r))
@@ -1083,9 +1378,12 @@ get.tau.bootstrap <- function(posmat,
   }
   
   if (data.frame == FALSE) {
+       class(rc) <- "taubstrap"
        return(rc)
   } else if (data.frame == TRUE) {
-       return(data.frame(r.low=r.low, r=r, t(rc)))
+       rc = data.frame(r.low=r.low, r=r, t(rc))
+       class(rc) <- "taubstrap"
+       return(rc)
   }
 }
 
@@ -1129,7 +1427,7 @@ get.tau.typed.bootstrap <- function(posmat,
   } else if (comparison.type == "independent") {
     comp.type.int <- 1
   } else {
-    stop("unkown comparison type specified")
+    stop("unknown comparison type specified")
   }
 
   rc <- matrix(nrow=boot.iter, ncol=length(r))
@@ -1280,7 +1578,7 @@ get.tau.typed.permute <- function(posmat,
   } else if (comparison.type == "independent") {
     comp.type.int <- 1
   } else {
-    stop("unkown comparison type specified")
+    stop("unknown comparison type specified")
   }
 
   rc <- matrix(nrow=permutations, ncol=length(r))
diff --git a/README.md b/README.md
index 8c35ccd..f6e754a 100644
--- a/README.md
+++ b/README.md
@@ -1,34 +1,64 @@
-## IDSpatialStats
+# IDSpatialStats
 
-This GitHub repository provides source code for the `IDSpatialStats` R package, which is designed to help epidemiologists assess the scale of spatial and temporal dependence in epidemic case occurrence data. 
+**Previous users: please read [news on the latest release](../master/NEWS.md "News on the latest release") to update you on changes to the tau statistic functions.**
 
-The current implementation of the package includes a function which simulates infectious disease spread as a spatial branching process, along with two novel spatial statistics that estimate: 1) the mean of the spatial transmission kernel, which is a measure of fine-scale spatial dependence between two cases, and 2) the tau-statistic, a measure of global clustering based on pathogen subtype.
+This GitHub repository provides source code for the `IDSpatialStats` R package, which helps epidemiologists assess the scale of spatial and temporal dependence in epidemic case occurrence data. This package can simulate infectious disease spread as a spatial branching process, along with two novel spatial statistics that estimate:
 
-Detailed description of the methods can be found here:
+1. the mean of the spatial transmission kernel, which is a measure of fine-scale spatial dependence between two cases, and
+2. the tau statistic &tau;, a measure of global clustering based on any/all of pathogen subtype; serotype; case onset time.
 
+This package is maintained by John Giles (GitHub: @[gilesjohnr](https://github.com/gilesjohnR)) and Justin Lessler (GitHub: @[jlessler](https://github.com/jlessler)) as part of the Johns Hopkins Bloomberg School of Public Health Infectious Disease Dynamics team (GitHub: @[HopkinsIDD](https://github.com/HopkinsIDD)).
+
+## Detailed description of the methods and relevant literature
+### Transmission distance estimation and tau statistic
 [The IDSpatialStats R package: Quantifying spatial dependence of infectious disease spread (Giles et al. 2019)](https://journal.r-project.org/archive/2019/RJ-2019-043/index.html)
 
 [Measuring spatial dependence for infectious disease epidemiology (Lessler et al. 2016)](http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0155249)
 
+### Transmission distance specific
 [Estimating infectious disease transmission distances using the overall distribution of cases (Salje et al. 2016)](http://www.sciencedirect.com/science/article/pii/S1755436516300317)
 
-### Installation
+### Tau statistic specific
+[Developments in statistical inference when assessing spatiotemporal disease clustering with the tau statistic (Pollington et al. 2020 pre-proof)](https://doi.org/10.1016/j.spasta.2020.100438)
+
+[The spatiotemporal tau statistic: a review (Pollington et al. preprint)](https://arxiv.org/abs/1911.11476)
+
+Three tau statistic plots are illustrated below for diagnostic, graphical hypothesis tests and clustering range parameter estimation purposes [(image source, CC-BY licence)](https://doi.org/10.1016/j.spasta.2020.100438):
+<img src="https://ars.els-cdn.com/content/image/1-s2.0-S2211675320300324-fx1_lrg.jpg" width="720" align="left" title="Click for larger image">
+<br clear="left">
 
-To install the offical release of the `IDSpatialStats` package, open `R` and type:
-```r
+## Installation
+
+To install the official release of the `IDSpatialStats` package, open `R` and type:
+```
 install.packages('IDSpatialStats')
 ```
-
-To install the install the development version, first install the `devtools` package and then install `IDSpatialStats` from source via GitHub:
-```r
+or for the development version, first install the `devtools` package and then `IDSpatialStats` from source via GitHub:
+```
 install.packages('devtools')
 devtools::install_github('HopkinsIDD/IDSpatialStats')
 ```
 
-### Troubleshooting
+## Troubleshooting and contributions
 
 For general questions, contact package maintainers Justin Lessler (jlessler@unc.edu) or John Giles (jrgiles@uw.edu).
 
-To report bugs or problems with documentation, please go to the [Issues](https://github.com/HopkinsIDD/IDSpatialStats/issues) page associated with this GitHub page and click *new issue*.
+To report bugs or problems with documentation, please go to the [Issues](https://github.com/HopkinsIDD/IDSpatialStats/issues) page associated with this GitHub page and click *New issue*. Bugs clearly reported using the `reprex` package are encouraged.
+
+If you wish to contribute to `IDSpatialStats`, please first get in touch via email. Then if we agree in principle please:
+
+1. Fork a copy of the current *development* version on GitHub
+2. Add your functions and edits to your forked copy
+  * pay attention to existing naming conventions and outputs
+  * add examples
+  * line comments are welcome for non-intuitive commands.
+  * commit to your own forked version:
+    * often
+    * describe what was done and why, but not how
+    * use the imperative
+    * &lt; 72 characters  
+      e.g. "*Replace percentile CI with BCa CI, as tau bootstrap distrib. non-symm*"
 
-If you wish to contribute to `IDSpatialStats`, please get in touch via email and then fork the latest version of the package. After committing your code to your own forked version, submit a pull request when you are ready to share.
+3. Any modified functions must return identical output as the current functions. Check that modified functions return the same output using package `testthat` and consider writing new test cases if appropriate. For new functions, write test cases should test that functions return expected values given expected inputs, and that they behave as expected in boundary conditions.
+4. Add conditional stops to functions so that they fail gracefully with unexpected inputs.
+5. Submit a pull request when you are ready to share. Thank you!
diff --git a/inst/CITATION b/inst/CITATION
new file mode 100644
index 0000000..be94c03
--- /dev/null
+++ b/inst/CITATION
@@ -0,0 +1,60 @@
+bibentry(
+   header = "To cite IDSpatialStats in publications please use:",
+   key = "Giles2019",
+   bibtype = "Unpublished",
+   title = "The IDSpatialStats R Package: Quantifying Spatial Dependence
+          of Infectious Disease Spread",
+   author = c(as.person("John R. Giles"),as.person("Henrik Salje"),as.person("Justin Lessler")),
+   year = 2019,
+   doi = "10.32614/RJ-2019-043",
+   journal = "The R Journal",
+   url = "https://journal.r-project.org/archive/2019/RJ-2019-043/index.html",
+   note = "Accepted article",
+)
+
+bibentry(
+   header = "If you calculate mean transmission distances using IDSpatialStats please also cite:",
+   key = "Salje2016",
+   bibtype = "Article",
+   title = "Estimating infectious disease transmission distances using the overall distribution of cases",
+   author = c(as.person("Henrik Salje"),as.person("Derek A.T. Cummings"),as.person("Justin Lessler")),
+   year = 2016,
+   doi = "10.1016/j.epidem.2016.10.001",
+   journal = "Epidemics",
+   url = "https://www.sciencedirect.com/science/article/pii/S1755436516300317"
+)
+
+c(bibentry(
+   header = "For the two foundational papers responsible for the tau statistic please also cite:",
+   key = "Salje2012",
+   bibtype = "Article",
+   title = "Revealing the microscale spatial signature of dengue transmission and immunity in an urban population",
+   author = c(as.person("Henrik Salje"),as.person("Justin Lessler"),as.person("Timothy P. Endy"),as.person("Frank C. Curriero"),as.person("Robert V. Gibbons"),as.person("Ananda Nisalak"),as.person("Suchitra Nimmannitya"),as.person("Siripen Kalayanarooj"),as.person("Richard G. Jarman"),as.person("Stephen J. Thomas"),as.person("Donald S. Burke"),as.person("Derek A. T. Cummings")),
+   year = 2012,
+   doi = "doi.org/10.1073/pnas.1120621109",
+   journal = "PNAS",
+   url = "https://www.pnas.org/content/109/24/9535"
+   ),
+   bibentry(
+   key = "Lessler2016",
+   bibtype = "Article",
+   title = "Measuring Spatial Dependence for Infectious Disease Epidemiology",
+   author = c(as.person("Justin Lessler"),as.person("Henrik Salje"),as.person("M. Kate Grabowski"),as.person("Derek A. T. Cummings")),
+   year = 2016,
+   doi = "10.1371/journal.pone.0155249",
+   journal = "PLoS ONE",
+   url = "https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0155249"
+   )
+)
+
+bibentry(
+   header = "For the method and rationale of graphical hypothesis testing or clustering range estimation for the tau statistic please also cite:",
+   key = "Pollington2020",
+   bibtype = "Article",
+   title = "Developments in statistical inference when assessing spatiotemporal disease clustering with the tau statistic",
+   author = c(as.person("Timothy M. Pollington"),as.person("Michael J. Tildesley"),as.person("T. Déirdre Hollingsworth"),as.person("Lloyd A. C. Chapman")),
+   year = 2020,
+   journal = "Spatial Statistics",
+   url = "https://doi.org/10.1016/j.spasta.2020.100438",
+   doi = "10.1016/j.spasta.2020.100438"
+)
\ No newline at end of file
diff --git a/inst/tests/test-getpi.r b/inst/tests/test-getpi.r
deleted file mode 100644
index 9c83f10..0000000
--- a/inst/tests/test-getpi.r
+++ /dev/null
@@ -1,277 +0,0 @@
-context("get.pi")
-test_that("get.pi returns 1 when labels are ignored", {
-
-    #generate a set of 100 random points even labeled between the two
-    x<-cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
-
-    colnames(x) <-c("type","x","y")
-
-    test <- function(a,b) {return(1)}
-
-    #with no lower limit
-    res <- get.pi(x,test,seq(10,100,10))$pi
-    expect_that(res,equals(rep(1,10)))
-
-    #with lower and upper limit
-    res <- get.pi(x,test,seq(10,100,10), seq(0,90,10))$pi
-    expect_that(res,equals(rep(1,10)))
-})
-
-test_that("get.pi returns appropriate values for cannonical test case 1 (equilateral triangle)", {
-
-    x <- rbind(c(1,0,0), c(1,1,0),c(2,.5,sqrt(.75)))
-    colnames(x) <-c("type","x","y")
-
-    test <- function(a,b) {
-        if (a[1] != 1) return(3)
-        if (b[1] == 2) return(1)
-        return(2)
-    }
-
-    #first no lower limit
-    res <- get.pi(x,test,1.5)$pi
-    res2 <- get.pi.typed(x,1,2,1.5)$pi
-
-    expect_that(res, equals(.5))
-    expect_that(res2, equals(.5))
-
-    #now with a lower limit
-
-    res <- get.pi(x,test,1.5,.5)$pi
-    res2 <- get.pi.typed(x,1,2,1.5,.5)$pi
-
-    expect_that(res, equals(.5))
-    expect_that(res2, equals(.5))
-
-})
-
-test_that("get.pi returns appropriate values cannonical test case 2 (points on a line)", {
-    x<-rbind(c(1,0,0), c(2,1,0), c(2,-1,0), c(3,2,0),
-             c(2,-2,0), c(3,3,0),c(3,-3,0))
-
-    colnames(x) <-c("type","x","y")
-
-    test <- function(a,b) {
-        if (a[1] != 1) return(3)
-        if (b[1] == 2) return(1)
-        return(2)
-    }
-
-    #pi 0,1.5 should be 1, 1.5-2.5 should be 0.5 and 2.5+ should be 0
-    res <- get.pi(x, test, c(1.5,2.5,Inf), c(0,1.5,2.5))$pi
-    res2 <- get.pi.typed(x, 1, 2, c(1.5,2.5,1000), c(0,1.5,2.5))$pi
-
-    expect_that(res,equals(c(1,0.5,0)))
-    expect_that(res2,equals(c(1,0.5,0)))
-
-})
-
-
-test_that("get.pi and get.pi.typed have same behavior on random data", {
-
-    #generate a set of 1000 random points even labeled between the two
-    x<-cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
-
-    colnames(x) <-c("type","x","y")
-
-    test <- function(a,b) {
-        if (a[1] != 1) return(3)
-        if (b[1] == 1) return(1)
-        return(2)
-    }
-
-    #no lower limit
-    res1 <- get.pi(x,test,seq(10,100,10))$pi
-    res2 <- get.pi.typed(x, 1,1, seq(10,100,10))$pi
-    expect_that(res1,equals(res2))
-
-    #lower limit
-    res1 <- get.pi(x,test,seq(10,100,10), seq(0,90,10))$pi
-    res2 <- get.pi.typed(x, 1,1, seq(10,100,10), seq(0,90,10))$pi
-    expect_that(res1,equals(res2))
-})
-
-test_that("get.pi returns identical results regardless of column order",
-          {
-              x<-cbind(rep(c(1,2),50), x=runif(100,0,100),
-                       y=runif(100,0,100))
-
-              colnames(x) <-c("type","x","y")
-
-              test <- function(a,b) {
-                  if (a[1] != 1) return(3)
-                  if (b[1] == 1) return(1)
-                  return(2)
-              }
-
-              res1 <- get.pi(x,test,seq(10,100,10), seq(0,90,10))$pi
-
-              test <- function(a,b) {
-                  if (a[3] != 1) return(3)
-                  if (b[3] == 1) return(1)
-                  return(2)
-              }
-
-              res2 <- get.pi(x[,c(3,2,1)],test,seq(10,100,10), seq(0,90,10))$pi
-
-              test <- function(a,b) {
-                  if (a[2] != 1) return(3)
-                  if (b[2] == 1) return(1)
-                  return(2)
-              }
-
-              res3 <- get.pi(x[,c(2,1,3)],test,seq(10,100,10), seq(0,90,10))$pi
-
-              expect_that(res1, equals(res2))
-              expect_that(res2, equals(res3))
-
-          })
-
-
-test_that ("get.pi fails nicely if x and y column names are not provided", {
-    x<-cbind(rep(c(1,2),500), a=runif(1000,0,100), b=runif(1000,0,100))
-
-    test <- function(a,b) {
-        if (a[1] != 2) return(3)
-        if (b[1] == 3) return(1)
-        return(2)
-    }
-
-    expect_that(get.pi(x,test,seq(10,50,10), seq(0,40,10))$pi,
-                throws_error("unique x and y columns must be defined"))
-
-})
-
-
-##################DEPRECATED TESTS...TAKE TO LONG...NOW USING SMALLER CANONICAL
-##################TESTS THAT HAVE VALUES THAT CAN BE WORKED OUT BY HAND
-## test_that("get.pi and get.pi.typed have same behavior and both return about .5 when they should", {
-
-##     set.seed(787)
-
-##     #generate a set of 1000 random points even labeled between the two
-##     x<-cbind(rep(c(1,2),500), x=runif(1000,0,100), y=runif(1000,0,100))
-
-##     colnames(x) <-c("type","x","y")
-
-##     test <- function(a,b) {
-##         if (a[1] != 1) return(3)
-##         if (b[1] == 1) return(1)
-##         return(2)
-##     }
-
-##     #no lower limit
-##     res1 <- get.pi(x,test,seq(10,100,10))
-##     res2 <- get.pi.typed(x, 1,1, seq(10,100,10))
-##     expect_that(res1,equals(res2))
-##     expect_that(round(res1,1),equals(rep(.5,10)))
-##     expect_that(round(res2,1),equals(rep(.5,10)))
-
-
-##     #lower limit
-##     res1 <- get.pi(x,test,seq(10,100,10), seq(0,90,10))
-##     res2 <- get.pi.typed(x, 1,1, seq(10,100,10), seq(0,90,10))
-##     expect_that(res1,equals(res2))
-##     expect_that(round(res1,1),equals(rep(.5,10)))
-##     expect_that(round(res2,1),equals(rep(.5,10)))
-## })
-
-
-## test_that("get.pi with equality test returns about .5 when points are uniform",
-##       {
-##           set.seed(787)
-
-##           #generate a set of 1000 random points even labeled between the two
-##           x<-cbind(rep(c(1,2),500), x=runif(1000,0,100), y=runif(1000,0,100))
-
-##           colnames(x) <-c("type","x","y")
-
-##           test <- function(a,b) {
-##               if (a[1]==b[1]) return(1)
-##               return(2)
-##           }
-
-##           #no lower limit
-##           res1 <- get.pi(x,test,seq(10,100,10))
-##           expect_that(round(res1,1),equals(rep(.5,10)))
-
-##           #lower limit
-##           res1 <- get.pi(x,test,seq(10,100,10), seq(0,90,10))
-##           expect_that(round(res1,1),equals(rep(.5,10)))
-##       })
-
-
-## test_that("get.pi is montonically decreasing for normally distributed clusters",
-##           {
-##               set.seed(787)
-
-##               #first generate 500 random uniform points
-##               x<-cbind(1, x=runif(500,0,100), y=runif(500,0,100))
-##               colnames(x) <-c("type","x","y")
-
-##               #add a seed point
-##               x<-rbind(x,c(2,50,50))
-
-##               #generate 500 normally distibuted points around this
-##               x<-rbind(x,cbind(3,rnorm(1000,50,20),rnorm(1000,50,20)))
-
-##               #check wit get.pi.typed
-##               res1 <- get.pi.typed(x,2,3,seq(10,50,10), seq(0,40,10))
-
-##               expect_that(res1[1]>res1[2] & res1[2]>res1[3] &
-##                           res1[3]>res1[4] & res1[4]>res1[5], is_true())
-
-
-##               #do test for not pi version
-##               test <- function(a,b) {
-##                   if (a[1] != 2) return(3)
-##                   if (b[1] == 3) return(1)
-##                   return(2)
-##               }
-
-##               res2 <- get.pi(x,test,seq(10,50,10), seq(0,40,10))
-
-
-##               expect_that(res2[1]>res2[2] & res2[2]>res2[3] &
-##                           res2[3]>res2[4] & res2[4]>res2[5], is_true())
-##               expect_that(res1,equals(res2))
-##           })
-
-## test_that("get.pi returns identical results regardless of column order",
-##           {
-##               set.seed(787)
-
-##               x<-cbind(rep(c(1,2),500), x=runif(1000,0,100), y=runif(1000,0,100))
-##               colnames(x) <-c("type","x","y")
-
-##               test <- function(a,b) {
-##                   if (a[1] != 1) return(3)
-##                   if (b[1] == 1) return(1)
-##                   return(2)
-##               }
-
-##               res1 <- get.pi(x,test,seq(10,100,10), seq(0,90,10))
-
-##               test <- function(a,b) {
-##                   if (a[3] != 1) return(3)
-##                   if (b[3] == 1) return(1)
-##                   return(2)
-##               }
-
-##               res2 <- get.pi(x[,c(3,2,1)],test,seq(10,100,10), seq(0,90,10))
-
-##               test <- function(a,b) {
-##                   if (a[2] != 1) return(3)
-##                   if (b[2] == 1) return(1)
-##                   return(2)
-##               }
-
-##               res3 <- get.pi(x[,c(2,1,3)],test,seq(10,100,10), seq(0,90,10))
-
-##               expect_that(res1, equals(res2))
-##               expect_that(res2, equals(res3))
-
-##           })
-
-
-
diff --git a/inst/tests/test-getpipermute.r b/inst/tests/test-getpipermute.r
deleted file mode 100644
index 2701617..0000000
--- a/inst/tests/test-getpipermute.r
+++ /dev/null
@@ -1,156 +0,0 @@
-context("get.pi.permute")
-
-test_that("get.pi.permute returns appropriate values for test case 1 (equilateral triangle)" ,{
-
-    x <- rbind(c(1,0,0), c(1,1,0),c(2,.5,sqrt(.75)))
-    colnames(x) <-c("type","x","y")
-
-    test <- function(a,b) {
-        if (a[1] != 1) return(3)
-        if (b[1] == 2) return(1)
-        return(2)
-    }
-
-
-    #should return .5 for every permutation
-    res <- get.pi.permute(x, test, 1.5, 0, 500)[,-(1:2)]
-    res2 <- get.pi.typed.permute(x, 1, 2, 1.5, 0, 500)[,-(1:2)]
-
-    expect_that(as.numeric(res), equals(rep(0.5,500)))
-    expect_that(as.numeric(res2), equals(rep(0.5,500)))
-
-})
-
-test_that("get.pi.permute returns appropriate values for test case 2 (points on a line)" ,{
-    x<-rbind(c(1,0,0), c(2,1,0), c(2,-1,0), c(3,2,0),
-             c(2,-2,0), c(3,3,0),c(3,-3,0))
-
-    colnames(x) <-c("type","x","y")
-
-    test <- function(a,b) {
-        if (a[1] != 1) return(3)
-        if (b[1] == 2) return(1)
-        return(2)
-    }
-
-    #the mean of the null distribution should be 0.5
-    #the 95% CI equals 0,1 with windows
-    res <- get.pi.permute(x, test, c(1.5,2.5,3.5), c(0,1.5,2.5), 500)[,-(1:2)]
-    res2 <- get.pi.typed.permute(x, 1, 2, c(1.5,2.5,3.5), c(0,1.5,2.5), 500)[,-(1:2)]
-
-    expect_that(rowMeans(res,na.rm=T), equals(rep(.5,3), tolerance=0.1))
-    expect_that(rowMeans(res2, na.rm=T), equals(rep(.5,3), tolerance=0.1))
-
-    for (i in 1:3) {
-        expect_that(as.numeric(quantile(res[i,], probs=c(.025,.975))),
-                    equals(c(0,1)))
-        expect_that(as.numeric(quantile(res2[i,], probs=c(.025,.975))),
-                    equals(c(0,1)))
-    }
-
-    #without windows the 95% CI should be around 2* 0.5+/- 1/sqrt(4) * 0.25
-    #since quantiles, that is 0.25 and 0.75
-    res <- get.pi.permute(x, test, 4,0, 500)
-    res2 <- get.pi.typed.permute(x, 1, 2, 4,0, 500)
-    expect_that(as.numeric(quantile(res[1,], probs=c(.025,.975))),
-                equals(c(0.25,0.75)))
-    expect_that(as.numeric(quantile(res2[1,], probs=c(.025,.975))),
-                equals(c(0.25,0.75)))
-})
-
-
-
-test_that ("fails nicely if x and y column names are not provided", {
-    x<-cbind(rep(c(1,2),500), a=runif(1000,0,100), b=runif(1000,0,100))
-
-    test <- function(a,b) {
-        if (a[1] != 2) return(3)
-        if (b[1] == 3) return(1)
-        return(2)
-    }
-
-    expect_that(get.pi.permute(x,test,seq(10,50,10), seq(0,40,10),100),
-                throws_error("unique x and y columns must be defined"))
-
-})
-
-
-
-##################DEPRECATED TESTS...TAKE TO LONG...NOW USING SMALLER CANONICAL
-##################TESTS THAT HAVE VALUES THAT CAN BE WORKED OUT BY HAND
-## test_that("get.pi.permute cis enclose get.pi when no clustering exists",
-##           {
-##               set.seed(787)
-
-##               x<-cbind(rep(c(1,2),250), x=runif(500,0,100), y=runif(500,0,100))
-
-##               colnames(x) <-c("type","x","y")
-
-##               test <- function(a,b) {
-##                   if (a[1] != 1) return(3)
-##                   if (b[1] == 1) return(1)
-##                   return(2)
-##               }
-
-##               #plot(x[,"x"],x[,"y"], col=x[,"type"])
-
-##               res <- get.pi.permute(x, test, seq(10,100,10), seq(0,90,10), 300)
-##               res2 <- get.pi(x, test, seq(10,100,10), seq(0,90,10))
-
-##               for (i in 1:10) {
-##                   tmp <- quantile(res[,i], probs=c(0.025, .975), na.rm=T)
-##                   print(res2[i])
-##                   print(tmp)
-##                   expect_that(res2[i]>=tmp[1], is_true())
-##                   expect_that(res2[i]<=tmp[2], is_true())
-##               }
-##           })
-
-## test_that("get.pi.permute cis do not enclose get.pi at extremes when no clustering exists",
-##           {
-##               set.seed(787)
-
-##               #first generate 200 random uniform points
-##               x<-cbind(1, x=runif(200,0,100), y=runif(200,0,100))
-##               colnames(x) <-c("type","x","y")
-
-##               #add a seed point
-##               x<-rbind(x,c(2,50,50))
-
-##               #generate 200 normally distibuted points around this
-##               x<-rbind(x,cbind(3,rnorm(200,50,20),rnorm(200,50,20)))
-
-##               test <- function(a,b) {
-##                   if (a[1] != 2) return(3)
-##                   if (b[1] == 3) return(1)
-##                   return(2)
-##               }
-
-##               res <- get.pi.permute(x,test,seq(10,50,10), seq(0,40,10), 200)
-##               res2 <- get.pi(x,test,seq(10,50,10), seq(0,40,10))
-
-##               #print(res)
-##               #print(res2)
-
-##               for (i in c(1,5)) {
-##                   tmp <- quantile(res[,i], probs=c(0.025, .975), na.rm=T)
-##                   expect_that((res2[i]>=tmp[1]) & (res2[i]<=tmp[2]) ,
-##                               is_false())
-##               }
-
-
-##               res <- get.pi.typed.permute(x,2,3,seq(10,50,10),
-##                                            seq(0,40,10), 100)
-
-##               for (i in c(1,5)) {
-##                   tmp <- quantile(res[,i], probs=c(0.025, .975), na.rm=T)
-##                   #print(res2[i])
-##                   #print(tmp)
-##                   expect_that((res2[i]>=tmp[1]) & (res2[i]<=tmp[2]) ,
-##                               is_false())
-##               }
-
-
-##           })
-
-
diff --git a/man/DengueSimR01.Rd b/man/DengueSimR01.Rd
index 24fa3a4..6c760e3 100644
--- a/man/DengueSimR01.Rd
+++ b/man/DengueSimR01.Rd
@@ -5,17 +5,17 @@
 \alias{DengueSimR01}
 \title{Simulated dataset of dengue transmission with basic reproductive number of 1}
 \format{
-Matrix with five columns representing the X and Y coordinates of infected individuals, the time of infection, 
+Matrix with five columns representing the X and Y coordinates of infected individuals, the time of infection,
 the genotype of the infecting pathogen and the serotype of the infecting pathogen.
 }
 \usage{
 DengueSimR01
 }
 \description{
-Dataset simulated using an agent based model with a spatially heterogeneous population structure. Infectious agents 
-were introduced resulting in agent to agent transmission. The distance between successive cases in a transmission chain 
-were randomly drawn from a uniform distribution U(0,100). Each infectious agent resulted in a single transmission to 
-another agent after a delay of 15 days, reflecting the generation time of dengue. There are 11 transmission chains, 
+Dataset simulated using an agent based model with a spatially heterogeneous population structure. Infectious agents
+were introduced resulting in agent to agent transmission. The distance between successive cases in a transmission chain
+were randomly drawn from a uniform distribution U(0,100). Each infectious agent resulted in a single transmission to
+another agent after a delay of 15 days, reflecting the generation time of dengue. There are 11 transmission chains,
 each with a different genotype. The genotypes are subdivided into four serotypes.
 }
 \author{
diff --git a/man/DengueSimR02.Rd b/man/DengueSimR02.Rd
index 40037bf..0f9480a 100644
--- a/man/DengueSimR02.Rd
+++ b/man/DengueSimR02.Rd
@@ -5,17 +5,17 @@
 \alias{DengueSimR02}
 \title{Simulated dataset of dengue cases with basic reproductive number of 2}
 \format{
-Matrix with five columns representing the X and Y coordinates of infected individuals, the time of infection, 
+Matrix with five columns representing the X and Y coordinates of infected individuals, the time of infection,
 the genotype of the infecting pathogen and the serotype of the infecting pathogen.
 }
 \usage{
 DengueSimR02
 }
 \description{
-Dataset simulated using an agent based model with a spatially heterogeneous population structure. Infectious agents 
-were introduced resulting in agent to agent transmission. The distance between successive cases in a transmission chain 
-were randomly drawn from a uniform distribution U(0,100). Each infectious agent resulted in transmissions to two other 
-agents after a delay of 15 days, reflecting the generation time of dengue. There are 11 transmission chains, each with 
+Dataset simulated using an agent based model with a spatially heterogeneous population structure. Infectious agents
+were introduced resulting in agent to agent transmission. The distance between successive cases in a transmission chain
+were randomly drawn from a uniform distribution U(0,100). Each infectious agent resulted in transmissions to two other
+agents after a delay of 15 days, reflecting the generation time of dengue. There are 11 transmission chains, each with
 a different genotype. The genotypes are subdivided into four serotypes.
 }
 \author{
diff --git a/man/DengueSimRepresentative.Rd b/man/DengueSimRepresentative.Rd
index 30cdf48..2c9672b 100644
--- a/man/DengueSimRepresentative.Rd
+++ b/man/DengueSimRepresentative.Rd
@@ -5,19 +5,19 @@
 \alias{DengueSimRepresentative}
 \title{Simulated dataset of dengue cases with representative underlying population}
 \format{
-Matrix with five columns representing the X and Y coordinates of infected individuals, the time of infection, the genotype 
-of the infecting pathogen and the serotype of the infecting pathogen. Individuals representative from the underlying population 
-have '-999'for time, genotype and serotype.
+Matrix with five columns representing the X and Y coordinates of infected individuals, the time of infection, the genotype
+of the infecting pathogen and the serotype of the infecting pathogen. Individuals representative from the underlying population
+have '-999' for time, genotype and serotype.
 }
 \usage{
 DengueSimRepresentative
 }
 \description{
-Dataset simulated using an agent based model with a spatially heterogeneous population structure. Infectious agents were 
-introduced resulting in agent to agent transmission. The distance between successive cases in a transmission chain were 
-randomly drawn from a uniform distribution U(0,100). Each infectious agent resulted in transmissions to two other agents 
-after a delay of 15 days, reflecting the generation time of dengue. There are 11 transmission chains, each with a different 
-genotype. The genotypes are subdivided into four serotypes. 500 randomly selected individuals from the underlying population 
+Dataset simulated using an agent based model with a spatially heterogeneous population structure. Infectious agents were
+introduced resulting in agent to agent transmission. The distance between successive cases in a transmission chain were
+randomly drawn from a uniform distribution U(0,100). Each infectious agent resulted in transmissions to two other agents
+after a delay of 15 days, reflecting the generation time of dengue. There are 11 transmission chains, each with a different
+genotype. The genotypes are subdivided into four serotypes. 500 randomly selected individuals from the underlying population
 also included.
 }
 \author{
diff --git a/man/get.pi.Rd b/man/get.pi.Rd
index f92336b..d60ea01 100644
--- a/man/get.pi.Rd
+++ b/man/get.pi.Rd
@@ -7,7 +7,7 @@
 get.pi(posmat, fun, r = 1, r.low = rep(0, length(r)), data.frame = TRUE)
 }
 \arguments{
-\item{posmat}{a matrix with columns x, y and any other named columns
+\item{posmat}{a matrix with columns x, y and any other named
 columns needed by \code{fun}}
 
 \item{fun}{a function that takes in two rows of \code{posmat} and returns:
@@ -15,11 +15,12 @@ columns needed by \code{fun}}
      \item for pairs included in the numerator and denominator
      \item for pairs that should only be included in the denominator
      \item for pairs that should be ignored all together}
-Note that names from \code{posmat} are not preserved in calls to \code{fun}, so the columns of the matrix should be
+Note that names from \code{posmat} are not preserved in calls to \code{fun}, so the columns of
+the matrix should be
 referenced numerically
 so this is not available to the \code{fun}}
 
-\item{r}{the series of spatial distances (or there maximums) we are
+\item{r}{the series of spatial distances (or their maximums) we are
 interested in}
 
 \item{r.low}{the low end of each range, 0 by default}
@@ -33,7 +34,8 @@ pi value for each distance range that we look at. Where:
 }
 \description{
 Generalized version of the \code{get.pi} function that takes in an arbitrary function and
-returns the probability that a point within a particular range of a point of interest shares the relationship
+returns the probability that a point within a particular range of a point of interest shares
+the relationship
 specified by the passed in function with that point.
 }
 \examples{
@@ -55,7 +57,7 @@ sero.pi<-get.pi(DengueSimR02,sero.type.func,r=r.max,r.low=r.min)
 }
 }
 \seealso{
-Other get.pi: 
+Other get.pi:
 \code{\link{get.pi.bootstrap}()},
 \code{\link{get.pi.ci}()},
 \code{\link{get.pi.permute}()},
@@ -63,9 +65,12 @@ Other get.pi:
 \code{\link{get.pi.typed.permute}()},
 \code{\link{get.pi.typed}()}
 
-Other spatialtau: 
+Other spatialtau:
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
 \code{\link{get.tau}()},
-\code{\link{get.theta}()}
+\code{\link{get.theta}()},
+\code{\link{plot.tau}()}
 }
 \author{
 Justin Lessler and Henrik Salje
diff --git a/man/get.pi.bootstrap.Rd b/man/get.pi.bootstrap.Rd
index 5beab5d..fdcbbe6 100644
--- a/man/get.pi.bootstrap.Rd
+++ b/man/get.pi.bootstrap.Rd
@@ -27,7 +27,9 @@ get.pi.bootstrap(
 \item{data.frame}{logical indicating whether to return results as a data frame (default = TRUE)}
 }
 \value{
-pi values for all the distances we looked at
+Values of pi for all distance bands. Return value dependent on data.frame argument.
+Asa matrix (rows = bootstrap samples, columns = increasing distance bands)
+or a data.frame (r.low, r and increasing distance bands)
 }
 \description{
 Runs \code{get.pi} on multiple bootstraps of the data. Is formulated
diff --git a/man/get.pi.ci.Rd b/man/get.pi.ci.Rd
index 4c47b41..efb1546 100644
--- a/man/get.pi.ci.Rd
+++ b/man/get.pi.ci.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/spatialfuncs.r
 \name{get.pi.ci}
 \alias{get.pi.ci}
-\title{Calculate bootstrapped confidence intervals for \code{get.pi} values.}
+\title{Calculate bootstrapped BCa confidence intervals from \code{get.pi} values.}
 \usage{
 get.pi.ci(
   posmat,
@@ -10,66 +10,74 @@ get.pi.ci(
   r = 1,
   r.low = rep(0, length(r)),
   boot.iter = 1000,
-  ci.low = 0.025,
-  ci.high = 0.975,
+  ci.level = 0.95,
   data.frame = TRUE
 )
 }
 \arguments{
-\item{posmat}{a matrix with columns type, x and y}
+\item{posmat}{a matrix with named columns x and y for 2D individual location}
 
-\item{fun}{the function to decide relationships}
+\item{fun}{the function to decide transmission-related pairs}
 
-\item{r}{the series of spatial distances wer are interested in}
+\item{r}{the upper end of each distance band}
 
-\item{r.low}{the low end of each range. 0 by default}
+\item{r.low}{the low end of each distance band (default: a vector of zeroes)}
 
-\item{boot.iter}{the number of bootstrap iterations}
+\item{boot.iter}{the number of bootstrap iterations (default = 1000)}
 
-\item{ci.low}{the low end of the ci...0.025 by default}
+\item{ci.level}{the level of the desired BCa CI (default = 0.95)}
 
-\item{ci.high}{the high end of the ci...0.975 by default}
-
-\item{data.frame}{logical indicating whether to return results as a data frame (default = TRUE)}
+\item{data.frame}{logical: indicating whether to return results as a data frame (default = TRUE)}
 }
 \value{
-a matrix with a row for the high and low values and
-    a column per distance
+If \code{data.frame = TRUE} then a data frame of 5 variables \code{r.low}, \code{r}, \code{pt.est} (the point estimate from \code{get.pi}), the confidence envelope as \code{ci.low} and \code{ci.high}, with the observations representing ascending distance bands. Else a matrix with first row \code{ci.low} and second row \code{ci.high} with columns representing ascending distance bands.
 }
 \description{
-Wrapper to \code{get.pi.bootstrap} that takes care of calculating the
-confidence intervals based on the bootstrapped values..
+Wrapper using \pkg{coxed} package to calculate the
+BCa (bias-corrected and accelerated) confidence interval (CI) for \eqn{\pi}(\code{r.low}, \code{r}), based on bootstrapped values from \code{get.pi.bootstrap}.
 }
-\examples{
-\donttest{
+\section{Depends on}{
 
-#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)
+coxed::bca()
 }
 
-r.max<-seq(10,100,10)
-r.min<-seq(0,90,10)
-r.mid <- (r.max+r.min)/2
+\examples{
+\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")
 
+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
+}
 
-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)
+# define distance band set
+r.max <- seq(10,100,10)
+r.min <- seq(0,90,10)
+r.mid <- (r.max + r.min)/2
 
-plot(r.mid, pi$pi, type="l")
-lines(r.mid, pi.ci[,2] , lty=2)
-lines(r.mid, pi.ci[,3] , lty=2)
+# 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$ci.low, lty=2)
+lines(r.mid, pi.ci$ci.high, lty=2)}
 }
+\references{
+\href{https://arxiv.org/pdf/1911.08022v4.pdf#page=12}{Rationale for BCa rather than percentile CIs} is described in Pollington et al. (2020)
+Developments in statistical inference when assessing
+spatiotemporal disease clustering with the tau statistic.
+\emph{arXiv/stat.ME: 1911.08022v4}.
 }
 \seealso{
-Other get.pi: 
+Other get.pi:
 \code{\link{get.pi.bootstrap}()},
 \code{\link{get.pi.permute}()},
 \code{\link{get.pi.typed.bootstrap}()},
@@ -78,6 +86,6 @@ Other get.pi:
 \code{\link{get.pi}()}
 }
 \author{
-Justin Lessler
+Justin Lessler and Timothy M Pollington
 }
 \concept{get.pi}
diff --git a/man/get.tau.D.param.est.Rd b/man/get.tau.D.param.est.Rd
new file mode 100644
index 0000000..43b3555
--- /dev/null
+++ b/man/get.tau.D.param.est.Rd
@@ -0,0 +1,96 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/spatialfuncs.r
+\name{get.tau.D.param.est}
+\alias{get.tau.D.param.est}
+\title{Cluster range estimation using \code{get.tau.D.param.est}}
+\usage{
+get.tau.D.param.est(r, tausim, GETres = NULL)
+}
+\arguments{
+\item{r}{the series of spatial distances (or their maximums) we are
+interested in}
+
+\item{tausim}{the set of spatially-bootstrapped simulations. Has to be \code{taubstrap} class; use \code{get.tau.bootstrap(..., data.frame = FALSE)} to obtain this.}
+
+\item{GETres}{is a required object and is obtained from a previous global hypothesis test using \code{get.tau.GET}. It ensures that the user has performed a graphical hypothesis test first and has considered there is evidence against H_0, before deciding to estimate the clustering range.}
+}
+\value{
+An object of class \code{tauparamest} which can then be plotted using \code{plot.tau()}. The object consists of:
+\itemize{
+  \item envelope the distribution of clustering range estimates
+}
+}
+\description{
+Estimates the range of spatiotemporal clustering. It records the place on the horizontal tau=1 line where each spatially bootstrapped simulation touches. This distribution then represents an empirical distribution for the clustering range and a confidence interval can be computed.
+}
+\section{Attributes}{
+
+\itemize{
+  \item BCaCI the BCa CI for the distribution of clustering range estimates
+}
+}
+
+\examples{
+\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)
+}
+}
+\seealso{
+Other get.tau: 
+\code{\link{get.tau.GET}()},
+\code{\link{get.tau.bootstrap}()},
+\code{\link{get.tau.ci}()},
+\code{\link{get.tau.permute}()},
+\code{\link{get.tau.typed.bootstrap}()},
+\code{\link{get.tau.typed.permute}()},
+\code{\link{get.tau.typed}()},
+\code{\link{get.tau}()},
+\code{\link{plot.tau}()}
+
+Other spatialtau: 
+\code{\link{get.pi}()},
+\code{\link{get.tau.GET}()},
+\code{\link{get.tau}()},
+\code{\link{get.theta}()},
+\code{\link{plot.tau}()}
+}
+\author{
+Timothy M Pollington
+}
+\concept{get.tau}
+\concept{spatialtau}
diff --git a/man/get.tau.GET.Rd b/man/get.tau.GET.Rd
new file mode 100644
index 0000000..e63f469
--- /dev/null
+++ b/man/get.tau.GET.Rd
@@ -0,0 +1,110 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/spatialfuncs.r
+\name{get.tau.GET}
+\alias{get.tau.GET}
+\title{Global hypothesis testing for the tau statistic}
+\usage{
+get.tau.GET(posmat, fun, r, r.low, permutations = 2500, comparison.type)
+}
+\arguments{
+\item{posmat}{a matrix with columns x, y and any other named
+columns needed by \code{fun}}
+
+\item{fun}{a function that takes in two rows of posmat and returns:
+\enumerate{
+     \item for pairs included in the numerator (and the denominator for independent data)
+     \item for pairs that should only be included in the denominator
+     \item for pairs that should be ignored all together}
+Note that names from \code{posmat} are not preserved in calls to
+\code{fun}, so the columns of the matrix should be referenced numerically
+so this is not available to fun}
+
+\item{r}{the series of spatial distances (or their maximums) we are
+interested in}
+
+\item{r.low}{the low end of each range, 0 by default}
+
+\item{permutations}{number of simulations of H_0. 2,500 is an optimal number according to Myllymäki et al. (2017).}
+
+\item{comparison.type}{what type of points are included in the comparison set.
+\itemize{
+  \item "representative" if comparison set is representative of the underlying population
+  \item "independent" if comparison set is cases/events coming from an indepedent process
+}}
+}
+\value{
+An object of class \code{tauGET} which can then be plotted using \code{plot.tau()} and an additional \code{tau} class object. The object consists of:
+\itemize{
+  \item r that inputted earlier
+  \item obs the tau point estimate computed internally using \code{get.tau()}
+  \item central the median estimate of all simulation curves that represent the null hypothesis. Comparing this to the tau=1 line indicates if it is reasonable to assume that H_0 was adequately simulated.
+  \item lo the lower bound of the global envelope
+  \item hi the upper bound of the global envelope
+  \item tau.permute the entire record of simulations of H_0, to plot with the global envelope using \code{plot.tau()}.
+}
+}
+\description{
+Performs a graphical hypothesis test to assess the evidence against the null hypothesis (H_0: tau = 1 i.e. no spatiotemporal clustering nor inhibition). A global envelope test from the \code{GET} package is used to see if any part of the point estimate connected line is outside the lower or upper bounds of the global envelope. The global envelope is formed on the tau estimator acting on time-permuted data to simulate H_0. The global envelope test is of 'extreme rank type' i.e. minimum of pointwise ranks with 95\% significance level.
+}
+\section{Attributes}{
+
+\itemize{
+  \item p_interval represents a range rather than a single p-value to assess the evidence against H_0. Accessed using \code{attr(x,"p_interval")}.
+}
+}
+
+\examples{
+\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)
+}
+}
+\seealso{
+Other get.tau: 
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.bootstrap}()},
+\code{\link{get.tau.ci}()},
+\code{\link{get.tau.permute}()},
+\code{\link{get.tau.typed.bootstrap}()},
+\code{\link{get.tau.typed.permute}()},
+\code{\link{get.tau.typed}()},
+\code{\link{get.tau}()},
+\code{\link{plot.tau}()}
+
+Other spatialtau: 
+\code{\link{get.pi}()},
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau}()},
+\code{\link{get.theta}()},
+\code{\link{plot.tau}()}
+}
+\author{
+Timothy M Pollington
+}
+\concept{get.tau}
+\concept{spatialtau}
diff --git a/man/get.tau.Rd b/man/get.tau.Rd
index 8f77683..e72e122 100644
--- a/man/get.tau.Rd
+++ b/man/get.tau.Rd
@@ -14,19 +14,19 @@ get.tau(
 )
 }
 \arguments{
-\item{posmat}{a matrix with columns x, y and any other named columns
-columns needed by fun}
+\item{posmat}{a matrix with columns x, y and any other named
+columns needed by \code{fun}}
 
 \item{fun}{a function that takes in two rows of posmat and returns:
 \enumerate{
-     \item for pairs included in the numerator (and the denominator for independent data)
+     \item for pairs included in the numerator (and the denominator for representative data)
      \item for pairs that should only be included in the denominator
      \item for pairs that should be ignored all together}
 Note that names from \code{posmat} are not preserved in calls to
 \code{fun}, so the columns of the matrix should be referenced numerically
 so this is not available to fun}
 
-\item{r}{the series of spatial distances (or there maximums) we are
+\item{r}{the series of spatial distances (or their maximums) we are
 interested in}
 
 \item{r.low}{the low end of each range, 0 by default}
@@ -37,10 +37,10 @@ interested in}
   \item "independent" if comparison set is cases/events coming from an indepedent process
 }}
 
-\item{data.frame}{logical indicating whether to return results as a data frame (default = TRUE)}
+\item{data.frame}{logical indicating whether to return results 'like' a data frame format (default = TRUE)}
 }
 \value{
-The tau value for each distance we look at. If \code{comparison.type} is "representative", this is:
+The tau value for each distance we look at as a tau class with a matrix or data frame style. If \code{comparison.type} is "representative", this is:
 
 \code{tau = get.pi(posmat, fun, r, r.low)/get.pi(posmat,fun,infinity,0)}
 
@@ -63,39 +63,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)
@@ -113,7 +98,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)
@@ -125,23 +116,47 @@ 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)
 }
 }
 \seealso{
-Other get.tau: 
+Other get.tau:
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
 \code{\link{get.tau.bootstrap}()},
 \code{\link{get.tau.ci}()},
 \code{\link{get.tau.permute}()},
 \code{\link{get.tau.typed.bootstrap}()},
 \code{\link{get.tau.typed.permute}()},
-\code{\link{get.tau.typed}()}
+\code{\link{get.tau.typed}()},
+\code{\link{plot.tau}()}
 
-Other spatialtau: 
+Other spatialtau:
 \code{\link{get.pi}()},
-\code{\link{get.theta}()}
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
+\code{\link{get.theta}()},
+\code{\link{plot.tau}()}
 }
 \author{
-Justin Lessler and Henrik Salje
+Justin Lessler, Timothy M Pollington and Henrik Salje
 }
 \concept{get.tau}
 \concept{spatialtau}
diff --git a/man/get.tau.bootstrap.Rd b/man/get.tau.bootstrap.Rd
index f0f77a4..9f67af7 100644
--- a/man/get.tau.bootstrap.Rd
+++ b/man/get.tau.bootstrap.Rd
@@ -40,7 +40,7 @@ calculated
 \examples{
 \donttest{
 
-#compare normally distributed with uniform points
+# 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")
@@ -55,26 +55,32 @@ r.max<-seq(10,100,10)
 r.min<-seq(0,90,10)
 r.mid <- (r.max+r.min)/2
 
-tau<-get.tau(x,fun,r=r.max,r.low=r.min)
-tau.boot<-get.tau.bootstrap(x,fun,r=r.max,r.low=r.min,boot.iter=50)
+tau<-get.tau(x,fun,r=r.max,r.low=r.min,"representative", data.frame = TRUE)
+tau.ci = get.tau.ci(x, fun, r.max, r.min, 50, "representative", 0.95, data.frame = TRUE)
 
-tau.ci<-apply(tau.boot[,-(1:2)],1,quantile,probs=c(0.25,0.75))
+## plot.tau() method
+plot.tau(tau, r.mid = TRUE, ptwise.CI = tau.ci)
 
-plot(r.mid, tau$tau ,ylim=c(min(tau.ci),max(tau.ci)), type="l", log="y")
+## previous plot() method using connected lines to join the top and bottoms of the pointwise CIs.
+#This may lead the user to perform graphical hypothesis testing using this plot without considering
+#the specific distance band of interest before plotting.
+plot(r.mid, tau$tau.pt.est ,ylim=c(min(tau.ci$ci.low),max(tau.ci$ci.high)), type="l", log="y")
 lines(c(0,100),c(1,1), lty=3, col="grey")
-lines(r.mid, tau.ci[1,] , lty=2)
-lines(r.mid, tau.ci[2,] , lty=2)
-
+lines(r.mid, tau.ci$ci.low, lty=2)
+lines(r.mid, tau.ci$ci.high, lty=2)
 }
 }
 \seealso{
-Other get.tau: 
+Other get.tau:
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
 \code{\link{get.tau.ci}()},
 \code{\link{get.tau.permute}()},
 \code{\link{get.tau.typed.bootstrap}()},
 \code{\link{get.tau.typed.permute}()},
 \code{\link{get.tau.typed}()},
-\code{\link{get.tau}()}
+\code{\link{get.tau}()},
+\code{\link{plot.tau}()}
 }
 \author{
 Justin Lessler and Henrik Salje
diff --git a/man/get.tau.ci.Rd b/man/get.tau.ci.Rd
index f6fbbcb..14c460f 100644
--- a/man/get.tau.ci.Rd
+++ b/man/get.tau.ci.Rd
@@ -11,8 +11,7 @@ get.tau.ci(
   r.low = rep(0, length(r)),
   boot.iter = 1000,
   comparison.type = "representative",
-  ci.low = 0.025,
-  ci.high = 0.975,
+  ci.level = 0.95,
   data.frame = TRUE
 )
 }
@@ -29,9 +28,7 @@ get.tau.ci(
 
 \item{comparison.type}{the comparison type to pass to get.tau}
 
-\item{ci.low}{the low end of the ci...0.025 by default}
-
-\item{ci.high}{the high end of the ci...0.975 by default}
+\item{ci.level}{significance level of the BCa CI, default = 0.95}
 
 \item{data.frame}{logical indicating whether to return results as a data frame (default = TRUE)}
 }
@@ -40,7 +37,7 @@ a data frame with the point estimate of tau and its low and high confidence inte
 }
 \description{
 Wrapper to \code{get.tau.bootstrap} that takes care of calulating
-the confidence intervals based on the bootstrapped values
+the confidence intervals based on the bootstrapped values.
 }
 \examples{
 \donttest{
@@ -60,24 +57,32 @@ r.max<-seq(10,100,10)
 r.min<-seq(0,90,10)
 r.mid <- (r.max+r.min)/2
 
-tau <- get.tau.ci(x,fun,r=r.max,r.low=r.min,boot.iter=50)
+tau.CI <- get.tau.ci(x,fun,r=r.max,r.low=r.min,boot.iter=50, comparison.type = "representative")
 
-plot(r.mid, tau$pt.est, ylim=c(1/max(tau[,3:5]), max(tau[,3:5])), type="l", log="y",
-     xlab="Distance", ylab="Tau")
-lines(r.mid, tau$ci.low , lty=2)
-lines(r.mid, tau$ci.high, lty=2)
-lines(c(0,100),c(1,1), lty=3, col="grey")
+## plot.tau() method
+tau = get.tau(x,fun,r=r.max,r.low=r.min, comparison.type = "representative")
+plot.tau(x = tau, ptwise.CI = tau.CI)
 
+## previous plot() method
+plot(r.mid, tau.CI$pt.est, ylim=c(min(tau.CI$pt.est,tau.CI$ci.low),
+                                  max(tau.CI$pt.est,tau.CI$ci.high)), type="l", xlab="Distance",
+     ylab="Tau")
+lines(r.mid, tau.CI$ci.low , lty=2)
+lines(r.mid, tau.CI$ci.high, lty=2)
+lines(c(0,100),c(1,1), lty=3, col="grey")
 }
 }
 \seealso{
-Other get.tau: 
+Other get.tau:
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
 \code{\link{get.tau.bootstrap}()},
 \code{\link{get.tau.permute}()},
 \code{\link{get.tau.typed.bootstrap}()},
 \code{\link{get.tau.typed.permute}()},
 \code{\link{get.tau.typed}()},
-\code{\link{get.tau}()}
+\code{\link{get.tau}()},
+\code{\link{plot.tau}()}
 }
 \author{
 Justin Lessler and Henrik Salje
diff --git a/man/get.tau.permute.Rd b/man/get.tau.permute.Rd
index 321cfac..7f0cce1 100644
--- a/man/get.tau.permute.Rd
+++ b/man/get.tau.permute.Rd
@@ -60,21 +60,26 @@ tau.null<-get.tau.permute(x,fun,r=r.max,r.low=r.min,permutations=50,comparison.t
 
 null.ci<-apply(tau.null[,-(1:2)],1,quantile,probs=c(0.25,0.75))
 
+# note these plots are only for illustrative purposes to show how get.tau.permute() can generate
+# the null distribution. These should not be used for graphical hypothesis testing nor parameter
+# estimation of the clustering endpoint.
 plot(r.mid, tau$tau, ylim=c(1/max(tau$tau),max(tau$tau)), type="l", log="y")
 lines(c(0,100),c(1,1), lty=3, col="grey")
 lines(r.mid, null.ci[1,] , lty=2)
 lines(r.mid, null.ci[2,] , lty=2)
-
 }
 }
 \seealso{
-Other get.tau: 
+Other get.tau:
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
 \code{\link{get.tau.bootstrap}()},
 \code{\link{get.tau.ci}()},
 \code{\link{get.tau.typed.bootstrap}()},
 \code{\link{get.tau.typed.permute}()},
 \code{\link{get.tau.typed}()},
-\code{\link{get.tau}()}
+\code{\link{get.tau}()},
+\code{\link{plot.tau}()}
 }
 \author{
 Justin Lessler and Henrik Salje
diff --git a/man/get.tau.typed.Rd b/man/get.tau.typed.Rd
index 0cd1079..3a2c2c8 100644
--- a/man/get.tau.typed.Rd
+++ b/man/get.tau.typed.Rd
@@ -63,13 +63,16 @@ abline(h=1,lty=2)
 }
 }
 \seealso{
-Other get.tau: 
+Other get.tau:
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
 \code{\link{get.tau.bootstrap}()},
 \code{\link{get.tau.ci}()},
 \code{\link{get.tau.permute}()},
 \code{\link{get.tau.typed.bootstrap}()},
 \code{\link{get.tau.typed.permute}()},
-\code{\link{get.tau}()}
+\code{\link{get.tau}()},
+\code{\link{plot.tau}()}
 }
 \author{
 Justin Lessler and Henrik Salje
diff --git a/man/get.tau.typed.bootstrap.Rd b/man/get.tau.typed.bootstrap.Rd
index cb350a8..6b1070d 100644
--- a/man/get.tau.typed.bootstrap.Rd
+++ b/man/get.tau.typed.bootstrap.Rd
@@ -51,16 +51,16 @@ r.max<-seq(20,1000,20)
 r.min<-seq(0,980,20)
 r.mid<-(r.max+r.min)/2
 
-# Lets see if there's a difference in spatial dependence between those that occurred 
+# Lets see if there's a difference in spatial dependence between those that occurred
 # late versus early in the outbreak
 
 type <- 2 - (DengueSimR02[,"time"] < 120)
 tmp <- cbind(DengueSimR02, type=type)
 
-typed.tau <- get.tau.typed(tmp, typeA=1, typeB=2, r=r.max, r.low=r.min, 
+typed.tau <- get.tau.typed(tmp, typeA=1, typeB=2, r=r.max, r.low=r.min,
                            comparison.type = "independent")
 
-typed.tau.type.bs <- get.tau.typed.bootstrap(tmp, typeA=1, typeB=2, r=r.max, r.low=r.min, 
+typed.tau.type.bs <- get.tau.typed.bootstrap(tmp, typeA=1, typeB=2, r=r.max, r.low=r.min,
                                              boot.iter=100, comparison.type = "independent")
 
 ci <- apply(typed.tau.type.bs[,-(1:2)], 1, quantile, probs=c(0.025,0.975))
@@ -77,13 +77,16 @@ lines(r.mid, ci[2,] , lty=2)
 }
 }
 \seealso{
-Other get.tau: 
+Other get.tau:
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
 \code{\link{get.tau.bootstrap}()},
 \code{\link{get.tau.ci}()},
 \code{\link{get.tau.permute}()},
 \code{\link{get.tau.typed.permute}()},
 \code{\link{get.tau.typed}()},
-\code{\link{get.tau}()}
+\code{\link{get.tau}()},
+\code{\link{plot.tau}()}
 }
 \author{
 Justin Lessler and Henrik Salje
diff --git a/man/get.tau.typed.permute.Rd b/man/get.tau.typed.permute.Rd
index 7a9d907..4b4cdce 100644
--- a/man/get.tau.typed.permute.Rd
+++ b/man/get.tau.typed.permute.Rd
@@ -55,15 +55,15 @@ r.mid<-(r.max+r.min)/2
 type <- 2 - (DengueSimR02[,"time"] < 120)
 tmp <- cbind(DengueSimR02, type=type)
 
-typed.tau <- get.tau.typed(tmp, typeA=1, typeB=2, r=r.max, r.low=r.min, 
+typed.tau <- get.tau.typed(tmp, typeA=1, typeB=2, r=r.max, r.low=r.min,
                            comparison.type = "independent")
 
-typed.tau.type.null<-get.tau.typed.permute(tmp, typeA=1, typeB=2, r=r.max, r.low=r.min, 
+typed.tau.type.null<-get.tau.typed.permute(tmp, typeA=1, typeB=2, r=r.max, r.low=r.min,
                                            permutations=100, comparison.type = "independent")
 
 null.ci <- apply(typed.tau.type.null[,-(1:2)], 1, quantile, probs=c(0.025,0.975))
 
-plot(r.mid, typed.tau$tau, ylim=c(0.3,4), log="y", cex.axis=1.25, 
+plot(r.mid, typed.tau$tau, ylim=c(0.3,4), log="y", cex.axis=1.25,
      xlab="Distance (m)", ylab="Tau", cex.main=0.9, lwd=2, type="n")
 abline(h=1,lty=1)
 lines(r.mid,typed.tau$tau,pch=20,col=1,lwd=3)
@@ -73,13 +73,16 @@ lines(r.mid, null.ci[2,] , lty=2)
 }
 }
 \seealso{
-Other get.tau: 
+Other get.tau:
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
 \code{\link{get.tau.bootstrap}()},
 \code{\link{get.tau.ci}()},
 \code{\link{get.tau.permute}()},
 \code{\link{get.tau.typed.bootstrap}()},
 \code{\link{get.tau.typed}()},
-\code{\link{get.tau}()}
+\code{\link{get.tau}()},
+\code{\link{plot.tau}()}
 }
 \author{
 Justin Lessler and Henrik Salje
diff --git a/man/get.theta.Rd b/man/get.theta.Rd
index cbea253..fb96490 100644
--- a/man/get.theta.Rd
+++ b/man/get.theta.Rd
@@ -7,8 +7,8 @@
 get.theta(posmat, fun, r = 1, r.low = rep(0, length(r)), data.frame = TRUE)
 }
 \arguments{
-\item{posmat}{a matrix with columns x, y and any other named columns
-columns needed by fun}
+\item{posmat}{a matrix with columns x, y and any other named
+columns needed by \code{fun}}
 
 \item{fun}{a function that takes in two rows of posmat and returns:
 \enumerate{
@@ -19,7 +19,7 @@ Note that names from \code{posmat} are not preserved in calls to \code{fun}, so
 referenced numerically
 so this is not available to the fun}
 
-\item{r}{the series of spatial distances (or there maximums) we are
+\item{r}{the series of spatial distances (or their maximums) we are
 interested in}
 
 \item{r.low}{the low end of each range, 0 by default}
@@ -55,7 +55,7 @@ sero.theta<-get.theta(DengueSimR02,sero.type.func,r=r.max,r.low=r.min)
 }
 }
 \seealso{
-Other get.theta: 
+Other get.theta:
 \code{\link{get.theta.bootstrap}()},
 \code{\link{get.theta.ci}()},
 \code{\link{get.theta.permute}()},
@@ -63,9 +63,12 @@ Other get.theta:
 \code{\link{get.theta.typed.permute}()},
 \code{\link{get.theta.typed}()}
 
-Other spatialtau: 
+Other spatialtau:
 \code{\link{get.pi}()},
-\code{\link{get.tau}()}
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
+\code{\link{get.tau}()},
+\code{\link{plot.tau}()}
 }
 \author{
 Justin Lessler and Henrik Salje
diff --git a/man/get.theta.ci.Rd b/man/get.theta.ci.Rd
index ac8f572..132a6dd 100644
--- a/man/get.theta.ci.Rd
+++ b/man/get.theta.ci.Rd
@@ -10,8 +10,7 @@ get.theta.ci(
   r = 1,
   r.low = rep(0, length(r)),
   boot.iter = 1000,
-  ci.low = 0.025,
-  ci.high = 0.975,
+  ci.level = 0.95,
   data.frame = TRUE
 )
 }
@@ -26,19 +25,16 @@ get.theta.ci(
 
 \item{boot.iter}{the number of bootstrap iterations}
 
-\item{ci.low}{the low end of the ci...0.025 by default}
-
-\item{ci.high}{the high end of the ci...0.975 by default}
+\item{ci.level}{significance level of the 95% BCa CI, default = 0.95}
 
 \item{data.frame}{logical indicating whether to return results as a data frame (default = TRUE)}
 }
 \value{
-a matrix with a row for the high and low values and
-    a column per distance
+a matrix with a row for the high and low values and a column per distance
 }
 \description{
 Wrapper to \code{get.theta.bootstrap} that takes care of calculating the
-confience intervals based on the bootstrapped values.
+confidence intervals based on the bootstrapped values.
 }
 \examples{
 \donttest{
@@ -68,7 +64,7 @@ lines(r.mid, theta.ci[,3] , lty=2)
 }
 }
 \seealso{
-Other get.theta: 
+Other get.theta:
 \code{\link{get.theta.bootstrap}()},
 \code{\link{get.theta.permute}()},
 \code{\link{get.theta.typed.bootstrap}()},
diff --git a/man/plot.tau.Rd b/man/plot.tau.Rd
new file mode 100644
index 0000000..48a32a6
--- /dev/null
+++ b/man/plot.tau.Rd
@@ -0,0 +1,63 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/spatialfuncs.r
+\name{plot.tau}
+\alias{plot.tau}
+\title{Plotting the results from tau functions}
+\usage{
+\method{plot}{tau}(
+  x,
+  r.mid = TRUE,
+  tausim = NULL,
+  ptwise.CI = NULL,
+  GET.res = NULL,
+  d.param.est = NULL,
+  ...
+)
+}
+\arguments{
+\item{x}{\code{tau} object; create using \code{get.tau(..., data.frame = TRUE)}. Required for all plots.}
+
+\item{r.mid}{If \code{TRUE}(default) then for each point the x-coordinate of the midpoint of a distance band is plotted and if \code{FALSE} the endpoint of the distance band is plotted.}
+
+\item{tausim}{the set of spatially-bootstrapped simulations of \code{taubstrap} class; use \code{get.tau.bootstrap()} to obtain this. Required for Estimation of the clustering range plot.}
+
+\item{ptwise.CI}{the set of pointwise CIs of \code{tauCI} class; create using \code{get.tau(..., data.frame = TRUE)}. Optional for the diagnostic plot but should not be supplied for the other plots.}
+
+\item{GET.res}{is a required object for the graphical hypothesis test plot but should not be supplied for the other plots. It is obtained from \code{get.tau.GET(..., data.frame = TRUE)}. It ensures that the user has performed a graphical hypothesis test first and has considered there is evidence against H_0, before deciding to estimate the clustering range.}
+
+\item{d.param.est}{a required object for Estimating the clustering range plot from \code{get.tau.D.param(..., data.frame = TRUE)}, but should not be supplied for the other plots. A \code{taubstrap} object will also be necessary.}
+
+\item{...}{other arguments which are standard for \code{plot()} for plot customisation}
+}
+\description{
+Three types of plots:
+\enumerate{
+\item Diagnostic plot to indicate the structure or magnitude of spatiotemporal clustering. Requires \code{tau} object; \code{tauCI} object optional to draw pointwise CIs. This plot is only suitable for the purpose of a graphical hypothesis test in the situation that a specific distance band is selected prior to graph creation.
+\item Graphical hypothesis test to assess the evidence against the null hypothesis (no spatiotemporal clustering nor inhibition). Requires \code{tau} and \code{tauGET} objects.
+\item Estimation of the clustering range (the distribution of the places on the horizontal tau=1 line, where decreasing bootstrap simulations first intercept). Requires \code{tau}, \code{tauparamest} and \code{taubstrap} objects.
+}
+}
+\seealso{
+Other get.tau: 
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
+\code{\link{get.tau.bootstrap}()},
+\code{\link{get.tau.ci}()},
+\code{\link{get.tau.permute}()},
+\code{\link{get.tau.typed.bootstrap}()},
+\code{\link{get.tau.typed.permute}()},
+\code{\link{get.tau.typed}()},
+\code{\link{get.tau}()}
+
+Other spatialtau: 
+\code{\link{get.pi}()},
+\code{\link{get.tau.D.param.est}()},
+\code{\link{get.tau.GET}()},
+\code{\link{get.tau}()},
+\code{\link{get.theta}()}
+}
+\author{
+Timothy M Pollington
+}
+\concept{get.tau}
+\concept{spatialtau}
diff --git a/tests/testthat/test-getpi.r b/tests/testthat/test-getpi.r
new file mode 100644
index 0000000..77f07d7
--- /dev/null
+++ b/tests/testthat/test-getpi.r
@@ -0,0 +1,143 @@
+test_that("get.pi returns 1 when labels are ignored", {
+
+    #generate a set of 100 random points even labeled between the two
+    x<-cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
+
+    colnames(x) <-c("type","x","y")
+
+    test <- function(a,b) {return(1)}
+
+    #with no lower limit
+    res <- get.pi(x,test,seq(10,100,10))$pi
+    expect_that(res,equals(rep(1,10)))
+
+    #with lower and upper limit
+    res <- get.pi(x,test,seq(10,100,10), seq(0,90,10))$pi
+    expect_that(res,equals(rep(1,10)))
+})
+
+test_that("get.pi returns appropriate values for cannonical test case 1 (equilateral triangle)", {
+
+    x <- rbind(c(1,0,0), c(1,1,0),c(2,.5,sqrt(.75)))
+    colnames(x) <-c("type","x","y")
+
+    test <- function(a,b) {
+        if (a[1] != 1) return(3)
+        if (b[1] == 2) return(1)
+        return(2)
+    }
+
+    #first no lower limit
+    res <- get.pi(x,test,1.5)$pi
+    res2 <- get.pi.typed(x,1,2,1.5)$pi
+
+    expect_that(res, equals(.5))
+    expect_that(res2, equals(.5))
+
+    #now with a lower limit
+
+    res <- get.pi(x,test,1.5,.5)$pi
+    res2 <- get.pi.typed(x,1,2,1.5,.5)$pi
+
+    expect_that(res, equals(.5))
+    expect_that(res2, equals(.5))
+
+})
+
+test_that("get.pi returns appropriate values cannonical test case 2 (points on a line)", {
+    x<-rbind(c(1,0,0), c(2,1,0), c(2,-1,0), c(3,2,0),
+             c(2,-2,0), c(3,3,0),c(3,-3,0))
+
+    colnames(x) <-c("type","x","y")
+
+    test <- function(a,b) {
+        if (a[1] != 1) return(3)
+        if (b[1] == 2) return(1)
+        return(2)
+    }
+
+    #pi 0,1.5 should be 1, 1.5-2.5 should be 0.5 and 2.5+ should be 0
+    res <- get.pi(x, test, c(1.5,2.5,Inf), c(0,1.5,2.5))$pi
+    res2 <- get.pi.typed(x, 1, 2, c(1.5,2.5,1000), c(0,1.5,2.5))$pi
+
+    expect_that(res,equals(c(1,0.5,0)))
+    expect_that(res2,equals(c(1,0.5,0)))
+
+})
+
+
+test_that("get.pi and get.pi.typed have same behavior on random data", {
+
+    #generate a set of 1000 random points even labeled between the two
+    x<-cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
+
+    colnames(x) <-c("type","x","y")
+
+    test <- function(a,b) {
+        if (a[1] != 1) return(3)
+        if (b[1] == 1) return(1)
+        return(2)
+    }
+
+    #no lower limit
+    res1 <- get.pi(x,test,seq(10,100,10))$pi
+    res2 <- get.pi.typed(x, 1,1, seq(10,100,10))$pi
+    expect_that(res1,equals(res2))
+
+    #lower limit
+    res1 <- get.pi(x,test,seq(10,100,10), seq(0,90,10))$pi
+    res2 <- get.pi.typed(x, 1,1, seq(10,100,10), seq(0,90,10))$pi
+    expect_that(res1,equals(res2))
+})
+
+test_that("get.pi returns identical results regardless of column order",
+          {
+              x<-cbind(rep(c(1,2),50), x=runif(100,0,100),
+                       y=runif(100,0,100))
+
+              colnames(x) <-c("type","x","y")
+
+              test <- function(a,b) {
+                  if (a[1] != 1) return(3)
+                  if (b[1] == 1) return(1)
+                  return(2)
+              }
+
+              res1 <- get.pi(x,test,seq(10,100,10), seq(0,90,10))$pi
+
+              test <- function(a,b) {
+                  if (a[3] != 1) return(3)
+                  if (b[3] == 1) return(1)
+                  return(2)
+              }
+
+              res2 <- get.pi(x[,c(3,2,1)],test,seq(10,100,10), seq(0,90,10))$pi
+
+              test <- function(a,b) {
+                  if (a[2] != 1) return(3)
+                  if (b[2] == 1) return(1)
+                  return(2)
+              }
+
+              res3 <- get.pi(x[,c(2,1,3)],test,seq(10,100,10), seq(0,90,10))$pi
+
+              expect_that(res1, equals(res2))
+              expect_that(res2, equals(res3))
+
+          })
+
+
+test_that ("get.pi fails nicely if x and y column names are not provided", {
+    x<-cbind(rep(c(1,2),500), a=runif(1000,0,100), b=runif(1000,0,100))
+
+    test <- function(a,b) {
+        if (a[1] != 2) return(3)
+        if (b[1] == 3) return(1)
+        return(2)
+    }
+
+    expect_that(get.pi(x,test,seq(10,50,10), seq(0,40,10))$pi,
+                throws_error("unique x and y columns must be defined"))
+
+})
+
diff --git a/inst/tests/test-getpibootstrap.r b/tests/testthat/test-getpibootstrap.r
similarity index 59%
rename from inst/tests/test-getpibootstrap.r
rename to tests/testthat/test-getpibootstrap.r
index 029f8cd..2f13e19 100644
--- a/inst/tests/test-getpibootstrap.r
+++ b/tests/testthat/test-getpibootstrap.r
@@ -1,26 +1,23 @@
-context("get.pi.bootstrap")
+test_that("get.pi.bootstrap runs and returns 1 when it should", {
 
-test_that("get.pi.boostrap runs and returns 1 when it should", {
+    x <- cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
 
-    x<-cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
-
-    colnames(x) <-c("type","x","y")
+    colnames(x) <- c("type","x","y")
 
     test <- function(a,b) {return(1)}
 
-    #should return a matrix of all ones
+    # should return a matrix of all ones
     res <- get.pi.bootstrap(x, test, seq(10,100,10), seq(0,90,10), 20)[,-(1:2)]
     expect_that(sum(res!=1),equals(0))
     expect_that(ncol(res),equals(20))
-
+    
 })
 
-
-test_that("get.pi.ci returns bootstrap cis when same seed", {
+test_that("get.pi.ci returns bootstrap CIs when same seed", {
      
-    x<-cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
+    x <- cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
 
-    colnames(x) <-c("type","x","y")
+    colnames(x) <- c("type","x","y")
 
     test <- function(a,b) {
         if (a[1] != 1) return(3)
@@ -32,27 +29,20 @@ test_that("get.pi.ci returns bootstrap cis when same seed", {
     res <- get.pi.bootstrap(x, test, seq(15,45,15), seq(0,30,15), 20)[,-(1:2)]
 
     set.seed(787)
-    ci1 <- get.pi.ci(x, test, seq(15,45,15), seq(0,30,15), 20)[,4:5]
+    ci1 <- get.pi.ci(x, test, seq(15,45,15), seq(0,30,15), 20, ci.level = 0.95)[,4:5]
 
-    expect_that(as.numeric(ci1[1,]),
-                equals(as.numeric(quantile(res[1,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+    expect_that(as.numeric(ci1[1,]), 
+                equals(coxed::bca(as.numeric(res[1,]),conf.level = 0.95)))
 
     expect_that(as.numeric(ci1[2,]),
-                equals(as.numeric(quantile(res[2,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+                equals(coxed::bca(as.numeric(res[2,]),conf.level = 0.95)))
 
     expect_that(as.numeric(ci1[3,]),
-                equals(as.numeric(quantile(res[3,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
-
+                equals(coxed::bca(as.numeric(res[3,]),conf.level = 0.95)))
+    
 })
 
 
-
 test_that("performs correctly for test case 1 (equilateral triangle)", {
     x <- rbind(c(1,0,0), c(1,1,0),c(2,.5,sqrt(.75)))
     colnames(x) <-c("type","x","y")
@@ -64,23 +54,23 @@ test_that("performs correctly for test case 1 (equilateral triangle)", {
     }
 
     res <- get.pi.bootstrap(x, test, 1.5, 0.1, 500)[,-(1:2)]
+    res = res[!is.na(res)]
     res2 <- get.pi.typed.bootstrap(x, 1,2, 1.5, 0.1, 500)[,-(1:2)]
-
+    res2 = res2[!is.na(res2)]
+    
     #should have 95% CI of 0,1 and mean/median of 0.5
-    expect_that(as.numeric(quantile(res[1,], probs=c(.025,.975), na.rm=T)),
+    expect_that(coxed::bca(res, conf.level = 0.95),
                 equals(c(0,1)))
-    expect_that(as.numeric(quantile(res2[1,], probs=c(.025,.975), na.rm=T)),
+    expect_that(coxed::bca(res2, conf.level = 0.95),
                 equals(c(0,1)))
-
-
 })
 
 test_that("performs correctly for test case 2 (points on a line)", {
 
-    x<-rbind(c(1,0,0), c(2,1,0), c(2,-1,0), c(3,2,0),
+    x <- rbind(c(1,0,0), c(2,1,0), c(2,-1,0), c(3,2,0),
              c(2,-2,0), c(3,3,0),c(3,-3,0))
 
-    colnames(x) <-c("type","x","y")
+    colnames(x) <- c("type","x","y")
 
     test <- function(a,b) {
         if (a[1] != 1) return(3)
@@ -100,28 +90,35 @@ test_that("performs correctly for test case 2 (points on a line)", {
     expect_that(median(as.numeric(res2[2,]), na.rm=T), equals(0.5))
     expect_that(median(as.numeric(res2[3,]), na.rm=T), equals(0))
 
-
+    # For the first and third ranges we use the quantile method as the BCa method breaks down for 
+    # constant vectors of ones or zeroes
+    
     #FIRST RANGE
     #deterministically 1
+    
     expect_that(as.numeric(quantile(res[1,], probs=c(.025,.975), na.rm=T)),
                 equals(c(1,1)))
     expect_that(as.numeric(quantile(res2[1,], probs=c(.025,.975), na.rm=T)),
                 equals(c(1,1)))
 
     #SECOND RANGE...should be 0 and 1 respectively a fairly large % of the time
-    expect_that(as.numeric(quantile(res[2,], probs=c(0.025,.975), na.rm=T)),
+    res1.2 = na.omit(as.numeric(res[2,]))
+    res2.2 = na.omit(as.numeric(res2[2,]))
+    
+    expect_that(coxed::bca(as.numeric(res1.2), conf.level = 0.95),
                 equals(c(0,1)))
-    expect_that(as.numeric(quantile(res2[2,], probs=c(.025,.975), na.rm=T)),
+    expect_that(coxed::bca(as.numeric(res2.2), conf.level = 0.95),
                 equals(c(0,1)))
+    
 
     #THIRD RANGE
     #deterministically 0
+    
     expect_that(as.numeric(quantile(res[3,], probs=c(.025,.975), na.rm=T)),
                 equals(c(0,0)))
     expect_that(as.numeric(quantile(res2[3,], probs=c(.025,.975), na.rm=T)),
                 equals(c(0,0)))
 
-
 })
 
 
@@ -142,31 +139,3 @@ test_that ("fails nicely if x and y column names are not provided", {
     expect_that(get.pi.ci(x,test,seq(10,50,10), seq(0,40,10),100),
                 throws_error("unique x and y columns must be defined"))
 })
-
-
-##################DEPRECATED TESTS...TAKE TO LONG...NOW USING SMALLER CANONICAL
-##################TESTS THAT HAVE VALUES THAT CAN BE WORKED OUT BY HAND
-
-
-
-## test_that("CIs calculated from get.pi.bootstrap include the true value", {
-##     set.seed(787)
-
-##     x<-cbind(rep(c(1,2),250), x=runif(500,0,100), y=runif(500,0,100))
-
-##     colnames(x) <-c("type","x","y")
-
-##     test <- function(a,b) {
-##         if (a[1] != 1) return(3)
-##         if (b[1] == 1) return(1)
-##         return(2)
-##     }
-
-##     res <- get.pi.ci(x, test, seq(10,100,10), seq(0,90,10), 100)
-
-##     expect_that(sum(!(res[1,]<res[2,])),equals(0))
-##     expect_that(sum(!(res[1,]<.5)),equals(0))
-##     expect_that(sum(!(res[2,]>.5)),equals(0))
-## })
-
-
diff --git a/tests/testthat/test-getpipermute.r b/tests/testthat/test-getpipermute.r
new file mode 100644
index 0000000..36b9e9b
--- /dev/null
+++ b/tests/testthat/test-getpipermute.r
@@ -0,0 +1,88 @@
+test_that("get.pi.permute returns appropriate values for test case 1 (equilateral triangle)" ,{
+
+    x <- rbind(c(1,0,0), c(1,1,0),c(2,.5,sqrt(.75)))
+    colnames(x) <-c("type","x","y")
+
+    test <- function(a,b) {
+        if (a[1] != 1) return(3)
+        if (b[1] == 2) return(1)
+        return(2)
+    }
+
+
+    #should return .5 for every permutation
+    res <- get.pi.permute(x, test, 1.5, 0, 500)[,-(1:2)]
+    res2 <- get.pi.typed.permute(x, 1, 2, 1.5, 0, 500)[,-(1:2)]
+
+    expect_that(as.numeric(res), equals(rep(0.5,500)))
+    expect_that(as.numeric(res2), equals(rep(0.5,500)))
+
+})
+
+test_that("get.pi.permute returns appropriate values for test case 2 (points on a line) with
+          windows" ,{
+    x<-rbind(c(1,0,0), c(2,1,0), c(2,-1,0), c(3,2,0),
+             c(2,-2,0), c(3,3,0),c(3,-3,0))
+
+    colnames(x) <-c("type","x","y")
+
+    test <- function(a,b) {
+        if (a[1] != 1) return(3)
+        if (b[1] == 2) return(1)
+        return(2)
+    }
+
+    #the mean of the null distribution should be 0.5
+    #the 95% CI equals 0,1 with windows
+    set.seed(seed = 1)
+    res <- get.pi.permute(x, test, c(1.5,2.5,3.5), c(0,1.5,2.5), 500)[,-(1:2)]
+    res2 <- get.pi.typed.permute(x, 1, 2, c(1.5,2.5,3.5), c(0,1.5,2.5), 500)[,-(1:2)]
+
+    expect_that(rowMeans(res,na.rm=T), equals(rep(.5,3), tolerance=0.1))
+    expect_that(rowMeans(res2, na.rm=T), equals(rep(.5,3), tolerance=0.1))
+
+    for (i in 1:3) {
+        expect_that(coxed::bca(as.numeric(res[i,]), conf.level = 0.95),
+                    equals(c(0,1)))
+        expect_that(coxed::bca(as.numeric(res2[i,]), conf.level = 0.95),
+                    equals(c(0,1)))
+    }
+
+})
+
+test_that("get.pi.permute returns appropriate values for test case 2 (points on a line),
+          no windows" ,{
+    x<-rbind(c(1,0,0), c(2,1,0), c(2,-1,0), c(3,2,0),
+             c(2,-2,0), c(3,3,0),c(3,-3,0))
+    
+    colnames(x) <-c("type","x","y")
+    
+    test <- function(a,b) {
+        if (a[1] != 1) return(3)
+        if (b[1] == 2) return(1)
+        return(2)
+    }
+    
+    #without windows the distributions is asymmetric and 95% BCa CI is rather than [0.25,0.75] for 
+    # percentile CIs
+    set.seed(seed = 1)
+    res3 <- get.pi.permute(x, test, 4,0, 500)
+    res4 <- get.pi.typed.permute(x, 1, 2, 4, 0, 500)
+    expect_that(coxed::bca(as.numeric(res3[1,]), conf.level = 0.95), equals(c(1/3,1)))
+    expect_that(coxed::bca(as.numeric(res4[1,]), conf.level = 0.95), equals(c(1/3,1)))
+})
+
+
+test_that ("fails nicely if x and y column names are not provided", {
+    x<-cbind(rep(c(1,2),500), a=runif(1000,0,100), b=runif(1000,0,100))
+
+    test <- function(a,b) {
+        if (a[1] != 2) return(3)
+        if (b[1] == 3) return(1)
+        return(2)
+    }
+
+    expect_that(get.pi.permute(x,test,seq(10,50,10), seq(0,40,10),100),
+                throws_error("unique x and y columns must be defined"))
+
+})
\ No newline at end of file
diff --git a/inst/tests/test-gettau.r b/tests/testthat/test-gettau.r
similarity index 60%
rename from inst/tests/test-gettau.r
rename to tests/testthat/test-gettau.r
index ec1c910..c0836f4 100644
--- a/inst/tests/test-gettau.r
+++ b/tests/testthat/test-gettau.r
@@ -1,5 +1,3 @@
-context("get.tau")
-
 test_that("get.tau returns 1 when labels are ignored", {
 
     #generate a set of 1000 random points even labeled between the two
@@ -11,7 +9,7 @@ test_that("get.tau returns 1 when labels are ignored", {
 
     #######FIRST WITH REPRESENTATIVE SAMPLE ASSUMED
     #with no lower limit
-    res <- get.tau(x,test,seq(10,100,10))$tau
+    res <- get.tau(x,test,seq(10,100,10))$tau.pt.est
     expect_that(res,equals(rep(1,10)))
 
     #with lower and upper limit
@@ -172,103 +170,8 @@ test_that ("selection of an invalid comparison type fails nicely", {
 
 
     expect_that( get.tau(x, test, c(1.5,2.5,Inf), c(0,1.5,2.5), comparison.type="foobar"),
-                throws_error("unkown comparison type specified"))
+                throws_error("unknown comparison.type specified"))
     expect_that( get.tau.typed(x, 1, 2, c(1.5,2.5,1000), c(0,1.5,2.5), comparison.type="foobar"),
-                throws_error("unkown comparison type specified"))
+                throws_error("unknown comparison.type specified"))
 })
 
-
-##################DEPRECATED TESTS...TAKE TO LONG...NOW USING SMALLER CANONICAL
-##################TESTS THAT HAVE VALUES THAT CAN BE WORKED OUT BY HAND
-
-
-## test_that("get.tau and get.tau.typed have same behavior and both return about 1 when they should", {
-##     set.seed(380)
-
-##     #generate a set of 1000 random points even labeled between the two
-##     x<-cbind(rep(c(1,2),500), x=runif(1000,0,100), y=runif(1000,0,100))
-
-##     colnames(x) <-c("type","x","y")
-
-##     test <- function(a,b) {
-##         if (a[1] != 1) return(3)
-##         if (b[1] == 1) return(1)
-##         return(2)
-##     }
-
-##     #no lower limit
-##     res1 <- get.tau(x,test,seq(10,100,10))
-##     res2 <- get.tau.typed(x, 1,1, seq(10,100,10))
-##     expect_that(res1,equals(res2))
-##     expect_that(round(res1,1),equals(rep(1,10)))
-##     expect_that(round(res2,1),equals(rep(1,10)))
-
-
-##     #lower limit
-##     res1 <- get.tau(x,test,seq(10,100,10), seq(0,90,10))
-##     res2 <- get.tau.typed(x, 1,1, seq(10,100,10), seq(0,90,10))
-##     expect_that(res1,equals(res2))
-##     expect_that(round(res1,1),equals(rep(1,10)))
-##     expect_that(round(res2,1),equals(rep(1,10)))
-## })
-
-## test_that("get.tau with equality test returns about .5 when points are uniform",
-##       {
-##           set.seed(787)
-
-##           #generate a set of 1000 random points even labeled between the two
-##           x<-cbind(rep(c(1,2),500), x=runif(1000,0,100), y=runif(1000,0,100))
-
-##           colnames(x) <-c("type","x","y")
-
-##           test <- function(a,b) {
-##               if (a[1]==b[1]) return(1)
-##               return(2)
-##           }
-
-##           #no lower limit
-##           res1 <- get.tau(x,test,seq(10,100,10))
-##           expect_that(round(res1,1),equals(rep(1,10)))
-
-##           #lower limit
-##           res1 <- get.tau(x,test,seq(10,100,10), seq(0,90,10))
-##           expect_that(round(res1,1),equals(rep(1,10)))
-##       })
-
-
-## test_that("get.tau is montonically decreasing for normally distributed clusters",
-##           {
-##               set.seed(787)
-
-##               #first generate 500 random uniform points
-##               x<-cbind(1, x=runif(500,0,100), y=runif(500,0,100))
-##               colnames(x) <-c("type","x","y")
-
-##               #add a seed point
-##               x<-rbind(x,c(2,50,50))
-
-##               #generate 500 normally distibuted points around this
-##               x<-rbind(x,cbind(3,rnorm(1000,50,20),rnorm(1000,50,20)))
-
-##               #check wit get.tau.typed
-##               res1 <- get.tau.typed(x,2,3,seq(10,50,10), seq(0,40,10))
-
-##               expect_that(res1[1]>res1[2] & res1[2]>res1[3] &
-##                           res1[3]>res1[4] & res1[4]>res1[5], is_true())
-
-
-##               #do
-##               test <- function(a,b) {
-##                   if (a[1] != 2) return(3)
-##                   if (b[1] == 3) return(1)
-##                   return(2)
-##               }
-
-##               res2 <- get.tau(x,test,seq(10,50,10), seq(0,40,10))
-
-
-##               expect_that(res2[1]>res2[2] & res2[2]>res2[3] &
-##                           res2[3]>res2[4] & res2[4]>res2[5], is_true())
-##               expect_that(res1,equals(res2))
-##           })
-
diff --git a/inst/tests/test-gettaubootstrap.r b/tests/testthat/test-gettaubootstrap.r
similarity index 55%
rename from inst/tests/test-gettaubootstrap.r
rename to tests/testthat/test-gettaubootstrap.r
index ba62690..983c250 100644
--- a/inst/tests/test-gettaubootstrap.r
+++ b/tests/testthat/test-gettaubootstrap.r
@@ -1,6 +1,4 @@
-context("get.tau.bootstrap")
-
-test_that("get.tau.bootstrap runs and returs 1 when it should", {
+test_that("get.tau.bootstrap runs and returns 1 when it should", {
 
     x<-cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
 
@@ -10,14 +8,15 @@ test_that("get.tau.bootstrap runs and returs 1 when it should", {
 
     ########### REPRESENTATIVE
     #should return a matrix of all ones
-    res <- get.tau.bootstrap(x, test, seq(10,100,10), seq(0,90,10), 10)[,-(1:2)]
+    set.seed(1)
+    res <- get.tau.bootstrap(x, test, seq(10,100,10), seq(0,90,10), 10, data.frame = FALSE)
     expect_that(sum(res!=1),equals(0))
     expect_that(nrow(res),equals(10))
 
 
     ########### INDEPENDENT
     res <- get.tau.bootstrap(x, test, seq(10,100,10), seq(0,90,10), 10,
-                             comparison.type="independent")[,-(1:2)]
+                             comparison.type="independent", data.frame = FALSE)
     expect_that(sum(res!=1),equals(0))
     expect_that(nrow(res),equals(10))
 
@@ -35,24 +34,27 @@ test_that("performs correctly for test case 1 (equilateral triangle)", {
     }
 
     ########### REPRESENTATIVE
-    res <- get.tau.bootstrap(x, test, 1.5, 0.1, 500)[,-(1:2)]
+    res <- get.tau.bootstrap(x, test, 1.5, 0.1, 500, data.frame = FALSE)
     res2 <- get.tau.typed.bootstrap(x, 1,2, 1.5, 0.1, 500)[,-(1:2)]
 
-    #should have 95% CI of 1,1
-    expect_that(as.numeric(quantile(res[1,], probs=c(.025,.975), na.rm=T)),
+    #should have 95% CI of 1,1. quantile() method used as coxed::bca() breaks
+    # down under Inf conditions
+    expect_that(as.numeric(quantile(res[,1], probs=c(.025,.975), na.rm=TRUE)),
                 equals(c(1,1)))
 
     expect_that(as.numeric(quantile(res2[1,], probs=c(.025,.975), na.rm=T)),
                 equals(c(1,1)))
 
     ########### INDEPENDENT
+    set.seed(1)
     res <- get.tau.bootstrap(x, test, 1.5, 0.1, 500,
-                             comparison.type="independent")[,-(1:2)]
+                             comparison.type="independent", data.frame = FALSE)
     res2 <- get.tau.typed.bootstrap(x, 1,2, 1.5, 0.1, 500,
                                     comparison.type="independent")[,-(1:2)]
 
-    #should have 95% CI of 1,1
-    expect_that(as.numeric(quantile(res[1,], probs=c(.025,.975), na.rm=T)),
+    #should have 95% CI of 1,1. quantile() method used as coxed::bca() breaks
+    # down under Inf conditions
+    expect_that(as.numeric(quantile(res[,1], probs=c(.025,.975), na.rm=T)),
                 equals(c(1,1)))
 
     expect_that(as.numeric(quantile(res2[1,], probs=c(.025,.975), na.rm=T)),
@@ -76,25 +78,25 @@ test_that("performs correctly for test case 2 (points on a line) - representativ
 
     ########### REPRESENTATIVE
     #the medians for the null distribution should be 2,1,0
-    res <- get.tau.bootstrap(x, test, c(1.5,2.5,3.5), c(0,1.5,2.5), 1500)[,-(1:2)]
+    set.seed(1)
+    res <- get.tau.bootstrap(x, test, c(1.5,2.5,3.5), c(0,1.5,2.5), 1500, data.frame = FALSE)
     res2 <- get.tau.typed.bootstrap(x, 1, 2, c(1.5,2.5,3.5), c(0,1.5,2.5), 1500)[,-(1:2)]
 
-    expect_that(median(as.numeric(res[1,]), na.rm=T), equals(2))
-    expect_that(median(as.numeric(res[2,]), na.rm=T), equals(1))
-    expect_that(median(as.numeric(res[3,]), na.rm=T), equals(0))
+    expect_that(median(as.numeric(res[,1]), na.rm=T), equals(2))
+    expect_that(median(as.numeric(res[,2]), na.rm=T), equals(1))
+    expect_that(median(as.numeric(res[,3]), na.rm=T), equals(0))
 
     expect_that(median(as.numeric(res2[1,]), na.rm=T), equals(2))
     expect_that(median(as.numeric(res2[2,]), na.rm=T), equals(1))
     expect_that(median(as.numeric(res2[3,]), na.rm=T), equals(0))
 
-
-
+    # quantile() used over coxed::bca() as latter breaks down under these toy conditions or cannot provide the interval required.
 
     #FIRST RANGE
     #max would be only 1 type 2 used and in range = 1/(1/6) = 6...should occur
     #more than 2.5% of time
     #min would be 1, occuring just over .01% of the time
-    expect_that(as.numeric(quantile(res[1,], probs=c(.001,.975), na.rm=T)),
+    expect_that(as.numeric(quantile(res[,1], probs=c(.025,.975), na.rm=T)),
                 equals(c(1,6)))
     expect_that(as.numeric(quantile(res2[1,], probs=c(.001,.975), na.rm=T)),
                 equals(c(1,6)))
@@ -102,21 +104,17 @@ test_that("performs correctly for test case 2 (points on a line) - representativ
     #SECOND RANGE
     #max would be 6, should occur less than 1% of the time
     #min should be 0, should occur 2.5% of the time
-    expect_that(as.numeric(quantile(res[2,], probs=c(.025), na.rm=T)),
+    expect_that(as.numeric(quantile(res[,2], probs=c(.025), na.rm=T)),
                 equals(0))
     expect_that(as.numeric(quantile(res2[2,], probs=c(.025), na.rm=T)),
                 equals(0))
 
-    expect_that(as.numeric(quantile(res[2,], probs=c(.99), na.rm=T))<6,
-                is_true())
-    expect_that(as.numeric(quantile(res2[2,], probs=c(.99), na.rm=T))<6,
-                is_true())
-
-
+    expect_true(as.numeric(quantile(res[,2], probs=c(.99), na.rm=T))<6)
+    expect_true(as.numeric(quantile(res2[2,], probs=c(.99), na.rm=T))<6)
 
     #THIRD RANGE
     #Should be determinsitically 0 or NaN
-    expect_that(as.numeric(quantile(res[3,], probs=c(.025,.975), na.rm=T)),
+    expect_that(as.numeric(quantile(res[,3], probs=c(.025,.975), na.rm=T)),
                 equals(c(0,0)))
     expect_that(as.numeric(quantile(res2[3,], probs=c(.025,.975), na.rm=T)),
                 equals(c(0,0)))
@@ -139,26 +137,28 @@ test_that("performs correctly for test case 2 (points on a line) - independent c
 
     ########### INDEPENDENT
     #the medians for the null distribution should be Inf,1,0
+    set.seed(1)
     res <- get.tau.bootstrap(x, test, c(1.5,2.5,3.5), c(0,1.5,2.5), 1500,
-                             comparison.type="independent")[,-(1:2)]
+                             comparison.type="independent", data.frame = FALSE)
     res2 <- get.tau.typed.bootstrap(x, 1, 2, c(1.5,2.5,3.5), c(0,1.5,2.5), 1500,
                                     comparison.type="independent")[,-(1:2)]
 
-    expect_that(median(as.numeric(res[1,]), na.rm=T), equals(Inf))
-    expect_that(median(as.numeric(res[2,]), na.rm=T), equals(1))
-    expect_that(median(as.numeric(res[3,]), na.rm=T), equals(0))
+    expect_that(median(as.numeric(res[,1]), na.rm=T), equals(Inf))
+    expect_that(median(as.numeric(res[,2]), na.rm=T), equals(1))
+    expect_that(median(as.numeric(res[,3]), na.rm=T), equals(0))
 
     expect_that(median(as.numeric(res2[1,]), na.rm=T), equals(Inf))
     expect_that(median(as.numeric(res2[2,]), na.rm=T), equals(1))
     expect_that(median(as.numeric(res2[3,]), na.rm=T), equals(0))
 
 
-
+    # quantile() used over coxed::bca() as latter breaks down under these toy conditions or cannot
+    # provide the interval required.
 
     #FIRST RANGE
-    #max would be Inf, occuring most of the time
-    #min would be 1, occuring just over .01% of the time
-    expect_that(as.numeric(quantile(res[1,], probs=c(.001,.975), na.rm=T)),
+    #max would be Inf, occurring most of the time
+    #min would be 1, occurring just over .01% of the time
+    expect_that(as.numeric(quantile(res[,1], probs=c(.001,.975), na.rm=T)),
                 equals(c(1,Inf)))
     expect_that(as.numeric(quantile(res2[1,], probs=c(.001,.975), na.rm=T)),
                 equals(c(1,Inf)))
@@ -167,21 +167,19 @@ test_that("performs correctly for test case 2 (points on a line) - independent c
     #max would be Inf, should occur around 25% of the time. .7 should be
     # reliably less than
     #min should be 0, should occur 2.5% of the time
-    expect_that(as.numeric(quantile(res[2,], probs=c(.025), na.rm=T)),
+    expect_that(as.numeric(quantile(res[,2], probs=c(.025), na.rm=T)),
                 equals(0))
     expect_that(as.numeric(quantile(res2[2,], probs=c(.025), na.rm=T)),
                 equals(0))
 
-    expect_that(as.numeric(quantile(res[2,], probs=c(.7), na.rm=T))!=Inf,
-                is_true())
-    expect_that(as.numeric(quantile(res2[2,], probs=c(.7), na.rm=T))!=Inf,
-                is_true())
+    expect_true(as.numeric(quantile(res[2,], probs=c(.7), na.rm=T))==Inf)
+    expect_true(as.numeric(quantile(res2[2,], probs=c(.7), na.rm=T))!=Inf)
 
 
 
     #THIRD RANGE
     #Should be determinsitically 0 or NaN
-    expect_that(as.numeric(quantile(res[3,], probs=c(.025,.975), na.rm=T)),
+    expect_that(as.numeric(quantile(res[,3], probs=c(.025,.975), na.rm=T)),
                 equals(c(0,0)))
     expect_that(as.numeric(quantile(res2[3,], probs=c(.025,.975), na.rm=T)),
                 equals(c(0,0)))
@@ -202,53 +200,41 @@ test_that("get.tau.ci returns bootstrap cis when same seed", {
 
     ####REPRESENTATIVE
     set.seed(787)
-    res <- get.tau.bootstrap(x, test, seq(15,45,15), seq(0,30,15), 20)[,-(1:2)]
+    res <- get.tau.bootstrap(x, test, seq(15,45,15), seq(0,30,15), 20, data.frame = FALSE)
 
     set.seed(787)
-    ci1 <- get.tau.ci(x, test, seq(15,45,15), seq(0,30,15), 20)[,-(1:3)]
+    ci1 <- get.tau.ci(x, test, seq(15,45,15), seq(0,30,15), 20, comparison.type = "representative",
+                      ci.level = 0.95, data.frame = FALSE)
 
-    expect_that(as.numeric(ci1[1,]),
-                equals(as.numeric(quantile(res[1,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+    expect_that(as.numeric(t(ci1)[1,]), 
+                equals(coxed::bca(as.numeric(res[,1]),conf.level = 0.95)))
 
-    expect_that(as.numeric(ci1[2,]),
-                equals(as.numeric(quantile(res[2,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+    expect_that(as.numeric(t(ci1)[2,]), 
+                equals(coxed::bca(as.numeric(res[,2]),conf.level = 0.95)))
 
-    expect_that(as.numeric(ci1[3,]),
-                equals(as.numeric(quantile(res[3,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+    expect_that(as.numeric(t(ci1)[3,]), 
+                equals(coxed::bca(as.numeric(res[,3]),conf.level = 0.95)))
 
     ### INDEPENDENT
     set.seed(787)
     res <- get.tau.bootstrap(x, test, seq(15,45,15), seq(0,30,15), 20,
-                             comparison.type="independent")[,-(1:2)]
+                             comparison.type="independent", data.frame = FALSE)
 
     set.seed(787)
-    ci1 <- get.tau.ci(x, test, seq(15,45,15), seq(0,30,15), 20,
-                      comparison.type="independent")[,-(1:3)]
+    ci1 <- get.tau.ci(x, test, seq(15,45,15), seq(0,30,15), 20, comparison.type="independent", 
+                      data.frame = FALSE)
 
-    expect_that(as.numeric(ci1[1,]),
-                equals(as.numeric(quantile(res[1,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+    expect_that(as.numeric(t(ci1)[1,]), 
+                equals(coxed::bca(as.numeric(res[,1]),conf.level = 0.95)))
 
-    expect_that(as.numeric(ci1[2,]),
-                equals(as.numeric(quantile(res[2,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+    expect_that(as.numeric(t(ci1)[2,]), 
+                equals(coxed::bca(as.numeric(res[,2]),conf.level = 0.95)))
 
-    expect_that(as.numeric(ci1[3,]),
-                equals(as.numeric(quantile(res[3,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+    expect_that(as.numeric(t(ci1)[3,]), 
+                equals(coxed::bca(as.numeric(res[,3]),conf.level = 0.95)))
 
 })
 
-
 test_that("fails nicely if x and y column names are not provided", {
 
     x<-cbind(rep(c(1,2),500), a=runif(1000,0,100), b=runif(1000,0,100))
@@ -266,47 +252,3 @@ test_that("fails nicely if x and y column names are not provided", {
                 throws_error("unique x and y columns must be defined"))
 })
 
-
-##################DEPRECATED TESTS...TAKE TO LONG...NOW USING SMALLER CANONICAL
-##################TESTS THAT HAVE VALUES THAT CAN BE WORKED OUT BY HAND
-
-
-## test_that("CIs calculated from get.tau.bootstrap include the true value", {
-##     set.seed(777)
-
-##     x<-cbind(rep(c(1,2),250), x=runif(500,0,100), y=runif(500,0,100))
-
-##     colnames(x) <-c("type","x","y")
-
-##     test <- function(a,b) {
-##         if (a[1] != 1) return(3)
-##         if (b[1] == 1) return(1)
-##         return(2)
-##     }
-
-##     res <- get.tau.ci(x, test, seq(10,100,10), seq(0,90,10), 200)
-
-##     #print(res)
-
-##     expect_that(sum(!(res[1,]<res[2,])),equals(0))
-##     expect_that(sum(!(res[1,]<1)),equals(0))
-##     expect_that(sum(!(res[2,]>1)),equals(0))
-
-##     #repeat for typed data
-##     res <- get.tau.typed.bootstrap(x, typeA=1, typeB=1,
-##                                    seq(10,100,10), seq(0,90,10), 200)
-
-##     ci <- matrix(nrow=2, ncol=ncol(res))
-
-##     for (i in 1:ncol(ci)) {
-##         ci[,i] <- quantile(res[,i], probs=c(0.025, 0.975))
-##     }
-
-##     res <- ci
-
-##     expect_that(sum(!(res[1,]<res[2,])),equals(0))
-##     expect_that(sum(!(res[1,]<1)),equals(0))
-##     expect_that(sum(!(res[2,]>1)),equals(0))
-
-## })
-
diff --git a/inst/tests/test-gettaupermute.r b/tests/testthat/test-gettaupermute.r
similarity index 52%
rename from inst/tests/test-gettaupermute.r
rename to tests/testthat/test-gettaupermute.r
index 1f1cccc..48adc6f 100644
--- a/inst/tests/test-gettaupermute.r
+++ b/tests/testthat/test-gettaupermute.r
@@ -1,5 +1,3 @@
-context("get.tau.permute")
-
 test_that("get.tau.permute returns appropriate values for test case 1 (equilateral triangle)" ,{
 
     x <- rbind(c(1,0,0), c(1,1,0),c(2,.5,sqrt(.75)))
@@ -97,84 +95,3 @@ test_that ("fails nicely if x and y column names are not provided", {
 
 })
 
-
-##################DEPRECATED TESTS...TAKE TO LONG...NOW USING SMALLER CANONICAL
-##################TESTS THAT HAVE VALUES THAT CAN BE WORKED OUT BY HAND
-## test_that("get.tau.permute cis enclose get.tau when no clustering exists",
-##           {
-##               set.seed(788)
-
-##               x<-cbind(rep(c(1,2),250), x=runif(500,0,100), y=runif(500,0,100))
-
-##               colnames(x) <-c("type","x","y")
-
-##               test <- function(a,b) {
-##                   if (a[1] != 1) return(3)
-##                   if (b[1] == 1) return(1)
-##                   return(2)
-##               }
-
-##               #plot(x[,"x"],x[,"y"], col=x[,"type"])
-
-##               res <- get.tau.permute(x, test, seq(10,100,10), seq(0,90,10), 200)
-##               res2 <- get.tau(x, test, seq(10,100,10), seq(0,90,10))
-
-##               for (i in 1:10) {
-##                   tmp <- quantile(res[,i], probs=c(0.025, .975), na.rm=T)
-##                   #print(res2[i])
-##                   #print(tmp)
-##                   expect_that(res2[i]>=tmp[1], is_true())
-##                   expect_that(res2[i]<=tmp[2], is_true())
-##               }
-##           })
-
-## test_that("get.tau.permute cis do not enclose get.tau at extremes when no clustering exists",
-##           {
-##               set.seed(788)
-
-##               #first generate 200 random uniform points
-##               x<-cbind(1, x=runif(300,0,100), y=runif(300,0,100))
-##               colnames(x) <-c("type","x","y")
-
-##               #add a seed point
-##               x<-rbind(x,c(2,50,50))
-
-##               #generate 200 normally distibuted points around this
-##               x<-rbind(x,cbind(3,rnorm(300,50,20),rnorm(300,50,20)))
-
-##               test <- function(a,b) {
-##                   if (a[1] != 2) return(3)
-##                   if (b[1] == 3) return(1)
-##                   return(2)
-##               }
-
-##               ## res <- get.tau.permute(x,test,seq(10,50,10), seq(0,40,10), 200)
-##               res2 <- get.tau(x,test,seq(10,50,10), seq(0,40,10))
-
-
-##               for (i in c(1,5)) {
-##                   tmp <- quantile(res[,i], probs=c(0.025, .975), na.rm=T)
-##                   #print(tmp)
-##                   #print(res2[i])
-##                   expect_that((res2[i]>=tmp[1]) & (res2[i]<=tmp[2]) ,
-##                               is_false())
-##               }
-
-
-
-##               res <- get.tau.typed.permute(x,2,3,seq(10,50,10),
-##                                            seq(0,40,10), 100)
-
-
-##               for (i in c(1,5)) {
-##                   tmp <- quantile(res[,i], probs=c(0.025, .975), na.rm=T)
-##                   #print(res2[i])
-##                   #print(tmp)
-##                   expect_that((res2[i]>=tmp[1]) & (res2[i]<=tmp[2]) ,
-##                               is_false())
-##               }
-
-
-##           })
-
-
diff --git a/inst/tests/test-gettheta.r b/tests/testthat/test-gettheta.r
similarity index 99%
rename from inst/tests/test-gettheta.r
rename to tests/testthat/test-gettheta.r
index 264ca80..03423df 100644
--- a/inst/tests/test-gettheta.r
+++ b/tests/testthat/test-gettheta.r
@@ -1,4 +1,3 @@
-context("get.theta")
 test_that("get.theta returns Inf when all relations are 1", { #Would throwing an error be better?
     #generate a set of 100 random points even labeled between the two
     x<-cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
diff --git a/inst/tests/test-getthetabootstrap.r b/tests/testthat/test-getthetabootstrap.r
similarity index 83%
rename from inst/tests/test-getthetabootstrap.r
rename to tests/testthat/test-getthetabootstrap.r
index a4372dd..1d35990 100644
--- a/inst/tests/test-getthetabootstrap.r
+++ b/tests/testthat/test-getthetabootstrap.r
@@ -1,6 +1,4 @@
-context("get.theta.bootstrap")
-
-test_that("get.theta.boostrap runs and returns Inf when all relations are 1", {
+test_that("get.theta.bootstrap runs and returns Inf when all relations are 1", {
 
     x<-cbind(rep(c(1,2),50), x=runif(100,0,100), y=runif(100,0,100))
 
@@ -31,22 +29,16 @@ test_that("get.theta.ci returns bootstrap cis when same seed", {
     res <- get.theta.bootstrap(x, test, seq(15,45,15), seq(0,30,15), 20)[,-(1:2)]
 
     set.seed(787)
-    ci1 <- get.theta.ci(x, test, seq(15,45,15), seq(0,30,15), 20)[,-(1:3)]
-
-    expect_that(as.numeric(ci1[1,]),
-                equals(as.numeric(quantile(res[1,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+    ci1 <- get.theta.ci(x, test, seq(15,45,15), seq(0,30,15), 20, ci.level = 0.95)[,-(1:3)]
 
+    expect_that(as.numeric(ci1[1,]), 
+                equals(coxed::bca(as.numeric(res[1,]),conf.level = 0.95)))
+    
     expect_that(as.numeric(ci1[2,]),
-                equals(as.numeric(quantile(res[2,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
-
+                equals(coxed::bca(as.numeric(res[2,]),conf.level = 0.95)))
+    
     expect_that(as.numeric(ci1[3,]),
-                equals(as.numeric(quantile(res[3,],
-                                           probs=c(.025,.975),
-                                           na.rm=T))))
+                equals(coxed::bca(as.numeric(res[3,]),conf.level = 0.95)))
 
 })
 
@@ -65,7 +57,8 @@ test_that("performs correctly for test case 1 (equilateral triangle)", {
     res2 <- get.theta.typed.bootstrap(x, 1,2, 1.5, 0.1, 500)[,-(1:3)]
 
 
-    #should have 95% CI of 0,1 and mean/median of 0.5
+    #should have 95% CI of 0,1 and mean/median of 0.5. quantile() method used as coxed::bca() breaks
+    # down under Inf conditions
     expect_that(as.numeric(quantile(res[1,], probs=c(.025,.975), na.rm=T)),
                 equals(c(0,Inf)))
     expect_that(as.numeric(quantile(res2[1,], probs=c(.025,.975), na.rm=T)),
@@ -88,7 +81,8 @@ test_that("performs correctly for test case 2 (points on a line)", {
         return(2)
     }
 
-    #the medians for the null distribution should be 1,0.5,0
+    #the medians for the null distribution should be 1,0.5,0. quantile() method used as 
+    #coxed::bca() breaks down under Inf conditions
     res <- get.theta.bootstrap(x, test, c(1.5,2.5,3.5), c(0,1.5,2.5), 500)[,-(1:3)]
     res2 <- get.theta.typed.bootstrap(x, 1, 2, c(1.5,2.5,3.5), c(0,1.5,2.5), 500)[,-(1:3)]
 
diff --git a/inst/tests/test-getthetapermute.r b/tests/testthat/test-getthetapermute.r
similarity index 85%
rename from inst/tests/test-getthetapermute.r
rename to tests/testthat/test-getthetapermute.r
index 7444711..e5bf9a8 100644
--- a/inst/tests/test-getthetapermute.r
+++ b/tests/testthat/test-getthetapermute.r
@@ -1,5 +1,3 @@
-context("get.theta.permute")
-
 test_that("get.theta.permute returns appropriate values for test case 1 (equilateral triangle)" ,{
 
     x <- rbind(c(1,0,0), c(1,1,0),c(2,.5,sqrt(.75)))
@@ -34,9 +32,9 @@ test_that("get.theta.permute returns appropriate values for test case 2 (points
         return(2)
     }
 
-    #the median of the null distribution should be 1 (includes infs so
-    #  mean does not work)
-    #the 95% CI equals 0,Inf with windows
+    #the median of the null distribution should be 1 (includes infs so mean does not work)
+    #the 95% CI equals 0,Inf with windows. As it contains Infs we use quantile() rather than bca()
+    #to compute the confidence interval
     res <- get.theta.permute(x, test, c(1.5,2.5,3.5), c(0,1.5,2.5), 500)[,-(1:2)]
     res2 <- get.theta.typed.permute(x, 1, 2, c(1.5,2.5,3.5), c(0,1.5,2.5), 500)[,-(1:2)]
 
@@ -50,9 +48,12 @@ test_that("get.theta.permute returns appropriate values for test case 2 (points
                     equals(c(0,Inf)))
     }
 
-    #without windows the 95% CI should be 1/3 and 3
+    #without windows the 95% CI should be 1/3 and 3. As it contains Infs we use quantile() rather 
+    # than coxed::bca() to compute the confidence interval
     res <- get.theta.permute(x, test, 4,0, 500)[,-(1:2)]
     res2 <- get.theta.typed.permute(x, 1, 2, 4,0, 500)[,-(1:2)]
+    
+    
     expect_that(as.numeric(quantile(res[1,], probs=c(.025,.975))),
                 equals(c(1/3,3)))
     expect_that(as.numeric(quantile(res2[1,], probs=c(.025,.975))),
diff --git a/inst/tests/test-simulateepidemic.r b/tests/testthat/test-simulateepidemic.r
similarity index 99%
rename from inst/tests/test-simulateepidemic.r
rename to tests/testthat/test-simulateepidemic.r
index a2232a1..ae1fdcd 100644
--- a/inst/tests/test-simulateepidemic.r
+++ b/tests/testthat/test-simulateepidemic.r
@@ -1,5 +1,3 @@
-context("sim.epidemic")
-
 test_that("Plausible parameter values produce simulations: R", {
      
      for(i in seq(1, 2, 0.25)) {
diff --git a/inst/tests/test-thetaweights.r b/tests/testthat/test-thetaweights.r
similarity index 99%
rename from inst/tests/test-thetaweights.r
rename to tests/testthat/test-thetaweights.r
index 66c20c7..9e2b6c0 100644
--- a/inst/tests/test-thetaweights.r
+++ b/tests/testthat/test-thetaweights.r
@@ -1,5 +1,3 @@
-context("theta weights")
-
 test_that("Correct array for theta values is returned", {
      
      case.times <- c(1,2,2,3,3)
diff --git a/inst/tests/test-transdist.r b/tests/testthat/test-transdist.r
similarity index 98%
rename from inst/tests/test-transdist.r
rename to tests/testthat/test-transdist.r
index 1803190..95beb34 100644
--- a/inst/tests/test-transdist.r
+++ b/tests/testthat/test-transdist.r
@@ -1,5 +1,3 @@
-context("estimate transdist")
-
 test_that("Data checks performed", {
   
   set.seed(1)
diff --git a/inst/tests/test-transdistbootstrapci.r b/tests/testthat/test-transdistbootstrapci.r
similarity index 98%
rename from inst/tests/test-transdistbootstrapci.r
rename to tests/testthat/test-transdistbootstrapci.r
index 0efde22..cf8b59a 100644
--- a/inst/tests/test-transdistbootstrapci.r
+++ b/tests/testthat/test-transdistbootstrapci.r
@@ -1,5 +1,3 @@
-context("estimate transdist bootstraps")
-
 test_that("Data checks performed", {
   
   set.seed(1)
diff --git a/inst/tests/test-transdisttemporal.r b/tests/testthat/test-transdisttemporal.r
similarity index 99%
rename from inst/tests/test-transdisttemporal.r
rename to tests/testthat/test-transdisttemporal.r
index 5175029..e6988b1 100644
--- a/inst/tests/test-transdisttemporal.r
+++ b/tests/testthat/test-transdisttemporal.r
@@ -1,5 +1,3 @@
-context("estimate transdist temporal")
-
 test_that("Data checks performed", {
   
   set.seed(1)
diff --git a/inst/tests/test-wallingateunis.r b/tests/testthat/test-wallingateunis.r
similarity index 99%
rename from inst/tests/test-wallingateunis.r
rename to tests/testthat/test-wallingateunis.r
index fdedf32..f82b64e 100644
--- a/inst/tests/test-wallingateunis.r
+++ b/tests/testthat/test-wallingateunis.r
@@ -1,5 +1,3 @@
-context("Wallinga-Teunis")
-
 test_that("Basic case time WT weights are calculated correctly", {
      
      # Toy example from paper
diff --git a/inst/tests/test-wrapperfuncs.r b/tests/testthat/test-wrapperfuncs.r
similarity index 98%
rename from inst/tests/test-wrapperfuncs.r
rename to tests/testthat/test-wrapperfuncs.r
index c09ea57..5d270d9 100644
--- a/inst/tests/test-wrapperfuncs.r
+++ b/tests/testthat/test-wrapperfuncs.r
@@ -1,5 +1,3 @@
-context("crossK and crossPCF wrapper functions")
-
 test_that("Data checks performed", {
      
      msg <- "epi.data must be a numeric matrix"