fmixed : A SAS Macro for Smoothing-Spline-Based Functional Mixed Eﬀects Models

In this article we implement the smoothing-spline-based functional mixed eﬀects models (Guo 2002) by a SAS macro by exploiting the connection between mixed eﬀects models and smoothing splines. The macro can handle ﬂexible design matrices and is easy to use. Input parameters and output results are described and explained. A numeric example and a real data example are used for illustration.


Introduction
Functional mixed effects models (see Guo 2002;Morris and Carroll 2006, among others) are extremely powerful tools in modeling longitudinal profiles. Among them, the smoothingspline-based functional mixed effects model (SSFMM; Guo 2002) has been gaining attentions in recent years. For example Stacey et al. (2009) used SSFMM to fit the subject specific curves of cytokines in human immunodeficiency virus infections, and Fagiolini et al. (2009) used SSFMM to evaluate the temporal differences of treatment effects in bipolar disorder patients. However, the lack of a corresponding software package still limits the general accessibility of SSFMM. In this article, we implement SSFMM in SAS and demonstrate how to use the resultant macro fmixed by a numeric example and a data example.
In literature, there are several approaches in functional mixed effects modeling, and some of them have public-accessible softwares implemented. Rice and Wu (2001) used B-spline basis and the mixed effects modeling is on the B-spline coefficients. James et al. (2000) used a combination of B-spline and functional principal components to reduce the number of parameters required for the mixed effect modeling of the B-spline coefficients. Guo (2002) modeled both the functional fixed effects and random effects as smoothing splines, which only requires limited number of parameters. Wu and Zhang (2002) used local polynomial, the mixed ef-fects modeling is on the polynomial coefficients and this method has been implemented in the companion MATLAB softwares to Wu and Zhang (2006). Morris and Carroll (2006) used wavelets which can handle abrupt and spiky signals, Bayesian approach was adopted for the mixed effects modeling on the wavelet coefficients and free software is available as indicated in their paper. Di et al. (2009) and Aston et al. (2010) used functional principal components as the basis and this approach has been implemented in Crainiceanu and Goldsmith (2010). The goal of this article is to implement Guo (2002) by a SAS macro, which adds one alternative among the large selection of different tools. A thorough comparison of those different approaches in functional mixed effects modeling, however, is beyond the scope of this article.
The key of the implementation is the equivalence between a smoothing spline and its Bayesian model (Wahba 1978). By utilizing this equivalence, Wang (1998) showed that a smoothing spline can be represented as a mixed effects model and consequently it can be implemented by SAS PROC MIXED. He also showed that the restricted maximum likelihood (REML) based on the Bayesian model outperforms generalized cross validation in selecting the smoothing parameters. Following this approach, Guo (2002) showed that SSFMM can also be represented as a mixed effects model and subsequently inference can be carried out based on REML such as likelihood ratio test. The implementation in this article follows the approach of Guo (2002). The SAS macro fmixed is written using SAS/STAT and SAS/IML (SAS Institute Inc. 2008a,b), and can handle flexible design matrices. The items needed for fitting the mixed effects model representation of SSFMM are calculated by SAS/IML, the model fitting is done by PROC MIXED, and the estimates and 95% confidence intervals of the fixed functions and the individual fittings are obtained by rearranging the PROC MIXED output using SAS data steps.
This article is organized as follows. In Section 2 we review the methodological background of SSFMM. In Section 3 we introduce the input parameters for fmixed and its output data sets. In Section 4 we illustrate fmixed by a numeric example. In Section 5 we apply fmixed to a data example. And in Section 6 we give some discussions on the issue of computational burdens.

Methodological background
In this section we briefly review SSFMM proposed by Guo (2002).

