# Robust Medical Test Evaluation Using Flexible Bayesian Semiparametric Regression Models

Authors: Adam J. Branscum; Wesley O. Johnson; Andre T. Baron
Journal: Epidemiology Research International
 (2013)

---

## Abstract

The application of Bayesian methods is increasing in modern epidemiology. Although parametric Bayesian analysis has penetrated the population health sciences, flexible nonparametric Bayesian methods have received less attention. A goal in nonparametric Bayesian analysis is to estimate unknown functions (e.g., density or distribution functions) rather than scalar parameters (e.g., means or proportions). For instance, ROC curves are obtained from the distribution functions corresponding to continuous biomarker data taken from healthy and diseased populations. Standard parametric approaches to Bayesian analysis involve distributions with a small number of parameters, where the prior specification is relatively straight forward. In the nonparametric Bayesian case, the prior is placed on an infinite dimensional space of all distributions, which requires special methods. A popular approach to nonparametric Bayesian analysis that involves Polya tree prior distributions is described. We provide example code to illustrate how models that contain Polya tree priors can be fit using SAS software. The methods are used to evaluate the covariate-specific accuracy of the biomarker, soluble epidermal growth factor receptor, for discerning lung cancer cases from controls using a flexible ROC regression modeling framework. The application highlights the usefulness of flexible models over a standard parametric method for estimating ROC curves.


---

## Body

## 1. Introduction

Bayesian analysis is often used in the support of epidemiologic research [1–7]. A contemporary area of research in the population health sciences involves the development and application of statistical methods for evaluating the accuracy of medical tests. With binary outcome data, statistical methods focus on estimating sensitivity and specificity, while with quantitative data, standard objects of interest are the receiver operating characteristic (ROC) curve and area under the curve (AUC).The ROC curve can be regarded as a graphical portrayal of the degree of separation between the distributions of test outcomes for “diseased” and nondiseased populations. The formula for an ROC curve depends onSe(k) and Sp(k), the sensitivity and specificity of the test at classification threshold k. We proceed under the innocuous assumption that test outcomes tend to be higher for individuals in the diseased population. Let y denote a continuous test outcome for a disease D, where disease status is labeled as D=0 for disease absent and D=1 for disease present. In general, y can be any continuously measured classifier that varies according to a cumulative distribution function G0 for nondiseased individuals and varies according to G1 for diseased individuals. The ROC curve is a plot of true positive probability versus false positive probability across all classification thresholds, that is, a plot of pairs (1-Sp(k),Se(k)) for all k. Since Se(k)=Pr(y>k∣D=1)=1-G1(k) and 1-Sp(k)=Pr(y>k∣D=0)=1-G0(k), accurate direct estimation of the ROC curve depends crucially on accurate estimation of the distribution functions G0 and G1.This motivates the use of flexible data-driven methods for estimating functions, which is a central goal of nonparametric Bayesian analysis. The standard parametric approach imposes the strong condition thatG0 and G1 be members of the normal family (i.e., the normal-normal model), while a flexible nonparametric approach treats them as arbitrary continuous distribution functions. A flexible method is desirable because G0 and/or G1 could exhibit unanticipated right or left skewness, multimodality, or they could be symmetric but not bell-shaped. As a result, the methodology we present has the flexibility to handle data from many different settings, and notably it will be designed to include the normal-normal model as a special case.The development of modern nonparametric and semiparametric Bayesian procedures for ROC analysis is an active area of research [8–17]. The main goals of this paper are (i) to use flexible Bayesian models and ROC curves to evaluate the accuracy of a biomarker that may depend on patient characteristics and (ii) to provide a nontechnical description of Polya tree priors and show how they can be implemented in epidemiological research using standard statistical software.


## 2. Materials and Methods

