UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Small sample improvement over Bayes prediction under model uncertainty 2000

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_2000-566463.pdf [ 4.52MB ]
Metadata
JSON: 1.0089709.json
JSON-LD: 1.0089709+ld.json
RDF/XML (Pretty): 1.0089709.xml
RDF/JSON: 1.0089709+rdf.json
Turtle: 1.0089709+rdf-turtle.txt
N-Triples: 1.0089709+rdf-ntriples.txt
Citation
1.0089709.ris

Full Text

Small Sample Improvement Over Bayes Prediction Under Model Uncertainty by Hubert Wong B.A.Sc, U B C 1992 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF Doctor of Philosophy in T H E F A C U L T Y OF G R A D U A T E STUDIES (Department of Statistics) we accept this thesis as conforming to the required standard The University of British Columbia August 2000 © Hubert Wong, 2000 In presenting t h i s thesis in p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It i s understood that copying or p u b l i c a t i o n of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Si Department of The University of B r i t i s h Columbia Vancouver, Canada Date Abstract Existing criteria for evaluating the adequacy of a predictive model are model-based (e.g. AIC, BIC, MSPE) or empirical (e.g. PRESS and other cross-validation type criteria). We introduce a new class of "mongrel" criteria for on-line prediction that evaluates candidate predictors based on both model information and past empirical performance. Simulation results showed that the mongrel procedure produced more accurate predictions than the standard Bayes procedure for small sample sizes. This improvement was observed over a wide range of data-generators for the problem of variable selection in normal linear models. ii Contents Abstract i i Contents iii List of Tables v List of Figures vii Acknowledgements xi 1 Introduction 1 2 Mongrel Risk 8 3 Application to Normal Linear Models 18 3.1 Notation 19 3.2 Formulae for Choice and Weighting Strategies 21 3.3 Some Intuition 23 3.4 A Simulation Example . . . 26 3.4.1 Model Averaging 28 3.4.2 Model Choice 30 3.4.3 Summary 30 in 4 Finite Samples: Adaptive Selection 51 4.1 Model Averaging 55 4.1.1 Minimax meta-risk 55 4.1.2 Minimum-weighted-average meta-risk 56 4.1.3 Bayes-near-minimum meta-risk 58 4.1.4 Bayes-factor-decisive meta-risk 59 4.2 Model Choice 60 4.3 Relationship between mWA and mM meta-risk 60 4.4 Discussion 62 5 Finite Samples: Global Selection 89 5.1 Closeness to the Bayes solution 90 5.1.1 Approximating (5.6) 92 5.1.2 Simulation Results 94 5.2 Equalizing meta-risk 95 5.3 Discussion 96 6 Asymptotics 118 6.1 Consistency of Model Weights 118 6.2 Conditional Predictive Distributions 128 7 Discussion 130 Bibliography 136 Appendix 138 iv List o f Tables 3.1 Key for comparing the mongrel procedures to the Bayes procedure. 31 3.2 Summary comparison of the naive mongrel averaging strategy with S" = R/, to the Bayes strategy (S° = Y ( n ) ) 31 3.3 Summary comparison of the naive mongrel choice strategy with S£ = Kh to the Bayes strategy (S° = Y ( n ) ) 32 4.1 Summary comparison of the mM mongrel averaging strategy with T£ = S£ to the Bayes strategy 64 4.2 Summary comparison of the m M mongrel averaging strategy with T£ = Y(n) to the Bayes strategy 64 4.3 Summary comparison of the mWA mongrel averaging strategy with (T£, T£) = (S£, S£) to the Bayes strategy 65 4.4 Summary comparison of the mWA mongrel averaging strategy with (T£, T£) = (S£, Y ( n ) ) to the Bayes strategy 65 4.5 Summary comparison of the B N M mongrel averaging strategy with (T£, T£) = (S° S°) to the Bayes strategy 66 4.6 Summary comparison of the B N M mongrel averaging strategy with (T£, T£) = (S£, Y ( n ) ) to the Bayes strategy 66 4.7 Summary comparison of the B F D mongrel averaging strategy with (T£, T£) = (S£, S%) to the Bayes strategy 67 4.8 Summary comparison of the B F D mongrel averaging strategy with (T°,T£) = (S£,Y ( n )) to the Bayes strategy 67 5.1 Summary comparison of the ROB mongrel averaging strategy to the Bayes strategy 97 5.2 Summary comparison of the R O B mongrel choice strategy to the Bayes strategy 97 vi List of Figures 3.1 Performance of naive averaging strategies: a 2 ] 0 = 0.2, 72 = 0. 33 3.2 Performance of naive averaging strategies: 02,0 = 0.2, 72 = 0.2. 34 3.3 Performance of naive averaging strategies: a 2 > 0 = 0.2, 72 = 0.4. 35 3.4 Performance of naive averaging strategies: a 2 ; 0 = 0.2, 72,= 0. 36 3.5 Performance of naive averaging strategies: a2,o = 0.5, 72 = 0.2. 37 3.6 Performance of naive averaging strategies: a 2 > 0 = 0-5, 72 = 0.4. 38 3.7 Performance of naive averaging strategies: a 2 , 0 = 0.8, 72 = 0. 39 3.8 Performance of naive averaging strategies: a 2 , 0 = 0.8, 72 = 0.2. 40 3.9 Performance of naive averaging strategies: a 2 > 0 = 0.8, 72 = 0.4. 41 3.10 Performance of naive choice strategies: o 2 ) 0 = 0.2, 72 = 0. . . 42 3.11 Performance of naive choice strategies: a 2, 0 = 0.2, 72 = 0.2. . 43 3.12 Performance of naive choice strategies: a 2 ) 0 = 0.2, 72 = 0.4. . 44 3.13 Performance of naive choice strategies: a 2 | 0 = 0.2, 72 = 0. . . 45 3.14 Performance of naive choice strategies: a2)0 = 0.5, 72 = 0.2. . 46 3.15 Performance of naive choice strategies: a 2 ) 0 = 0.5, 72 = 0.4. . 47 3.16 Performance of naive choice strategies: a 2 > 0 = 0.8, 72 = 0. . . 48 3.17 Performance of naive choice strategies: a 2, 0 = 0.8, 72 = 0.2. . 49 3.18 Performance of naive choice strategies: a 2, 0 = 0.8, 72 = 0.4. . 50 vii 4.1 Meta-risk profiles for the first 12 sequences from an averaging strategy with T£ = T£ = S£: a 2 ) 0 = 0.2, 72 = 0. The solid (dashed) curve assumes the full (reduced) model is true. The dotted curve is a weighted average 68 4.2 Meta-risk profiles for the first 12 sequences from an averaging strategy with T£ = T£ = S£: a 2 ; 0 = 0.2, 7 2 = 0. The solid (dashed) curve assumes the full (reduced) model is true. The dotted curve is a weighted average 69 4.3 Histograms of the number of predictuals to include in S" that was selected by the minimax meta-risk procedure with T£ = S°: a2,o = 0.2, 7 2 = 0 . 70 4.4 Performance of m M averaging strategies: a2jC, = 0.2, 72 = 0. . 71 4.5 Performance of m M averaging strategies: a 2 ) 0 = 0.2, 72 = 0.2. 72 4.6 Performance of mM averaging strategies: a 2 , 0 = 0.2, 72 = 0.4. 73 4.7 Performance of mM averaging strategies: a2>0 = 0.5, 72 = 0. . 74 4.8 Performance of mM averaging strategies: a2]0 = 0.5, 7 2 = 0.2. 75 4.9 Performance of m M averaging strategies: a 2 i 0 = 0.5, 72 = 0.4. 76 4.10 Performance of mM averaging strategies: a 2, 0 = 0.8, 72 = 0. . 77 4.11 Performance of m M averaging strategies: a 2, 0 = 0.8, 72 = 0.2. 78 4.12 Performance of mM averaging strategies: a 2 ) 0 = 0.8, 72 = 0.4. 79 4.13 Performance of mWA averaging strategies: a 2 , 0 = 0.2, 72 = 0. 80 4.14 Performance of mWA averaging strategies: a 2 ] 0 = 0.2, 7 2 = 0.2. 81 4.15 Performance of mWA averaging strategies: a 2 > 0 = 0.2, 72 = 0.4. 82 4.16 Performance of mWA averaging strategies: a 2, 0 = 0.5, 72 = 0. 83 4.17 Performance of mWA averaging strategies: a 2 ) 0 = 0.5, 72 = 0.2. 84 4.18 Performance of mWA averaging strategies: a 2 , 0 = 0.5, 7 2 = 0.4. 85 4.19 Performance of mWA averaging strategies: a2,0 = 0.8, 72 = 0. 86 viii 4.20 Performance of mWA averaging strategies: a 2, 0 = 0.8, 72 = 0.2. 87 4.21 Performance of mWA averaging strategies: a2>0 = 0.8, 72 = 0.4. 88 5.1 Robustness profiles as a function of the number of predictuals included in S£ for ROB strategy: a2>0 = 0.5, 72 = 0.2 98 5.2 Histograms of the number of predictuals to include in S£ as selected by ROB strategy: a2o = 0.5, 72 = 0 2 99 5.3 Performance of ROB averaging strategy: a 2, 0 = 0.2 72 = 0. 100 5.4 Performance of ROB averaging strategy: a2>0 = 0.2 72 = 0.2. 101 5.5 Performance of ROB averaging strategy: a 2 ; 0 = 0.2 72 = 0.4. 102 5.6 Performance of ROB averaging strategy: a 2 , 0 = 0.5 72 = 0. . 103 5.7 Performance of ROB averaging strategy: a 2 ] 0 = 0.5, 72 = 0.2. 104 5.8 Performance of ROB averaging strategy: a 2 ; 0 = 0.5, 72 = 0.4. 105 5.9 Performance of ROB averaging strategy: a 2 , 0 = 0.8, 72 = 0. . 106 5.10 Performance of ROB averaging strategy: a 2 ] 0 = 0.8, 72 = 0.2. 107 5.11 Performance of ROB averaging strategy: a 2, 0 = 0.8, 72 = 0.4. 108 5.12 Performance of ROB choice strategy: a 2, 0 = 0.2 72 = 0. 109 5.13 Performance of ROB choice strategy: a 2 ) 0 = 0.2 72 = 0.2 110 5.14 Performance of ROB choice strategy: a 2 ] 0 = 0.2 72 = 0.4 111 5.15 Performance of ROB choice strategy: a2<0 = 0.5 72 = 0. 112 5.16 Performance of ROB choice strategy: a 2, 0 = 0.5 72 = 0.2 113 5.17 Performance of ROB choice strategy: a 2 ) 0 = 0.5 72 = 0.4 114 5.18 Performance of ROB choice strategy: a2>0 = 0.8, 72 = 0. 115 5.19 Performance of ROB choice strategy: a 2 ] 0 = 0.8, 72 = 0.2. 116 5.20 Performance of ROB choice strategy: a 2 ) 0 = 0.8, 72 = 0.4. 117 ix Distribution of MSPE(aff) - MSPE(a2xfmM) for predicting Y1Q with simulation parameters a 2 , 0 = 0.5, 72 = 0.2. (Bottom panel displays the smaller frequency range on an expanded scale.) . Acknowledgements I would like to thank Professors Jim Zidek, John Petkau, and Paul Gustafson for their invaluable advice and support all of these past years, even when it was not clear what this thesis was about. Special thanks go to my supervisor, Professor Bertrand Clarke, for being even crazier than me. Partial funding for this project was provided by NSERC Operating Grant OGP-0138122. HUBERT W O N G The University of British Columbia August 2000 xi Chapter 1 Introduction This dissertation describes a new approach to on-line prediction in the pres- ence of model uncertainty. Here, a prediction, or forecast - we use the two terms interchangeably - is a statement about the outcome of an as yet unob- served random variable. The term "on-line" indicates that we are predicting a sequence Y = (Yi, Y" 2 , . . . ) wherein the prediction for the outcome Yn+i, is made at time-point n using only the information that is available at that time. The term "model uncertainty" means that we have a collection of models (i.e. parametric families) M = {Mek : 0k 6 Qk,k € JC} as candidates for the distribution of Y . Each Mgk generates it's own candidate predictor Yk>n+i for Yn+\. However, because of model uncertainty, we require a criterion for evaluating the worth of each model in order to obtain the prediction that will be used, either by choosing one model from or averaging over the models in Ai. In existing literature, the criteria for evaluating the candidate models can be classified into one of two main types. We label them "model-based" and "empirical". The value of a model-based criterion depends the structure of an as- 1 siimed probability model. This dependence typically manifests as a measure of fit of the data to this model. For instance, a likelihood or an expected risk of the predictor computed conditionally on this model and the data would be model-based. In general, two different models will generate a different values for the criterion even if both models had generated the same sequence of pre- dictions in the past. Well-known examples of model-based criteria used for model choice include the Akaike Information Criterion (AIC) and variants on it such as the Bayes Information Criterion (BIC). These criteria are computed as minus two times the maximized log-likelihood plus a penalty term that de- pends on the number of fitted parameters in the model and sample size. For model averaging, the most active research area currently is in Bayesian model averaging wherein the adequacy of a model is judged on it's posterior probabil- ity. Recent references include Raftery et al (1997), Clyde (1999), Hoeting et al (1999). One criticism of these criteria is that they essentially measure the fit of the data to the model rather than evaluate the expected accuracy of the cur- rent prediction. To overcome this limitation, we can use a (decision-theoretic) risk criterion instead. Here, the risk Pi(Yk,n+i) of candidate predictor Yk,n+u assuming model i is true and conditional on the data, is evaluated for each i in turn. The value of the criterion is taken to be the average risk J2i &iPi(Yk,n+i) where a* is the weight assigned to model i. However, whichever of these model- based criteria is used, all of them can be criticized on the basis that they do not consider the past predictive performance of the candidate models. In contrast, an empirical criterion assesses the worth of a candidate model strictly on its observed predictive performance. That is, the worth of the predictor Ykj from model k is given by a loss function L(Ykj,Yj) that depends only on the observed values of Ykj and Yj. A l l other information is ignored. Thus if two models had generated the same set of predictions, 2 the value of the empirical criterion would be the same for both models re- garless of their underlying structure. The paradigmatic empirical approach in non-sequential settings is "leave-one-out" cross-validation. Here, individual observations Yj are omitted one at a time and Ykj is obtained from fitting the model using the remaining data. The criterion for assessing model k is the total Tfc = J2j L(Yk,j, Yj)- A smaller value of Tk indicates a better model. For squared error loss this criterion is the well-known predicted residual sum of squares, or PRESS, statistic (Allen, 1974). In a sequential setting, the predic- tion at a given time-point must be issued without knowledge of later data and hence the PRESS criterion is artificial. Dawid (1984) suggested that an appro- priate modification is to base the criterion on the sequence of one-step ahead prediction losses L(Yk>n+i, Yn+i) that already have been incurred where the forecast Yk,n+\ is based on only data known at time-point n. Subsequent work by Dawid and others (e.g., Dawid 1992, Sellier-Moiseiwitsch and Dawid 1993, Skouras 1996) developed the asymptotic theory underlying this "prequential approach" to forecasting. The empirical approach has several attractive features. Consider, for example, two sequences of on-line forecasts for the next-day maximum tem- perature where one of the sequences was generated by a meteorologist's sta- tistical model and the other sequence was generated by an old man based on "the feel in his bones". The performance of these two sequences would be non-comparable using a model-based approach since no model is available for the old man's sequence. But this comparison is easy to make using an em- pirical approach; simply define the criterion as Y,j L(Yk}j,Yj), say. While the model-free nature of an empirical criterion allows for comparison of arbitrary forecasting procedures, we do not view this property as an advantage necessar- ily. Indeed, we will indicate shortly why it is a disadvantage in small samples. 3 Instead, we feel that the most attractive feature of an empirical criterion is that it provides an evaluation based directly on past predictive performance and we seek to preserve this feature in our new approach. Initially, our work was motivated by the possibility that judicious use of the observed one-step ahead losses might allow for better predictions in small- samples while retaining the desirable asymptotic properties. The intuition is that losses incurred early in the sequence are likely to have less bearing on the quality of a candidate predictor than losses later in the sequence. Early on, the models are poorly estimated and therefore yield less precise information about performance. Dawid (1992) recognized this possibility and omitted the losses from early time-points in his simulations when computing the total loss. We, however, were interested specifically in such early losses since our focus was on small-sample performance. Through simulations, we investigated the influence of early losses in ad hoc fashion by downweighting the earlier losses when computing the total loss, that is, rather than using £ " = 1 L(Ykji, Yi) as the criterion when predicting for time-point n + 1, we used J2i=iWiL(Yk>i,Y) where the weights satisfied Wi < Wi+\. Some weighting choices we considered were: (1) tUj = i and (2) Wi — 0 if i < n/2, Wi = 1 if n/2 < i < n. Many of our choices yielded statistically significant improvements over using the simple total loss. However the magnitude of the improvement seemed small. Typically the reduction in the squared error prediction loss was around 0.5%. We felt that there were two major limitations to this initial approach. First, by taking a purely empirical view, we were left without a probabilistic framework. This meant we could not quantify the relative importance of resid- uals at different time-points and thus could not determine optimal weighting strategies. Second, the initial specification of candidate models and prior be- 4 liefs about their plausibility represented information that may not have been used fully by looking only at the incurred losses. This information would be particularly valuable early in the sequence when little data has accumulated. Hence, we concluded that a purely empirical approach was not suitable for evaluating candidate models in small samples. In Chapter 2, we describe a novel approach that combines aspects of both the model-based and the empirical approaches: we assume a probability framework for computing an expected risk but the expectation is computed conditional on a statistic S„ that reflects the observed predictive performance of the candidate models rather than conditional on the data values. To emphasize the dual aspects of the approach, we give the label "mongrel risk" to the resulting criteria. Examples of S„ that reflect empirical performance include past losses L(Yk,u Y), or past "predictuals" Yi — Yktt (the residuals that would arise from using the predictions from model k), from the candidate models. Different choices for S n generate different members in the class of mongrel risk criteria. The task is to determine good strategies for selecting S„ at any given time-point. A simple strategy is to set S n equal to the last, say, n/2 losses or predictuals always. But such a rule may be too naive as it ignores both information intrinsic to the model structure and supplied by the data. More sophisticated strategies, which we label as "global", would incorporate model and covariate information but does not use the outcomes of the response Y(„) in selecting S„. If a strategy also uses the outcomes of Y(„) to select S n , we call it "adaptive". In adaptive selection, we define a "meta-risk" for assessing the adequacy of each candidate S„. We describe two examples of such meta-risks: a minimax meta-risk and a weighted average meta-risk. In Chapter 3, we apply the mongrel risk criterion to prediction based on normal linear models. Explicit formulae for computing the mongrel risk 5 are given for the class S„ = U T Y ( n ) + c where U and c do not depend on Y( n )j that is, Sn is affine in the response vector. This class of S n includes past predictuals. We illustrate the implementation of the approach using a simulation study which uses naive choices of S„. We observed that setting S„ equal to a few recent predictuals resulted in more accurate predictions than the Bayes procedure (obtained by setting S n = Y( n )) in small samples for many but not all data-generating models. We investigate the adaptive selection of S„ for small-samples in Chapter 4. Specifically, we seek the optimal number of predictuals to include in S n . We argue that adaptive selection maximizes the improvement in predictive accu- racy. We implemented both the minimax and the weighted average versions of the meta-risk based on the simulated data used in Chapter 3. Our simulation results showed that in a model averaging context, the minimax meta-risk crite- rion produced more accurate predictions than the Bayes procedure uniformly over all of the scenarios tested. Moreover, the magnitude of the improvement was substantially greater than that seen in Chapter 3 where global choices of S n were used. In Chapter 5, we use robustness considerations to suggest global strate- gies. Unfortunately, these strategies are difficult to implement since the com- putations require integration over the distribution of Y ( n + 1 ) that are not tractable typically. We implemented one strategy for selecting S n in our sim- ulation study by approximating the needed quantities. Our results indicated that the performance of this strategy beat out the Bayes strategy in many scenarios but also lost badly in a few cases. The asymptotic theory for S„ is developed in Chapter 6. We character- ize the sub-class of affine S n that will yield asymptotic consistency of model weights. This consistency condition ensures that if one of the candidate mod- 6 els is true then as n —> oo, the predictor from the true model will always be chosen (in a model choice approach) or that the weight assigned to the true model tends to 1 (in a model averaging approach). In Chapter 7, we discuss additional simulation results that support the use of the mongrel procedure. In addition, we indicate areas of current and future work. 7 Chapter 2 Mongrel Risk Let Y = (Yi,Y2,...) be the sequence of random variables that is to be pre- dicted. At each time point n, we must issue a prediction concerning the value of Yn+l. To aid us in constructing the prediction for time point n + 1, typically the following information is available: 1. a p-vector of covariates X n + i whose elements may be related to Yn+i, 2. the outcomes and covariates already observed up to time point n, i.e., Y ( n ) = ( F i , . . . , Yn) and X ( „ ) , the n x p matrix with row i equal to X j , 3. prior information about the structure, which we call the model, that describes the probabilistic dependence of the outcomes on the covariates and the set of unknown parameters 6 that indexes the model, and 4. prior distributions on the values of the unknown 9. When the model posited in 3. is uncertain, we often entertain a collec- tion of candidate models M = {M6k : 6k € 0^, k G JC}. Each candidate model 8 k is assigned a prior probability akfi that reflects the plausibility that it is the true model. We say that model k is true if it contains the true distribution of Y and no sub-model (still in M) of model k contains this distribution. For each model k, we can proceed in a variety of ways to obtain a function of the observed data which we will use as the point prediction of Yn+i. If the loss function L(a, y n+i) describes the loss incurred by predicting using the value a = a(Y(„)) when Yn+\ = yn+\ obtains, the Bayes predictor Yk<n+\, conditional on model k being true, is the value of a minimizing the posterior risk, i.e., Yk,n+i = a rgminE M Y ( n ) I / (a , Yn+1) (2.1) Here, Efc|s indicates the conditional expectation given a statistic S assuming model k, marginalized with respect to the unknown parameters, is true. That is, if we let T be a minimal extension of S to Y ( „ + 1 ) , then for any function 0 (Y( n + i ) ) , E*| S <7(Y ( n + 1 ) ) = J 9{y{n+i))p(y(n+i)\0k,S)p{ek\S)dekdt (2.2) where p(-) denotes the appropriate conditional density and dt represents inte- gration over the sigma-field generated by T. Analogously, the notations C f c |s and Vfc|s will indicate the covariance and variance operators respectively. Each candidate model generates a forecast but we must give a single forecast that will be used. This problem is usually solved either by choosing one of the candidate forecasts or by using the forecast generated from a mixture over the candidate models. We use the term "model choice" to describe the first approach and the term "model averaging" for the second. To implement a model choice strategy, we require a criterion to assess the risk of the forecast derived from each candidate model. If model k were 9 true, the posterior risk of using Yk,n+i is Pk{Yk,n+V, Y ( N ) ) = Ek\Y(n)L{Yktn+i, Yn+i). (2.3) But because the true model is uncertain, we should also consider the posterior risk of using Yk>n+i when a different candidate model i / k is true. That is, we also consider Pi{Yk,n+i; Y ( „ ) ) = ~Ei\Y(n)L(Ykin+i,Yn+i). (2.4) In the Bayes decision approach, the overall assessment of risk for each candidate forecast is given by the weighted average over the collection of pos- terior risks that would be incurred by this forecast under different true models. That is, the adequacy of the predictor YktTl+i would be the average posterior risk p{Yk,n+l, Y ( n ) ) = Y;ai(Y(n))pi(Yk ,n+li ~Y(n) ) (2-5) ieK. where the model weights are a I ( Y ( N ) ) = P ( M I | Y ( N ) ) 1 (2.6) i.e., CVJ(Y(„)) is the posterior probability that model i is true. Based on the criterion p ( Y f c i n + 1 ; Y ( N ) ) , the best choice for the predictor of Yn+i is Yk*,n+i where k* satisfies k* = axgminp(Y f c n + 1 ; Y ( N ) ) . (2.7) For the special case of squared error loss, k* is the index of the model with the highest posterior probability, i.e., (2.7) reduces to k* = axgmaxdjfe(Y(N)). (2.8) 10 In the case of a model averaging strategy, the Bayes decision approach is to construct the predictive mixture distribution for Yn+i given by |̂v(B) = E a i ( Y w ) ^ | Y w (2-9) ieic where i * i | Y ( n ) is the posterior distribution of Yn+i given Y ( n ) assuming model i is true and the weight o i j ( Y ( n ) ) is determined as in (2.6), the same as for choice strategies. (Here, the subscript rn refers to the mixture rather than a candidate model.) The risk of using a = a(Y( n ) ) to forecast Yn+\ under the mixture is p ( a ; Y ( n ) ) = E m | Y ( n ) L ( a , y n + 1 ) (2.10) = £ a i ( Y ( n ) ) A ( a , Y ( n ) ) (2.11) i€K and the optimal forecast is now Ym,n+i = a igminp(a ; Y ( n ) ) . (2.12) For the special case of squared error loss, (2.12) reduces to Ym,n+l = Y ak(Y(n))Yk,n+\- (2-13) keic We wi l l refer to the solution defined by (2.1) through (2.13) as the Bayes procedure (for model choice or model averaging, as appropriate). The problem with the Bayes procedure is that none of the assessments of risk it uses reflect the true risk. First , the evaluation of Pi(-) assumes a distribution for Y ( n + i ) obtained by averaging over the distribution of the parameters whereas the true parameter values are fixed numbers which almost certainly do not equal those implied in the averaging. The Bayesian's position is that this averaging represents the best that one can do according to the rules 11 of probability. However this does not make the assessed risk true. Moreover, since only one of the models can in fact be true, at least some of the risks Pi(-) are based on an incorrect model. These observations raise two issues: (i) Are these risk assessments valid? and (ii) Can the criteria for assessing risk be modified to produce better predictors? As to (i), it is clear that we cannot avoid making some "incorrect" assessments of risk since the true distribution is unknown. So pending a better resolution we must be satisfied with the process of weighting the risks or the models according to beliefs about the merit of each model as a reasonable means to handle the model uncertainty. This procedure is no different than what we use to handle parameter uncertainty in the strictly parametric case. However, we answer (ii) by showing that better predictors can be obtained by changing the way we calculate « j ( - ) and pi(-). To put our proposed approach in context, let S n = S n (Y( n ) ) be any statistic and consider generalizing the Bayes decision procedure by replacing occurences of Y(„) by S„ in (2.4) through (2.13). The choice for S n need not be the same for each instance. We let S" and S£ denote the choices used when evaluating ai(-) and Pi{-), respectively. (More generally, we are free to choose a different S n for every instance of Y ( n ) . For example, the choice of S£ in (2.6) could be different for «i(S") than for av(S").) Our thesis is that better predictors can be obtained by choosing and/or S£ to be different from Y ( n ) . In this more general framework, the formula for the weights becomes Oi(SZ) = P(M, |S : n a (2.14) For a choice strategy, equations (2.4), (2.5), and (2.7) become (2.15) (2.16) 12 r = a rgnmj5 (n , B + 1 ;S« ,S£) . (2.17) If the loss is squared error and S£ = Y ( n ) , then it can easily be shown that (2.17) reduces to A:* = argmaaa i(SS). (2.18) For averaging strategies, (2.9) through (2.12) become Fm\S^n = I > i ( S M | S £ > (2.19) p(a;S£,S£) = E m | S s , s , L ( a , Yn+1), (2.20) y m, n + 1 = argminp(a ;S£ ,S£) . (2.21) For squared error loss (2.13) becomes Ym,n+l = £ a*(S£)Efc|Sp Y„+i. (2.22) keK When S£ = Y(„), the expectations in (2.22) are posterior means and so (2.22) further reduces to W i = E t t * ( S n ) V ' (2-23) The collection (2.14) through (2.23) defines a new class of risk cri- teria, indexed by the choice of a pair of cr-fields (S£,S£), for obtaining a predictor Yk*tTl+i (model choice) or Fm,n+i (model averaging). We will use Yn+X = F„+i(S£, S£) to refer to a predictor of either form {Yk*>n+1 or ^m,n+l)- We conjecture that suitable choices for S£ and S£ are vectors of statistics that reflect how well each of the candidate models have performed in predicting earlier data points. Such choices of S° and S£ motivate the following defini- tions. Definition 2.1 The mongrel risk of the predictor Yn+\ when model i is true is Pi{Yn+l-S£) = Eils?iL(Yn+1, Yn+l). (2.24) 13 Defini t ion 2.2 The average mongrel risk of the predictor Yn+X is j 5 (y n + 1 ;SS ,S£) = ^ ( S S t o f o + i l S S ) . (2.25) ieic The label "mongrel" is intended to reflect the hybridization of a model- based framework wi th empirical performance. The Bayes criterion is a special case in which we set S" = S£ = Y(n). The losses L(Ykj; Yj) from previous time-points j that would have been incurred had the predictor Ykj been used are natural candidates for statistics to be included in S° or S£. The following statistics are of particular interest and so we name them specifically. Defini t ion 2.3 The predictual resulting from using the predictor Ykj to pre- dict Yj is Rk,j =Yj — Yjfcj. (2.26) (The index k in the predictuals or the losses need not be the same as the index for the model that is under evaluation). B y choosing S£ and/or S£ of this form and conditioning on them, we obtain predictors that are functions of the actual performance of the candidate models rather than simply on data values. It is not obvious which predictuals or losses should be included and one of our main goals is to find good choices for (S", S£). We wi l l focus on using predictuals rather than losses for two reasons. The pragmatic reason is that in the normal linear models with which we wi l l be working, predictuals (more generally, any affine functions of Y(n)) allow us to evaluate the risks in (2.14) through (2.23) analytically. The conceptual reason is that losses, typically, do not distinguish between the bias and variance aspects in the error and this information may be relevant to assessing the quality of the candidates. 14 Since it is intuitively reasonable to expect that more recent predict- uals contain more information (if only because recently fitted models are more stable), we prefer the inclusion of more recent predictuals in (S°, S£). There- fore, we wi l l consider only collections of predictuals wherein the inclusion of a predictual from time-point j implies that the predictuals from all time-points greater than j are also included. Note that the special case of using the pre- dictuals from all past time-points is equivalent to using the Bayes procedure since the a-field generated by al l predictuals and the a-field generated by the data are equivalent. (The value of Y i can be recovered from the predictual at the first time-point. A t any time-point n > 1, the value of Yn can be recovered given the n-th predictual and Y ( n _ X ) . B y induction, Y(„) can be recovered from the set of all past predictuals.) The formulae (2.14) through (2.23) serve as a means for obtaining the predictor for a given choice of (S",S£) but say nothing about what (S°,S£) to use. (This step is unnecessary in the Bayes approach since in that case S° = S£ = Y ( n ) always.) One simple specification for (S°,S£) is to use always, say, only the n/2 (or some other pre-specified function of n) most recent predictuals when predicting for time-point n. This "naive" specification does not consider the structure of the underlying models or the observed data. In Chapter 3, we w i l l see that this choice often yields better predictions than those obtained using the Bayes procedure but does not do so consistently. A more sophisticated selection procedure could take into account the model structure and the covariate values X ( n ) but not the outcomes of Y ( n ) . This "global" approach, discussed in Chapter 5, is difficult to implement because the computations require an integration over the distribution of Y( n +i ) that is not tractable typically. However, we place less importance on the global approach because we believe that to obtain the greatest improvement, (S°, S£) must be 15 chosen adaptively, that is, rather than using the same ( S ° , S £ ) regardless of the outcomes of Y ( n ) , a different ( S ° , S £ ) is chosen for each sequence. The reasons for this belief will be developed more fully in Chapter 4. For now, we describe only the basic mathematical framework. In order to compare different (S", S£)'s, we need to assess the adequacy of each ( S « , S £ ) . D e f i n i t i o n 2.4 The meta-risk o / ( S " , S £ ) is a number that assesses the ade- quacyof(S%,Stt. The optimal ( S ° , S £ ) minimizes this meta-risk. There are two ways of interpreting meta-risk. The first way involves noting that any given ( S ° , S £ ) determines completely Yn+i. Thus, it seems reasonable that the meta-risk should be computed conditional on that choice of Yn+i. In this case, potential definitions for the meta-risk include the maximum risk of Yn+\ over different true models p v ( S £ , S £ ; T £ ) = m a x P j ( F n + 1 ; T £ ) , (2.27) or the weighted average risk of Yn+i over different true models p + ( S « , S £ ; T £ , T £ ) = ^al(T-)pi(Yn+l;T^ (2.28) i where we have introduced the statistics T° and T £ to emphasize that these statistics, which are used to evaluate ( S ° , S £ ) , are distinct from the statistics ( S £ , S £ ) , which are used to evaluate the candidate predictors. Taking T £ = S" and T £ = S £ would be natural because it would mean that we are using the same model weights and the same mongrel risks in both deriving and assessing Yn+i. However, we are not constrained to do so. Note that the dependence on (S", S £ ) is implicit to the construction of F n + 1 . 16 The second way of interpreting the meta-risk takes the view that the choice of (S£, S£) precedes consideration of how the predictor is subsequently obtained for the selected ( S " , S£). Thus, the meta-risk should reflect the risk of different candidate predictors that might arise (as opposed to does arise) from using ( S ° , S£) . This type of assessment is straightforward in the model choice context. For example, we could define the meta-risk as, say, the maximum average mongrel risk over candidate predictors p v ( S £ , S £ ) = m a x p ( l \ n + 1 ;S£ , S a (2.29) k or the weighted average of the average mongrel risk over candidate predictors ^ v ( S : , S : ; T : ) ^ a 1 ( T : ) f c 1 ; S : , S : ) . (2.30) k Again, it would be natural, but not necessary, to take T " = S " . Unfortunately, the analog to (2.29) or (2.30) in model averaging is not obvious. The class of "candidate predictors" in model averaging is the space of all functions of ( S ° , S£) . This class is problematic because of its large size and there does not appear to be a natural restriction. We w i l l use the first interpretation of meta-risk, in part because it is simpler to compute, but also because we find the second one odd in giving weight to the risk of a candidate predictor when it is known that that predictor wi l l not be used. 17 Chapter 3 Application to Normal Linear Models In this chapter, we derive the computational formulae needed to implement the criteria defined in (2.14) through (2.20) in the context of normal linear models. The candidate models are subset regression models. For expository simplicity, we will use the symbol S n in statements that are applicable to both S° and S£. Except where noted, in this chapter and for the remainder of thesis, we assume that S„ is a J = J{n) vector obtained as an affine transformation of Y(„), i.e., S n = U T ( Y ( n ) + c ) , (3.1) where the n x J matrix U and the n-vector c do not depend on Y(n). With- out loss of generality, we can assume that U is of full rank; if it is not, we simply remove linearly dependent rows until it is full rank. Ultimately, we are interested in S n only for the a-field it generates. Choices for Sn that satisfy (3.1) include the special cases where S n is a constant, S n = Y(n), and S„ is a vector of past predictuals. The inclusion of past predictuals in this group 18 follows from the fact that Bayes predictors in linear models are affine in Y ( „ ) . The manner in which predictuals impact mongrel risk assessments can be rather subtle. So we will provide some basic results in order to develop an intuitive understanding. We conclude the chapter by presenting simulation results for some simple forms of S n. Our primary intent here is to illustrate the procedure; the identification of "good" choices for S n is taken up in the subsequent chapters. 3.1 Notation We consider a collection of subset regression models of the form Y ( n ) | X ( n ) , & ~ A / - (x ( n ) D f c &,a 2 l ) (3.2) as candidate models where Dk is a p x pk matrix that "picks out" the pk covariates associated with model k. For simplicity, assume that a 2 is known. The parameter vector, Bk, is assumed to have a prior distribution irk given by pk~M(bk,rk). (3.3) For notational convenience, let Zfc,(n) = X ( n ) D j t . (3.4) Then the marginal density for Y ( n ) after mixing over the prior is J\f (ukjn, ^ktn) where Vktn = Zfci(n)bfc (3.5) $k,n = o2l + Zfc)(„)rfcZ^(n). (3.6) For each model k and given data Y ( n ) at time n, the Bayes rule with respect to 7Tjt and under squared error loss is 19 • for estimating Bk: J3k = arg ininEjfc iY^) (Bk - a ) 2 = Efc|Y (n )/?fc = C * , n * f c , n Y ( n ) + (bk - C k , n ^ n Z k t { n ) b k ) (3.7) where Cki7l — TkZT^ is the covariance between Bk and Y ( n ) . • for predicting Yn+i. Yfc,n+1 = Zl>n+1Pk = Zr ,n+lCfc,n*fc,nY( n) + Z £ n + 1 (bfc - C f c ) n*^Z f c ] ( r l)bfc) = U fc ,n+l Y (n) + (Zfc,n+1 - U f c ^ Z ^ ) ) bfc (3.8) where Ufc,n+1 = Zfc^+iCfc^fcjj (3.9) can be recognized as (marginalized with respect to fik) the covariance of Yn+i and Y ( n ) multiplied by the inverse of the variance of Y ( n ) . The Bayes predictors YktTl+1 arising from considering different models will constitute the collection of candidate forecasts. The predictual arising from using Yk<n+i, the predictor from model k, to predict Yn+i is Rk,n+1 = Yn+i — Y/c,n+l = Yn+i - u J ] 7 i + 1 Y ( n ) - (Zkjn+i - u^ n + 1 Z f c i ( n )) bk = u*kTn+i (Y(n+l) ~ Zfc>(„+i)bfc) (3.10) where u*kTn+l = ( -u^ n + 1 , l ) . 20 Given data up to time-point n, it can be useful to express the predictuals for all time points less than n in terms of the full data set (Y(„), Z f c j( n)), i.e., we write the predictual from time-point j as Rk,j = u£ j (Y( n ) - Z fc](n)b fcj (3-11) where u£j = (u£j, 0 , . . . , 0). The difference between (3.10) and (3.11) is purely notational. However, when S n is composed of predictuals only, the matrix U in (3.1) is easy to contruct if the predictuals are in the form (3.11); simply set each row in U to be with the desired choice of k and j. When model i is true, Rk,n+\ is normally distributed with mean and variance given by, respectively, Ej-ftfej = u°kTj (Zi^n)bi — Zk:(n)bkJ (3.12) ViRkJ = ufj%,nu°kij. (3.13) 3.2 Formulae for Choice and Weighting Strate- gies To implement the criteria (2.14) through (2.20), we require explicit formulae for c*j(S"), Pi(Yktn+i, S£) and F^. Since S n is affine in Y(„), these formu- lae can be obtained in a straightforward manner using the properties of the multivariate normal distribution. If model i is assumed to be true, then S n is distributed as a Af(pi, £j) with mean and variance Hi = U T Z t , ( n ) b t + U T c , (3.14) = a 2 U T U + U r Z , i W r , Z j ( n ) U . (3.15) 21 Then the mongrel risk of the predictor Yk,n+i is Pi(Yfc > n+i;S£) = EllSPnRl>n+1 + ( E , i i f c ) r i + 1 + Hfc.iEr1 (Spn - /ii)) 2 (3.16) where the covariance between the predictual Rk,n+i and S£ is H f c ) i = Ci(Rk,n+1, S£) = Dl^i>nU (3.17) and Dk,i = ^"^Zi^^riZ^+x - ^^Zfc^^TfcZfc^+i. (3.18) The predictive distribution under model i conditional on S£ is given by, ^ | S S ~ ^ , T 0 (3.19) where A = EiisP^+x = z£B+1b, + Z ^ . r . Z ^ U E - 1 (S£ - M l ) , (3.20) T i = VjispYn+x = ° 2 + z ^ n + 1 r j Z i ] n + x - z^„ + 1 r i z^ ( n ) uE~ 1 u T Zj ) ( n ) r jZj i „ + 1 = * 2 + Z j n + 1 ( r r 1 + a - X w U ^ U ) " 1 ^ ) ) - 1 Z i > B + 1.(3.21) The model weights a i (S£) are obtained using Bayes theorem which yields <*(SZ) = ^'°mfc(S5aV (3-22) Efca*,om*( sS) where a i ) 0 denotes the prior weight given to model i and m,(S*) = (2 7 r) J / 2 |E J | - 1 / 2 exp [-\(San - ^f^1^ - M l)} (3-23) is the marginal density of S" when model i is true. 22 3.3 Some Intuition It is worthwhile at this point to develop some intuition as to how past pre- dictuals represent information. The key point is that predictuals are relevant only when model uncertainty is present. Indeed, in Corollary 3.1 below, we see that when a single model is taken to be true a priori, the assessed quality of the predictor from this model does not depend on any affine function of Y ( N ) . One important notion relevant to choosing S£ and S£ is risk-sufficiency. Def ini t ion 3.1 (a) A a-field S£ is risk-sufficient for a i / a i ( S £ ) = OJJ(Y(„)) for every i. (b) A a-field S£ is risk-sufficient for p if pi(Yk,n+1; S£) = pi(Yk,n+l; Y ( N ) ) for every pair (k,i). In taking a mongrel risk approach, we need to avoid choosing a pair (S" ,S£) that is risk-sufficient since then our procedure devolves to the Bayes procedure that we are trying to beat. The following Lemma can be useful for verifying risk-sufficiency of a given S£. L e m m a 3.1 pi(Ykin+i;S^) = Pi(YktJl+x\Y(N)) •<=>• at least one of the follow- ing holds: (i) Ekii = 0, (ii) rank(U) = n, or (iii) Dk>i lies in the column space of U , i.e., there exists a vector d such Dk,i = U d . Proof The sufficiency of condition (i) is obvious from inspection of (3.16). If condition (i) does not hold, we need Ap = pi(YktTl+1; S£) - Pi(Yktn+1; Y ( N ) ) = 0. 23 Let Xk,i = C j ( i ? f c ) T l + 1 , Y ( n ) ) . Applying (3.16) and simplifying, we obtain Ap = - E ^ E " 1 - ^ + (EiRk>n+1 + S f c i j E t _ 1 ( S ^ - //0) = £ £ ^ ( 1 - P)£> f c l i + c ( Y ( n ) - !/0 r(I - P)£>fc,i (3.24) where P = U ( U T * i U ) - 1 U T * i (3.25) is a projection matrix onto the column space of U and c = E ; / 4 ! n + i + (Y( n) — ui)T(I + P)Dkti. Clearly A , = 0 iff (I - P)Dkii = 0. (3.26) Since rank(U) = n implies I — P = 0, condition (ii) is sufficient. Otherwise, we need the null space of I — P to contain Dk>i. But the null space of I — P is equal to the column space of U and hence we need Dkj to lie in the column space of U . This is condition (iii). | As an example consider that, conditional on a given model i, the min- imal sufficient statistic for the parameter is Z ^ n ) V f ~ * Y ( n ) . This suggests that a natural candidate for risk-sufficiency (for p) in a collection of K models is S£ = U T Y ( n ) where U = [ ^ n Z U n ) | • • • | ^ ^ n Z ^ ( n ) ) . Indeed, by taking components of d to be 0 or of the form r^Zi^+i , it is clear from (3.18) that this U can generate any Dktki for any choice of k and k'. This shows that S £ is risk-sufficient. Note that Lemma 3.1 also suggests that risk-sufficiency is a weaker notion of sufficiency than parametric sufficiency. If we have, say, only two candidate models, then we can obtain one-dimensional risk-sufficient statistic simply by setting U = Dkj. Moreover, this choice can be made ir- respective of the number of parameters in the models. This result contrasts 24 with parametric sufficiency where the minimal sufficient statistic usually has dimension equal to the number of parameters. (This phenomena may also explain what's happening in Lemma 3.2 below where only a few predictuals is needed to get risk-sufficiency.) This observation leads us to conjecture that given K candidate models, it should be possible to construct a K — 1 dimen- sional risk-sufficient statistic. Corollary 3.1 The risk pk(Yktn+i, S£) is constant in S£. Proof Setting i = k in (3.17), (3.18) gives Ektk = 0 and the result follows from condition (i) in Lemma 3.1. | This result is not surprising really. Conditional on a fixed model, the Bayes predictor is the optimal predictor so the corresponding predictual must be uncorrelated with any affine function of Y ( n ) . (Otherwise one can con- struct a better predictor and thereby contradict the optimality of the Bayes predictor.) The characterization of risk-sufficiency for a is more complicated and we are unable to provide a simple test for general choices of S£ and S£. However, for our purposes, we are more concerned with the behaviour of S" and S£ that are comprised of predictuals. In our computational work, we have established that when model i is true, the mongrel risk of the predictor from a different model, k say, can depend nontrivially on the past predictuals generated by model i or model k, that is, Fact 3.1 Suppose model i is true and let S" and S£ each consist of a set of past predictuals generated by model i. Then in general, Pl(Yktn+l]S^^Pl(Yktn+l]0). (3.27) 25 Thus, our mongrel risk procedure indeed generates predictions that are different from the Bayes procedure. But our computational results also in- dicate that care must be exercised to ensure that we do not choose a set of risk-sufficient predictuals. Fact 3.2 Suppose k and k! are a pair of nested models where model k con- tains p additional predictors. Let S% and S£ each consist of the p most recent predictuals from each model. Then Pk(Xk',n+i] S£) = Pk(Yk',n+u Y ( n ) ) , (3.28) pk,(Yk,!n+1;Spn) = pk,(Yk,tn+1;Y{n)), (3.29) and c*i(S£) = ai(Yin)) for both i = k,k' (3.30) That is, this choice of S° and S£ is risk-sufficient. These results also suggest that predictive risk assessments using pre- dictuals have a Markov property; which of two models is better is determined by only the most recent predictual(s) from each model. It can also be shown that including predictuals from only the larger model has no impact on any of the risk assessments. The combination of these two facts suggests that we include only the predictuals from the smaller model. This is what we do in all of our simulations. 3.4 A Simulation Example Through simulation we assessed the forecasting performance for several choices of (S°, S£). Data sequences of length forty were randomly generated according 26 to the model Yn = 70 + 7 l * l , n + l2X2,n + Cn (3.31) where Xi,n, X2tTl, tn were all independent standard normal variables. We fixed 7o = 1 and 71 = 0.8, but varied the value of 72 to be 0, 0.2, or 0.4 in different simulations. We feel that these three choices for 72 cover a large enough range for testing our methodology. With unit variances for X2<n and en, the correlation between Y and X2, conditional on Xx, is 72/1 /1 4- j2. When 72 = 0.4, this correlation is ~ 0.37, a reasonably strong association. We assume that any stronger association typically would have been evident, or at least suspected, prior to the analysis. Consequently, there would have been no doubt that such a covariate should be included in all of the candidate models, i.e., such a covariate would not be subject to our predictor selection procedures. The collection of candidate models contained two models: • Model 1 (the "reduced model"): containing the intercept and Xx only, i.e., Yn = p* + p*Xhn + en (3.32) • Model 2 (the "full" model"): containing the intercept and both X\ and X2, i.e., y„ = A> + PiXi,n + P2X2,n + en. (3.33) For prior distributions on the parameters, we assumed that (/?Q, f3{) ~ Af((l, 0.8), I) and (Po,Pi,P2) ~ j V ( ( l , 0 .8 ,0 .2) , I). We considered three choices for the a2j0, the prior probability of Model 2: 0.2, 0.5, and 0.8. Hence a total of 9 sce- narios (3 choices for 72 x 3 choices for a2fi) were considered. The number of sequences used in each scenario was m — 5000. 27 In all of the simulations presented here, we set S£ = Y ( n ) . The choices for S° were • R i (the most recent predictual from Model 1) • R 5 (the most recent 5 predictuals from Model 1) • R/j (the most recent half of (i.e., n/2, rounded down) predictuals from Model 1) • R n = Y( n ) (all past predictuals = full data). 3.4.1 Model Averaging Figures 3.1 to 3.9 plot the results from taking a model averaging approach. The top panel in each page plots of the mean squared prediction error -i m M S P E = - £ ( y n + 1 - F n + 1 ) 2 (3.34) m i-l incurred by using each specified choice of S". The second panel from the top shows the average weight that was assigned to the full model. The standard which we are trying to beat is the Bayes procedure in which S" = Y( n ) (labeled 'a2ff' in the plots). The bottom two panels compare the difference in M S P E between the Bayes procedure and the choice S" = R/j ('a2hf') or the choice S" = R 5 ('a25f'). On these two plots, a curve lying above zero indicates that our mongrel approach, 'a2hf or la25f, is beating out the Bayes procedure. To facilitate comparison across different plots, we classified the perfor- mance of 'a2hf relative to 'a2ff' into the groups shown in Table 3.1. Table 3.2 summarizes the results of this classification in 3 time intervals: 10 to 20, 20 to 30, and 30 to 40. 28 For small n, the mongrel strategy = Rh performed no worse and in most scenarios better than the Bayes strategy. The magnitude of the improve- ment decreased as the length of the sequence increased. When n is large and 72 = 0.4, the mongrel strategy performed worse than the Bayes strategy. At this point, one could question whether the better performance by the mongrel strategies is coincidental. That is, could a mongrel strategy be doing better because it gives higher weight to the "right" model than does the Bayes strategy just by chance? Figure 3.4 suggests that this is not the case. Here, the true model has 72 = 0. So a strategy that gives greater weight to the reduced model on average ought to perform better than one that gives less weight to the reduced model. We see in the second panel that the mongrel strategies ('a2hf and 'a25f') both give less weight to the reduced model than does the Bayes strategy so, on average, the mongrel strategies are at a disadvantage with respect to giving high weight to the right model. Yet, from the third and fourth panels, we see that the mongrel strategies are beating the Bayes strategy. This suggests that the mongrel strategies are more intelligent. The results for model averaging are qualitatively very similar to the results for model choice. Once again the mongrel strategy is no worse and sometimes better than the Bayes strategy for small n. Also, the mongrel strategy is worse when n is large and 72 = 0.4. In general, the mongrel strategy beats the Bayes strategy in more scenarios and by a greater degree than was seen in the model choice approach. The same counter-intuitive phenomenom seen in the choice approach occurs here. When 72 = 0, the mongrel procedures on average give less weight to the reduced model (the true model) than the Bayes strategy and yet tend to perform better early in the sequence. 29 3.4.2 Model Choice Figures 3.10 to 3.18 plot the results from taking a model choice approach. The panels contain the same information as in the model choice results except that the second panel now shows the proportion of times that the full model was selected. Table 3.3 summarizes these plots. The results for model choice are qualitatively very similar to the results for model averaging. Once again the mongrel strategy is never worse and sometimes better than the Bayes strategy for small n. Also, the mongrel strategy is worse when n is large and 72 = 0.4. In general, the difference in performance between the mongrel strategies beats the Bayes strategy in the model choice approach was smaller than what was observed in model averaging. The intelligent behaviour of mongrel strategies that was seen in the model averaging approach manifests here as well. When 72 = 0, the mongrel procedures on average choose the reduced model (the true model) less often than the Bayes strategy and yet tend to perform better (early in the sequence). 3.4.3 Summary The simulation results presented here provide some evidence that taking a mongrel risk approach often beats out the Bayes procedure, particularly at early time points. The advantange gained here, using simple choices for the mongrel risk criteria, decreased as time progressed and, when 72 = 0.4 became a disadvantage. Moreover, the magnitude of the gains was relatively small. However, these choices were only intended to illustrate of the technique. In the next chapter, we show that larger and more lasting gains can be obtained by optimizing the choice of ( S " , S £ ) . 30 Table 3.1: Key for comparing the mongrel procedures to the Bayes procedure. MSPE(Bayes) - MSPE(mongrel) +++ > 0.02 ++ between 0.01 and 0.02 + between 0 and 0.01 0 no clear difference - between 0 and —0.01 — between -0.01 and -0.02 < -0.02 Table 3.2: Summary comparison of the naive mongrel averaging strategy with S" = Hh to the Bayes strategy (S° = Y(n)). n r 2 ^2,0 72 10 to 20 20 to 30 30 to 40 0 0 0 0 0.2 0.2 0.4 ++ 0 + 0 i 0 ++ + + 0.5 0.2 ++ + + 0.4 ++ 0 — 0 + + + 0.8 0.2 ++ ++ ++ 0.4 ++ + 0 31 Table 3.3: Summary comparison of the naive mongrel choice strategy with S£ = Rh to the Bayes strategy (S° = Y ( n ) ) . n r 2 0:2,0 72 10 to 20 20 to 30 30 to 40 0 0 0 0 0.2 0.2 ++ + + 0.4 0 - — i 0 ++ + + 0.5 0.2 ++ 0 0 0.4 0 — 0 0 0 0 0.8 0.2 0 0 + 0.4 0 0 - 32 # of simulations = 5000 true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) I 10 l 20 30 l - r ~ 40 I I 10 I - r - 30 I 20 40 I Mean SE I 30 l 10 20 l 40 I s= o Mean SE 10 20 I 30 40 Time Figure 3.1: Performance of naive averaging strategies: a2fi = 0.2, 72 = 0. 33 # of simulations = 5000 true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) 10 20 30 _ : 40 l Time Figure 3.2: Performance of naive averaging strategies: a 2 ) 0 = 0.2, 72 = 0.2. 34 # of simulations = 5000 true coefficient values: mean = (1,0.8,0.4) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) Figure 3.3: Performance of naive averaging strategies: a 2, 0 = 0.2, 72 = 0.4. 35 CL CO E g> J D oo o d o o CL CO CO o o Mean SE I 10 I 20 l 30 _ J 40 l OJ ns I 5 CL CO O O CO o o Mean SE I 10 I 20 I 30 I 40 Time Figure 3.4: Performance of naive averaging strategies: a2>0 = 0.2, j2 = 0. 36 # of simulations = 5000 true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) 0.5 Mean SE 10 _ l 20 i — r ~ 30 40 Mean SE I 10 20 30 40 Time Figure 3.5: Performance of naive averaging strategies: a 2 ) 0 = 0.5, 72 = 0.2. 37 38 39 Figure 3.8: Performance of naive averaging strategies: a 2 , 0 = 0.8, 72 = 0.2. 40 Figure 3.9: Performance of naive averaging strategies: a 2 ] 0 = 0.8, 72 = 0.4. 41 # of simulations = 5000 true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(1,l,1) prior for small model: mean = (1,0.8) sd = diag(1,1) 10 20 _ J 30 I 40 l Figure 3.10: Performance of naive choice strategies: a 2 ) 0 = 0.2, 72 = 0. 42 # of simulations = 5000 true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) Time Figure 3.11: Performance of naive choice strategies: a 2 ] 0 = 0.2, 72 = 0.2. 43 44 # of simulations = 5000 true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.5 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) "i i i 1 r 0 10 20 30 40 J 1 I I L_ c2ff c2hf c25f - - - c21f ^ "i 1 1 1 r 0 10 20 30 40 Mean SE - r - 10 20 30 40 Time igure 3.13: Performance of naive choice strategies: a2y0 — 0.2, 72 = 0. 45 0 10 20 30 0 10 20 30 Time Figure 3.14: Performance of naive choice strategies: a 2 , 0 = 0.5, 72 = 0.2. 46 Figure 3.15: Performance of naive choice strategies: a 2 j 0 = 0.5, 72 = 0.4. 47 Figure 3.16: Performance of naive choice strategies: a2,0 = 0.8, 72 = 0. 48 Time Figure 3.17: Performance of naive choice strategies: a 2, 0 = 0.8, 72 = 0.2. 49 # of simulations = 5000 true coefficient values: mean = (1,0.8,0.4) sd = (0,0,0) prior prob. on big model: prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) Time Figure 3.18: Performance of naive choice strategies: a 2 , 0 = 0.8, 72 = 0.4. 50 Chapter 4 Finite Samples: Adaptive Selection We classify the methods for choosing (S°, S£) as either global or adaptive. In a global method, the choice is made without regard to the past outcomes Y( n). In contrast, an adaptive method selects a different (S°, S£) for each sequence based on Y(n). The reason for classifying in this way is to differentiate between methods according to whether the method uses past predictive performance. Note that a choice of (S", S£) that depends on X(„) is considered to be global in this classification scheme. We conjecture that adaptive methods yield better results than global methods. The intuition is that using all of the data, i.e., setting S° = S£ = Y( n), provides the most accurate risk assessments for the 'typical' sequence but generates misleading assessments in other sequences. An adaptive method identifies and compensates for the misleading results by conditioning on past predictive performance. To implement an adaptive selection method, we must specify a meta- 51 risk criterion, such as given by (2.27) or (2.28), that is to be used to evaluate each choice of (S", S£). We will use the simulation results from the framework described in Chapter 3 to help illustrate the method. As in Chapter 3, we set S£ = Y( n ) and vary the choice for S". Here, we optimize the choice over the class S" G { R j : J = 0 , 1 , . . . , n}. The index J represents the number of most recent predictuals from the reduced model that are included in S°. Recall that if we take J = n, then the optimal choice corresponds to that from the Bayes procedure. The intuition underlying the meta-risks p v and p+ is seen most easily by considering meta-risk profiles. Defin i t ion 4.1 The meta-risk profile assuming model % is true, is the collec- tion of meta-risks pi(Yn+i;T^) generated by varying (S",S£). Thus, in our simulation framework, the meta-risk profile under a given model for predicting Yn+X contains the n+1 points corresponding to (S°, Sn) = (R j , Y( n ) ) , J = 0 , 1 , . . . , n. Figure 4.1 displays graphically the meta-risk pro- files as a function of J for each of the first 12 simulated sequences of data. The value being predicted here is Yw. The prediction Yn+X was generated using an averaging strategy and we have set T£ = S". In this particular scenario, 72 = 0.4 in the data-generator. The solid curve connects points in the meta- risk profile assuming the big model true the dashed curve connects points in the meta-risk profile assuming the small model is true. In specifying a meta-risk criterion, we must consider the values on both these curves since the true model is unknown. On average, we expect that meta-risks computed conditional on more predictuals (larger J) tend to be more accurate. But, for the reasons discussed in Chapter 2, none of these meta-risks is the true risk so it cannot be argued that conditioning on full 52 data ( J = n) is necessarily optimal. Instead, we should feel free to use fewer predictuals in S° if the resulting predictor is perceived to have more desirable characteristics. For example, if we find that the meta-risk of the predictor based on full data, say J = 10, is large under one of the models, call it k, but the meta-risk of the predictor based on J = 9 is (relatively) small under all of the models, then we might prefer to select J — 9 to avoid the large meta-risk that would be incurred if in fact model k was true. There are a variety of reasonable ways of reducing the meta-risk profiles to a choice of J . If for each J , we take the maximum of the meta-risks over all of the models, we get the maximum meta-risk criterion p v described by (2.28). Then, we obtain a "minimax" strategy by choosing the number of predictuals in S£ to be JmM = a rgmjnp v (Rj , Y ( n ) ; T £ ) (4.1) = a rgminmaxp f c (y n + i (R j , Y ( n ) ) , T £ ) . (4.2) (Note that in contrast to a standard risk function plot in which the parameter appears along the horizontal axis, the meta-risk profile plot places the elements in the decision space, i.e., the choices for S", on this axis.) Alternatively, if for each J , we average the values over all of the models using weights aj(T"), we get the weighted average meta-risk criterion p + described by (2.28). Then we obtain a "minimum-weighted-average" strategy by choosing the number of predictuals in S" to be JmWA = a r g m i n p + ( R J , Y { n ) ; T ^ , T ^ ) (4.3) = a r g m i n ^ a ^ T ^ p ^ y ^ ^ R ^ Y ^ ) , ^ ) . (4.4) J k Additional strategies are discussed later in this chapter. 53 In Figure 4.1, the points in the weighted average meta-risk profile, with T£ = S£, are connected by the dotted curve. Let to denote the panel in the i-tb. row and j-th column. We observe that for the five sequences in [1,2], [2,1], [2,2], [4,3], and [5,3] the maximum meta-risk at each J corresponds to the big model and that the minimum over these maxima appears to occur when J = 9, that is, when S" is equivalent to the full data. In contrast, for the sequence in, say, [1,1] the meta-risk profiles cross and the minimum of the maxima appears to occur when J = 6. The weighted average meta- risk profiles exhibit similar patterns; the minimum is achieved with J = 9 for the sequences in [2,1], [2,3], and [3,3] but with J < 9 in all the remaining sequences. We conjecture that when the minimum meta-risk is achieved with J < 9, the predictor based on the minimizing S" is less sensitive to model mis-specification than the Bayes predictor (for which J = 9). The general pattern seems to be that the meta-risk profile under the big model tends to decrease as J increases, whereas the profile under the small model remains fairly flat for all J or increases as J increases. The meta-risk profiles for the choice Tpn = Y(n) (see Figure 4.2 are somewhat different in that the profile under the big model does not tend to decrease but rather remains relatively flat or increases slightly with J. The profile under the small model behaves in much the same way as when T£ = . The difference in behaviour for the profile under the big model results in much different values of J M M - To assess how often the minimum was achieved for each value of J , we constructed histograms at time-points n = 10, 25, and 40. These histograms are shown in the first column of Figure 4.3 for a minimax strategy with T£ = S". The last bin in each histogram represents the proportion of times that full data was selected for S". For this particular scenario, we see that the 54 modal choice was full data, but that the choices are spread over a very wide range with a tendency to picking a large number of predictuals. The second column in Figure 4.3 shows the histograms for when T £ = Y ( „ ) . In this case, the modal choice is to use no predictuals and the tendency is to select a smaller number of predictuals. These histograms will be seen to be useful for explaining the characteristics of different strategies; we defer further discussion to the following sections. We will discuss the results from taking a model averaging approach in detail and only comment briefly on the model choice approach. The reason is that we are able to implement the needed computation for model averaging but not for model choice. 4.1 Model Averaging Through simulations, we first examined the performance of the "minimax" and "minimum-weighted-average" strategies from a model averaging approach. Be- cause the results suggested deficiencies in the mWA strategies, we subsequently also considered two modified strategies which we label "Bayes-near-minimum" and "Bayes-factor-decisive". We look at the performance of each of the four strategies in turn. 4.1.1 Minimax meta-risk Figures 4.4 through 4.12 show the performance obtained by using a mini- max (mM) meta-risk strategy (see (4.2)). The M S P E obtained using the two strategies ( T £ , T £ ) = ( S « , S « ) (labeled 'a2xxmM') and ( T « , T £ ) = ( S « , Y ( n ) ) ('a2xfmM') are shown in the top panel. For comparison, we have also plot- ted the M S P E obtained by the Bayes strategy ('a2ff') and, for reference, the 55 M S P E that would have been obtained had the full model ('big') or the reduced model ('small') been used at all times. The second panel from the top shows the average weight given to the big model. The third and fourth panels show the difference in M S P E between the Bayes strategy and each of the two mWA strategies. Positive values in these two plots indicate that the mWA strategy is better than the Bayes strategy. The most important feature of 'a2xxmM' is that it never performs worse and often performs substantially better than Bayes. The performance of 'a2xfmM' is generally similar to 'a2xxmM'. However, 'a2xfmM' performs slightly worse than Bayes when 72 = 0 and 0:2,0 is 0.2 or 0.5. This loss is mitigated somewhat in that the size of the gain, when present, tends to be substantially greater than the gain seen for 'a2xxmM'. Both m M strategies ex- hibit greatest improvement early in the sequence. Tables 4.1 and 4.2 summa- rize the comparison of the Bayes strategy to the 'a2xxmM' and the 'a2xfmM' strategies, respectively, over all of the scenarios. The similarity in performance is rather surprising since the histograms of J M M (e.g. Figure 4.3) are radically different. Whereas 'a2xxmM' tends to pick all of the data or nearly all of the data, 'a2xfmM' tends to use none or nearly none of the data. At present, we do not have a good explanation for this difference. 4.1.2 Minimum-weighted-average meta-risk Figures 4.13 through 4.21 show the performance obtained by using a minimum- weighted-average strategy (see (4.4). The M S P E obtained using the two strategies (T° T£) = (S»,S°) (labeled 'a2xxwa') and (T£,T£) = ( S « , Y ( n ) ) ('a2xfwa') are shown in the top panel. For comparison, we have also plot- 56 ted the M S P E obtained by the Bayes strategy ('a2ff') and, for reference, the M S P E that would have been obtained had the full model ('big') or the reduced model ('small') been used at all times. The second panel from the top shows the average weight given to the big model. The third and fourth panels show the difference in M S P E between the Bayes strategy and each of the two mWA strategies. Positive values in these two plots indicate that the mWA strategy is better than the Bayes strategy. Early in the sequence, the 'a2xxwa' strategy beat the Bayes strategy substantially in 6 of the 9 scenarios, was slightly better in 1 scenario, was slightly worse in the scenario with 72 = 0.4 and a2,0 — 0.5, and substantially worse in the scenario with 72 = 0.4 and a2,0 = 0.8 As time increased, the performance of 'a2xxwa' relative to Bayes decreased in every scenario. The trend was consistent in every scenario. By time-point 40, 'a2xxwa' beat Bayes substantially in only 1 scenario. In all 3 scenarios where 72 = 0.4, 'a2xxwa' was substantially worse. In the remaining 5 scenarios, 'a2xxwa' was slightly better in three of them, about the same in one, and slightly worse in one. The 'a2xfwa' strategy performed similarly to the 'a2xxwa' strategy when 72 = 0. But for most of the other cases, 'a2xfwa' performed worse than 'a2xxwa'. Generally, 'a2xfwa' behaved more like Bayes than 'a2xxwa' did, that is, while the differences in performance between 'a2xxwa' and Bayes were often large the differences between 'a2xfwa' and Bayes were relatively smaller. Table 4.3 (4.4) summarizes the comparisons of the Bayes strategy with the 'a2xxwa' ('a2xfwa') strategy for all of the figures. In general, the performance of the mWA strategies appears to deteri- orate as 72 increases. Additionally, the problem seems to be enhanced when a 2 , 0 is small. We conjecture that using only some of the predictuals to update the model weights inhibits the identification of when the big model is true and 57 this identification is needed for avoiding bias. This conjecture is supported by the plots of the mean posterior weight assigned to the big model (second panel). When 72 = 0.4, a 2 ) 0 increases relatively rapidly for the Bayes strategy as data accumulates but not for the mWA strategies. The mWA strategies behave oppositely to the m M strategies in that while mWA strategies work better when 72 is small the m M strategies work better when 72 is large. 4.1.3 Bayes-near-minimum meta-risk Because averaging generally is considered a good way of handling model un- certainty, we thought that it would be nice if the minimum-weighted-average meta-risk strategy could be "patched up" to perform well even when 72 is large. In a large proportion of the sequences, the average meta-risk profile often dropped to its minimum p m i n = p+(Yn+1(RJmWA, Y ( n ) ) ; T ° , T£) once a small number of predictuals was included in S" and then remained at relatively flat for all larger J . If we take the view that a small difference between p m ; n and the meta-risk ps = p+(Y„+i(R n, Y( n)); T", T£) for the full data (Bayes) procedure is not important, then we might want to default to using the full data rather than the number of predictuals corresponding to the minimum risk. That is, the number of predictuals to include in S ° is for a suitable cut-off value c. We call this the "Bayes-near-minimum" (BNM) strategy. Tables 4.5 and 4.6 summarize the results for B N M strategies with T£ = S" and T£ = Y( n ) , respectively. (We have not included the Figures that (4.5) 58 would be analogous to Figures 4.13 through 4.21 due to length considerations.) The cut-off value was c = 1.05. While this approach managed to reduce the magnitude of the poor performance when 72 = 0.4, it did not eliminate it completely. Moreover, it also tended to reduce the size of the gains. Hence, we did not pursue this modification believing it to be ineffective. This approach yielded another interesting fact. The given cut-off re- sulted in over 95% of the predictions, typical across all scenarios, being based on using all of the data. Yet, the differences in performance between the Bayes strategy and the N M strategies were sometimes nearly as large as those between the Bayes strategy and the mWA strategies. This result provides fur- ther evidence that a small proportion of the sequences generates much of the differences in performance (for better or for worse). 4.1.4 Bayes-factor-decisive meta-risk Another way of assessing whether to default to the full data is to use the Bayes factor (BF) of the small model with respect to the big model, say. The idea here is similar to "Occam's window" - a model that has little support from the data (indicated by a small Bayes factor) ought to be viewed as discredited and discarded from consideration. Since we think that, on average, mongrel criteria give a less precise measure than the Bayes criterion for the true risk, we ought to avoid use of the mongrel criteria when we have confidence in the Bayes criterion (which we assume is reflected in a sufficiently large or small BF) . Hence, the number of predictuals to include in S" is n if B F > c or B F < 1/c o.w. (4.6) 59 for some cut-off value c. We call this the "Bayes-factor-decisive" (BFD) strat- egy- Tables 4.7 and 4.8 summarize the results for B F D strategies with T£ = S° and T n = Y ( n ) , respectively. (Again, we have not included the Figures that would be analogous to Figures 4.13 through 4.21 due to length consider- ations.) The cut-off Bayes factor was c = 4. Overall, the B F D strategies were not successful in improving on the mWA strategies. The modification more often reduced the gains than reduced the losses seen originally in the mWA strategies. 4.2 Model Choice Unfortunately, we are unable to optimize meta-risk when using model choice approach because the meta-risks cannot be evaluated in closed form. Specifi- cally, Yn+i takes the value Y i , n + i if Y( n ) is in the set Si — {p (Y 1 ) I l + 1 ; S", S n ) < p(y 2 ,n+i; S£,S£)} or the value Y2>n+i if Y( n ) is in the complement, <S2, of Si. Hence the meta-risk is ft(yB+1; T n ) = E i | T { . £ ( y n + 1 - YwfXSk (4.7) k where XA is the indicator of the set A. The integral in (4.7) is not tractable analytically. 4.3 Relationship between mWA and mM meta- risk A standard result in estimation theory is that, typically, the minimax deci- sion is equivalent to the Bayes decision under a least favourable prior on the 60 parameter. A natural question to ask is whether the m M meta-risk choice of ( S " , S £ ) is equivalent to the mWA meta-risk choice. The following result shows the m M and mWA procedures are equivalent only in a trivial case. L e m m a 4.1 Let TT* = argmax m i n i V ^ f c + i ; T n ) (4-8) n (Sn>Sn) k where TT = (TTI, ..., TTK) is a probability on the model space. That is, TT* is a least favourable prior. Then the mM and m WA meta-risks yield the same choice for (S", S n ) (almost surely in Y^) iffVk, ak(Tan) = TT*k. (4.9) Proof The problem can be viewed as a finite zero-sum game in which the choice of ( S £ , S £ ) represents the decision to be made and the model space is the parameter space. Then the Minimax Theorem (Berger, p. 345; 1985) applies so that the game has value V = inf suP5>*p f c(rn + 1;Tn) = sup inf Y,*kPk(Yn+l; T£) (4.10) s n fc w s fc where the infimum is over all randomized choices of ( S ° , S n ) and the supremum is over all probabilities on the model space. By Lemma 1 (Berger, p.318; 1985), the value of the game can also be expressed as V = , s L n L 8 U P * & + i ; T S ) = inf . E ^ A + i j T S ) . (4.11) Finally, the infimums and supremums can be replaced by minima and maxima since both the model and decision spaces are finite and therefore closed, so we have , c m i c np N " \ A X p k ; T ^ = m a x ^ S nkPk (Yn+i; T£). (4.12) (S",S n ) * (Sg,S n ) 61 The choice of (S",S£) on the L.H.S. of (4.12) corresponds to the choice ob- tained using the m M strategy. The choice of (S£, S£) on the R.H.S. of (4.12) will correspond to the choice obtained using the mWA strategy iff Vfc, ak(T%) = <• I However, since 7r£ is constant in T", it is clear that (4.9) is possible only if T" is held constant for all v. We are primarily interested in the case where T" varies with (S",S£) and so in general the results we obtain using the minimax meta-risk are not equivalent to the results obtained using the minimum weighted average meta-risk for any such choice of T£. 4.4 Discussion At first glance, some of the differences in behaviour of these strategies seem natural while others seem counterintuitive. The results suggest that the mM strategies are more robust than the WA strategies. Whereas 'a2xxmM' never performs worse than Bayes and 'a2xfmM' at most performs slightly worse than Bayes, both 'a2xxwa' and 'a2xfwa' perform much worse than Bayes when 72 = 0.4. This result agrees with the usual sense in which minimax strategies are robust. On the other hand, weighted averaging strategies are expected to perform better overall. This does not seem to occur here - the gains from the mM strategies appear to exceed those seen for the WA strategies in general. A partial explanation may be that the optimization over the meta-risk is not intended to generate an optimal predictor directly but merely to decide on a good (S£,S£) on which the optimal predictor will then be computed. In this case, the most desirable property of (S",S£) is that it minimizes the chance that our subsequent choice of a predictor will be poor because it derives 62 from a bad model. That is, minimizing meta-risk is inherently a robustness issue and that may be why the mM meta-risks generate better results. A comparison of strategies using T £ = S £ to T £ = Y(„) shows a clear pattern. For all four methods of selecting the optimal value of J , the choice T £ = S " generated more accurate predictions in nearly all of the scenarios. This pattern would be expected if our reasoning for the advantages of using mongrel risk is correct, i.e., we get better assessments of risk by conditioning on empirical performance rather than the raw data. However, this is only a conjecture at this time. Note that if the choices for and S £ are varied independently of each other, then we would have n2 potential choices for ( S " , S £ ) since each of the two statistics can range from 0 up to n — 1 predictuals. We have not optimized over all n 2 choices for two reasons. First, the expectations in the meta-risks (2.27) or (2.28) can be evaluated analytically only when t r ( S ° ) C CT(T£). (SO that the weights ak(S") can be taken outside the expectation E J | T P ) . In gen- eral, T £ is not set to be equivalent to the full data so some choices of S " cannot be handled. The second reason is that we wish to dissociate the effects due to varying S " from the effects due to varying S £ . We have examined only the impact of varying S £ because this situation was easiest to implement compu- tationally. The case with S £ varied requires additional computational effort and has been left for future investigation. 63 Table 4.1: Summary comparison of the m M mongrel averaging strategy with T £ = S ° to the Bayes strategy. n r 2 72 10 to 20 20 to 30 30 to 40 0 + 0 0 0.2 0.2 + 0 0 0.4 + + + i 0 ++ + 0 0.5 0.2 ++ + + 0.4 +++ ++ + 0 ++ ++ + 0.8 0.2 +++ ++ ++ 0.4 +++ ++ + Table 4.2: Summary comparison of the mM mongrel averaging strategy with T £ = Y ( n ) to the Bayes strategy. n r 2 «2,o 72 10 to 20 20 to 30 30 to 40 0 - - - 0.2 0.2 ++ ++ ++ 0.4 +++ +++ +++ i 0 0 - - 0.5 0.2 ++ ++ ++ 0.4 +++ +++ +++ 0 ++ ++ + 0.8 0.2 +++ ++ ++ 0.4 +++ ++ + 64 Table 4.3: Summary comparison of the mWA mongrel averaging strategy with ( T ° T£) = (S«, S«) to the Bayes strategy. n r 2 ^2 ,0 72 10 to 20 20 to 30 30 to 40 0 ++ + + 0.2 0.2 + 0 - 0.4 — i 0 +++ ++ + 0.5 0.2 ++ + 0 0.4 - 0 +++ ++ ++ 0.8 0.2 +++ ++ + 0.4 ++ Table 4.4: Summary comparison of the mWA mongrel averaging strategy with (T£, T£) = (S*, Y ( n ) ) to the Bayes strategy. n r 2 ^2 ,0 72 10 to 20 20 to 30 30 to 40 0.2 0 0.2 0.4 ++ + + 0 + i 0.5 0 0.2 0.4 ++ 0 + + 0.8 0 0.2 0.4 0 65 Table 4.5: Summary comparison of the B N M mongrel averaging strategy with i ' Jn) to the Bayes strategy. n r 2 72 10 to 20 20 to 30 30 to 40 0 ++ + + 0.2 0.2 ++ + 0 0.4 — i 0 ++ ++ + 0.5 0.2 ++ 0 0 0.4 0 — 0 ++ + + 0.8 0.2 ++ + 0 0.4 ++ - - Table 4.6: Summary comparison of the B N M mongrel averaging strategy with (T° T£) = (S£, Y ( n ) ) to the Bayes strategy. n r 2 «2,0 72 10 to 20 20 to 30 30 to 40 0 ++ + + 0.2 0.2 ++ + 0 0.4 — i 0 + • + + 0.5 0.2 0 - 0 0.4 - — — 0 0 0 0 0.8 0.2 0 0 0 0.4 0 0 0 66 Table 4.7: Summary comparison of the B F D mongrel averaging strategy with ( T S , T B ) = (S«,SS) to the Bayes strategy. n r 2 0*2,0 72 10 to 20 20 to 30 30 to 40 0 + + + 0.2 0.2 0 - — 0.4 i 0 ++ ++ + 0.5 0.2 0 0 - 0.4 — 0 ++ ++ + 0.8 0.2 ++ + + 0.4 0 — — Table 4.8: Summary comparison of the B F D mongrel averaging strategy with (T° T£) = (S£, Y ( n ) ) to the Bayes strategy. n r 2 Ot-2,0 72 10 to 20 20 to 30 30 to 40 0.2 0 0.2 0.4 + 0 + + i 0.5 0 0.2 0.4 ++ 0 ++ 0 + 0.8 0 0.2 0.4 — - + 0 67 Meta-risk profiles of sample sequences - a2xx Figure 4.1: Meta-risk profiles for the first 12 sequences from an averaging strategy with T£ = T£ = S£: a 2 , 0 = 0.2, 7 2 = 0. The solid (dashed) curve assumes the full (reduced) model is true. The dotted curve is a weighted average. 68 Meta-risk profiles of sample sequences - a2xf 00 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 J (= # of predictuals in S) J (= # of predictuals in S) J (= # of predictuals in S) Figure 4.2: Meta-risk profiles for the first 12 sequences from an averaging strategy with T £ = T£ = S£: a2<0 = 0.2, 72 = 0. The solid (dashed) curve assumes the full (reduced) model is true. The dotted curve is a weighted average. 69 Frequency of selecting given # of predictuals - a2**mM a2xxmM a2xfmM 0 1 2 3 4 5 6 7 8 All # of predictuals in S _ 3 6 9 12 15 18 21 All # of predictuals in S 1A 0 4 8 12 17 22 27 32 37 # of predictuals in S m a cc 8" J5 CD CC C4 O O OJ O O O o OJ o o o d 0 1 2 3 4 5 6 7 8 All # of predictuals in S TTTTTTTTTTTTiTrT^ 0 3 6 9 12 15 18 21 All # of predictuals in S TlTTnTTTTTTTlTMTTT TTTTTnTrrlTfl 0 4 8 12 17 22 27 32 37 # of predictuals in S Figure 4.3: Histograms of the number of predictuals to include in S£ that was selected by the minimax meta-risk procedure with T£ = S": a 2 j 0 = 0.2, 72 = 0. 70 a2ff a2xxmM a2xfmM big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) I 10 I 20 30 I 40 a2ff a2xxmM a2xfmM big small Time Figure 4.4: Performance of mM averaging strategies: a 2 ] 0 = 0.2, 72 = 0. 71 a2ff a2xxmM a2xfmM big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) —r~ 10 _J 20 l I 30 l 40 I a2ff a2xxmM a2xfmM big small -r— 20 - r - 30 L _ 10 _l 40 l Mean SE - r - 20 i 10 30 i 40 L _ Mean SE 10 20 - 1 - 30 40 Time Figure 4.5: Performance of mM averaging strategies: a 2 , 0 = 0.2, 72 = 0.2. 72 a2ff a2xxmM a2xfmM big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0.4) sd = (0,0,0) prior prob. on big model • prior for big model: mean = (1,0.8,0.4) sd = diag(l,l,1) •prior for small model: mean = (1,0.8) .sd = diag(1-,1) • I 10 I 20 _ l 30 — T - 40 I a2ff a2xxmM a2xfmM - - - big "- small Mean SE I 10 I 20 30 40 Time Figure 4.6: Performance of mM averaging strategies: a2fi = 0.2, 72 = 0.4. 73 a2ff a2xxmM a2xfmM big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.5 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) 20 l 10 _ l 30 l 40 a2ff a2xxmM a2xfmM big small Time Figure 4.7: Performance of mM averaging strategies: a 2 i 0 = 0.5, 72 = 0. 74 a2ff a2xxmM a2xfmM big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model • prior for big model: mean = (1,0.8,0.4) sd = diag(1,l,1) prior for small model: mean = (1,0.8) sd = diag(1,1) 10 _i 20 _i 30 40 I a2ff a2xxmM a2xfmM t big small Time Figure 4.8: Performance of mM averaging strategies: a 2, 0 = 0.5, 72 = 0.2. 75 a2ff a2xxmM a2xfmM big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0.4) sd = (0,0,0) prior prob. on big model prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) •prior.for small model: mean = (-1,0.8) sd = diag(1,1) 10 20 -r~ 30 - r— 40 i a2ff a2xxmM a2xfmM big small Figure 4.9: Performance of mM averaging strategies: a 2 , 0 = 0.5, 72 = 0.4. 76 a2ff a2xxmM a2xfmM big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.8 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) I 20 10 30 l —r~ 40 _ J _ a2ff a2xxmM a2xfmM big small Mean SE 10 20 30 40 Time Figure 4.10: Performance of mM averaging strategies: a 2 ) 0 = 0.8, 72 = 0. 77 a2ff a2xxmM a2xfmM big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) - r - 10 20 l 30 40 l a2ff a2xxmM a2xfmM big small ~ T - 20 _ l 10 I 30 _ 1 _ 40 Mean SE I 30 i 10 —T- 20 i 40 Mean SE 10 20 30 40 Time Figure 4.11: Performance of m M averaging strategies: a^0 = 0.8, 72 = 0.2. 78 a. cj T3 a2ff a2xxmM a2xfmM big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0.4) sd = (0,0,0) prior prob. on big model: prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean'=. (1,0.8) / s d = diag(1,1) —I— 20 10 _l 30 i I 40 oo d a2ff a2xxmM a2xfmM big small Time Figure 4.12: Performance of mM averaging strategies: a2t0 = 0.8, 72 = 0.4. 79 \ a2ff a2xxwa a2xfwa big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) I 30 L _ 10 I 20 l 40 a2ff a2xxwa a2xfwa big small Mean SE -r— 10 - r - 30 20 Time 40 Figure 4.13: Performance of mWA averaging strategies: a 2 , 0 = 0.2, 72 = 0. 80 Figure 4.14: Performance of mWA averaging strategies: a2>0 = 0.2, j2 = 0.2. 81  # of simulations = 5000 true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.5 prior for big model: mean = (1,0.8,0.4) sd = diag(l,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) I 10 20 30 40 l a2ff a2xxwa a2xfwa big small Time Figure 4.16: Performance of mWA averaging strategies: a2>0 = 0.5, 72 = 0. 83 a2ff a2xxwa a2xfwa big small # of simulations = 5000 true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model = 0.5 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) 10 20 30 l I 40 a2ff a2xxwa a2xfwa big small 30 40 i Mean SE 10 20 l 30 _ l 40 — Mean - SE 10 20 Time 30 40 Figure 4.17: Performance of mWA averaging strategies: a 2 j 0 — 0.5, 72 = 0.2. 84 # of simulations = 5000 true coefficient values: mean = (1,0.8,0.4) sd = (0,0,0) prior prob. on big model = 0.5 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) •prior.for small model: mean = (1;0.8) sd ='diag(1,'l> 10 _i 20 I 30 40 I a2ff a2xxwa a2xfwa big small -r— 20 l 10 30 40 1 20 Time 10 30 40 Figure 4.18: Performance of mWA averaging strategies: a 2 ) 0 = 0.5, 72 = 0.4. 85 # of simulations = 5000 true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.8 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) 10 - r ~ 30 i 20 40 l a2ff a2xxwa a2xfwa big small Figure 4.19: Performance of mWA averaging strategies: a 2 ) 0 = 0.8, 72 = 0. 86 # of simulations = 5000 true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model = 0.8 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag{1,1) —I— 10 _! 20 30 I 40 I a2ff a2xxwa a2xfwa big small Mean S E I 10 20 ~~r~ 30 40 Time Figure 4.20: Performance of mWA averaging strategies: a2,0 = 0.8, 72 = 0.2. 87 # of simulations = 5000 true coefficient values: mean = (1,0.8,0.4) sd = (0,0,0) prior prob. on big model: prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean'=.(1,0.8) . 'sd = diag(1,1 )...••'' 10 20 _ J 30 i - r - 40 L _ 00 ci c i a2ff a2xxwa a2xfwa big small CO - o _ 2x  d CO _ i o o . CM c > 3- U J CL C O 00 ~~ 2 p Time Figure 4.21: Performance of mWA averaging strategies: a 2 ] 0 = 0.8, 72 = 0.4. 88 Chapter 5 Finite Samples: Global Selection The specification of a global selection of ( S ° , S n ) can be approached in a variety of ways. The basic principle that we assume here is that we hope to capture sufficient information to obtain an accurate assessment of risk but simultaneously we do not want the the assessment to be responsive to "bad data". For example, consider the selection of S n . Intuitively, as the size of the sigma field cr(Sn) increases, the risk assessment becomes evermore sen- sitive to the data. Treating /9j(F f e )„ + 1; S n ) as a random variable and using ViPi(Yk,n+i \ S£) as a measure of the sensitivity, the following result formalizes this intuition. L e m m a 5.1 Let S n , T n be two statistics such that T n = #(S n), or equiva- lent^, cr(T n) C o"(Sn). Then pi(Yk >n+i;T^) is not more sensitive to the data than pi(Yktn+1; S£), that is, yiPi(Yk>n+1; T B ) < V i f t ( n > B + 1 ; S B ) . (5.1) 89 Proof Iterating in the definition of pi we see V i P i ( n i n + 1 ; T n ) < since ViPi(Yk,n+1; S n) = V ^ E ^ , P i ( Y k t n + 1 ; S„) + EiVilTppi(Yk,n+1; S£). | This result suggests that the choice S n = Y( n) yields a mongrel risk that is maximally sensitive to the data, while the choice S n = 0 yields a minimally sensitive mongrel risk. In the case where S n is a collection of past predicted predictuals, the inclusion of more predictuals increases sensitivity. 5.1 Closeness to the Bayes solution One type of robustness is that the risk criterion should be insensitive to in- correctly specified models. We illustrate for the case of model choice with two candidate models. In choosing between the two models, we base our decision on the difference in our assessments of the risk of the predictors from the two models. The essential idea is that we would like this difference to be "close" to the true difference in risk irrespective of the data-generator. But because the data-generator is unknown, we use instead the posterior distributions from the Bayes procedure as surrogates. For two models indexed by k and k', we set the target for the difference in risks between the two models to be h = Pk(Yk,n+U Y(n)) - Pk(Yk',n+l', Y(n)) (5.2) 90 when model k is true, or &k> = Pk<(Yk,n+l, Y(„)) - Pfc'(Vfc',n+l') Y( n ) ) (5.3) when model is true. Meanwhile, we select (S°, S£) based on our assessment of the difference in average mongrel risk 5 = p(Yk>n+1; S«, S£) - P{Yk,,n+l-S£, S£). (5.4) Ideally, we would like to select (S°, S£) such that 6 to be close to both Sk and If we are willing to weight the models according to numbers wk and wki, then we might try selecting (S°,S£) to minimize, say, wk{5-8k)2 + wk<(5-5k,)2 (5.5) (pointwise in Y( n ) ) , or, p = wkEk(6 - 6k)2 + wk,Ekl(5 - 8klf. (5.6) (The obvious choice for the weights would be io, = ajj(y(n))-) It can be shown that the minimum for (5.5) can be obtained by setting S" = S£ = Y( n ) so this criterion is not useful for selecting (S£, S^). The criterion (5.6) can be thought of as a robustness criterion. We are trying to select the (S",S£) that keeps 5 close to both 5k and 5k< (in a weighted average sense) regardless of the value of Y( n ) that obtains. Unfortunately, the evaluation of the expectations in p are not tractable analytically and so the criterion is not easily implemented. We attempt to circumvent this problem by finding an approximation p to p on which we base our selection of (S£, S£). 91 5.1.1 Approximating (5.6) Note that both 5 — 8k and 6 — 5ki is a weighted sum of six mongrel risk terms which we can write as, for 5 — 5k say, 6 6 ~ <** = H K3r3 (5-7) i = i where Kj G { ± 1 , ±a ic (S£) , ±aki(S°)} is a weight and Tj is a mongrel risk Pi(-;-). In the normal case, rj can be evaluated using (3.12) through (3.16) and these yield expressions that are quadratic in S n . By substituting for S n by its function of Y( n ) these expressions can be written in the form r; = G + (7; + # Y ( B ) ) 2 (5-8) where Q and jj are numbers and 5j is a vector. L e m m a 5.2 Suppose the weights ai j(S") in (5.7) are replaced by their observed values a j ( s " ) . Then an approximation to ~Ek(5 — 5k)2 is given by 2 Ek = 2tr(A**A* f c)+411^111 + 8 ^ * ^ + 411̂ 111, + (tr(A9k) + + 2i%Z + £ Kj(Q + yj)J (5.9) where £ = Y^jljKjbj, A = A A r , and A is the matrix with columns ,JKJ$J- P r o o f Treating the Kj as numbers instead of random variables and applying (A7), we have r 2>jrj = ( Y + b ) r A ( Y + b) + c (5.10) where b and c are defined in (A9) and (A10). Applying (A4) and (A5), we obtain 92 = 2 t r ( A * f c A * f c ) + 4||i/ f c + b | | i ^ A + (tr(AVk)+ \\vk + m2A + c)2 • The lemma follows from expanding the terms and applying the relation A b = £• I Thus, we can approximate p by p = wkEk + wk'Ek/, (5.11) and then the optimal choice for (S", S£) is (S°,S£)* = arg min p. (5.12) The computations in Lemma 5.2 can be can be avoided in special cases such as the following: Lemma 5.3 Suppose S£ = Y(n) and that the weights aij(S") in (5.7) are replaced by their observed values Oii(s"). Then the optimal choice for S° min- imizes p = wkal,(s«)EkB4 + wk,a2{s«)Ek,B* (5.13) where B2 = ( Y f e > n + i - Yfc', n+i) 2- P r o o f Expanding 5 — 5k using the definition of p with S£ = Y(n), we get 5 — 6k = akpk(Yktn+i,y'(„)) + ak'Pk'(Yk,n+l \ Y ( n ) ) -{®kPk{Yk> i Y ( n ) ) + a*'Afe'(**',n+i; Y ( „ ) ) ) +Pk(Yk,n+l', Y ( n ) ) - Pfc(Vfc',n+i; Y ( n ) ) = aA:Pfc(Vfc,n+i; Y(„)) + ^/(pfc^Yfc'.n+i; Y(„)) + B2) ,n-f l i Y ( n ) ) +Pfc(Tfc,„+i; Y ( n ) ) - (pfc(F/fc,n+i; Y ( n ) ) + B 2 ) = 2ak,B2 (5.14) 93 where the second equality follows from the property pk(Yk',n+i) — Pk(Yk,n+i) + (Xk,n+i — ^ f c ' , n + i ) 2 - Similarly, 8 — 5^ = — 2akB2. The Lemma then follows by treating the weights ak, ak' as numbers so that they can be taken outside the expectations in (5.6). | 5 .1 .2 S i m u l a t i o n R e s u l t s We applied the approximation provided by (5.9) to optimize the choice for Jg, the number of predictuals to be included in S£. (The choice Jg = 0 was omitted in order to avoid certain additional computations not readily available.) We set S£ = Y ( n ) . We call this method of choosing (S£, S n ) the "ROB" strategy. (Note that the optimization is over choices of S n rather than S" as was the case in earlier chapters; the reason is that the needed computations for optimizing over were not readily available.) Figure 5.1 plots each of the two expectations in (5.11) and their weighted sum p as a function of Jg for each of the first 12 sequences at time-point 10 from the scenario 72 = 0.2 and a2t0 = 0.5. The expectation with respect to the full model is much larger by than the one with respect to small model always. One possible explanation for these differences is that, intuitively, large models are more sensitive to data than small models. But we are uncertain whether this explanation can justify such dramatic differences or if there is an error in our computations. At this time, we have not been able to find any coding errors so we shall proceed on the assumption that the computations are correct but are subject to verification. The choice histogram for the scenario with 72 = 0.2 and 012,0 = 0.5 (Figure 5.2) is bimodal and concentrates at the extremes, i.e., there was a tendency to use either very few or almost all of the predictuals. This shape 94 was observed all of the scenarios, albeit with different splits over the two extremes; typically, a larger value of 72 was associated with a greater chance that few predictuals would be selected. Figures 5.3 through 5.11 compare the performance of the ROB strategy based on model averaging (labeled 'a2robwa') relative to the Bayes strategy. Table 5.1 summarizes these comparisons for all of the scenarios. Figures 5.12 through 5.20 compare the performance the the ROB strategy based on model choice (labeled 'c2robwa') relative to the Bayes strategy. Table 5.2 summarizes the comparisons for all of the scenarios. The performance of both the ROB averaging and choice strategies are similar in that the R O B strategy tends to do better than the Bayes strategy when 72 is small (0 or 0.2) but substantially worse when 72 = 0.4 and a2,o = 0.2 or 0.5. Also, the performance of the ROB strategies relative to Bayes is greatest when n is small. 5.2 Equalizing meta-risk Rather than focusing on robustness to model misspecification, we might focus instead on the robustness of the predictors F n + 1 ( S " , S£). One measure of the robustness of Y n + i (S£ , S£) is that its meta-risk should be relatively constant across candidate predictors in an averaged (over Y(„)) sense. This constancy property should hold regardless of which model is true. For instance, when there are only two candidate predictors, we might select the (S" ,S n ) that minimizes maxV* ( P l ( y n + 1 ; T n ) - p2(Yn+1;Tpn).) (5.15) This criterion leads to a solution that conceptually is similar to an averaged (over Y( n)) version of the adaptive minimax criterion. The plots of the meta- 95 risk profiles from Chapter 4 suggest that the minimax solution typically is found where the profile under the big model crosses the profile under the small model (if such a crossing occurs) or at the point where the two profiles are nearest to each other (if no crossing occurs). If these profiles reflected average (over Y( nj) rather than per-sequence properties, then these solutions would correspond typically to the locations where the minimum of (5.15) is achieved. Once again, the difficulty with implementing a criterion based on (5.15) is that the needed expectations cannot be evaluated simply. 5.3 Discussion In both the previous and the current chapters, we have assumed that the same S£ is used in evaluating p f c(Y" f e ) n + 1; S£) for all i. However, it may be desirable to use a different S£ for each model depending on, say, its posterior weight. For example, suppose we are taking a model averaging approach and consider the possible influence of a model that is larger than the true model. The chance of obtaining a poor predictive distribution is relatively high since the predictive distribution would be quite sensitive to the data for this model. The sensitivity would be transferred to the mixture distribution and might result in a poor predictor. Hence, as the posterior weight of a model decreases, it may be beneficial to reduce the number of predictuals in S£ in order to limit the impact of data. Conversely, if the posterior weight increases we feel evermore certain that we have the correct model and so it may be beneficial to increase the number of predictuals in S„ to take more advantage of the data. 96 Table 5.1: Summary comparison of the R O B mongrel averaging strategy to the Bayes strategy. n r 2 0:2,0 72 10 to 20 20 to 30 30 to 40 0 + + + 0.2 0.2 ++ + + 0.4 + - i 0 ++ + + 0.5 0.2 ++ ++ + 0.4 +++ ++ 0 0 + 0 - 0.8 0.2 ++ ++ ++ 0.4 +++ ++ ++ Table 5.2: Summary comparison of the R O B mongrel choice strategy to the Bayes strategy. n r 2 0:2,0 72 10 to 20 20 to 30 30 to 40 0 +++ ++ + 0.2 0.2 ++ ++ + 0.4 0 1 0 +++ ++ + 0.5 0.2 +++ ++ + 0.4 — 0 0 - - 0.8 0.2 0 + + 0.4 0 + + 97 s H - i 1 1 J - 2 4 6 8 J (= # of predictuals in S) J (= # of predictuals in S) 4 6 J (= # of predictuals in S) 2 4 6 8 J (= # of predictuals in S) J (= # of predictuals in S) 6 8 J (= # of predictuals in S) 2 4 6 8 J (= # of predictuals in S) J (= # of predictuals in S) J (= # of predictuals in S) J (= # of predictuals in S) J (= # of predictuals in S) J (= # of predictuals in S) Figure 5.1: Robustness profiles as a function of the number of predictuals included in S£ for ROB strategy: a 2 ) 0 = 0.5, 72 = 0.2. 98 Frequency of selecting given # of predictuals - rob C M d o o d 4 5 6 # of predictuals in S All in d d co 3 C T CO CO C M CO cc o d 1 2 3 4 5 6 7 8 9 10 12 14 16 # of predictuals in S 18 20 22 All 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 All # of predictuals in S Figure 5.2: Histograms of the number of predictuals to include in S£ as selected by ROB strategy: a 2 ] 0 = 0.5, 72 = 0.2. 99 Time Figure 5.3: Performance of ROB averaging strategy: a2,0 = 0.2, 72 = 0. 100 # of simulations = 5000 a2ff a2robwa big small true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(l,1,l) prior for small model: mean = (1,0.8) sd = diag(1,1) I 10 I - r — 30 _ L _ 20 _ L _ 40 _ L _ a2ff a2robwa big small Figure 5.4: Performance of ROB averaging strategy: a2,0 — 0.2, 72 = 0.2. 101 # of simulations = 5000 Time Figure 5.5: Performance of ROB averaging strategy: a 2 > 0 = 0.2, 72 = 0.4. 102 Figure 5.6: Performance of ROB averaging strategy: a2,0 = 0.5, 72 = 0. 103 Figure 5.7: Performance of ROB averaging strategy: a 2 ] 0 = 0.5, 72 = 0.2. 104 Figure 5.8: Performance of ROB averaging strategy: a 2 , 0 = 0.5, 72 = 0.4. 105 Figure 5.9: Performance of ROB averaging strategy: a 2 , 0 = 0.8, 72 = 0. 106 # of simulations = 5000 a2ff a2robwa - - big — small true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model = 0.8 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(l,1) Figure 5.10: Performance of ROB averaging strategy: a2,0 = 0.8, 72 = 0.2. 107 Figure 5.11: Performance of ROB averaging strategy: a 2 > 0 = 0.8, 72 = 0.4. 108 # of simulations = 5000 c2ff c2robwa big small true coefficient values: mean = (1,0.8,0) sd = (0,0,0) prior prob. on big model = 0.2 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) prior for small model: mean = (1,0.8) sd = diag(1,1) 10 20 _ l _ 30 _ 1 _ 40 _ l _ c2ff c2robwa big small m  •- P - T o P -I o —\— 10 20 30 40 Time Figure 5.12: Performance of ROB choice strategy: a 2 , 0 = 0.2, 72 = 0. 109 # of simulations = 5000 true coefficient values: mean = (1,0.8,0.2) sd = (0,0,0) prior prob. on big model = 0.2 10 I 20 I 30 _ 1 _ 40 _ L _ C2ff c2robwa big small 110 Figure 5.14: Performance of ROB choice strategy: a 2, 0 = 0.2, 72 = 0.4. I l l Figure 5.15: Performance of ROB choice strategy: a 2 > 0 = 0.5, 72 = 0. 112 Figure 5.16: Performance of ROB choice strategy: a 2 , 0 = 0.5, 72 = 0.2. 113 # of simulations = 5000 - c2ff - - c2robwa • -• big - - small A , true coefficient values: mean = (1,0.8,0.4) sd = (0,0,0) prior prob. on big model = 0.5 prior for big model: mean = (1,0.8,0.4) sd = diag(1,1,1) /prior for small model: ^\ ; ' / V' 'me , an = ifSifi) 10 20 I 30 I 40 I C2ff c2robwa big small I 10 20 _L_ I 30 I 40 Figure 5.17: Performance of ROB choice strategy: a 2 , 0 = 0.5, 72 = 0.4. 114 Figure 5.18: Performance of ROB choice strategy: a 2 ] 0 = 0.8, 72 = 0. 115 Figure 5.19: Performance of ROB choice strategy: a 2 , 0 = 0.8, 72 = 0.2. 116 Figure 5.20: Performance of ROB choice strategy: a 2 , 0 = 0.8, 72 = 0.4. 117 C h a p t e r 6 A s y m p t o t i c s In this chapter, we establish guidelines for ensuring that the conditioning statistics S£ and S n have good asymptotic properties. Once again the context will be a pair of normal linear models with S° and S n restricted to affine func- tions of the response. We characterize the class of sequences that guarantee that the prediction generated by a model choice or averaging strategy will be asymptotically equivalent to the Bayes predictor from the true model. 6.1 Consistency of Model Weights Regardless of whether a model choice or a model averaging approach is taken, we may wish to require that the model weights, ak, be consistently estimated, that is, if model k, say, were true, then we would like ak —> 1 (weakly or strongly) as more data become available. Otherwise, in the model choice approach, the components of the risk that are computed under the wrong model (which, as such, are suspect) continue to influence the overall assessment even asymptotically. In the model averaging context, consistent estimation of 118 ak ensures that the predictor is derived ultimately from only the correct model. We give a complete proof for weak convergence using Chebyshev's inequality and then outline a Wald-type proof for strong convergence. Let k, k' index normal linear models with model k' nested within model k. Partition the design matrix as Zk,(n) = (Zfc',(n) | Z) where Z consists of the p = Pk — Pk1 covariates that are present in model k but not in model k'. Suppose that Tk is block diagonal with respect to the partitioned parameter vector, i.e., We will establish conditions on U that yield consistent estimation of CKJ for i = k and i = k'. Let Aj(A), A m i n ( A ) , A m a x ( A ) denote respectively, the z-th, the mini- mum, and the maximum eigenvalues of the matrix A . Let | |x | |^ = x T A x . For any non-negative definite matrices A and B of the same dimension, we write A < B iff | | x | | i < | |x| |^ for all vectors x. It will be useful to re-express (3.22) as MK) = n ^ M (6-2) where M = (6.3) is the ratio of the marginal densities for S° from the two models and to define - 2 log M = A - log D (6.4) where D = M = | E - 1 / 2 E , S - 1 / 2 | (6.5) 2-> k' 119 A = | | S S - / i t f | £ - i - | | S « - / i f c (6.6) = WK ~ /**'IIE-I_S-I - 2(S n - p^f^1^ - pk) - \\pk, - /ifc|gj(fi.7) Denote the true parameter value by /30 = (/3*,/3) where /3 = 0 when model k' is true. Let Po = U TZ f c ]( n)/3 0 - U T (zk,,(„)/?: + z )̂ (6.8) denote the expected value of S£ Applying (A4) to (6.6), E j A tr ( E ^ E , ) + | | P o - / ^ l l ' - ^ - J - [tr ( E ^ E . ) + \\p0 - fc' fc' fc 1 fc = t r - S f c x Ei ) + ||/x0 - P f c / l ^ o - ^ | | | - i £ . s ^ i . (6.9) To obtain a bound on the variance of —2 log M , we apply the relation V(A + B) < 2(V(A) + VLB)) to (6.7) and then use (A5) to obtain V i A < 2 = 4 V i ( | | S n - pk>|g-i_s-i + Vi (2 (S n - pk)T^(pkl - MA;)) fc' k t r ((^fc'1 _ ^ f c 1 ) ^ ) 2 + 2\\Po - / i A . / | | ( S - i _ E - i ) E . ( s - . i _ s - i j (6.10) Lemma 6.1 TTie quantities \\p0 - /ifc||2 i , ||yu0 - and | | P f c/ - p,k\\t- i fc -^fc ^ / b are bounded irrespective of whether model k or model k' is true. P r o o f Substituting for E ^ 1 using (3.15) and applying (A3), \Po~Pk\W-i = | | / V - b f c | | Z T n S - i u T Z fc K,(n) fc < Wo - bfc||r-i. fc,(n) (6.11) The expression on the final line involves only constants and hence is bounded. The proofs for the other two quantities follow analogously. | 120 L e m m a 6.2 When model k' is true, \\LL0 — Pk'Wt-1 is bounded. P r o o f Making use of the fact that B = 0, substituting for EjT,1 using (3.6) and applying (A3), l l / / 0 - / i * ' l l ! - i = | | u T z f c , ( n ) ( ^ * - b r o i l s - ! k' k' < l l ^ - b f c H r - L (6.12) The expression on the last line involves only constants and hence is bounded. I L e m m a 6.3 Let G = E ^ , 1 / 2 U T Z = ( U r $ f c , U ) ~ 1 / 2 U T Z and let H = G f G T . The following are equivalent: (*) A m a x ( G T G ) oo (6.13) (ii) A m a x ( H ) - 4 oo (6.14) (m) t r (H)->oo . (6.15) P r o o f (i) (ii): Observe that A m i n ( f ) G G r < H < A m a x ( f ) G G T implies Amin(r') Amax ( G G T ) < A m a x ( H ) < A m a x ( r ) A m a x ( G G T ) . The result then follows from the fact that the non-zero eigenvalues of A A r and A T A are equal for any matrix A . (ii) (iii): Obvious since H has at most p non-zero eigenvalues and tr (H) equals the sum of its eigenvalues. | Theorem 6.1 / / model k! is true, then a necessary and sufficient condition for a.k' —> 1 is that U satisfies (6.13). P r o o f Clearly, ak> ->• 1 ak -> 0 (-21ogM) ->• -oo . It is sufficient to show that D —> oo <̂=>- (6.13) holds while both E f c / A and "VVA are bounded irrespective of whether (6.13) holds. 121 Observe that E f c = Efci + U r Z f Z T U (6.16) implies D = |I + H | . Clearly, I + H has at most p eigenvalues greater than one and the remaining eigenvalues are equal to one. Since the determinant of a matrix equals the product of its eigenvalues, it follows that (1) D > 1, and (2) by Lemma 6.3 D -> oo (6.13) holds. Setting i = k! in (6.9) and simplifying, Ek,A = tr (I - £ f c - ,%) + ||MO - Vk>\\y - \\(io ~ ^ l l | - i E f c # E - i - (6-17) The first term in (6.17) is bounded since tr ( I - E ^ E ^ ) = tr (EK 1 U T Z f Z T U ) = t r ( f 1 / 2 Z T U E f c 1 U T Z f 1 / 2 ) < t r ( f 1 / 2 f f 1 / 2 ) = p where the inequality follows from applying (A3) after substituting for E f c using (6.16). The second term is bounded by Lemma 6.2 and the third term is bounded by Lemma 6.1 since E f c/ < E f c implies \\LI0 - p f c | | 2 1 i , < \\LL0 - Mfcllsr1' k k' k fc Setting i = A;' in (6.10) and simplifying, W A < 4 t r [ ( ! - E f c ' % ) 2 ] + 2||// 0 - A**' 1  CE- 1—s^- 1 >s-1 (x:fc, — 1 ) (6.18) • II M2 The first term is bounded since tr (A2) < (tr (A))2. Since ( E ^ - E ^ E ^ E ^ - E^1) = E,T, 1 / 2(I - (I + H)" 1)E j t", 1 / 2 < E ^ 1 , applying Lemma 6.2 shows that the second term is bounded. Finally, the third term is bounded by Lemma 6.1 and the fact E ^ < E f c . | Theorem 6.2 If model k is true, then (6.13) is a necessary condition for ak -¥ 1. P r o o f Since ak —>• 1 •<=>• (—21ogM) —> oo, it is sufficient to show that both E f c ( -21ogM) and V f c ( - 2 1 o g M ) are bounded when (6.13) fails to hold. This 122 task reduces to showing that both E f c A and VfcA are bounded since it has already been seen that log(D) is bounded if (6.13) fails to hold. So suppose (6.13) does not hold. Setting i = k in (6.9) and simplifying, we have E f c A = tr (H) + ||// 0 - ^ f e , | | | l E t S _ ! - ||/z0 - /Zfclg-i . (6.19) The first term is bounded by supposition and the third term is bounded according to Lemma 6.1. To show that the second term is bounded, let c = 1 + A m a x ( H ) and observe that E ^ E ^ , 1 = E^ 1 / 2 (I + H ) E ^ , 1 / 2 < cE*, 1. Then \\l*>o ~ ^ f c ' l l i r / E k S - 1 < c l l M o - A**'HE- 1 fc' k fc' *:' = c||Z f c/ i ( n)(^* - bfc») + Z / ? | | u s - i u T k' < c\\Zk>,(n)(/3o - bfcOHus-iur + c|| Z / 3 | | u s - i u T k' Kr = - b v | | Z j ; I ( B ) U E 4 - / u r z t , i ( B ) + c | | ^ | | G r G < cll^-b^llr-.+cll^llGTG. (6.20) The second inequality follows from Cauchy-Schwarz and the third one from an application of (A3). Both of the terms in (6.21) are bounded by supposition. Setting i = k in (6.10) and simplifying, we have V f c A < 4 < 4 tr (H 2 ) + 2\\n0 - ̂ ' l l ' s - i ^ - i ^ ^ - i . s - i ) + IK' - tr (H 2 ) + 2\\fi0 - ^ * ' | | E - i S t £ - i + \W ~ Pk\\l 2 - i fc (6.21) where the last inequality follows from the fact (E^,1 - E f c 1 )E f c (E f e , 1 - E/71) < E^EfcEfc, 1. A l l of the terms in in (6.21) are bounded by supposition or have already have been shown to be bounded. | 123 Theorem 6.3 If model k is true, then (6.13) together with the condition that A m a x ( G r G ) < # A m i n ( G r G ) (6.22) for all n and some constant K is sufficient for ak —> 1 • P r o o f It is sufficient to show that (6.13) and (6.22) together imply that E f c ( -21ogM) -> oo and V f c ( - 2 1 o g M ) / E ^ ( - 2 1 o g M ) -> 0 since Chebyshev's inequality then implies (—21ogM) —> oo. Clearly, from (6.19), E f c ( -21ogM) = t r ^ + l l / X o - A t f c ' l l 2 . - ! E t E - i - \\fi0 - »k\\l-i - log |I + H | k' k fc' k > I I ^ - ^ ' I I I - I ^ E - I (6.23) fc' " k' since tr (H) - log |I + H | = £ (Aj(H) - log(l + A {(H)) > 0. We now show that 1  Mo — A t * 'H E - i s f c s - i increases at a rate of at least O ( A m i n ( ( G T G ) 2 ) ) . Since E^EfcSfc,1 = E^, 1 + E f c / U ^ f ^ U E f c , 1 , we have ll _ II2 •> ll _ l l 2 11 Mo Mfe ' l ls -ZEuE-, 1 — H ^ 0 M f c ' l l s - i u ^ z f z ^ u s - , 1 k' " k' fc' fc' = l l ^ l l z T ( i i ) U E - / / = H E - 1 / 2 U ^ Z „ , ( n ) + H ^ l l G T G f G ^ G - 2 ^ T Z f e ' , ( n ) E f c , 1 U T Z f G T G ^ . (6.24) The second term in (6.24) is at least O (XMIN((GTG)2)^ since H ^ l l G T G f G T G - ^min(r)||/3 | |(G TG) 2 > A m i n ( f ) A m i n ( ( G r G ) 2 ) | | / 3 | | 2 . (6.25) In contrast, the first term in (6.24) increases at a rate of at most O ( A m j n ( G T G ) ) since I I A J I 7 T T T r - 1 / 2 H V - 1 / 2 T 7 T 7 — ^ m a x ( G f G T ) 11/3*11̂  T J „ _ i T r T / Zfc,(n) U Lfc' H i V U Z M » ) Z ,fc',(n) U 2 jfc' U Z «=' , (n) < A m a x ( r ) A m a x ( G : r G ) i : C ' < A m a x ( f ) ia m i n (G T G ) i r (6.26) 124 where condition (6.22) has been used in the last inequality. Hence the second term in (6.24) dominates the first term. By the Cauchy-Schwarz inequality, the the first term also dominates the third term. Therefore Efc(—21ogM) —>• oo at a rate of at least O ( A m i n ( ( G T G ) 2 ) ) . To show that Vfc(—21ogM) increases at arate of at most 0 ( A m i n ( ( G T G ) 2 ) , it is sufficient to show that the first term in (6.21) is at most 0 ( A m ; n ( ( G T G ) 2 ) since it has been shown already that the third term in (6.21) is bounded and that the second term is at most O ( A m i n ( ( G T G ) 2 ) ) . But tr (H 2 ) < (tr (H)) 2 < ( A m a x ( f ) t r ( G T G ) ) 2 < ( A m a x ( f ) p A m a x ( G r G ) ) 2 < ( A m a x ( f ) p K \ m m ( G T G ) ) 2 is clearly 0 ( A m i n ( ( G T G ) 2 ) ) . | We conjecture that an alternative simpler condition to (6.22) is that J > p since we believe this condition implies (6.22). (Clearly J > p is a necessary condition for (6.22) since otherwise G T G has a zero eigenvalue.) Another perhaps even more attractive alternative would be to remove (6.22) and instead impose a condition on ( Z ^ , ^ ) , Z ) directly. So long as G T G is full rank, the lower bound obtained in (6.25) obtains only for exceptional sequences of (Zfc'^n), Z ) . For typical sequences, the lower bound is 0 ( A m a x ( ( G r G ) 2 ) so if we exclude the exceptional sequences from consideration, no additional conditions need be added (beyond requiring G T G to be full rank). Condition (6.22) aside, Theorems 6.1 through 6.3 essentially state that weak consistency is characterized by whether or not U satisfies condition (6.13). The remaining material in this subsection gives a Wald-type argument that gives sufficient conditions for strong consistency. Conjecture 6.1 A sufficient condition for oti -4 1 a.s. when model i is true for both i — k and i — k' is that lim rank(U) ->• oo. (6.27) 125 Sketch of proof when model k true: Let Q,ki and £lk denote the parame- ter spaces under models k' and k respectively. We need to show that (6.27) implies that for any e > 0, L = Pp0{M > e) —>• 0 when B0 6 Q,k but j30 g £lki. The density of S£ given the parameter value Bi under model i is P i = (27r)(-^ 2)|a 2U TU|- 1/2exp {~\\s{n) - U T Z , ( n ) A | | 2 ( T 2 U r u ) - 1 } . (6.28) Let the pa denote the density Pi evaluated at the true parameter value B0. Let B((30, 8) denote the open ball centred at B0 with radius 8 such that B(B0,8) (1 Q,k> = 0 and let BC(B0,8) denote its complement. Let Vi((3i) denote the prior density on Bi. Then fo vkt(8k/)Pk>d8ki M = *' (6 29) Snk Vk(Pk)Pkdf3k and L = P0o(M>e) < p f^PBc(MPk PQ \ - 0 O \ p 0 SB(Mvk{Pk)PkdBk ) ( ^ { M P k > e _ n a \ + f t V k { p ) d 0 k < . V Po J \JB(p„,6) e J (6.30) The first inequality follows from recognizing that s u p B c ^ o ^ p k > supQ f c / p k = supnfc)Pfc' > Jnk,uk'{/3k')Pk'dBkl and SB^0,8)uk{Bk)pkdBk < Jnkvk(Bk)Pkd/3k. The second inequality follows from first using the relation E = (E C\ A) U (E n Ac) C A U (E n Ac) where E = {M > e} and and then applying the relation { x y > p } Pi {x < q) C {qy > p ) to E f l Ac and re-arranging terms. Now let W „ = A _ 1 / / 2 S ° where A = a 2 U T U . Then the 126 elements of W „ are independent normals with unit variance and hence pk as a function of W n is given by pk = ( 2 7 r ) - J / V 2 U r U | - V 2 e x P { - i | | w n - A - 1 / 2 U r Z f c ; ( n ) ^ | | 2 } = ( 2 7 r ) - ^ | a 2 U r U | - 1 / 2 e x p | - i f : ( ^ - f t ( / 3 f c ) ) 2 | (6.32) where Wi and gi((3k) are the z-th elements of w„ and A~1^2UTZk^n'j/3k respec- tively. As a function of w„, the expression for pk is not a density due to the additional factor | a 2 U T U | _ 1 / 2 . However this factor appears for each instance of pk and cancel in each term in (6.30). Hence therein pk can be treated as a true density. If gi(Bk) were independent and identically distributed (i.i.d.) for all i , then as J —> oo, the first term of (6.30) converges to 0 by Wolfowitz's (1949) result and the second term is 0(1/J) by a result of Clarke and Barron (1990). Hence if these two results were extended to cover the case of indepen- dent and not identically distributed (i.n.i.d.) data as arising in this problem, then the theorem is proved. The extension from i.i.d. to i.n.i.d. is usually routine although often tedious and involving many conditions. For example, Hoadley (1971) extends Wolfowitz's result by generalizing standard regularity conditions. The extension of the Clarke and Barron work is expected to be similar. | Sketch of proof when mode l k' t rue: Let B(80,81,82) denote an open rect- angular neighborhood of B0 formed as the Cartesian product of B(j30,5\) and B(P0, $2) where B(B0,81) and B(B0,82) are open rectangular neighborhoods of B0 in the subspaces of By and /3 respectively. Let 81 (82) denote the length of each side of of the rectangle. Let v(B) denote the volume of the neighborhood B. Then v(B((30,8U 82)) = v(B({30,81))v(B(P0, <S2)). (6.33) 127 Suppose for the moment that V<5i, 62 there exists a sufficiently large N such that / VkPkdBh > / VkPkdBk (6.34) for all n > N. Then for sufficiently large n, M = > Snk Vk(Pk)PkdBk 1 fB(p0,5lt52) V(B(ZA))D^K' (6.35) v(B(B0,52))fB{,M2) <B£Xh))dpk- As 5i and <52 go to 0, the integrals in the numerator and denominator of (6.35) both converge to pk(/30) and hence M —> oo. So the theorem is proved if it is shown that (6.27) implies (6.34). Roughly speaking, 6.34 is requiring that the density concentrates around B0 as n gets large. If we were dealing with i.i.d. observations then (6.27) would be sufficient so again the key is an extension to i.n.i.d. data. | For a model choice strategy, consistency of the model weights guarantees that the final prediction matches the Bayes predictor from the true model asymptotically. Roughly speaking, if model i is correct (and so —>• 1), then asymptotically p*{k; S £ ) « p(k; i, S £ ) and obviously p(k; i, S n ) is minimized by taking k = i. 6.2 Conditional Predictive Distributions Consistency of model weights in a model averaging procedure does not guar- antee that the final prediction matches the Bayes predictor from the correct model in general. In model averaging the final prediction is derived from 128 (when model i is true) whereas the Bayes predictive distribution is -Fi|Y { n )- For point prediction, these two distributions will generate the same predictor asymptotically if and only if the means p(FI^SPI) and p(Fi\Y(n)) for the two distributions are equal asymptotically, i.e., A ^ M ^ ^ - M ^ s ^ ^ O . (6.36) A generalization of Lemma 3.1 provides one characterization of the solution to 6.36. Theorem 6.4 Asymptotically, the distributions F^SPN andFi\y(n) have the same mean iff the vector ^ ^ Z j ^ T i Z i ^ + i lies in the null space o / U . P r o o f Substituting in (6.36) using (3.8) and (3.20), and simplifying, we get A , = ( Y ( n ) - Z i , ( n ) b i ) r [ * r 1 - \J(\JT%V)-1\JT]Zi,{n)riZi,n+l = ( Y ( n ) - Z i ) ( n ) b 2 ) r ( I - P ) * r i z i i ( B ) r i Z i , n + 1 (6.37) where P is given by (3.25). By the argument used in the proof of Lemma 3.1, A ^ = 0 iff the condition in the statement of the theorem holds. | 129 C h a p t e r 7 D i s c u s s i o n We have proposed a new class of criteria, the mongrel risk, for selecting on- line predictors. The mongrel risk combines both model information and past empirical performance to evaluate the candidate predictors. The application of the mongrel risk requires a rule for selecting the conditioning statistic (S£, S„). The selection can be made in an adaptive or global approach. Our simulations show that an adaptive mongrel approach beats out the Bayes procedure uniformly over a wide range of data-generators in small sam- ples. An analytic proof of this result would, of course, be desirable but it is unlikely to be obtainable given that the expressions needed to assess the per- formance of a mongrel procedure cannot be evaluated analytically. Moreover, because we are dealing strictly with small sample performance, we cannot appeal to the approximation techniques that are often employed in proving asymptotic results. We believe that we have investigated a sufficient variety of simulation parameters to conclude that our results hold in enough generality to be com- pelling. The choices for the coefficient, 72, of X2 cover a practically meaningful 130 range of values and the choices for the prior probability, a 2 , 0 on the full model span a range that seems reasonable for practical applications. We also investigated the impact of the choice of prior distributions on the parameters. For our main simulation work, we set the prior variances, T i & r 2 , on the regression parameters to be identity matrices because we felt that such values reflected the typical (mild to moderate) amount of prior information available in practice. In simulation results not presented here, the use of very weak priors (Ti = 251) gave qualitatively the same results. Our choice of prior means on the parameters may also seem fortuitous in that when 72 is 0 or 0.4, one of the models will have prior means that completely match the data-generator. In practice, we would expect that neither of the prior means would match the means in the data-generator. Hence, we also conducted simulations in which a different set of coefficients for the data-generator was generated for each sequence randomly. That is, we used a random effects model for generating the data. The distribution of the random coefficients was normal with means equal to the values from the fixed coefficients model and 0.21 as the variance. Once again, we found that the results (not presented here) were qualitatively similar to what we found using the fixed coefficients model. As an alternative to using the M S P E (3.34), we evaluated the perfor- mance of the naive mongrel strategies 'a2hf' and 'c2hf' using the difference in relative entropies AD = D(pr\\pB)-D{pr\\pM), (7.1) where pr is the true predictive density (from the data-generator), pB is the density from taking a Bayes strategy, and PM is the density from the mongrel strategy. We limited the investigation to scenarios with a i 2 ) 0 = 0.5 and 72 = 0, 131 0.2, or 0.4. The plots (not shown) of the AD in all three cases exhibited patterns very similar to their analogs under M S P E loss (the third panel in Figures 3.4, 3.5, 3.6 (model averaging) or 3.13, 3.14, 3.15 (model choice)). The mongrel strategies beat out or lost out to the Bayes strategy in the same scenarios and across roughly the same time intervals regardless of whether relative entropy or M S P E loss was used. These results suggest that the gains seen from taking a mongrel approach generalize to loss functions other than M S P E (though this claim is perhaps mitigated by the use of normal models for which relative entropy is closely related to squared error loss). Relative entropy can also be used to select (S", S n ) . Consider the case of model averaging. The candidate predictive distributions are of the form of the mixture distribution (2.19) for different choices of (S", S n ) . Let pm|(s«,s&) denote the density of the mixture distribution. The adequacy of p ( S ° , S n ) as a predictive distribution in relative entropy distance is £)(pr|bm|(sg,s£))- We cannot use this measure to compare different choices of (S°,S£) because the true density px is unknown. Instead, we can compare them based on the relative entropy distance with respect to the posterior density Pi\YM of each candidate model, i.e., for each candidate model i. As in Chapter 4, we might summarize over all models by taking the maximum or taking a weighted average over the A 's . This leads to select (S£, S£) as A = £>(Pi|Y ( n )|bm|(Sg,S{; (7.2) S m M  = arg min max D, (SS,S£) i (7.3) or SmWA = arg min V c ^ A (S2,S£) Y (7.4) 132 for a minimax or weighted average strategy (with weights Ui) respectively. The choice of predictive distribution is then given by p m | s m M o r P m | s m W A - I*1 general, Di cannot be evaluated analytically. However, it's evaluation requires only a one-dimensional integral so numerical integration is feasible. Relative entropy is generally viewed as a "natural" measure of loss so it may be preferred over squared error loss outside of normal models. Both practical and conceptual difficulties have limited our work to only two candidate models. On a practical level, the computational time required to complete a simulation increases roughly as the square of the number of can- didate models so that going from two to three models would more than double the time needed. More importantly, we have yet to develop fully the idea of risk sufficiency so that it is unclear which predictuals should be considered for inclusion in (S", S£) . For instance, should we still use only the predictuals from the smallest model? a combination of predictuals from all but the largest model? or some entirely different set of predictuals? So long as ( S ° , S£) is not risk-sufficient, the generalization of most of our techniques to three or more models is straightforward. But based on our experience with two models, it seems not too difficult to unintentionally achieve risk-sufficiency. We view the development of a greater understanding of risk-sufficiency as important future work. Our arguments for why appropriate use of the mongrel risk criterion generates more accurate predictions than the Bayes procedure have been, for the most part, heuristic. Much additional work is needed to explain fully the properties of the mongrel procedure that give it the advantage. We are currently examining the results for individual sequences to identify the cir- cumstances under which the mongrel approach is better. For example, Figure 7.1 displays a histogram of the difference MSPE(a2ff) - MSPE(a2xfmM) in 133 performance between the Bayes averaging strategy (a2ff) and the minimax mongrel averaging strategy for the 5000 individual sequences. (The predic- tions are for y 1 0 with simulation parameters 72 = 0.2 and a2y0 = 0.5.) The distribution is tightly concentrated about zero but nearly all of the large de- viations are positive. Thus it appears that the mongrel and Bayes strategies typically perform about equally well, but in a small fraction of sequences the mongrel procedure performs much better. This result agrees well with our intuition that the mongrel approach is good at identifying the exceptional se- quences where using all of the data produces misleading risk assessments. Our efforts in this area are ongoing. 134 MSPE(a2ff) - MSPE(a2xxmM) Figure 7.1: Distribution of MSPE(aff) - MSPE(a2xfmM) for predicting Yw with simulation parameters a2,0 = 0.5, 72 = 0.2. (Bottom panel displays the smaller frequency range on an expanded scale.) 135 B i b l i o g r a p h y Allen, D . M . (1974). The relationship between variable selection and predic- tion. Technometrics, 16, 125-127. Berger, J.O. (1985), Statistical Decision Theory and Bayesian Analysis (2nd ed.), New York, Springer-Verlag. Clarke, B.S. and Barron, A .R . (1990). Information-theoretic asymptotics of Bayes methods. IEEE Trans. Inform. Theory, 36, 453-471. Clyde, M . A . (1999). Bayesian model averaging and model search strategies. In J . M . Bernardo, J.O. Berger, A.P. Dawid and A . F . M . Smith (eds), Bayesian Statistics 6, 157-185. Oxford University Press. Dawid, A.P. (1984). Statistical theory: the prequential approach (with discus- sion). J. Roy. Statist. Soc. A, 147, 278-292. Dawid, A.P . (1992). Prequential Data Analysis. In M . Ghosh and P.K. Pathak (eds), Current Issues in Statistical Inference: Essays in Honor of D. Basu, 113-126. IMS Lecture Notes — Monograph Ser. 17. Hayward, C A , Insti- tute of Mathematical Statistics. 136 Hoadley, B. (1971). Asymptotic properties of maximum likelihood estimators for the independent not identically distributed case. Ann. Math. Statist., 42, 1977-1991. Hoeting, J.A., Madigan, D., Raftery, A . E . and Volinsky, C T . (1999). Bayesian model averaging: A tutorial (with discussion). Statist. Sci., 14, 382-417. Raftery, A . E . , Madigan, D. and Hoeting, J .A. (1997). Bayesian model aver- aging for linear regression models. J. Am. Statist. Assoc., 92, 179-191. Seillier-Moiseiwitsch, F. and Dawid, A.P. (1993). On testing the validity of sequential probability forecasts. J. Am. Statist. Assoc., 88, 355-359. Skouras, K . (1998). On the Optimal Performace of Forecasting Systems: The Prequential Approach. Ph.D. Thesis, Department of statistical Science, University College London. Wolfowitz, J. (1949). On Wald's proof of the consistency of the maximum likelihood estimate. Ann. Math. Statist., 20, 601-602 137 Appendix A . 1 Let A and D be nonsingular matrices. Then ( A + B D B T ) _ 1 = A - 1 - A - 1 B ( D - 1 + B r A _ 1 B ) _ 1 B r A _ 1 . (Al) Rearranging this equality and left and right multiplying by A , we obtain B ( D - 1 + B r A _ 1 B ) _ 1 B T + A ( A " 1 + B D B r ) _ 1 A = A (A2) so that, if A and D are p.d., then each term in (A2) is n.n.d. and hence B ( D " 1 + B T A - 1 B ) ~ 1 B r < A . (A3) A . 2 Let Y ~ jV(m, S). Then for any n.n.d. matrix A E ( Y T A Y ) = tr (AS) + m T A m , (A4) V ( Y r A Y ) = 2tr ( A S A S ) + 4 m r A S A m . (A5) Note that for any n.n.d. matrix B , t r ( B 2 ) = YJIUJ^IJ < HiUjhi^jj = ( t r B ) 2 where is the i, j-th element of B . Setting B = A 1 / 2 S A 1 / 2 , we have t r ( A S A S ) = t r ( B 2 ) < ( t r B ) 2 = (tr A S ) 2 and hence V ( Y r A Y ) < 2(tr A S ) 2 + 4 m T A S A m . (A6) 138 A.3 Let pi = Ci+ (7* + 8fYy, where Ci, 7i are scalars and r5i5 Y are n-vectors. If r < n, then for any scalars i/j, £ ^ft = (Y + b) T A (Y + b) + c (A7) i=i where A = ^UiSiSf, (A8) i=l r b = A-J2nvrfi (A9) i=i c = E ^ ( C i + 7 ' ) - b r A b (A10) i=l and A~ is a generalized inverse of A. Note also that Ab = £ [ = 1 7^(5;. 139

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 2 1
China 1 57
City Views Downloads
Ashburn 2 0
Beijing 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items