5

4

5

2

# Nonlinear regression models in Julia

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> 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.

09/12/2013

11 months ago

71 commits