################################################################################
### Appendix 3. Beale et al. 2007.                                           ###
###                                                                          ###
### R code to run the simulations described in Beale et al. 2007.            ###
###                                                                          ###
### To run this code you will first need to download and install R from a    ###
### suitablemirror via http://www.r-project.org/ and will also need to       ###
### download and install the packages RandomFields and nlme. Functions       ###
### written specially for this analysis are provided separately and are      ###
### sourced within this code.                                                ###
################################################################################

### First, make the specially written functions and those within RandomFields
## and nlme available:
library(RandomFields)
library(nlme)
source("D:/My~path/SpatialACFunctions.R")                                        # Change the path to an appropriate location.

### Set the random number seed:
set.seed(1)

### Generate the first set of Gaussian Random Fields for analysis:
x1 <- GaussRF(x=1:32, y=1:32, model="stable", grid=TRUE,                         # fine scale AC.
              param=c(0, 10, 0, 10, 1))
x2 <- GaussRF(x=1:32, y=1:32, model="stable", grid=TRUE,                         # moderate scale AC.
              param=c(0, 4, 0, 50, 1))
x3 <- GaussRF(x=1:32, y=1:32, model="stable", grid=TRUE,                         # large scale AC.
              param=c(0, 10, 0, 200, 1))
y1 <- GaussRF(x=1:32, y=1:32, model="stable", grid=TRUE,                         # dependant with AC, no true relationship with
              param=c(0, 4, 0, 50, 1))                                           #    x1, x2, or x3.

### Fit the first OLS model using all the data:
full.lm1 <- lm(as.vector(y1) ~ as.vector(x1) + as.vector(x2) + as.vector(x3))

### Generate the setup parameters for calculating Moran's I using the same
## number of bins as Hawkins et al. 2007:
setup <- matrixMoranSetUp(mat = x1, equal.lag.dist = FALSE, nbins = 21)

### Calculate Moran's I for residuals for the residuals of the full OLS model:
MM1 <- matrixMoran(mat = matrix(resid(full.lm1), ncol = 32),  nbins = 21,
              plot.it = FALSE, resids = TRUE, useStartUp = setup)

### Identify at which distance Moran's I first becomes negative - the distance
## at which there remains no spatial autocorrelation in the residuals.
first0 <- MM1[[1]][,2] < 0                                                       # identify negative Moran's I values.
first0 <- cumsum(first0)                                                         # calculate how many negative I values have
                                                                                 #     occurred by each distance class.
first0 <- min(which(first0 == 1))                                                # find the first negative Moran's I value.
AC1 <- ceiling(32*MM1[[1]][first0,1])                                            # estimate the minimum distance needed to avoid
                                                                                 #     autocorrelation in the residuals.
### Do the subsampling routine and fit OLS models for 500 subsamples.
sub.results1 <- matrix(0, 500,4)                                                 # create matrix for subsample results.
for (i in 1:500) {
    hg1 <- hexGrid(ACsize = AC1)                                                 # identify a random sample of squares at least
                                                                                 #    AC1 distance apart.
    sub.results1[i,] <- coef(lm(as.vector(y1[hg1]) ~ as.vector(x1[hg1]) +        # fit subsampled OLS model.
           as.vector(x2[hg1]) + as.vector(x3[hg1])))
}
sub.lm1 <- colMeans(sub.results1)                                                # calculate mean of subsample coefficients.

### Build GLS model for full dataset:
x <- rep(1:32, each = 32)                                                        # generate x coordinates
y <- rep(1:32, 32)                                                               # generate y coordiantes
dat1 <- as.data.frame( cbind (as.vector(y1), as.vector(x1), as.vector(x2),       # create data.frame of covariates, explanatory
          as.vector(x3), x, y))                                                  #    variable and coordinates.
names(dat1)[1:4] <- c("y1", "x1", "x2", "x3")
gls1 <- gls(y1 ~ x1 + x2 + x3, correlation = corExp(form = ~ x + y), data =      # fit gls model with exponential spatial
          dat1)                                                                  #    autocorrelation structure.

save.image("D:/SpatialAC1.Rdata")                                                # Save results to a file.

### Repeat above procedure a further 999 times ommiting the first (slow)
## calculations of the setup parameters for the Moran's I test:
for(ii in 2:1000) {

    set.seed(ii)
    x1 <- GaussRF(x=1:32, y=1:32, model="stable", grid=TRUE,
                  param=c(0, 4, 0, 50, 1))
    x2 <- GaussRF(x=1:32, y=1:32, model="stable", grid=TRUE,
                  param=c(0, 10, 0, 10, 1))
    x3 <- GaussRF(x=1:32, y=1:32, model="stable", grid=TRUE,
                  param=c(0, 10, 0, 200, 1))
    y1 <- GaussRF(x=1:32, y=1:32, model="stable", grid=TRUE,
                  param=c(0, 4, 0, 50, 1))

    full.lm1 <- lm(as.vector(y1) ~ as.vector(x1) + as.vector(x2) +
                      as.vector(x3))

    MM1 <- matrixMoran2(mat = matrix(resid(full.lm1), ncol = 32),  nbins = 21,
                  plot.it = FALSE, resids = TRUE, useStartUp = setup)

    first0 <- MM1[[1]][,2] < 0
    first0 <- cumsum(first0)
    first0 <- min(which(first0 == 1))
    AC1 <- ceiling(32*MM1[[1]][first0,1])

    sub.results1 <- matrix(0, 500,4)
    for (i in 1:500) {
        hg1 <- hexGrid(ACsize = AC1)
        sub.results1[i,] <- coef(lm(as.vector(y1[hg1]) ~ as.vector(x1[hg1]) +
               as.vector(x2[hg1]) + as.vector(x3[hg1])))
    }
    sub.lm1 <- colMeans(sub.results1)
    dat1 <- as.data.frame(cbind(as.vector(y1), as.vector(x1), as.vector(x2),
        as.vector(x3), x, y))
    names(dat1)[1:4] <- c("y1", "x1", "x2", "x3")
    gls1 <- try(gls(y1 ~ x1 + x2 + x3, correlation = corExp(form = ~ x + y),
        data = dat1), TRUE)
    if(class(gls1) == 'try-error') gls1 <- try(gls(y1 ~ x1 + x2 + x3,            # if REML model doesn't converge, use
        correlation = corExp(form = ~ x + y), data = dat1, method = "ML"), TRUE) #    Maximum Likelihood instead.

    fil <- paste("D:/SpatialAC", ii, ".Rdata", sep = "")
    save.image(file = fil)

    print(paste(ii, "completed"))
}