I have a dataset that contains two sites (i.e., sites 1 and 2). The data in each site is arranged in a row-column layout (e.g., 3 rows by 4 columns). Below is an example for my dataset, where trt is the main factor of interest, and row, column, and site are blocking factors.
nsite <- 2
nrow <- 3
ncol <- 4
nperSite <- nrow*ncol
site <- factor(rep(1:nsite, each = nperSite))
d1 <- c(1, 2, 3, 4, 2, 3, 4, 1, 3, 4, 1, 2)
d2 <- c(1, 3, 4, 2, 3, 4, 2, 1, 4, 2, 1, 3)
trt <- factor(c(d1,d2))
row <- c(rep(1:nrow,each=ncol),rep(1:nrow,each=ncol))
col <- c(rep(1:ncol,nrow),rep(1:ncol,nrow))
y <- matrix(0,nrow=nsite*nperSite,1)
data <- data.frame(trt, row, col, site, y)
I would like to build a spatial model for my data, where each site has its own random error (epsilon1 for site 1 and epsilon2 for site 2). Within each site, I want to use the first-order auto-regressive model to model a spatial correlation for rows and columns (see Section 3.2.1 in https://idahoagstats.github.io/guide-to-field-trial-spatial-analysis/background.html). This means for site 1, we have a spatial correction: AR(r1) and AR(c1), and for site 2, we have a spatial correlation: AR(r2) and AR(c2). Does anyone know whether it is possible to fit such a spatial model in spaMM (https://gitlab.mbb.univ-montp2.fr/francois/spamm-ref), glmmTMB, or R-INLA?
Thank you in advance.