smspline: Smoothing splines in NLME

In lmeSplines: Add smoothing spline modelling capability to nlme.

Description

Functions to generate matrices for a smoothing spline covariance structure, to enable fitting smoothing spline terms in LME/NLME.

Usage

  smspline(formula, data)
  smspline.v(time)

Arguments

formula	model formula with right hand side giving the spline covariate
data	optional data frame
time	spline ‘time’ covariate to smooth over

Details

A smoothing spline can be represented as a mixed model (Speed 1991, Verbyla 1999). The generated Z-matrix from smspline() can be incorporated in the users's dataframe, then used in model formulae for LME random effects terms at any level of grouping (see examples). The spline random terms are fitted in LME using an identity 'pdMat' structure of the form pdIdent(~Z - 1). The model formulation for a spline in time (t) is as follows (Verbyla 1999):

y = X_s β_s + Z_s u_s + e

where X_s = [1 | t] , Z_s = Q(t(Q) Q)^{-1} , and u_s ~ N(0,G_s), is a set of random effects. We transform the set of random effects u_s to independence with u_s = L v_s, where

v_s ~ N(0,I σ^2_s)

is a set of independent random effects. The Z-matrix is transformed accordingly to Z = Z_s L, where L is the lower triangle of the Choleski decomposition of G_s.

The function smspline.v() is called by smspline(), and can be used to access the matrices X_s, Q, G_s. See Verbyla (1999) for further information.

Value

For smspline(), a Z-matrix with the same number of rows as the data frame. After fitting, the LME model output gives a standard deviation parameter for the random effects, estimating σ_s. The smoothing parameter from the penalised likelihood formulation is

λ = σ^2/σ^2_s

For smspline.v(), a list of the form

Xs	X-matrix for fixed effects part of the model
Zs	Z-matrix for random effects part of the model
Q,Gs,R	Matrices Q, G_s, R associated to the mixed-model form of the smoothing spline.

Note

The time points for the smoothing spline basis are, by default, the unique values of the time covariate. This is the easiest approach, and model predictions at the fitted data points, can be obtained using predict.lme. By interpolation, using approx.Z, the Z-matrix can be obtained for any set of time points and can be used for fitting and/or prediction. (See examples). Synopsis:data$Z <- smspline(formula1, data); fit <-lme(formula2, data, random= ...)

Author(s)

Rod Ball [e-mail] www.scionresearch.com

References

The correspondence between penalized likelihood formulations of smoothing splines and mixed models was pointed out by Speed (1991). The formulation used here for the mixed smoothing spline matrices are given in Verbyla (1999). LME/NLME modelling is introduced in Pinheiro and Bates (2000).

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. Statist. Sci., 6, 42–44.

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

See Also

approx.Z nlme

Examples

# smoothing spline curve fit
data(smSplineEx1)
# variable `all' for top level grouping
smSplineEx1$all <- rep(1,nrow(smSplineEx1))
# setup spline Z-matrix
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 cut down 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)))
# note: virtually identical df, loglik.
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))
# overall spline term, random plot and Tree effects
spruce.fit1 <- lme(logSize ~ days, data=Spruce,
                   random=list(all= pdIdent(~Zday -1),
                     plot=~1, Tree=~1))
# try overall spline term plus plot level linear + spline term
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)

Example output

Loading required package: nlme
Linear mixed-effects model fit by REML
 Data: smSplineEx1 
       AIC     BIC   logLik
  281.9241 292.264 -136.962

Random effects:
 Formula: ~Zt - 1 | all
 Structure: Multiple of an Identity
               Zt1        Zt2        Zt3        Zt4        Zt5        Zt6
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
               Zt7        Zt8        Zt9       Zt10       Zt11       Zt12
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt13       Zt14       Zt15       Zt16       Zt17       Zt18
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt19       Zt20       Zt21       Zt22       Zt23       Zt24
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt25       Zt26       Zt27       Zt28       Zt29       Zt30
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt31       Zt32       Zt33       Zt34       Zt35       Zt36
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt37       Zt38       Zt39       Zt40       Zt41       Zt42
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt43       Zt44       Zt45       Zt46       Zt47       Zt48
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt49       Zt50       Zt51       Zt52       Zt53       Zt54
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt55       Zt56       Zt57       Zt58       Zt59       Zt60
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt61       Zt62       Zt63       Zt64       Zt65       Zt66
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt67       Zt68       Zt69       Zt70       Zt71       Zt72
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt73       Zt74       Zt75       Zt76       Zt77       Zt78
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt79       Zt80       Zt81       Zt82       Zt83       Zt84
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt85       Zt86       Zt87       Zt88       Zt89       Zt90
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt91       Zt92       Zt93       Zt94       Zt95       Zt96
StdDev: 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436 0.01627436
              Zt97       Zt98  Residual
