,Choice of Parameters of Restricted Cubic Spline Functions
,!In Logistic Regression


,Frank E. Harrell Jr., Ph.D.
,"Division of Biometry
,Duke University Medical Center


,%April 25, 1990















,"Order No. 0000930153
,National Center for Health Statistics
,*			     

,

Choice of Parameters of Restricted Cubic Spline Functions in Logistic 
,
Regression

Introduction
	The logistic regression model is a powerful tool for analyzingdata from medical, epidemiologic, and other studies.  The logisticmodel for a binary response is stated as follows.  Let Y denote theresponse and X1, X2,..., Xp denote predictor variables.  Let B1, B2,..., Bp denote weights or regression coefficients corresponding to theX's , B0 an intercept parameter, and XB = B0 + B1 X1 + ... + BpXp,where X is the vector of Xs.  The probability of a positive responseis assumed to be (Walker and Duncan, 1967)

		Prob{Y=1| X} = 1 / [1 + expXB].

In its simplest form the logistic model assumes that a given predictorXj   is linearly related to the log odds that Y=1, i.e.,

		logit Prob{Y=1|X} = B0 + BjXj + linear combination ofother variables,

where logit P = log {P/(1P)}, the log odds of response.
	Linearity of  logit P with respect to Xj  is often assessed bystratifying the data by levels of Xj and plotting logit of proportionsof responses by mean Xj in each stratum.  Other methods such aspolynomials or nonparametric smoothing methods are also used.  Thecubic spline function or piecewise cubic polynomial is a flexible toolwith advantages over other methods.  Spline functions allow formalassessments of linearity, yield graphical estimates oftransformations, and flexibly parameterize a relationship between oneor more predictors and the logit of response.
	Stone and Koo (1985) advocated a modification of cubic splines torestrict the function to be linear in the tails, i.e., the intervalsbeyond the other join points or knots.  Let the knots for a predictorX be t1, t2, ..., tk.    The restricted cubic spline (RCS) function ina predictor X is obtained by expanding X into k1 component variablesX1, X2, ..., Xk1 (Devlin and Weeks, 1986), where X1=X and

		Xj+1 = (Xtj)+3  (Xtk1)+3(tktj)/(tktk1)

		                       + (Xtk)+3(tk1tj)/(tktk1), j= 1, ..., k2,  and 
		z+ = z,  z>0
	       	        0,  z0 . 

For k=5 and X varying from 0 to 1, the RCS component variables X2, X3,and X4 are shown below, with t1t5 being .05, .275, .5, .725, .95. The left graph magnifies the lower portion of the curves.















The RCS logistic model allowing, for example, one variable to be asmooth nonlinear function, is
			               k
		logit Prob{Y=1|X} = B0 +   BiXi + other terms.
		                                            i
Stone and Koo (1985) suggest using k=5 knots for many applications,and suggest setting the knots at the 5th lowest  value of  thepredictor X, the 5th largest value of X, and making the interior knotsequally spaced.  Stone (1986) states that the choice of k is much moreimportant than the knot locations.  Durrleman and Simon (1989) provideexamples in which knots are allowed to be estimated as freeparameters, jointly with the regression coefficients.  They found thateven though the "optimal" knots were often far from apriori knotlocations, the model fits were virtually identical.
	Harrell, Lee, and Pollock (1988) provide more reasons for usingRCS regression, and they advocate making knots equally spaced on thebasis of the density of the distribution of X rather than on actualvalues of X.  They found that Stone and Koo's recommendation fordefault knot placement does not work well when X is very skewed.  Thepurpose of this investigation is to study the effect of k and the knotlocations in a variety of situations to provide guidance in the choiceof these parameters.

Choice of k in Fitting a Complex Function
	In this section I present an example using simulated data thatdepicts how the flexibility of RCS functions increases with increasingk.  This example also shows that when the underlying function iscomplex, k will have to be greater than the usual default of 5.  
	Assume an underlying population logistic model given by

		logit Prob{Y=1|X} = f(X),
where
		f(X) = 2/[1 + {8 |X.45|}1.8]  1.

Draw a sample of n=1500 random uniform U(0,1) values of X, X1, ...,Xn,  and let Y1, ..., Yn denote realizations of Bernoulli randomvariables with probabilities given by

		Pi = 1/[1+expf(Xi)] .

