Distance

The first step is to calculate the distances among pixels as a distance matrix D. Here, we’ll calculate relative distances with distm_scaled() from our coordinate matrix.

D <- distm_scaled(coords)
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)

distm_scaled() scales distances across the spatial domain so that the greatest distance between two pixels is 1. Note that because distm_scaled() uses geosphere::distm(), it treats coordinates as degrees longitude and latitude on a sphere and calculates distances accordingly.

Covariance

Next, we need to estimate the expected correlation between the random errors of our spatial response variable (estimates of \(\beta_1\)) based on their distances. To do so, we need a spatial covariance function. In this example, we will use an exponential covariance function to estimate correlations: \(V = \exp\big(\frac{-D}{r}\big)\) where \(r\) is a parameter that dictates the range of spatial autocorrelation. The function covar_exp() corresponds to this covariance function.

V represents the correlation among points if all variation is accounted for. However, it is safest to assume that there is some additional source of unexplained and unmeasured variance (a nugget \(\eta\)). Therefore, we assume that the covariance structure among pixels is given by \(\Sigma = \eta I + (1-\eta)V\) where \(I\) is the identity matrix.

If we know the range parameter \(r\), we can calculate V from D with covar_exp():

r <- 0.1
V <- covar_exp(D, r)

And we could add a known nugget to V to obtain Sigma:

nugget <- 0.2 
I <- diag(nrow(V)) # identity matrix
Sigma <- nugget*I + (1-nugget)*V

See ?covar_exp() for a description of the covariance functions provided by remotePARTS and for more information regarding their use.

Known parameters

If the parameters that govern the spatial autocorrelation \(\Sigma_\gamma\) are known, a GLS can be fit with fitGLS(). Here, we will fit the GLS by providing (i) a model formula, (ii) a data source, (iii) our V matrix, which was pre-calculated with a spatial parameter \(r = 0.1\) and (iv) a nugget of \(\eta = 0.2\).

The specific task in this examples is to estimate the effect of land-cover class on our time trend. Because land is a factor, we’ll also specify a no-intercept model.

Note that the GLS fitting process requires an inversion of V. This means that even with only the 3000-pixel subset, it will take a few minutes to finish the computations on most computers.

GLS.0 <- fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = V, nugget = nugget)

Note also that fitGLS adds the nugget to V internally. If we wanted to do this ourselves, we could pass the covariance matrix Sigma which already contains a nugget component and then set the nugget argument of fitGLS to 0:

fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = Sigma, nugget = 0) # equivalent

The estimates for our land class effects can be extracted with coefficients().

coefficients(GLS.0)
## landShrubland   landSavanna landGrassland 
##  0.0008646565  0.0002554862  0.0011294576

The full model fit is given by

GLS.0
## 
## Call:
## fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = V, 
##     nugget = nugget)
## 
## t-tests:
##                     Est t.stat pval.t
## landShrubland 0.0008647 0.5065 0.6126
## landSavanna   0.0002555 0.1493 0.8813
## landGrassland 0.0011295 0.6648 0.5062
## 
## F-test:
##                model df_F    SSE       MSE logLik Fstat  pval_F
## 1 AR_coef ~ 0 + land    2 0.1128 3.764e-05  12808 4.461 0.01163
## 2        AR_coef ~ 1 2997 0.1131 3.775e-05  12802    NA      NA

Note that, although none of the three land-cover classes shows a statistically significant time trend in (proportional) NDVI, the F-test shows that there is a statistically significant difference among land-cover classes, because the model including land-cover classes gives a statistically significantly better fit to the data than the intercept-only model.

Parameter estimation

In practice, we rarely know the values of the parameters that govern spatial autocorrelation (e.g., the range and nugget) in advance. Therefore, these parameters will need to be estimated for most data.

corfit <- fitCor(resids = residuals(ARfit), coords = coords, covar_FUN = "covar_exp", 
                 start = list(range = 0.1), fit.n = 3000)
(range.opt = corfit$spcor)
##     range 
## 0.1083729

max.dist <- max(distm_km(coords))
corfit.km <- fitCor(resids = residuals(ARfit), coords = coords, covar_FUN = "covar_exp", 
                    start = list(range = max.dist*0.1), distm_FUN = "distm_km", fit.n = 3000)

V.opt <- covar_exp(D, range.opt)

GLS.opt <- fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = V.opt, nugget = NA, no.F = FALSE)
(nug.opt = GLS.opt$nugget)
## [1] 0.1342314
coefficients(GLS.opt)
## landShrubland   landSavanna landGrassland 
##  0.0009201230  0.0003899361  0.0011467324

