[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