[R] creating a new corClass for lme()

[Viechtbauer Wolfgang (STAT)]

No idea if this helps, but if the optimizer (nlminb) is giving you problems, you could try switching to another optimizer (nlm) with:

control=list(opt="nlm")

Best,

--
Wolfgang Viechtbauer                        http://www.wvbauer.com/
Department of Methodology and Statistics    Tel: +31 (43) 388-2277
School for Public Health and Primary Care   Office Location:
Maastricht University, P.O. Box 616         Room B2.01 (second floor)
6200 MD Maastricht, The Netherlands         Debyeplein 1 (Randwyck)


----Original Message----
From: r-help-bounces at r-project.org
[mailto:r-help-bounces at r-project.org] On Behalf Of Michael Steven
Rooney Sent: Wednesday, April 14, 2010 07:07 To: r-help at r-project.org
Subject: [R] creating a new corClass for lme()

> 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 at 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.