[R-sig-Geo] Linear mixed model with correlated random effects in R

[ONKELINX, Thierry]

Dear Timo,

Have a look at the INLA package (http://www.r-inla.org/)

The correlation structures in nlme work only on the residuals, not on the random effects.

Best regards,

Thierry

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx at inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey


-----Oorspronkelijk bericht-----
Van: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org] Namens Timo Schmid
Verzonden: woensdag 1 augustus 2012 11:13
Aan: r-sig-geo at r-project.org
Onderwerp: [R-sig-Geo] Linear mixed model with correlated random effects in R


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
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.