Flexible Regression Models for Rate Differences, Risk Differences and Relative Risks

3.1 Definition and properties

The B-splines are a series of basis functions for polynomial splines, defined by a grid of qc fixed turning points ξc1<⋯<ξcqc, where ξc1=wc(0) and ξcqc=wc(1) are positioned at the extremes of the domain. The grid is expanded to form a sequence τc of knots that determine the degree and continuity of the resulting curve [12].

We will restrict our focus to the case of B-splines of order 3, where each basis function is made up of a series of quadratic curves between each pair of turning points, with adjacent curves constrained to have equal gradient at their boundaries. The sequence of qc+4 knots is defined such that there are three knots at the lower boundary, three knots at the upper boundary and one at each of the internal turning points. That is,

Given a knot sequence, the Dc=qc+1 B-splines can be calculated recursively [10]. We use Bcd(w) to represent the d th basis function Bd(w;τc) on knot sequence τc, and Bc=B(τc) to denote the associated function space for fˆc.

The B-splines are normalized such that

for all w, meaning that a constraint must be applied to the parameters for each c to ensure that they are identifiable. We do this by setting γctc=0 for some choice of tc∈{1,…,Dc}.

Each basis function Bcd(w) is positive for all w∈(τcd,τc(d+3)) and zero elsewhere, taking its maximum value for some w∈(τc(d+1),τc(d+2)). This means that each parameter γcd has only local influence on the smooth function fˆc.

As an illustration we present a graphical representation in Figure 1. The four B-spline basis functions of order 3 with a single internal turning point, and hence 7 knots, are shown in Figure 1(A). Figure 1(C) shows the curve fˆ that results from taking a particular linear combination of these basis functions:

which can be equivalently expressed as

by using eq. (3). Figure 1(B) and (D) illustrates the use of B-splines to restrict fˆ to be monotonic, and will be discussed in Section 3.3.

Figure 1:

(A) B-spline basis functions and (C) example resulting curve with one internal knot, and (B) the equivalent monotonic B-spline basis functions and (D) resulting curve based on the same coefficients.

3.2 Model constraints

The models considered in this paper have constraints on their parameter spaces due to the restricted range of the response variable: rates must be non-negative and risks must lie in [0,1]. These constraints are often the source of instability in standard fitting algorithms, particularly when the MLE is close to the boundary of the parameter space [2].

The CEM algorithm is ideal for fitting these models, as it applies the parameter space constraints while guaranteeing stable convergence to the MLE. A crucial step in the definition of a CEM algorithm is that the parameter space is partitioned into a sequence of subspaces, each of which corresponds to a particular set of constraints on individual parameters. The properties of the B-splines as discussed in Section 3.1 allow us to extend the existing methods to include these semi-parametric functions.

As an example, we restrict our attention to fˆc∈Bc+, the space of strictly non-negative curves in Bc. This can be done without loss of generality due to eq. (3), such that the range of fˆc is determined by the intercept α0, and its shape by the B-spline coefficients.

The function space Bc+ can be partitioned into subspaces defined by the index of the smallest coefficient. That is, if we define

it is easy to see that

For a particular identifiability constraint γctc=0, if the remaining coefficients are restricted to be non-negative, the resulting fˆc will be a strictly non-negative curve. Furthermore, γctc will be the smallest of the coefficients, that is, fˆc∈Bc+(tc).

One characterization of the constrained function space is that any curve fˆc(w)∈Bc+(tc) will have a local minimum for w∈[τc(tc+1),τc(tc+2)]. The special case of tc=1 corresponds to the family of non-negative curves that take their minimum at the lower limit of the domain, wc(0), and likewise tc=Dc will constrain the curves to take their minimum at wc(1). This is demonstrated with an example in Section 4.2.

An analogous result can be achieved by fixing γctc=0 and restricting the remaining parameters to be non-positive. Then fˆc(w) will be a non-positive curve with γctc being the largest coefficient, a function space we denote by Bc−(tc). Similarly to the non-negative case discussed above, Bc−(tc) can be characterized by its restriction on the location of local maxima.

