Title:	Estimation of Natural Regeneration Potential
Version:	1.0.1
Description:	Functions for estimating the potential dispersal of tree species using regeneration densities and dispersal distances to nearest seed trees. A quantile regression is implemented to determine the dispersal potential. Spatial prediction can be used to identify natural regeneration potential for forest restoration as described in Axer et al (2021) <doi:10.1016/j.foreco.2020.118802>.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
URL:	https://github.com/MaximilianAxer/quaxnat
BugReports:	https://github.com/MaximilianAxer/quaxnat/issues
Encoding:	UTF-8
Depends:	R (≥ 3.5.0)
Imports:	stats, terra
RoxygenNote:	7.3.2
NeedsCompilation:	no
Packaged:	2024-10-28 15:52:57 UTC; maxer
Author:	Maximilian Axer [aut, cre], Robert Schlicht [aut], Robert Nuske [ctb], Nordwestdeutsche Forstliche Versuchsanstalt (NW-FVA) [fnd], Staatsbetrieb Sachsenforst [fnd], Technische Universität Dresden [fnd, cph]
Maintainer:	Maximilian Axer <maximilian.axer@nw-fva.de>
Repository:	CRAN
Date/Publication:	2024-10-28 16:50:07 UTC

Extracting distances to nearest seed source for point data

Description

Extracts a distance for the inventory plots. The distance to the nearest seed source is used for the analysis of the regeneration potential.

Usage

extract_dist(raster, geom, species)

Arguments

Value

Numeric vector with distances of every inventory plot to the nearest seed source of a specific tree species.

Examples

## Create raster data set
set.seed(2023)
rr <- terra::rast(
matrix(sample(0:10, 20 * 20, replace = TRUE),
        nrow = 20, ncol = 20))

## Create vector data set
vec <- terra::vect(rbind(c(5,10), c(5,15)))

## Extract distance for the inventory plot
extract_dist(raster=rr, geom=vec, species="10")

Dispersal kernels from exponential power family

Description

k_exponential_power computes the value, multiplied by N, of a dispersal kernel from the exponential power family that includes, as special cases, Gaussian kernels and kernels that follow an exponential function of the distance.

Usage

k_exponential_power(x, par, N = 1, d = NCOL(x))

Arguments

Details

The dispersal kernel, i.e. spatial probability density of the location of a seed relative to its source, is here given by

k(x)={b\Gamma (d/2) \over 2\pi ^{d/2}a^{d}\Gamma (d/b)} e^{-(\left\|{x}\right\|/a)^{b}},

which corresponds to a probability density of the distance given by

p(r)={b \over a^{d}\Gamma (d/b)}r^{d-1}e^{-(r/a)^{b}},

where d is the spatial dimension, \left\|{\,}\right\| denotes the Euclidean norm and the normalizing constants involve the gamma function; see Bateman (1947), Clark et al. (1998), Austerlitz et al. (2004), Nathan et al. (2012) for the planar case. This means the bth power of the distance has a gamma distribution with shape parameter d/b and scale parameter a^{b}.

The kernel has its maximum at zero and represents a rather flexible family that includes, for b=2 the classical Gaussian kernels and for b=1, kernels decreasing exponentially with the distance. For b<1 the distance distribution is fat-tailed in the sense of Kot et al. (1996). Such kernels have consequently been applied in a number of theoretical studies that address dispersal (Ribbens et al. 1994, Bullock et al. 2017).

Value

Numeric vector of function values k(x) multiplied by N.

References

Bateman, A. (1947). Contamination in seed crops: III. relation with isolation distance. Heredity 1, 303–336. doi:10.1038/hdy.1947.20

Kot, M., Lewis, M.A., van den Driessche, P. (1996). Dispersal Data and the Spread of Invading Organisms. Ecology 77(7), 2027–2042. doi:10.2307/2265698

Ribbens, E., Silander Jr, J.A., Pacala, S.W. (1994). Seedling recruitment in forests: calibrating models to predict patterns of tree seedling dispersion. Ecology 75, 1794–1806. doi:10.2307/1939638

Clark, J.S., Macklin, E., Wood, L. (1998). Stages and spatial scales of recruitment limitation in southern Appalachian forests. Ecological Monographs 68(2), 213–235. doi:10.2307/2657201

