Model Averaging: Beyond Model Uncertainty in Risk Analysis *

1. Model Averaging: Beyond Model Uncertainty in Risk Analysis * Matthew W. Wheeler NIOSH MWheeler@cdc.gov *The findings and conclusions in this report are those of the authors and do not necessarily represent the views of the National Institute for Occupational Safety and Health

2. Acknowledgements • A. John Bailer Miami University/NIOSH

3. Outline • Introduction/Motivation. • Model Averaging 101. • Validation Simulation Study. • MA Software for dichotomous response.

4. Introduction/Motivation • Model choice is frequently often a point of contention among risk assessors/risk managers. • Frequently multiple models describe the data “equivalently.” • Two model’s estimates of risk, especially at the lower bound, can differ dramatically. • Model uncertainty is inherent in most risk estimation, though practically ignored in most situations.

5. Introduction/Motivation (cont) • Consider the problem of estimating a benchmark dose (BMD) from dichotomous dose response data. • Here we seek to estimate the BMD from a “plausible” model, given experimental data. • In these experiments: • Animals are exposed to some potential hazard. • The adverse response is assumed to be distributed binomially. • Risk (i.e, probability of adverse response) is estimated using regression modeling. • Multiple dose-response models can be used to estimate risk.

6. Common Dose-Response Models Used: • logistic model: (1) • log-logistic model: (2) • gamma: (3) • multistage (4) • probit (5)

7. Common Dose-Response Models Used: • log-probit (6) • quantal-linear (7) • quantal-quadratic (8) • Weibull (9) • where Γ(α)= gamma function evaluated at α, for Ф(x) = CDF N(0,1) and πi = γ when di=0 for models (2) and (7).

8. Benchmark dose estimation • BMD is the dose associated with the a specified increase in response relative to the control response (BMR); e.g.,: dose d such that: BMR=[πd- π0]/[1- π0] or BMR= πd- π0 • The BMR is commonly set at values of 1%, 5%, 10%. • BMDL = 100(1-α)% lower confidence limit on the BMD. NOTE: As πd is dependent on a model thus the BMD is model dependent!!

9. Typical Risk Estimation Process • Given data (in absence of mechanistic information), a typical analyst will: • Estimate the regression coefficients for models (1)-(9). • Estimate the BMD/BMDL given the model. • Pick the “best model.”

10. Example of this: • Consider TiO2 lung tumor data which has been combined from the studies of Heinrich et al. (1995) Muhle et al. (1991) and Lee et al. (1985). • Here the benchmark dose (BMD) as well as its lower bound (BMDL) are estimated at BMRs of 10 and 1%.

11. ** All fits were obtained using the US EPA’s BMDS ‡ Calculated using the number of parameters vs. the number of non-bounded parameters.

13. Model Choice • If we pick best AIC the quantal-quadratic model risk estimates would be chosen. • If we pick the best Pearson Χ2 test statistic the 3-degree multistage model would be chosen. • All of the models are reasonable based on these statistics.

14. Estimates are “reasonably similar” at the 10% estimate. Estimates vary by a factor of 5 for the 1% estimate. If a BMR of 0.1% is used (results not shown) the Ti02 BMD/BMDL estimates differ by a factor of 35. This heterogeneity in model estimates exists even though the model fit statistics are very similar. Model uncertainty (and heated debate) result when any one of the above models is chosen. Model Choice (cont)

15. Model Averaging • A better way would be to find an adequate way to combine all estimates, and thus describe/account for model uncertainty. • Model Averaging (MA) is one such method that may satisfactorily account for model uncertainty. • Instead of focusing on a single model it allows researchers to focus on “plausible behavior.”

16. Model Averaging (cont.) • We can think of any model contributing information (including possible bias) to an analysis. • Picking any one model ignores other plausible information, and possibly introduces bias into the analysis. • Model averaging is a method that attempts to synthesize all of the information available.

17. Model Averaging (cont) • Kang et al. (Regulatory Toxicology and Pharmacology, 2000) and Bailer et al. (Risk Analysis,2005) proposed model averaging for risk assessment. • They used an “Average-BMD” methodology.

18. In this situation the benchmark dose is computed as : Here the BMDk is estimated from one of K models individually The BMDL is computed as: Weights are based upon the formula: Where Ii=AIC, Ii=KIC , or Ii=BIC. Other weights are possible. “Average Dose” MA

19. “Average Dose” MA • This method was found to perform unreliably in many situations in terms of nominal coverage and bias. (Wheeler and Bailer, Environmental and Ecological Statistics, In Press). • The results suggested a different BMD estimation technique may perform better. We call this the “average model” MA.

20. “Average-Model” MA • Given the fits of models (1)-(9) MA: • Calculates the dose-response based upon a weighted average of dose-responses Raftery et al. (1997), Buckland et al. (1997) • Estimates the MA dose-response curve as: • Weights are formed as above.

22. “Average-Model” MA Benchmark Dose • Given this “Average-model,” the benchmark dose is then computed by finding the dose that satisfies the equation BMR = [πMA(d)i- πMA(0)]/[1- πMA(0)]. • The BMDL is computed through a parametric bootstrap. Here the 5th percentile of the bootstrap distribution is used to compute the 95% lower tailed confidence limit estimate on the BMD.

23. “Average-Model” vs. “Average-Dose” • This is substantially different from the “average-dose” method. • “Average-Dose:” Model Fits → Individual BMD estimation → BMD MA estimate • “Average-Model:” Model Fits → MA Model Estimate → MA-BMD estimation

