[R] - help with the predict function

Ben Bolker bbolker at gmail.com
Tue Jun 5 16:52:39 CEST 2012


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)



More information about the R-help mailing list