A natural polynomial regression spline takes a form nearly identical to the piecewise regression spline seen above. It is different, however, in that it restricts the function to be linear at the boundaries. This sets the last term on the right-hand side of (2) to 0. A natural polynomial regression spline of order p thus takes the form If f (x) were modeled as such, (1) could then be re-written as where and β = (β₀, β₁, β₂, …, β_1+k), and for m observations collected at (x₁, x₂, …, x_m), it would be written as

The natural polynomial regression spline written in Equation (3) illustrates the concept behind splines, but the design matrix composed of the row vectors c (x₁, t_k), …, c (x_m, t_k) is unstable and thus rarely used. B-splines, which are computationally more stable, are preferred (Hastie and Tibshirani, 1990; Zhou and Shen, 2001). The functional form of a b-spline is more complex than that of the natural polynomial regression spline, but is readily available (deBoor, 1978). We write it as where b(x, t_k) is the design vector associated with a b-spline evaluated at x with k knots at t_k = (t₁, t₂, …, t_k), and α is a fixed set of parameters. For m observations collected at (x₁, …, x_m), the model would be written as With this model, estimating f(x) becomes a problem of estimating the number of knots, k, the locations of those knots, t_k, and α.

A variety of methods have been developed in estimating k and t_k. Halpern (1973) approached this problem using Bayesian methods. Allowing knots to only be placed at the design points in the experiment, he considered all of the subsets of the design points. Halpern assigned prior probabilities to all of these subsets, and calculated the corresponding posterior probabilities of these subsets. His estimator of the function was based on these posterior probabilities. Denison, Mallick, and Smith (1998) placed priors on the number of knots, k, and their locations, t_k. With these priors, they calculated a joint posterior distribution which included the variables k and t_k and then sampled from this posterior distribution using reversible jump MCMC methods. They too restricted the knots to be located only at the design points of the experiment. DiMatteo, Genovese, and Kass (2001) proposed a method similar to that of Denison, et al. They did not restrict the knots to be located only at the design points of the experiment and they penalized models with unnecessarily large values of k by averaging over the fixed effect coefficients, α.

Friedman (1991) and Stone et al. (1997) try to solve the problem of knot selection using backward and/or forward knot selection. This process continues until the “best” model has been identified. Zhou and Shen (2001) used an alternative method in identifying the locations of the knots. They constructed an algorithm which favored adding knots in locations where several knots had already been added. Lindstrom (1999) used similar methods when selecting knot locations.

Splines have also been used to model curves that vary within and between subjects (Shi et al., 1996, Rice and Wu, 2001, Behseta et al., 2005, Crainiceanu et al., 2005). Brumback and Rice (1998) use splines to model nested functions which vary between subjects and between groups of subjects. Functional models which only account for variation between and within subjects take the general form where Y_i (x_j) is the observation of the i^th individual at x_j, f (x) can be thought of as a population curve, and G_i (x) is a random curve specific to subject i. While these functions can be modeled in a variety of ways, we model these two functions as free-knot b-splines with k and k′ knots located at t_k and t_k′, respectively. Setting f(x) = b(x; t_k)α and G_i(x) = b(x; t_k′)γ_i, where γ_i are independent random vectors such that γ_i ~ N(0, Σ_γ), equation (4) becomes Bigelow and Dunson (2007) consider a variation of this model and identify the location of the knots by sampling from the full posterior distribution. They lessen the computational burden of this sampling procedure by setting t_k′ = t_k and restricting Σ_γ to be diagonal. In this paper, we consider a more flexible model and explore the computational issues associated with it. To be more specific, we eliminate the restrictions imposed by Bigelow and Dunson and then try to identify the number and location of the knots by sampling from p(k, t_k, k′, t_k′|y). This posterior can not be sampled from exactly, however, because eliminating the diagonal restriction on Σ_γ makes the likelihood p (y|k, t_k, k′, t_k′) intractable. We thus sample from two approximations to this posterior; each approximation corresponds to a different approximation to the likelihood p (y|k, t_k, k′, t_k′). To assess the accuracy of these approximations, we compare the posterior distribution of knots in each case to the marginal distribution of knots when sampling from p(k, t_k, k′, t_k′, ψ, σ²|y), where ψ is a set of covariance parameters and σ² is the within-subject variability. The marginal distribution of knots when sampling from this higher-dimensional posterior is exact as its corresponding likelihood is available in closed form. Sampling from this larger posterior, however, is less efficient because of the additional parameters that need to be sampled.

We place prior distributions on k, k′, t_k, t_k′, α, Σ_γ, and σ². Poisson priors are given to k and k′. Additionally, since there is no reason to favor knots at any particular location on the domain of f and G, flat priors are placed on both t_k and t_k′. These priors are written as where t_q(j) is the j^th smallest knot in the vector t_q, and (a, b) is the domain of f and G.

