Effect of Smoothing in Generalized Linear Mixed Models on the Estimation of Covariance Parameters for Longitudinal Data

2.1 Parametric nonlinear function

By parametric nonlinear function we refer to a regression model whose mean is linear in parameters but has higher order of covariates. The simplest and most common way to represent curvature in regression models is using polynomials of the covariates, typically quadratics. Adding a parametric nonlinear term (e.g., quadratic), f(⋅) in (1) can be modeled as

One advantage of this parametric approach is that it is more efficient if the model is correct. However, conventional low order polynomials do not always fit the data very well. High order polynomials follow the data more closely but often fit poorly at extreme values of the covariates [14].

2.2 Fractional polynomials

Fractional polynomials (FP), proposed by Royston and Altman [14], are a family of curves whose power terms are restricted to a small predefined set of integer and non-integer values. For one covariate x, a FP of degree d with powers p1,…,pd is given by

where powers p1,…,pd are taken from S,

Here, x0 is defined as logex. Usually, d=1 or d=2 is sufficient for a good fit as models with degree higher than two are rarely required in practice [14]. Ambler and Royston [15] suggested an iterative algorithm for covariate selection and model fitting when several covariates are available. For d=2, their algorithm first selects the best-fitting second-degree FP which has the lowest deviance (−2× log likelihood among all possible models with two powers, (p1,p2)∈S. The best selected model is then tested against the null model, a straight line, and the best-fitting first-degree FP (with lowest deviance among all models with power p1∈S) in an attempt at further simplification. All significance tests are performed by calculating an approximate P-value based on the difference in deviances having a chi-squared or F distribution, depending on the regression in use. This approach has been implemented in most popular statistical software (e.g. the mfp package in R, macros in SAS (www.imbi.uni-freiburg.de/biom/mfp), and the fp function in Stata).

2.3 Semiparametric mixed models (SPMMs)

The nonlinear association between an outcome and covariates can be modeled using penalized splines within the framework of a mixed effects model. This approach is a trade-off between regression splines (which rely on the number and position of knots) and smoothing splines (which require intensive computation for larger datasets) [16]. To capture the structure of f(⋅) in (1), we define K distinct knots t1,…,tK in the range of x and consider the regression spline

where β0,…,βp are the regression coefficients, uk denotes the spline coefficient at knot tk, p is the degree of polynomial used (e.g., two for quadratic, three for cubic), and the truncated basis function component (x−tk)+p=max0,(x−tk)p. Note that when p=1 the knot coefficients uk represent changes in slope from one segment to the next. Unpenalized estimation of uk would lead to a bumpy fit due to the large number of truncated polynomials. To avoid overfitting we assume the uk are independently distributed as

This constraint shrinks the uk towards zero to reduce the magnitude of slope changes, leading to a smoother fit.

Denoting w=w11…wmnmT, x=x11…xmnmT, z=z1…zmT with zi=zi1…ziniT, the vector of fixed effects β=(γT,β0,…,βp)T, vector of random effects u=(u1,…,uK,b1,…,bm)T, the associated design matrices X=[w 1 x x2 … xp] and Z=(x−t11)+p…(x−tK1)+pblockdiagonal1≤i≤mzi, where 1 is the vector of all ones, and the response vector Y=Y11…YmnmT, we can rewrite (1), (4) and (5) as the generalized linear mixed model

where

This mixed model representation of penalized splines allows us to take full advantage of existing methods and software for GLMMs. The resulting models are called semiparametric mixed models (SPMMs) [6]. We discuss maximum likelihood (ML) estimation of the SPMMs in the subsequent subsection.

An important special case of (6) is the Gaussian mixed model with random intercepts bi∼N(0,σb2). Fitting such a SPMM is equivalent to minimizing the penalized least squares (PLS)

where the smoothing parameter λ=σε2/σu2. Here, σε2 is the variance of the random error, and for a general vector v, v≡vTv is the usual Euclidean norm of the vector v. One could take the approach of minimizing the PLS directly, but obtaining ML estimates from the SPMM is advantageous because the smoothing parameter, λ, is estimated directly from the data as λˆ=σˆε2/σˆu2, though it is possible to specify it directly if desired [4]. Moreover, using the SPMM approach to estimate the smoothing parameter works better (in terms of mean squared error performance and computational stability) than generalized cross validation (GCV), Akaike information criterion (AIC), etc. [17].

SPMMs require the number and position of the knots to be specified in advance in addition to specifying the degree (p) of polynomial for the curve segments. In practice, the simplest choice for p is 1 (see, for example, Gurrin et al. [4], Wand [5] and Durban et al. [16]) so that the truncated line basis for the knots t1,…,tK is

These bases are preferred because of their simple mathematical form, which is very useful when formulating complicated models [5]. Greater smoothness can be achieved using higher degree truncated polynomial bases such as truncated cubic (p=3) basis functions, however the use of truncated polynomial (p≥2) bases has been criticized by many authors because of their sub-optimal numerical properties [5] and poor behaviour in the tails [18]. One solution to this is to use natural cubic splines. More complex bases such as B-splines [19] or radial basis functions [6] that have better numerical properties could also be easily used to fit models under this framework.

2.4 Estimation

Once f is specified by any of the methods discussed above, model (1) can be written in a general form

where u∗∼ fu∗(u∗|G)=N(0,G). Let ρ be the vector of components in CorrYij,Yij′|ui∗. The likelihood function is then given by

Except in the case of normally distributed outcomes, the likelihood in (7) does not have a closed form as the integral is intractable. In these cases, numerical approximations are required. While various approximation procedures and Bayesian methods using Gibbs sampling and EM algorithm have been proposed in the literature (see, Demidenko [20] for details), the penalized quasi-likelihood (PQL) method of Breslow and Clayton [3] is the most frequently used estimation technique because of its simplicity and relatively small computational effort. Lin and Zhang [7] proposed the double penalized quasi-likelihood to make approximate inference in generalized additive mixed models using smoothing splines. For the linear mixed effects model, obtaining maximum likelihood estimates from (7) is quite straightforward. However, to obtain less biased estimates for the covariance parameters, restricted maximum likelihood estimation (REML) is often used.

3.1 Data generation

3.2 Analysis of simulated data

Each simulated dataset was analysed by fitting a generalized linear mixed model that included a random intercept and AR(1) serial correlation structure. In some scenarios we generated data with either the serial correlation (ρ) or the random effect variance (σb2) set to zero, but fit models that included both random intercept and AR(1) correlation. We call these “misspecified” although they are just special cases with extreme values of the model parameters. Model parameters were estimated via restricted maximum likelihood estimation (REML) for normally distributed data and the penalized quasi likelihood (PQL) approach of Breslow and Clayton [3] for Poisson and binary data. In each case, we adopted the following approaches to estimate the smooth curve, f(⋅):

1. Semiparametric mixed model (SPMM) using truncated line bases with number and position of knots specified by the simple and reasonable default rule given by Wand [5] as

where 1≤k≤K, and K=min(n/4,35). For the parameters of all the fixed effects and the smooth term, we obtained the approximate (frequentist) covariance matrix following the approach described in Lin and Zhang [7].

1. Fractional polynomials (FP) of degree 2 to identify the best curve by ignoring the correlated nature of the data. To select the best FP, we followed an algorithm proposed by Ambler and Royston [15] (as implemented in mfp package in R) that has been described in Section 2.2. We then fit the GLMM to accommodate the correlation structure of the data using the best selected curve from the FP fit.

2. Quadratic polynomials as the simplest parametric non-linear function.

3. Linear function.

As a sensitivity analysis we also fit a model using the true data-generating curve. For the linear, quadratic polynomial and FP, curve f(⋅) was estimated as a fixed effects component. The approximate standard error of the estimated curve fˆ(⋅) was obtained conditional on the estimates of the random-effects parameters and therefore this only accounted for the uncertainty of the fixed effects. All simulations and analyses were carried out in R software employing lme and glmmPQL functions for REML and PQL methods, respectively. For both functions default procedure values were used unless otherwise noted. We allowed 200 iterations for convergence in all models fit.

3.3 Performance indicators

For the covariance parameters estimators ρˆ, σˆb2, σˆε2, and fixed effect estimator βˆ, we computed percentage relative biases (PRBs), simulated mean squared errors (MSEs), and empirical coverage probabilities (CPs). PRB was computed as

while the empirical coverage probability (CP) of an estimator αˆ of α was obtained as

where αˆ∈{ρˆ,σˆb2,σˆε2,βˆ}, and αˆL and αˆU, respectively, are the lower and upper limits of the approximate confidence intervals for α. The approximate confidence intervals for the parameters were obtained using a normal approximation to the distribution of the (restricted) maximum likelihood estimators following Pinheiro and Bates [23].

For the smoothing curve estimators, fˆ, we computed pointwise mean average squared distance/error (MASE) from the true curve, and the mean average coverage probability (MACP). The pointwise MASE was defined as the mean over the 1,000 replicated datasets of the pointwise average squared error,

The 95% pointwise MACP was obtained as the mean of the 1,000 pointwise average coverage probabilities,

where SE(fˆ(xij)) is the standard error of fˆ(xij).

To compare the fits graphically, we plotted the mean fitted values and 95% pointwise coverage probabilities of the true function. At each observed value of x, the mean fitted value was obtained by taking the average over the 1,000 replications. Pointwise coverage probability was computed as the proportion of replications in which the confidence interval contained the true curve.

3.4 Simulation results

We analyzed the simulated data for 100 and 500 clusters with homogeneous (ni=5) and heterogeneous (1≤ni≤10) cluster sizes. Since the results showed consistent patterns, we report only the case of ni=5 for m=500. Results from one scenario with varying cluster sizes (1≤ni≤10) are also presented in the supplementary materials (Supplementary Table 4).

Simulation results for normally distributed data reported in Table 2 indicate that for complex functional forms (e.g., double hump), the SPMM method yielded very good estimates for the correlation, variance, and fixed effect parameters in terms of bias and MSE. It also performed very well in estimating the nonparametric function (with respect to MASE) and had the correct average coverage probability. On the contrary, FP, quadratic, and linear fits performed very poorly in estimating all parameters except for the fixed effect (β). For slightly nonlinear associations (e.g., linear quadratic threshold), the SPMM approach estimated all the parameters extremely well and produced very good coverage probability for the nonparametric function. The FP and quadratic fits performed similarly to SPMM although they had slightly lower coverage probabilities for nonparametric function estimates. Using a linear functional form, however, performed poorly. When the true association was linear, results from all the methods were very similar with almost unbiased estimates for all parameters and small MSEs. The mean coverage probabilities for the estimated function were exactly 95% for all methods. Similar results were observed when ρ=0.75,σb2=σε2=0.1 and ρ=0.25,σb2=σε2=0.1 (data not shown).

In the two slightly misspecified cases, that is, when the correlation was due to (i) only random intercepts but no serial correlation (ρ=0,σb2=σε2=0.1) or (ii) serial correlation but no random intercepts (ρ=0.75,σb2=0,σε2=0.1), the results were similar to those that have been reported above except that the estimates of σε2 were slightly biased in case (ii) for all smoothing methods and shapes (see supplementary Table 1).

Table 3 shows results when the responses were Poisson distributed. The SPMM approach worked well in estimating all parameters and had average coverage probability for f(x2) close to the nominal value when the true functional form was complex (e.g., double hump). FP, quadratic, and linear fits, on the other hand, performed badly in estimating the variance of the random effects and f(x2), although the auto-correlation and fixed effect (β) were reasonably estimated. For the quadratic threshold shape, SPMM, FP, and quadratic fits estimated all the parameters reasonably well. As expected, the linear fit estimated the random effects variance and f(x2) poorly. For the linear functional form, all methods performed similarly and produced good estimates for all parameters. Similar results were observed for all other scenarios investigated: ρ=0.75,σb2=0.25; ρ=0.25,σb2=0.25; ρ=0.75,σb2=0; ρ=0,σb2=0.25 (see supplementary Table 2 for the results from last two scenarios).

Table 4 presents the simulation results for binary response data. Results suggest that for the double hump shape, the SPMM technique outperformed all other methods and estimated all the parameters fairly well. For the linear quadratic threshold and linear shapes of association, the SPMM, FP, and quadratic approaches estimated all parameters well except that the σb2 estimates were positively biased. As ρ increased, the bias in estimating σb2 also increased for all shapes of association (see supplementary Table 3). The estimates of σb2 were found to be biased for all methods in the slightly misspecified cases, i.e., when ρ=0,σb2=0.25 and ρ=0.5,σb2=0 (supplementary Table 3). While the results are not shown, we note that biases similar to those given by the SPMM method existed even when the true data-generating model was used for estimation.

The estimated (model based) standard errors (SEs) of the parameter estimates (βˆ,ρˆ,σˆb2 and σˆε2) agreed well with the empirical (simulation based) standard errors for all methods in all cases except for σˆb2 in the binary data, where empirical SEs were 3% to 28% larger than estimated SEs for all methods (results not shown).

The coverage probability (CP) was generally near nominal level for β for all methods in all cases. For ρ, σb2 and σε2, the CPs obtained from the SPMM approach were close to the nominal level for all shapes of association considered while for other methods CPs varied depending of the true shape of the association (results not shown).

The left panel of Figure 1 illustrates the ability of the four estimation methods to recover the true functions by comparing the four estimated curves based on 1,000 replications for normally distributed responses. The SPMM reconstructed the true curves extremely well for all the shapes. FP and quadratic fits performed well in capturing the true curves for linear quadratic threshold and linear shapes, whereas both the methods yielded unsatisfactory fits for the double hump function.

The right panel of Figure 1 compares the empirical pointwise coverage probabilities of the 95% confidence intervals of three functions calculated using four different methods. For normal data, the coverage probabilities of all confidence intervals were very close to 95% at well estimated values of x2. For the double hump curve, the SPMM gave undercoverage for small values of x2≤0.1, i.e., the region where the signal-to-noise ratio was low. Afterwards, the coverage probability remained close to 95%. For quadratic threshold shape, the coverage probability of the SPMM confidence interval agreed slightly better with the nominal value throughout the range of x2 and performed better than other methods. When the true functional form was linear, confidence intervals from all methods provided correct coverage probabilities although SPMM had slightly inferior performance.

Similar results in estimating nonparametric functions were obtained for Poisson and binary responses. However, the SPMM estimated curves were slightly biased and had slightly lower coverage probabilities (see supplementary Figures 1 and 2).

Figure 1:

True and estimated nonparametric functions based on the mean fitted value at each observed value of x2 (left) and estimated 95% pointwise coverage probabilities of the true functions (right) for each of the four fits. These results are for normally distributed data with m=500, ni=5, ρ=0.5,σb2=σε2=0.1 and from 1,000 replications.

Figure 2:

CD4 cell counts with parametric, fractional polynomial, and semiparametric mixed model estimates for the mean trend time obtain from (a) LMM fits, (b) Poisson GLMM fits and (c) Logistic GLMM fits.

Convergence or singularity problems for the SPMM method were relatively minor (<0.5%) and only occurred for the case when the true functional form was linear. Results from cases that did not converge were omitted.

4.1 Data and variables

The Multicenter AIDS Cohort study (MACS) followed 369 HIV infected men aged between 23 and 64 in the USA in the mid 1980s (see, Kaslow et al. [9], and Zeger and Diggle [10] for details). Repeated measurements were obtained for each subject on CD4 cell counts, an important immunological marker which measures the body’s ability to fight off infections [24]. Time, measured in years since the date of seroconversion, was known approximately for each subject from several years before seroconversion and up to 6 years after seroconversion. Observations were taken every six months on average for a total of 2,376 measurements (on average 6.5 measures per patient, varying between 1 and 12). Several other measures such as age at study entry, number of cigarette packs smoked per day, number of sexual partners, a measure of depression status (CESD > 0 indicates possible depression), etc. were also collected. This dataset has been analyzed previously (e.g., by Diggle et al. [25]).

4.2 Data analysis

The relationship between time since seroconversion and CD4 count is not immediately apparent from visual inspection of the scatter plot (see Figure 2(a)), apart from a sparsity of subjects with higher CD4 count after 2 years of seroconversion, suggesting that CD4 count decreases with time after seroconversion. With the goal of identifying the progression of mean CD4 counts as a function of time since seroconversion, we used smoothing techniques within a linear mixed effects model (LMM), considering the square root of CD4 count as continuous outcome. We also considered CD4 count as Poisson outcome and as a binary measure, where Y=1 indicated low CD4 count, i.e., CD4 count ≤ 350 cells/mm³ and Y=0 indicates CD4 count > 350 cells/mm3. Note that as of March, 2015, the guidelines from the U.S. Department of Health and Human Services suggest starting HIV treatment when one’s CD4 count falls to 350 cells/mm³ or below (www.aids.gov).

To adjust for the effects of potential confounders, smoking (number of packs per day), recreational drug use (yes/no), number of sexual partners, and depression symptoms (as measured by the CESD scale) were also included in all the models as covariates. We centered the number of packs smoked per day, and number of sexual partners at their mean values and included them as linear covariates.

Considering the square root of the CD4 count as a continuous response, we constructed the empirical correlation matrix and trajectory plot. Based on these, we adopted a random intercept model with an autoregressive AR(1) correlation structure for the errors. We fit the model

where f is some smooth function of time, i=1,…,369, j=1,…,ni, 1≤ni≤12, βpack,…,βcesd are regression parameters, and bi∼N(0,σb2).

Assuming the square root of the CD4 cell counts follow a normal distribution, we considered gEYij|bi=ECD4ij|bi and assumed the random error εij∼N(0,σε2V(ρ)) such that the correlation between two observations measured at time tij and tik for ith subject was Corr(Yij,Yik)=ρ|tij−tik| for j≠k=1,…ni.

Representing CD4 count as Poisson distributed outcome, we fit the GLMM (14) with gEYij|bi=logECD4ij|bi and considering AR(1) serial correlation for responses. For binary GLMM, we fit model (14), where gEYij|bi=logitPCD4ij≤350|bi and responses were AR(1) correlated. In all cases, we considered the continuous specification of the AR(1) correlation structure that took into account the actual distance between time points.

For each data distribution, we considered four approaches to estimate the nonlinear effect, f of time on CD4 count.

1. A parametric approach, as suggested by Diggle et al. [25]. Following their approach, we specified the mean time-trend for CD4 count as constant prior to seroconversion and quadratic in time thereafter.

2. Using second degree fractional polynomials, where we first identify the best curve within the fractional polynomials framework by ignoring the correlated nature of the data, and then fit the GLMM to respect the correlation structure of the data using the best selected curve.

3. Applying the semiparametric mixed model to smooth the regression spline implementation of the time-CD4 association. We considered linear splines with knots specified as in (13). We also computed the degrees of freedom (df) of the smooth fit following Ruppert et al. [6] (Sections 8.3 and 11.4) to quantify the amount of smoothing used for estimating f. A linear function uses 1 df whereas a nonlinear function uses some number greater than 1, depending on the bumpiness of the function.

4. A linear function.

To estimate the model parameters we used the marginal likelihood approach where estimates were obtained by REML for LMM. The parameters of the GLMMs were estimated via PQL method. All analyses were performed in R.

4.3 Results

Table 5 shows the summaries of the baseline characteristics for 369 included men. Mean and standard deviation (SD) are reported for quantitative variables whereas for categorical variables counts and percentages are provided. The mean patient age was 37.4 years.

Results from all models fits are presented in Table 6. The first panel of Table 6 summarizes results obtained from the mixed effects models fit considering square root of CD4 cell count as a continuous outcome and adopting four different smoothing approaches. The variance and correlation parameter estimates were very similar for all the smoothing methods. The estimated within-subject AR(1) correlation was small (around 0.15). The variance of the random intercept σb2 was estimated as high as around 15 for each of the considered models. The estimated degrees of freedom (df) of the smooth fit was 11.8. The fitted CD4-time trends are displayed in 2 (a) at the mean value of the confounders. The SPMM fit shows that the mean CD4 cell counts remained relatively stable until the time of seroconversion but decreased sharply over the first year after seroconversion. The counts then decreased gradually for the remaining observed time. The FP and parametric (following Diggle et al. [25]) fits did not closely reproduce the behaviour of the SPMM. Indeed, for these two fits there was an increase in CD4 counts after 4 years of the seroconversion. This may be a consequence of the parametric forms of the regression functions or due to a treatment or survivor effect.

The second panel of Table 6 shows results obtained from the models fit considering CD4 cell count as a Poisson outcome. The estimates of the correlation and variance parameters were very similar for all approaches. Fitted curves are shown in Figure 2(b) which are visually indistinguishable from those obtained from models considering CD4 cell counts as continuous.

Results from model fits considering CD4 cell count as a binary response with P(CD4 cell ≤350)=1 are shown in the third panel of Table 6. The estimates of the between-subject variance component σb2 were similar for all the models and quite large on the logit scale for binary outcomes. The estimates of ρ from all the models were near zero indicating no evidence of AR(1) correlation between observations within a subject. The estimated degrees of freedom was 3.13. Figure 2(c) displays the fitted probability of CD4 cell count ≤350, considering three smooth curves. For the SPMM fit, it is apparent that the risk of having a CD4 cell count ≤350 was zero up until the time of seroconversion, and it increased over time after seroconversion. The models fit using FP and parametric (following Diggle et al. [25]) largely reproduced the trend of SPMM except that they showed a slight decrease in risk after five years of seroconversion.

Simulation results suggested that the correlation and variance parameters were not well-estimated from the linear fit in the presence of curvature. Nevertheless, in our data analysis, the linear fit yielded similar estimates for the correlation and variance parameters to other smoothing approaches. This may be attributed to the fact that the CD4-time association was only slightly nonlinear. However, the shape of the association was better captured by the smoothing methods.