Comparing two spatial variables with the probability of agreement

Computing the agreement between two continuous sequences is of great interest in statistics when comparing two instruments or one instrument with a gold standard. The probability of agreement (PA) quantifies the similarity between two variables of interest, and it is useful for accounting what constitutes a practically important difference. In this article we introduce a generalization of the PA for the treatment of spatial variables. Our proposal makes the PA dependent on the spatial lag. As a consequence, for isotropic stationary and nonstationary spatial processes, the conditions for which the PA decays as a function of the distance lag are established. Estimation is addressed through a first-order approximation that guarantees the asymptotic normality of the sample version of the PA. The sensitivity of the PA is studied for finite sample size, with respect to the covariance parameters. The new method is described and illustrated with real data involving autumnal changes in the green chromatic coordinate (Gcc), an index of"greenness"that captures the phenological stage of tree leaves, is associated with carbon flux from ecosystems, and is estimated from repeated images of forest canopies.


Introduction
The comparison of two sequences is a fundamental problem in several scientific disciplines, and it can be addressed in many different ways.For example, the Student's t-test, the correlation coefficient, and the Wilcoxon two-sample rank test are three techniques that, under different assumptions, provide information about some aspects of comparisons between two independent populations (Lin et al., 2012).When the goal is to measure the agreement between two variables to validate an assay, a process, or a newly developed instrument, it may be relevant to evaluate whether its performance is concordant with other existing ones or a "gold standard" (Lin et al., 2002).
The probability of agreement ("PA") measures the level of agreement between two continuous sequences.It was first introduced in a series of papers by Nathaniel Stevens et al. (Stevens and Anderson-Cook, 2017;Stevens et al., 2017Stevens et al., , 2018Stevens et al., , 2020;;Stevens and Lu, 2020) as an alternative approach to some existing agreement measures, including the concordance correlation coefficient between measurements generated by two different methods (Lin, 1989) and its extensions.Subsequently, Leal et al. (2019) studied the PA in a context of local influence.De Castro and Galea (2021) developed Bayesian PA methods to compare measurement systems with either homoscedastic or heteroscedastic measurement errors.We note that the literature on PA is extensive and many relevant uses for it in different scenarios have been published, but there is as yet no single definition for PA.A salient example is by Ponnet et al. (2021), who adapted the probability of concordance or C-index to the specific needs of a discrete frequency and severity model that typically is used during the technical pricing of a non-life insurance product.
In this paper, we generalize the PA (sensu Lin et al., 2002;Stevens and Anderson-Cook, 2017;Stevens et al., 2017) to the case of two georeferenced sequences in the plane.This generalization makes the PA dependent on a spatial lag, similar to how the variogram and covariance functions are used in spatial statistics.Specific conditions for the variance of the difference between the two sequences are imposed so that the PA is a decreasing function when the norm of the spatial lag increases.This monotonic feature of the PA can be established easily for the bivariate Matérn and Wendland covariance functions.We then extend the PA for the case of spatiotemporal processes so that two images of the same scene taken at different times can be compared.This extension is motivated by the need to track temporal changes in spatial patterns.In our spatiotemporal extension of the PA, we consider that the process includes a temporal and, perhaps, a spatial trend, and random noise.The resulting process is nonstationary in the mean and is flexible enough to account for a number of different trends in time.
Estimation of the PA is addressed via plug-ins and the delta method, assuming that the estimates of the parameters of the covariance function exist and are asymptotically normal (Mardia and Marshall, 1984;Acosta and Vallejos, 2018).A simple expression for the asymptotic variance of the sample PA is also derived.To assess the properties of the PA for finite sample sizes, we carried out two numerical experiments: a sensitivity study that illustrates that PA decreases as a function of the norm of the spatial lag, and a Monte Carlo simulation study to estimate the PA and the relevant parameters of two spatial covariance models.Finally, we apply our spatial PA to a temporal sequence of images of a forest canopy and estimate the PA of two images taken of the same scene at different times.Images like these are used routinely to identify seasonal changes in the unfolding, maturation, and senescence (with accompanying fall colors) of individual trees and entire forest canopies, and to estimate fluxes of carbon, water vapor, and other gases between the forest and the atmosphere (e.g., Richardson et al., 2018).
In Section 2 we give some additional, albeit brief, background on PA.In Sections 3 and 4 we introduce the idea of PA for spatial processes and establish the main theoretical results for stationary processes (Section 3) and spatiotemporal processes (Section 4).Section 5 discusses estimation of PA, which is then illustrated with numerical simulation experiments (Section 6) and an empirical example (Section 7).We conclude with an outline for future research and developments in this area (Section 8).Proofs of the four theorems and two lemmas used in the paper are given in the Appendix.

