IPEC | R Documentation |
Calculates the RMS intrinsic and parameter-effects curvatures of a nonlinear regression model. The curvatures are global measures of assessing whether a model/data set combination is close-to-linear or not. See Bates and Watts (1980) and Ratkowsky and Reddy (2017) for details.
The DESCRIPTION file:
This package was not yet installed at build time.
Index: This package was not yet installed at build time.
We are deeply thankful to Drs. Paul Gilbert and Jinlong Zhang for their invaluable help during creating this package.
Peijian Shi [aut, cre], Peter Ridland [aut], David A. Ratkowsky [aut], Yang Li [aut]
Maintainer: Peijian Shi <pjshi@njfu.edu.cn>
Bates, D.M and Watts, D.G. (1988) Nonlinear Regression Analysis and its Applications. Wiley, New York. doi: 10.1002/9780470316757
Ratkowsky, D.A. (1983) Nonlinear Regression Modeling: A Unified Practical Approach. Marcel Dekker, New York.
Ratkowsky, D.A. (1990) Handbook of Nonlinear Regression Models, Marcel Dekker, New York.
Ratkowsky, D.A. & 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. doi: 10.1093/aesa/saw098
hessian
in package numDeriv, jacobian
in package numDeriv, rms.curv
in package MASS
#### 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 res1 <- derivIPEC( MM, theta=par1, z=x1[1], method="Richardson", method.args=list(eps=1e-4, d=0.11, zero.tol=sqrt(.Machine$double.eps/7e-7), r=6, v=2) ) res1 # To calculate curvatures res2 <- curvIPEC( MM, theta=par1, x=x1, y=y1, alpha=0.05, method="Richardson", method.args=list(eps=1e-4, d=0.11, zero.tol=sqrt(.Machine$double.eps/7e-7), r=6, v=2) ) res2 # To calculate bias res3 <- biasIPEC(MM, theta=par1, x=x1, y=y1, tol= 1e-20) res3 set.seed(123) res4 <- bootIPEC( MM, x=x1, y=y1, ini.val=par1, control=list(reltol=1e-20, maxit=40000), nboot=2000, CI=0.95, fig.opt=TRUE ) res4 set.seed(NULL) # To calculate skewness res5 <- skewIPEC(MM, theta=par1, x=x1, y=y1, tol= 1e-20) res5 ################################################################################################# #### 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 resu1 <- derivIPEC( myfun, theta=par2, z=x2[1], method="Richardson", method.args=list(eps=1e-4, d=0.11, zero.tol=sqrt(.Machine$double.eps/7e-7), r=6, v=2) ) resu1 # To calculate curvatures resu2 <- curvIPEC( myfun, theta=par2, x=x2, y=y2, alpha=0.05, method="Richardson", method.args=list(eps=1e-4, d=0.11, zero.tol=sqrt(.Machine$double.eps/7e-7), r=6, v=2) ) resu2 # To calculate bias resu3 <- biasIPEC(myfun, theta=par2, x=x2, y=y2, tol= 1e-20) resu3 set.seed(123) resu4 <- bootIPEC( myfun, x=x2, y=y2, ini.val=ini.val1, nboot=2000, CI=0.95, fig.opt=TRUE ) resu4 set.seed(NULL) # To calculate skewness resu5 <- skewIPEC(myfun, theta=par2, x=x2, y=y2, tol= 1e-20) resu5 ################################################################################################# #### 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)) } # Define the simplified beta sigmoid model (simp.bsm) simp.bsm <- function(P, x, tmin=0){ P <- cbind(P) ropt <- P[1] topt <- P[2] tmax <- P[3] 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)) } # For the original beta sigmoid model 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 re1 <- derivIPEC( bsm, theta=par3, x3[15], method="Richardson", method.args=list(eps=1e-4, d=0.11, zero.tol= sqrt(.Machine$double.eps/7e-7), r=6, v=2) ) re1 re2 <- curvIPEC( bsm, theta=par3, x=x3, y=y3, alpha=0.05, method="Richardson", method.args=list(eps=1e-4, d=0.11, zero.tol= sqrt(.Machine$double.eps/7e-7), r=6, v=2) ) re2 re3 <- biasIPEC( bsm, theta=par3, x=x3, y=y3, tol= 1e-20 ) re3 re4 <- bootIPEC( bsm, x=x3, y=y3, ini.val=ini.val2, control=list(trace=FALSE, reltol=1e-20, maxit=50000), nboot=2000, CI=0.95, fig.opt=TRUE, fold=3.5 ) re4 re5 <- skewIPEC( bsm, theta=par3, x=x3, y=y3, tol= 1e-20 ) re5 # For the simplified beta sigmoid model # (in comparison with the original beta sigmoid model) ini.val7 <- c(40, 30, 50) RESU0 <- fitIPEC( simp.bsm, x=x3, y=y3, ini.val=ini.val7, xlim=NULL, ylim=NULL, xlab=xlab2, ylab=ylab2, fig.opt=TRUE, control=list(trace=FALSE, reltol=1e-20, maxit=50000) ) par7 <- RESU0$par par7 RESU1 <- derivIPEC( simp.bsm, theta=par7, x3[15], method="Richardson", method.args=list(eps=1e-4, d=0.11, zero.tol=sqrt(.Machine$double.eps/7e-7), r=6, v=2) ) RESU1 RESU2 <- curvIPEC( simp.bsm, theta=par7, x=x3, y=y3, alpha=0.05, method="Richardson", method.args=list(eps=1e-4, d=0.11, zero.tol=sqrt(.Machine$double.eps/7e-7), r=6, v=2) ) RESU2 RESU3 <- biasIPEC( simp.bsm, theta=par7, x=x3, y=y3, tol= 1e-20 ) RESU3 set.seed(123) RESU4 <- bootIPEC( simp.bsm, x=x3, y=y3, ini.val=ini.val7, control=list(trace=FALSE, reltol=1e-20, maxit=50000), nboot=2000, CI=0.95, fig.opt=TRUE, fold=3.5 ) RESU4 set.seed(NULL) RESU5 <- skewIPEC( simp.bsm, theta=par7, x=x3, y=y3, tol= 1e-20 ) RESU5 ################################################################################################## #### 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. graphics.off() 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 cons2 <- curvIPEC( isom.fun, theta=par8, x=X, y=Y, alpha=0.05, method="Richardson", method.args=list(eps=1e-4, d=0.11, zero.tol=sqrt(.Machine$double.eps/7e-7), r=6, v=2)) cons2 cons3 <- biasIPEC( isom.fun, theta=par8, x=X, y=Y, tol= 1e-20 ) cons3 set.seed(123) cons4 <- bootIPEC( isom.fun, x=X, y=Y, ini.val=ini.val8, control=list(trace=FALSE, reltol=1e-20, maxit=50000), nboot=2000, CI=0.95, fig.opt=TRUE, fold=10000 ) cons4 set.seed(NULL) cons5 <- skewIPEC( isom.fun, theta=par8, x=X, y=Y, tol= 1e-20 ) cons5 ##################################################################################################
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.