
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 to model between-group correlation.

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.  		 	   		  
	[[alternative HTML version deleted]]

