## Goal of this document: Study how ASReml-R splits up the error variance between spatially correlated error variance and spatially uncorrelated error variance (i.e. "nugget variance" in geostatistics).
## Load packages
library(RandomFields)
library(LaplacesDemon)
library(asreml)

## Set random seeed
RFoptions(seed=0)

## Create dataset with 50 row and 100 columns
data <- expand.grid(row=1:50,col=1:100)
data$obs <- as.factor(1:nrow(data))

## Define parameters
intercept <- 3
siguncorres <- sqrt(9)
sigcorres <- sqrt(25)
rho=0.5

## Initialize AR1 simulation model
model <- RMexp(scale=-1/log(rho), var=sigcorres^2)

## Add spatially correlated and spatially uncorrelated error
autocorrerror <- RFsimulate(model,x=data$row,y=data$col)@data[[1]]
uncorrerror <- rnorm(nrow(data),0,siguncorres)
data$y <- intercept+autocorrerror+uncorrerror

## Convert as factor for asreml
data$row <- as.factor(data$row)
data$col <- as.factor(data$col)

## Dataset should be with columns nested within rows
data <- data[order(data$row),]
m0 <- asreml(y ~ 1,rcov=~ar1(row):ar1(col),random=~obs,na.method.X="include",data=data,trace=F)
m0 <- update(m0)
## Equivalent model with the exponential function
#m1 <- asreml(y ~ 1,rcov=~exp(row):exp(col),random=~obs,na.method.X="include",data=data,trace=F)

m <- m0
m$coefficients$fixed
summary(m)$varcomp
## Check 95% confidence intervals
c(summary(m)$varcomp[1,2]-qnorm(0.975)*summary(m)$varcomp[1,3],summary(m)$varcomp[1,2]+qnorm(0.975)*summary(m)$varcomp[1,3])
c(summary(m)$varcomp[2,2]-qnorm(0.975)*summary(m)$varcomp[2,3],summary(m)$varcomp[2,2]+qnorm(0.975)*summary(m)$varcomp[2,3])

## Check variogram
plot.asrVariogram(variogram(m))