Clark, J.S. (1998). Why trees migrate so fast: confronting theory with dispersal biology and the paleorecord. The American Naturalist 152(2), 204–224. doi:10.1086/286162

Austerlitz, F., Dick, C.W., Dutech, C., Klein, E.K., Oddou-Muratorio, S., Smouse, P.E., Sork, V.L. (2004). Using genetic markers to estimate the pollen dispersal curve. Molecular Ecology 13, 937–954. doi:10.1111/j.1365-294X.2004.02100.x

Bullock, J. M., Mallada González, L., Tamme, R., Götzenberger, L., White, S.M., Pärtel, M., Hooftman, D.A. (2017). A synthesis of empirical plant dispersal kernels. Journal of Ecology 105, 6–19. doi:10.1111/1365-2745.12666

Nathan, R., Klein, E., Robledo‐Arnuncio, J.J., Revilla, E. (2012). Dispersal kernels: review, in Clobert, J., Baguette, M., Benton, T.G., Bullock, J.M. (eds.), Dispersal ecology and evolution, 186–210. doi:10.1093/acprof:oso/9780199608898.003.0015

Examples

k_exponential_power(2:5, par=c(0,0), d=2)

Dispersal kernels for log-normal distance distributions

Description

k_lognormal computes the value, multiplied by N, of a dispersal kernel based on seeds having a distance with a log-normal distribution from the their source.

Usage

k_lognormal(x, par, N = 1, d = NCOL(x))

Arguments

Details

The dispersal kernel, i.e. spatial probability density of the location of a seed relative to its source, is here given by

k(x)={\Gamma (d/2) \over 2\pi ^{d/2}\left\|{x}\right\|^{d}\sqrt{2\pi b}} e^{-{1 \over 2b}(\log (\left\|{x}\right\|/a))^{2}} ={\Gamma (d/2)e^{d^{2}b/2} \over 2\pi ^{d/2}a^{d}\sqrt{2\pi b}} e^{-{1 \over 2b}(\log {\left\|{x}\right\| \over a}+db)^{2}},

which corresponds to a probability density of the distance given by

p(r)={1 \over r\sqrt{2\pi b}}e^{-{1 \over 2b}(\log (r/a))^{2}} ={e^{b/2} \over a\sqrt{2\pi b}} e^{-{1 \over 2b}(\log {r \over a}+b)^{2}},

where d is the spatial dimension, \left\|{\,}\right\| denotes the Euclidean norm and the normalizing constant of the kernel involves the gamma function; see Greene and Johnson (1989), Stoyan and Wagner (2001) for the planar case. Thus, the distance is assumed to have the log-normal distribution such that the log-distance has a normal distribution with mean \log a and variance b. Here \log k(x) is a quadratic function of \log \left\|{x}\right\| with a maximum at \log a-db, while \log p(r) is a quadratic function of \log r with a maximum at \log a-b.

This kernel is particularly suitable if the maximum regeneration density is not directly at the seed source (e.g. Janzen–Connell effect), cf. Nathan et al. (2012).

Value

Numeric vector of function values k(x) multiplied by N.

References

Greene, D.F., Johnson, E.A. (1989). A model of wind dispersal of winged or plumed seeds. Ecology 70(2), 339–347. doi:10.2307/1937538

Stoyan, D., Wagner, S. (2001). Estimating the fruit dispersion of anemochorous forest trees. Ecol. Modell. 145, 35–47. doi:10.1016/S0304-3800(01)00385-4

Nathan, R., Klein, E., Robledo‐Arnuncio, J.J., Revilla, E. (2012). Dispersal kernels: review, in Clobert, J., Baguette, M., Benton, T.G., Bullock, J.M. (eds.), Dispersal ecology and evolution, 186–210. doi:10.1093/acprof:oso/9780199608898.003.0015

Examples

k_lognormal(2:5, par=c(0,0), d=2)

Power-law dispersal kernels

Description

k_power computes the value, multiplied by N, of a dispersal kernel that follows a power law of a constant a plus the distance.

Usage

k_power(x, par, N = 1, d = NCOL(x))

Arguments

Details

The dispersal kernel, i.e. spatial probability density of the location of a seed relative to its source, is here given by

