dummy-link

NLreg

Nonlinear regression in Julia

Readme

Nonlinear regression models in Julia

Build Status

In this Julia package nonlinear regression models are formulated as Julia types that inherit from NLregMod. A simple example is the predicted concentration in a 1 compartment model with a single bolus dose at time 0.

conc = exp(logV) * exp(-exp(logK)*t)

where logV and logK are the logarithms of the volume of distribution and K is the elimination rate constant and t is the time of the measurement.

The logsd1 type represents this model and the data to which it is to be fit. The fields of this type are t, the vector of times at which samples are drawn, y, the vector of measured concentrations, mu the mean responses at the current parameter values, resid the residuals at the current parameter values, and tgrad, the transpose of the gradient matrix. The external constructors for this model allow it to be specified from t and y or in a Formula/Data specification.

A nonlinear regression model must provide methods for pnames, the parameter names, updtmu, update the mean response, residuals and tgrad from new parameter values, and initpars, determine initial parameter estimates from the data.

julia> using DataFrames, NLreg

julia> const sd1 = readtable(Pkg.dir("NLreg","data","sd1.csv.gz"));

julia> nl = fit(BolusSD1(CONC ~ TIME, sd1))
Nonlinear least squares fit to 580 observations

     Estimate Std.Error t value Pr(>|t|)
V     1.14296 0.0495656 23.0595  < eps()
K    0.245688 0.0202414 12.1379  < eps()

Residual sum of squares at estimates: 110.597
Residual standard error = 0.43743 on 578 degrees of freedom

Plans for the near future

  • Nonlinear mixed-effects models fit using the Laplace approximation to the log-likelihood

  • Specification of partially linear models

  • Composite models consisting of a parameter transformation and a nonlinear model.

Partially linear models (this used to work but is now broken)

Partially linear models (those models with some parameters that occur linearly in the model expression) are expressed as types that inherit from the PLregMod abstract type. A instance of a model type is created from the values of any covariates in the model.

Example - a Michaelis-Menten fit

In the Michaelis-Menten model for enzyme kinetics,

v = Vm * c / (K + c)

the relationship between the velocity, v, of a reaction and the concentration, c, of the substrate depends on two parameters; Vm, the maximum velocity and K, the Michaelis parameter. The Vm parameter occurs linearly in this expression whereas K is a nonlinear parameter.

To fit such a model we create a MicMen object from the vector of observed concentrations and a PLregFit object from this model and the responses.

julia> pur = dataset("datasets","Puromycin");

julia> purtrt = sub(pur, pur[:State] .== "treated")
12x3 SubDataFrame{Array{Int64,1}}
|-------|------|------|-----------|
| Row # | Conc | Rate | State     |
| 1     | 0.02 | 76   | "treated" |
| 2     | 0.02 | 47   | "treated" |
| 3     | 0.06 | 97   | "treated" |
| 4     | 0.06 | 107  | "treated" |
| 5     | 0.11 | 123  | "treated" |
| 6     | 0.11 | 139  | "treated" |
| 7     | 0.22 | 159  | "treated" |
| 8     | 0.22 | 152  | "treated" |
| 9     | 0.56 | 191  | "treated" |
| 10    | 0.56 | 201  | "treated" |
| 11    | 1.1  | 207  | "treated" |
| 12    | 1.1  | 200  | "treated" |

julia> pm1 = fit(MicMen(Rate ~ Conc, purtrt), true)
Iteration:  1, rss = 1679.58, cvg = 0.257787 at [201.837,0.0484065]
   Incr: [0.012547]  f = 1.0, rss = 1211.92
Iteration:  2, rss = 1211.92, cvg = 0.0122645 at [210.623,0.0609536]
   Incr: [0.00280048]  f = 1.0, rss = 1195.66
Iteration:  3, rss = 1195.66, cvg = 0.000161307 at [212.448,0.063754]
   Incr: [0.000331041]  f = 1.0, rss = 1195.45
Iteration:  4, rss = 1195.45, cvg = 1.56158e-6 at [212.661,0.0640851]
   Incr: [3.27084e-5]  f = 1.0, rss = 1195.45
Iteration:  5, rss = 1195.45, cvg = 1.4557e-8 at [212.681,0.0641178]
   Incr: [3.15934e-6]  f = 1.0, rss = 1195.45
Iteration:  6, rss = 1195.45, cvg = 1.35192e-10 at [212.684,0.0641209]
   Incr: [3.04477e-7]  f = 1.0, rss = 1195.45

