nlr function

Nonlinear Regression in R

It performs nonlinear regression usually for pharmacokinetic and pharmacodynamic models.

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.

Returns

• 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 k@acr.kr

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)