Generalized Smoothing Spline Models

Suppose that data have the form $(y_i, \mbox{\boldmath$t$}_i), i=1, \cdots, n$, where $y_i$'s are independent observations and $\mbox{\boldmath$t$}_i=(t_{1i}, \cdots, t_{di})$. The distribution of $y_i$ is from an exponential family with density function

$$g(y_i; f_i, \phi)={\mbox exp}\{y_ih(f_i)-b(f_i)/a(\phi) + c(y_i, \phi)\},$$ (25)

where $f_i=f(\mbox{\boldmath$t$}_i)$, $h(f_i)$ is a monotone transformation of $f_i$ known as the canonical link, and $\phi$ is a dispersion parameter. Assume that $f\in{\cal H}$ where ${\cal H}$ is given in (). The penalized likelihood estimate of $f$ is the minimizer of

$$-\sum_{i=1}^nl_i(f_i) + \frac{n\lambda}{2}\sum_{k=1}^p \theta_{k}^{-1} \vert\vert P_{k}f\vert\vert^2,$$ (26)

where $l_i$ is the log-likelihood of $y_i$. Again, the solution to () has the form () (Wahba et al., 1995), and $\mbox{\boldmath$c$}$ and $\mbox{\boldmath$d$}$ are solved by minimizing (). Usually the coefficients cannot be solved directly. If all $l_i(f_i)$'s are strictly concave, the Newton-Raphson iterative procedure can be used to calculate $\mbox{\boldmath$c$}$ and $\mbox{\boldmath$d$}$ for fixed smoothing parameters. The smoothing parameters $\lambda/\theta_1,\cdots,\lambda/\theta_p$ can be estimated at each iteration using GCV, GML and UBR methods (Gu, 1990; Wahba et al., 1995; Gu, 1992). It was found that when the dispersion parameter is known, the UBR method works better than the GCV and GML methods (Wang et al., 1995). For binary, binomial, Poisson and gamma data, this procedure was implemented in GRKPACK (Wang, 1997). In our ASSIST package, the functions gdsidr and gdmudr serve as intermediate interface between S and several drivers in GRKPACK.

The argument family in ssr specifies the distribution of $y$ as in glm. Families supported are ``binary'', ``binomial'', ``poisson'', ``gamma'' and ``gaussian'' for Bernoulli, binomial, Poisson, gamma and Gaussian distributions respectively. The default is Gaussian.

Laplace approximations to the posterior mean and variance can be calculated by the predict function (Wahba et al., 1995). Then Bayesian confidence intervals can be constructed.

For example, one may fit a cubic spline to binary data with the UBR choice of the smoothing parameter and compute approximate posterior means and variances by

    a <- ssr(y~t, rk=cubic(t), family=``binary'', spar=``u'', varht=1)
    predict(a)

where varht specifies fixed variance (dispersion) parameter as 1 for the UBR function.