Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- ## 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))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement