Data-driven linear mixed effect model spline modelling


Function that models a linear or limear mixed model depending on the best fit. Alternatively, the function can return THE derivation information of the fitted models for the fixed (original) times points and a chosen basis.


lmmSpline(data, time, sampleID, timePredict, deri, basis, knots, keepModels,numCores)



data.frame or matrix containing the samples as rows and features as columns


numeric vector containing the sample time point information.


character, numeric or factor vector containing information about the unique identity of each sample


numeric vector containing the time points to be predicted. By default set to the original time points observed in the experiment.


logical value. If TRUE returns the predicted derivative information on the observed time points.By default set to FALSE.


character string. What type of basis to use, matching one of "cubic", "p-spline" or "cubic p-spline". The "cubic" basis (default) is the cubic smoothing spline as defined by Verbyla et al. 1999, the "p-spline" is the truncated p-spline basis as defined by Durban et al. 2005.


Alternatively an integer, the number of knots used for the "p-spline" or "cubic p-spline" basis calculation. Otherwise calculated as proposed by Ruppert 2002. Not used for the "cubic" smoothing spline basis as it used the inner design points.


alternative logical value if you want to keep the model output. Default value is FALSE


Alternative numeric value indicating the number of CPU cores to be used. Default value is automatically estimated.


The first model (modelsUsed=0) assumes the response is a straight line not affected by individual variation.

Let y_{ij}(t_{ij}) be the expression of a feature for individual (or biological replicate) i at time t_{ij}, where i=1,2,...,n, j=1,2,...,m_i, n is the sample size and m_i is the number of observations for individual i for the given feature. We fit a simple linear regression of expression y_{ij}(t_{ij}) on time t_{ij}. The intercept β_0 and slope β_1 are estimated via ordinary least squares: y_{ij}(t_{ij})= β_0 + β_1 t_{ij} + ε_{ij}, where ε_{ij} ~ N(0,σ^2_{ε}). The second model (modelsUsed=1) is nonlinear where the straight line in regression replaced with a curve modelled using here for example a spline truncated line basis (basis="p-spline") as proposed Durban et al. 2005:

y_{ij}(t_{ij})= f(t_{ij}) +ε_{ij},

where ε_{ij}~ N(0,σ_{ε}^2).

The penalized spline is represented by f, which depends on a set of knot positions κ_1,...,κ_K in the range of {t_{ij}}, some unknown coefficients u_k, an intercept β_0 and a slope β_1. The first term in the above equation can therefore be expanded as:

f(t_{ij})= β_0+ β_1t_{ij}+∑\limits_{k=1}^{K}u_k(t_{ij}-κ_k)_+,

with (t_{ij}-κ_k)_+=t_{ij}-κ_k, if t_{ij}-κ_k > 0, 0 otherwise.

The choice of the number of knots K and their positions influences the flexibility of the curve. If the argument knots=missing, we use a method proposed by Ruppert 2002 to estimate the number of knots given the measured number of time points T, so that the knots κ_1 … κ_K are placed at quantiles of the time interval of interest:

K= max(5,min(floor(\frac{T}{4}) , 40)).

In order to account for individual variation, our third model (modelsUsed=2) adds a subject-specific random effect U_i to the mean response f(t_{ij}). Assuming f(t_{ij}) to be a fixed (yet unknown) population curve, U_i is treated as a random realization of an underlying Gaussian process with zero-mean and variance σ_U^2 and is independent from the random error ε_{ij}:

y_{ij}(t_{ij}) = f(t_{ij}) + U_i + ε_{ij}

with U_{i} ~ N(0,σ_U^2) and ε_{ij} ~ N(0,σ_{ε}^2). In the equation above, the individual curves are expected to be parallel to the mean curve as we assume the individual expression curves to be constant over time. A simple extension to this model is to assume individual deviations are straight lines. The fourth model (modelsUsed=3) therefore fits individual-specific random intercepts a_{i0} and slopes a_{i1}:

y_{ij}(t_{ij}) = f(t_{ij}) + a_{i0} + a_{i1}t_{ij} + ε_{ij}

with ε_{ij} ~ N(0,σ_ε^2) and (a_{i0},a_{i1})^T ~ N(0,Σ). We assume independence between the random intercept and slope. @return lmmSpline returns an object of class lmmspline containing the following components:

  • predSplinedata.frame containing predicted values based on linear model object or linear mixed effect model object.

  • modelsUsednumeric vector indicating the model used to fit the data. 0 = linear model, 1=linear mixed effect model spline (LMMS) with defined basis ('cubic' by default) 2 = LMMS taking subject-specific random intercept, 3 = LMMS with subject specific intercept and slope.

  • modellist of models used to model time profiles.

  • derivativelogical value indicating if the predicted values are the derivative information.


Durban, M., Harezlak, J., Wand, M. P., & Carroll, R. J. (2005). Simple fitting of subject-specific curves for longitudinal data. Stat. Med., 24(8), 1153-67.

Ruppert, D. (2002). Selecting the number of knots for penalized splines. J. Comp. Graph. Stat. 11, 735-757

Verbyla, A. P., Cullis, B. R., & Kenward, M. G. (1999). The analysis of designed experiments and longitudinal data by using smoothing splines. Appl.Statist, 18(3), 269-311.

Straube J., Gorse A.-D., Huang B.E., Le Cao K.-A. (2015). A linear mixed model spline framework for analyzing time course 'omics' data PLOSONE, 10(8), e0134540.

See Also

summary.lmmspline, plot.lmmspline, predict.lmmspline, deriv.lmmspline


## Not run: 
# running for samples in group 1
G1 <- which(kidneySimTimeGroup$group=="G1")
testLMMSpline<- lmmSpline(data=kidneySimTimeGroup$data[G1,],time=kidneySimTimeGroup$time[G1],
DerivTestLMMSplineTG<- lmmSpline($data[G1,]),
## End(Not run)

Want to suggest features or report bugs for Use the GitHub issue tracker. Vote for new features on Trello.

comments powered by Disqus