Background and preliminaries
Following Leal et al. (2019), we assume that {(X 11 , X 21 ), . . ., (X 1n , X 2n )} is a random sample from a bivariate normal distribution with mean vector µ and covariance matrix Σ.Then, a method to quantify the degree of agreement between the variables X 1 and X 2 relies on the differences between their corresponding values: The probability of agreement is defined as where c denotes the maximum acceptable difference from a practical perspective.In such a case the interval (−c, c) is often called the "clinically acceptable difference" (CAD).
Because of the normality assumption, the PA as given in Equation ( 1) takes the form where Φ(•) denotes the cumulative distribution function of the standard normal, µ D = µ 1 − µ 2 , and σ 2 D = σ 11 + σ 22 − 2σ 12 .How large the PA should be to consider the variables interchangeable is up to the practitioner; Stevens et al. (2017) suggest using ψ c ≥ 0.95 as a guideline.
Under the assumption of normality, inference for µ and Σ can be addressed via maximum likelihood (ML) (Anderson, 2003, §3.2).By substituting such estimates into Equation (2), an ML estimate for ψ c , denoted by ψ c , can be obtained.Under mild assumptions, Leal et al. (2019) established the asymptotic normality of ψ c , which relies on the asymptotic distribution of the ML under normality and the delta method.It is then straightforward to estimate approximate confidence intervals and test hypotheses about ψ c .

