[R-sig-ME] error in nlme()

[Ben Bolker]

amelie pinet <amelie.pinet at ...> writes:

> 
> Hello,
> 
> *Here it is a copy of my data:*
> 
  [snip]

> 
> *I also post my data here (be careful, it's a .xls file)*
> http://www.fichier-xls.fr/2012/03/06/7gvwshi/*
> 
> The code used is as follow:*
> 

  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.