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]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