For Σ_γ and σ², we place prior distributions on transformed values of these parameters. The transformations we consider create a set of unconstrained parameters which are preferable to work with since their sampling distributions can be well approximated with the normal density. For Σ_γ, we applied the modified Cholesky parameterization proposed by Pourhamadi (1999). Pourhamadi decomposes the inverse of the covariance matrix Σ_γ as , where These parameters can be calculated directly from the elements of the covariance matrix using the functions g (·;·) and h (·;·), where for 2 ≤ j ≤ p , ϕj = (ϕ_j,1, ϕ_j,2, …, ϕ_j,j−1)′ = h (Σ_γ;j) = , Σ_γ,j is the minor matrix within Σ_γ composed of its first j − 1 columns and rows, and σ_j is the vector within Σ_γ composed of the first j – 1 elements of the j^th column. To ease our notational burden, we will refer to these modified Cholesky parameters collectively as ψ, where , and we will replace the symbol Σ_γ with ψ (or, to illustrate that Σ_γ is a function of ψ, with Σ_γ(ψ)). All of the elements in the vector ψ are unconstrained, and a multivariate normal unit-information prior is assigned to them. The mean of the unit-information prior on ψ is the maximum likelihood estimate of the parameters, , and the variance is inversely proportional to approximately one unit of information on ψ. The prior can be written as where p(y|k, t_k, k′, t_k′, α, ψ, σ²) is the distribution of all the observed data in the study, n is the number of subjects in the study, and , , and are the maximum likelihood estimates of α, ψ, and σ, respectively. The exact formula for I_ψ is given in the Appendix.

For σ², a unit-information prior was placed on log (σ²). The prior can be written as where and m_i is the number of observations for subject i. The exact formula for I_log(σ²), i is given in the Appendix.

The prior placed on α is another unit-information prior and can be written as where The exact formula for I_α,i is given in the Appendix.

Kass and Wasserman (1994) and Pauler (1998) defend the use of unit-information priors, especially in the context of model selection. They show that when using such priors, the resulting Bayes factors can be well approximated with the Bayesian information criterion. Although centering unit-information priors at the maximum likelihood estimates has been suggested in Bayesian adaptive regression splines (Paciorek, 2006), the final results presented in this paper were not sensitive to the means selected for these three prior distributions.

The integral in (7) can not be calculated analytically. We thus explore two different approximations to p (y|k, t_k, k′, t_k′). The first simply plugs in the maximum likelihood estimate of ψ and σ². The second approximates (7) using a Laplace approximation.

The first approximation that we consider estimates p (y|k, t_k, k′, t_k′) with where (t_k, t_k′) and ² (t_k, t_k′) are the maximum likelihood estimators of ψ and σ² corresponding to the model with fixed effect knots at t_k and random effect knots at t_k′. The dependence of these maximum likelihood estimators on (t_k, t_k′) will be suppressed, and the estimators will now simply be denoted as and ². We refer to this likelihood approximation as the “plugged-in” approximation and denote it as _{Plugged In} (y|k, t_k, k′, t_k′). Although this approximation ignores the penalty from increasing the dimension of ψ, it is computationally inexpensive and is similar to the methods employed by Taplin (1993), Draper (1995), Raftery et al. (1996), and Dominici et al. (2002). They approximate marginal likelihoods by plugging in the maximum likelihood estimates of the nuisance parameters that can not be integrated out.

Rather than sampling the fixed and random effect knots together in one MCMC iteration (which results in very low acceptance rates), we sample from the posterior (k, t_k, k′, t_k′|y) using a Gibbs Sampling type algorithm. In the algorithm, one value of the pair (k, t_k) is sampled from (k, t_k|k′, t_k′, y) using the Metropolis-Hastings algorithm. Conditioned on this sampled pair of (k, t_k), a pair (k′, t_k′) is sampled from (k′, t_k′|k, t_k, y) using identical Metropolis-Hastings methods. This algorithm is given below, and will be referred to as Algorithm I.

The jump probabilities and details of this algorithm are given in the Appendix.

While LP will take longer to implement than PLG (as calculating _Lp (k, t_k, k′, t_k′|y) will take longer than calculating _Plg(k, t_k, k′, t_k′|y)), it is expected that PLG will result in unnecessarily large values of k′. We anticipate this positive bias in the number of random effect knots since the plugged-in approximation to p(y|k, t_k, k′, t_k′) does not average over the random effect covariance parameters, ψ; it simply plugs in the maximum likelihood estimate of ψ. The tendency of PLG to estimate more random effect knots than LP becomes clear after close examination of the ratio calculated in step 3(a) of Algorithm I. In PLG, this ratio is calculated as , and in LP, this ratio is calculated as . In general, we expect when . If this inequality holds, then it does follow that PLG will over estimate the number of random effect knots. In Theorem 1, it is proven that when is the true set of knots, the inequality given in (10) asymptotically holds.

Theorem 1. Let the true values of k, t_k, k′ and t_k′ be k_new, t_knew, , respectively. If the four conditions below are met, then as n (the total number of subjects) goes to ∞, The proof of this theorem is given in the Appendix.

