Bootstrapping Methods for Confidence Intervals in Statistics, Lecture notes of Statistics

Download Bootstrapping Methods for Confidence Intervals in Statistics and more Lecture notes Statistics in PDF only on Docsity! Paper PK02 A Practical Introduction to the Bootstrap Using the SAS System Nancy Barker, Oxford Pharmaceutical Sciences, Wallingford, UK ABSTRACT Discovering new medications is a field populated by many unknowns. Even when the human physiology is well understood it can be difficult to predict the way in which the body will react to a new drug in order to assess what the effects will be. Moreover, knowledge of the corresponding distributions for any derived measurements may be limited. In today's world of clinical trials it is not only important to know how a novel product performed, it is necessary to give some indication of the accuracy of any estimate of the performance. Without knowledge of the distribution, standard parametric techniques cannot be reliably executed and so an alternative is required. Many conventional statistical methods of analysis make assumptions about normality, including correlation, regression, t tests, and analysis of variance. When these assumptions are violated, such methods may fail. Bootstrapping, a data-based simulation method, is steadily becom ing more popular as a statistical methodology. It is intended to simplify the calculation of statistical inferences, sometimes in situations where no analytical answer can be obtained. As computer processors become faster and more powerful, the time and effort required for bootstrapping decreases to levels where it becomes a viable alternative to standard parametric techniques. Although the SAS/STAT® software does not have any specific bootstrapping procedures, the SAS® system may be used to perform bootstrap methodology. Even with the speed of modern computers, careful use of efficient programming techniques is required in order to keep processing time to a minimum. The author will attempt to address programming techniques for bootstrapping methodology. This will include an introduction to the techniques of bootstrapping – including the calculation of standard errors and confidence intervals with the associated SAS code. This paper also compares and contrasts three different methods of calculating bootstrap confidence intervals. INTRODUCTION Discovering new medications is a field populated by many unknowns. Even when the human physiology is well understood it can be difficult to predict the way in which the body will react to a new drug in order to assess what the effects will be. Moreover, knowledge of the corresponding distributions for any derived measurements may be limited. In the field of statistics, there are lots of methods that are practically guaranteed to work well if the data are approximately normally distributed and if all we are interested in are linear combinations of these normally distributed variables. In fact, if our sample sizes are large enough we can use the central limit theorem which tells us that we would expect means to converge on normality so we do not even need to have samples from a normal distribution as N increases. So if we have two groups of say 100 subjects each and we are interested in mean change from baseline of a variable then we have no need to worry and can apply standard statistical methods with only the most basic of checks for statistical validity. However, what happens if this is not the case? Suppose we want to make inferences about the data when one of the following is true: - Small sample sizes where the assumption of normality does not hold - A non-linear combination of variables (e.g. a ratio) - A location statistic other than the mean Bootstrapping, a data-based simulation method for assigning measures of accuracy to statistical estimates, can be used to produce inferences such as confidence intervals without knowing the type of distribution from which a sample has been taken. The method is very computationally intensive, so it is only with the age of modern computers that it has been a viable technique. The use of statistics in pharmaceutical research is becoming more and more sophisticated. It is increasingly common for proposed methodology to go beyond standard parametric analyses. In addition, cost data – with its extremely non-normal distribution is regularly collected. Bootstrapping methodology has become a recognised technique for dealing with these issues. In a recent CPMP guidance document “Guidance on Clinical Trials in Small Populations” (released for consultation in March 20051), it is stated that “…some forms of Bootstrap methods make no assumptions about data distributions and so can be considered a ‘safe’ option when there are too few data to test or verify model assumptions … they prove particularly useful where very limited sample data are available…”. The aim of this paper is to introduce the reader to the process of bootstrapping and to look in a little more detail at two of the more common applications of bootstrapping: estimating standard error (SE) and estimating confidence intervals (CI) for statistics of interest. It should be noted that while all of the necessary calculations and SAS code has been included, a great deal of the statistical theory has been glossed over. Readers interested in understanding fully the statistical theory involved should read Efron and Tibshirani(1993)2. WHAT IS BOOTSTRAPPING? The method of bootstrapping was first introduced by Efron as a method to derive the estimate of standard error of an arbitrary estimator. Finding the standard error of an estimate is an important activity for every statistician as it is rarely enough to find a point estimate; we always want to know how reasonable is the estimate – what is the variability of the estimator? In fact, sometimes statisticians are really greedy and not only want to know the point estimate and its standard error but also things like the bias or even the complete distribution of the estimator. If available these can be used to create confidence intervals or to test statistical hypotheses around our estimator. The use of the term 'bootstrap' comes from the phrase “To pull oneself up by one's bootstraps ” - generally interpreted as succeeding in spite of limited resources. This phrase comes from the adventures of Baron Muchausen - Raspe (1786)3 In one of his many adventures, Baron Munchausen had fallen to the bottom of a lake and just as he was about to succumb to his fate he thought to pick himself up by his own bootstraps! The method is extremely powerful and Efron once mentioned that he considered calling it ‘The Shotgun’ since it can “… blow the head of any problem if the statistician can stand the resulting mess”. The quotation relates to the bootstrap’s wide applicability in combination with the large amount of data that results from its application together with the large volume of numerical computation that is required. Figure 1 : Illustration of the Bootstrap Method 202 202 -206 100 -206 188 280 -206 280 280 Original sample 196, -12, 280, 212, 52, 100, -206, 188, -100, 202 Mean = 91.2 100 100 188 280 -100 -100 188 52 188 196 … … 188 202 100 100 212 52 212 -12 -12 202 -100 -206 52 202 280 196 -206 -12 188 -100 Sample 1 Sample 2 Sample 3 Sample B Mean=109.2 Mean=124.4 Mean=91.4 Mean=29.4 with the same number of records. We can then go on to perform the analysis required. Note that it is important to use a different seed each time you create a sample in the RANUNI( ) function or you will create the same sample each time. EFFICIENT PROGRAMMING When performing bootstrapping operations it is extremely important to use efficient programming techniques. If we were to create each bootstrap sample individually, then perform the analysis required (PROC MEANS in the case of the example above), and append the results to the previous samples we require a lot longer for our programming. Whenever possible , it is much better to use BY processing. If we were to create all our bootstrap samples simultaneously and then used PROC MEANS with a BY or CLASS statement, the difference in times can be quite startling. Creating the bootstrap samples simultaneously is simply a matter of using an additional DO statement in the DATA STEP as follows: data bootsamp; do sampnum = 1 to 1000; /* To create 1000 bootstrap replications */ do i = 1 to nobs; x = round(ranuni(0) * nobs); set original nobs = nobs point = x; output; end; end; stop; run; This will create a data set with 1000*no. observations in ORIGINAL. Each bootstrap sample will be identified using the variable SAMPNUM which can be used to perform later BY processing. Using the first method of creating each sample and finding its mean individually and appending the result to the data set containing the means from the previously calculated samples took the author’s home computer a total of 2 minutes and 41 seconds. Using the second method of creating all samples simultaneously then using PROC MEANS with a CLASS statement took approximately 1 second. From this we can see that the optimal method of running code is to create all samples within a single data set and then use by processing (with the sample number as the BY variable). However, this is not always possible. You may be running an analysis that requires some form of coding that has no BY processing (e.g. SAS/IML) it may be that your initial data set is so large that your computer could not deal with running everything at once. In such a case the author would recommend closing your log before running (ODS LOG CLOSE;) as this gives a slight time saving. Running the code in this manner reduced the time from 2 minutes and 41 seconds to 1 minute and 50 seconds. BOOTSTRAPPING FOR TREATMENT GROUP COMPARISONS In the majority of clinical work, we are looking to compare two or more treatments. To do this we often use randomisation to assign the subjects to each treatment. Suppose we have a trial with two treatment groups. If we were to use the bootstrapping methods to randomly select subjects we would have many bootstrap samples where the balance of subjects in each treatment group is different from the original sample. Instead, we need to create bootstrap samples which have the same number of observations in each group as are found in the original. To do this we re-sample independently for each treatment group and then reset the observations together. Consider our original example: let us suppose that we had another 10 subjects in the study who were on a different treatment. Suppose that the first set of data came from subjects treated with ‘active’ (referred to as Active) medication while a further 10 subjects came from a group on ‘placebo’ (referred to as Placebo) medication (with values = 120, -80, -63, 200, 23, 54, -198, 165, -8, 19). This time, we wish to find the difference between the means of the two groups. As discussed above, we create bootstrap samples which each have 10 Active values and 10 Placebo values, by independently sampling from each group. We then calculate the mean for each group within each sample and find the difference between the means (our statistic of interest) for each sample. We can then look at the distribution of this statistic in exactly the same way as we did for the ‘single group’ mean above. To do this with SAS we need to split the data before re-sampling. Unfortunately, the POINT= option does not work when the WHERE= option is also used in the set statement, so two DATA steps are required: the first to divide the datasets for each treatment group and the second to perform the bootstrap sampling. If the original data set (ORIGINAL) contains two variables called CHANGE (for change from Baseline value) and TRT (treatment: 1=Active, 2=Placebo) – the following code will perform the desired analysis. %macro bootse(b); data orig1 (where=(trt=1)) orig2 (where=(trt=2)); /* Create one data set for each treatment */ set original; run; data boot; %do t=1 %to 2; /* Create independent sets of replications */ do sample=1 to &b; do i = 1 to nobs; pt = round(ranuni(&t)*nobs) ; set orig&t nobs = nobs point=pt; /*use treatment-specific data */ output; end; end; %end; stop; run; proc means /* find the mean, by sample and treatment */ data=boot noprint nway; class sample trt; var change; output out=x mean=mean; run; data diffmean; merge x (where=(trt=1) rename=(mean=mean1)) x (where=(trt=2) rename=(mean=mean2)); by sample; diffmean=mean1-mean2; /* calculate the difference between the means */ run; proc means data=diffmean std; var diffmean; output out=bootse std=bootse; run; %mend; CALCULATING CONFIDENCE INTERVALS USING THE BOOTSTRAP So far we have looked at the computation of bootstrap standard errors. Standard errors are often used to assign approximate confidence intervals to a parameter θ . For example, given an estimate θ̂ and an estimated standard error esˆ , the usual 95% confidence interval is: Mean ± 1.96*SE. = esˆ*96.1ˆ ±θ This is based on the assumption that: Z = es) θθ −ˆ ∼N(0, 1) (i.e. is approximately distributed to the Normal distribution with mean=0 and sd=1) This is valid as n→ ∞, but is only an approximation for finite samples. A better approximation would perhaps use the t distribution and assume that (for a sample size of n): Z = es) θθ −ˆ ∼ 1−nt (i.e. is approximately distributed to the t-distribution with n-1 degrees of freedom) For small samples this is a better approximation, but the use of the t distribution doesn’t adjust the confidence interval to account for the skewness in the underlying population – or other errors that can result when θ̂ is not the sample mean. Through the use of bootstrapping we can obtain accurate intervals without having to make normal theory assumptions. There are a variety of methods available which compute confidence intervals and we will go into a little detail on each. In general, the variability between estimators for confidence intervals is greater than for standard errors and so it is better to have a larger number of bootstrap replications. The 100 replications recommended to calculate standard errors is not really acceptable for confidence limits. This author would generally recommend a set of 10,000 replications. The sampling error and number of replications required will not be discussed here. Efron and Tibshirani1 (Chapter 19) gives an excellent discussion of this for those wishing for more detail. For each of the following types of confidence intervals, we will demonstrate the SAS code required, assuming that we are trying to calculate a 95% confidence interval for the difference between the means for the example data above. In all cases, the %BOOTSE macro will be run and will create the data sets BOOT, DIFFMEANS and BOOTSE are in the section above. THE BOOTSTRAP-T INTERVAL This method allows us to estimate the distribution of our statistic Z = es) θθ −ˆ directly from the data. Suppose we create a set of bootstrap samples. For each, we estimate a value of Z as follows: ( ) ( ) ( )bes b b *ˆ ˆ*ˆ * θθ − =Ζ where ( )b*θ̂ is the estimate of θ̂ for the b th bootstrap sample and ( )bes *ˆ is the estimated standard error of *θ̂ for the bth bootstrap sample. The α th percenti le of ( )b*Ζ is estimated as the value ( )αt̂ such that: ( ) ( ){ } α α = Β ≤Ζ tb ˆ*# where #{x} denotes the number of observations which meet condition x The confidence interval is then calculated as: ( ) ( )( )estest ˆˆˆ,ˆˆˆ 1 ×−×− − αα θθ data ci_perc; set diffmean end=eof; /* pull out 250th and 9750th values as confidence limits */ retain conf_lo conf_hi; if _n_= 250 then conf_lo=diffmean; if _n_=9750 then conf_hi=diffmean; if eof then output; keep conf_lo conf_hi; run; A print out of the final data set ci_t would give us the following: conf_lo conf_hi -26 208.7 Giving us a confidence limit of (-26.0, 208.7) for the difference in the means between the groups. BETTER CONFIDENCE INTERVALS One of the most important goals of bootstrap methodology is that we produce confidence intervals which closely match ‘exact’ confidence intervals in the situations where statistical theory yields an exact answer. Unfortunately, neither the bootstrap-t method, nor the percentile method defined above meet this goal. The bootstrap-t intervals have good theoretical coverage – but tend to be erratic in practice. In addition, if no statistical formulae are available to estimate the standard error for the statistic of interest, a nested bootstrap procedure is required which can increase the required number of samples (and hence required computational time) exponentially. The percentile intervals are simpler to use and less erratic but comparison with exact statistical methods have demonstrated these to have less than satisfactory coverage in some cases. Part of the problem with these confidence intervals is that they are unable to deal with the issues related to: 1. the bootstrap estimates are biased with respect to the original estimate 2. the standard error varies with the value of the estimate An improved version of the percentile method has been formulated to help deal with these problems. The bias -corrected and accelerated (BCa) intervals have been shown to have a substantial improvement over the percentile method in both theory and practice. The method does however require a little knowledge of bias and jackknifing techniques before the computation can be performed. BIAS In statistics, bias is defined as “any systematic failure of a sample to represent its population base”. In the case of bootstrapping, we want to see if there is a median bias – i.e. a systematic underestimation or overestimation of our sample statistic when calculated by the bootstrap replications. To do this we look at the proportion of samples which calculate a sample statistic which is below the statistic calculated by the original sample. Let θ̂ = the actual statistic calculated on the original records ( )b*θ̂ = the statistic calculated from the b th bootstrap sample B = the number of bootstrap samples Then the bias correction 0ẑ = ( ){ }         <Φ − B b θθ ˆˆ# * 1 where { }condition# = number of samples to meet the specified condition and ( )x1−Φ indicates the inverse function of a standard normal distribution (PROBIT function in SAS) e.g. ( )975.01−Φ = 1.96. To calculate this in SAS for our example for difference between the means we would use the following code: %bootse(10000); data bootorig; /* set original and sample data together so can calculate */ set original (in=a) /* mean values in one procedure */ boot; if a then sample=0; run; proc means data=bootorig noprint nway; class sample trt; var change; output out=bootmean /* mean for each sample */ mean=mean; run; data diffmean (keep=sample diffmean) /* data set containing bootstrap values */ bias (keep=bias); /* data set containing bias correction value */ merge bootmean (where=(trt=1) rename=(mean=mean1)) bootmean (where=(trt=2) rename=(mean=mean2)) end=eof; by sample; diffmean=mean1-mean2; /* calculate the difference between the means */ retain origdiff; /* add value of original estimate */ if sample=0 then origdiff=diffmean; /* difference to each record */ if diffmean lt origdiff then lessthan=1; /*flag if bootstrap sample gives lower */ else lessthan=0; /*value than original sample */ retain nless 0; if sample gt 0 then nless=nless+lessthan; /* count samples with flag lessthan */ if sample ne 0 then output diffmean; /* output only bootstrap sample statistics */ if eof then do; /* for the last value calculate: */ propless=nless/sample; /* 1. proportion of values below original estimate */ bias=probit(propless); /* 2. inverse normal of that proportion */ output bias; /* 3. output only that record to new data set */ end; run; For our data, 54.5% of bootstrap samples gave differences between the means which were less than that found by the original sample. This gives us a value for the bias of 0.1135. JACKKNIFING AND ACCELERATION The acceleration statistic refers to the rate of change of the standard error of the estimate of our sample statistic with respect to the true value of the statistic. Efron and Tibshirani proposed that this be calculated using jackknife values of the estimate of the statistic . Jackknifing looks at how much each individual record influences the estimate. To do this we recalculate the estimate of the LED with all but one record - removing each record in turn. For our example we would create a set of (10+10)=20 estimates for the difference between the mean which are each based on 19 of the records with a different record being removed each time. Let ( )iθ̂ = the statistic of interest calculated on the jackknifed sample with the ith record removed n = the number of jackknifed samples i.e. the number of records across all treatment groups ( )•θ̂ = the mean of the n jackknife samples = ∑ = n i n 1 θ̂ Then the acceleration â = ( ) ( )( ) ( ) ( )( ) 2 3 1 2 3 1 ˆˆ6 ˆˆ       − − ∑ ∑ = • = • n i i n i i θθ θθ This statistic may look a little daunting, but is easy to calculate. We simply create our jackknife samples, find the difference between the mean of the samples and each sample and calculate both a squared value and a cubed value of our difference. Then we sum these squared and cubed differences and plug the values into the equation. In SAS, for our example, we would use the following code. data origjack; /* create a new data set which contains observation */ set original end=eof; /* numbers 1 to &nobs (no. obs in data set) */ obsnum=_n_; if eof then call symput('nobs', put(obsnum, 2.)); run; %macro jackdata; /* use macro for %do processing utility */ data jackdata; set %do i=1 %to &nobs; /* do loop to create all samples */ origjack (in = in&i where=(obsnum ne &i)) /* remove a different value each time */ %end;; %do i=1 %to &nobs; if in&i then repeat=&i; /* add repeat number for each sample */ %end; run; %mend; %jackdata; proc means data=jackdata noprint nway; class repeat trt; var change; output out=jacksum /* mean for each sample */ mean=mean; run; A print out of the final data set ci_bca would give us the following: conf_lo conf_hi -17.5 216.4 Giving us a confidence limit of (-17.5, 216.4 ) for the difference in the means between the groups. COMPARISON OF CONFIDENCE INTERVALS We now have three different ways of calculating bootstrap confidence intervals. But how do they compare? The author has stated that the BCa method is thought to address some of the identified issues with the bootstrap-t method and the percentile method – but how can we be sure that this is the better method? To show the differences, the author created a simulated data set with 100 data points from a log-normal distribution, X. This distribution has the property that if we took the Y=log(X), then the values would come from a normal distribution. In this case, we chose the distribution such that log(X) = Y ∼N(2, 4) i.e. so that the log of the values has mean(µ)=2 and variance(σ2)=4 (and so standard deviation(σ)=2). The X values take the variable name CHANGE (to simulate a change from baseline value) and so the Y values are LOGCHG. The simulation was performed using the RANNOR() function which is generates random values from a N(0, 1) distribution. To do the random sampling within SAS you would use the following code: data original; do pid=1 to 100; /* create 100 observations */ logchg=(rannor(0)*2 + 2); /* simulate LOGCHG as normal with mean=2 and sd=2 */ change = exp(logchg); /* if LOG(CHANGE)=LOGCHG then EXP(LOGCHG)=CHANGE */ output; end; run; The log-normal distribution is a nice one to use as an example as it is very highly skewed and a confidence interval calculated with normal theory methods would be highly suspect. However, some statistical theory has been formulated to give a reliable estimate for the confidence interval for the mean (modified-Cox method - see Zhao and Gou (1997)2). The modified-Cox method for calculating a confidence interval for the mean was then used to compare the three methods for bootstrap confidence intervals (10,000 iterations) with the following results: Table 2: Comparison of bootstrap-confidence intervals for mean of log-normal data Method Confidence Interval Modified-Cox (25.5, 92.3) Bootstrap-t (23.9, 137.2) Percentile (18.8, 71.6) BCa (24.4, 101.7) If we take the modified-Cox as the ‘true’ value for the confidence interval, we can see that the BCa method is by far the closest to this true value. The bootstrap-t method is relatively close on the lower limit of the confidence interval but is almost 50% larger at the upper limit than the modified-Cox value. The percentile interval gives limits which are too low at both ends – most importantly at the upper end showing that it does not cover the true interval at all. In comparison, the BCa method covers the interval given by the modified-Cox method and is only slightly more conservative at both ends. A CAUTIONARY NOTE In the non-bootstrap case, a parametric estimate accepts the possibility that a sample can be non-representative by random occurrence. A bootstrap depends on the sample being representative. If the sample is not representative, the conclusions will be completely inappropriate. It is very important to be confident that your sample is a good representation of the whole population from which you are sampling. CONCLUSION The bootstrap method can be generalised into these three steps 1. Create the bootstrap replications by re-sampling with replacement (replicating independently with two (or more) samples , such as for different treatment groups) 2. Calculate the statistic of interest for each sample 3. Use the distribution for the replicated statistic to make inferences. Despite being a relatively new methodology, the bootstrap is quickly becoming a regular part of the clinical statistician’s arsenal. It allows the statistician to assess the statistical accuracy of complicated procedures, by exploiting the power of the computer. Its aim is to carry out familiar statistical calculations such as standard errors, biases and confidence interval in an unfamiliar way – by purely computational means rather than through the use of statistical formulae. The bootstrap relieves the statistician from having to do complex mathematical derivations, or in some instances , provides an answer where no analytical solution is possible. It allows us to assess the properties of a statistical procedure for the data at hand. The bootstrap is a fairly crude form of inference but it can be used when the statistician is unable or unwilling to carry out extensive modelling. For large samples, it gives accurate answers no matter what the underlying population. Using bootstrapping techniques within the SAS system is relatively easy. The RANUNI function together with the POINT= option within the SET statement of a DATA STEP makes the creation of bootstrap samples a simple process. The BY processing allowed within most of the procedures available in SAS allows us to utilise efficient programming techniques, which can cut down on the computational time required. Although we have only shown a brief introduction here, the world of the bootstrap is large and expanding. As well as calculating standard errors and confidence intervals, bootstrapping methodology is available for regression, hypothesis testing, adaptive estimation, calibration and bioequivalence techniques. REFERENCES 1. Committee for Medicinal products for Human Use (CPMP), Guideline for Clinical Trials in Small Populations (Draft). Release for Consultation: March 2005. http://www.emea.eu.int/pdfs/human/ewp/8356105en.pdf 2. Efron and Tibshirani (1993), An Introduction to the Bootstrap, Chapman &Hall/CRC (IBSN: 0-412-04231-2). 3. Rudolf Erich Ras pe (1786), The Singular Travels, Campaigns and Adventures of Baron Munchausen. 4. Zhou, X-H., and Gao, S. (1997), “Confidence intervals for the log-normal mean,” Statistics in Medicine , 16, 783-790. ACKNOWLEDGEMENTS With many thanks to Katherine Hutchinson at OPS for time, encouragement and a thorough review, RECOMMENDED READING An Introduction to the Bootstrap: Efron and Tibshirani, Chapman &Hall/CRC 1993 (IBSN: 0 -412-04231-2). This book gives an excellent introduction to the methodology of the bootstrap. It is very readable and gives good examples throughout. CONTACT INFORMATION Your comments and questions are valued and encouraged. Contact the author at: Nancy Barker Oxford Pharmaceutical Sciences Ltd The Stables, 114 Preston Crowmarsh, Wallingford OXON OX10 6SL UK Work Phone: +44 (0)1491 821680 Fax: +44 (0)8704 580729 Email: Nancy.Barker@ops-web.com SAS and all other SAS Institute Inc. product or service names are registered trademarks or trademarks of SAS Institute Inc. in the USA and other countries. ® indicates USA registration. Other brand and product names are trademarks of their respective companies.