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.