LP is not only preferred to PLG, but also to EX. The Markov chain resulting from LP converges to the stationary distribution more quickly than that corresponding to EX. This is the case because, in LP, transitions are only made in the space of the fixed and random effect knots. In EX, transitions are made in the space of the fixed knots, the random effect knots, and the covariance parameters. The acceptance probability of the chain corresponding to LP is thus expected to be larger than that corresponding to EX.

The population curve and five randomly sampled subject specific curves corresponding to each model are given in Figure 1. The simulation experiments were performed at two different values of n (n = 25 and n = 40, where n is the number of subjects). The subjects in all simulations each had 20 observations, 5 of which were randomly selected (m_i = 5 for all i). With such a small sample size per subject, it was not feasible to smooth observations separately for each individual as has been done in previous work (Behseta et al., 2005).

For each simulated data set, we sampled from _{Plugged In} (k, t_k, k′, t_k′| y) and _Laplace(k, t_k, k′, t_k′|y) using Algorithm I, and p(k, t_k, k′, t_k′, ψ, σ²|y) using Algorithm II. We refer to the posterior sample from p(k, t_k, k′, t_k′, ψ, σ²|y) as the “truth” or the “true” posterior sample. The number of iterations for each chain of _{Plugged In}(k, t_k, k′, t_k′|y) and _Laplace(k, t_k, k′, t_k′|y) was at least 10,000, and for p(k, t_k, k′, t_k′, ψ, σ²|y) the chains were twice as long. All chains were run using R 2.5.0, and the average time it took to complete 1000 iterations at n = 25 was 18 minutes for _{Plugged In}(k, t_k, k′, t_k′|y), and 54 minutes for _Laplace(k, t_k, k′, t_k′|y) and p(k, t_k, k′, t_k′, ψ, σ²|y). At n = 40 the average time it took to complete 1000 iterations was 30 minutes for _{Plugged In}(k, t_k, k′, t_k′| y), and 78 minutes for _Laplace(k, t_k, k′, t_k′|y) and p(k, t_k, k′, t_k′, ψ, σ²|y). Convergence to a stationary distribution did not seem to be an issue with any of the chains run. The distributions of k and k′ corresponding to each approximation and p(k, t_k, k′, t_k′, ψ, σ²|y) are shown in Figures (2) – (5). The captions of each figure also give the average acceptance rates for the plugged-in, Laplace, and “true” chains. The acceptance rate of the Laplace chains are, in each case, higher than those for the “true” or plugged-in chains. Figure 6 gives the posterior distributions of t_k and t_k′ for data simulated from model 1 at n = 25. The accuracy of the Laplace and plugged-in as shown in Figure 6 is similar to the accuracy seen in data simulated from model 1 at n = 25 and from data simulated from model 2 at n = 25 and n = 40 (not shown).

The marginal distributions of k′ show that the Laplace method is more accurate than the plugged-in method. In each case, the marginal distributions of k′ for the Laplace method essentially matches the “true” marginal distributions of k′. The plugged-in estimator greatly overestimates the true number of random effect knots (see Figures 2 – 5). This is to be expected since both ψ and σ² were not averaged over (averaging over the parameters would penalize vectors of ψ with large dimensions). The marginal distributions of t_k′ in Figure 6 also show that the Laplace method is more accurate than the plugged-in method. The posterior distributions of t_k′ corresponding to _Laplace(k, t_k, k′, t_k′|y) and p(k, t_k, k′, t_k′, ψ, σ²|y) are essentially the same. The posterior distribution of t_k′ corresponding to _Plugged–In(k, t_k, k′, t_k′|y) is far from the true posterior distribution of t_k′; the plugged-in method fits far too many random effect knots and thus can not correctly identify their true locations. Thus, in terms of accuracy, the Laplace approximation is definitely preferred to the plugged-in approximation.

In addition to examining the two posterior distributions of knots, we also examine the two corresponding posterior distributions of one randomly selected cow’s curve. Recall that the curve of one cow is the sum of the fixed curve and that cow’s random effect curve. To sample one cow’s curve from the chain corresponding to _Laplace(k, t_k, k′, t_k′|y), values of ψ and σ² were sampled using an independence sampler. Conditioned on these sampled values, a fixed effect curve was sampled, and conditioned on this sampled fixed effect curve, a random effect curve was sampled. The algorithm used to do this requires a slight addition to Algorithm I. This addition is noted below: This addition to Algorithm I slows it down only marginally.

We are currently investigating other approaches to the knot selection problem. Rather than sampling from the posterior p (k, t_k, k′, t_k′|y), we set t_k = t_k′ and then try to reduce the number of principal component curves associated with the random effects by sampling from the posterior p (k, t_k, r|y) where r = the number of random effect principal components retained (Botts, 2005). This is similar in spirit to the methods employed by Shi et al. (1996), and James et al. (2000).