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 [email protected] 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: [email protected] [mailto:[email protected]] Namens Timo Schmid Verzonden: woensdag 1 augustus 2012 11:13 Aan: [email protected] 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 [email protected] 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. _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
