#' Compute Range Shift Index
#'
#'
#' The range shift index is a dimensionless measure of the distance of the centroids of two ranges divided by the diameter of the 95\% area. This function uses the 95\% confidence intervals from a range shift fit to calculate a point estimate and 95\% confidence intervals of the RSI.
#'
#' @param FIT a rnage shift object, outputted by \code{\link{estimate_shift}}
#' @param n1 the indices of the ranges to estimate from and to, i.e., for single shift, 1 and 2. For three ranges (two shifts) it can be 1 and 2, 2 and 3, or 1 and 3 - if the ultimate shift is the one of interest.
#' @param nboot number of bootstrap simulation
#' @param n2 see n1
#'
#' @return returns a data frame reporting the distance traveled, the RSI and respective bootstrapped confidence intervals.
#' @export
getRSI <- function(FIT, n1=1, n2=2, nboot = 1e3){
p.hat <- FIT$p.hat
par.CI <- FIT$p.CI
getMCsample <- function(par, n=1e3)
{
p.CI <- (par.CI %>% t %>% as.data.frame)[par][,1]
p.sd <- mean(p.CI[1]-p.CI[2], p.CI[3]-p.CI[1])/2
list(hat = p.CI[1], MC = rnorm(n, p.CI[1], p.sd))
}
x1 <- paste0("x",n1)
x2 <- paste0("x",n2)
y1 <- paste0("y",n1)
y2 <- paste0("y",n2)
if(is.null(par.CI)){
z1 <- p.hat[x1] + 1i*p.hat[y1]
z2 <- p.hat[x2] + 1i*p.hat[y2]
D.hat <- Mod(z2 - z1)
RSI.hat <- D.hat / sqrt(4*p.hat["A"]/pi)
return(c(D = D.hat, RSI = RSI.hat))
} else {
A.mc <- getMCsample("A", nboot)$MC
Z1.mc <- getMCsample(x1, nboot)$MC + 1i * getMCsample(y1, nboot)$MC
Z2.mc <- getMCsample(x2, nboot)$MC + 1i * getMCsample(y2, nboot)$MC
A.hat <- getMCsample("A", nboot)$hat
Z1.hat <- getMCsample(x1, nboot)$hat + 1i * getMCsample(y1, nboot)$hat
Z2.hat <- getMCsample(x2, nboot)$hat + 1i * getMCsample(y2, nboot)$hat
D.mc <- Mod(Z2.mc - Z1.mc)
RSI.mc <- D.mc / sqrt(4*A.mc/pi)
RSI.hat <- Mod(Z2.hat - Z1.hat) / sqrt(4*A.hat/pi)
return(data.frame(
D = c(p.hat = Mod(Z2.hat - Z1.hat), CI.low = as.numeric(quantile(D.mc, 0.025)), CI.high = as.numeric(quantile(D.mc, 0.975))),
RSI = c(p.hat = RSI.hat, CI.low = as.numeric(quantile(RSI.mc, 0.025)), CI.high = as.numeric(quantile(RSI.mc, 0.975)))) %>% t
)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.