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]]