Nonlinear least squares fit to 12 observations

      Estimate  Std.Error t value Pr(>|t|)
Vm     212.684    6.94715 30.6145  3.2e-11
K    0.0641212 0.00828095 7.74323   1.6e-5

Residual sum of squares at estimates: 1195.45
Residual standard error = 10.9337 on 10 degrees of freedom
julia> using RDatasets, NLreg

julia> purtrt = sub(dataset("datasets","Puromycin"),:(State .== "treated"));

julia> pm1 = fit(MicMen(Conc ~ Rate, purtrt),true)
Iteration:  1, rss = 0.188744, cvg = 0.0888416 at [-0.0786133,-220.728]
   Incr: [-4.66305]  f = 1.0, rss = 0.173277
Iteration:  2, rss = 0.173277, cvg = 0.00102418 at [-0.0995101,-225.391]
   Incr: [-0.688546]  f = 1.0, rss = 0.173117
Iteration:  3, rss = 0.173117, cvg = 6.54049e-6 at [-0.10249,-226.08]
   Incr: [0.0574836]  f = 1.0, rss = 0.173116
Iteration:  4, rss = 0.173116, cvg = 7.53229e-8 at [-0.102242,-226.022]
   Incr: [-0.00614653]  f = 1.0, rss = 0.173116
Iteration:  5, rss = 0.173116, cvg = 8.30647e-10 at [-0.102269,-226.028]
   Incr: [0.000645718]  f = 1.0, rss = 0.173116

Nonlinear least squares fit to 12 observations

      Estimate Std.Error  t value Pr(>|t|)
Vm   -0.102266 0.0315309 -3.24335   0.0088
K     -226.028   7.08463 -31.9039  2.2e-11

Residual sum of squares at estimates: 0.173116
Residual standard error = 0.131574 on 10 degrees of freedom

We can also use parameter transformations

julia> pm2 = fit([LogTr] * MicMen(Rate ~ Conc, purtrt), true)
Iteration:  1, rss = 1679.58, cvg = 0.257787 at [201.837,-3.02812]
   Incr: [0.259201]  f = 1.0, rss = 1198.55
Iteration:  2, rss = 1198.55, cvg = 0.00234006 at [211.785,-2.76892]
   Incr: [0.0198568]  f = 1.0, rss = 1195.48
Iteration:  3, rss = 1195.48, cvg = 2.12492e-5 at [212.598,-2.74906]
   Incr: [0.00188326]  f = 1.0, rss = 1195.45
Iteration:  4, rss = 1195.45, cvg = 1.96812e-7 at [212.675,-2.74718]
   Incr: [0.000181183]  f = 1.0, rss = 1195.45
Iteration:  5, rss = 1195.45, cvg = 1.82666e-9 at [212.683,-2.747]
   Incr: [1.74545e-5]  f = 1.0, rss = 1195.45

Nonlinear least squares fit to 12 observations

        Estimate Std.Error  t value Pr(>|t|)
Vm       212.684   6.94715  30.6145  3.2e-11
log(K)  -2.74698  0.129145 -21.2705   1.2e-9

Residual sum of squares at estimates: 1195.45
Residual standard error = 10.9337 on 10 degrees of freedom

Creating a PLregMod type

A PLregMod type contains the transposed gradient, usually called tgrad with the conditionally linear parameters first, the three-dimensional Jacobian array, usually called MMD, with each face corresponding to the partial derivative of the conditionally linear rows of tgrad with respect to the nonlinear parameters, and the values of any covariates needed to evaluate the model. The model-matrix function, mmf, is a function of two read-only arguments; nlp, the nonlinear parameters in the model function as a vector and x the covariates for a single observation, also as a vector, and two arguments, tg and MMD that are updated in the function. The tg vector is updated with the model matrix for the conditionally linear parameters and the MMD slice, considered as a matrix, is updated with the gradient.

For the Michaelis-Menten model the model-matrix function is

function MicMenmmf(nlp,x,tg,MMD)
    x1 = x[1]
    denom = nlp[1] + x1
    MMD[1,1] = -(tg[1] =  x1/denom)/denom
end

The arguments are untyped, to allow for submatrices or views of matrices and arrays, but they should be treated as three vectors and a matrix, for the purposes of indexing.

First Commit

09/12/2013

Last Touched

11 months ago

Commits

71 commits

Used By: