update vignette sources

vignettes/main1_pkgintro.Rmd

@@ -2,8 +2,8 @@
 title: "Introduction to the spEDM package"
 author: "Wenbo Lv"
 date: |
-  | Last update: 2025-12-01
-  | Last run: 2025-12-20
+  | Last update: 2026-01-01
+  | Last run: 2026-01-01
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteIndexEntry{1. Introduction to the spEDM package}
@@ -46,7 +46,7 @@ The *spEDM* package includes three illustrative datasets with well-known causati
 2. **Population density and its potential drivers in mainland China** — A county-level dataset capturing population density alongside relevant socio-environmental drivers.
    *Data files*: `popd.csv`, `popd_nb.gal`
 
-3. **Farmland NPP and related variables** — A raster dataset capturing net primary productivity (NPP) of farmland, key climatic variables, elevation, and human activity footprints across mainland China, suitable for analyzing the interactions between agricultural productivity, environmental conditions, and human activities.
+3. **Farmland NPP and related variables** — A raster dataset capturing farmland net primary productivity (NPP), key climatic variables, elevation, and human activity footprints across mainland China, suitable for analyzing the interactions between agricultural productivity, environmental conditions, and human activities.
    *Data files*: `npp.tif`
 
 These datasets can be loaded as shown below:

vignettes/main1_pkgintro.Rmd.orig

@@ -2,7 +2,7 @@
 title: "Introduction to the spEDM package"
 author: "Wenbo Lv"
 date: |
-  | Last update: 2025-12-01
+  | Last update: 2026-01-01
   | Last run: `r Sys.Date()`
 output: rmarkdown::html_vignette
 vignette: >
@@ -52,7 +52,7 @@ The *spEDM* package includes three illustrative datasets with well-known causati
 2. **Population density and its potential drivers in mainland China** — A county-level dataset capturing population density alongside relevant socio-environmental drivers.
    *Data files*: `popd.csv`, `popd_nb.gal`
 
-3. **Farmland NPP and related variables** — A raster dataset capturing net primary productivity (NPP) of farmland, key climatic variables, elevation, and human activity footprints across mainland China, suitable for analyzing the interactions between agricultural productivity, environmental conditions, and human activities.
+3. **Farmland NPP and related variables** — A raster dataset capturing farmland net primary productivity (NPP), key climatic variables, elevation, and human activity footprints across mainland China, suitable for analyzing the interactions between agricultural productivity, environmental conditions, and human activities.
    *Data files*: `npp.tif`
 
 These datasets can be loaded as shown below:

vignettes/main2_ssr.Rmd

@@ -2,8 +2,8 @@
 title: "State Space Reconstruction"
 author: "Wenbo Lv"
 date: |
-  | Last update: 2025-12-15
-  | Last run: 2025-12-20
+  | Last update: 2026-01-01
+  | Last run: 2026-01-01
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteIndexEntry{2. State Space Reconstruction}
@@ -19,7 +19,7 @@ Takens theory proves that for a dynamical system $\phi$, if its trajectory conve
 
 According to the generalized embedding theorem, for a compact $d$-dimensional manifold $M$ and a set of observation functions $\langle h_1, h_2, \ldots, h_L \rangle$, the mapping $\psi_{\phi,h}(x) = \langle h_1(x), h_2(x), \ldots, h_L(x) \rangle$ is an embedding of $M$ when $L \geq 2d + 1$. Here, *embedding* refers to a one-to-one map that resolves all singularities of the original manifold. The observation functions $h_i$ can take the form of time-lagged values from a single time series, lags from multiple time series, or even completely distinct measurements. The former two are simply special cases of the third.
 
-This embedding framework can be extended to *spatial cross-sectional data*, which lack temporal ordering but are observed over a spatial domain. In this context, the observation functions can be defined using the values of a variable at a focal spatial unit and its surrounding neighbors (known as *spatial lags* in spatial statistics). Specifically, for a spatial location $s$, the embedding can be written as:
+This embedding framework can be extended to *spatial cross-sectional data*, which lack temporal ordering but are observed over a spatial domain. In this context, the observation functions can be defined using the values of a variable at a focal spatial unit and its surrounding neighbors (known as *spatial lags* in spatial statistics). Specifically, for a spatial unit $s$, the embedding can be written as:
 
 $$
 \psi(x, s) = \langle h_s(x), h_{s(1)}(x), \ldots, h_{s(L-1)}(x) \rangle,
@@ -92,7 +92,7 @@ plot3D::scatter3D(v[,4], v[,5], v[,10], colvar = NULL, pch = 19,
                   bty = "f", clab = NA, tickmarks = FALSE)
 ```
 
-![**Figure 1**. The reconstructed shadow manifolds for the variable `popd`.](../man/figures/ssr/fig1-1.png)
+![**Figure 1**. The reconstructed attractors for the variable `popd`.](../man/figures/ssr/fig1-1.png)
 
 <br>
 
@@ -138,4 +138,4 @@ plot3D::scatter3D(r[,1], r[,2], r[,3], colvar = NULL, pch = 19,
                   bty = "f", clab = NA, tickmarks = FALSE)
 ```
 
-![**Figure 2**. The reconstructed shadow manifolds for the variable `npp`.](../man/figures/ssr/fig2-1.png)
+![**Figure 2**. The reconstructed attractors for the variable `npp`.](../man/figures/ssr/fig2-1.png)

vignettes/main2_ssr.Rmd.orig

@@ -2,7 +2,7 @@
 title: "State Space Reconstruction"
 author: "Wenbo Lv"
 date: |
-  | Last update: 2025-12-15
+  | Last update: 2026-01-01
   | Last run: `r Sys.Date()`
 output: rmarkdown::html_vignette
 vignette: >
@@ -25,7 +25,7 @@ Takens theory proves that for a dynamical system $\phi$, if its trajectory conve
 
 According to the generalized embedding theorem, for a compact $d$-dimensional manifold $M$ and a set of observation functions $\langle h_1, h_2, \ldots, h_L \rangle$, the mapping $\psi_{\phi,h}(x) = \langle h_1(x), h_2(x), \ldots, h_L(x) \rangle$ is an embedding of $M$ when $L \geq 2d + 1$. Here, *embedding* refers to a one-to-one map that resolves all singularities of the original manifold. The observation functions $h_i$ can take the form of time-lagged values from a single time series, lags from multiple time series, or even completely distinct measurements. The former two are simply special cases of the third.
 
-This embedding framework can be extended to *spatial cross-sectional data*, which lack temporal ordering but are observed over a spatial domain. In this context, the observation functions can be defined using the values of a variable at a focal spatial unit and its surrounding neighbors (known as *spatial lags* in spatial statistics). Specifically, for a spatial location $s$, the embedding can be written as:
+This embedding framework can be extended to *spatial cross-sectional data*, which lack temporal ordering but are observed over a spatial domain. In this context, the observation functions can be defined using the values of a variable at a focal spatial unit and its surrounding neighbors (known as *spatial lags* in spatial statistics). Specifically, for a spatial unit $s$, the embedding can be written as:
 
 $$
 \psi(x, s) = \langle h_s(x), h_{s(1)}(x), \ldots, h_{s(L-1)}(x) \rangle,
@@ -55,7 +55,7 @@ v = spEDM::embedded(popd_sf,"popd",E = 10)
 v[1:5,c(4,5,10)]
 ```
 
-```{r fig1,fig.width=4.25,fig.height=4.5,fig.dpi=100,fig.cap=knitr::asis_output("**Figure 1**. The reconstructed shadow manifolds for the variable `popd`.")}
+```{r fig1,fig.width=4.25,fig.height=4.5,fig.dpi=100,fig.cap=knitr::asis_output("**Figure 1**. The reconstructed attractors for the variable `popd`.")}
 plot3D::scatter3D(v[,4], v[,5], v[,10], colvar = NULL, pch = 19,
                   col = "red", theta = 45, phi = 10, cex = 0.35,
                   bty = "f", clab = NA, tickmarks = FALSE)
@@ -81,7 +81,7 @@ r = spEDM::embedded(npp,"npp",E = 5,tau = 20)
 r[which(!is.na(r),arr.ind = T)[1:5],1:3]
 ```
 
-```{r fig2,fig.width=4.25,fig.height=4.5,fig.dpi=100,fig.cap=knitr::asis_output("**Figure 2**. The reconstructed shadow manifolds for the variable `npp`.")}
+```{r fig2,fig.width=4.25,fig.height=4.5,fig.dpi=100,fig.cap=knitr::asis_output("**Figure 2**. The reconstructed attractors for the variable `npp`.")}
 plot3D::scatter3D(r[,1], r[,2], r[,3], colvar = NULL, pch = 19,
                   col = "#e77854", theta = 45, phi = 10, cex = 0.01,
                   bty = "f", clab = NA, tickmarks = FALSE)

vignettes/main3_gccm.Rmd

@@ -2,8 +2,8 @@
 title: "Geographical Convergent Cross Mapping"
 author: "Wenbo Lv"
 date: |
-  | Last update: 2025-12-15
-  | Last run: 2025-12-20
+  | Last update: 2026-01-01
+  | Last run: 2026-01-01
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteIndexEntry{3. Geographical Convergent Cross Mapping}
@@ -23,58 +23,58 @@ S_{(\tau)}(x_1) & S_{(2\tau)}(x_1) & \cdots & S_{(E\tau)}(x_1) \\
 S_{(\tau)}(x_2) & S_{(2\tau)}(x_2) & \cdots & S_{(E\tau)}(x_2) \\
 \vdots & \vdots & \ddots & \vdots \\
 S_{(\tau)}(x_n) & S_{(2\tau)}(x_n) & \cdots & S_{(E\tau)}(x_n)
-\end{bmatrix}
+\end{bmatrix}.
 $$
 
 Here, $S_{(j)}(x_i)$ denotes the $j$ th-order spatial lag value of spatial unit $i$, $\tau$ is the step size for the spatial lag order, and $E$ is the embedding dimension and $M_{x}$ corresponds to the shadow manifolds of $X$.
 
 With the reconstructed shadow manifolds $M_{x}$, the state of Y can be predicted with the state of X through
 
 $$
-\hat{Y}_s \mid M_x = \sum\limits_{i=1}^k \left(\omega_{si}Y_{si} \mid M_x \right)
+\hat{Y}_s \mid M_x = \sum\limits_{i=1}^k \left(\omega_{si}Y_{si} \mid M_x \right),
 $$
 
 where $s$ represents a spatial unit at which the value of $Y$ needs to be predicted, $\hat{Y}_s$ is the prediction result, $k$ is the number of nearest neighbors used for prediction, $si$ is the spatial unit used in the prediction, $Y_{si}$ is the observation value of $Y$ at $si$. $\omega_{si}$ is the corresponding weight defined as:
 
 $$
-\omega_{si} \mid M_x = \frac{weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}{\sum_{i=1}^{L+1}weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}
+\omega_{si} \mid M_x = \frac{weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}{\sum_{i=1}^{L+1}weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)},
 $$
 
 where $\psi(M_x, s_i)$ is the state vector of spatial unit $s_i$ in the shadow manifold $M_x$, and $weight (\ast, \ast)$ is the weight function between two states in the shadow manifold, defined as:
 
 $$
 weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right) =
-\exp \left(- \frac{dis \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}{dis \left(\psi\left(M_x,s_1\right),\psi\left(M_x,s\right)\right)} \right)
+\exp \left(- \frac{dis \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}{dis \left(\psi\left(M_x,s_1\right),\psi\left(M_x,s\right)\right)} \right),
 $$
 
 where $\exp$ is the exponential function and $dis \left(\ast,\ast\right)$ represents the distance function between two states in the shadow manifold defined as:
 
 $$