Robust inference for ROC curves stems from allowingG0 and G1 to be members of broad classes of distributions. Our approach is to model them using Mixtures of Finite Polya Trees (MFPT) priors, which have been carefully discussed in the statistics literature [18, 19]. We avoid a full technical description because it requires considerable mathematical detail. Instead, a conceptual introduction appears in Appendix A.Briefly, a major advantage beyond flexibility is that MFPT priors forG0 and G1 can be “centered” on separate parametric families. In part, this means that the flexible model generalizes the standard parametric model since the normal-normal model is allowed as a special case. Additionally, it means that the expected value of Gℓ under the prior will be the distribution function that corresponds to the centering family. In our illustrations, the centering families will be Fℓ≡{N(μℓ,σℓ2):μℓ∈ℝ,σℓ2∈ℝ+},ℓ=0,1. The MFPT priors will also have positive weight parameters, c0 and c1, and positive integers J0 and J1 that define the (finite) length of the trees. Large values for the weights indicate higher prior confidence in the centering distributions, while values near zero allow for considerable flexibility around the centering family. To aid the selection of a weight parameter, consider that the parametric model Fℓ is the special case obtained when cℓ→∞. Hence, using a large weight value is similar to fitting the parametric centering model to the data (possible oversmoothing). The opposite extreme of setting cℓ≐0 produces empirical estimates (possible undersmoothing). In our experience, a weight parameter between 0.5 and 1 is often a positive compromise between these extremes. Tree lengths generally range between 4 and 10, depending on the sample size. Small values for Jℓ tend to oversmooth the data. A fully nonparametric prior is obtained as Jℓ→∞, but large values of Jℓ can substantially increase computing time. Finally, there will be priors on (μℓ,σℓ); write them as pℓ(μℓ,σℓ), for ℓ=0,1. We assume G0 and G1 are independent, with priors denoted by Gℓ~MFPT(Fℓ,pℓ,cℓ), for ℓ=0,1.


### 2.1. Models

In the absence of covariates, consider two independent samples of continuous biomarker measurements, wherey1,y2,…,yn0 constitute a random sample from an unknown distribution G0, and z1,z2,…,zn1 are realizations from G1. We regard the yi’s as biomarker outcomes from n0 nondiseased individuals and the zj’s as outcomes from n1 diseased individuals. Our model is represented as yi∣G0~indepG0 with G0~MFPT(F0,p0,c0) and similarly for the zj’s. The prior pℓ(μℓ,σℓ) in our illustrations involves a normal distribution on μℓ and an independent uniform distribution on σℓ.There are many possible extensions of this two-group model depending on the complexity of the data. We describe a semiparametric regression model that can be easily adapted to handle a variety of scenarios. The model specifies separate linear regressions with arbitrary residual distributions for the data from the nondiseased and diseased populations:(1)yi=β1xi1+β2xi2+⋯+βgxig+ϵi≡xiβ+ϵi,ϵi∣G0~G0,G0~MFPT(N(μ0,σ02),p0,c0);zj=α1x~j1+α2x~j2+⋯+αhx~jh+ϵ~j≡x~jα+ϵ~j,ϵ~j∣G1~G1,G1~MFPT(N(μ1,σ12),p1,c1).There are no intercepts in the regression part of the models since each residual distribution has an arbitrary unknown mean. The covariate vectorsxi and x~j can be distinct or have some overlap (e.g., both might include age) or complete overlap (xi=x~j). In the absence of available prior information, independent normal prior distributions with mean 0 and large variance can be used for the regression coefficients. Alternatively, when prior information is available, we recommend conditional means priors [20]. Special cases include a model with no covariates for the nondiseased population (no xir variables), a model with no covariates for the diseased population (no x~js variables) and a model without covariates (the two-group model).Our model can be used to estimate covariate-specific ROC curves and AUC. Letx*=(x,x~) define two particular subpopulations of nondiseased and diseased people. Then, ROC(x*) and its corresponding AUC(x*) are measures of the tests’ accuracy when applied to nondiseased people with covariate vector x and diseased people with covariate vector x~ (often, the researcher will set x equal to x~). Approximate posterior inference for ROC curves can be made by calculating posterior summaries of (1-G0(k-xβ),1-G1(k-x~α)) across a fine grid of values for k that spans the sample space.


## 2.1. Models

In the absence of covariates, consider two independent samples of continuous biomarker measurements, wherey1,y2,…,yn0 constitute a random sample from an unknown distribution G0, and z1,z2,…,zn1 are realizations from G1. We regard the yi’s as biomarker outcomes from n0 nondiseased individuals and the zj’s as outcomes from n1 diseased individuals. Our model is represented as yi∣G0~indepG0 with G0~MFPT(F0,p0,c0) and similarly for the zj’s. The prior pℓ(μℓ,σℓ) in our illustrations involves a normal distribution on μℓ and an independent uniform distribution on σℓ.There are many possible extensions of this two-group model depending on the complexity of the data. We describe a semiparametric regression model that can be easily adapted to handle a variety of scenarios. The model specifies separate linear regressions with arbitrary residual distributions for the data from the nondiseased and diseased populations:(1)yi=β1xi1+β2xi2+⋯+βgxig+ϵi≡xiβ+ϵi,ϵi∣G0~G0,G0~MFPT(N(μ0,σ02),p0,c0);zj=α1x~j1+α2x~j2+⋯+αhx~jh+ϵ~j≡x~jα+ϵ~j,ϵ~j∣G1~G1,G1~MFPT(N(μ1,σ12),p1,c1).There are no intercepts in the regression part of the models since each residual distribution has an arbitrary unknown mean. The covariate vectorsxi and x~j can be distinct or have some overlap (e.g., both might include age) or complete overlap (xi=x~j). In the absence of available prior information, independent normal prior distributions with mean 0 and large variance can be used for the regression coefficients. Alternatively, when prior information is available, we recommend conditional means priors [20]. Special cases include a model with no covariates for the nondiseased population (no xir variables), a model with no covariates for the diseased population (no x~js variables) and a model without covariates (the two-group model).Our model can be used to estimate covariate-specific ROC curves and AUC. Letx*=(x,x~) define two particular subpopulations of nondiseased and diseased people. Then, ROC(x*) and its corresponding AUC(x*) are measures of the tests’ accuracy when applied to nondiseased people with covariate vector x and diseased people with covariate vector x~ (often, the researcher will set x equal to x~). Approximate posterior inference for ROC curves can be made by calculating posterior summaries of (1-G0(k-xβ),1-G1(k-x~α)) across a fine grid of values for k that spans the sample space.


## 3. Results

### 3.1. Fitting MFPT Models in SAS

The SAS 9.3 (SAS Institute, Inc., Cary, North Carolina) code that was used to fit a one-sample MFPT model to the simulated data described in AppendixA is in Appendix B.


### 3.2.  ROC Analysis of sEGFR Data

