Description Usage Arguments Details Value Note Author(s) References See Also Examples
Functions to generate matrices for a smoothing spline covariance structure, to enable fitting smoothing spline terms in LME/NLME.
1 2 |
formula |
model formula with right hand side giving the spline covariate |
data |
optional data frame |
time |
spline ‘time’ covariate to smooth over |
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.
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. |
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= ...)
Rod Ball rod.ball@scionresearch.com www.scionresearch.com
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 | # 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)
|
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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.