In the sample generated, the frequency of Y=1 was d=614.  Figure 1shows the true logit regression function f(X) along with logisticregression models fitted using maximum likelihood estimation fork=3,...,8.  For k<7, knots were placed at the .05 and .95 quantiles ofX and at k2 equally spaced quantiles in between.  For k7, outerquantiles were set at .025 and .975.  Note that the generated samplecontained usual sampling variation and was not constrained toprecisely reflect the true regression function f(X).
	The function f(X) rises to a narrow peak and thus is not easilyfitted using polynomials.  The spline fits for k=3 and k=4 are toosmooth to be able to "climb" this peak.  Fits for k>4 becomeincreasingly flexible, although fits for large k begin to display badtail behavior from estimation of too many parameters.   Eight knotswere required to fit the peak well.

Effect of Knot Placement for k=5
	To verify Stone's contention that the RCS fit is relativelyinsensitive to the choice of knot positions (see also Durrleman andSimon, 1989), I let the knots vary at random.  For the regressionfunction f(X) used in the previous section, a sample of n=500 X's wasdrawn, and 500 Y's were generated (228 had Y=1).  The X's and Y's werethen fixed and 12 sets of t1, ..., t5  were generated from the U(0,1)distribution.  Figure 2 depicts the 12 fits, with the true f(X) shownin bold lines.  When most knots are far from the center, the fits aredamped.  When the knots are clustered, the fits appear worse. However, the fits do not vary wildly over the sets of random knots. 
	Figure 2 provides evidence that the knots should be restricted tobe roughly uniformly spaced.  I suggest that for k=5, default knotlocations be at the .05, .275, .5, .725, and .95 quantiles of X.  Forn<100, the 5th smallest and largest values replace the outerquantiles.  To study this choice of knots, I fitted the same data usedin Figure 2 with knots that were moderately constrained to not crossmidpoints between these quantiles.  The following table provides theuniform distributions used for generating the knots.

