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!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.