k(x)={\Gamma (d/2) \over 2\pi ^{d/2}a^{d}B(d,b)} (1+{\left\|{x}\right\| \over a})^{-(b+d)},

which corresponds to a probability density of the distance given by

p(r)={1 \over a^{d}B(d,b)}r^{d-1}(1+{r \over a})^{-(b+d)},

where d is the spatial dimension, \left\|{\,}\right\| denotes the Euclidean norm and the normalizing constants involve the beta and gamma functions; see Nathan et al. (2012) for the planar case (with b replaced by b-d). This means the distance is da \over b times a random variable having an F distribution with 2d and 2b degrees of freedom. This is a fat-tailed distribution for all choices of the parameter b.

Value

Numeric vector of function values k(x) multiplied by N.

References

Nathan, R., Klein, E., Robledo‐Arnuncio, J.J., Revilla, E. (2012). Dispersal kernels: review, in Clobert, J., Baguette, M., Benton, T.G., Bullock, J.M. (eds.), Dispersal ecology and evolution, 186–210. doi:10.1093/acprof:oso/9780199608898.003.0015

Austerlitz, F., Dick, C.W., Dutech, C., Klein, E.K., Oddou-Muratorio, S., Smouse, P.E., Sork, V.L. (2004). Using genetic markers to estimate the pollen dispersal curve. Molecular Ecology 13, 937–954. doi:10.1111/j.1365-294X.2004.02100.x

Examples

k_power(2:5, par=c(0,0), d=2)

Dispersal kernels from spatial t distribution

Description

k_t computes the value, multiplied by N, of the dispersal kernel from Clark et al. (1999) that represents a multivariate t distribution.

Usage

k_t(x, par, N = 1, d = NCOL(x))

Arguments

Details

The dispersal kernel, i.e. spatial probability density of the location of a seed relative to its source, is here given by

k(x)={\Gamma ((b+d)/2) \over \pi ^{d/2}a^{d}\Gamma (b/2)} (1+{\left\|{x}\right\|^{2} \over a^{2}})^{-(b+d)/2},

which corresponds to a probability density of the distance given by

p(r)={2 \over a^{d}B(d/2,b/2)}r^{d-1} (1+{r^{2} \over a^{2}})^{-(b+d)/2},

where d is the spatial dimension, \left\|{\,}\right\| denotes the Euclidean norm and the normalizing constants involve the beta and gamma functions; see Clark et al. (1999) and Austerlitz et al. (2004) for the planar case (with a,b replaced by \sqrt{u},2p and a,2b-d, respectively). This means the position is a \over \sqrt{b} times a random vector having a standard d-variate t distribution with b degrees of freedom (a standard Gaussian vector divided by \sqrt{z/b}, where z is independent and chi-squared distributed with b degrees of freedom), and the squared distance is da^{2} \over b times a random variable having an F distribution with d and b degrees of freedom.

This results from the kernel being defined as a mixture of Gaussian kernels with an inverse variance having a gamma distribution with shape parameter b\over 2 and inverse scale parameter a^{2}\over 2, which for a=1 is a chi-squared distribution with b degrees of freedom.

The dispersal kernel always has its maximum at zero, and the distance has a fat-tailed distribution for all choices of b.

Value

Numeric vector of function values k(x) multiplied by N.

References

Clark, J.S., Silman, M., Kern, R., Macklin, E., HilleRisLambers, J. (1999). Seed dispersal near and far: patterns across temperate and tropical forests. Ecology 80, 1475–1494. doi:10.1890/0012-9658(1999)080[1475:SDNAFP]2.0.CO;2

Austerlitz, F., Dick, C.W., Dutech, C., Klein, E.K., Oddou-Muratorio, S., Smouse, P.E., Sork, V.L. (2004). Using genetic markers to estimate the pollen dispersal curve. Molecular Ecology 13, 937–954. doi:10.1111/j.1365-294X.2004.02100.x

Examples

k_t(2:5, par=c(0,0), d=2)

Dispersal kernels for Weibull distance distributions

Description

k_weibull computes the value, multiplied by N, of the dispersal kernel from Tufto et al. (1997) based on seeds having a distance with a Weibull distribution from their source.

Usage

k_weibull(x, par, N = 1, d = NCOL(x))

Arguments

Details

The dispersal kernel, i.e. spatial probability density of the location of a seed relative to its source, is here given by