We investigated a soluble isoform of the epidermal growth factor receptor (sEGFR) as a biomarker for lung cancer in premenopausal and postmenopausal women. sEGFR also has been associated with ovarian cancer [21, 22]. The data were collected from a case-control study conducted at the Mayo Clinic in Rochester, Minnesota between 1998 and 2003. There were 140 controls and 101 lung cancer cases, and 92 premenopausal and 149 postmenopausal women enrolled in the study.In a preliminary analysis, we found no clear evidence of a difference in the distributions ofy=sEGFR for premenopausal versus postmenopausal lung cancer cases (Wilcoxon P=0.15). However, there was evidence for a statistical difference in the distribution of y for controls based on menopausal status (Wilcoxon P<0.01). We thus included menopausal status as a covariate in models for the control data. A flexible approach is supported over the standard parametric approach because normal quantile plots (not shown) indicated a right skew in the distributions of y for cases and postmenopausal controls.We compared two analyses. The first analysis had menopausal status as a covariate for the control group in an ROC regression model. The second analysis modeled outcomes among pre- and postmenopausal controls with completely separate (flexible) distributions, without regression structure. We compared these models using the log pseudomarginal likelihood (LPML) [23], [24, Section 4.9].For all models that used Polya trees, we set weight parameters (thecℓ’s) equal to 0.5 and tree lengths (the Jℓ’s) equal to 5. In this setting, values of y tended to be higher among controls instead of cases, so we reversed the roles for cases and controls. Let x~=1 for premenopausal women and x~=0 for postmenopausal women. For both models under consideration, the data from cases were modeled as independently distributed according to an unknown distribution G0. The prior was G0~MFPT(N(μ0,σ02),p0,c), with p0(μ0,σ0) selected to be μ0~N(50,400) independent of σ0~Uniform(0,100). The prior mean of 50 for μ0 was chosen because that was approximately the mean of y for male lung cancer cases in a similar study, but we allowed for a high degree of uncertainty about the value of μ0 through a large prior variance.For controls, model 1 for the data was the regression(2)zj=α1x~j+ϵ~j,ϵ~j∣G1~G1,G1~MFPT(N(μ1,σ12),p1,c),
where, independently, α1~N(0,400), μ1~N(50,400), and σ1~Uniform(0,100).Model 2 specified separate distributions for data from premenopausal and postmenopausal controls. Denote the distributions byG2 and G3, respectively, where Gw was assigned an MFPT prior that was centered at N(μw,σw2), with μw~N(50,400) and σw~Uniform(0,100), for w=2,3. Using the same priors, we also considered the standard parametric analysis (model 3) in which G0, G2, and G3 were normal distributions with distinct means and variances.The LPML statistics were similar for models 1 and 2 (LPML≐-1187). We therefore selected the more flexible model 2 that had separate distributions across x~. For the standard parametric model 3, LPML=-1201, indicating a strong preference for the more flexible models (pseudo-Bayes factor =e14).Estimated ROC curves were obtained from model 2, together with estimates from the parametric normal model 3 and from nonparametric empirical distribution functions. The estimated ROC curve from the parametric analysis differs drastically from estimates obtained from the MFPT and empirical distribution functions for postmenopausal women (Figure1). Not surprisingly, the parametric model fails to track either of the data-driven estimates. Note that although the parametric model is its prior expectation, the MFPT model has the flexibility to allow the data to reallocate probability mass away from a bell-curve to produce an appropriately smooth and flexible estimate of the ROC curve. For premenopausal women, the parametric and Polya tree models also give different inferences, including at clinically important cutoffs corresponding to false positive probabilities up to 10% (Figure 2). The empirical partial AUC over that region was 0.026, while the estimate from the parametric analysis was increased by 37% to 0.041. A 90% posterior interval for the partial AUC from the parametric model was (0.027, 0.055), which does not contain the empirical estimate. The Polya tree analysis is a middle ground between these extremes, with an estimated ROC curve that is essentially a smoothed version of the empirical curve and an estimated partial AUC of 0.019 (0.009, 0.036). Similar results were obtained from a sensitivity analysis that placed diffuse Gamma priors with mean 1 and variance 10000 on the μℓ’s and σℓ’s.Figure 1
Empirical ROC curve for postmenopausal women and estimates of the ROC curve from a parametric normal model (dashed line) and a mixture of finite Polya trees analysis (smooth solid line).Figure 2
Empirical ROC curve for premenopausal women and estimates of the ROC curve from a parametric normal model (dashed line) and a mixture of finite Polya trees analysis (smooth solid line).


## 3.1. Fitting MFPT Models in SAS

The SAS 9.3 (SAS Institute, Inc., Cary, North Carolina) code that was used to fit a one-sample MFPT model to the simulated data described in AppendixA is in Appendix B.


## 3.2.  ROC Analysis of sEGFR Data

