Hi Uwe,
I downloaded R 2.11.0 and the latest version of nlme. I am still getting the
problem. I should also note that I got a message from Viechtbauer Wolfgang
saying that he had experienced a similar problem using nlme.
I have pasted the code below. (The data is generated in the first ~15
lines). I have bunch of debug() statements, which may or may not be useful.
The code is not as elegant as it could be since I have been constantly
editing it to try to track down the bug (I replaced the corMatrix C routine
with R code). FYI, the objective of the code is to simultaneously fit
temporal and spatial correlation (adding one corExp matrix plus another)
where the nuggets are constrained to add to 1.
I really hope you have a chance to look at this. Thanks so much!!
## BEGIN
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)
df <- df[-c(22,23,67),]
library(nlme)
corSPT <- function(a,b) {
object <- list("time"=a,"space"=b)
class(object) <- c("corSPT","corSpatial","corStruct")
attr(object, "formula") <- formula(object)
attr(object, "nugget") <- attr(object[["time"]], "nugget")
attr(object, "metric") <- attr(object[["time"]], "metric")
attr(object, "fixed") <- attr(object[["time"]], "fixed")
return(object)
}
coef.corSPT <- function(object,...) {
c("time"=coef(object[["time"]],...),"space"=coef(object[["space"]],...)[1])
}
"coef<-.corSPT" <- function (object, ..., value) {
value <- as.numeric(value)
if (length(value) != 3) {
stop("Cannot change the length of the parameter of a corSPT object")
}
object[["time"]][] <- value[1:2]
object[["space"]][] <- c(value[3],-value[2]) # nugget for space is
determined by the nugget for time
attr(object, "factor") <- NULL
attr(object, "factor") <- corFactor(object)
attr(object, "logDet") <- NULL
attr(object, "logDet") <- logDet(object)
object
}
Initialize.corSPT <- function (object, data, ...) {
if (!is.null(attr(object[["time"]], "minD"))) {
return(object)
}
object[["time"]] <-
getS3method("Initialize","corStruct")(object[["time"]], data)
object[["space"]] <-
getS3method("Initialize","corStruct")(object[["space"]], data)
object[["time"]][1] <- log(object[["time"]][1])
object[["space"]][1] <- log(object[["space"]][1])
object[["time"]][2] <- log(object[["time"]][2]/(1 -
object[["time"]][2]))
object[["space"]][2] <- -object[["time"]][2]
attr(object,"groups") <- attr(object[["time"]],"groups")
attr(object,"Dim") <- attr(object[["time"]],"Dim")
attr(object,"covariate") <- NULL
attr(object, "nugget") <- attr(object[["time"]], "nugget")
attr(object[["time"]],"minD") <- min(unlist(attr(object[["time"]],
"covariate")))
attr(object[["space"]],"minD") <- min(unlist(attr(object[["space"]],
"covariate")))
attr(object, "minD") <- NULL
###
corD <- Dim(object[["time"]])
attr(object[["time"]],"covariate") <- lapply(seq(corD$M),function(M) {
x <- attr(object[["time"]],"covariate")[[M]]
l <- corD$len[M]
g <- rep(seq(l),(l-1):0)
mat <- do.call(rbind,lapply(seq(l),function(G) {c(rep(0,G),x[g==G])}))
mat <- mat + t(mat)
})
attr(object[["space"]],"covariate") <- lapply(seq(corD$M),function(M) {
x <- attr(object[["space"]],"covariate")[[M]]
l <- corD$len[M]
g <- rep(seq(l),(l-1):0)
mat <- do.call(rbind,lapply(seq(l),function(G) {c(rep(0,G),x[g==G])}))
mat <- mat + t(mat)
})
###
attr(object, "factor") <- corFactor(object)
attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
object
}
corMatrix.corSPT <- function(object,covariate = NULL,corr=TRUE, ...) {
if (corr) {
temp <- coef(object[["time"]],unconstrained=FALSE)
cat(temp)
time.rng <- temp[1]
time.nug <- temp[2]
temp <- coef(object[["space"]],unconstrained=FALSE)
space.rng <- temp[1]
space.nug <- temp[2]
cor.time <- lapply(attr(object[["time"]],"covariate"),function(mat)
(1-time.nug)*exp(-mat/time.rng))
cor.space <- lapply(attr(object[["space"]],"covariate"),function(mat)
(1-space.nug)*exp(-mat/space.rng))
x <- lapply(seq(length(cor.time)),function(mat)
round(cor.time[[mat]]+cor.space[[mat]],4))
attr(x, "logDet") <- NULL
} else {
d <- Dim(object)
gg <- rep(seq(d$M),each=d$len^2)
cf <- corFactor(object)
lD <- attr(cf, "logDet")
x <- lapply(seq(d$M),function(mm) {mat <- cf[gg==mm];
matrix(mat,nrow=d$len[mm])})
attr(x,"logDet") <- lD
}
return(x)
}
corFactor.corSPT <- getS3method("corFactor","corStruct")
logDet.corSPT <- getS3method("logDet","corStruct")
recalc.corSPT <- getS3method("recalc","corStruct")
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"]],...)
}
debug(corSPT)
debug(getS3method("recalc","corSpatial"))
debug(getS3method("recalc","corStruct"))
debug(getS3method("Initialize","corSPT"))
debug(getS3method("coef<-","corSPT"))
debug(getS3method("coef","corSPT"))
debug(getS3method("corFactor","corSPT"))
debug(getS3method("logDet","corStruct"))
debug(getS3method("logLik","lmeStruct"))
debug(getS3method("recalc","modelStruct"))
debug(getS3method("recalc","corStruct"))
debug(getS3method("corMatrix","corSPT"))
debug(getS3method("coef<-","reStruct"))
debug(getS3method("coef<-","modelStruct"))
debug(getS3method("corMatrix","corSpatial"))
debug(getS3method("Initialize","corSpatial"))
debug(getS3method("Initialize","corStruct"))
mymodel <- lme(fixed = e ~ a,random= ~ 1 |
g,data=df,correlation=corSPT(corExp(c(1,.5),form = ~ t |
g,nugget=TRUE),corExp(c(1,.5),form= ~ x + y | g,
nugget=TRUE)),control=list(msVerbose=TRUE,opt="nlminb"))
## End
2010/4/21 Uwe Ligges <lig...@statistik.tu-dortmund.de>
>
>
> On 21.04.2010 18:19, Michael Steven Rooney wrote:
>
>> The R version is 2.9.2.
>>
>
> Please try with R-2.11.0 RC (release candidate) to be released tomorrow as
> well as a recent version of nlme.
>
> If its still fails, please send code (and data, if you do not generate them
> in the code) that reproduces the problem.
>
> Uwe Ligges
>
>
>
>
>
> nlme version is 3.1-93
>> The OS is Windows XP.
>> This is for my work computer but I get the same behavior at home using
>> Vista.
>> ---------- Forwarded message ----------
>> From: Sarah Goslee<sarah.gos...@gmail.com>
>> Date: Wed, Apr 21, 2010 at 11:50 AM
>> Subject: Re: [R] R crashing oddly
>> To: Michael Steven Rooney<michael.s.roo...@gmail.com>
>>
>>
>> You definitely need to provide the information requested in the posting
>> guide,
>> specifically OS, R version, package version (and update to the newest if
>> you
>> haven't yet).
>>
>> Sarah
>>
>> On Wed, Apr 21, 2010 at 11:41 AM, Michael Steven Rooney
>> <michael.s.roo...@gmail.com> wrote:
>>
>>> Hi,
>>>
>>> I am working with the package nlme, and I tried creating a new
>>> correlation
>>> class (which, according to the help pages, is possible if you write a few
>>> new method functions). Anyways, I think I am 99% of the way there, but I
>>> have a recurring problem with R crashing on seemingly innocuous
>>> statements
>>> (I have set debug() on nearly every function, so I can see where it is
>>> failing).
>>>
>>> --
>> Sarah Goslee
>> http://www.functionaldiversity.org
>>
>> [[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.
>>
>
[[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.
