Hi,
I have been using the function lme() of the package nlme to model grouped
data that is auto-correlated in time and in space (the data was collected on
different days via a moving monitor). I am aware that I can use the
correlation classes corCAR1 and corExp (among other options) to model the
temporal and spatial components of the auto-correlation. However, as far as
I can tell, I can only model using one correlation class or the other. I
would like my model to account for both temporal and spatial autocorrelation
simultaneously.
The ?corClasses help page provides some advice on how to create your own
correlation class. I was able to create a new class that I called corSPT: I
add the spatial and temporal autocorrelation matrices to produce the corSPT
correlation matrix. The code for this is pasted below:
################## BEGIN ###############
# Create some data
N <- 100
x <- round(sin(rep(1:23/2,length.out=N)),digits=2)+1:N*2/N
y <- round(cos(rep(1:23/2,length.out=N)),digits=2)+1:N*2/N
g <- rep(1:5,each=N/5)
a <- round(runif(N,0,10))
t <- 1:N
r <- runif(N,0,5)
e <- 5*sin(4*x) +
5*cos(4*y) +
5*sin(t) +
2*g +
a +
r
e <- round(e)
df <- data.frame(x,y,g,a,t,r,e)
library(nlme)
# Define constructor function and methods
corSPT <- function(a,b) {
object <- list("time"=a,"space"=b)
class(object) <- c("corSPT","corStruct")
return(object)
}
coef.corSPT <- function(object,...) {
c(coef(object[["time"]],...)[1],coef(object[["space"]],...))
}
"coef<-.corSPT" <- function(object,...,value) {
coef(object[["time"]]) <- value[1]
coef(object[["space"]]) <- value[2:3]
class(object) <- c("corSPT","corStruct")
return(object)
}
Initialize.corSPT <- function(object,...) {
object <-
list("time"=Initialize(object[["time"]],...),"space"=Initialize(object[["space"]],...))
class(object) <- c("corSPT","corStruct")
return(object)
}
corMatrix.corSPT <- function(object,covariate = NULL, ...) {
a <-
corMatrix(object[["time"]],covariate=getCovariate(object[["time"]]),...)
b <-
corMatrix(object[["space"]],covariate=getCovariate(object[["space"]]),...)
lapply(seq(length(a)),function(n) pmax(-1,pmin(1,a[[n]]+b[[n]])))
}
formula.corSPT <- function(object,...) {
a <- as.character(formula(object[["time"]]))
b <- as.character(formula(object[["space"]]))
a <- strsplit(a[2],split="|",fixed=TRUE)[[1]]
b <- strsplit(b[2],split="|",fixed=TRUE)[[1]]
as.formula(paste("~",a[1],"+",b[1],"|",a[2]))
}
Dim.corSPT <- function(object,...) {
Dim(object[["time"]],...)
}
############ END ###############
When I run this model on:
mymodel <- lme(fixed = e ~ a,random= ~ 1 |
g,data=df,correlation=corSPT(corCAR1(.2,form = ~ t | g),corExp(c(1,.5),form=
~ x + y | g, nugget=TRUE)),control=list(msVerbose=TRUE))
I get sensible results. However if I change the way the temporal data is
modeled:
mymodel <- lme(fixed = e ~ a,random= ~ 1 |
g,data=df,correlation=corSPT(corExp(.2,form = ~ t | g),corExp(c(1,.5),form=
~ x + y | g, nugget=TRUE)),control=list(msVerbose=TRUE))
I get a C runtime error. I have been trying to track this down using
browser() and debug(), but all I have been able to tell is that the problem
arises during
.Call(R_port_nlminb,obj,grad,hess,rho,low,upp,d=rep(as.double(scale),length.out=length(par)),iv,v)
in the nlminb() function.
Please let me know if you know what might be going on. I am using R version
2.7.2 on Windows XP.
Thanks,
Mike
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