Probability of agreement for stationary processes
In this section we introduce the PA in the context of georeferenced variables.
Let Z(s) = (X(s), Y (s)) be a bivariate second-order stationary random field with s, h ∈ R 2 , mean (µ X , µ Y ) , and covariance function where and (•) means transposition.Define the difference This difference measures the discrepancy between the processes when there is a separation vector equal to h between them.Assume that Then, the PA between processes X(s) and which assumes the form where Φ(•) is as in Equation ( 2).If we also assume isotropic processes, the probability (Equation ( 4)) can be plotted as a function of h (where • is the Euclidean norm in R 2 ) in a similar way as the covariance function is plotted for several parametric processes.In that case, the PA (Equation (1)) is obtained as a particular case of Equation (4).In fact, ψ c = ψ c (0).We first consider the Matérn covariance function (Matérn, 1986) to illustrate Equation (5).This function is widely used in spatial statistics because of its theoretical properties and its flexibility for modeling local behavior of spatial correlations (Stein, 1999;Guttorp and Gneiting, 2006).It also is used in machine learning and with neural networks (Rasmussen and Williams, 2006).The Matérn covariance function takes the form where K ν (•) is a modified Bessel function of the second kind, a > 0 is a parameter that controls the rate of decay of the correlation, and ν > 0 is the smoothing parameter that is related to the behavior of the correlation near the origin.A special case of the Matérn function is when By choosing m = 0, we have the simplest form M (h, 1/2, a) = exp(−a h ).In the sequel, when dealing with isotropic models, we will use h for h ∈ R 2 .
For a bivariate Gaussian random field the Matérn covariance function has been extended (Gneiting et al., 2010) and defined as where σ 2 X > 0, σ 2 Y > 0, and ρ XY is the co-located correlation coefficient between X(•) and Y (•).Specific conditions for the parameters are required so that the covariance model described in equations ( 8)-( 10) is positive definite.In this model we have that For illustrative purposes, consider σ X = 1, σ Y = 2, σ XY = 1.8, a XY = 2, ρ XY = 0.9 and ν = ν XY = {0.5, 1.5, 2.5}.We plot ψ c (h) versus h for h ∈ {0, 1, . . ., 15}, c = {1.5, 2, 2.5}, and using the Matérn covariance function (Figure 1 in Supplementary Material).In all cases, ψ c (h) decreases as a function of h and the curves decay more rapidly to zero as ν decreases.This is a consequence of the monotonic property of the Matérn covariance as shown in the following example.
Example 1.Let Z(s) be a bivariate second-order stationary Gaussian process with the Matérn covariance function given in equations ( 8)-( 10).Assume that ν XY = m + 1 2 and, without loss of generality, assume that the Gaussian process has mean zero and a = 1 in Equation ( 7).Then, This implies that where ϕ(•) is the probability density function of a standard normal random variable.Because the right hand side of Equation ( 11) is positive, to prove that ψ c (h) is decreasing as a function of h, it is enough to get the sign of the derivative of M (h, ν XY , a XY ) with respect to h.We write (Acosta and Vallejos, 2018) and Q is a polynomial with negative coefficients.Therefore, we have explicitly shown that ψ c (h) is a decreasing function of h.
It should be noted that ψ c ( h ) will not necessarily be a decreasing function of h .As an example, consider a bivariate process with mean (µ, µ) and a separable covariance function However, for certain parametric models, as in Example 1, ψ c ( h ) is a monotonic function.A sufficient condition for a parametric covariance model that ensures that ψ c ( h ) is a decreasing monotonic function is given by Theorem 1 (the proof of this and subsequent theorems are in the Appendix).
Theorem 2 shows that for the bivariate Matérn covariance function defined in Equations ( 8)-( 10), σ 2 D ( h ) is an increasing function of h .Theorem 2. Suppose that σ 2 D ( h ) is obtained using the Matérn covariance model and assume that ρ XY ≥ 0. Then σ D ( h ) is an increasing function of h .In consequence, the conditions of Theorem 1 are satisfied and ψ c ( h ) is a decreasing function of h .
Example 1, therefore, is a direct consequence of Theorems 1 and 2.Moreover, in Theorem 2 it is not necessary to restrict the smoothness parameter of bivariate Matérn covariance model, The Generalized Wendland family of covariance functions (Gneiting, 2002) is defined, for an integer κ > 0, as where B denotes the beta function and µ must be positive with a lower bound as given in Gneiting (2002).By continuity, for κ = 0, we have Specific conditions for µ can also be found in Bevilacqua et al. (2019).
Lemma 1.The function in Equation ( 12) is decreasing in h for all κ ≥ 0.
For a bivariate Gaussian random field the Wendland-Gneiting covariance function has been extended (Daley et al., 2015) and is defined as where σ 2 X > 0, σ 2 Y > 0, and ρ XY is the co-located correlation coefficient between X(•) and Y (•).Specific conditions for the parameters c ii , b ii , γ ii , i = 1, 2, and b 12 , c 12 , γ 12 , ν, and κ are needed so that the covariance model described in equations ( 13)-( 15) is positive definite (Daley et al., 2015).In this model we have that Theorem 3. Suppose that σ 2 D (h) is as in Equation ( 16), and assume that ρ XY ≥ 0 and c 12 ≥ 0.Then, σ D (h) is an increasing function of h.In consequence, the conditions of Theorem 1 are satisfied and ψ c ( h ) is a decreasing function of h .

