A simple example of fitting a nonlinear regression and nonlinear mixed-effects model

First we load the NLreg package and and the DataFrames package. The readtable function in creates a DataFrame from a comma-separated or tab-separated file's contents. The file can be compressed, as shown here.

In []: 
using DataFrames,NLreg
const sd1 = readtable(Pkg.dir("NLreg","data","sd1.csv.gz"));

Nonlinear regression models are represented as types in the NLreg package. The type can implement several constructors. Here we generate a BolusSD1 model from a formula expression and a data frame.

In []: 
m = BolusSD1(:(CONC ~ TIME), sd1);

The formula indicates that the response is CONC and the first (and only) covariate is TIME.

Fitting a nonlinear regression model

A nonlinear least squares fit is represented by the type NonlinearLS, which is constructed from a nonlinear regression model and, optionally, initial values for the parameters.

If initial values for the parameters are not given the model's initpars method is used to generate initial values.

In [1]: 
show(initpars(m))

[0.9265302850601753,0.2793108696162248]

The gnfit function fits the model using the Gauss-Newton algorithm. An optional second argument determines if verbose output is generated.

In [2]: 
nl = gnfit(NonlinearLS(m),true)

Iteration:  1, rss = 122.661, cvg = 0.0980845 at [0.92653,0.279311]
   Incr: [0.21725,-0.0428423]  f = 1.0, rss = 110.69
Iteration:  2, rss = 110.69, cvg = 0.00081764 at [1.14378,0.236469]
   Incr: [-0.00211657,0.00838402]  f = 1.0, rss = 110.598
Iteration:  3, rss = 110.598, cvg = 2.64837e-6 at [1.14166,0.244853]
   Incr: [0.00119367,0.000779098]  f = 1.0, rss = 110.597
Iteration:  4, rss = 110.597, cvg = 1.086e-8 at [1.14286,0.245632]
   Incr: [9.1389e-5,5.07036e-5]  f = 1.0, rss = 110.597
Iteration:  5, rss = 110.597, cvg = 4.41408e-11 at [1.14295,0.245682]
   Incr: [5.83282e-6,3.23308e-6]  f = 1.0, rss = 110.597
Out[2]:
Nonlinear least squares fit to 580 observations

     Estimate Std.Error t value Pr(>|t|)
V     1.14295  0.049565 23.0596  < eps()
K    0.245682  0.020241 12.1378  < eps()

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

There are several extractor methods such as coef, coeftable, deviance, pnames, stderr and vcov for this type.

In [3]: 
coeftable(nl)
Out[3]:
     Estimate Std.Error t value Pr(>|t|)
V     1.14295  0.049565 23.0596  < eps()
K    0.245682  0.020241 12.1378  < eps()
In [4]: 
deviance(nl)  # residual sum-of-squares
Out[4]:
110.59728633992346

Fitting a simple nonlinear mixed-effects model

A simple nonlinear mixed-effects model has a random effect for each parameter for each group (e.g. subject). Uncorrelated random effects are indicated by specifying a Diagonal \(\Lambda\) matrix. Correlated random effects user a Triangular \(\Lambda\) matrix. The model form and initial estimates of the fixed-effects parameters can be taken from a fitted NonlinearLS object.

In [6]: 
nm = NLreg.fit(SimpleNLMM(nl,vector(sd1[:ID]),Diagonal))
Out[6]:
Simple, nonlinear mixed-effects model fit by Laplace approximation
 logLik: -102.527190, deviance: 205.054379

 Variance components:
                Variance    Std.Dev.
 V              0.494279    0.703050
 K              0.013330    0.115457
 Residual       0.023797    0.154264
 Number of obs: 580; levels of grouping factor: 200

Fixed effects parameters:
     Estimate
V     1.19333
K    0.294698
In [7]: 
nm1 = NLreg.fit(SimpleNLMM(nl,vector(sd1[:ID]),Triangular))
Out[7]:
Simple, nonlinear mixed-effects model fit by Laplace approximation
 logLik: -102.167404, deviance: 204.334808

 Variance components:
                Variance    Std.Dev.
 V              0.484640    0.696161
 K              0.013369    0.115622
 Residual       0.023999    0.154916
 Number of obs: 580; levels of grouping factor: 200

Fixed effects parameters:
     Estimate
V     1.19672
K    0.299901

Although not currently shown in the output, the within-subject correlation of the random effects for this model fit is -0.091. A likelihood ratio test comparting nm versus nm1 has a p-value of 39.6% from which we would conclude that the more complex model, nm1, does not provide a sufficiently better fit to justify the additional parameter.