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 [11]:
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 [12]:
m = BolusSD1(CONC ~ TIME, sd1);
typeof(m)
Out[12]:
This model is a concrete type inheriting from the abstract type PLregModF
, a partially-linear regression model function.
In [13]:
typeof(m) <: PLregModF
Out[13]:
formula
methodsThe formula indicates that the response is CONC
and the first (and only) covariate is TIME
. The names of the parameters for the model are provided by pnames
.
In [14]:
pnames(m)
Out[14]:
The first parameter, V
, representing the volume of distribution, is conditionally linear. The second parameter, K
, the rate constant is nonlinear.
A partially linear least squares fit is represented by the type PLinearLS
, 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 [15]:
show(initpars(m))
Note that the initial estimates are for the nonlinear parameters only.
The fit
function fits the model using the Golub-Pereyra algorithm. An optional second argument determines if verbose output is generated.
In [16]:
pl = fit(m,true)
Out[16]:
There are several extractor methods such as coef
, coeftable
, deviance
, pnames
, stderr
and vcov
for this type.
In [17]:
coeftable(pl) # the coefficient table shown above
Out[17]:
In [18]:
deviance(pl) # residual sum-of-squares
Out[18]:
In this model we used the elimination rate constant, K
, as a parameter. Suppose that instead we wish to use the logarithm of the rate constant as a parameter. We could create a new model type or we could combine the existing model with a parameter transformation.
In [19]:
pl2 = fit([LogTr]*BolusSD1(CONC ~ TIME, sd1))
Out[19]:
The two fits are equivalent in that they produce the same sum of squared residuals and the parameter estimates can be transformed from one scale to the other.
In [20]:
using Base.Test
@test_approx_eq_eps coef(pl) [coef(pl2)[1], exp(coef(pl2)[2])] 1e-5
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 [21]:
nm = fit(SimpleNLMM(pl2,array(sd1[:ID]),Diagonal))
In [22]:
nm1 = fit(SimpleNLMM(pl2,array(sd1[:ID]),Triangular))
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.