A useful property of the existing CEM algorithms is that they can easily accommodate non-negativity or non-positivity constraints on individual coefficients. By applying such constraints, we can find the MLE of the GAM under the restriction that each fˆc∈Bc+(tc) or Bc−(tc), and by repeating this process for all possible choices of the identifiability constraint, we can find the overall MLE. This process is explained in more detail for the specific models of interest in Section 4.

3.3 Monotonic B-splines

In some contexts, it may be sensible to restrict fˆc(w) to increase with increasing w. A sufficient condition for fˆc to be monotonically non-decreasing is that the coefficients of the B-splines are themselves strictly non-decreasing [13], that is, γc1≤⋯≤γcDc.

Assuming non-negative coefficients, the identifiability constraint γc1=0 ensures that fˆc(wc(0))=0 is the minimum of the smooth function. We can then introduce Dc−1 new parameters that represent the increments between successive B-spline coefficients:

for d=2,…,Dc. Restricting these increments to be non-negative will then ensure that the original coefficients are monotonically non-decreasing, as desired.

This can be achieved in the same manner as the unrestricted case by re-expressing eq. (2) as

Here, Bcdm+, d=2,…,Dc, denotes a new series of non-negative monotonic basis functions, which we will refer to as monotonic B-splines. In fact, for the order 3 B-splines used here, these new basis functions are equivalent to the piecewise quadratic integrated splines (I-splines) introduced by Ramsay [12] and used by Tutz and Leitenstorfer [14] for generalized smooth monotonic regression.

Continuing with the illustration provided in Section 3.1, the monotonic B-spline basis functions associated with the B-spline bases in Figure 1(A) are shown in Figure 1(B). Because the B-spline coefficients in eq. (4) are monotonically non-decreasing, the resulting curve in Figure 1(C) can be expressed in terms of the monotonic B-splines

and this is demonstrated in Figure 1(D), where the intercept term is shown as a horizontal line.

For the alternate case in which we wish for the curve fˆc to be non-positive and monotonically non-decreasing with its maximum value at wc(1), a similar process applies. We fix γcDc=0, and define Dc−1 new parameters, δcd=γcd−γc(d+1) for d=1,…,Dc−1. The associated monotonic basis functions are defined as

and their coefficients δcd are constrained to be non-positive.

3.4 Knot selection

The number and placement of the turning points will influence the shape of the resulting function by determining the space in which our estimated fˆc is constrained to lie. Ideally, external information would be used to determine the turning points; however, in many applications this will not be available and we are forced to depend on our data to guide this decision.

In some situations, such as OLS, the positioning of the turning points is crucially important [15], and a large number of knot selection methods have been derived (e.g. [10], pp. 174–196, [16]). Many of these methods can be integrated into the approach we describe, but they often lead to a large increase in the computational burden. In a general setting, Ramsay [12] noted that the shape of the resulting estimate is not particularly sensitive to knot placement.

We apply an adaptation of cardinal splines discussed by Hastie and Tibshirani ([7], p. 24), placing the qc−2 internal turning points at evenly spaced quantiles of the observed covariate values wic. This is often the standard approach used with regression splines (e.g. [17]). Alternatively, our implementation in R discussed in Section 6 optionally allows the user to specify their own list of knots, so other knot-placement regimes could be used, such as equally spaced knots, nested knot structures or the approach proposed by Yao and Lee [18], in which knots are placed at the local minima and maxima of an initial estimate of the smooth curve, fitted using a basic knot structure.

Of more importance is the choice of the number of turning points. With too few, we may fail to detect important features of the relationship, but with too many, we are at risk of over-fitting. One common way to resolve this trade-off is to use a sufficiently large number of turning points to broaden the function space for fˆc, but add a penalty term to the likelihood function such that spurious fluctuations in the smooth function are avoided [19].

