SAS Global Forum 2007 Posters

1 Paper 143- 2007 regression Splines with Longitudinal DATA Chan Hee Jo, University of Arkansas Medical Sciences, Little Rock, AR Jeff Gossett, University of Arkansas Medical Sciences, Little Rock, AR Pippa Simpson, University of Arkansas Medical Sciences, Little Rock, AR ABSTRACT In many clinical trial studies, patients are observed and/or measured on multiple occasions. To account for the longitudinal nature of the data, a mixed model analysis implemented using SAS PROC MIXED is commonly used. It is typical to make comparisons between dose or treatment groups, possibly controlling for demographic variables. When the dose response curve is non-linear, it can be difficult to find an appropriate curve, and model misspecification is common. We will present an alternative to non-linear mixed models that avoids the problem of specifying a parametric model.

2 We demonstrate the use of a semi-parametric regression in the mixed model framework. This use of PROC MIXED allows for nonlinear fitting with a covariance structure for the correlated observations. We illustrate its use from the longitudinal data with various time points. Keywords: mixed models, regression splines, random coefficients, longitudinal. INTRODUCTION In many clinical trial studies, patients are observed and/or measured on multiple occasions. To account for the longitudinal nature of the data, a mixed model analysis implemented using SAS PROC MIXED is commonly used. For non-linear curves, it can be difficult to find an appropriate curve. We present an alternative to non-linear mixed models namely spline regression . We use a semi-parametric regression in the mixed model framework. This use of PROC MIXED allows for nonlinear fitting with a covariance structure for the correlated observations.

3 We illustrate its use with simulated data. SMOOTHING SPLINES AND regression SPLINES BACKGROUND Assume that the responses Y1, .., Yn are obtained at X1, .. Xn of an independent variable X. A standard nonparametric regression model is given as Yi = m(Xi) + i _ i =1,.. n, where m is an unknown regression function and the i s are zero mean, uncorrelated random errors with common variance 2 . There are a number of estimators that can be used for the regression function m with spline smoothing providing one popular choice. We use m to denote the regression function and g as the spline model function. For a given >0 , a smoothing spline is usually derived as the minimizer of the penalized sum of squares criterion over all functions g in the space Our smoothing spline method becomes less practical when the sample size increases, because the smoothing splines use all observations as knots.

4 On the other hand, a regression spline uses a fixed number of knots. To estimate m one can use a regression spline and the general form of a regression spline is: 1 SASG lobalForum2007 Posters where p( 1) is an integer, ( 1+p .. l+p ) is a vector of regression coefficients, zp+= zpI* I(z>0). In the above equation, 1 .. l represent fixed knots, which are fixed points in the range of the independent variable x. In general, the knots for a regression spline are on a grid of equally-spaced sample quantiles and they are used in the ( )p+ function. Ruppert and Carroll (2000) suggested a penalty approach like smoothing splines but by putting a roughness penalty on the squares of 1+p .. l+p. Using McCulloch and Searle s (2001) equation, let D be the l l diagonal matrix with p+ 1 zeroes followed by l ones along the diagonal.

5 Then we define the estimator of 1+p .. l+p to be the minimizer of where is a smoothing parameter. The larger the value of the smoothing value, the more the spline fit is shrunk towards smoothing curve fit. Simple calculations show that the estimator is given by where B is the n p+1 design matrix with ith row equal to The estimator above is called a ridge regression estimator since it moves the regression spline towards the least squares fit to a polynomial model (See Hastie and Tibshirani, 1990). It is quick to compute especially if one uses a algorithm by Ruppert (2002). The average squared residual (ASR) at given is and the hat matrix S is Then to select a smoothing parameter , one can use the generalized cross-validation (GCV) criterion: and choose the solution which makes that criterion smallest over a grid of values.

6 From a Bayesian point of view, the smoothing parameter is generated by Bayesian machinery by itself and thus it justifies the use of regression P splines. Assume the prior density for vector is a normal with mean 0 and variance -1 TD then the resulting posterior distribution of will be It turns out that the regression spline estimation is equivalent to a mixed effects models. We use this method to set up the model using PROC MIXED. NONPARAMETRIC BASIS TERM To utilize (x- ) p+ as a random effect for the PROC MIXED framework, we will build the basis function. Firstly, the number of knot points l must be determined to build the function. Ruppert (2002) found that the choice of the number of knots, l is not important, provided that l is at least a certain minimum value, generally 20 for most situations in clinical trials setting.

7 In our case, set of l=10 equally spaced knots are selected in this paper. Using SAS Macro, one can find sample 10th quantiles, xknots, of x for knots. Using array function in SAS, one can build the pth regression spline basis, xknot. We used p=2 so quadratic spline basis functions, and thus x and x*x for a linear and quadratic fixed terms in our regression spline set up. Proc Mixed data=LMM method=reml; 2 SASG lobalForum2007 Posters Model y=x+x*x/solution; Random xknots/solution type=toep(1); Type specifies the covariance structure of nonparametric basis term, the random effects, xknots, and we can compare the results for the different covariance structures. We use a diagonal variance covariance matrix for the random effects. LONGITUDINAL DEPENDENT VARIABLE WITH VARYING COEFFICIENT MODELS For many longitudinal studies, the covariate variables are independent of time, so that their observations are cross-sectional.

8 But the effects of covariates on the dependent may not be constant over time, in particular, the group effect may not be the same. So to incorporate the effects of time and covariates, one can consider the time varying coefficient model as below: Yi (t) = 0 (t)*Groupi+ f(Xi) + i(t) _ i =1,.. n, t=1,..,mi, where f(Xi) = quadratic regression spline funtion with 10 knots. Note that for simplicity, we consider only 0 is time varying but the rest of the s are constant over time. Very often, failure to take into account correlations between repeatedly measured variables will result in biased estimation of betas. One of the major advantages using PROC MIXED is to allow both splines and autocorrelations at the same time. You can add the REPEATED statement with specified covariance structure to the PROC Mixed.

9 In addition, we consider the time design points to be different from subject to subject. For this paper, we consider some of measurements to be missing at random. SIMULATION To evaluate the performance of Proc MIXED for our nonlinear fitting with a covariance structure for the correlated observations, we consider the following model: Yi (t) =b1 (t)*Groupi+ m(Xi) + i(t) _ i =1,.. n, t=1,..,mi where Group and X are independent random variables. Two categories are set up for group variable but it can be extended to many. X is assumed to be normally distributed with mean 0 and variance 4. For the simplicity, we assume that we have nonlinear function of X which is time invariant. The simulation process randomly generated n patients and their corresponding covariate vectors of (groupi, Xi), i=1.

10 N. The errors i(t) s were generated from multivariate normal distribution with the covariance structure: 1 2 .. m-1 1 .. m-2 Correlation within subjects: .. 1 .. m-1 .. 1 Correlation between subjects = 0 Using SAS/IML which contains variable matrix operations and built-in functions for statistical computation, we are able to generate auto-correlated multivariate normal random samples for errors as shown in the SAS code below: l=t(root(eps)); /* calculate cholesky root of cov matrix eps */ z=normal(j(1,n,seed)); /* generate n samplesize normals */ 3 SASG lobalForum2007 Posters xx=l*z*v; /* multiply by cholesky root and variance */ tx=t(xx); do i=1 to Group[i] = round(ranuni(int(time()))); x[i] = rannor(int(time()))*2 ; do r=1 to b0[r]= *cos((r-& )*/ ); yy[i,r]=b0[r]*Group[i]+ *sin(2* *x[i]/5); end; end; yy=yy+tx.