[R] SAR structure in linear mixed model in R

[Timo Schmid]

Hey,

I want to estimate a spatial linear mixed model y=X\beta+Zv+e with e\sim
N(0,4) and v\sim N(0,G). 
G is a covariance matrix of a simultaneously autoregressive model (SAR) and
is given by 
G=\sigma_v^2((I-pW)(I-pW^T))^{-1} where p is a spatial correlation parameter
and W is a matrix which describes the neighbourhood structure between the
random effects v. 

I have already seen some R-code when the SAR model is included in the error
term without area effects. This could be done with the errorsarlm-function
in the spdep package. I am looking for a "sarlme" function or something
similar.
Has anybody an idea?

I have also tried to incorporate the SAR model using the lme function. Here
is a short example:

The first code is for a linear mixed model without spatial correlation of
the random effects v. v\sim N(0,1).

/library (nlme)
library(mnormt)
library(spdep)
set.seed(10)
nclusters=12 #Number of area effects
areasize=20  #Number of units in each area
test.x <- runif(areasize*nclusters) #Covariates
test.e <- rnorm(areasize*nclusters) #Error term
ClusterID <- sort(sample(rep(1:nclusters,areasize)))
v <- numeric(nclusters)
for(i in 1:nclusters){v[((i-1)*areasize+1):(areasize+(i-1)*areasize)] <-
rep(rnorm(1, 0, 1), areasize)}
test.y <- 100+2*test.x+v+test.e 
test.fit <- lme(test.y~test.x,random = ~1 | ClusterID) #Fitting the model/

The second code should be for a linear mixed model with spatial correlation
in the random effects v. v\sim N(0,G). Therefore we have to generate the
matrix G.
/
library (nlme)
library(mnormt)
library(spdep)
set.seed(10)
nclusters=12 #Number of area effects
areasize=20  #Number of units in each area
test.x <- runif(areasize*nclusters) #Covariates
test.e <- rnorm(areasize*nclusters) #Error term
ClusterID <- sort(sample(rep(1:nclusters,areasize)))

# Generation of the neighbourhood matrix W
nb_structure <- cell2nb(3,4,type="rook") 
xyccc <- attr(nb_structure, "region.id")
xycc <- matrix(as.integer(unlist(strsplit(xyccc, ":"))), ncol=2, byrow=TRUE)
plot(nb_structure, xycc)
W_list<-nb2listw(nb_structure)
W <- nb2mat(nb_structure, glist=NULL, style="W") # Spatial Matrix W

spatial_p<-0.6 #spatial correlation parameter
G<-diag(1,nclusters)%*%solve((diag(1,nclusters)-diag(spatial_p,nclusters)%*%W)%*%(diag(1,nclusters)-diag(spatial_p,nclusters)%*%t(W)))
v1<-rmnorm(1, mean=rep(0,nclusters ), G)
v2<-v1/(sd(v1[1:nclusters]))*1
v<-rep(v2,each=areasize)
test.y <- 100+2*test.x+v+test.e/

Now the problem is the estimation. I have tried the lme and I wanted to
specify the correlation structure in the following way
/test.fit <- lme(test.y~test.x,random = ~1 | ClusterID, correlation=G)/
The point is that the correlation structure is not correct specfied with the
the matrix G. I know the matrix G must be a corStruct object, but I do not
know how to generate a SAR model as a corStruct object. In the corClasses
there is no SAR-model available. 
How can I incorporate the matrix G directly into the formula?

Thank you very much for your help.

Cheers,
Timo 






--
View this message in context: http://r.789695.n4.nabble.com/SAR-structure-in-linear-mixed-model-in-R-tp4637941.html
Sent from the R help mailing list archive at Nabble.com.