This dataset comes from a study of cold tolerance of a grass
species, *Echinochloa crus-galli*, described in Potvin et al. (1990)
and Pinheiro and Bates (2000). A total of twelve four-week-old
`plants` were used in the study. There were two `types` of
plants: six from Quebec and six from Mississippi. Two `treatments`,
nonchilling and chilling, were assigned to three plants of each type.
Nonchilled plants were kept at C and chilled plants were
subject to 14 hours of chilling at C. After 10 hours of recovery
at C, CO `uptake` rates (in ) were measured
for each plant at seven `concentrations` of ambient CO in increasing,
consecutive order. Plots of observations are shown in Figure
as circles. The objective of the experiment was to evaluate the effect
of plant `type` and chilling `treatment` on the CO
`uptake`.

Pinheiro and Bates (2000) gave detailed analyses of
this dataset based on NLMMs using their software `nlme`.
They reached the following model:

where denotes the CO uptake rate of

> options(contrasts=rep("contr.treatment", 2)) > co2.fit1 <- nlme(uptake~exp(a1)*(1-exp(-exp(a2)*(conc-a3))), fixed=list(a1+a2~Type*Treatment,a3~1), random=a1~1, groups=~Plant, data=CO2, start=c(log(30),0,0,0,log(0.01),0,0,0,50)) > co2.fit1 Nonlinear mixed-effects model fit by maximum likelihood Model: uptake ~ exp(a1) * (1 - exp(-exp(a2) * (conc - a3))) Data: CO2 AIC BIC logLik 393.2869 420.0259 -185.6434 Random effects: Formula: a1 ~ 1 | Plant a1.(Intercept) Residual StdDev: 0.08221494 1.857658 Fixed effects: list(a1 + a2 ~ Type * Treatment, a3 ~ 1) Value Std.Error DF t-value p-value a1.(Intercept) 3.73338 0.052536 64 71.06301 <.0001 a1.TypeMississippi -0.29080 0.075535 64 -3.84990 0.0003 a1.Treatmentchilled -0.07274 0.074633 64 -0.97459 0.3334 a1.TypeMississippi:Treatmentchilled -0.51321 0.109733 64 -4.67688 <.0001 a2.(Intercept) -4.57570 0.086116 64 -53.13387 <.0001 a2.TypeMississippi -0.09635 0.103950 64 -0.92687 0.3575 a2.Treatmentchilled -0.17048 0.091377 64 -1.86565 0.0667 a2.TypeMississippi:Treatmentchilled 0.70555 0.205851 64 3.42748 0.0011 a3 49.98833 4.576255 64 10.92341 <.0001 ...

Fits of model () are plotted in Figure as dotted lines. Based on model (), one may conclude that the CO2 uptake is higher for plants from Quebec and that chilling, in general, results in lower uptake, and its effect on Mississippi plants is much larger than on Quebec plants. These conclusions are comparable to the results in Potvin et al. (1990).

We aim to use this dataset to demonstrate how to fit SNM models with
covariates, and how to check if an NLMM is appropriate. As an extension of
(), we fit the following SNM model

where for some fixed and the second stage model is paralleled to (). In order to test if the parametric model () is appropriate, we construct the following -spline to model . The hypothesis is : . Let , where denotes the th derivative operator. The kernel space of is . Define a linear operator such that . Let . Define inner products on , and as , and . Then it is easy to check that . The reproducing kernel of is given in ().

Note that in () is excluded from () to make free of constraint on the vertical scale. We need the side conditions that and for to separate from . The first condition reduces to and is satisfied by all functions in . We do not enforce the second condition because it is satisfied by all reasonable estimates. Thus the null space for becomes and the model space is . The penalty is still .

With the initial values chosen from the NLMM fit, our program converged after 5 iterations.

> M <- model.matrix(~Type*Treatment, data=CO2)[,-1] > co2.fit2 <- snm(uptake~exp(a1)*f(exp(a2)*(conc-a3)), func=f(u)~list(~I(1-exp(-u))-1,lspline(u, type="exp")), fixed=list(a1~M-1,a3~1,a2~Type*Treatment), random=list(a1~1), group=~Plant, verbose=T, start=co2.fit1$coe$fixed[c(2:4,9,5:8)], data=CO2) > summary(co2.fit2) Semi-parametric Nonlinear Mixed Effects Model fit Model: uptake ~ exp(a1) * f(exp(a2) * (conc - a3)) Data: CO2 AIC BIC logLik 406.4864 441.625 -188.3760 Random effects: Formula: a1 ~ 1 | Plant a1.(Intercept) Residual StdDev: 0.09304172 1.816200 Fixed effects: list(a1 ~ M - 1, a3 ~ 1, a2 ~ Type * Treatment) Value Std.Error DF t-value p-value a1.MTypeMississippi -0.28569 0.055952 65 -5.10591 <.0001 a1.MTreatmentchilled -0.07212 0.054741 65 -1.31739 0.1923 a1.MTypeMississippi:Treatmentchilled -0.54879 0.099184 65 -5.53309 <.0001 a3 50.67565 4.226094 65 11.99113 <.0001 a2.(Intercept) -4.56698 0.085130 65 -53.64708 <.0001 a2.TypeMississippi -0.12162 0.098575 65 -1.23377 0.2217 a2.Treatmentchilled -0.16161 0.087009 65 -1.85740 0.0678 a2.TypeMississippi:Treatmentchilled 0.81924 0.211587 65 3.87187 0.0003 ... GCV estimate(s) of smoothing parameter(s): 1.864811 Equivalent Degrees of Freedom (DF): 4.867178 Converged after 5 iterations

Fits of model () are plotted in Figure as solid lines. Since the data set is small, different initial values may lead to different estimates. However, the overall fits are similar. We also fitted models with AR(1) within-subject correlations and covariate effects on . None of these models improve fits significantly. The estimates are comparable to the nonlinear fits and the conclusion is similar to that based on ().

To check if the parametric NLMM () is appropriate, we
calculated approximate posterior means and variances using the
function `intervals`. Then we plotted the estimated (overall), its
projection onto (the parametric part) and
(the smooth part) in Figure .

> co2.grid2 <- data.frame(u=seq(0.3, 11, len=50)) > co2.ci <- intervals(co2.fit2, newdata=co2.grid2, terms=matrix(c(1,1,1,0,0,1), ncol=2,byrow=T)) > plot.bCI(co2.ci,x=co2.grid2$u,layout=c(3,1), type.name=c("overall","parametric","smooth"))

For -splines, Heckman and Ramsay (2000) showed that selecting the right form of penalty function via a differential operator can reduce the bias. An -spline allows us to incorporate prior knowledge on the main features of into the penalty. Usually the form of is known with coefficients depending on some unknown parameters. Heckman and Ramsay (2000) proposed several methods to estimate these parameters. When the whole model can be written in the form of an SNM model as in this example, our methods can also be used to estimate the penalty. Our approach is the same as the PL (penalized likelihood) approach in Heckman and Ramsay (2000). They commented that the PL method is ``philosophically appealing''. However it is also ``time-consuming'' because a grid search was used. Our method estimates the coefficients in and the functions iteratively, thus making it less computationally intensive. In addition, we allow these coefficients to depend on covariates.