!@@!@@"""  i#Distribution of ti
1#U(.01,.1625)
2#U(.1625,.3875)
3#U(.3875,.6125)
4#U(.6125,.8375)
5#U(.8375,.99)

The results are shown in Figure 3.  As expected, the fits are moreuniform.  When the center knot was near .45, the fit was more peaked. Overall, the fits are relatively insensitive to the knot locations.

Choice of k as a Function of n and d
In binary logistic regression, the number of variables that can be modeledreliably is principally a function of the sample size n and the numberof responses d (if d>n/2, substitute the number of nonresponses ndfor d).  It  is reasonable to expect that the maximum number of knotsthat can be modeled without inflating the variance of predictions willalso be a function of n and d.  Simulation studies are currentlyunderway to develop databased prescriptions for k.

Examples using Real Data##--77AAKK
A dataset of patients from intensive care units at several hospitals wasprovided courtesy of Dr. Douglas Wagner of George WashingtonUniversity's ICU Research Unit.  This dataset consisted of n=1000patients of whom 234 died in the hospital.  For each patient, variousphysiologic variables were measured on the first day of the hospitalstay.  Figures 410 display various spline logistic fits for the logodds of hospital death.  In each figure, the two upper plots representfits with k=4 and k=5 default knots, respectively (i.e., usingquantiles equally spaced between .05 and .95).  The left plot in thebottom row depicts Stone and Koo's knot placement for k=5: knots atthe 5th smallest and largest values of X and equally spaced knots onthe Xscale.  The last plot corresponds to k=6 using default knotplacement.  Knots are indicated with arrows on the xaxes.
Various statistics appear within each plot.  L.R. refers to the likelihoodratio {{chi}}2 for the logistic model.  The Wald test statistic fortotal association is a {{chi}}2 statistic testing the same hypothesisas L.R.  Below it is found the Wald {{chi}}2 for testing linearity ofthe fitted regression function.  AIC refers to Akaike's informationcriterion, L.R.  2 x (number of regression parameters).  AIC ispenalized for complexity of the model.  In many situations the modelwith the highest AIC has the greatest chance of predicting accuratelyin future samples.
The relation between temperature (0C) and log odds of hospital death is shownin Figure 4.  The fit with 4 default knots has the largest AIC andslightly smaller confidence bands than the other 3 models.  Fits donot vary greatly.
Figure 5 depicts relationships with  heart rate (beats/minute).  The k=6 modelhad the best AIC and a flatter relationship around the normal range ofheart rate.  Stone and Koo's knot placements had the next best AIC.
The analysis of respiratory rate (/minute) is in Figure 6.  The k=6 model hadthe best AIC, but k=4 was better than k=5.  Stone and Koo's methodapparently placed knots too far from the nadir.
Analyses of serum creatinine are shown in Figure 7.  The k=4 model had the bestAIC, followed by k=5.  Stone and Koo's knots were too far in the righttail of this highly skewed variable.  Fits of the other 3 models wereremarkably similar.
Mean blood pressure plots are in Figure 8.  The k=6 model had the highest AICbecause it picked up a sharper dip as well as curvature in the leftand right portions of the curve.  The k=5 and Stone and Koo modelswere not far behind, but the k=4 model was too smooth.
Urine output analyses are in Figure 9.  Again, Stone and Koo's knots were toofar to the right in this skewed variable.  The k=5 model had the bestAIC, and fits for the 3 quantilebased methods are very similar.
Figure 10 shows analyses for BUN.  Stone and Koo's knot placement yielded a fitthat is too smooth due to insufficient knots in the left half of theskewed distribution of BUN.  The k=6 model had by far the best AICbecause it was able to fit an additional hump around BUN=13.  
A different way to study the effect of k and the knot locations is to study theglobal influence on the fit of radically changing data in oneinterval.  For the BUN example of Figure 10, 29 of the 66 patientswith BUN>60 died.  If we pretend that all patients with BUN>60 lived,we obtain the fits in Figure 11.  The right portion of each curve isrightfully pulled down.  The far left tail of the k=4 curve turned upslightly, and that of the k=5 curve turned down slightly.  However,the left portion of the k=6 fit remained the same.  Apparently, 6knots allowed enough flexibility in the left and middle portions ofthe curve to overcome a drastic change in the right segment.

Software Implementation of RCS Regression
Once the k2 nonlinear component variables are computed, RCS regressioninvolves no extra computations.  Thus ordinary least squares ormaximum likelihood estimation programs can be used to obtain RCS fitsand associated test statistics.  RCS tools have been implemented asfunctions in the New S language as well as SAS macros.  These toolswork well in conjunction with a variety of S and SAS regressionmodeling procedures.  For standalone software systems, new moduleswill have to be added to automatically expand chosen predictors intospline components.  Automatic knot selection (as a default) isessential, so the software will have to compute quantiles on  thepredictors and eliminate nonunique knots in the case of ties in thedata.
Graphs such as those in Figure 11 were obtained in S by statements such as

rcspline.plot(x,y,nk=no. knots),

which uses the S functions glim and coxreg written by Dr. TerryTherneau.  The SAS statement

%PSPLINET(x,y,nk=k),

which uses the LOGIST and PHGLM procedures, produces similar output(though without the analysis of loglikelihood table).

Summary
The examples and simulations presented support the idea that k=4 or 5 knots areusually sufficient unless the underlying regression function is verycomplex.  RCS regression is not very sensitive to knot locations aslong as the locations are crudely equally spaced where there are datapoints.  Knots can be selected without referencing Y unless the curvetakes a sudden turn.  Since RCS regression involves only an enhanceddesign matrix, its implementation in a variety of software environmentis simple.

References
!@@"""  1@@"""
Walker SH, Duncan DB (1967).  Estimation of the probability of anevent as a function of several independent variables.  Biometrika54:16779.N

Stone CJ, Koo CY (1985).  Additive splines in statistics.  Proc StatComputing Sect ASA 458.N

Devlin TF, Weeks BJ (1986).  Spline functions for logistic regressionmodeling.  Proc Eleventh Ann SAS Users Group Intnl Conf.  Cary NC:SAS Institute, Inc., 64651.N

Stone CJ (1986).  Comment: Generalized additive models.  Stat Sci1:31214.N

Durrleman S, Simon R (1989).  Flexible regression models with cubicsplines.  Stat in Med 8:55161.N

Harrell FE, Lee KL, Pollock BG (1988).  Regression models in clinicalstudies: determining relationships between predictors and response. J Natl Cancer Inst 80:11981202.N