However, with a penalty term, the CEM algorithm central to these stable methods cannot be used directly. We discuss this further, and propose a possible solution in Section 7. In general, we will use the Akaike information criterion AIC=2K−2ℓ [20] to choose between models with different numbers of knots, where ℓ is the log-likelihood of the fitted model. This similarly includes a penalty for model complexity, and in the case of unpenalized maximum likelihood estimation, the effective degrees of freedom K is simply the number of estimated parameters in the model [8]. The AICc is a bias-corrected version of the AIC for small samples [21], which becomes virtually identical as n increases. Figure 4 illustrates the use of the AIC in determining the optimal number of knots for the example in Section 5.

As a criterion for choosing the optimal smoothing parameter in a wide range of scenarios, the AIC and AICc have been shown to be generally superior to other classical approaches such as BIC and GCV in the context of general non-parametric regression by Hurvich, Simonoff and Tsai [22], semi-parametric Cox regression by Malloy, Spiegelman and Eisen [23] and likelihood-based boosting by Leitenstorfer and Tutz [13].

Figure 2 shows the results from a simulation study in which we examined the performance of the AIC value in selecting the optimal number of knots. We simulated 500 datasets of 750 binomial observations, each with true risk functions that were based on B-splines with 1, 2 or 3 internal knots. For each dataset, we found the MLE corresponding to a binomial model in which the continuous variable was included as a linear term or as a B-spline with between 0 and 5 internal knots. Of these, the model with the smallest AIC was selected as the optimal model. This approach selected the correct number of knots in the vast majority of samples, and the mean number of knots selected was close to the true value. As has been observed in other contexts, the BIC was biased towards a low number of knots (oversmoothing), and the GCV criterion tended to undersmooth by choosing a higher number of knots more often (data not shown).

Figure 2:

Results from 500 simulations of the performance of the AIC for selecting the optimal number of knots. The panels on the left-hand side show the true risk function p(w) with (A) 1, (B) 2 and (C) 3 internal knots. The panels on the right-hand side show a histogram of the model selected by AIC in each simulation, with the bar corresponding to the true model shaded.

The computational effort required by full cross-validation renders it infeasible in this situation. Nevertheless, the focus of this paper is on providing a method of estimation for a single model, which may be applied within any scheme for determining the optimal number of knots.

4.1 CEM algorithm

Each of the methods that we use to fit the models in the subsequent sections is an application of CEM algorithms [6]. A CEM algorithm is a general approach in which we consider a finite family of complete data models, indexed by t∈T, each of which has a parameter space Θ(t) that is a subset of the parameter space Θ for the model of interest, such that

The complete data models are defined such that an expectation–maximization (EM) algorithm [24] can be used to find the constrained maximum of the likelihood θˆ(t), within each Θ(t). Then, due to eq. (8), the θˆ(t) that attains the highest likelihood is the MLE θˆ.

When the model of interest is a GLM, the complete data models typically impose some constraint on the individual model parameters, such as non-negativity. In the case of GAMs as described below, we use augmented complete data models, where the effect of such constraints on the spline coefficients is to constrain the shape of the smooth curves fˆc. The overall MLE will then correspond to one of these constrained maxima. We now describe the details for each of the specific models of interest.

4.2 Adjusted rate differences

Adjusted rate differences can be estimated by fitting an identity link Poisson GAM. Specifically, if Ni is the period of time over which Yi events were observed, and we assume Yi∼Poisson(Niλi), then λi is the event rate for an individual with covariate vector (ui,vi,wi).

With link function g(λ)=λ, the absolute rate difference is simply the difference between two linear functions. Keeping all other covariates constant, the adjusted rate difference associated with changing the ath categorical covariate from level uia to uja is αa(uja)−αa(uia). Likewise, the adjusted rate difference for a 1-unit increase in the bth linear covariate is βb.

The Poisson means must be non-negative, and so the parameter space within which the MLE θˆ must lie is

where U, V and W denote the covariate spaces defined by the Cartesian products of the observed ranges of covariate values:

For C=0, Marschner [3] has described a stable CEM algorithm to find the MLE θˆ∈Θ. This is achieved by partitioning the parameter space into a sequence of distinct constrained parameter spaces Θ′(r,s), in which (r,s)=(r1,…,rA,s1,…,sB) is the covariate vector associated with the minimum fitted Poisson mean.

