[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