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.