smspline: Smoothing Splines in NLME

In lmeSplines: Add Smoothing Spline Modelling Capability to `nlme`

Smoothing Splines in NLME

Description

Functions to generate matrices for a smoothing spline covariance structure, enabling the fitting of smoothing spline terms in linear mixed-effects models (LME) or nonlinear mixed-effects models (NLME). A smoothing spline can be represented as a mixed model, as described by Speed (1991) and Verbyla (1999). The generated Z-matrix can be incorporated into a data frame and used in LME random effects terms with an identity covariance structure (pdIdent(~Z - 1)).

The model formulation for a spline in time (t) is:

y = X_s \beta_s + Z_s u_s + e

where X_s = [1 | t], Z_s = Q (t(Q) Q)^{-1}, and u_s \sim N(0, G_s) is a set of random effects. The random effects are transformed to independence via u_s = L v_s, where v_s \sim N(0, I \sigma^2_s) and L is the lower triangle of the Cholesky decomposition of G_s. The Z-matrix is transformed to Z = Z_s L.

Usage

smspline(formula, data)

smspline.v(time)

Arguments

formula	Model formula with the right-hand side specifying the spline covariate (e.g., ~ time). Must contain exactly one variable.
data	Optional data frame containing the variable specified in formula. If not provided, the formula is evaluated in the current environment.
time	Numeric vector of spline time covariate values to smooth over.

Value

For smspline, a Z-matrix with the same number of rows as the input data frame or vector, representing the random effects design matrix for the smoothing spline. After fitting an LME model, the standard deviation parameter for the random effects estimates \sigma_s, and the smoothing parameter is \lambda = \sigma^2 / \sigma^2_s.

For smspline.v, a list containing:

Xs

Matrix for fixed effects, with columns [1 | t].

Zs

Matrix for random effects, computed as Q (t(Q) %*% Q)^-1 L.

Q

Matrix used in the spline formulation.

Gs

Covariance matrix for the random effects.

R

Cholesky factor of Gs.

Note

The time points for the smoothing spline basis are, by default, the unique values of the time covariate. Model predictions at the fitted data points can be obtained using predict.lme. For predictions at different time points, use approx.Z to interpolate the Z-matrix.

Author(s)

Rod Ball <[e-mail]>

References

Pinheiro, J. and Bates, D. (2000) Mixed-Effects Models in S and S-PLUS. Springer-Verlag, New York.

Speed, T. (1991) Discussion of "That BLUP is a good thing: the estimation of random effects" by G. Robinson. Statistical Science, 6, 42–44.

Verbyla, A. (1999) Mixed Models for Practitioners. Biometrics SA, Adelaide.

See Also

approx.Z, lme

Examples

# Smoothing spline curve fit
data(smSplineEx1)
smSplineEx1$all <- rep(1, nrow(smSplineEx1))
smSplineEx1$Zt <- smspline(~ time, data = smSplineEx1)
fit1s <- lme(y ~ time, data = smSplineEx1,
             random = list(all = pdIdent(~ Zt - 1)))
summary(fit1s)
plot(smSplineEx1$time, smSplineEx1$y, pch = "o", type = "n",
     main = "Spline fits: lme(y ~ time, random = list(all = pdIdent(~ Zt - 1)))",
     xlab = "time", ylab = "y")
points(smSplineEx1$time, smSplineEx1$y, col = 1)
lines(smSplineEx1$time, smSplineEx1$y.true, col = 1)
lines(smSplineEx1$time, fitted(fit1s), col = 2)

# Fit model with reduced number of spline points
times20 <- seq(1, 100, length = 20)
Zt20 <- smspline(times20)
smSplineEx1$Zt20 <- approx.Z(Zt20, times20, smSplineEx1$time)
fit1s20 <- lme(y ~ time, data = smSplineEx1,
               random = list(all = pdIdent(~ Zt20 - 1)))
anova(fit1s, fit1s20)
summary(fit1s20)

# Model predictions on a finer grid
times200 <- seq(1, 100, by = 0.5)
pred.df <- data.frame(all = rep(1, length(times200)), time = times200)
pred.df$Zt20 <- approx.Z(Zt20, times20, times200)
yp20.200 <- predict(fit1s20, newdata = pred.df)
lines(times200, yp20.200 + 0.02, col = 4)

# Mixed model spline terms at multiple levels of grouping
data(Spruce)
Spruce$Zday <- smspline(~ days, data = Spruce)
Spruce$all <- rep(1, nrow(Spruce))
spruce.fit1 <- lme(logSize ~ days, data = Spruce,
                   random = list(all = pdIdent(~ Zday - 1),
                                 plot = ~ 1, Tree = ~ 1))
spruce.fit2 <- lme(logSize ~ days, data = Spruce,
                   random = list(all = pdIdent(~ Zday - 1),
                                 plot = pdBlocked(list(~ days, pdIdent(~ Zday - 1))),
                                 Tree = ~ 1))
anova(spruce.fit1, spruce.fit2)
summary(spruce.fit1)

lmeSplines documentation built on Aug. 31, 2025, 5:08 p.m.