confcurves: Wald Confidence Curves and the Likelihood Confidence Curves

Description

Calculates the Wald confidence curves and the likelihood confidence curves of model parameters.

Usage

confcurves( expr, x, y, ini.val, weights = NULL, control=list(), 
            fig.opt = TRUE, fold = 5, np = 20, alpha = seq(1, 0.001, by=-0.001), 
            show.CI = NULL, method = "Richardson", method.args = 
            list(eps = 1e-04, d = 0.11, zero.tol = sqrt(.Machine$double.eps/7e-07), 
            r = 6, v = 2, show.details = FALSE), side = NULL )

Value

partab: The estimates, standard errors and confidence intervals of model parameters; also see the parinfo function
parlist: A list for the estimate, Wald interval curves and likelihood interval curves of each model parameter.

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 confidence curves of each parameter
fold: The fold of $ \mbox{SE}\left(\hat{\theta_{i}}\right) $ for controlling the width of the confidence interval of $\hat{\theta_{i}}$ that represents the estimate of the $i$th parameter
np: The number of data points for forming the lower or upper bounds of a likelihood confidence interval of $ \hat{\theta_{i}} $, which controlls the step size (i.e., $\delta$) in the $y$ coordinates of the likelihood confidence curves
alpha: The significance level(s) for calculating the $x$ coordinate(s) of the $(1-\alpha)100\%$ Wald confidence curves, which equals to $ t_{\alpha/2}\left(n-p\right) $
show.CI: The $ t_{\alpha/2}\left(n-p\right) $ value(s) associated with the confidence level(s) of each parameter to be showed, i.e., c(0.80, 0.90, 0.95, 0.99)
method: It is the same as the input argument of method of the hessian function in package numDeriv
method.args: It is the same as the input argument of method.args of the hessian function in package numDeriv
side: It is the same as the input argument of side of the jacobian function in package numDeriv

Author

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.

Details

For the $(1-\alpha)100\%$ Wald confidence curves, the corresponding $x$ and $y$ coordinates are: $$ x = t_{\alpha/2}\left(n-p\right), $$ and $$ y = \hat{\theta_{i}} \pm t_{\alpha/2}\left(n-p\right)\,\mbox{SE}\left(\hat{\theta_{i}}\right), $$ where $n$ denotes the number of the observations, $p$ denotes the number of model parameters, and $\mbox{SE}\left(\hat{\theta_{i}}\right)$ denotes the standard error of the $i$th model parameter's estimate.

$\quad$ For the likelihood confidence curves (Cook and Weisberg, 1990), the corresponding $x$ and $y$ coordinates are: $$ x = \sqrt{\frac{\mbox{RSS}\left(\hat{\theta}^{\,\left(-i\right)}\right)-\mbox{RSS}\left(\hat{\theta}\right)}{\mbox{RSS}\left(\hat{\theta}\right)/(n-p)}}, $$ where $\mbox{RSS}\left(\hat{\theta}\right)$ represents the residual sum of squares for fitting the model with all model parameters; $\mbox{RSS}\left(\hat{\theta}^{\,\left(-i\right)}\right)$ represents the residual sum of squares for fitting the model with the $i$th model parameter being fixed to be $\hat{\theta_{i}} \pm k\,\delta$. Here, $k$ denotes the $k$th iteration, and $\delta$ denotes the step size, which equals $$ \delta = \frac{\hat{\theta_{i}} \pm \mbox{fold}\,\mbox{SE}\left(\hat{\theta_{i}}\right)}{\mbox{np}}. $$ $$ y = \hat{\theta_{i}} \pm k\,\delta. $$ Here, fold and np are defined by the user in the arguments.

$\quad$ For other arguments, please see the fitIPEC and parinfo functions for details.

References

Cook, R.D. and Weisberg, S. (1990) Confidence curves in nonlinear regression. J. Am. Statist. Assoc. 82, 221$-$230. tools:::Rd_expr_doi("10.1080/01621459.1990.10476233")

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")

Ratkowsky, D.A. (1990) Handbook of Nonlinear Regression Models, Marcel Dekker, New York.

Examples

Run this code

#### Example 1 ###################################################################################
# Weight of cut grass data (Pattinson 1981)
# References:
#   Clarke, G.P.Y. (1987) Approximate confidence limits for a parameter function in nonlinear 
#       regression. J. Am. Stat. Assoc. 82, 221-230.
#   Gebremariam, B. (2014) Is nonlinear regression throwing you a curve? 
#       New diagnostic and inference tools in the NLIN Procedure. Paper SAS384-2014.
#       http://support.sas.com/resources/papers/proceedings14/SAS384-2014.pdf
#   Pattinson, N.B. (1981) Dry Matter Intake: An Estimate of the Animal
#       Response to Herbage on Offer. unpublished M.Sc. thesis, University
#       of Natal, Pietermaritzburg, South Africa, Department of Grassland Science.

# 'x4' is the vector of weeks after commencement of grazing in a pasture
# 'y4' is the vector of weight of cut grass from 10 randomly sited quadrants

x4 <- 1:13
y4 <- c(3.183, 3.059, 2.871, 2.622, 2.541, 2.184, 
        2.110, 2.075, 2.018, 1.903, 1.770, 1.762, 1.550)

# Define the first case of Mitscherlich equation
MitA <- function(P, x){    
    P[3] + P[2]*exp(P[1]*x)
}

# Define the second case of Mitscherlich equation
MitB <- function(P2, x){
    if(P2[3] <= 0)
        temp <- mean(y4)
    if(P2[3] > 0)
       temp <- log(P2[3]) + exp(P2[2] + P2[1]*x)
    return( temp )
}

# Define the third case of Mitscherlich equation
MitC <- function(P3, x, x1=1, x2=13){
    theta1 <- P3[1]
    beta2  <- P3[2]
    beta3  <- P3[3]
    theta2 <- (beta3 - beta2)/(exp(theta1*x2)-exp(theta1*x1))
    theta3 <- beta2/(1-exp(theta1*(x1-x2))) - beta3/(exp(theta1*(x2-x1))-1)
    theta3 + theta2*exp(theta1*x)
}

ini.val3 <- c(-0.1, 2.5, 1)
RESU1    <- confcurves( MitA, x=x4, y=y4, ini.val=ini.val3, fig.opt = TRUE, 
                        fold=5, np=20, alpha=seq(1, 0.001, by=-0.001), 
                        show.CI=c(0.8, 0.9, 0.95, 0.99) )

ini.val4 <- c(-0.10, 0.90, 2.5)
RESU2    <- confcurves( MitB, x=x4, y=y4, ini.val=ini.val4, fig.opt = TRUE, 
                        fold=5, np=200, alpha=seq(1, 0.001, by=-0.001), 
                        show.CI=c(0.8, 0.9, 0.95, 0.99) )

ini.val6 <- c(-0.15, 2.5, 1)
RESU3    <- confcurves( MitC, x=x4, y=y4, ini.val=ini.val6, fig.opt = TRUE, 
                        fold=5, np=20, alpha=seq(1, 0.001, by=-0.001), 
                        show.CI=c(0.8, 0.9, 0.95, 0.99) )
##################################################################################################

graphics.off()

Run the code above in your browser using DataLab