[R-sig-ME] Problem using pdSymm matrix to fit spatial models
Virgilio Gomez Rubio
Virgilio.Gomez at uclm.es
Mon Dec 8 19:27:51 CET 2008
Dear all,
I am trying to use a Spatially Autoregressive specification for the
variance of the random effects to fit a model with spatial structure but
I am having some trouble with that. Basically, I want to define the
follwing variance matrix for the random effects:
V=sigma^2*((I-rho*W)*(I-rho*W))^{-1}
where I is the identity matrix, rho is the autocorrelation coefficient
and W is an adjacency matrix that measures dependence between areas.
The model is then
y_{ij}=mu_i+u_i+e_{ij}
(u_1, ..., u_n)~N(0, V)
e_{ij}~N(0, kappa^2)
The following is an example of what I am trying to do:
#Example with simulated data
library(nlme)
ssize<-rep(10, 42)
set.seed(1)
mu<-rnorm(42)*10
y<-rnorm(sum(ssize), mean=rep(mu, each=ssize), sd=.5)
Sdf<-data.frame(y=y, REGION=rep(1:42, each=ssize))
rho<-.2
n<-42
#W is an adjacency matrix
W<-matrix(0, ncol=42, nrow=42)
for(i in 1:41){W[i, i+1]<-.5;W[i+1,i]<-.5}
m<-pdSymm(solve( (diag(n)-rho*W)%*%(diag(n)-rho*t(W)) ),
nam=as.character(1:42),
form=~REGION-1)
model2<-lme(y~1, random=list(REGION=m), data=Sdf)
But I keep getting the following error:
Error in `Names<-.pdMat`(`*tmp*`, value = "REGION") :
Length of names should be 42
I think that I have used argument 'nam' in the right way when defining
the matrix, but it seems that I am missing something here.
I would appreciate any comments on this, in special if you think that
another approach would be better. The next step is to provide a way of
estimating rho, sigma^2 and kappa^2 (probably, by calling lme several
times) so that I can fit a wide range of spatial models.
Many thanks,
Virgilio
More information about the R-sig-mixed-models
mailing list