-dis \left( \psi(M_x, s_i), \psi(M_x, s) \right) = \frac{1}{E} \sum_{m=1}^{E} \left| \psi_m(M_x, s_i) - \psi_m(M_x, s) \right|
+dis \left( \psi(M_x, s_i), \psi(M_x, s) \right) = \frac{1}{E} \sum_{m=1}^{E} \left| \psi_m(M_x, s_i) - \psi_m(M_x, s) \right|,
 $$
 
 where $dis \left( \psi(M_x, s_i), \psi(M_x, s) \right)$ denotes the average absolute difference between corresponding elements of the two state vectors in the shadow manifold $M_x$, $E$ is the embedding dimension, and $\psi_m(M_x, s_i)$ is the $m$-th element of the state vector $\psi(M_x, s_i)$.
 
 The skill of cross-mapping prediction $\rho$ is measured by the Pearson correlation coefficient between the true observations and corresponding predictions, and the confidence interval of $\rho$ can be estimated based the $z$-statistics with the normal distribution:
 
 $$
-\rho = \frac{Cov\left(Y,\hat{Y}\mid M_x\right)}{\sqrt{Var\left(Y\right) Var\left(\hat{Y}\mid M_x\right)}}
+\rho = \frac{Cov\left(Y,\hat{Y}\mid M_x\right)}{\sqrt{Var\left(Y\right) Var\left(\hat{Y}\mid M_x\right)}},
 $$
 
 $$
-t = \rho \sqrt{\frac{n-2}{1-\rho^2}}
+t = \rho \sqrt{\frac{n-2}{1-\rho^2}},
 $$
 
 where $n$ is the number of observations to be predicted, and
 
 $$
-z = \frac{1}{2} \ln \left(\frac{1+\rho}{1-\rho}\right)
+z = \frac{1}{2} \ln \left(\frac{1+\rho}{1-\rho}\right).
 $$
 
 The prediction skill $\rho$ varies by setting different sizes of libraries, which means the quantity of observations used in reconstruction of the shadow manifold. We can use the convergence of $\rho$ to infer the causal associations. For GCCM, the convergence means that $\rho$ increases with the size of libraries and is statistically significant when the library becomes largest.
 
 $$
-\rho_{x \to y} = \lim_{L \to \infty} cor \left( Y,\hat{Y}\mid M_x \right)
+\rho_{x \to y} = \lim_{L \to \infty} cor \left( Y,\hat{Y}\mid M_x \right),
 $$
 
 where $\rho_{x \to y}$ is the correlation after convergence, used to measure the causation effect from $Y$ to $X$, despite the notation suggesting the reverse direction.
@@ -143,7 +143,7 @@ pd_res = spEDM::gccm(data = popd_sf, cause = "pre", effect = "popd",
                      E = c(3,9), k = 12, nb = popd_nb, progressbar = FALSE)
 endTime = Sys.time()
 print(difftime(endTime,startTime, units ="mins"))
-## Time difference of 1.867319 mins
+## Time difference of 1.816188 mins
 pd_res
 ##    libsizes pre->popd  popd->pre
 ## 1       100 0.1199174 0.03313697
@@ -170,7 +170,7 @@ plot(pd_res, xlimits = c(0, 2800), draw_ci = TRUE) +
   ggplot2::theme(legend.justification = c(0.65,1))
 ```
 
-![**Figure 1**. The cross-mapping prediction outputs between population density and county-level precipitation.](../man/figures/gccm/fig1-1.png)
+![**Figure 1**. The cross-mapping results between population density and county-level precipitation.](../man/figures/gccm/fig1-1.png)
 
 <br>
 
@@ -239,7 +239,7 @@ npp_res = spEDM::gccm(data = npp, cause = "pre", effect = "npp",
                       progressbar = FALSE)
 endTime = Sys.time()
 print(difftime(endTime,startTime, units ="mins"))
-## Time difference of 0.9425755 mins
+## Time difference of 0.9297233 mins
 npp_res
 ##   libsizes  pre->npp  npp->pre
 ## 1       10 0.1235069 0.1061684
@@ -259,4 +259,4 @@ plot(npp_res,xlimits = c(5, 135),ylimits = c(0.05,1), draw_ci = TRUE) +
   ggplot2::theme(legend.justification = c(0.25,1))
 ```
 
