rst | R Documentation |
This function simulates spatial and temporal random effects with mean zero. The method is described in Algorithm 3.1 of Rue & Held 2015.
rst( n = 1, type = c("s", "t", "st")[1], type.s = "ICAR", type.t = c("RW1", "RW2")[2], Amat = NULL, n.t = NULL, scale.model = TRUE )
n |
sample size |
type |
type of random effects: temporal (t), spatial (s), or spatial-temporal (st) |
type.s |
type of spatial random effect, currently only ICAR is available |
type.t |
type of temporal random effect, currently only RW1 and RW2 are available |
Amat |
adjacency matrix for the spatial regions |
n.t |
number of time points for the temporal random effect |
scale.model |
logical indicator of whether to scale the random effects to have unit generalized variance. See Sørbye 2013 for more details |
a matrix (for spatial or temporal) or a three-dimensional array (for spatial-temporal) of the random effects.
Zehang Richard Li
Rue, H., & Held, L. (2005). Gaussian Markov random fields: theory and applications. CRC press.
Sørbye, S. H. (2013). Tutorial: Scaling IGMRF-models in R-INLA. Department of Mathematics and Statistics, University of Tromsø.
## Not run: data(DemoMap) ## Spatial random effects out <- rst(n=10000, type = "s", Amat = DemoMap$Amat) # To verify the mean under the conditional specification mean(out[,1] - apply(out[,c(2,3,4)], 1, mean)) mean(out[,2] - apply(out[,c(1,3)], 1, mean)) mean(out[,3] - apply(out[,c(1,2,4)], 1, mean)) mean(out[,4] - apply(out[,c(1,3)], 1, mean)) ## Temporal random effects (RW1) out <- rst(n=1, type = "t", type.t = "RW1", n.t = 200, scale.model = FALSE) par(mfrow = c(1,2)) plot(1:dim(out)[2], out, col = 1, type = "l", xlab = "Time", ylab = "Random effects") # verify the first order difference is normally distributed first_diff <- diff(as.numeric(out[1,])) qqnorm(first_diff ) abline(c(0,1)) ## Temporal random effects (RW2) out <- rst(n=1, type = "t", type.t = "RW2", n.t = 200, scale.model = FALSE) par(mfrow = c(1,2)) plot(1:dim(out)[2], out, col = 1, type = "l", xlab = "Time", ylab = "Random effects") # verify the second order difference is normally distributed first_diff <- diff(as.numeric(out[1,])) second_diff <- diff(first_diff) qqnorm(second_diff) abline(c(0,1)) ## Spatial-temporal random effects out <- rst(n=1, type = "st", type.t = "RW2", Amat = DemoMap$Amat, n.t = 50) dimnames(out) par(mfrow = c(1,1)) plot(1:dim(out)[3], out[1,1,], col = 1, type = "l", ylim = range(out), xlab = "Time", ylab = "Random effects") for(i in 2:4) lines(1:dim(out)[3], out[1,i,], col = i) legend("bottomright", colnames(DemoMap$Amat), col = c(1:4), lty = rep(1,4)) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.