We investigated a soluble isoform of the epidermal growth factor receptor (sEGFR) as a biomarker for lung cancer in premenopausal and postmenopausal women. sEGFR also has been associated with ovarian cancer [21, 22]. The data were collected from a case-control study conducted at the Mayo Clinic in Rochester, Minnesota between 1998 and 2003. There were 140 controls and 101 lung cancer cases, and 92 premenopausal and 149 postmenopausal women enrolled in the study.In a preliminary analysis, we found no clear evidence of a difference in the distributions ofy=sEGFR for premenopausal versus postmenopausal lung cancer cases (Wilcoxon P=0.15). However, there was evidence for a statistical difference in the distribution of y for controls based on menopausal status (Wilcoxon P<0.01). We thus included menopausal status as a covariate in models for the control data. A flexible approach is supported over the standard parametric approach because normal quantile plots (not shown) indicated a right skew in the distributions of y for cases and postmenopausal controls.We compared two analyses. The first analysis had menopausal status as a covariate for the control group in an ROC regression model. The second analysis modeled outcomes among pre- and postmenopausal controls with completely separate (flexible) distributions, without regression structure. We compared these models using the log pseudomarginal likelihood (LPML) [23], [24, Section 4.9].For all models that used Polya trees, we set weight parameters (thecℓ’s) equal to 0.5 and tree lengths (the Jℓ’s) equal to 5. In this setting, values of y tended to be higher among controls instead of cases, so we reversed the roles for cases and controls. Let x~=1 for premenopausal women and x~=0 for postmenopausal women. For both models under consideration, the data from cases were modeled as independently distributed according to an unknown distribution G0. The prior was G0~MFPT(N(μ0,σ02),p0,c), with p0(μ0,σ0) selected to be μ0~N(50,400) independent of σ0~Uniform(0,100). The prior mean of 50 for μ0 was chosen because that was approximately the mean of y for male lung cancer cases in a similar study, but we allowed for a high degree of uncertainty about the value of μ0 through a large prior variance.For controls, model 1 for the data was the regression(2)zj=α1x~j+ϵ~j,ϵ~j∣G1~G1,G1~MFPT(N(μ1,σ12),p1,c),
where, independently, α1~N(0,400), μ1~N(50,400), and σ1~Uniform(0,100).Model 2 specified separate distributions for data from premenopausal and postmenopausal controls. Denote the distributions byG2 and G3, respectively, where Gw was assigned an MFPT prior that was centered at N(μw,σw2), with μw~N(50,400) and σw~Uniform(0,100), for w=2,3. Using the same priors, we also considered the standard parametric analysis (model 3) in which G0, G2, and G3 were normal distributions with distinct means and variances.The LPML statistics were similar for models 1 and 2 (LPML≐-1187). We therefore selected the more flexible model 2 that had separate distributions across x~. For the standard parametric model 3, LPML=-1201, indicating a strong preference for the more flexible models (pseudo-Bayes factor =e14).Estimated ROC curves were obtained from model 2, together with estimates from the parametric normal model 3 and from nonparametric empirical distribution functions. The estimated ROC curve from the parametric analysis differs drastically from estimates obtained from the MFPT and empirical distribution functions for postmenopausal women (Figure1). Not surprisingly, the parametric model fails to track either of the data-driven estimates. Note that although the parametric model is its prior expectation, the MFPT model has the flexibility to allow the data to reallocate probability mass away from a bell-curve to produce an appropriately smooth and flexible estimate of the ROC curve. For premenopausal women, the parametric and Polya tree models also give different inferences, including at clinically important cutoffs corresponding to false positive probabilities up to 10% (Figure 2). The empirical partial AUC over that region was 0.026, while the estimate from the parametric analysis was increased by 37% to 0.041. A 90% posterior interval for the partial AUC from the parametric model was (0.027, 0.055), which does not contain the empirical estimate. The Polya tree analysis is a middle ground between these extremes, with an estimated ROC curve that is essentially a smoothed version of the empirical curve and an estimated partial AUC of 0.019 (0.009, 0.036). Similar results were obtained from a sensitivity analysis that placed diffuse Gamma priors with mean 1 and variance 10000 on the μℓ’s and σℓ’s.Figure 1
Empirical ROC curve for postmenopausal women and estimates of the ROC curve from a parametric normal model (dashed line) and a mixture of finite Polya trees analysis (smooth solid line).Figure 2
Empirical ROC curve for premenopausal women and estimates of the ROC curve from a parametric normal model (dashed line) and a mixture of finite Polya trees analysis (smooth solid line).


## 4. Conclusions

We have described the use of MFPT priors in Bayesian analysis. Even when a parametric model is thought to hold, a Polya tree analysis offers a way to assess parametric assumptions and to perform a sensitivity analysis to address deviations from them. For example, a simple approach to sensitivity analysis involves comparing estimated ROC curves and AUC from the normal-normal and semiparametric models.A finite Polya tree prior is not technically nonparametric. Bayesian nonparametric models are characterized by having an infinite number of parameters. Such models often use finite approximations, as in our case. It is thus more accurate to view models that contain finite Polya tree priors as parametric, with a large number of parameters. In fact, Bayesian statistical analysis with finite Polya tree priors has been called “parametric nonparametric statistics,” because it uses parametric models (with many parameters) that maintain flexibility [25]. In our experience, finite approximations of Polya trees are sufficiently flexible and scalable to perform well for a wide variety of complex settings.



---

*Source: 131232-2013-12-11.xml*