k(x)={b\Gamma (d/2) \over 2\pi ^{d/2}a^{b}}\left\|{x}\right\|^{b-d} e^{-(\left\|{x}\right\|/a)^{b}},

which corresponds to a probability density of the distance given by

p(r)={b \over a^{b}}r^{b-1}e^{-(r/a)^{b}},

where d is the spatial dimension, \left\|{\,}\right\| denotes the Euclidean norm and the normalizing constants involve the gamma function; see Tufto et al. (1997) for the planar case. Thus, the distance is assumed to have the Weibull distribution with scale parameter a and shape parameter b. Equivalently, the bth power of the distance has an exponential distribution with scale parameter a^{b}.

Consequently, if and only if b<1, the distance distribution has a heavier tail than an exponential distribution, although with tail probabilities still decreasing faster than any power law; it is a fat-tailed distribution in the sense of Kot et al. (1996). The kernel coincides with a Gaussian kernel in the special case b=d=2.

Value

Numeric vector of function values k(x) multiplied by N.

References

Tufto, J., Engen, S., Hindar, K. (1997). Stochastic dispersal processes in plant populations, Theoretical Population Biology 52(1), 16–26. doi:10.1006/tpbi.1997.1306

Austerlitz, F., Dick, C.W., Dutech, C., Klein, E.K., Oddou-Muratorio, S., Smouse, P.E., Sork, V.L. (2004). Using genetic markers to estimate the pollen dispersal curve. Molecular Ecology 13, 937–954. doi:10.1111/j.1365-294X.2004.02100.x

Kot, M., Lewis, M.A., van den Driessche, P. (1996). Dispersal Data and the Spread of Invading Organisms. Ecology 77(7), 2027–2042. doi:10.2307/2265698

Nathan, R., Klein, E., Robledo‐Arnuncio, J.J., Revilla, E. (2012). Dispersal kernels: review, in Clobert, J., Baguette, M., Benton, T.G., Bullock, J.M. (eds.), Dispersal ecology and evolution, 186–210. doi:10.1093/acprof:oso/9780199608898.003.0015

Examples

k_weibull(2:5, par=c(0,0), d=2)

Prediction of potential regeneration densities

Description

Prediction of the potential regeneration density as a function of the distance to the nearest seed tree.

Usage

predict_quax(distmap, quax)

Arguments

Details

, defined by the study area. The potential regeneration density is calculated and given for each raster cell.

Value

A SpatRaster with the same resolution as the input raster containing the regeneration density on the same scale (e.g. numbers per hectare) as in the input data.

Examples

## Prepare artificial data:
set.seed(0)
r <- rgamma(200, shape=2, scale=150)
simulated.data <- data.frame(distance = r, density = rpois(length(r),
k_lognormal(r, par=c(6,0), N=1000000, d=2)))

## Run quax function:
f1 <- quax(x = simulated.data$distance, y = simulated.data$density,
          tau = 0.9, fun = k_lognormal)

## Create raster data set
 rr <- terra::rast(
 matrix(sample(0:10, 20 * 20, replace = TRUE),
        nrow = 20, ncol = 20))

## Compute distance for prediction area
distance <- seed_tree_distmap(raster = rr, species = "10")

## Prediction
p <- predict_quax(distmap = distance, quax = f1)
terra::plot(p)

Estimating potential regeneration densities by quantile regression

Description

quax estimates parameters of a spatial dispersal kernel that describes the regeneration potential as the \tau th quantile of the regeneration density. Here \tau is between 0 and 1, with typical values close to 1 representing the situation that the full regeneration potential is realized only at a small fraction of all sites.

Usage

quax(...)

## Default S3 method:
quax(
  ...,
  y,
  tau,
  fun = k_lognormal,
  weights = 1,
  dim = 2,
  par = c(log.a = 8, log.b = 1)
)

## S3 method for class 'formula'
quax(
  formula,
  data,
  tau,
  fun = k_lognormal,
  subset,
  weights,
  na.action,
  offset,
  ...
)

Arguments

Details

The function estimates the parameters N and \theta of the regeneration potential Nk_{\theta } by minimizing

\displaystyle \sum _{i=1}^{n}w_{i}\rho _{\tau }(y_{i}-Nk_{\theta } (x_{i})),

