Table Of ContentModel Diagnostics for Smoothing Spline
ANOVA Models
Chong Gu(cid:3)
Abstract
Functional ANOVA decompositions can be incorporated in multivariate function estimation
through the penalized likelihood method. In this article, we propose some simple diagnostics
for the \testing" of selected model terms in the decomposition; the elimination of practically
insigni(cid:12)cant terms generally enhances the interpretability of the estimates, and sometimes may
also have inferential implications. What we try to achieve are the tasks of the traditional
likelihood ratio tests, but in the absence of sampling distributions due to the typically in(cid:12)nite
dimensional nulls in nonparametric settings. The diagnostics are illustrated in the settings of
regression, probability density estimation, and hazard rate estimation.
KEY WORDS: ANOVA decomposition; Diagnostic; Kullback-Leibler projection; Penalized
likelihood estimate.
1 Introduction
Penalizedlikelihoodmethodestimatesafunctionofinterest, say(cid:17), bytheminimizerofafunctional
of the form
L((cid:17)jdata)+(cid:21)J((cid:17)) (1)
where L((cid:17)jdata), usually taken as the minus log likelihood, measures the goodness-of-(cid:12)t of (cid:17) to
the data, and J((cid:17)), usually taken as a quadratic roughness functional, measures the smoothness
(regularity)of(cid:17); thesmoothingparameter(cid:21)controlsthetradeo(cid:11)betweenthetwoconflictinggoals.
Comprehensive treatments of penalized likelihood estimation can be found in, e.g., Wahba (1990)
and Gu (2002).
On a product domain, say X = X (cid:2)X , one may decompose the function as
1 2
(cid:17)(x) = (cid:17)(xh1i;xh2i) = (cid:17);+(cid:17)1(xh1i)+(cid:17)2(xh2i)+(cid:17)1;2(xh1i;xh2i);
(cid:3)Chong Gu is Associate Professor, Department of Statistics, Purdue University, West Lafayette, IN 47907. The
authoracknowledgestheassistanceofJingyuanWangintheearlystageofthisresearch. Thisresearchwassupported
by National Institutes of Health under grant R33HL68515.
1
where x = (x ;x ) 2 X with x 2 X and x 2 X and (cid:17) , (cid:17) , and (cid:17) satisfy certain
h1i h2i h1i 1 h2i 2 1 2 1;2
side conditions to ensure identi(cid:12)ability; extensions to higher dimensions are straightforward. Such
decompositionextendstheclassicalANOVAdecompositionondiscretedomainstogenericdomains,
and is simply referred to as (functional) ANOVA decomposition. ANOVA decompositions can be
builtinpenalizedlikelihoodestimationusingthetensorproductsplinetechnique; seethereferences
givenaboveandx2.2. AswithclassicalANOVAmodels,regressionmodelsconsistingofmaine(cid:11)ects
and selected lower order interactions are the most useful in practical applications.
For the estimation of a probability density f(x;y) of random variables (X;Y) on domain
R
X (cid:2) Y, one may employ the logistic density transform f = e(cid:17)= e(cid:17) and estimate (cid:17)(x;y) =
X(cid:2)Y
(cid:17) (x)+(cid:17) (y)+(cid:17) (x;y), where the constant is eliminated for a one-to-one transform. The absence
x y x;y
of the interaction (cid:17) characterizes the independence of the random variables X and Y. More
x;y
generally, the exclusions of selected interactions in higher dimensions may characterize various
(conditional) independence structures among random variables; see, e.g., Gu (2002, x1.3.3) for a
few more examples.
For the estimation of hazard rate based on censored lifetime data, one may employ the log
transform and estimate the log hazard (cid:17)(t;u) = (cid:17); + (cid:17)t + (cid:17)u + (cid:17)t;u as a function of time t and
covariate u; the absence of (cid:17) characterizes a proportional hazard model.
t;u
The purpose of this article is to devise and illustrate some simple diagnostics for practically
insigni(cid:12)canttermsinanANOVAdecompositionestimatedthroughthepenalizedlikelihoodmethod.
The diagnostics are based on Kullback-Leibler geometry, and are implemented and illustrated in
the settings of regression, density estimation, and hazard estimation.
The rest of the article is organized as follows. In x2, detailed formulation of penalized likelihood
estimation is briefly reviewed. The diagnostics are derived in x3 and illustrated in x4; x5 reports
some supplementary simulation results. Brief discussion in x6 concludes the article.
2 Penalized Likelihood Estimation
The likelihood term L((cid:17)jdata) in (1) reflects the sampling structure of the observed data, and the
penalty term J((cid:17)) largely characterizes the model to be (cid:12)tted to the data. Some details concerning
the background formulations are given below.
2.1 Likelihood
First consider regression problems with data from exponential families,
Yjx (cid:24) expf(y(cid:18)((cid:17)(x))−b((cid:18)((cid:17)(x))))=a((cid:30))+c(y;(cid:30))g;
2
where (cid:17) is a monotone transform of the canonical parameter (cid:18), taking values in (−1;1), and
a((cid:30)) is a common dispersion parameter; a unrestricted range of (cid:17) prevents the complication of
constrained minimization of (1). Observing (x ;Y ), i = 1;:::;n, one may take
i i
Xn
1
L((cid:17)jdata) = − fY (cid:18)((cid:17)(x ))−b((cid:18)((cid:17)(x )))g; (2)
i i i
n
i=1
where the term c(y;(cid:30)), independent of (cid:17), is dropped and the dispersion a((cid:30)) is to be absorbed
into the smoothing parameter (cid:21). Examples of this setting include least squares regression with (cid:17)
the normal mean, logistic regression with (cid:17) the logit, Poisson regression with (cid:17) the log intensity,
etc.; see, e.g., Gu (2002, x5.4) for further details. The general formulation of penalized likelihood
regression with exponential family responses was due to O’Sullivan, Yandell, and Raynor (1986).
See, e.g., Gu (2002, Ch. 5) for a comprehensive treatment and further references. See also Gu and
Kim (2002).
For probability density estimation on a bounded domain X with independent samples X ,
i
i = 1;:::;n, the minus log likelihood is simply
Xn n Z o
1
L((cid:17)jdata) = − (cid:17)(X )−log e(cid:17)(x)dx ; (3)
i
n
X
i=1
R
wherethelogisticdensitytransformf = e(cid:17)= e(cid:17) isemployed; aboundeddomainpermitsaproper
X
uniform distribution with (cid:17) = 0. Penalized likelihood density estimation was pioneered by Good
and Gaskins (1971); see also Silverman (1982) and Gu and Qiu (1993). A comprehensive treatment
and further references can be found in Gu (2002, Ch. 6).
Let T be the lifetime of an item with survival function S(tju) = P(T > tju), possibly dependent
on a covariate U. The hazard function is given by e(cid:17)(t;u) = −@logS(tju)=@t. Let Z be the left-
truncation time and C be the right-censoring time, independent of T and of each other. Observing
(Ui;Zi;Xi;(cid:14)i), i = 1;:::;n, whereX = min(T;C), (cid:14) = I[T(cid:20)C], andZ < X, theminusloglikelihood
is seen to be 1 Xn n Z Xi o
L((cid:17)jdata) = − (cid:14) (cid:17)(X ;U )− e(cid:17)(t;Ui)dt : (4)
i i i
n
i=1 Zi
The formulation through (4) was found in Gu (1996), which covers various earlier models as special
cases; see, e.g., Gu (2002, Ch. 7) for a comprehensive treatment and further references.
Whenparametricmodelsareassumedonthetimeaxis, sayinthecaseofacceleratedlifemodels
with T~ = logT following a location-scale family distribution with density f((t~−(cid:22))=(cid:27)), the task is
to estimate a parameter of the parametric model as a function of the covariate u; it is essentially a
regression problem. For example, for T Weibull, the hazard is given by ((cid:11)=t)(t=(cid:12))(cid:11) = e(t~−(cid:22))=(cid:27)=(cid:27)t
3
where (cid:11) = 1=(cid:27) is the shape parameter and (cid:12) = e(cid:22) is the scale parameter, and the minus log
likelihood, upon setting (cid:17) = (cid:22), is seen to be
Xn (cid:8) (cid:9)
1
L((cid:17)jdata) = − (cid:14) ((cid:11)(logX −(cid:17)(U ))+log(cid:11))−(X(cid:11)−Z(cid:11))e−(cid:11)(cid:17)(Ui) : (5)
n i i i i i
i=1
Other examples of accelerated life models include the log normal and the log logistic distributions;
see, e.g., Gu (2002, x7.5) for further details.
2.2 Penalty
We now discuss the construction of J((cid:17)). To keep the notation simple, the domain of (cid:17) is to be
denoted by X regardless of the stochastic setting.
The minimization of (1) is in a space H (cid:18) f(cid:17) : J((cid:17)) < 1g in which J((cid:17)) is a squre semi-
norm. The evaluation functional [x](cid:17) = (cid:17)(x) appears in the term L((cid:17)jdata), and is assumed to be
continuous in H. A space H in which the evaluation is continuous is called a reproducing kernel
Hilbertspace(RKHS)possessingareproducingkernel(RK)R((cid:1);(cid:1)),anon-negativede(cid:12)nitefunction
satisfying R ((cid:1)) = R(x;(cid:1)) 2 H, 8x 2 X, and hR(x;(cid:1));f((cid:1))i = f(x), 8f 2 H, where h(cid:1);(cid:1)i is the inner
x
product in H; the RK R((cid:1);(cid:1)) and the space (H;h(cid:1);(cid:1)i) determine each other uniquely. Typically,
h(cid:1);(cid:1)i = J((cid:1);(cid:1))+J~((cid:1);(cid:1)), where J((cid:1);(cid:1)) is the semi inner product associated with J((cid:1)) and J~((cid:1);(cid:1)) is an
inner product in the null space N = f(cid:17) : J((cid:17)) = 0g of J((cid:17)) when restricted therein. There exists
J
a tensor sum decomposition H = N (cid:8)H , where the space H has J((cid:17)) as its square norm and
J J J
an RK R satisfying J(R (x;(cid:1));f((cid:1))) = f(x), 8f 2 H . See, e.g., Gu (2002, x2.1).
J J J
R
An example on X = [0;1] is J((cid:17)) = 1(cid:17)¨2dx with N = f(cid:17) : (cid:17)(x) = (cid:12) + (cid:12) xg; this yields
0 R J R R 0 R 1
the popular cubic splines. A choice of J~(f;g) is ( 1fdx)( 1gdx) + ( 1f_dx)( 1g_dx), yielding
R R 0 0 0 0
H = f(cid:17) : 1(cid:17)dx = 1(cid:17)_dx = 0; J((cid:17)) < 1g and the RK R (x;y) = k (x)k (y)−k (x−y), where
J 0 0 J 2 2 4
k = B =(cid:23)! are scaled Bernoulli polynomials. See, e.g., Gu (2002, x2.3.3).
(cid:23) (cid:23)
For X a product domain, tensor product RKHS can be constructed from RKHS on marginal
domains, and ANOVA decompositions can be induced from one-way decompositions on marginal
domains; the construction of J((cid:17)) is indirect through that of the RK R . The construction can be
J
P
denoted by H = (cid:8)g H and J((cid:17)) = g (cid:18)−1J ((cid:17) ), where (cid:17) 2 H , 0 < (cid:18) < 1, and J ((cid:1)) is
(cid:12)=0 (cid:12) (cid:12)=1 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)P (cid:12)
the square norm in H , (cid:12) > 0. One has N = H , H = (cid:8)g H , and R = g (cid:18) R , where
(cid:12) J 0 J (cid:12)=1 (cid:12) J (cid:12)=1 (cid:12) (cid:12)
R is the RK in H . The (cid:18) ’s are an extra set of smoothing parameters, often suppressed in the
(cid:12) (cid:12) (cid:12)
notation, that adjust the relative penalties on the roughness of di(cid:11)erent components. See, e.g., Gu
(2002, x2.4).
As an example, consider the tensor product cubic splines on X = [0;1]2. On the marginal
4
domain [0;1], one has
(cid:8) R (cid:9)
(cid:17) : 1(cid:17)¨2dx < 1 = H (cid:8)H (cid:8)H
0 00 01 1
(cid:8) R R R (cid:9)
= spanf1g(cid:8)spanfk (x)g(cid:8) (cid:17) : 1(cid:17)dx = 1(cid:17)_dx = 0; 1(cid:17)¨2dx < 1 ;
1 0 0 0
where k (x) = x−0:5, with RKs R (x;y) = 1, R (x;y) = k (x)k (y), and R = k (x)k (y)−
1 00 01 1 1 1 2 2
R
k (x−y); anone-wayANOVAisbuiltinwithf 2 H (cid:8)H satisfyingthesidecondition 1fdx = 0.
4 01 1 0
Taking tensor product, one has 9 tensor sum terms H = H ⊗H on X = [0;1]2, (cid:23);(cid:22) =
(cid:23);(cid:22) (cid:23)h1i (cid:22)h2i
00;01;1, with RKs R (x;y) = R (x ;y )R (x ;y ). The 4 subspaces with (cid:23);(cid:22) = 00;01
(cid:23);(cid:22) (cid:23)h1i h1i h1i (cid:22)h2i h2i h2i
are of one-dimension each, and can be lumped together as H = N . The other 5 subspaces can
0 J
be indexed as H , (cid:12) = 1;:::;5. For the ANOVA decomposition, H contains the constant,
(cid:12) 00;00
H (cid:8)H contains the x main e(cid:11)ect, H (cid:8)H contains the x main e(cid:11)ect, and the
01;00 1;00 h1i 00;01 00;1 h2i
other 4 contain the interaction. Excluding the 4 subspaces corresponding to the interaction, one
gets an additive model. See, e.g., Gu (2002, Example 2.8).
2.3 Computation
When the term L((cid:17)jdata) involves only evaluations such as in (2) and (5), it can be shown that the
P P
minimizer of (1) has an expression (cid:17)(x) = m d (cid:30) (x)+ n c R (x ;x), where f(cid:30) g is a basis
(cid:23)=1 (cid:23) (cid:23) i=1 i J i (cid:23)
of N . When L((cid:17)jdata) involves more than evaluations such as in (3) and (4), the exact minimizer
J
of (1) does not appear to be computable. In either case, it can be shown that the minimizer of (1)
in a space N (cid:8)spanfR (z ;(cid:1));j = 1;:::;qg shares the same asymptotic convergence rates as the
J J j
minimizer in H, where fz g is a random subset of fx g and q can be as small as O(n2=9) for (tensor
j i
product) cubic splines; see, e.g., Gu and Wang (2002) and Gu and Kim (2002). Plugging into (1)
the expression
Xm Xq (cid:16)Xg (cid:17)
(cid:17)(x) = d (cid:30) (x)+ c (cid:18) R (z ;x) (6)
(cid:23) (cid:23) j (cid:12) (cid:12) j
(cid:23)=1 j=1 (cid:12)=1
and(cid:12)xing(cid:18) and(cid:21),onesolvesforthecoe(cid:14)cientsd andc . TheANOVAtermsareeasilyextracted
(cid:12) (cid:23) j
by regrouping the terms of (6). For the selection of the smoothing parameters (cid:18) and (cid:21), various
(cid:12)
cross-validation scores have been developed to target the Kullback-Leibler losses derived in the
various stochastic settings; see, e.g., Craven and Wahba (1979), Xiang and Wahba (1996), and Gu
and Wang (2002).
Most of the regression, density and hazard models discussed here are implemented in the R
package gss by the author; simple illustrations of the usage of gss can be found in the examples
of x4.
5
3 Diagnostics
We are now ready to derive the diagnostics in the settings of regression, density estimation, and
hazard estimation.
3.1 General Method
Suppose the estimation of (cid:17) has been done in a space H, but in fact (cid:17) 2 H(cid:3) (cid:26) H. The task is to
devise tools to \test" the null hypothesis that (cid:17) 2 H(cid:3). In a parametric analysis, the likelihood ratio
test or the like serve as the primary tool for the purpose. In a nonparametric setting, the null H(cid:3)
is typically in(cid:12)nite-dimensional, and consequently the sampling distribution of the likelihood ratio
statistic is no longer available, hence one has to look elsewhere.
We propose some heuristic diagnostics based on Kullback-Leibler geometry. Denote by (cid:17)^ the
estimate of (cid:17) obtained by minimizing (1) in H, with the smoothing parameters selected by cross-
validation. Let (cid:17)~ be the Kullback-Leibler projection of (cid:17)^ in H(cid:3), the minimizer of KL((cid:17)^;(cid:17)) for
(cid:17) 2 H(cid:3), and let (cid:17)c be the constant model; (cid:17)c is the maximum likelihood estimate of (cid:17) = (cid:17); for
the regressions of (2) and (5) and for the hazard estimation of (4), and is (cid:17) = 0 for the density
estimation of (3). It typically holds that
KL((cid:17)^;(cid:17) ) = KL((cid:17)^;(cid:17)~)+KL((cid:17)~;(cid:17) ); (7)
c c
where KL((cid:17)^;(cid:17) ) is the \total entropy" of (cid:17)^ and KL((cid:17)~;(cid:17) ) is the \predicted entropy" by the null
c c
P P
H(cid:3), very much like the total variation (Y −Y(cid:22))2 and the predicted variation (Y^ −Y(cid:22))2 in
i i i i
standardleastsquaresregression. Theratio(cid:26) = KL((cid:17)~;(cid:17) )=KL((cid:17)^;(cid:17) ),verymuchlikeanR2 statistic,
c c
indicates how much of (cid:17)^ actually sits in the null H(cid:3), and shall be used to diagnose the feasibility
of the null hypothesis (cid:17) 2 H(cid:3).
3.2 Exponential Family Regression
For regression with exponential family responses as in (2), one has
Xn
1
KL((cid:17)^;(cid:17)) = f(cid:22)^(x )((cid:18)((cid:17)^(x ))−(cid:18)((cid:17)(x )))−(b((cid:18)((cid:17)^(x )))−b((cid:18)((cid:17)(x ))))g; (8)
i i i i i
n
i=1
where (cid:22)(x) = (db=d(cid:18))(x) = E[Yjx]. Substituting in the expression (6) with m and g the versions
associated with H(cid:3), the minimization of (8) can be solved via Newton iteration. The resulting
KL((cid:17)^−(cid:17)~) depends on the (cid:18) ’s hidden in (6), and an outer loop of optimization shall be performed
(cid:12)
to minimize KL((cid:17)^−(cid:17)~) with respect to the (cid:18) ’s.
(cid:12)
6
Note that after substituting in (6) one is no longer projecting into H(cid:3) but only into some
(m+q)-dimensional subspace, indexed by (cid:18) , therein. The outer loop optimization serves to locate
(cid:12)
the nearest (to (cid:17)^) of such (m+q)-dimensional subspaces.
Setting (cid:17) = (cid:17)~+(cid:11)((cid:17)~−(cid:17) ) in (8) for (cid:11) real, di(cid:11)erentiating with respect to (cid:11), and evaluating at
c
(cid:11) = 0, one has
Xn
1
((cid:22)~(x )−(cid:22)^(x ))h~(x )((cid:17)~(x )−(cid:17) (x )) = 0; (9)
i i i i c i
n
i=1
where h = d(cid:18)=d(cid:17); note that (cid:17)~ minimizes (8) in H(cid:3) and (cid:17)~−(cid:17) 2 H(cid:3). When (cid:17) = (cid:18) is the canonical
c
parameter, (9) yields (7) for the Kullback-Leibler distance de(cid:12)ned in (8); this covers the cases of
Gaussian regression, logistic regression, and Poisson regression.
3.3 Density Estimation
For density estimation as in (3), one has
Z Z
KL((cid:17)^;(cid:17)) = (cid:22) ((cid:17)^−(cid:17))−log e(cid:17)^(x)dx+log e(cid:17)(x)dx; (10)
(cid:17)^
X X
R R
where (cid:22) (h) = h(x)eg(x)dx= eg(x)dx. One can again substitute the expression of (6) into (10)
g X X
and perform an inner loop Newton iteration nested under an outer loop optimization with respect
tothe(cid:18) ’s. Similarto(9), onehas(cid:22) ((cid:17)~−(cid:17) ) = (cid:22) ((cid:17)~−(cid:17) ), whichyields(7)fortheKullback-Leibler
(cid:12) (cid:17)~ c (cid:17)^ c
distance de(cid:12)ned in (10).
3.4 Hazard Estimation
For hazard estimation as in (4), one has
Z
1 Xn Xi(cid:8) (cid:9)
KL((cid:17)^;(cid:17)) = e(cid:17)^(t;Ui)((cid:17)^(t;U )−(cid:17)(t;U ))−(e(cid:17)^(t;Ui)−e(cid:17)(t;Ui)) dt; (11)
i i
n
i=1 Zi
see, e.g., Gu (2002, x7.2.1). The same procedure can be followed to locate (cid:17)~ that minimizes (11)
in the nearest (m+q)-dimensional subspace characterized by (6). Similar to (9), one has
Z
1 Xn Xi
(e(cid:17)~(t;Ui)−e(cid:17)^(t;Ui))((cid:17)~(t;U )−(cid:17) (t;U ))dt = 0;
i c i
n
i=1 Zi
which yields (7) for the Kullback-Leibler distance de(cid:12)ned in (11).
For Weibull regression as in (5), with the hazard ((cid:11)=t)e(cid:11)(logt−(cid:17)(u)) parametric in t and the
7
location parameter (cid:17)(u) of logT as a function of the covariate U, (11) becomes
Xn
1
KL((cid:17)^;(cid:17)) = (X(cid:11)−Z(cid:11))f(e−(cid:11)(cid:17)^(Ui)−e−(cid:11)(cid:17)(Ui))−(cid:11)e−(cid:11)(cid:17)^(Ui)((cid:17)^(U )−(cid:17)(U ))g: (12)
n i i i i
i=1
Similar to (9), one has
Xn
1
(X(cid:11)−Z(cid:11))(e−(cid:11)(cid:17)~(Ui)−e−(cid:11)(cid:17)^(Ui))((cid:17)~(U )−(cid:17) (U )) = 0;
n i i i c i
i=1
which yields (7) for the Kullback-Leibler distance given in (12).
4 Examples
We now illustrate the diagnostics through a few examples. The computation can be done using the
R package gss.
4.1 AIDS Incubation
TostudytheAIDSincubationtime,avaluablesourceofinformationisintherecordsofpatientswho
were infected with the HIV virus through blood transfusion, of which the date can be ascertained
retrospectively. A data set collected by the Centers for Disease Control and Prevention (CDC) is
listed in Wang (1989), which includes the time X from the transfusion to the diagnosis of AIDS,
and the time Y from the transfusion to the end of study (July 1986), both in months, for 295
individuals. It is clear that X (cid:20) Y, i.e., the data are truncated.
With the data in an R data frame aids with components incu and infe, a tensor product
cubic spline model can be (cid:12)tted using the following commands.
aids.fit <- ssden(~incu*infe, data=aids, quad=list(pt=quad.pt,wt=quad.wt),
domain=data.frame(incu=c(0,100),infe=c(0,100)), seed=2375)
The seed option was needed for the results to be reproducible; note that z in (6) is a random
j
subset. The quadrature in quad.pt and quad.wt is a 40(cid:2)40 uniform grid truncated to fx (cid:20) yg,
with half weights on the diagonal; the code for the generation of quad.pt and quad.wt can be
found in the documentation of ssden.
The model has three terms, "incu", "infe", and "incu:infe". To \test" for the practical
signi(cid:12)cance of the term "incu:infe", use
project(aids.fit, include=c("incu","infe"))
8
which yields 1 − (cid:26) = KL((cid:17)^;(cid:17)~)=KL((cid:17)^;(cid:17) ) = 0:0144, suggesting the adequacy of pre-truncation
c
independence,
4.2 Progression of Diabetic Retinopathy
The Wisconsin Epidemiological Study of Diabetic Retinopathy (WESDR) is an on-going epidemi-
ological study of a cohort of patients receiving their medical care in an 11-county area in southern
Wisconsin, who were (cid:12)rst examined in 1980-1982, then again in 1984-1986, 1990-1992, and 1994-
1996. AsubsetderivedfromtheWESDRdataisdistributedinGRKPACK(Wang1997),consisting
of the baseline measures of duration of diabetes in years, percent of glycosylated hemoglobin, body
mass index, and a binary indicator of retinopathy progression at the (cid:12)rst follow-up, of 669 patients.
With the data in an R data frame wesdr with components dur, gly, bmi, and ret, a tensor
product linear spline model can be (cid:12)tted and the additivity \tested."
wesdr.fit0 <- gssanova1(ret~dur*gly*bmi, family="binomial", type="linear",
data=wesdr, nbasis=135, seed=2375)
project(wesdr.fit0, include=c("dur","gly","bmi"))
Tensor product linear splines are used so the number of smoothing parameters (cid:18) is kept man-
(cid:12)
ageable, and nbasis is speci(cid:12)ed to override the default q = 10n2=9 set for cubic splines. The
Kullback-Leibler projection yields 1 − (cid:26) = KL((cid:17)^;(cid:17)~)=KL((cid:17)^;(cid:17) ) = 0:0000, so an additive model
c
appears adequate. A cubic spline additive model can then be (cid:12)tted and evaluated.
wesdr.fit <- gssanova1(ret~dur+bmi+gly, "binomial", , wesdr, seed=2375)
project(wesdr.fit, include=c("gly","bmi"))
project(wesdr.fit, include=c("dur","bmi"))
project(wesdr.fit, include=c("dur","gly"))
The Kullback-Leibler projections yield KL((cid:17)^;(cid:17)~)=KL((cid:17)^;(cid:17) ) = 0:0808, 0:7354, 0:0781; all the terms
c
appear practically signi(cid:12)cant.
4.3 Survival After Heart Transplant
One of the most demonstrated survival data is the Stanford heart transplant data. The data listed
in Miller and Halpern (1982) contain survival or censoring times of 184 patients after their ((cid:12)rst)
heart transplant, in days, and their ages at transplant. To spread out the observations for adequate
nonparametric modeling, a square root transform shall be applied to the time axis.
With the data in an R data frame stan with components futime, status, and age, a hazard
model can be (cid:12)tted and evaluated as follows.
9
stan.fit <- sshzd(Surv(futime,status)~futime*age, data=stan, seed=2375)
project(stan.fit, include=c("futime","age"))
The Kullback-Leibler projection yields KL((cid:17)^;(cid:17)~)=KL((cid:17)^;(cid:17) ) = 0:1245, indicating that a proportional
c
hazard model may not be adequate. The time axis transformation does not a(cid:11)ect hazard propor-
tionality or the lack of it.
4.4 Treatments of Gastric Cancer
Thesurvivaltimesof90gastriccancerpatientsarelistedinMoreau,O’Quigley,andMesbah(1985).
Half of the patients were treated by chemotherapy, the other half by chemotherapy combined with
radiotherapy.
WiththedatainanRdataframegastricwithcomponentsfutime,status,andtrt,ahazard
model can be (cid:12)tted and evaluated as follows; trt is a factor variable.
gastric.fit <- sshzd(Surv(futime,status)~futime*trt, data=gastric, seed=2375)
project(stan.fit, include=c("futime","trt"))
The Kullback-Leibler projection yields KL((cid:17)^;(cid:17)~)=KL((cid:17)^;(cid:17) ) = 0:7307, suggesting that hazard pro-
c
portionality is strongly violated.
5 Simulations
Togainsomeperspectiveconcerningthecalibrationofthediagnostics,simplesimulationsoflimited
scales are presented in this section. The purpose of the simulations is to seek practical insights,
but not quite to establish universal decision rules.
5.1 Density Estimation
Samples X , i = 1;:::;n, were generated from bivariate normal distributions with mean 0, unit
i
variance, and correlation r. For each of the six combinations of sample size n = 100;200 and
correlation r = 0;:3;:5, twenty replicates were generated. For each replicate, the domain X is
rectangular with the marginals extending 5% on each end beyond the data ranges. The domain is
mapped to [0;1]2 and the tensor product cubic spline of x2.2 is used in (3). Five di(cid:11)erent random
subsets fz g (cid:26) fX g were used to calculate (cid:12)ve di(cid:11)erent (cid:12)ts with interaction, and Kullback-Leibler
j i
projectionsintotheadditivemodelspacewerecalculated. Thesizesoffz gwereq = 30forn = 100
j
and q = 33 for n = 200.
The simulation results are summarized in Figure 1. The left frames illustrate the boxplots of
KL((cid:17)^;(cid:17)~)=KL((cid:17)^;(cid:17) ) (of 5 (cid:12)ts) for each of the 6(cid:2)20 replicates. The right frames plot the median
c
10
Description:Model Diagnostics for Smoothing Spline. ANOVA Models. Chong Gu∗. Abstract. Functional ANOVA decompositions can be incorporated in multivariate function