Probability of agreement for spatiotemporal processes
0 is a stationary Gaussian spatiotemporal process with mean 0 and covariance function C(h, u) = Cov(Z(s, t), Z(s where the mean function µ(s, t) = F (s, t)β, with F (s, t) known and β is a vector of unknown parameters.Now, let us define the difference The quantity defined in Equation ( 17) measures the discrepancy between the process and itself for a spatial separation h and temporal separation u.Then, under the Gaussian assumption where and γ(h, u) is the spatiotemporal semivariogram (Sherman, 2011).Consequently, a natural extension of the probability of agreement between Y (s, t) and Y (s + h, t + u) is Therefore, the PA takes the form where Φ(•) is as in Equation ( 2).As an illustration, if µ(s, t) = β 0 +β 1 t (linear trend in time) and C(h, u) = σ 2 exp(− h /φ s ) exp(−|u|/φ t ) (separable exponential covariance model) the mean and variance are Thus, Equation ( 19) can be written as No matter the sign of β 1 , ψ c (h, u) decreases as u increases, which is in agreement with the fact that the PA become smaller when separation over time is enlarged.
Note that we can also use Theorems 1-3 and Equation ( 19) to consider two different spatial processes in time.In this case, the covariance function has the same structure, but The proof is virtually identical to that of Theorem 1.

Estimation
The purpose of this section is to describe the estimation of the PA defined in Equation ( 5).Here we emphasize that the variance of the difference depends on the parameters of the correlation structure, which we denote as σ D (h) = σ D (h, θ), where θ ∈ R q , q ∈ N, is a parameter vector associated with the covariance function.The next definition stresses the dependence of the PA on θ.
Definition 1. Suppose that (X(s), Y (s)) is a bivariate second-order stationary random field with s, h ∈ R 2 , mean (µ X , µ Y ) , and parametric covariance function C(h; θ).The probability of agreement between processes X(s) and Y (s + h) is defined through If µ D , and θ are estimators of µ D , and θ, respectively, obtained from the sample (X(s 1 ), is consistent and asymptotically Gaussian with Theorem 4. Suppose that (X(s), Y (s)) is a bivariate stationary Gaussian random field.As- and µ D is independent of θ.Then, ψ c (h) is consistent and asymptotically Gaussian with As a consequence of the limiting distribution established in Theorem 4, an approximate hypothesis test for the PA can be constructed.Consider the null hypothesis versus one of the following three alternative hypotheses H 1 : this hypothesis test is relevant because it compares the PA with the nominal value suggested by Stevens et al. (2017) for a fixed h.In practice, if under the conditions of Theorem 1, the test of H 0 : ψ c (0, µ D , θ) = 0.95 versus H 1 : ψ c (0, µ D , θ) < 0.95 can be considered; if H 0 is rejected, then H 0 : ψ c ( h , µ D , θ) = 0.95 is rejected for all h , because of the monotone property of the PA.
In a spatiotemporal context, denote the covariance function parameterized by θ as C(h, u; θ).Then, the PA in this case is ψ c (h, u; β, θ), similar to that defined in Equation ( 19), with

Numerical Experiments
We carried out two numerical experiments to gain more insights into the properties of the PA for finite sample sizes.The first was a sensitivity analysis that examined how variation in key parameters affected the estimation of the PA for Gaussian random fields with a specific covariance function.The second was a Monte Carlo simulation study that considered spatiotemporal processes with a linear trend and either separable or non-separable covariance structures.In all cases, the estimates were obtained using a pairwise maximum-likelihood method implemented in the R software system version 4.0.5 (R Core Team, 2022).Code is available at https://github.com/JAcosta-Hub/Comparing-two-spatial-variables-with-the-probability-of-agreement.

Gaussian spatiotemporal process
The Monte Carlo simulation study considered a spatiotemporal process as defined in §4.For the purposes of this study, we assumed that µ(s, t) = a 0 + a 1 t, and for the correlation structure, we considered both separable and non-separable cases: In the absence of a nugget, the covariance function is C(h, u) = σ 2 R(h, u).We used the GeoModels library version 1.0.0 (Beivlacqua et al., 2022) in R for the simulations, as it allowed us to simulate spatiotemporal processes with linear mean and our defined correlation structures.A regular grid of size N S × N S was considered for the spatial coordinates, N T points in time from 1 to N T with step 1; examples are shown in Figs. 2 and 3 of the Supplementary Material.
The estimates of the spatial scale (extent) parameter φ s in the covariance terms had very large ranges, extending from < 0 to greater than the maximum size of the grid.Thus, we only estimated PA for values of φ s such that 0 < φ s < the maximum size of the grid (in this case, 50).
The parameter estimates for the spatiotemporal process with a separable covariance structure (Fig. 2 in the Supplementary Material) and two different values for the trend parameter a 1 (−0.1 and 0.1) are given in Table 1; the estimators were practically unbiased and consistent.The corresponding estimates of PA for different values of h, u and c, and fixed parameters given in Table 1 are shown in Figs.4-7 of the Supplementary Material.The corresponding results from the simulations of a spatiotemporal process with a non-separable covariance structure and identical values for the trend parameter a 1 (−0.1 and 0.1) are given, respectively, in Table 1 and Figs.8-10 of the Supplementary Material, and Table 1 and Fig. 2. As with the separable case, the estimators were practically unbiased and consistent.The model with the Iacocesare covariance had a slightly better performance than the model with the exponential-separable covariance and a higher percent of valid cases used for estimating the parameters.

Motivation
Near-Earth remote sensing provides a great deal of information about ongoing environmental change and its effects on the Earth's climate (e.g., Richardson et al., 2007Richardson et al., , 2018;;Yang et al, 2013).Of particular interest are the times of year when trees in the northern hemisphere emerge from dormancy and produce new leaves ("spring green-up") and when the leaves of these same trees senesce in the fall before the trees go dormant for the winter.The growing season is the time between the spring green-up and leaf senescence in the fall, and is that time of year when forests in the northern hemisphere remove a substantial amount of carbon dioxide from the atmosphere (Barichivich et al., 2012).There is substantial evidence that because of ongoing, anthropogenically-driven climate change, on average spring green-up is occurring earlier in the year (e.g., Richardson et al., 2007;Keenan et al, 2014) and leaf senescence is occurring later in the fall (e.g., Moon, 2022).
Table 1: Parameter estimates for the spatiotemporal process with an exponential separable and a non-separable Iacocesare covariance structure.For the time trend, the cases of positive and negative slope are included."Percent valid" is the percentage of simulations for which estimates of φs were greater than zero and less than the maximum size of the grid.   1, and for Ns = 50.Note differences in range limits of the y-axis among the nine panels.

Imagery
We analyzed a series of 15 annual images of the same scene taken in mid-October from 2008-2022 (Fig. 3A).This is the time of year when leaves of deciduous trees are senescing, which leads to the spectacular display of fall colors in New England (USA), Japan, and other parts of the northern hemisphere.We chose images taken in a 3-day window (12-15 October) each year, as this represents the approximate historical "peak" of fall colors in New England.As the regional climate has warmed, however, this peak has begun to show a shift towards later dates, and identifying the rate and spatial patterning of this shift is of interest to ecologists, foresters, tourism boards, and economists (e.g., Moon, 2022).The 15 images we used were taken from the database of the PhenoCam Network,1 a network of more than 700 fixed observation sites across North America and elsewhere in the world that since 2008 has collected high-frequency and high-resolution imagery with networked digital cameras to track the timing of vegetation change (phenology) in a range of ecosystems (Seyednasrollah, 2019).We used images from the Harvard Forest, where they have been captured at 30-60-minute intervals since 2008 with a 2048 × 1636-pixel CMOS sensor in an outdoor StarDot NetCam XL 3MP camera (Richardson, 2021).

A C B D
To focus attention on the forest canopy, we clipped a 570 × 660-pixel rectangular section from each image Fig. 3B).The edges of the clipped image (black rectangle in Fig. 3A) were chosen to maximize the size of the clipped image while avoiding sky (top), wires (left) and other instrumentation on the extended boom (lower right).The high-resolution clipped image (Fig. 3B) was then downscaled ≈100-fold (to 44 × 58 pixels; Fig. 3C)) for further analysis and estimation of spatiotemporal PA.Downscaling was done by rasterizing the image using adjoining 15 × 15-pixel windows; we used the mean RGB value from these windows in the rasterized image (Fig. 3C).This downscaling was done for two reasons.First, reasonable values of h (i.e., < 15) would have been within a single leaf of the high-resolution image, and it is of more interest to look at changes among leaves and among entire trees.Second, we estimated that estimating the PA of the two high-resolution clipped images would require ≈ 6 PB of RAM, whereas the estimation of the downscaled images, which were close to the same size as those used for our simulation studies (Figs. 2 and 3 in the Supplementary Material), preserved sufficient visual differences among trees while being computationally more manageable.
Last, for each pixel, we calculated its green chromatic coordinate (G cc ), an index of "greenness" that captures the phenological stage of tree leaves, is associated with carbon flux from ecosystems, and is estimated from repeated images of forest canopies (Richardson et al., 2018).For a pixel in an RGB image, G cc is calculated as where • DN is the digital number of the Green (G), Red (R), and Blue (B) channels, respectively, assigned to each pixel in a digital image.The PhenoCam network calculates G cc for each pixel and in their "provisional" data products reports its mean, and the 50 th , 75 th , and 90 th percentiles of the G cc for the ROI of each image, and their 1-day and 3-day running means and percentiles (Richardson et al., 2018).Original, clipped, and rasterized images, and code used for rasterizing and analyzing these images are all available on GitHub https://github.com/JAcosta-Hub/Comparing-two-spatial-variables-with-the-probability-of-agreement.
We note that our clipped rectangle is different in shape, but approximately the same size, as the "region of interest" (ROI) defined and analyzed by researchers who use these images for phenological studies (the unmasked area in Fig. 3D).The ROI for each PhenoCam site is identified as the area of the image that maximizes the amount of vegetation of interest (i.e., deciduous forest at Harvard Forest) while avoiding sky, topographic features, instrumentation, other human artefacts (e.g., buildings, wires), and other areas of the image that could give seasonally biased results (e.g., soil covered by snow) (Richardson et al., 2007).Our estimates of the mean G cc calculated by the PhenoCam network for the corresponding ROI fell within the range of G cc values for each pixel of our clipped and rasterized images (Fig. 4).

Estimates
The empirical variogram of the original data showed a strong temporal dependence, so we used linear regression independent of spatial covariance to remove the effect of a deterministic trend (slope = 0.00235, intercept = 0.3854; P < 0.0001 for both).Fig. 11 in the Supplementary Material shows the marginal empirical variograms after this trend had been removed (i.e., the empirical variogram of the residuals of the simple linear regression).
To model G cc , we considered a linear temporal trend µ(s, t) = µ(t) = a 0 + a 1 t and different covariance models (separable and non-separable), each with fixed nugget effect equal to 0. The value of the objective function (log composite likelihood) for the Exponential and Iacocesare models, respectively, were 19668351.95 and 19682852.64,and the values of the associated pseudo-Akaike information criterion were, respectively, -39336693.90and -39365689.29.These results suggested a better fit to the data when using the Iacocesare non-separable covariance model.The parameter estimates (and the values we used to initialize the estimation routine) for the temporal trend and the covariance model are given in Table 2.
The practical spatial range, defined as the distance at which 95% of the sill (i.e., σ 2 in this case) is reached, will depend on the temporal separation.For the Iacocesare covariance model, this distance can be obtained by using the estimates of the parameters in this equation: Using the estimates given in Table 2 and setting u = 0 (no time lag), the estimated practical range is 12.7 pixels.initial 0.50000 0.01000 1.00000 1.00000 2.00000 6.67616 1.00000 0.00100 estimate 0.38453 0.00258 1.06297 0.95371 1.94867 6.53170 2.08316 0.00044 Finally, Fig. 5 illustrates estimates of the PA of the G cc between images as a function of spatial lag h for four different time lags (u) and four different maximum acceptable differences c.Regardless of the values of u and c, PA decreases with increasing h , but the rate of decrease declines rapidly with increasing u.For u > 1, PA is practically independent of h (and c), which we interpret to mean that the spatiotemporal processes separated by two or more years are independent and are essentially Markovian in time.

Discussion and future work
The probability of agreement has been generalized for the analysis of concordance between two georeferenced variables.This extension possesses monotonic properties-the PA declines as a function of the norm of the spatial lag-and thus quantifies the effective (or practical) spatial range.Our spatial PA is meaningful for isotropic processes.The hypothesis testing developed in Section 5 allows the estimation of the PA as a function of the spatial range, and provides a way to rule out spatial agreement when the null hypothesis is rejected for all h .Our theoretical extension of the PA also works for spatiotemporal processes with a trend, allowing for the analysis of nonstationary processes in the mean.
The monotonic properties of the PA for finite sample size were supported by Monte Carlo simulation experiments, which also showed that the parameter estimates using the composite likelihood have small bias and low variance.The value of c plays a crucial role in the estimation and the PA is sensitive to the choice of it.In practice, the value of c needs to be scaled by the square root of the sill in order to account for the scale of the data.
The application presented in Section 7 illustrated that the PA can describe the change of a spatiotemporal variable in time while accounting for spatial and temporal dependence.In this particular case, the PA also provided information about the dependence of the trend on the spatial information of the past realizations of the process; such information may be of value in modeling the trend and developing such models should be a focus of future work.Although parameter estimates were not very sensitive to the choice of the covariance function, identification of the best covariance model (e.g., separable or non-separable, and types of each) can still be improved.If the trend cannot be modeled easily with a well-known function, prior exploration will be needed to characterize a parametric function that accurately identifies observed patterns.Nonparametric trend estimation (e.g., Strandberg et al., 2019) could also be used.
The computational efficiency of the composite likelihood method used in Sections 6 and 7 deserves attention.The computational implementations used in this article worked well for images of the order of size 50 × 50; we had to rasterize and downscale the original images used in Section 7 by two orders of magnitude before we could estimate the relevant parameters.Overcoming memory limitations to enable parameter estimation from much larger images (more than 10 6 pixels) will require innovative parallelizable algorithms.
We also note that the Monte Carlo simulations and the application were developed and illustrated using spatial processes (images) defined on regular grids.However, there is no apparent reason that the proposed methods could not also be applied to irregularly spaced spatial data.
Finally, a natural but unexplored extension of the present work is the definition of the PA for spatiotemporal marked point processes.Because the randomness in point processes is in the location (as well as in the marks) and repeated observations of a given point process will yield a new set of locations, it is of theoretical interest to estimate the PA of point patterns that evolve through time.Such estimation would have immediate applicability to ecological spatial datasets, which are predominantly samples of point processes, not rasters (e.g., Plant, 2019)

Proof of Theorem 3
Without loss of generality, we assume b 12 = 1 and note that σ 2 D (h) is an increasing function of h if and only if GW(h; κ, µ) is a decreasing function in h for all κ.Therefore, the result holds by Lemma 1, since ν + γ 12 + 1 > 0.

Proof of Lemma 2
Let σ D (h; θ) = g( θ).Applying a Taylor expansion of order 1 for g( θ) around θ, we have that σ Then, Applying expected value and variance in both sides of Equation ( 21), the result follows.Figure 4: Estimates of the probability of agreement as a function of h , u and c for a spatiotemporal Gaussian process with a negative linear trend and an exponential separable covariance structure with fixed parameters given in Table 1 (in the manuscript), and for NS = 20.Note differences in range limits of the y-axis among the nine panels.
Figure 5: Estimates of the probability of agreement as a function of h , u and c for a spatiotemporal Gaussian process with a positive linear trend and an exponential separable covariance structure with fixed parameters given in Table 1 (in the manuscript), and for NS = 20.Note differences in range limits of the y-axis among the nine panels.
Figure 6: Estimates of the probability of agreement as a function of h , u and c for a spatiotemporal Gaussian process with a negative linear trend and an exponential separable covariance structure with fixed parameters given in Table 1 (in the manuscript), and for NS = 50.Note differences in range limits of the y-axis among the nine panels.
Figure 7: Estimates of the probability of agreement as a function of h , u and c for a spatiotemporal Gaussian process with a positive linear trend and an exponential separable covariance structure with fixed parameters given in Table 1 (in the manuscript), and for NS = 50.Note differences in range limits of the y-axis among the nine panels.
Figure 8: Estimates of the probability of agreement as a function of h, u and c for a spatiotemporal Gaussian process with a negative linear trend and a non-separable Iacocesare covariance structure with fixed parameters given in Table 1 (in the manuscript), and for NS = 20.Note differences in range limits of the y-axis among the nine panels.
Figure 9: Estimates of the probability of agreement as a function of h , u and c for a spatiotemporal Gaussian process with a negative linear trend and a non-separable Iacocesare covariance structure with fixed parameters given in Table 1 (in the manuscript), and for NS = 20.Note differences in range limits of the y-axis among the nine panels.
Figure 10: Estimates of the probability of agreement as a function of h , u and c for a spatiotemporal Gaussian process with a negative linear trend and a non-separable Iacocesare covariance structure with fixed parameters given in Table 1 (in the manuscript), and for NS = 50.Note differences in range limits of the y-axis among the nine panels.

Figure 1 :
Figure 1: The effect of variation in ρXY (top left), µD (top right), φXY (bottom left), and σY (bottom right), on the behavior of PA as a function of lag h in a bivariate Gaussian random field with a Matérn covariance structure.Note differences in range limits of the y-axis among the four panels.

Figure 2 :
Figure2: Estimates of the probability of agreement as a function of h , u and c for a spatiotemporal Gaussian process with a positive linear trend and a non-separable Iacocesare covariance structure with fixed parameters given in Table1, and for Ns = 50.Note differences in range limits of the y-axis among the nine panels.

Figure 3 :
Figure 3: A. The PhenoCam image taken by a stationary camera from the EMS tower at Harvard Forest, Massachusetts, USA on 15 October 2008 at 13:31 (UTC −4).A rectangular section (570 × 660 pixels) of the image was clipped (B; black outline in A) and then down-scaled (to 44 × 58) pixels (C) for estimation of the spatial PA.Note that our rectangular image is different from the clipped "region of interest" (ROI; unmasked area in D) analyzed by the PhenoCam network.

Figure 4 :
Figure 4: Green chromatic indices (Gcc) of each of the rasters in the 15 downscaled images (grey symbols) and the mean Gcc estimated for the entire region of interest (ROI; masked area in Fig. 3D) by the PhenoCam network.Also shown are the slopes and intercepts of the temporal trend in the Gcc assuming independence (red line), separable (green), and non-separable (blue) spatial covariances.

Figure 5 :
Figure 5: Probability of Agreement as a function of h for four time lags u (years between images (individual panels), each with different of maximum acceptable differences c (black, red, blue and green curves within each graph).

Figure 2 :
Figure 2: Simulated realizations of a spatiotemporal process defined by a Gaussian random field with an exponential separable covariance function for NS = 50 and NT = 6.

Figure 3 :
Figure 3: Simulated realization of a spatiotemporal process defined by a Gaussian random field with an Iacosecare non-separable covariance function for NS = 50 and NT = 6.

Figure 11 :
Figure 11: Empirical description of the data set.

Table 2 :
Parameter estimates for the linear temporal trend and the Iacocesare covariance structure.