For a particular choice of r and s, a complete data model is defined, consisting of A+B+1 independent latent random variables underlying each observed random variable Yi. These latent variables have Poisson distributions with means that each depend on just one of the model parameters. Using this separation of parameters, an EM algorithm is then defined to find the constrained MLE θˆ(r,s)∈Θ′(r,s).

By considering all possible reference vectors (r,s) that could result in the minimum fitted mean, and finding the MLE within each corresponding constrained parameter space, we can be sure that one of these constrained MLEs will be the overall MLE θˆ∈Θ.

For C>0, it is straightforward to extend this idea to handle smooth semi-parametric components. The reference vector is expanded to include t=(t1,…,tC), where each tc∈{1,…,Dc}. For a particular choice of t, the identifiability constraint γctc=0 is applied for each c=1,…,C.

The complete data model is augmented with ∑c(Dc−1) latent Poisson random variables Yi(cd), each having an expected value that depends on one of the remaining parameters, that is,

If the basis function values Bcd(wic) are viewed as additional continuous covariates, this augmented complete data model has the same form as that used by Marschner [3]. Thus the EM algorithm can be applied directly to find its MLE θˆ(r,s,t), noting that due to eq. (10), each γcd is constrained to be non-negative.

As discussed in Section 3.2, the non-negativity constraints on the γcd force each estimated fˆc to belong to Bc+(tc) as defined in eq. (5). The constrained MLE will then belong to the parameter space

where

Let R denote the set of all possible choices for r, S the set of possible choices for s and T the set of possible choices for t. Then due to eq. (6), we have that

We have thus defined a CEM algorithm for finding the overall MLE θˆ∈Θ, which will be the constrained estimate θˆ(r,s,t) associated with the highest likelihood.

This approach can be implemented directly using the existing CEM algorithm for identity link Poisson GLMs. Using this algorithm, it is straightforward to apply non-negativity constraints to the coefficients of some continuous covariates by considering only one of the two possible reference levels for that covariate.

For a particular choice of t, if we apply the CEM algorithm to the categorical and linear covariates as usual, and include the basis function values Bcd(wic) (for d≠tc) as linear covariates with non-negativity constraints, we will find a constrained MLE θˆ(t)∈Θ(t), where

Repeating this for all t∈T will result in a collection of constrained MLEs θˆ(t), one of which will be the overall MLE θˆ∈Θ. There are a total of ∏c=1CDc elements in T, so by applying the CEM algorithm for each, the maximum total number of applications of the EM algorithm required to find the MLE is

The process can be halted if one such application converges to a point in the interior of the constrained parameter space, as we can be sure that this is the overall MLE.

The process of cycling through the parameter space partition is illustrated in Figure 3 for a Poisson model with C=1 covariate. Here we have simulated 500 observations with covariate values wi and means λi=50(α0+fˆ(wi)), based on a B-spline with two internal knots

where the usual dependence on c has been suppressed since C=1. By virtue of the constraint (3), one of the γd parameters must be set to zero while all others must be non-negative. The specific B-spline used to simulate the data for Figure 3 is

which has γ3=0 and is represented by the red dotted line in each plot. In this example, cycling through the parameter space partition involves computing the constrained MLE for each of the five models corresponding to the constraint γd=0, d=1,…,5. Thus, we would expect the constrained MLE with γ3=0 to yield the overall MLE. This is shown in Figure 3(A), while the constrained MLEs obtained by using the other constraints are shown in Figure 3(B)–(E). As expected, the log-likelihood for the curve with γ3=0 is higher than that for the other constrained MLEs (−1378 compared to −1437 through −1548) so it is the overall MLE.

Figure 3:

Illustration of the effect of parameter constraints on the smooth curve in a simulated Poisson dataset. The red dotted line represents the true underlying rate and the grey dotted lines denote the knot locations. The black solid line in each panel is the MLE under the identifiability constraint shown in the heading, with all other parameters constrained to be non-negative. The MLE is obtained with the constraint γ3=0.

