[R-sig-ME] mysterious syntax error in call to nlmer
Ben Bolker
bbolker at gmail.com
Sun Mar 27 23:36:32 CEST 2011
Jacob Wegelin <jacobwegelin at ...> writes:
Thank you for your very clear and reproducible example!
> I'm trying to fit my first nlmer model. I follow Doug Bates's code
> at
> http://lme4.r-forge.r-project.org/slides/2009-07-01-Lausanne/8NLMMD.pdf
> and use lme4_0.999375-39. A call to nls works fine, but when I turn
> it into a call to nlmer and add the random component, errors are
> returned. Providing start as a named vector returns one kind of
> error; providing start as a list returns a different error. My data
> are simulated but they contain noise (per the warning on the nlmer
> help page). You can see a plot of the data at
> http://jacobwegelin.net/r-help/2011_03_27_1311_simulated.pdf. The
> error messages are as follows. Code to reproduce the entire problem
> is below under Details.
set.seed(1)
N<-100; nIDs<-20
X<-seq(1,19, length=N)
Y<-sin(X) + rnorm(N)/10
IDs<-1:nIDs
## I have a silly tangential question. Why put the commas
## on the next line? I haven't seen this convention much before,
## but it has started appearing on this list. In general it
## makes much more sense to me (and occasionally makes a
## difference) to put connectors (commas, arithmetic operators etc.)
## at the end of the previous line rather than the beginning of
## the next ...
trueRandomEffects<-data.frame(
xintercept=rnorm(nIDs)/3
, amplitude=1+ rnorm(nIDs)/9
, frequency =1+rnorm(nIDs)/10
, xshift =runif(nIDs) * pi
)
rownames( trueRandomEffects)<- IDs
mdoot<- DAT<-data.frame(ID=IDs[1], X=X, Y=Y)
for(thisperson in IDs[-1]) {
nextOne<-DAT
nextOne$Y<- (
trueRandomEffects[thisperson,"xintercept"]
+ trueRandomEffects[thisperson,"amplitude"]
* sin(
trueRandomEffects[thisperson,"frequency"] * X
- trueRandomEffects[thisperson,"xshift"]
)
+ rnorm(N)/10
)
nextOne$ID<- thisperson
mdoot<-rbind( mdoot, nextOne)
}
attr(mdoot, "ranef")<-trueRandomEffects
require(ggplot2)
print(
ggplot(data=mdoot, mapping=aes(x=X, y=Y, group=ID, color=ID))
+ geom_path())
whump<-nls(
Y ~ xintercept + amplitude * sin( frequency * X - xshift),
data=mdoot,
start=list(xintercept=1, amplitude=1.5, frequency=1, xshift=pi/2)
)
print(whump)
print(summary(whump))
mdoot$nlsfit<-predict(whump)
print(head(mdoot))
print(
ggplot(data=mdoot, mapping=aes(x=X, y=Y, group=ID, color=ID)) +
geom_path() +
geom_point(mapping=aes(x=X, y=nlsfit))
)
require(lme4)
## at the moment nlmer assumes that the objective function has
## a gradient attribute associated. *If* you have a simple
## objective function (i.e. one that R's built-in deriv()
## function can handle), this is simple to address
parnames <- c("xintercept","amplitude","frequency","xshift")
dd <- deriv(expression(xintercept + amplitude *
sin( frequency * X - xshift)),
namevec=parnames, function.arg=c(parnames,"X"))
roop <-nlmer(
Y ~ dd(xintercept,amplitude,frequency,xshift,X)~xshift|ID,
data=mdoot,
start=list(fixef=c(xintercept=1, amplitude=1.5,
frequency=1, xshift=pi/2))
)
summary(roop)
detach("package:lme4")
## however, you can also do this (possibly better?)
## with nlme:
library(nlme)
roop2 <-nlme( Y ~ dd(xintercept,amplitude,frequency,xshift,X),
data=mdoot,
fixed = xintercept + amplitude + frequency + xshift ~ 1,
random = xshift~1|ID,
start=c(xintercept=1, amplitude=1.5,
frequency=1, xshift=pi/2)
)
summary(roop2)
## it's a little bit alarming that nlme finds a non-zero xshift
## variance while nlmer doesn't ...
More information about the R-sig-mixed-models
mailing list