[R] Complex nonlinear model

Douglas Bates bates at stat.wisc.edu
Sat Dec 8 15:06:48 CET 2001


I have been away from email for a week and may be coming late to a
discussion here.  The short answer is that you must have names on
your list or vector of starting estimates so that the arguments and
covariates in Y.model can be distinguished.

Try using
 simparj.st <- c(gN = gN, MnmN = MnmN, OptN = OptN, DIs = DIs,
      beta = beta, eta1 = eta1, eta2 = eta2)

Murray Jorgensen <maj at waikato.ac.nz> writes:

> I am running 1.3.1 on a Windows (NT 4.0) machine.
> 
> I am trying to fit a nonlinear model intended to predict crop yield from
> nutrient information.
> 
> I am troubled by an error message:
> 
> > library(nls)
> > simparj.fm <- nls(Y ~ Y.model(gN, MnmN, OptN, DIs,  beta, eta1, eta2,
> +       Popn, Dmax, AWC, SumEp, PotYield3, Nsupply),
> +       start = simparj.st, trace =T)
> Error in numericDeriv(form[[3]], names(ind), env) : 
>         theta should be of type character
> 
> I suspect that I am doing something stupid. My code in full is
> 
> # parameter assignments
> Nmin      <- 0.8852
> Nopt1     <- 16.78
> gN <- 0.5511
> E.nfert1 <- 0.3271
> E.nfert2 <- 0.6132
> Beta <- 0.8902
> Dls <- 0.5378
> eta1 <- 0.3791
> eta2 <- 0.6332
> PopStd <- 90468
> beta <- Beta
> DIs <- Dls
> MnmN <- Nmin
> OptN <- Nopt1
> 
> # simulate experimental data for predictors
> nsim <- 70
> Popn <- rnorm(nsim,PopStd,0.1*PopStd)
> Dmax <- rnorm(nsim,140.68,47.45)
> AWC <- rnorm(nsim,186.86,47.41)
> SumEp <- rnorm(nsim,318.54,32.53)
> PotYield3 <- rnorm(nsim,0.16180,0.01167)
> Nsoil <- rnorm(nsim,94.07,34.06)
> Bdfield <- rnorm(nsim,1.0590,0.1420)
> Bdlab <- rnorm(nsim,0.7876,0.1169)
> 
> Nfert.broad <- runif(nsim,95.3,576.5)
> Nfert.band <- runif(nsim,122,250)
> broad <- rbinom(nsim,1,0.2)
> Nfert.broad <- Nfert.broad*broad
> Nfert.band <- Nfert.band*(1 - broad)
> Nsupply<-Nsoil*Bdfield/Bdlab + Nfert.broad*E.nfert1 + Nfert.band*E.nfert2
> 
> # define model function
> 
> Y.model <- function(gN, MnmN, OptN, DIs,  beta, eta1, eta2,
>       Popn, Dmax, AWC, SumEp, PotYield3, Nsupply)
>       {
>       Ymax<- 1-ifelse(Popn<=PopStd, eta1, eta2)*log(Popn/PopStd)
>       Ymax <- Ymax*PotYield3*Popn/1000
>       Ymax <- Ymax*ifelse(Dmax<=DIs*AWC, 1, 1 - beta*(Dmax -DIs*AWC)/SumEp)
>       Nstar <- (Nsupply- MnmN*Ymax) / (OptN*Ymax - MnmN*Ymax)
>       Nstar<-pmax ( 0,Nstar)
>       Ystar<-ifelse(Nstar<1, (1 + gN*(1 - Nstar))* Nstar^(1+gN ), 1)
>       Ystar<-pmax ( 0, Ystar)
>       Y.model <- Ystar*Ymax
>       }
> 
> # generate response variable from model
> 
> Y <- Y.model(gN, MnmN, OptN, DIs,  beta, eta1, eta2,
>       Popn, Dmax, AWC, SumEp, PotYield3, Nsupply)
> Y <- Y + rnorm(nsim,0,1)
> 
> # attempt to fit starting from true parameters
> library(nls)
> simparj.st <- c(gN, MnmN, OptN, DIs,  beta, eta1, eta2)
> simparj.fm <- nls(Y ~ Y.model(gN, MnmN, OptN, DIs,  beta, eta1, eta2,
>       Popn, Dmax, AWC, SumEp, PotYield3, Nsupply),
>       start = simparj.st, trace =T)
> 
> 
> 
> 
> Dr Murray Jorgensen      http://www.stats.waikato.ac.nz/Staff/maj.html 
> Department of Statistics, University of Waikato, Hamilton, New Zealand 
> Email: maj at waikato.ac.nz                            Fax +64-7 838 4155
> Phone +64-7 838 4773 home phone +64-7 856 6705  Mobile +64-21 139 5862
> 
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
> 

-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list