gls_spatial: Run spatial GLS models with optimal spatial correlation...

In natedemaagd/LandCover: Model the Spread of Landcovers and Effect on Water Yield

Description

This function runs a GLS model on a spatial dataset. It chooses among five types of correlation structure, and returns the model with the lowest AIC.

This function runs a GLS model on a spatial dataset. It chooses among five types of correlation structure, and returns the model with the lowest AIC.

Usage

gls_spatial(
  data,
  landcover_varname,
  landcover_vec,
  reg_formula,
  error_formula,
  num_cores = parallel::detectCores() - 1,
  maxIter = 50,
  msMaxIter = 50,
  tolerance = 1e-06,
  silent = FALSE
)

gls_spatial(
  data,
  landcover_varname,
  landcover_vec,
  reg_formula,
  error_formula,
  num_cores = parallel::detectCores() - 1,
  maxIter = 50,
  msMaxIter = 50,
  tolerance = 1e-06,
  silent = FALSE
)

Arguments

data	data.frame with spatial data
landcover_varname	character string specifying the landcover variable from data
landcover_vec	vector of all landcover types in data$landcover_varname to be analyzed. This should include the landcover that is spreading along with all landcovers that are susceptible to the spread.
reg_formula	regression formula to be used, as in lm(), with columns from data (e.g. c ~ a + b)
error_formula	one-sided formula specifying coordinate columns from data (e.g. ~ x + y)
num_cores	numerical. Number of cores for parallel processing, max 5
maxIter	numerical. Max iterations for GLS optimization, see ?nlme::glsControl()
msMaxIter	numerical. Max iterations for optimization step inside GLS optimization, see ?nlme::glsControl()
tolerance	numerical. Tolerance for the convergence criterion, see ?nlme::glsControl()
silent	logical. Hides notifications for completion of each landcover type. Default is FALSE

Details

Landcovers with large numbers of observations will take a long time to run. It is strongly recommended to use 5 cores in parallel. The code runs the landcover types sequentially, and the spatial GLS models for each landcover in parallel. The function will output the landcover and time of completion as it progresses, unless silent = TRUE. If there are too many observations for your system, you must take a representative subset. Each list element in the output corresponds to a landcover type. The objects are of class gls and will report the correlation structure that resulted in the lowest model AIC. If an element in the output list is NA, the log likelihood model failed to converge using all correlation structures. Try increasing the number of iterations or increasing the tolerance.

Landcovers with large numbers of observations will take a long time to run. It is strongly recommended to use 5 cores in parallel. The code runs the landcover types sequentially, and the spatial GLS models for each landcover in parallel. The function will output the landcover and time of completion as it progresses, unless silent = TRUE. If there are too many observations for your system, you must take a representative subset. Each list element in the output corresponds to a landcover type. The objects are of class gls and will report the correlation structure that resulted in the lowest model AIC. If an element in the output list is NA, the log likelihood model failed to converge using all correlation structures. Try increasing the number of iterations or increasing the tolerance.

Value

A list of two objects. The first list element is a list of regression results for each landcover type. The second list element is a residual density plot, separated by landcover type. See details.

A list of two objects. The first list element is a list of regression results for each landcover type. The second list element is a residual density plot, separated by landcover type. See details.

Examples

# load packages
library(LandCover); library(foreach)

# initialize data.frame with coordinates
dat <- expand.grid(x = 1:20, y = 1:20, KEEP.OUT.ATTRS = FALSE)


# create some data: elevation, landcover, and temp/ET dependent on elevation and landcover
dat$elevation <- with(dat, 50 + 2*x + 5*y + rnorm(nrow(dat), sd = 7))
dat$landcover <- ifelse(dat$elevation < median(dat$elevation), 1, 2)
dat[dat$x < median(dat$x) & dat$landcover == 2, 'landcover'] <- 3
dat$temp      <- with(dat, (120-0.7*(0.5*elevation + 0.3*y - 0.5*x + ifelse(landcover == 'lc1', -30, 0) + rnorm(nrow(dat)))))
dat$ET        <- with(dat, (   -0.4*(-2*temp       + 0.5*y - 1.0*x + ifelse(landcover == 'lc1', +20, 0) + rnorm(nrow(dat)))))

# run the gls model
regression_results <- gls_spatial(data = dat, landcover_varname = 'landcover', landcover_vec = c(1,2),
                                  reg_formula = ET ~ elevation + temp, error_formula = ~ x + y)

# load packages
library(LandCover); library(foreach)

# initialize data.frame with coordinates
dat <- expand.grid(x = 1:20, y = 1:20, KEEP.OUT.ATTRS = FALSE)


# create some data: elevation, landcover, and temp/ET dependent on elevation and landcover
dat$elevation <- with(dat, 50 + 2*x + 5*y + rnorm(nrow(dat), sd = 7))
dat$landcover <- ifelse(dat$elevation < median(dat$elevation), 1, 2)
dat[dat$x < median(dat$x) & dat$landcover == 2, 'landcover'] <- 3
dat$temp      <- with(dat, (120-0.7*(0.5*elevation + 0.3*y - 0.5*x + ifelse(landcover == 'lc1', -30, 0) + rnorm(nrow(dat)))))
dat$ET        <- with(dat, (   -0.4*(-2*temp       + 0.5*y - 1.0*x + ifelse(landcover == 'lc1', +20, 0) + rnorm(nrow(dat)))))

# run the gls model
regression_results <- gls_spatial(data = dat, landcover_varname = 'landcover', landcover_vec = c(1,2),
                                  reg_formula = ET ~ elevation + temp, error_formula = ~ x + y)

natedemaagd/LandCover documentation built on April 1, 2021, 4:14 p.m.