Skip to content

error in nlme()

5 messages · amelie pinet, Ben Bolker

Original
Ben Bolker
# 
amelie pinet <amelie.pinet at ...> writes:
It's hard to answer without a reproducible example:
http://tinyurl.com/reproducible-000 ... the problem *could* be inside
your dbleseg1 function.  Does it work outside of nls/nlme?  Are there
any simplified versions that you have gotten to work?
It's not necessarily *the* problem, but it would generally be advised
to update to a more recent version of nlme -- especially if we have
trouble reproducing your problem with more recent versions.

  You can use traceback() to try to see where the error came from,
although it's probably pretty deep within a nested set of function
calls ...
3 days later
Ben Bolker
# 
amelie pinet <amelie.pinet at ...> writes:
[snip]
The main problem with your code is that you used DATA$Traitement
within the functions; when na.omit() scrubbed NAs from the input
data, it led to a mismatch between the data length and the length
of the 'Traitement' variable.  In general it's often simpler/more
robust to scrub NAs explicitly beforehand (using na.omit()) unless
you have particular reasons to need to retain those data throughout
the modeling process.  You can also (as I have done below) add
'Traitement' as a variable in your functions; in general it is 
bad practice to refer directly to global variables in your
nls functions, for exactly this reason ...

library(gdata)
## http://www.fichier-xls.fr/2012/03/06/7gvwshi/Phyllocrone_AP.xls
DATA <- read.xls("Phyllocrone_AP.xls",na.strings=c("NA","#VALUE!"),
comment.char="")

## for some reason I couldn't get #VALUE! read in as NA ...
DATA <- transform(DATA,Nb_phyto=as.numeric(as.character(Nb_phyto)))

# broken line with no effet of traitement
dbleseg0 <- function (JourJulien,a0,b0,a1,T) (a0*JourJulien+b0) +
  a1*(JourJulien >T)*(T-JourJulien)

# broken line with effet of traitement on a0 parameter
dbleseg1 <- function (JourJulien,a0,b0,a1,T,a0a,Traitement) {
    (((a0+a0a*(Traitement!="Control"))*JourJulien)+b0) +
a1*(JourJulien >T)*(T-JourJulien)
}

# broken line with effet of traitement on b0 parameter
dbleseg2 <- function (JourJulien,a0,b0,a1,T,b0a,Traitement) {
    ((a0*JourJulien)+(b0+b0a*(Traitement!="Control"))) +
a1*(JourJulien >T)*(T-JourJulien)
                }

# broken line with effet of traitement on a0 and b0 parameters
dbleseg3 <- function (JourJulien,a0,b0,a1,T,a0a,b0a,Traitement) {

  (((a0+a0a*(Traitement!="Control"))*JourJulien)+
    (b0+b0a*(Traitement!="Control")))+ a1*(JourJulien >T)*(T-JourJulien)
}

a0init <- 0.49
b0init <- -65
a1init <- 0.30
Tinit <- 239
a0ainit <- -0.40
b0ainit <- 6
a1ainit <- -0.1

library(nlme)

start1 <- c(a0=a0init,b0=b0init,a1=a1init,T=Tinit)
with(c(DATA,as.list(start1)),dbleseg0(JourJulien,a0,b0,a1,T))
Mod1 <- nls(Nb_phyto ~dbleseg0(JourJulien,a0,b0,a1,T),
data =DATA,start=start1,na.action=na.omit)

Mod2 <- nls(Nb_phyto ~dbleseg1(JourJulien,a0,b0,a1,T,a0a,Traitement),
data =DATA,start=start1,
na.action=na.omit)

Mod3 <- nls(Nb_phyto ~dbleseg2(JourJulien,a0,b0,a1,T,b0a,Traitement),
data =DATA,start=c(a0=a0init,b0=b0init,a1=a1init,T=Tinit,b0a=b0ainit),
na.action=na.omit)

Mod4 <- nls(Nb_phyto ~dbleseg3(JourJulien,a0,b0,a1,T,a0a,b0a,Traitement),
data = DATA,start1,na.action=na.omit)

*Mod2, Mod3, Mod4 ran but I obtained this warnings:
Il y a eu 50 avis ou plus (utilisez warnings() pour voir les 50 premiers)

  Did you look at the warnings????

  I had a problem with Mod4 ("singular gradient matrix
at initial parameter estimates") -- you should look at your
initial values and see whether they give rise to a sensible
estimate or not ...

Then I tried to use nlme()* (with the same start values) :

dataAna <- groupedData(Nb_phyto ~1|Plante,data=DATA)

Mod5 <- nlme(Nb_phyto ~dbleseg0(JourJulien,a0,b0,a1,T),data =
dataAna,fixed=a0+b0+a1+T~1,random=pdDiag(b0~1),
start=c(a0=a0init,b0=b0init,a1=a1init,T=Tinit),na.action=na.omit)

Mod6 <- nlme(Nb_phyto ~dbleseg1(JourJulien,a0,b0,a1,T,a0a,Traitement),data =
dataAna,fixed=a0+b0+a1+T+a0a~1,random=pdDiag(a0~1),
start=c(a0=a0init,b0=b0init,a1=a1init,T=Tinit,a0a=a0ainit),na.action=na.omit)

  these worked for me.