version 0.1.0

still a lot of work to do, but code is working and complete for now

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: ecotoner
 Title: ecotoner - ecotones / ecological boundaries across environmental gradients
-Version: 0.0.0.9051
+Version: 0.1.0
 Authors@R: person("Daniel", "Schlaepfer", email = "daniel.schlaepfer@unibas.ch", role = c("aut", "cre"))
 Description: ecotoner locates ecotones and extracts geographic and geometric measures of ecotones
 Depends:

R/classes_input.R

@@ -18,7 +18,7 @@ GridInfo <- setClass("GridInfo",
 	prototype = list(crs = sp::CRS(NA_character_), res_m = NA_real_, longlat = NA, origin = c(NA_real_, NA_real_), rotation = NA)
 )
 
-# TODO(drs): validify that units of slot 'res_m' are correct
+# TODO(drs): validate that units of slot 'res_m' are correct
 
 #' An S4-class to represent the information of the input raster grids
 #'

R/locate_flowpath.R

@@ -260,7 +260,6 @@ calc.Identify_GoodvsBadMigration <- function(patches4, tally, paths, paths_succe
 	patchIDs_crossedByPaths <- lapply(paths, FUN = function(p) unique(na.omit(raster::extract(patches4, p)))) #entries in order of downstream patches4 crossed
 	#patches4 that connect directly with the x-border
 	patchIDs_BadMigration <- tally[tally[, "PatchCoversGoodBorder_TF"] == 1 | (tally[, "SuccessfulFlowpaths_N"] > 0 & tally[, "GoodMigration_N"] > 0 & tally[, "BadMigration_N"] == 0), "PatchID"]
-#TODO(drs): shouldn't this be 'patchIDs_GoodMigration'??
 
 	if (end_toLeft) {
 		#Check network of flowpaths (originating from every y-border cell) and decide for each patch whether to be excluded or not

R/measure_eppinga.R

@@ -1,6 +1,6 @@
 #------Eppinga, M.B., Pucko, C.A., Baudena, M., Beckage, B. & Molofsky, J. (2013) A new method to infer vegetation boundary movement from 'snapshot' data. Ecography, 36, 622-635.
 
-version_Eppinga2013Ecography <- function() numeric_version("0.2.3")
+version_Eppinga2013Ecography <- function() numeric_version("0.2.4")
 
 
 #---Eppinga et al. 2013: 'Analytical analysis of vegetation boundary movement' (Fig. 1)
@@ -63,110 +63,108 @@ calc_Eppinga2013_stats <- function(FR_dist_T17_veg1, FR_dist_mean_T17_veg1, FR_d
 	res[["FrontsAdvBeyondOptBoundary"]] <- all(FR_dist_mean_T17_veg1 > 0, FR_dist_mean_T17_veg2 < 0)
 
 
-	if (res[["FrontsAdvBeyondOptBoundary"]]) {
-		if (!is.na(seed)) set.seed(seed)
-
-		#---Test if vegetation types have advanced comparably
-		#	i.e., test if difference between front runner distances (FD) of veg1 vs veg2 is 0: eq. 18
-		data_T17 <- cbind(FR_dist_T17_veg1, -FR_dist_T17_veg2)
-		delta_T17 <- apply(data_T17, 1,  function(x) x[1] - x[2])
-		res[["FD_mean_T17"]] = mean(delta_T17, na.rm = TRUE)
-		res[["FD_sd_T17"]] = sd(delta_T17, na.rm = TRUE)
-		res[["FD_mean_m"]] = backtransformation17(res[["FD_mean_T17"]])
-		res[["FD_boots_R"]] <- 1e5L # Eppinga et al. 2013: R = 1e5
-
-		bmds <- list(iid = list("boot" = NULL, "ci" = NULL),
-					 dep = list("boot" = NULL, "ci" = NULL))
+	if (!is.na(seed)) set.seed(seed)
+
+	#---Test if vegetation types have advanced comparably
+	#	i.e., test if difference between front runner distances (FD) of veg1 vs veg2 is 0: eq. 18
+	data_T17 <- cbind(FR_dist_T17_veg1, -FR_dist_T17_veg2)
+	delta_T17 <- apply(data_T17, 1,  function(x) x[1] - x[2])
+	res[["FD_mean_T17"]] = mean(delta_T17, na.rm = TRUE)
+	res[["FD_sd_T17"]] = sd(delta_T17, na.rm = TRUE)
+	res[["FD_mean_m"]] = backtransformation17(res[["FD_mean_T17"]])
+	res[["FD_boots_R"]] <- 1e5L # Eppinga et al. 2013: R = 1e5
+
+	bmds <- list(iid = list("boot" = NULL, "ci" = NULL),
+				 dep = list("boot" = NULL, "ci" = NULL))
+
+	if (requireNamespace("boot", quietly = TRUE)) {
+		#- Bootstrap approach assuming independent data (as used by Eppinga et al. 2013)
+		bmds[["iid"]][["boot"]] <- boot::boot(data = data_T17,
+									  statistic = indexed_mean_of_diffs,
+									  R = res[["FD_boots_R"]],
+									  sim = "ordinary", stype = "i",
+									  parallel = "no")
 
-		if (requireNamespace("boot", quietly = TRUE)) {
-			#- Bootstrap approach assuming independent data (as used by Eppinga et al. 2013)
-			bmds[["iid"]][["boot"]] <- boot::boot(data = data_T17,
-										  statistic = indexed_mean_of_diffs,
-										  R = res[["FD_boots_R"]],
-										  sim = "ordinary", stype = "i",
-										  parallel = "no")
+		# bias corrected and accelerated bootstrap (BCa) interval
+		bmds[["iid"]][["ci"]] <- boot::boot.ci(bmds[["iid"]][["boot"]],
+									conf = c(0.95, 0.99, 0.999), type = "bca")
+	} else {
+		warning("Package 'boot' not installed: 'calc_Eppinga2013_stats' will not estimate iid bootstrap")
+	}
+
+	if (requireNamespace("boot", quietly = TRUE) && requireNamespace("np", quietly = TRUE)) {
+		# Stationary block bootstrap (Politis, D. N., and J. P. Romano. 1994. The Stationary Bootstrap. Journal of the American Statistical Association 89:1303-1313.)
+		#	with optimal mean block length (Patton, A., D. N. Politis, and H. White. 2009. Correction to "Automatic Block-Length Selection for the Dependent Bootstrap" by D. Politis and H. White (vol 23, pg 53, 2004). Econometric Reviews 28:372-375.)
+
+		if (anyNA(delta_T17)) {
+			# Multiple imputations of time-series data
+			# 'np::b.star' calls acf() and ccf() with default value for 'na.action', i.e., 'na.fail'
+			# 'np::b.star' fails if anyNA(data); data contain NAs if a transect row has no cells of a vegetation type
+			if (requireNamespace("Amelia", quietly = TRUE)) {
+				im_T17 <- cbind(seq_len(nrow(data_T17)), data_T17)
+				am_T17 <- Amelia::amelia(im_T17, m = 10, ts = 1, splinetime = 6,
+										p2s = 0, parallel = "no")
 
-			# bias corrected and accelerated bootstrap (BCa) interval
-			bmds[["iid"]][["ci"]] <- boot::boot.ci(bmds[["iid"]][["boot"]],
-										conf = c(0.95, 0.99, 0.999), type = "bca")
-		} else {
-			warning("Package 'boot' not installed: 'calc_Eppinga2013_stats' will not estimate iid bootstrap")
-		}
-
-		if (requireNamespace("boot", quietly = TRUE) && requireNamespace("np", quietly = TRUE)) {
-			# Stationary block bootstrap (Politis, D. N., and J. P. Romano. 1994. The Stationary Bootstrap. Journal of the American Statistical Association 89:1303-1313.)
-			#	with optimal mean block length (Patton, A., D. N. Politis, and H. White. 2009. Correction to "Automatic Block-Length Selection for the Dependent Bootstrap" by D. Politis and H. White (vol 23, pg 53, 2004). Econometric Reviews 28:372-375.)
-
-			if (anyNA(delta_T17)) {
-				# Multiple imputations of time-series data
-				# 'np::b.star' calls acf() and ccf() with default value for 'na.action', i.e., 'na.fail'
-				# 'np::b.star' fails if anyNA(data); data contain NAs if a transect row has no cells of a vegetation type
-				if (requireNamespace("Amelia", quietly = TRUE)) {
-					im_T17 <- cbind(seq_len(nrow(data_T17)), data_T17)
-					am_T17 <- Amelia::amelia(im_T17, m = 10, ts = 1, splinetime = 6,
-											p2s = 0, parallel = "no")
-
-					est_BstarSB <- mean(sapply(am_T17$imputations, function(x)
-											np::b.star(apply(x, 1,  function(x) x[2] - x[3]), round = TRUE)[, "BstarSB"]))
-				} else {
-					warning("Package 'Amelia' not installed: 'calc_Eppinga2013_stats' cannot estimate optimal block-length for dependent bootstrap with missing data; poor approximation based on complete-cases used instead")
-					est_BstarSB <- np::b.star(delta_T17[complete.cases(delta_T17)], round = TRUE)[, "BstarSB"]
-				}
+				est_BstarSB <- mean(sapply(am_T17$imputations, function(x)
+										np::b.star(apply(x, 1,  function(x) x[2] - x[3]), round = TRUE)[, "BstarSB"]))
 			} else {
-				est_BstarSB <- np::b.star(delta_T17, round = TRUE)[, "BstarSB"]
+				warning("Package 'Amelia' not installed: 'calc_Eppinga2013_stats' cannot estimate optimal block-length for dependent bootstrap with missing data; poor approximation based on complete-cases used instead")
+				est_BstarSB <- np::b.star(delta_T17[complete.cases(delta_T17)], round = TRUE)[, "BstarSB"]
 			}
-
-			res[["FD_depboot_bstar"]] <- min(nrow(data_T17), est_BstarSB)
-			bmds[["dep"]][["boot"]] <- boot::tsboot(tseries = data_T17,
-										 statistic = indexed_mean_of_diffs,
-										 R = res[["FD_boots_R"]],
-										 sim = "geom", l = res[["FD_depboot_bstar"]], endcorr = TRUE, n.sim = nrow(data_T17),
-										 orig.t = TRUE, parallel = "no")
-
-			# BCa and studentized CI don't apply for tsboot objects; use instead percentile method
-			bmds[["dep"]][["ci"]] <- boot::boot.ci(bmds[["dep"]][["boot"]],
-										conf = c(0.95, 0.99, 0.999), type = "perc")
 		} else {
-			warning("Package 'boot' and/or 'np' not installed: 'calc_Eppinga2013_stats' will not estimate dependent bootstrap")
+			est_BstarSB <- np::b.star(delta_T17, round = TRUE)[, "BstarSB"]
 		}
+
+		res[["FD_depboot_bstar"]] <- min(nrow(data_T17), est_BstarSB)
+		bmds[["dep"]][["boot"]] <- boot::tsboot(tseries = data_T17,
+									 statistic = indexed_mean_of_diffs,
+									 R = res[["FD_boots_R"]],
+									 sim = "geom", l = res[["FD_depboot_bstar"]], endcorr = TRUE, n.sim = nrow(data_T17),
+									 orig.t = TRUE, parallel = "no")
+
+		# BCa and studentized CI don't apply for tsboot objects; use instead percentile method
+		bmds[["dep"]][["ci"]] <- boot::boot.ci(bmds[["dep"]][["boot"]],
+									conf = c(0.95, 0.99, 0.999), type = "perc")
+	} else {
+		warning("Package 'boot' and/or 'np' not installed: 'calc_Eppinga2013_stats' will not estimate dependent bootstrap")
+	}
 
-		# Extract bootstrap data
-		ptol <- sqrt(.Machine$double.eps)
-		ntol <- -sqrt(.Machine$double.neg.eps)
-		for (ib in names(bmds)) if (!is.null(bmds[[ib]][["boot"]])) {
-			res[[paste0("FD_", ib, "boot_mean")]] = mean(bmds[[ib]][["boot"]]$t, na.rm = TRUE)
-			res[[paste0("FD_", ib, "boot_bias")]] <- res[[paste0("FD_", ib, "boot_mean")]] - res[["FD_mean_T17"]]
-			res[[paste0("FD_", ib, "boot_se")]] <- as.numeric(sqrt(var(bmds[[ib]][["boot"]]$t, na.rm = TRUE)))
-
-			# Test approach 1a: Is 0 contained in ci?
-			i_item <- 4 # the i-th item in the list that is returned by boot::boot.ci
-			res[[paste0("FD_", ib, "boot_ci_type")]] <- which(names(bmds[[ib]][["ci"]])[i_item] == boot_ci_types)
-			conf <- bmds[[ib]][["ci"]][[i_item]]
-			pid <- !(as.integer(apply(conf[, 4:5], 1, function(x) sum(x > ptol))) == 1)
-			res[[paste0("FD_", ib, "boot_ci0_p")]] <- 1 - if (sum(pid) > 0) max(conf[pid, "conf"]) else 0 # this represents steps of 1, 0.05, 0.01, and 0.001 as upper bound of the p-value
-			
-			# Test approach 1b: Calculate p-value (for H0: diff = 0); eq. 19
-			# Direct interpretation of eq. 19: sum(Heaviside(abs(bmdiid$t) + abs(bmdiid$t0) - abs(bmdiid$t + bmdiid$t0))) / bmdiid$R
-			res[[paste0("FD_", ib, "boot_freq_p")]] <- (if (bmds[[ib]][["boot"]]$t0 > ptol) sum(bmds[[ib]][["boot"]]$t <= ptol) else if (bmds[[ib]][["boot"]]$t0 < ntol) sum(bmds[[ib]][["boot"]]$t >= ntol) else bmds[[ib]][["boot"]]$R) / bmds[[ib]][["boot"]]$R
-		}
+	# Extract bootstrap data
+	ptol <- sqrt(.Machine$double.eps)
+	ntol <- -sqrt(.Machine$double.neg.eps)
+	for (ib in names(bmds)) if (!is.null(bmds[[ib]][["boot"]])) {
+		res[[paste0("FD_", ib, "boot_mean")]] = mean(bmds[[ib]][["boot"]]$t, na.rm = TRUE)
+		res[[paste0("FD_", ib, "boot_bias")]] <- res[[paste0("FD_", ib, "boot_mean")]] - res[["FD_mean_T17"]]
+		res[[paste0("FD_", ib, "boot_se")]] <- as.numeric(sqrt(var(bmds[[ib]][["boot"]]$t, na.rm = TRUE)))
+
+		# Test approach 1a: Is 0 contained in ci?
+		i_item <- 4 # the i-th item in the list that is returned by boot::boot.ci
+		res[[paste0("FD_", ib, "boot_ci_type")]] <- which(names(bmds[[ib]][["ci"]])[i_item] == boot_ci_types)
+		conf <- bmds[[ib]][["ci"]][[i_item]]
+		pid <- !(as.integer(apply(conf[, 4:5], 1, function(x) sum(x > ptol))) == 1)
+		res[[paste0("FD_", ib, "boot_ci0_p")]] <- 1 - if (sum(pid) > 0) max(conf[pid, "conf"]) else 0 # this represents steps of 1, 0.05, 0.01, and 0.001 as upper bound of the p-value
+
+		# Test approach 1b: Calculate p-value (for H0: diff = 0); eq. 19
+		# Direct interpretation of eq. 19: sum(Heaviside(abs(bmdiid$t) + abs(bmdiid$t0) - abs(bmdiid$t + bmdiid$t0))) / bmdiid$R
+		res[[paste0("FD_", ib, "boot_freq_p")]] <- (if (bmds[[ib]][["boot"]]$t0 > ptol) sum(bmds[[ib]][["boot"]]$t <= ptol) else if (bmds[[ib]][["boot"]]$t0 < ntol) sum(bmds[[ib]][["boot"]]$t >= ntol) else bmds[[ib]][["boot"]]$R) / bmds[[ib]][["boot"]]$R
+	}
 
-
-		if (requireNamespace("coin", quietly = TRUE)) {
-			# Test approach 2: exact Wilcoxon signed rank test (with Pratt correction of zeros)
-			wsrt <- coin::wilcoxsign_test(FR_dist_T17_veg1 ~ FR_dist_T17_veg2, distribution = "exact", alternative = "two.sided")
-			res[["FD_WSRT_Z"]] <- as.numeric(coin::statistic(wsrt, type = "test"))
-			res[["FD_WSRT_p"]] <- coin::pvalue(wsrt)
-			res[["FD_WSRT_midp"]] <- coin::midpvalue(wsrt)
-		} else {
-			warning("Package 'coin' not installed: 'calc_Eppinga2013_stats' will not calculate Wilcoxon signed rank test")
-		}
-
-		#---Retrospective power: Eppinga et al. 2013: eq. 20
-		tau <- 0.2 #effect size
-		n <- sum(complete.cases(data_T17))
-		tcrit <- qt(0.95, df = n) #95% confidence
-		res[["FD_retro_power"]] <- 1 - (1/2 * (1 + erf((tcrit - tau * sqrt(n) / res[["FD_sd_T17"]]) / (sqrt(2) * res[["FD_sd_T17"]]))))
+
+	if (requireNamespace("coin", quietly = TRUE)) {
+		# Test approach 2: exact Wilcoxon signed rank test (with Pratt correction of zeros)
+		wsrt <- coin::wilcoxsign_test(FR_dist_T17_veg1 ~ FR_dist_T17_veg2, distribution = "exact", alternative = "two.sided")
+		res[["FD_WSRT_Z"]] <- as.numeric(coin::statistic(wsrt, type = "test"))
+		res[["FD_WSRT_p"]] <- coin::pvalue(wsrt)
+		res[["FD_WSRT_midp"]] <- coin::midpvalue(wsrt)
+	} else {
+		warning("Package 'coin' not installed: 'calc_Eppinga2013_stats' will not calculate Wilcoxon signed rank test")
 	}
+
+	#---Retrospective power: Eppinga et al. 2013: eq. 20
+	tau <- 0.2 #effect size
+	n <- sum(complete.cases(data_T17))
+	tcrit <- qt(0.95, df = n) #95% confidence
+	res[["FD_retro_power"]] <- 1 - (1/2 * (1 + erf((tcrit - tau * sqrt(n) / res[["FD_sd_T17"]]) / (sqrt(2) * res[["FD_sd_T17"]]))))
 
 	res
 }

README.md

@@ -1,4 +1,4 @@
-# Ecotoner
+# Ecotoner -- locates ecotones and extracts geographic and geometric measures of ecotones
 
 We haven’t really published the code yet nor prepared it for sharing (though through our use of github made it openly accessible), it is actively and gradually being developed by the Schlaepfer lab, and there is no manual - we cannot give you individual support in setting up and running the code except if we agreed on a collaboration or similar agreement.
 
@@ -9,8 +9,11 @@ There is no graphical user interface, help pages and available documentation may
 If you make use of this model, please cite appropriate references, and we would like to hear about your particular study (especially a copy of any published paper).
 
 
-Some recent references
-* 
+Some references of the implemented methods
+* Danz, N.P., Frelich, L.E., Reich, P.B. & Niemi, G.J. (2012) Do vegetation boundaries display smooth or abrupt spatial transitions along environmental gradients? Evidence from the prairie–forest biome boundary of historic Minnesota, USA. Journal of Vegetation Science, 24, 1129-1140.
+* Eppinga, M.B., Pucko, C.A., Baudena, M., Beckage, B. & Molofsky, J. (2013) A new method to infer vegetation boundary movement from ‘snapshot’ data. Ecography, 36, 622-635.
+* Gastner, M., Oborny, B., Zimmermann, D.K. & Pruessner, G. (2009) Transition from connected to fragmented vegetation across an environmental gradient: scaling laws in ecotone geometry. The American Naturalist, 174, E23-E39.
+
 
 
 
@@ -28,12 +31,12 @@ help(package = "ecotoner")	# The index page with some of the functions (which ar
 ### Using the package
 I added the main code which I use to locate and measure my transects (producing the data for later analysis), as demo to the package. I also added data so that it can run  as a small contained demo ‘example’ (if the flag do.demo is set to TRUE). You could run this code directly with the demo() function, but this is probably not convenient:
 ```
-demo("BSE-TF_EcotonesAtSlope_v10", package = "ecotoner")
+demo("BSE-TF_ecotoner_LocateAndMeasure", package = "ecotoner")
 ```
 
 Instead, the following command should open my main code in your text editor (at least on unix-alike systems) for easier inspection:
 ```
-system2("open", file.path(system.file("demo", package = "ecotoner"), "BSE-TF_EcotonesAtSlope_v10.R"))
+system2("open", file.path(system.file("demo", package = "ecotoner"), "BSE-TF_ecotoner_LocateAndMeasure.R"))
 ```
 
 

demo/00Index

@@ -1 +1 @@
-BSE-TF_EcotonesAtSlope_v10	Code to locate, measure, and analyse big sagebrush -- temperate forest ecotones in the western USA.
+BSE-TF_ecotoner_LocateAndMeasure	Code to locate and measure big sagebrush--temperate forest ecotones in the western USA.