# Nonlinear Fitting Function Estimates the parameters of a given parametric model using the `optim` function in package `stats`. ```r fitIPEC( expr, x, y, ini.val, weights = NULL, control = list(), fig.opt = TRUE, xlim = NULL, ylim = NULL, xlab = NULL, ylab = NULL ) ``` ## Arguments - `expr`: A given parametric model - `x`: A vector or matrix of observations of independent variable(s) - `y`: A vector of observations of response variable - `ini.val`: A vector or list of initial values of model parameters - `weights`: An optional vector of weights to be used in the fitting process. `weights` should be NULL or a numeric vector. If non-NULL, weighted least squares is used with weights `weights`; otherwise ordinary least squares is used. - `control`: A list of control parameters for using the `optim` function in package `stats` - `fig.opt`: An option to determine whether to draw the fitted curve - `xlim`: The shown range of the **x**-axis - `ylim`: The shown range of the **y**-axis - `xlab`: The label of the **x**-axis - `ylab`: The label of the **y**-axis ## Details The Nelder-Mead algorithm is the default in the `optim` function in package `stats`. The user can accurately estimate the model parameters by setting smaller relative convergence tolerance and larger maximum number of iterations in the input argument of `control`, e.g. `control=list(trace=FALSE, reltol=1e-20, maxit=50000)`, at the expense of the running speed. `ini.val` can be a vector or a list that has saved initial values for model parameters, e.g. `y = beta0 + beta1 * x + beta2 * x^2`, `ini.val = list(beta0=seq(5, 15, len=2), beta1=seq(0.1, 1, len=9),beta2=seq(0.01, 0.05, len=5))`, which is similar to the usage of the input argument of `start` of `nls` in package `stats`. In the `weights` argument option, the default is `weights = NULL`. In that case, ordinary least squares is used. The residual sum of squares (RSS) between the observed and predicted $y$ values is minimized to estimate a model's parameters, i.e., $$ \mbox{RSS} = \sum_{i=1}^{n}\left(y_i-\hat{y}_i\right)^{2} $$ where $y_i$ and $\hat{y}_i$ represent the observed and predicted $y$ values, respectively; and $n$ represents the sample size. If `weights` is a numeric vector, the weighted residual sum of squares is minimized, i.e., $$ \mbox{RSS} = \sum_{i=1}^{n}w_i\left(y_i-\hat{y}_i\right)^{2} $$ where $w_i$ is the $i$ elements of `weights`. ## Returns - **expr**: The formula used - **par**: The vector of estimates of parameters - **RSS**: The residual sum of squares or the weighted residual sum of squares - **R.sq**: The coefficient of determination or the weighted coefficient of determination - **n**: The number of data points, namely the sample size ## References Nelder, J.A. and Mead, R. (1965) A simplex method for function minimization. **Comput. J.** 7, 308$-$313. tools:::Rd_expr_doi("10.1093/comjnl/7.4.308") ## Author(s) Peijian Shi pjshi@njfu.edu.cn , Peter M. Ridland p.ridland@unimelb.edu.au , David A. Ratkowsky d.ratkowsky@utas.edu.au , Yang Li yangli@fau.edu . ## Note This function can be applicable to a nonlinear parametric model with a single independent variable or with multiple independent variables. `R.sq` is only used to help users intuitively judge whether the fitted curve seriously deviates from the actual observations. However, it should NOT be used to decide which of several competing models is the most appropriate (Pages 44$-$45 in Ratkowsky 1990). `RSS` and curvatures are among the suitable candidates to answer such a question. ## See Also `bootIPEC`, `optim` in package `stats` ## Examples ```r #### Example 1 ################################################################################### graphics.off() # The velocity of the reaction (counts/min^2) under different substrate concentrations # in parts per million (ppm) (Page 269 of Bates and Watts 1988) x1 <- c(0.02, 0.02, 0.06, 0.06, 0.11, 0.11, 0.22, 0.22, 0.56, 0.56, 1.10, 1.10) y1 <- c(76, 47, 97, 107, 123, 139, 159, 152, 191, 201, 207, 200) # Define the Michaelis-Menten model MM <- function(theta, x){ theta[1]*x / ( theta[2] + x ) } res0 <- fitIPEC(MM, x=x1, y=y1, ini.val=c(200, 0.05), xlim=c(0, 1.5), ylim=c(0, 250), fig.opt=TRUE) par1 <- res0$par par1 res0 # The input names of parameters will not affect the fitted results. # We can use other names to replace theta1 and theta2. iv.list1 <- list( theta1=seq(100, 300, by=50), theta2=seq(10, 100, by=10) ) result0 <- fitIPEC( MM, x=x1, y=y1, ini.val=iv.list1, xlim=c(0, 1.5), ylim=c(0, 250), fig.opt=FALSE, control=list(trace=FALSE, reltol=1e-20, maxit=50000) ) param1 <- result0$par param1 ################################################################################################## #### Example 2 ################################################################################### graphics.off() # Development data of female pupae of cotton bollworm (Wu et al. 2009) # References: # Ratkowsky, D.A. and Reddy, G.V.P. (2017) Empirical model with excellent statistical # properties for describing temperature-dependent developmental rates of insects # and mites. Ann. Entomol. Soc. Am. 110, 302-309. # Wu, K., Gong, P. and Ruan, Y. (2009) Estimating developmental rates of # Helicoverpa armigera (Lepidoptera: Noctuidae) pupae at constant and # alternating temperature by nonlinear models. Acta Entomol. Sin. 52, 640-650. # 'x2' is the vector of temperature (in degrees Celsius) # 'D2' is the vector of developmental duration (in d) # 'y2' is the vector of the square root of developmental rate (in 1/d) x2 <- seq(15, 37, by=1) D2 <- c(41.24,37.16,32.47,26.22,22.71,19.01,16.79,15.63,14.27,12.48, 11.3,10.56,9.69,9.14,8.24,8.02,7.43,7.27,7.35,7.49,7.63,7.9,10.03) y2 <- 1/D2 y2 <- sqrt( y2 ) ini.val1 <- c(0.14, 30, 10, 40) # Define the square root function of the Lobry-Rosso-Flandrois (LRF) model sqrt.LRF <- function(P, x){ ropt <- P[1] Topt <- P[2] Tmin <- P[3] Tmax <- P[4] fun0 <- function(z){ z[z < Tmin] <- Tmin z[z > Tmax] <- Tmax return(z) } x <- fun0(x) if (Tmin >= Tmax | ropt <= 0 | Topt <= Tmin | Topt >= Tmax) temp <- Inf if (Tmax > Tmin & ropt > 0 & Topt > Tmin & Topt < Tmax){ temp <- sqrt( ropt*(x-Tmax)*(x-Tmin)^2/((Topt-Tmin)*((Topt-Tmin )*(x-Topt)-(Topt-Tmax)*(Topt+Tmin-2*x))) ) } return( temp ) } myfun <- sqrt.LRF xlab1 <- expression( paste("Temperature (", degree, "C)", sep="" ) ) ylab1 <- expression( paste("Developmental rate"^{1/2}, " (", d^{"-1"}, ")", sep="") ) resu0 <- fitIPEC( myfun, x=x2, y=y2, ini.val=ini.val1, xlim=NULL, ylim=NULL, xlab=xlab1, ylab=ylab1, fig.opt=TRUE, control=list(trace=FALSE, reltol=1e-20, maxit=50000) ) par2 <- resu0$par par2 resu0 ################################################################################################## #### Example 3 ################################################################################### graphics.off() # Height growth data of four species of bamboo (Gramineae: Bambusoideae) # Reference(s): # Shi, P., Fan, M., Ratkowsky, D.A., Huang, J., Wu, H., Chen, L., # Fang, S. and Zhang, C. (2017) Comparison of two ontogenetic # growth equations for animals and plants. Ecol. Model. 349, 1-10. data(shoots) # Choose a species # 1: Phyllostachys iridescens; 2: Phyllostachys mannii; # 3: Pleioblastus maculatus; 4: Sinobambusa tootsik # 'x3' is the vector of the investigation times from a specific starting time of growth # 'y3' is the vector of the aboveground height values of bamboo shoots at 'x3' ind <- 4 x3 <- shoots$x[shoots$Code == ind] y3 <- shoots$y[shoots$Code == ind] # Define the beta sigmoid model (bsm) bsm <- function(P, x){ P <- cbind(P) if(length(P) !=4 ) {stop(" The number of parameters should be 4!")} ropt <- P[1] topt <- P[2] tmin <- P[3] tmax <- P[4] tailor.fun <- function(x){ x[x < tmin] <- tmin x[x > tmax] <- tmax return(x) } x <- tailor.fun(x) ropt*(x-tmin)*(x-2*tmax+topt)/(topt+tmin-2*tmax)*( (x-tmin)/(topt-tmin))^((topt-tmin)/(tmax-topt)) } ini.val2 <- c(40, 30, 5, 50) xlab2 <- "Time (d)" ylab2 <- "Height (cm)" re0 <- fitIPEC( bsm, x=x3, y=y3, ini.val=ini.val2, xlim=NULL, ylim=NULL, xlab=xlab2, ylab=ylab2, fig.opt=TRUE, control=list(trace=FALSE, reltol=1e-20, maxit=50000) ) par3 <- re0$par par3 ################################################################################################## #### Example 4 ################################################################################### # Data on biochemical oxygen demand (BOD; Marske 1967) # References: # Pages 56, 255 and 271 in Bates and Watts (1988) # Carr, N.L. (1960) Kinetics of catalytic isomerization of n-pentane. Ind. Eng. Chem. # 52, 391-396. data(isom) Y <- isom[,1] X <- isom[,2:4] # There are three independent variables saved in matrix 'X' and one response variable (Y) # The first column of 'X' is the vector of partial pressure of hydrogen # The second column of 'X' is the vector of partial pressure of n-pentane # The third column of 'X' is the vector of partial pressure of isopentane # Y is the vector of experimental reaction rate (in 1/hr) isom.fun <- function(theta, x){ x1 <- x[,1] x2 <- x[,2] x3 <- x[,3] theta1 <- theta[1] theta2 <- theta[2] theta3 <- theta[3] theta4 <- theta[4] theta1*theta3*(x2-x3/1.632) / ( 1 + theta2*x1 + theta3*x2 + theta4*x3 ) } ini.val8 <- c(35, 0.1, 0.05, 0.2) cons1 <- fitIPEC( isom.fun, x=X, y=Y, ini.val=ini.val8, control=list( trace=FALSE, reltol=1e-20, maxit=50000) ) par8 <- cons1$par ################################################################################################## ```

fitIPEC function