We propose an automatic Bayesian approach to the selection of covariates and penalised splines transformations thereof in generalised additive models. Specification of a hyper-g prior for the model parameters and a multiplicity-correction prior for the models themselves is crucial for this task. We introduce the methodology in the normal model and illustrate it with an application to diabetes data. Extension to non-normal exponential families is finally discussed

Held, L

Kauermann, G

Sabanés Bové, D

ZORA

Student Oral PresentationHyper-g Priors for Generalised AdditiveModel SelectionDaniel Sabane´s Bove´1, Leonhard Held1 and Go¨ranKauermann21 Division of Biostatistics, Institute of Social and Preventive Medicine, Uni-versity of Zurich, Hirschengraben 84, CH-8001 Zurich, Switzerland, Email:daniel.sabanesbove@ifspm.uzh.ch, leonhard.held@ifspm.uzh.ch2 Centre for Statistics, Department of Economics and Business Administration,University Bielefeld, Postfach 300131, D-33501 Bielefeld, Germany, Email:gkauermann@uni-bielefeld.deAbstract: We propose an automatic Bayesian approach to the selection of covari-ates and penalised splines transformations thereof in generalised additive models.Specification of a hyper-g prior for the model parameters and a multiplicity-correction prior for the models themselves is crucial for this task. We introducethe methodology in the normal model and illustrate it with an application todiabetes data. Extension to non-normal exponential families is finally discussed.Keywords: Penalised splines; Bayesian variable selection; Shrinkage.1 IntroductionSuppose we have pmetrical covariates x1, . . . , xp and use the additive modely = β0 +∑pj=1mj(xj) + , where  ∼ N(0, σ2). When xj is included non-linearly in the model, we assumemj(xj) = xjβj +Zj(xj)Tujwhere Zj(xj) is the K × 1 spline basis vector at position xj and uj ∼N(0, σ2ρjI) is the corresponding coefficients vector. In order to combinen observations, we stack these to the n× 1 vector xj and the n×K basismatrix Zj , both modified to be zero-centred and orthogonal to each other.We then translate the variance parameter ρj into the corresponding degreeof freedom (Aerts, Claeskens and Wand, 2002, section 2.2)dj(ρj) = tr{(ZTj Zj + ρ−1j I)−1ZTj Zj}+ 1 ∈ (1,K + 1). (1)A larger ρj (or a larger dj) leads to a weaker penalty on the non-linearcomponent of the function mj . If xj is excluded from or linearly includedin the model we have mj(xj) ≡ 0 or mj(xj) = xjβj and set dj = 0 ordj = 1, respectively. Thus, the function mj is exactly defined by dj , whichwe may restrict to a finite set of values, say dj ∈ {0, 1} ∪ {2, 3, . . . ,K}.2 Hyper-g Priors for Generalised Additive Model SelectionAs default prior for the parameters β0, β = (βj : dj ≥ 1) and σ2 in a givenmodel specified via d = (d1, . . . , dp),y |β0,β,u, σ2 ∼ N(1β0 +Xβ +Zu, σ2I) (2)with X = (xj : dj ≥ 1), Z = (Zj : dj > 1) and u = (uTj : dj > 1)T ,we propose the hyper-g prior (Liang et al., 2008) described in Section 2.For the models we propose a multiplicity-correction prior in Section 3.The methodology is applied to diabetes data in Section 4 and extended togeneralised additive models in Section 5.2 Hyper-g Priors for Additive ModelsIntegrating out the spline coefficients vector u ∼ N(0, σ2D), where D =diag{ρjI : dj > 1}, from the conditional model (2) yields the marginalmodely |β0,β, σ2 ∼ N(1β0 +Xβ, σ2V)(3)with V = I + ZDZT having Cholesky decomposition V = RTR. Thetransformed response vector y˜ = R−Ty follows a linear model with sim-ilarly transformed design matrix X˜ and diagonal covariance matrix σ2I.It turns out that we can use the hyper-g prior (Liang et al., 2008) for thistransformed model, i. e. a locally uniform prior p(β0) ∝ 1 on the intercept,Jeffreys’ prior p(σ2) ∝ (σ2)−1 on the variance and the g-prior (Zellner,1986)β | g, σ2 ∼ N(0, gσ2(X˜TX˜)−1)(4)on the coefficients are combined with a uniform prior on the shrinkagecoefficient g/(1 + g). Note that σ−2X˜TX˜ = σ−2XTV −1X is the Fisherinformation matrix of β in the marginal model (3). The hyper-g prior leadsto a closed form of the marginal likelihood, which we can compute on theoriginal response scale via the change of variables formula:p(y |d) ∝ ‖y˜ − ˜¯y‖−(n−1) (ld + 2)−12F1(n− 12; 1;ld + 42; R˜2)|R|−1,where ld is the dimension of β, 2F1 is the Gaussian hypergeometric functionand R˜2 is the classical coefficient of determination in model (3).3 Model PriorWe propose a prior p(d) on the model space which explicitly corrects forthe multiplicity of testing inherent in the simultaneous analysis of manycovariates (see Scott and Berger, 2010): A priori, the number of covari-ates included in the model (ld) is uniformly distributed on {0, 1, . . . , p}.Daniel Sabane´s Bove´ et al. 3TABLE 1. Marginal posterior probabilities (x1: age, x2: systolic blood pressure,x3: cholesterol/HLD ratio, x4: BMI, x5: waist/hip ratio, x6: gender).x1 x2 x3 x4 x5 x6not included (dj = 0) 0.00 0.65 0.00 0.14 0.50 0.65linear (dj = 1) 0.71 0.33 0.93 0.81 0.48 0.35non-linear (dj > 1) 0.29 0.03 0.07 0.05 0.02 —Then the number of non-linearly included covariates (sd) is uniformly dis-tributed on {0, 1, . . . , ld}. The respective choice of the ld and sd covariatesis uniformly distributed on all possible configurations. Finally, the degreesof freedom of the non-linearly modelled covariates are independent anduniformly distributed on {2, 3, . . . ,K}. Altogether, this gives1/p(d) =(pld)(p+ 1)(ldsd)(ld + 1)(K − 1)sdand leads to marginal prior probabilities Pr(dj = 0) = 1/2, Pr(dj = 1) =Pr(dj > 1) = 1/4.4 ApplicationWe illustrate our modelling approach with the diabetes data from Harrell(2001). We study the association of (the negative reciprocal of) glycosolatedhaemoglobin of n = 377 study participants with the continuous covariatesage (in years), systolic blood pressure (in mmHg), cholesterol/HDL ratio,body mass index (BMI, in kg/m2) and waist/hip ratio as well as the bi-nary covariate gender. As the computational complexity is quadratic inthe spline basis dimension K, we want to use splines with few quantile-based knots. Therefore, we choose cubic O’Sullivan splines (Wand andOrmerod, 2008). Here, we get basis matrices Zj with K = 9 columnsfrom 7 knots. The exhaustive evaluation of the posterior model probabili-ties p(d |y) ∝ p(y |d)p(d) of all (K + 1)5 · 2 = 200 000 models takes only585 seconds due to an efficient C++ implementation which is available in anR-package from the first author. In Table 1 the marginal posterior proba-bilities for linear and non-linear inclusion of the six covariates are shown.There is strong evidence for linear inclusion of cholesterol/HDL ratio andBMI, while the posterior probability for inclusion of systolic blood pressureor gender is only 35%. There is overwhelming evidence for (non-linear) in-clusion of age, and the posterior odds for (linear) inclusion of waist/hipratio are around 1. The maximum a posteriori model includes age, choles-terol/HDL ratio and BMI all linearly. Note that these are the covariateswhich have inclusion probabilities larger than 50%, thus defining the set ofmedian probability models (Barbieri and Berger, 2004) d with d1, d3, d4 ≥ 14 Hyper-g Priors for Generalised Additive Model Selection20 40 60 80−0.04−0.020.000.020.040.06x1 (age in years)m1(x1)5 10 15 200.000.050.100.15x3 (cholesterol/HDL ratio)m3(x3)20 30 40 50−0.020.000.020.04x4 (BMI in kg/m2)m4(x4)FIGURE 1. Estimated covariate effects in the median probability model average,based on 10 000 samples: Posterior means (solid lines), pointwise (dashed lines)and simultaneous (dotted lines) 95%-credible intervals as well as positions of datapoints (ticks above x-axes) are shown.and d2 = d5 = d6 = 0. Figure 1 shows the estimated covariate effects fromthe resulting model average. While the age effect is slightly non-linear (with38% probability in the median probability models), both other covariateshave essentially linear effect estimates.5 Extension to Generalised Additive ModelsNow we assume more generally that the covariate effects mj(xj) enter addi-tively into the linear predictor η = β0+∑pj=1mj(xj) of an exponential fam-ily distribution with canonical parameter θ, mean E(y) = h(η) = db(θ)/dθand variance Var(y) = φ/w · v(µ) = φ/w · d2b(θ)/dθ2 (see McCullaghand Nelder, 1989). We restrict our attention to non-normal distributionswith fixed dispersion φ (as φ = 1 for the Bernoulli and Poisson distribu-tion) and known weight w. For n observations, the linear predictor vectorη = (η1, . . . , ηn)T is η = 1β0 + Xβ + Zu, where u ∼ N(0,D), and thelikelihood isp(y |β0,β,u) ∝ exp{n∑i=1yiθi − b(θi)φ/wi}. (5)A reasonable generalisation of (1) is (see Ruppert, Wand and Carroll, 2009,section 11.4)dj(ρj) = tr{(ZTj ŴZj + ρ−1j I)−1ZTj ŴZj}+ 1, (6)which uses a fixed weight matrix Ŵ = W (1β̂0) for all models, whereW (η) = diag{(dh(ηi)/dη)2v(h(ηi))−1φ−1wi}ni=1 is the usual generalisedlinear model (GLM) weight matrix and β̂0 is the intercept estimate fromthe null model. Therefore, we now arrange 1, xj and the columns of Zj tobe orthogonal with respect to the inner product in terms of Ŵ , so that (6)Daniel Sabane´s Bove´ et al. 5correctly captures the degrees of freedom associated with the non-linearpart of mj .In order to derive a generalised g-prior for β, we will use the iterativeweighted least squares (IWLS) approximation to (5) to come back to anormal model and then derive the resulting g-prior (4). So letz0 = η0 + diag{dh(η0)/dη}−1(y − h(η0))be the adjusted response vector resulting from a first-order approximationto h−1(y) around h(η0), such thatz0 |β0,β,u ∼ N(1β0 +Xβ +Zu, W (η0)−1)is the working normal model. This can be rewritten toz˜0 |β0,β,u ∼ N(1˜β0 + X˜β + Z˜u, I) (7)by setting z˜0 = W (η0)1/2z0 etc. Since (7) is analogous to (2), our proposalfor a generalised g-prior isβ | g ∼ N(0, gJ−1), (8)where J is the Fisher information for β in (7) with η0 = 0:J = X˜T(I + Z˜DZ˜T)−1X˜= XTW1/20 (I +W1/20 ZDZTW1/20 )−1W 1/20 X,abbreviating W 0 = W (0). Note that this prior directly generalises theprior proposed by Sabane´s Bove´ and Held (2011) for GLMs, to which itreduces when there are no spline effects in the model.The generalised hyper-g prior then consists of the improper prior p(β0) ∝ 1on the intercept β0, the g-prior (8) on the linear effects vector β, thepenalty prior u ∼ N(0,D) on the spline coefficients vector u and someproper hyper-prior p(g) on the hyper-parameter g in the g-prior. For theimplementation of posterior inference we can easily extend the approachof Sabane´s Bove´ and Held (2011, section 3). Let Xa = (1,X,Z) andβa = (β0,βT ,uT )T , such that η = Xaβa. The prior for βa conditionalon g has Gaussian form with mean zero and singular precision Ra =diag{0, g−1J(0),D−1}. Thus, the Laplace approximation of p(y | g,d), whichis based on a Gaussian approximation to the conditional posterior p(βa |y, g),can be obtained by the Bayesian IWLS algorithm (West, 1985). Afterwards,the marginal likelihoodp(y |d) =∫ ∞0p(y | g,d)p(g) dg,can be approximated by numerical integration of the Laplace approxima-tion p˜(y | g,d). Note that this strategy of integrated Laplace approxima-tions was proposed more generally by Rue, Martino and Chopin (2009).Finally, for sampling from the posterior of βa and g in a specific model dwe can use a tuning-free Metropolis-Hastings algorithm.6 Hyper-g Priors for Generalised Additive Model SelectionReferencesAerts, M., Claeskens, G., and Wand, M.P. (2002). Some theory for penal-ized spline generalized additive models. Journal of Statistical Plan-ning and Inference, 103, 455-470.Barbieri, M.M., and Berger, J.O. (2004). Optimal predictive model selec-tion. Annals of Statistics, 32, 870-897.Harrell, Jr., F.E. (2001). Regression Modeling Strategies. New York: Springer.Liang, F., Paulo, R., Molina, G., Clyde, M.A., and Berger, J.O. (2008). Mix-tures of g priors for Bayesian variable selection. Journal of the Amer-ican Statistical Association, 103, 410-423.McCullagh, P., and Nelder, J.A. (1989). Generalized Linear Models. NewYork: Chapman and Hall.Rue, H., Martino, S., and Chopin, N. (2009). Approximate Bayesian infer-ence for latent Gaussian models by using integrated nested Laplaceapproximations. Journal of the Royal Statistical Society. Series B(Methodological), 71, 319-392.Ruppert, D., Wand, M.P., and Carroll, R.J. (2003). Semiparametric Regres-sion. Cambridge: Cambridge University Press.Sabane´s Bove´, D., and Held, L. (2011). Hyper-g Priors for Generalized Lin-ear Models. Bayesian Analysis, 6, forthcoming article.URL: http://ba.stat.cmu.edu/abstracts/Sabanes.phpScott, J.G., and Berger, J.O. (2010). Bayes and empirical-Bayes multiplic-ity adjustment in the variable-selection problem. Annals of Statistics,38, 2587-2619.Wand, M.P., and Ormerod, J.T. (2008). On semiparametric regression withO’Sullivan penalized splines. Australian & New Zealand Journal ofStatistics, 50, 179-198.West, M. (1985). Generalized linear models: scale parameters, outlier ac-commodation and prior distributions. In: Bayesian Statistics 2: Pro-ceedings of the Second Valencia International Meeting, 531-558. Am-sterdam: North-Holland.Zellner, A. (1986). On assessing prior distributions and Bayesian regres-sion analysis with g-prior distributions. In: Bayesian Inference andDecision Techniques: Essays in Honor of Bruno de Finetti, 233-243.Amsterdam: North-Holland.

Hyper-g priors for generalised additive model selection

http://www.zora.uzh.ch/53039/1/sabanesbove.etal2011a.pdf

Hyper-g priors for generalised additive model selection

Abstract

Similar works

Full text

Available Versions

ZORA