Many thanks for your answer !
It works also for me and I did not have any warnings messages for nls()
fittings.

Thanks !

Amélie

2012/3/6 Ben Bolker <bbolker@gmail.com>

> amelie pinet <amelie.pinet@...> 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.
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



-- 
Dr. Pinet Amélie
INRA, Plantes et Systèmes de culture Horticoles
Domaine St Paul, Site Agroparc, 84941 Avignon Cedex 9, France

	[[alternative HTML version deleted]]


