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

time 

sampleID 

timePredict 

deri 

basis 

knots 
Alternatively an 
keepModels 
alternative 
numCores 
Alternative 
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
="pspline") 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 subjectspecific 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 zeromean 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 individualspecific 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 subjectspecific 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 subjectspecific curves for longitudinal data. Stat. Med., 24(8), 115367.
Ruppert, D. (2002). Selecting the number of knots for penalized splines. J. Comp. Graph. Stat. 11, 735757
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), 269311.
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.
summary.lmmspline
, plot.lmmspline
, predict.lmmspline
, deriv.lmmspline
1 2 3 4 5 6 7 8 9 10 11 12  ## Not run:
data(kidneySimTimeGroup)
# running for samples in group 1
G1 < which(kidneySimTimeGroup$group=="G1")
testLMMSpline< lmmSpline(data=kidneySimTimeGroup$data[G1,],time=kidneySimTimeGroup$time[G1],
sampleID=kidneySimTimeGroup$sampleID[G1])
summary(testLMMSpline)
DerivTestLMMSplineTG< lmmSpline(data=as.data.frame(kidneySimTimeGroup$data[G1,]),
time=kidneySimTimeGroup$time[G1],sampleID=kidneySimTimeGroup$sampleID[G1],
deri=TRUE,basis="pspline")
summary(DerivTestLMMSplineTG)
## End(Not run)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.