# How to implement a Gaussian Markov Random field (GMRF) model in R?

I would like to model a (conditional) GMRF using a linear mixed effects model without having grid Data but only a neighbourhood matrix $W$. My model is given by $$Y=X\beta+ \epsilon$$ and the error term $\epsilon$ can be decomposed additively as $\epsilon=s+e$, while $e$ is just an i.i.d. random error term, $s$ represents a spatial component that follows a GMRF process, i.e. $s_{i}|(s_{j}, j \in N_i)\sim \mathcal{N}(|N_i|^{-1}\sum_{j \in N_j}s_j,|N_i|^{-1}\sigma^2_{s} )$. Where $N_i$ represents the neighbourhood of observation $i$. I would like to estimate this approach in R. However, I mostly see packages that are constructed for spatial econometrics using models like the SAR model $$Y=X\beta+ \epsilon\text{, with }\epsilon = \rho W\epsilon+e$$ or the CAR model $$Y=X\beta+ \rho WY+ \epsilon,$$ where $W$ containes the spatial weights or represents a binary adjacency matrix. Because I want to estimate the spatial effects I have tried to implement the model in a linear mixed model framework with $$Y=X\beta+ \Omega^{1/2}\tilde{s} +\epsilon\text{, with } \tilde{s} \sim \mathcal{N}(0,\sigma^2_sI_N)\text{, and } \epsilon \sim \mathcal{N}(0,\sigma^2_\epsilon I_N)\text{ (1)}$$ where $\Omega$ is defined analogously to Kauermann, Haupt and Kaufmann (2012) as $$\Omega \equiv (I_N-W)^{-}\text{diag}(|N_i|^{-1})$$, where $(I_N-W)^{-}$ denotes the Moore-Penrose inverse of $I_N-W$. However, I am not able to fit model $(1)$ with standard linear mixed models procedures in R, as nlme, lm4 or varComp. The main difficulties arise as I do not have Grid Data, i.e. longitude and lattitude, but only a neighbourhood matrix $W$. Hence I can easily fit spatial econometric models as the SAR or the CAR model, but not a model with GMRF random effects. Any help is appreciated!