The general model
For subject i = 1, . . . , n and time point j = 1, . . . , m i , the observed value y ij is modeled as where β(t) = {β 1 (t), . . . , β p (t)} is a p×1 vector of fixed functions, α i (t) = {α 1i (t), . . . , α qi (t)} is a q × 1 vector of random functions that are modeled as realizations from the same stochastic processes, X ij and Z ij are design matrices which contain the information of experimental designs and covariates, ande ij is the independent and identically distributed error term.
Each element of β(t) is modeled as a cubic smoothing spline with the equivalent Bayesian model as where {b 1k , b 2k } ∼ N(0, τ I), τ → ∞ and W k (s) is a standard Wiener process. Each element of α i (t) is modeled as independent across subjects and as realizations from the same stochastic process where {a 1l , a 2l } ∼ N(0, D) and W l (s) is also a standard Wiener process. And β k (t) and α il (t) are mutually independent and are independent across k and l.
The reproducing kernels corresponding to (2) and (3) are usually difficult to calculate and we will use an alternative Hilbert space on [0, 1], where the reproducing kernels are specified in terms of Bernoulli polynomials (Gu 2002, Section 2.3.3). These kernels are easier to calculate and more numerically stable in the equivalent mixed effects model estimation.
We define the corresponding reproducing kernel through direct sum: R = R 0 ⊕ R 1 (Gu 2002, Section 2.3.3). We first define the Bernoulli polynomials and then the reproducing kernels: A cubic spline f (t) can be expanded as the following (Gu 2002, Section 2.3.3) Using this expansion, we can remodel the functional fixed effects and random effects in (2) and (3) as where W k (s), W l (s) and other coefficients are the same as before.

The mixed effects model representation
Similar to Wang (1998), Equation 6 and Equation 7 can be represented as equivalent mixed effects model at the design points t i = {t i1 , . . . , t i,m i }. We define T i as an n × 2 matrix with the first column as one and the second as . Each element of β(t) and α i (t) can be represented as mixed effects models Collecting items in β(t) and α i (t), we have the mixed effects model representation for Equation 1 where and ⊗ denotes the Kronecker product. In this representation,d are the fixed effects andα i are the random effects as in the conventional sense, and this model can be fitted with existing softwares such as SAS PROC MIXED.

Estimation and inference
The smoothing parameters and other variance components are estimated by REML. This approach is justified by two results: for correlated errors Wang (1998) showed that REML outperforms other methods such as generalized cross validation in selecting the smoothing parameters and in estimating other variance components, for independent errors Kohn et al. (1991) had a similar result. Based on REML, inference can be performed by likelihood ratio test and the null distribution of the test statistic can be approximated by simulation.
The spline estimates on the design points are the empirical best linear unbiased predictions (EBLUPs) conditional on the REML estimates. The function values other than at the design points can also be estimated as EBLUPŝ Their confidence intervals can be constructed similarly.
3. The input parameters and the output data sets

Input parameters
The following is a list of the input parameters for fmixed.
response_var: Specifies the dependent variable, or the "y" variable.
t_var: Specifies the independent variable, or the "time" variable.
sid: Specifies the variable that identifies the subjects in the model. fixed_effects = %bquote(): Specifies the design matrix for the functional fixed effects, or the "X " matrix. There are three scenarios: -fixed_effects = %bquote (0): A fixed linear trend will be fitted.
-random_effects = %bquote(x1,x2): One random cubic smoothing spline will be fitted for effect x1 and one for x2. x1 and x2 are separated by a comma. random_int_slope_type = %bquote(): Specifies the covariance structure for the random intercept and slope of the random splines. This parameter can be either "vc", which stands for variance component, or "un", which stands for unstructured. If this parameter is left blank, the macro will use the default setting as "vc". Note that by using "vc", the random intercept and the random slope will be independent and each will have its own variance. save_location = %bquote(): Specifies a library name to save the output data sets. If this parameter is left blank, the macro will save the output data sets to the default library "work".

