[R-sig-ME] mysterious syntax error in call to nlmer

[Ben Bolker]

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