rbind(GLS.0 = c(range = r, nugget = GLS.0$nugget, logLik = GLS.0$logLik, MSE = GLS.0$MSE),
      GLS.opt = c(range = range.opt, nugget = GLS.opt$nugget, logLik = GLS.opt$logLik, MSE = GLS.opt$MSE))
##             range    nugget   logLik          MSE
## GLS.0   0.1000000 0.2000000 12807.91 3.763878e-05
## GLS.opt 0.1083729 0.1342314 12809.46 5.003165e-05

fitopt <- fitGLS_opt(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, 
                     coords = ndvi_AK3000[, c("lng", "lat")], 
                     covar_FUN = "covar_exp", 
                     start = c(range = .1, nugget = .2),
                     method = "BFGS", # use BFGS algorightm (see ?stats::optim())
                     control = list(reltol = 1e-5)  # lower the convergence tolerance (see ?stats::optim()) 
                     )
fitopt$opt$par
#      range     nugget 
# 0.02497874 0.17914929 
fitopt$GLS$logLik
# [1] 12824.77
fitopt$GLS$MSE
# [1] 2.475972e-05

Intercept-only model

If we want to test the hypothesis that “there was a trend in Alaska NDVI Alaska from 1982-2013”, we can regress the AR coefficient on an intercept-only GLS model:

(GLS.int <- fitGLS(AR_coef ~ 1, data = ndvi_AK3000, V = V.opt, nugget = nug.opt, no.F = TRUE))
## 
## Call:
## fitGLS(formula = AR_coef ~ 1, data = ndvi_AK3000, V = V.opt, 
##     nugget = nug.opt, no.F = TRUE)
## 
## t-tests:
##                  Est t.stat pval.t
## (Intercept) 0.000972 0.4564 0.6481
## 
## Model statistics:
##      SSE       MSE logLik
## 1 0.1503 5.011e-05  12806

We can see from the t-test that the intercept is not statistically different from zero. In other words, there is no map-level temporal trend in NDVI across the entire data set. We have not performed an F-test, because the full model is the intercept-only model and is therefore the same as the reduced model.

Land-cover effects

If we want to test the hypothesis that “trends in Alaskan NDVI differ by land- cover class”, we can use GLS.opt() from earlier:

GLS.opt
## 
## Call:
## fitGLS(formula = AR_coef ~ 0 + land, data = ndvi_AK3000, V = V.opt, 
##     nugget = NA, no.F = FALSE)
## 
## t-tests:
##                     Est t.stat pval.t
## landShrubland 0.0009201 0.4306 0.6668
## landSavanna   0.0003899 0.1822 0.8554
## landGrassland 0.0011467 0.5385 0.5903
## 
## F-test:
##                model df_F    SSE       MSE logLik Fstat  pval_F
## 1 AR_coef ~ 0 + land    2 0.1499 5.003e-05  12809 3.263 0.03841
## 2        AR_coef ~ 1 2997 0.1503 5.014e-05  12805    NA      NA

The t-tests show that the trend in NDVI, for all land-cover classes, was not statistically different from zero, meaning that NDVI did not show a statistically significant trend in any land-cover class. The F-test (ANOVA table), however, shows that time trends in NDVI differ among the land-cover classes. The better fit of the model with land-cover classes can also be seen in the increase in the likelihood (logLik) compared to the intercept-only model.

Latitude effects

Finally, to test the hypothesis that “temporal trends in NDVI differ with latitude”, we can regress the AR coefficient on latitude in our GLS model:

(GLS.lat <- fitGLS(AR_coef ~ 1 + lat, data = ndvi_AK3000, V = V.opt, nugget = nug.opt, no.F = FALSE))
## 
## Call:
## fitGLS(formula = AR_coef ~ 1 + lat, data = ndvi_AK3000, V = V.opt, 
##     nugget = nug.opt, no.F = FALSE)
## 
## t-tests:
##                    Est   t.stat pval.t
## (Intercept) -8.617e-04 -0.04362 0.9652
## lat          2.986e-05  0.09337 0.9256
## 
## F-test:
##               model df_F    SSE       MSE logLik    Fstat pval_F
## 1 AR_coef ~ 1 + lat    1 0.1503 5.012e-05  12806 0.008719 0.9256
## 2       AR_coef ~ 1 2998 0.1503 5.012e-05  12806       NA     NA