Output data sets
Besides the default output from PROC MIXED and the intermediate data sets, macro fmixed generates four data sets containing the results. In the following names of the data sets, "response_var" will be replaced by the actual name of the dependent variable. Those output data sets will be stored at the directory specified by save_location = %bquote() or at the "work" library if no save location is given.
fit_stat_response_var: Contains the fit statistics such as -2 Res Log Likelihood.
cov_parms_response_var: Contains the estimated covariance parameters as follows.
-Fixed Spline:$xk$: The estimated smoothing parameter 1/λ bk for the fixed spline corresponding to variable xk.
-Random Spline:intercept: The estimated variance for the intercept part in the random spline.
-Random Spline:slope: The estimated variance for the slope part in the random spline.
-Random Spline:spline: The estimated variance, 1/λ α , for the smooth part in the random spline.
-Residual: The estimated variance for the error terms.
blup_response_var: Contains BLUPs for all the data points and the 95% confidence intervals with variables pred: Point estimates.
fun_fixed_response_var: Contains the estimated functional fixed effects. Corresponding to the three scenarios in specifying the functional fixed effects, this data set can have variables as the following.
-fixed_effects = %bquote(0): Data fun_fixed_response_var is not created. The results about the linear trend are available in the default output.
* t_var: The independent variable, or the "time" variable. * fun1: The estimated function values. * StdErr_fun1: The standard errors of the estimated function values.
* t_var: The independent variable, or the "time" variable.

A numeric example
In this section we demonstrate how to use fmixed by a numeric example. For subject i = 1, . . . , 20, and equally spaced time points j = 1, . . . , 30 on [0, 1], we generate data as where t = {t 1 , . . . , t 30 }. For each subject, x i1 is assigned by a draw from a uniform distribution on [−0.5, 0.5] and x i2 is assigned as 0 or 1 equally likely, r i (t) is a draw from N(0, σ 2 1 Σ(t))  with σ 2 1 = 900 and Σ = R 1 (t, t), r i (t)s are independent across different subjects, and e i is a draw from N(0, σ 2 2 I) with σ 2 2 = 1. To mimic real applications that subjects usually do not have the same number of observations, the generated data points are then subject to random missing with a probability of 0.2.