where \rho _{\tau }(u)=\int _{0}^{u}\tau -\mathbf{1}_{s<0}\;ds=\bigl\{ \begin{smallmatrix}u\tau & \text{if }u\geq 0\\u(\tau -1)&\text{if }u<0 \end{smallmatrix} (Koenker and Bassett 1978, Chapter 6.6 in Koenker 2005). The preceding line, after subtracting the same expression for N=0 and substituting s=y_{i}-tk_{\theta }(x_{i}) in the integral, becomes \int _{0}^{N}\sum _{i=1}^{n}w_{i}k_{\theta }(x_{i})( \mathbf{1}_{y_{i}<tk_{\theta }(x_{i})}-\tau )\;dt, and any N such that the last integrand is \leq 0 for t<N and \geq 0 for t>N, which can always be found as the integrand is increasing in t, minimizes this integral. The integrand being the difference of the sum of w_{i}k_{\theta }(x_{i}) over the i with y_{i}<tk_{\theta }(x_{i}) and \tau times the sum over all i, with relevant terms for nonzero k_{\theta }(x_{i}), this means that the estimate of N for a given vector \theta can be computed as a \tau th quantile. This is implemented as an inner, nested minimization, the result of which is minimized in \theta using optim.

This is a rather naive approach to quantile regression that appears to work reasonably well for scaled dispersal kernels Nk_{\theta } as considered here, see Appendix A in Axer et al. (2021). For general quantile regression problems the more sophisticated procedure nlrq in the package quantreg, based on Koenker and Park (1996), is expected to provide better results.

In particular, quax is subject to the usual numerical issues inherent in optimization: It can get stuck in a local minimum or altogether miss a minimum if the initial values (as specified by the argument par) are too far off or if the objective function exhibits bad behavior. Problems can further arise in the dispersal kernels if parameter values passed on a log scale become too large. It is therefore recommended to visually check the results (see Examples). Also, the optim arguments method and control can be added in ... to select and tune the optimization algorithm, but note that the objective function is usually not differentiable.

See Koenker (2005) for a detailed exposition of quantile regression.

Value

An object of class quax containing the estimated function, including an attribute o containing the results of optim. Generic functions with methods defined for quax objects invoke these methods; see summary.quax for an example.

References

Koenker, R., Bassett, G. (1978). Regression quantiles. Econometrica 46(1), 33–50. doi:10.2307/1913643

Axer, M., Schlicht, R., Wagner, S. (2021). Modelling potential density of natural regeneration of European oak species (Quercus robur L., Quercus petraea (Matt.) Liebl.) depending on the distance to the potential seed source: Methodological approach for modelling dispersal from inventory data at forest enterprise level. Forest Ecology and Management 482, 118802. doi:10.1016/j.foreco.2020.118802

Koenker, R., Park, B.J. (1996). An interior point algorithm for nonlinear quantile regression. Journal of Econometrics 71(1–2), 265–283. doi:10.1016/0304-4076(96)84507-6

Koenker, R. (2005). Quantile regression. Cambridge University Press. doi:10.1017/CBO9780511754098

See Also

Function nlrq in the package quantreg.

Examples

## Prepare artificial data:
set.seed(0)
r <- rgamma(200, shape=2, scale=150)
simulated.data <- data.frame(distance = r, density =
  rpois(length(r), k_lognormal(r, par=c(6,0), N=1000000, d=2)))
plot(density ~ distance, simulated.data)

## Run quax function:
f1 <- quax(x = simulated.data$distance, y = simulated.data$density,
  tau = 0.9, fun = k_lognormal)
summary(f1)
curve(f1(x), add=TRUE)

## Do the same using formula interface:
f1 <- quax(density ~ distance, simulated.data,
  tau = 0.9, fun = k_lognormal)
summary(f1)
#quantreg::nlrq(density ~ k_lognormal(distance,c(log.a,log.b),N=N,d=2),
#  simulated.data, start = c(log.a=6,log.b=0,N=1e6), tau = 0.9)  # similar

## Use another quantile:
f2 <- quax(density ~ distance, simulated.data,
  tau = 0.99, fun = k_lognormal)
summary(f2)
curve(f2(x), add=TRUE, lwd=0)

## Show effect of weights:
f3 <- quax(density ~ distance, simulated.data,
  tau = 0.9, fun = k_lognormal, weights = distance)
summary(f3)
curve(f3(x), add=TRUE, lty=3)

## Compare various dispersal models:
fun <- c("k_lognormal","k_t","k_weibull","k_power","k_exponential_power")
for (i in seq_along(fun))
  curve(quax(density ~ distance, simulated.data,
    tau = 0.9, fun = get(fun[i]), weights = distance)(x),
    add=TRUE, col=i, lty=3)
legend("topright", fun, col=seq_along(fun), lty=3)

## Use positions in computation:
simulated.data$position <- r *
  (\(a) cbind(cos(a),sin(a))) (rnorm(length(r)))
f3 <- quax(density ~ position, simulated.data,
  tau = 0.9, fun = k_lognormal, weights = distance)
summary(f3)

## Show problems with bad initial values and try another parameterization:
curve(quax(density ~ distance, simulated.data, par = c(log.a=0,log.b=0),
  tau = 0.99, fun = k_lognormal)(x), add=TRUE, lty=2)
curve(quax(density ~ distance, simulated.data, par = c(a=1,b=1),
  tau = 0.99, fun = function(x,par,N,d) if (any(par<=0)) rep(NA,NROW(x))
    else k_lognormal(x,log(par),N,d))(x), add=TRUE, lty=2)

## Use custom variant of lognormal model that includes a shift:
plot(simulated.data$position)
f4 <- quax(density ~ position, simulated.data,
  tau = 0.9, par = c(8, 1, 0, 0),
  fun = function(x, par, N, d)
    k_lognormal(x - rep(par[-(1:2)],each=NROW(x)), par[1:2], N, d)
)
summary(f4)

Regeneration densities at inventory plots and potential dispersal distances to nearest seed trees

Description

A dataset containing the regeneration densities of beech, oak and Douglas fir of the inventory plots and the distance to the nearest conspecific nearest seed tree.

Usage

data(regeneration)

Format

A data frame with 484 rows and 7 variables

• id. An identifier for each inventory plot as an integer

• distance_beech. Distance in m from the plot to the nearest beech (0–3206.57)

• distance_oak. Distance in m from the plot to the nearest oak (0–1481.2)

• distance_dgl. Distance in m from the plot to the nearest Douglas fir (0–1807)

• oak_regen. Regeneration density of oak (N/ha) of the plot (0–30)

• beech_regen. Regeneration density of beech (N/ha) of the plot (0–30)

• douglas_regen. Regeneration density of Douglas fir (N/ha) of the plot (0–30)

Distance map to nearest seed source

Description

Creation of a distance map for the study area. The distance to the nearest seed source is calculated for every raster cell.

Usage

seed_tree_distmap(raster, species)

Arguments

Value

A SpatRaster object containing the distances to seed source. The object has the same resolution and extent as the input raster.

Examples

## Create raster data set
 rr <- terra::rast(
 matrix(sample(0:10, 20 * 20, replace = TRUE),
        nrow = 20, ncol = 20))

## Compute distance for study area
distance <- seed_tree_distmap(raster = rr, species = "10")

## Plot the seed_tree_distmap
terra::plot(distance)

Summarizing quantile regression fits of potential regeneration densities

Description

This function is the summary method for class quax objects as returned by quax.

Usage

## S3 method for class 'quax'
summary(object, ...)

Arguments

Details

The value component of the result can be used to compare the quality of the fit of different dispersal kernels for the same quantile to the same data.

Value

A list with the following components:

coefficients

The parameters of the estimated dispersal kernel.

value

The attained value of the objective function that is minimised in the quantile regression.

Examples

## Prepare artificial data:
set.seed(0)
r <- rgamma(200, shape=2, scale=150)
simulated.data <- data.frame(distance = r, density =
  rpois(length(r), k_lognormal(r, par=c(6,0), N=1000000, d=2)))
plot(density ~ distance, simulated.data)

## Fit a log-normal and a power-law dispersal kernel to the data:
f1 <- quax(density ~ distance, simulated.data,
  tau = 0.9, fun = k_lognormal)
f2 <- quax(density ~ distance, simulated.data,
  tau = 0.9, fun = k_power)

## Compare both fits:
summary(f1)
summary(f2)