The t-tests show that temporal trends in NDVI did not differ with latitude. Note that the p-value from the F-test is equivalent to that of the t-test p-value for effect of latitude.

Conclusions (ndvi_AK3000)

We can see from these hypothesis tests that, at least among the 3000-pixel sub-sample of Alaska, the answer to all three questions that we posed is no: there is no statistical evidence for an overall greening in Alaska, nor differences among land-cover classes or latitude.

Partitioned GLS

fitGLS_parition() conducts GLS analyses on partitions of the data and then combines the results from the partitions to give overall statistical results. Specifically, this process (1) breaks the data into smaller and more manageable pieces (partitions), (2) conducts GLS on each partition, (3) calculates cross-partition statistics from pairs of partitions, and (4) summarizes the results with statistical tests that account for correlations among partitions. We will use the full ndvi_AK10000 data set to demonstrate fitGLS_parition().

We have already performed the time-series analyses on the full data set so you don’t have to. These are in the AR_coef column of ndvi_AK10000. However, we used a complete data set, so you will need to remove rare land-cover classes.

df = ndvi_AK10000

Step (1) is to divide pixels up into partitions, which is done with the function sample_parition(). Passing sample_partitions, the number of pixels in our map, and the argument partsize = 1500 will result in a partition matrix with 1,500 rows and 20 columns. Columns of the resulting partition matrix pm each contain a random sample of 1,500 pixels. Each of these 20 samples (partitions) are non-overlapping, containing no repeats. Setting npart = NA will automatically give the maximum number of partitions possible (i.e., 20).

pm <- sample_partitions(npix = nrow(df), partsize = 1500, npart = NA)
dim(pm)
## [1] 1500    6

Once we have our partition matrix, fitGLS_parititon() performs steps (2)-(4) of the analyses. The input is similar to fitGLS. For this example, we specify (i) a formula, (ii) the data as a data.frame (df), (iii) the partition matrix (partmat pm), (iv) the covariance function (covar_FUN), (v) a list of spatial parameters covar.pars including our optimized range parameter, and a nugget. If nugget is specified, this value is used for the calculations, while if nugget = NA (the default) it is estimated for each partition separately.

Note that, although the compute time is much faster than if we needed to invert the full covariance matrix V, this example still takes a long time to fit. Therefore, we have saved the output of this code as an R object partGLS.ndviAK so that you can look at its output without having to execute the function.

The model was fit with this code:

partGLS_ndviAK <- fitGLS_partition(formula = AR_coef ~ 0 + land, data = df, 
                                   partmat = pm, covar_FUN = "covar_exp", 
                                   covar.pars = list(range = range.opt),
                                   nugget = nug.opt, ncores = 4)

and can be loaded with

data(partGLS_ndviAK)

Here are the t-tests, that show that land cover class does not significantly affect NDVI trend:

partGLS_ndviAK$overall$t.test
## $p.t
##                        Est          SE    t.stat    pval.t
## landShrubland 0.0013094359 0.001899589 0.6893260 0.4906359
## landSavanna   0.0004379666 0.001902064 0.2302586 0.8178960
## landGrassland 0.0016210956 0.001894683 0.8556027 0.3922404
## 
## $covar_coef
##              [,1]         [,2]         [,3]
## [1,] 3.623978e-06 3.569499e-06 3.561642e-06
## [2,] 3.566697e-06 3.633439e-06 3.537596e-06
## [3,] 3.562415e-06 3.540957e-06 3.605232e-06

It is highly recommended that users read the full documentation (?fitGLS_parition()) before using fitGLS_partition to analyze any data.

Here is the p-value for the chisqr test of the partitioned GLS

partGLS_ndviAK$overall$pval.chisqr
## pval.chisqr 
##       1e-06 
## attr(,"rankdef.MSR")
## [1] 0

This again, indicates that the model which includes land cover classes better explains variation in NDVI trends than the intercept-only model. Note that the p-value is much lower than the p-value from the F-test conducted by fitGLS(). This is likely due to outliers in the data, which should be removed before conducting any real analysis. One simple way to filter potential trend outliers would be to remove any pixels whose time-series coefficient standard errors are unusually large (e.g., SE > 4*mean(SE)).

Conclusions

Our parititoned GLS has returned the same conclusions as our standard GLS analysis. No map-level trend in NDVI was found within any of the land cover classes.

This was only one of our three hypotheses tested earlier. The remaining two can easily be tested with fitGLS_partition(), as they were with fitGLS, by changing the formula argument.