One point to note is that the non-negativity of the coefficients γcd is a sufficient, but not necessary, condition for the non-negativity of the estimated fˆc. Thus, the parameter space over which we search for the MLE is not guaranteed to include all fˆc∈Bc+. However, in practice it would be highly unlikely that the true relationship we are modelling is exactly a quadratic spline on the selected knot sequence τc, and so if we find a MLE that is restricted by this deficiency, we can modify τc in order to achieve a better fit. In testing, we found that simply including an additional knot close to the minimum of the curve will generally resolve this issue.

4.3 Adjusted relative risks

With binomial data we have Yi∼Bin(Ni,λi), where Ni is the number of independent trials over which Yi binary events were observed, and λi is the constant event probability, or risk, at each trial. The adjusted relative risk is the ratio of event probabilities associated with a change in one covariate, keeping the others constant.

With g(λ)=log(λ), the risk is λi=exp{Λ(ui,vi,wi;θ)}, and so the relative risk is a ratio of two exponential functions. Keeping the other covariates constant, the adjusted relative risk associated with a change in the ath categorical covariate from uia to uja is then exp{αa(uja)−αa(uia)}. Similarly, the adjusted relative risk for a 1-unit increase in the bth linear covariate is exp(βb).

The fitted risks must lie within [0,1], and so the fitted linear predictors must be strictly non-positive. That is, the parameter space for this model is

where U, V and W are as defined in eq. (9).

With C=0, Marschner and Gillett [5] provide a CEM algorithm for finding the MLE θ∈Θ. The approach is similar to that used for the identity link Poisson GLM: for a given reference vector (r,s), the model is reparameterized and a complete data model is defined, which imposes non-positivity constraints on each of the transformed parameters. An EM algorithm is used to find the MLE θˆ(r,s) of the complete data model, which is constrained to lie in the parameter space

That is, (r,s) is the covariate vector associated with the maximum linear predictor, which is constrained to be non-positive.

The union of these constrained parameter spaces over all possible choices of (r,s) is Θ, and so the constrained MLE with the highest likelihood will be the overall MLE. It is straightforward to apply non-positivity constraints to individual parameters.

Extension to include smooth terms follows an analogous approach to that used in Section 4.2. For a particular t=(t1,…,tC), we set γctc=0 and apply the CEM algorithm, including the remaining basis function values Bcd(wic) as continuous covariates with non-positivity constraints. Then each fˆc∈Bc−(tc), the space of non-positive B-splines that have their shape constrained by the choice of γctc as the largest coefficient. The resulting estimate θˆ(t) is the constrained MLE over the parameter space

By considering each possible t∈T, and hence all possible locations for the maximum of each fˆc, we find a collection of constrained MLEs θˆ(t), and the one with the highest likelihood is then the overall MLE θˆ∈Θ.

Again, the sequence may stop early if a maximum in the interior of the parameter space is identified. Furthermore, the same caution also applies here as it did for the identity link Poisson model: non-positivity of the coefficients is only a sufficient condition for non-positivity of the linear predictor, so we are not searching the entire parameter space where fˆc∈Bc−. As before, this may be remedied by adjusting the choice of τc if necessary.

4.4 Adjusted risk differences

Adjusted risk differences can be estimated using a binomial GAM with identity link, that is, Yi∼Bin(Ni,λi) and g(λ)=λ. The adjusted risk difference associated with a change in the ath categorical covariate from uia to uja is αa(uja)−αa(uia). Similarly, βb represents the adjusted risk difference for a 1-unit increase in the bth linear covariate.

Again we require that the probabilities λ lie within [0,1], but now the parameter space Θ simultaneously imposes both lower and upper boundaries on the linear predictors, that is,

With C=0, Donoghoe and Marschner [4] have presented an approach for finding the MLE θˆ∈Θ. It exploits the multinomial–Poisson transformation [25] to convert the model into an equivalent identity link Poisson fit, and the CEM algorithm of Marschner [3] can then be applied to the transformed data to find the MLE. As with the other CEM algorithms, non-negativity constraints can be applied to some of the parameters, which is achieved by imposing such constraints when employing the identity link Poisson CEM algorithm.

