Joachim Audenaert <Joachim.Audenaert <at> pcsierteelt.be> writes:
> > Hi all, > > I would like to predict some values for an nls regression function > (functional response model Rogers type II). This is an asymptotic function > of which I would like to predict the asymptotic value > I estimated the paramters with nls, but can't seem to get predictions for > values of m choice...... > This is my script: > > RogersII_N <- > nls(FR~N0-lambertW(attackR3_N*Th3_N*N0* > exp(-attackR3_N*(24-Th3_N*N0)))/(attackR3_N*Th3_N), > start=list(attackR3_N=0.04,Th3_N=1.46),control=list(maxiter=10000)) > lista <- c(1,2,100,1000) > predict(RogersII_N,newdata=lista) > > I created a list (lista) with some values of which I would like the > predict function to give me function values > > What am I doing wrong? > If you want to use the predict() function, it is best (necessary?) to use the form of nls() where you pass data explicitly within a data frame, as the data= argument, instead of using variables in the global workspace. Also: for this model, the asymptotic value is (time/(handling time)) = 24/Th3_N ... so you don't really need predict() (although it is still useful) library(emdbook) rogers.pred <- function(N0,a,h,T) { N0 - lambertW(a*h*N0*exp(-a*(T-h*N0)))/(a*h) } params0 <- list(attackR3_N=0.05,Th3_N=1.5) params1 <- list(attackR3_N=0.04,Th3_N=1.46) dat <- expand.grid(N0=1:20,rep=1:10) dat$FR0 <- with(c(dat,params0), rogers.pred(N0,a=attackR3_N,h=Th3_N,T=24)) par(las=1,bty="l") with(dat,plot(N0,FR0,ylim=c(0,12))) set.seed(101) dat$FR <- rnorm(nrow(dat),dat$FR0,sd=1) with(dat,points(N0,FR,col=rgb(1,0,0,alpha=0.5))) RogersII_N <- nls(FR~rogers.pred(N0,attackR3_N,Th3_N,T=24), start=params1,data=dat,control=list(maxiter=10000)) lista <- c(1,2,100,1000,10000) predict(RogersII_N,newdata=data.frame(N0=lista)) N0vec <- 1:20 lines(N0vec,predict(RogersII_N,newdata=data.frame(N0=N0vec)), col=4) with(as.list(coef(RogersII_N)),24/Th3_N) ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.