24. TiO2 Analysis Revisited

25. Validation Study • MA seems like a good idea, however we need to know if it works well in practice. • A simulation study was conducted to investigate the behavior of MA.

26. Validation Study • A Monte Carlo simulation study was conducted to investigate the behavior of “average model” model averaging. • 54 true model conditions, using models (1) – (9) were used in the simulation. • Full study described in Wheeler and Bailer (Risk Analysis, in Press)

27. Step 1: Decide upon “TRUTH” Truth was determined by using various response patterns. These were chosen to reflect shallow, medium, and steep curvature as well as low and high background response rates. In all six response patterns were chosen.

28. Step 1 (cont) • Using these 6 dose-response patterns the parameters associated with models (1)-(9) (e.g. γ,α,β,θ,θ2 ) were determined. For example: • These fits were used as the true underlying model in the simulation. • Three points allowed a unique specification of all models. • A variety of dose-response patterns were represented by the 54(6 conditions x 9 models) true dose response curves.

29. Step 2: Given truth set up simulation conditions. • The simulation proceeded by generating hypothetical toxicology experiments with response probability π(d). • With π(d) specified by one of the 54 true dose-response curves estimated in step 1. • These experiments consisted of 4 dose group design with doses of 0, 0.25, 0.50, and 1.0. • n=50 for all dose groups. • 2000 experiments were generated per true dose-response curve. • Bias as well as coverage [i.e., Pr(BMDL ≤ BMDtrue)] was estimated. • Coverage is reported here.

30. Step 3: Gather appropriate statistics. • In each experiment the “average-model” BMD as well as the BMDL was estimated. • BMRs of 1% and 10% were used to estimate the BMD. • Two model spaces for averaging were considered. • One space consisted of three flexible models: the multistage, Weibull and the log-probit model. • The second space had seven models that added the probit, logistic, quantal-linear, and quantal-quadratic to the three model space. • The simulation took approximately 1 CPU year of computation.

31. Coverage • The BMDL was computed [with α = 0.05] and compared to the true benchmark dose at the specified BMR. • Coverage probability [i.e., Pr(BMDL ≤ BMDtrue)] was estimated across 2000 simulations.

32. Coverage BMR = 10%

33. Coverage BMR = 1%

34. Coverage (Summary) • Nominal coverage is reached for most simulation conditions. • MA fails to reach nominal coverage in the quantal-linear and similar cases.

35. Quantal-Linear Problems • The quantal-linear case represents a cause of concern if MA is to be used. • It is important to understand when (and why) it doesn’t work, and compare the method to current practice. • We first compare it to the common practice of picking the “best model.” • The “best-model” is chosen by a chi-squared goodness of fit statistic.

36. Coverage Comparison MA vs. Best Model *True model included in the model set.

37. MA vs. “Best Model” • As a decision rule, in terms of coverage, “Average-model” 3-model averaging performs better than picking the best fitting model. NOTE: The true model was included in the set of models the “best model” was chosen from.

38. Further Analysis of Quantal Linear Model • Despite the above it is disconcerting that a set of flexible models like the 3-model space model average fails to describe quantal-linear behavior. • We look at the sampling distribution of the experiment for the quantal-linear model, in attempt to explain this behavior.

39. Sampling distribution for the quantal-linear model

40. Average fit for 3-model MA models

41. Further Analysis of Quantal Linear Model • The flexibility of the models combined with the sampling distribution introduces bias into the estimation of the dose-response curve. • The bias carries through in BMD estimation. • This also may be the cause of the conservative behavior (i.e. coverage > 99%) seen in the quantal-quadratic case.

42. Bias Corrected Intervals • The evident bias suggests a search for a corrective mechanism. • Percentile based bootstrap confidence intervals are not optimal. • Biased Corrected and Accelerated (BCa) confidence intervals (CI) are preferred (Efron, 1994) • The BCa CIs are computationally more difficult. (i.e., full simulation becomes longer if these were being computed in all conditions) • Only the quantal-linear conditions were analyzed using the BCa confidence interval.

44. BCa CIs fail to be appropriate in shallow dose response cases. • In these cases an infinite BMD can be estimated through the bootstrap re-samples. • This causes a problem for the BCa technique, as it violates its underlying assumptions. • These problems are identifiable through histogram plots of the bootstrap resample.

45. Sampling distribution where BCa BMDLfails to cover the true BMD when the percentile based method does.

46. Simulation Conclusions • Model averaging recovers the central estimate, as well as the lower bound in most cases. (Some bias is however evident) • BMDs and BMDLs with BMRs of 1% are most problematic. • Problems exhibited with true quantal-linear dose-response mechanisms. • BCa confidence intervals can help remediate some of these problems.

47. “Average-model” model averaging is superior to picking the best model. • It recovers the true BMD uncertainty (in terms of correct coverage probabilities) when the true model is not in the MA model space. • The results show MA is not a panacea, it is however a step in the right direction.

49. Model Averaging Software • The results are promising but implementation of this approach is difficult. • The simulation code has been repackaged to allow users to implement dichotomous dose-response model averaging. • This is done in a simple MS Windows command prompt program.

50. Software is described that will implement Model Averaging for Dichotomous dose-response data (affectionately named MADr for short). • The software, as well as the code, is freely available under the GNU public license. • It will be available, in the near future, on NIOSH’s web site. A beta version is distributed with this talk.