However, the inclusion of the semi-parametric terms cannot proceed in exactly the same way as in Sections 4.2 and 4.3. The covariate space considered for continuous covariates entered into the identity link binomial CEM algorithm is the Cartesian product of the observed ranges of those covariates, that is, V in eq. (9). So if we include the observed B-spline values Bcd(wic) as continuous covariates, the covariate space will include vectors in which more than three of the basis functions for a given c are non-zero. By definition, this cannot correspond to the B-spline values for any w. Since the parameter space constraints are imposed for all points in the covariate space, using a larger covariate space than is necessary means that the parameter space is overly restrictive, and may not include the overall MLE.

Instead we must use a slightly different approach, based on the ordering of the B-spline coefficients. We begin by choosing t=(t1,…,tC), where each tc=(tc1,…,tcDc) is now a vector containing some permutation of {1,…,Dc}.

Recall that if we define new basis functions Bcdm+ as described in Section 3.3, non-negativity constraints on the coefficients δcd will impose monotonicity on the coefficients of the original basis functions, that is, γc1≤⋯≤γcDc. Similarly, for a particular permutation tc, we can define

and non-negativity constraints on the associated coefficients will impose an order restriction on the original coefficients, that is, 0=γctc1≤⋯≤γctcDc.

We denote by Bc(tc)+ the subspace of Bc+ in which the coefficients are ordered in this way, and note that the overall MLE must correspond to one such permutation. If for a particular t we enter the new basis function values Bcd(tc)(wic) as continuous covariates into the identity link binomial CEM algorithm, and impose non-negativity constraints on their coefficients, the resulting estimate θˆ(t) will be the MLE for the parameter space

There are Dc! possible choices for each tc, and hence ∏c=1CDc! possible choices for the vector t. In order to find the overall MLE θˆ∈Θ, we find the constrained MLE θˆ(t) for each possible choice of t, and choose the one with the highest likelihood.

As with the other models, we may stop early if we find a constrained MLE in the interior of the parameter space, and there are various strategies to identify the parameterizations that are more likely to contain this MLE in order to potentially reduce computing time [6].

4.5 Monotonic smooth regression

Imposing a monotonicity restriction on one or more of the smooth curves is straightforward for all of the models examined so far. As discussed in Section 3.3, a sufficient condition for fˆc to be monotonically non-decreasing is that the B-spline coefficients themselves are monotonically non-decreasing, that is, γc1≤⋯≤γcDc. Because only one possible choice for the value tc or the vector tc needs to be considered, adding a smooth monotonic covariate to a model requires no additional applications of the EM algorithm.

For the identity link binomial model in Section 4.4, applying this constraint is trivial. In iterating through the possible choices of t, we need only to consider tc=(1,…,Dc) rather than all possible permutations, giving 0=γc1≤⋯≤γcDc.

For the other models, we employ the methods discussed in Section 3.3. In an identity link Poisson model, we replace the B-spline function values Bcd(wic) by their monotonic B-spline equivalents, Bcdm+(wic), and enter these into the CEM algorithm as continuous covariates with non-negativity constraints on the associated coefficients δcd. The resulting estimate will have 0=γc1≤⋯≤γcDc, as desired.

For the log link binomial model we proceed similarly, replacing Bcd(wic) by Bcdm−(wic) as defined in eq. (7). These are entered into the CEM algorithm as continuous covariates with non-positivity constraints, ensuring γc1≤⋯≤γcDc=0.

It is important to note that the monotonicity of the coefficients is only a sufficient condition for the monotonicity of the resulting smooth function. Hence the function space over which we search for the MLE fˆc is only a subset of the space of monotonic functions on a B-spline basis. However, when Tutz and Leitenstorfer [14] used boosting techniques to impose constraints on the monotonicity of the function rather than only the coefficients, they found “much higher computational costs without much effect on performance”.