Hormone dynamic example
In this section we illustrate fmixed by modeling the dynamics of cortisol in patients with both chronic fatigue syndrome and fibromyalgia syndrome (see Crofford et al. 2004, and references therein). Chronic fatigue syndrome is characterized with unexplained fatigue which leads to substantial reduction in daily performance, and fibromyalgia syndrome is characterized with persistent chronic pain conditions. Both syndromes are stress related and therefore by studying the dynamics of cortisol, a stress related hormone, we may gain a better understanding of the pathophysiological pathways. Guo (2002) compared the cortisol levels between fibromyalgia patients and healthy controls and he found that the fibromyalgia patients have significantly higher cortisol levels. In this section we concentrate on patients diagnosed with both syndromes compared to healthy controls, and we are interested in comparing the group averages and in obtaining subject specific fittings. Figure 3 displays the raw data. There are 12 subjects for each of the patient group and the control group, and each subject has hourly cortisol measurements for a 24-hour period. As shown in the figure, both groups exhibit 24-hour cycles which are difficult to capture by parametric models, and the subject specific deviations also need to be flexibly accounted. Therefore we adopt the same approach as in Guo (2002) and fit the following SSFMM to the data: where y ij is the cortisol concentration for the ith subject and jth hour and i = 1, . . . , 24, j = 1, . . . , 24; X ij = {x ij1 , x ij2 } = 1 0 for the patient group and 0 1 for the control group; β(t) = {β 1 (t), β 2 (t)} where β 1 (t) is the functional group-average for the patient group and β 2 (t) for the control group; α i (t) is the functional subject specific deviation where the random intercept and slope are modeled with a diagonal covariance matrix; and e ij is the error term. The model fitting is done by calling fmixed as %fmixed(data = both, response_var = cort, t_var = time, sid = subject, fixed_effects = %bquote(x1,x2), random_effects = %bquote(1), random_int_slope_type = %bquote(), save_location = %bquote());  The estimated group average curves are displayed in Figure 4. The figure shows that the patient group has higher cortisol levels from 3pm to 10pm and lower cortisol levels from 11pm to 8am compared to the control group. This suggests that patients tend to have higher cortisol levels from middle afternoon until early night and lower cortisol levels during sleep.
To test if the group averages are the same, we reformulate the models so that the reduced model is nested in the full model. The full model is y ij = f 1 (t ij ) + α i (t ij ) + e ij , for controls, where both f 1 (t) and f 2 (t) are smoothing splines. The reduced model is thus that f 2 (t) ≡ 0.
Note that under the mixed effects model representations the working covariance structure of the model defined in Equation 13, which has zero blocks for the covariance between the patient group and the control group, does not include this reduced model's covariance structure, which has nonzero blocks for the covariance between the patient group and the control group, as a reduced form. With the reformulation, however, the reduced model is special case of the full model defined in Equation 14 and Equation 15 by setting d 2 = 0 and λ −1 b2 = 0. The full model is fitted as %fmixed(data = both, response_var = cort, t_var = time, sid = subject, fixed_effects = %bquote(x1,x3), random_effects = %bquote(1), random_int_slope_type = %bquote(), save_location = %bquote()); where x3 ≡ 1 for all the subjects, and the reduced model is fitted as %fmixed(data = both, response_var = cort, t_var = time, sid = subject, fixed_effects = %bquote(x3), random_effects = %bquote(1), random_int_slope_type = %bquote(), save_location = %bquote()); The difference of -2 log REML between the full model and the reduced model is 5.92. For the null distribution of the test statistic, Crainiceanu and Ruppert (2004) showed that more weight should be put on the Chi-square distribution with smaller degrees of freedom rather than the conventional 50 : 50 mixture of χ 2 2 and χ 2 3 . But even with a 100 : 0 mixture, the pvalue 0.0518 is still not significant at the 0.05 level. The null distributions of the test statistics, on the other hand, can be approximated by simulation. Implementation of the simulation based approach, however, is beyond the scope of this article.
The individual fittings are displayed in Figure 5 which shows that the individual profiles are captured pretty well. To test if a functional random deviation is necessary, we set λ −1 a = 0 by %fmixed(data = both, response_var = cort, t_var = time, sid = subject, fixed_effects = %bquote(x1,x2), random_effects = %bquote(0), random_int_slope_type = %bquote(), save_location = %bquote()); to fit a reduced model. The difference of -2 log REML between the full model and the reduced model is 17.18, which is highly significant even under the conservative 50 : 50 mixture of χ 2 1 and zero. Therefore the smooth random component cannot be dropped from the model.

Discussion
By utilizing the mixed effects model representations, we have implemented SSFMM in a SAS macro. This macro is easy to use and it can handle complex experimental designs and include covariates. The computational burden of this approach, however, is heavy. The mixed effects model representations of the functional fixed effects essentially introduce correlations among all the data points, both within and across subjects. Therefore we need to work with matrices with dimension nm for n subjects and m time points. This requires O(n 2 m 2 ) for storage and O(n 3 m 3 ) for calculating the inverse and the determinant. When the number of total observations is in the hundreds such as the numeric example and the hormone data example, it only takes 10 to 15 minutes for the model fitting on a personal computer with 2.40 GHz Intel i5 CPU and 3.42 GB RAM. But the computational time increases quickly as the number of total observations increases. For example it takes about 8 hours for nm = 1500 for the same model as the numeric example. Grid computing or mainframe computer are able to handle larger data sets, but the computing resources will still be quickly consumed because of the cubic relationship. An alternative is to adopt the state space approach as suggested in Guo (2002). This approach only requires O(n 2 m) for storage and O(n 2 m) for calculating the inverse and the determinant. Qin and Guo (2006) implemented this approach for periodic SSFMM without covariates. Unfortunately a general software package for SSFMM based on the state space approach is still unavailable.