[R] making self-starting function for nls

Douglas Bates dmbates at gmail.com
Thu Sep 1 17:27:47 CEST 2005


On 9/1/05, Bill Shipley <bill.shipley at usherbrooke.ca> wrote:
> Hello.  Following pages 342-347 of Pinheiro & Bates, I am trying to
> write a self-starting nonlinear function (a non-rectagular hyperbola) to
> be used in nonlinear least squares regression (and eventually for a
> mixed model).  When I use the getInitial function for my self-starting
> function I get the following error message:

> > getInitial(photo~NRhyperbola(Irr,theta,Am,alpha,Rd),dat)
> 
> Error in tapply(y, x, mean, na.rm = TRUE) :
> 
>         arguments must have same length
> 
> Since I do not explicitly call tapply in my function that makes
> NRhyperbola a self-starting function (called NRhyperbolaInit, see
> below), I assume that the error is coming from within the mCall function
> but I can't figure out where or how.

My guess is that it is occuring in the call to sortedXyData but I
won't be able to tell for sure without test data.

One of the things that sortedXyData does is to average the y values
for replicated x values.  It seems that in the call the lengths of the
x and y arguments are different.

> 
> 
> 
> Would someone who has successfully done this be willing to look at my
> code and see where the problem arises?
> 
> 
> 
> > NRhyperbolaInit
> 
> function(mCall,LHS,data)
> 
> {
> 
> xy<-sortedXyData(mCall[["x"]],LHS,data)
> 
> if(nrow(xy)<3){
> 
>  stop("Too few unique irradiance values")
> 
> }
> 
> theta<-0.75
> 
> Rd<-min(xy[,"y"])
> 
> Am<-max(xy[,"y"]) + abs(Rd)
> 
> if(sum(xy[,"x"]<50)>3)alpha<-coef(lm(y~x,data=xy,subset=x<50))[2]
> 
> if(sum(xy[,"x"]<50)<=3)alpha<-0.07
> 
> value<-c(theta,Am,alpha,Rd)
> 
> names(value)<-mCall[c("theta","Am","alpha","Rd")]
> 
> value
> 
> }
> 
> 
> 
> Bill Shipley
> 
> Bill.Shipley at USherbrooke.ca
> 
>  <http://callisto.si.usherb.ca:8080/bshipley/>
> http://callisto.si.usherb.ca:8080/bshipley/
> 
> 
> 
> 
>         [[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>




More information about the R-help mailing list