nlr {wnl}

Nonlinear Regression in R

Description

It performs nonlinear regression usually for pharmacokinetic and pharmacodynamic models.

Usage

nlr(Fx, Data, pNames, IE, LB, UB, Error="A", ObjFx=ObjDef, SecNames, SecForms, 
    Method="L-BFGS-B", Sx, conf.level=0.95, k, fix=0)

Arguments

Fx

Function for structural model. It should return a vector of the same length to observations.

Data

Data table which will be used in Fx. Fx should access this with e$DATA.

pNames

Parameter names in the order of Fx arguments

IE

Initial estimates of parameters

LB

Lower bound for optim function. The default value is 0.

UB

Upper bound for optim function. The default value is 1e+06.

Error

Error model. One of "A" for additive error, "POIS" for Poisson error, "P" for proportional error, "C" for combined error model, "S" for general error model. With Error="S", Sx should be provieded.

ObjFx

Objective function to be minimized. The default is maximum likelihood estimation function(-2 log likelihood).

SecNames

Names of secondary parameter estimates

SecForms

Formula to calculate the secondary parameter estimates

Method

"L-BFGS-B" is default. See optim for more detail.

Sx

Scale function. This is usually the inverse of weight. It should return the same length(nrow) of Y. When Error="S", Scale function should be provided as Sx.

conf.level

Confidence level for confidence interval

k

1/k likelihood interval(LI) will be provided. Currently recommended value is exp(qf(1 - alpha, 1, nRec-nPara)/2) + 1.

fix

indices of parameters to fix

Details

This uses scaled transformed parameters and environment e internally.

Value

Est

Point estimate(PE) with standard error(SE) and relative standard error(RSE)

LI

1/k likelihood interval, at which likelihood drops to 1/k of maximum likelihood. This reflects asymmetry better than confidence interval. This is estimated likelihood interval, not profile likelihood interval.

Skewness

Hougaard's skewness measure. This is printed only with additive error model. See also hSkew

Cov

Variance-covariance matrix of the objective function at the value of point estimates

run$m

Count of positive residuals

run$n

Count of negative residuals

run$run

Count of runs of residuals

run$p.value

P value of run test with excluding zero points

Objective Function Value

Minimum value of the objective function

-2LL

-2 times log likelihood

AIC

Akaike Information Criterion

AICc

Corrected Akaike Information Criterion

BIC

Schwarz Bayesian Information Criterion

Convergence

Convergence code from optim

Message

Message from optim.

Prediction

Fitted(predicted) values

Residuals

Residuals

Scale

Scales with Error="S". Variances for each points are scale vector multiplied by ScaleErrVar in Est.

Elapsed Time

Consumed time by minimization

Author(s)

Kyun-Seop Bae <[e-mail]>

Examples

  tData = Theoph
  colnames(tData) = c("ID", "BWT", "DOSE", "TIME", "DV")

  fPK = function(THETA) # Prediction function
  {
    DOSE = 320000 # in microgram
    TIME = e$DATA[, "TIME"] # use data in e$DATA

    K    = THETA[1]
    Ka   = THETA[2]
    V    = THETA[3]

    P  = DOSE/V*Ka/(Ka - K) * (exp(-K*TIME) - exp(-Ka*TIME))
    return(P)
  }

  IDs = unique(tData[,"ID"])
  nID = length(IDs)
  for (i in 1:nID) {
    Data = tData[tData$ID == IDs[i],]
    Res = nlr(fPK, Data, pNames=c("k", "ka", "V"), IE=c(0.1, 3, 500), 
              SecNames=c("CL", "Thalf", "MRT"), SecForms=c(~V*k, ~log(2)/k, ~1/k))
    print(paste("## ID =", i, "##"))
    print(Res)
  }

  # Another example from radioimmunoassay(RIA)
  d1 = data.frame(conc = c(200, 100, 50, 25, 12.5, 6.25, 3.125, 0),
                  DV = c(1.78, 1.5, 1.17, 0.74, 0.51, 0.31, 0.19, 0.04))

  PRED = function(TH) TH[1] + TH[2]*d1$conc^TH[4]/(TH[3]^TH[4] + d1$conc^TH[4])
  Scale = function(TH) 1/(PRED(TH) - (TH[1] + TH[2])/2)^2

  nlr(PRED, d1, pNames=c("R0", "Rmax", "RC50", "Hill"), IE=c(0.1, 3, 50, 1), 
      Error="S", Sx=Scale)