StdDev: 0.01627436 0.01627436 0.8603059

Fixed effects: y ~ time 
               Value  Std.Error DF  t-value p-value
(Intercept) 6.498800 0.17335977 98 37.48736       0
time        0.038723 0.00298034 98 12.99296       0
 Correlation: 
     (Intr)
time -0.868

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.8893014 -0.5811348  0.1164587  0.6593132  1.7839082 

Number of Observations: 100
Number of Groups: 1 
        Model df      AIC      BIC    logLik
fit1s       1  4 281.9241 292.2640 -136.9621
fit1s20     2  4 281.9491 292.2889 -136.9745
Linear mixed-effects model fit by REML
 Data: smSplineEx1 
       AIC      BIC    logLik
  281.9491 292.2889 -136.9745

Random effects:
 Formula: ~Zt20 - 1 | all
 Structure: Multiple of an Identity
             Zt201      Zt202      Zt203      Zt204      Zt205      Zt206
StdDev: 0.01613004 0.01613004 0.01613004 0.01613004 0.01613004 0.01613004
             Zt207      Zt208      Zt209     Zt2010     Zt2011     Zt2012
StdDev: 0.01613004 0.01613004 0.01613004 0.01613004 0.01613004 0.01613004
            Zt2013     Zt2014     Zt2015     Zt2016     Zt2017     Zt2018
StdDev: 0.01613004 0.01613004 0.01613004 0.01613004 0.01613004 0.01613004
         Residual
StdDev: 0.8618988

Fixed effects: y ~ time 
               Value  Std.Error DF  t-value p-value
(Intercept) 6.392423 0.17557148 98 36.40923       0
time        0.039741 0.00302297 98 13.14632       0
 Correlation: 
     (Intr)
time -0.87 

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.8961970 -0.5735952  0.1249489  0.6647801  1.7871074 

Number of Observations: 100
Number of Groups: 1 
            Model df       AIC       BIC    logLik   Test  L.Ratio p-value
spruce.fit1     1  6 -178.6251 -149.0304  95.31255                        
spruce.fit2     2  9 -230.3999 -186.0079 124.19997 1 vs 2 57.77485  <.0001
Linear mixed-effects model fit by REML
 Data: Spruce 
        AIC       BIC   logLik
  -178.6251 -149.0304 95.31255

Random effects:
 Formula: ~Zday - 1 | all
 Structure: Multiple of an Identity
               Zday1        Zday2        Zday3        Zday4        Zday5
StdDev: 0.0007578212 0.0007578212 0.0007578212 0.0007578212 0.0007578212
               Zday6        Zday7        Zday8        Zday9       Zday10
StdDev: 0.0007578212 0.0007578212 0.0007578212 0.0007578212 0.0007578212
              Zday11
StdDev: 0.0007578212

 Formula: ~1 | plot %in% all
        (Intercept)
StdDev:  0.08906169

 Formula: ~1 | Tree %in% plot %in% all
        (Intercept)  Residual
StdDev:   0.6217286 0.1765053

Fixed effects: logSize ~ days 
               Value  Std.Error  DF   t-value p-value
(Intercept) 4.140124 0.08522671 947  48.57778       0
days        0.003322 0.00002940 947 112.99157       0
 Correlation: 
     (Intr)
days -0.148

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-4.433907659 -0.482001385 -0.007197313  0.535760074  6.674697150 

Number of Observations: 1027
Number of Groups: 
                    all           plot %in% all Tree %in% plot %in% all 
                      1                       4                      79 

lmeSplines documentation built on May 2, 2019, 6:48 p.m.