-![**Figure 2**. The cross-mapping prediction outputs between farmland NPP and precipitation.](../man/figures/gccm/fig2-1.png)
+![**Figure 2**. The cross-mapping results between farmland NPP and precipitation.](../man/figures/gccm/fig2-1.png)

vignettes/main3_gccm.Rmd.orig

@@ -2,7 +2,7 @@
 title: "Geographical Convergent Cross Mapping"
 author: "Wenbo Lv"
 date: |
-  | Last update: 2025-12-15
+  | Last update: 2026-01-01
   | Last run: `r Sys.Date()`
 output: rmarkdown::html_vignette
 vignette: >
@@ -29,58 +29,58 @@ S_{(\tau)}(x_1) & S_{(2\tau)}(x_1) & \cdots & S_{(E\tau)}(x_1) \\
 S_{(\tau)}(x_2) & S_{(2\tau)}(x_2) & \cdots & S_{(E\tau)}(x_2) \\
 \vdots & \vdots & \ddots & \vdots \\
 S_{(\tau)}(x_n) & S_{(2\tau)}(x_n) & \cdots & S_{(E\tau)}(x_n)
-\end{bmatrix}
+\end{bmatrix}.
 $$
 
 Here, $S_{(j)}(x_i)$ denotes the $j$ th-order spatial lag value of spatial unit $i$, $\tau$ is the step size for the spatial lag order, and $E$ is the embedding dimension and $M_{x}$ corresponds to the shadow manifolds of $X$.
 
 With the reconstructed shadow manifolds $M_{x}$, the state of Y can be predicted with the state of X through
 
 $$
-\hat{Y}_s \mid M_x = \sum\limits_{i=1}^k \left(\omega_{si}Y_{si} \mid M_x \right)
+\hat{Y}_s \mid M_x = \sum\limits_{i=1}^k \left(\omega_{si}Y_{si} \mid M_x \right),
 $$
 
 where $s$ represents a spatial unit at which the value of $Y$ needs to be predicted, $\hat{Y}_s$ is the prediction result, $k$ is the number of nearest neighbors used for prediction, $si$ is the spatial unit used in the prediction, $Y_{si}$ is the observation value of $Y$ at $si$. $\omega_{si}$ is the corresponding weight defined as:
 
 $$
-\omega_{si} \mid M_x = \frac{weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}{\sum_{i=1}^{L+1}weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}
+\omega_{si} \mid M_x = \frac{weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}{\sum_{i=1}^{L+1}weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)},
 $$
 
 where $\psi(M_x, s_i)$ is the state vector of spatial unit $s_i$ in the shadow manifold $M_x$, and $weight (\ast, \ast)$ is the weight function between two states in the shadow manifold, defined as:
 
 $$
 weight \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right) =
-\exp \left(- \frac{dis \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}{dis \left(\psi\left(M_x,s_1\right),\psi\left(M_x,s\right)\right)} \right)
+\exp \left(- \frac{dis \left(\psi\left(M_x,s_i\right),\psi\left(M_x,s\right)\right)}{dis \left(\psi\left(M_x,s_1\right),\psi\left(M_x,s\right)\right)} \right),
 $$
 
 where $\exp$ is the exponential function and $dis \left(\ast,\ast\right)$ represents the distance function between two states in the shadow manifold defined as:
 
 $$
-dis \left( \psi(M_x, s_i), \psi(M_x, s) \right) = \frac{1}{E} \sum_{m=1}^{E} \left| \psi_m(M_x, s_i) - \psi_m(M_x, s) \right|
+dis \left( \psi(M_x, s_i), \psi(M_x, s) \right) = \frac{1}{E} \sum_{m=1}^{E} \left| \psi_m(M_x, s_i) - \psi_m(M_x, s) \right|,
 $$
 
 where $dis \left( \psi(M_x, s_i), \psi(M_x, s) \right)$ denotes the average absolute difference between corresponding elements of the two state vectors in the shadow manifold $M_x$, $E$ is the embedding dimension, and $\psi_m(M_x, s_i)$ is the $m$-th element of the state vector $\psi(M_x, s_i)$.
 
 The skill of cross-mapping prediction $\rho$ is measured by the Pearson correlation coefficient between the true observations and corresponding predictions, and the confidence interval of $\rho$ can be estimated based the $z$-statistics with the normal distribution:
 
 $$
-\rho = \frac{Cov\left(Y,\hat{Y}\mid M_x\right)}{\sqrt{Var\left(Y\right) Var\left(\hat{Y}\mid M_x\right)}}
+\rho = \frac{Cov\left(Y,\hat{Y}\mid M_x\right)}{\sqrt{Var\left(Y\right) Var\left(\hat{Y}\mid M_x\right)}},
 $$
 
 $$
-t = \rho \sqrt{\frac{n-2}{1-\rho^2}}
+t = \rho \sqrt{\frac{n-2}{1-\rho^2}},
 $$
 
 where $n$ is the number of observations to be predicted, and
 
 $$
-z = \frac{1}{2} \ln \left(\frac{1+\rho}{1-\rho}\right)
+z = \frac{1}{2} \ln \left(\frac{1+\rho}{1-\rho}\right).
 $$
 
 The prediction skill $\rho$ varies by setting different sizes of libraries, which means the quantity of observations used in reconstruction of the shadow manifold. We can use the convergence of $\rho$ to infer the causal associations. For GCCM, the convergence means that $\rho$ increases with the size of libraries and is statistically significant when the library becomes largest.
 
 $$
-\rho_{x \to y} = \lim_{L \to \infty} cor \left( Y,\hat{Y}\mid M_x \right)
+\rho_{x \to y} = \lim_{L \to \infty} cor \left( Y,\hat{Y}\mid M_x \right),
 $$
 
 where $\rho_{x \to y}$ is the correlation after convergence, used to measure the causation effect from $Y$ to $X$, despite the notation suggesting the reverse direction.
@@ -121,7 +121,7 @@ pd_res
 
 Visualize the result:
 
-```{r fig1,fig.width=4.5,fig.height=3.65,fig.dpi=100,fig.cap=knitr::asis_output("**Figure 1**. The cross-mapping prediction outputs between population density and county-level precipitation.")}
+```{r fig1,fig.width=4.5,fig.height=3.65,fig.dpi=100,fig.cap=knitr::asis_output("**Figure 1**. The cross-mapping results between population density and county-level precipitation.")}
 plot(pd_res, xlimits = c(0, 2800), draw_ci = TRUE) +
   ggplot2::theme(legend.justification = c(0.65,1))
 ```
@@ -176,7 +176,7 @@ npp_res
 
 Visualize the result:
 
-```{r fig2,fig.width=4.5,fig.height=3.5,fig.dpi=100,fig.cap=knitr::asis_output("**Figure 2**. The cross-mapping prediction outputs between farmland NPP and precipitation.")}
+```{r fig2,fig.width=4.5,fig.height=3.5,fig.dpi=100,fig.cap=knitr::asis_output("**Figure 2**. The cross-mapping results between farmland NPP and precipitation.")}
 plot(npp_res,xlimits = c(5, 135),ylimits = c(0.05,1), draw_ci = TRUE) +
   ggplot2::theme(legend.justification = c(0.25,1))
 ```