Document Sample

```					                                     Lecture 13

1. SAS version 9.2

2. Sampling distribution for the average

3. Arguments and output of bootstrap macros

4. Tuberculosis testing example

5. Bootstrap of correlation

6. Bootstrap of kappa

7. Bootstrap of proportions of positive and negative agreement

1

SAS version 9.2

See http://support.sas.com/documentation/whatsnew/index.html
for a description of changes: essentially none in statistical procedures.

Update when you have time to solve any problems that arise.

2
Review: Sampling distribution for the average

The bootstrap is based on the idea of taking many samples and computing a
statistic from each sample, in order to estimate the distribution of the statistic.

One of the confusing ideas in the introductory statistics course is the sampling
distribution of the average.

Confusing because a distribution means a histogram of sample averages,

while in real life, we have just one sample, we calculate just one average.

That’s not enough to make a histogram.

3

If we have a ﬁnite population, such as a deck of cards, then it’s easy to take more
samples.

Suppose this is the population histogram, with µ = 6 2 and σ = 2.7:
3

For a sample of 5 cards, should we sample with or without replacement?

4
Histograms of 100 averages of 100 simple random samples of sizes 2, 5, and 10:

Sample size = 2         Sample size = 5                    Sample size = 10

2   :   00              2   :                              2   :
2   :   5               2   :                              2   :
3   :   00              3   :   2                          3   :
3   :   55              3   :   68                         3   :
4   :   00000           4   :   0004                       4   :
4   :   5555            4   :   66668                      4   :   67
5   :   00000           5   :   0024444444                 5   :   344
5   :   55555555555     5   :   6688888                    5   :   55555555666677788889
6   :   000000000000    6   :   00000222222444444          6   :   0000011112224
6   :   55555           6   :   66666666688888             6   :   5555555666677778888888999
7   :   0000000000      7   :   0000002222222444           7   :   001111122222333444444
7   :   55555555555     7   :   6666666688                 7   :   566668899
8   :   0000000000      8   :   00002222222                8   :   0012223
8   :   555555          8   :   68                         8   :
9   :   0000000000      9   :   0                          9   :
9   :   5555            9   :                              9   :

5

¯
How does the size of the sample n affect the distribution of the average x n ?

Sample Size n           2          5             10     population

2
Mean of histogram           6.6        6.5            6.6       µ = 63

SD of histogram         1.8        1.24          0.85       σ = 2.7

σ   n     1.91       1.21          0.854

6
The answer is the Central Limit Theorem: as the sample size n gets large,

¯
the sampling distribution of the average x n gets closer to a Normal distribution

with mean equal to the population-mean µ and SD equal to σ               n,

where the population-SD is σ.

How is the bootstrap different from this procedure?

7

Apply bootstrap to draw samples of n = 10 and ﬁnd the average of each.

Approximate    Approximate
Approximate          Lower          Upper
Observed   Bootstrap     Standard         Confidence     Confidence    Confidence
Name   Statistic      Mean        Error             Limit          Limit        Level (%)
mean    6.66667      6.6593      0.84264             5.01513        8.31821          95

Sample Size n         2            5       10        population

Mean of histogram         6.6          6.5      6.6         µ = 62
3

SD of histogram        1.8          1.24    0.85         σ = 2.7

σ      n   1.91          1.21    0.854

8
%macro analyze(data=,out=);
proc means mean noprint data=&data; compute sample average
var x;
output out=&out mean=mean;
%bystmt;
run;
data &out;
set &out;
drop _type_ _freq_ ; remove extra stuff
run;
%mend analyze;

%boot(data=pop, size=10, samples = 2000, biascorr=0,chart=0)
%bootci(percentile, alpha=.05)

9

Bootstrap algorithm

2. Select B independent bootstrap samples x∗1, x∗2, . . . , x∗B , each a simple random
sample of size n drawn with replacement from the data x.

3. Calculate the statistic of interest θ (theta) from each bootstrap sample:
ˆ∗
θb = s(x∗b ),

for b = 1, . . . , B. These are called bootstrap replications.
ˆ∗
4. Use the B bootstrap replications {θb } as an estimate of the sampling
ˆ
distribution of θ for the population that gave the original data.
ˆ
Estimate the standard error of θ by the sample standard deviation of the
bootstrap replications.
Use the 2.5 and 97.5 percentiles to estimate a 95% conﬁdence interval.

10
Bootstrap macros

The main macro is %boot which calls the starting data set and calls %analyze
(which we must write) to calculate the statistic(s) of interest.

%boot produces two output temporary datasets in the WORK directory:

• bootdist contains the statistics calculated from each bootstrap sample. It is
good practice to make a histogram of these.

• bootdata contains all the bootstrap samples, indexed by the variable
_sample_

The second macro is %bootci which uses bootdist to compute conﬁdence
interval by different methods.

11

Arguments of the BOOT macro

%******************************* BOOT *******************************;
%macro boot(      /* Bootstrap resampling analysis */
data =,         /* Input data set, not a view or a tape file. */
the starting data set
samples =200,      /* Number of resamples to generate. */
200 is default, but use 1000 or more for conﬁdence intervals
residual=,         /* Name of variable in the input data set that
contains residuals; may not be used with SIZE= */
equation=,         /* Equation (in the form of an assignment statement)
for computing the response variable */
size=,             /* Size of each resample; default is size of the
input data set. The SIZE= argument may not be
used with BALANCED=1 or with a nonblank value
for RESIDUAL= */
balanced=,         /* 1 for balanced resampling; 0 for uniform
resampling. By default, balanced resampling
is used unless the SIZE= argument is specified,
in which case uniform resampling is used. */
random =0,           /* Seed for pseudorandom numbers. */
seed allows exact duplication of bootstrap results
stat=_numeric_,/* Numeric variables in the OUT= data set created
by the %ANALYZE macro that contain the values
of statistics for which you want to compute
12
bootstrap distributions. */
id=,              /* One or more numeric or character variables that
uniquely identify the observations of the OUT=
data set within each BY group. No ID variables
are needed if the OUT= data set has only one
observation per BY group.
The ID variables may not be named _TYPE_, _NAME_,
or _STAT_ */
biascorr =1,       /* 1 for bias correction; 0 otherwise */
alpha=.05,        /* significance (i.e., one minus confidence) level
for confidence intervals; blank to suppress normal
confidence intervals */
print=1,          /* 1 to print the bootstrap estimates;
0 otherwise. */
chart =1          /* 1 to chart the bootstrap resampling distributions;
0 otherwise. */
);

13

biascorr = 0

ˆ
Don’t use bias correction. If bias is the bootstrap estimate of bias in the estimate θ,
then the obvious bias-corrected estimator is
ˆ    ˆ
θc = θ − bias

ˆ                      ˆ
“Bias correction can be dangerous in practice. Even if θc is less biased than θ, it
may have substantially greater standard error.

ˆ                           ˆ      ˆ
If bias is small compared to the estimated SE(θ), then it is safer to use θ than θc .

ˆ
If bias is large compared to the estimated SE(θ), then it may be an indication that
ˆ
the statistic s(x) = θ is not an appropriate estimate of the parameter θ.”

(from Efron and Tibshirani, p 138)

14
chart=0           Proc Chart is an old procedure to make histograms.

Frequency
|                                      **
|                                      **
|                                      **
|                                      **
350   +                                      **
|                                      **
|                                      **
|                                      **
|                                      **
300   +                                      **      **
|                                      **      **
|                                      **      **
|                                      **      **
|                                      **      **
250   +                                      ** ** **
|                              **      ** ** **
|                              **      ** ** **
|                              **      ** ** **
|                              ** ** ** ** **
200   +                              ** ** ** ** **
|                              ** ** ** ** **
|                              ** ** ** ** **
|                              ** ** ** ** **
|                              ** ** ** ** ** **
150   +                              ** ** ** ** ** **
|                              ** ** ** ** ** **
|                              ** ** ** ** ** ** **
|                              ** ** ** ** ** ** **
|                              ** ** ** ** ** ** **
100   +                          ** ** ** ** ** ** ** **
|                      ** ** ** ** ** ** ** ** **
|                      ** ** ** ** ** ** ** ** **
|                      ** ** ** ** ** ** ** ** **
|                      ** ** ** ** ** ** ** ** ** **
50    +                      ** ** ** ** ** ** ** ** ** **
|                      ** ** ** ** ** ** ** ** ** ** **
|                      ** ** ** ** ** ** ** ** ** ** **
|              ** ** ** ** ** ** ** ** ** ** ** ** **
|              ** ** ** ** ** ** ** ** ** ** ** ** ** **
--------------------------------------------------------------------------
4   5   5   5   5   6   6   6   6   7   7   7   7   8   8   8   8   9
.   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
8   1   3   6   8   1   3   6   8   1   3   6   8   1   3   6   8   1
7   2   7   2   7   2   7   2   7   2   7   2   7   2   7   2   7   2
5   5   5   5   5   5   5   5   5   5   5   5   5   5   5   5   5   5

15

Instead, use Proc Univariate to make a good histogram:

Goptions reset=all vsize=5in ftext=simplex ctext=black lfactor=1.5;
proc univariate data = bootdist ;
var statistic ;
histogram   statistic / cfill=graye03;
inset q1 median mean q3 / position=NW noframe;

16
Arguments of the BOOTCI macro

%******************************* BOOTCI *******************************;
%macro bootci(    /* Bootstrap percentile-based confidence intervals.
Creates output data set BOOTCI. */
method,        /* One of the following methods must be specified:
PERCENTILE or PCTL
HYBRID
T
BC
BCA       Requires the %JACK macro
*/
stat=,         /* Numeric variables in the OUT= data set created
by the %ANALYZE macro that contain the values
of statistics for which you want to compute
bootstrap distributions. */
student=,      /* For the T method only, numeric variables in the
OUT= data set created by the %ANALYZE macro that
contain the standard errors of the statistics for which
you want to compute bootstrap distributions.
There must be a one-to-one between the VAR=
variables and the STUDENT= variables */
id=,           /* One or more numeric or character variables that
uniquely identify the observations of the OUT=
data set within each BY group. No ID variables
17

are needed if the OUT= data set has only one
observation per BY group.
The ID variables may not be named _TYPE_, _NAME_,
or _STAT_ */
alpha=.05 ,      /* significance (i.e., one minus confidence) level
for confidence intervals */
print=1);      /* 1 to print the bootstrap confidence intervals;
0 otherwise. */

Note that the BCA (bias-corrected accelerated) method calls another macro,
%JACK (for jackknife). This is included in the ﬁle of bootstrap macros.

18
Example: testing for tuberculosis

The standard method for testing a person for exposure to tuberculosis is the
tuberculin skin test (TST). A small amount of a test protein is injected under the
skin and 48 hours later the diameter of the red bump (induration) is measured.
A diameter > 10mm is considered a positive test, indicating the person has been
exposed to tuberculosis.

A second test, Quantiferon, uses a blood sample that is tested by exposure to
tuberulosis antigens. The response is the amount of INF-γ produced above a
control level, and responses ≥ 0.35 IU/ml are considered a positive test.

As part of an MPH project, the Quantiferon test was used in addition to TST in
screening 195 immigrants, in order to compare the two tests.
(Data from C. Baker.)

19

Correlation: Pearson or Spearman?
20
Bootstrap conﬁdence interval for the correlation coefﬁcient

We will start by looking at the correlation between the Quantiferon measurement
and the TST induration diameter.

Spearman Correlation Coefficients
Prob > |r| under H0: Rho=0
Number of Observations

Quantiferon                    TST

Quantiferon             1.00000                0.66569
<.0001
198                 195

For correlation, the nearly universal practice is to give a point estimate r with a
p-value, because that is what packages produce.
The p-value comes from H0 : ρ = 0 which is rarely relevant—we want to know how
large the correlation is, not whether it’s zero.

We want a conﬁdence interval.
21

Writing %analyze for a correlation

We need to get Proc Corr to produce a data set. This old procedure has a special
way to code the output statistics, with these options to the Proc Corr statement:

Output data set with Pearson correlation statistics               OUTP=
Output data set with Spearman correlation statistics               OUTS=

To get the Spearman correlation coefﬁcient out:

proc corr noprint data=a outs =check;
var Quantiferon TST;

proc print data=check;

Obs     _TYPE_     _NAME_             Quantiferon              TST

1       MEAN                                4.518         10.308
2       STD                                 6.778          9.652
3       N                                 198.000        195.000
4       CORR      Quantiferon               1.000          0.666
5       CORR      TST                       0.666          1.000

22
We need to get rid of the extra variables and the extra observations, so that
%analyze makes a dataset with only one observation: the Spearman correlation
between TST and Quantiferon.

%macro analyze(data=,out=);
proc corr noprint data=&data outs=&out;
var Quantiferon TST;
%bystmt;
run;
data &out;
set &out;
if ( _type_=’CORR’ and _name_=’TST’);
spearman_corr=Quantiferon;
drop    _type_ _name_ Quantiferon TST;
run;
%mend analyze;

Use IF to choose rows (observations), DROP to choose columns.

23

%include   "C: ... bootstrapmacros.sas";
%boot(data=a, samples = 1500, biascorr=0, chart=0)
%bootci(percentile, alpha=.05)

24
from BOOT                                            Approximate Approximate
Approximate             Lower       Upper
Observed Bootstrap   Standard            Confidence Confidence Confidence
Name     Statistic    Mean      Error                Limit       Limit    Level (%)
spearman_corr 0.66569   0.66367    0.045927             0.57568     0.75571      95

Method for              Minimum       Maximum
Confidence             Resampled     Resampled      Number of
Name          Interval               Estimate      Estimate      Resamples
spearman_corr Bootstrap Normal           0.48481       0.79199          1500

from BOOTCI(percentile)
Approximate         Approximate
Lower               Upper
Observed    Confidence          Confidence   Confidence
Name               Statistic      Limit               Limit       Level (%)
spearman_corr           0.66569      0.56874             0.75010          95

Method for         Number of
Name               Confidence Interval     Resamples
spearman_corr          Bootstrap percentile        1500

Spearman correlation between Quantiferon measurement and the TST induration
diameter was r = .67 (bootstrap 95% conﬁdence interval: .57–.75).

25

These are diagnostic tests, and correlation does not measure how well they agree
in identifying people exposed to tuberculosis. Lines show diagnostic cut-off values
(above cut-off is positive).

26
TST
Positive          Negative

Quantiferon    Positive        A                 B               R1

Negative         C                 D               R2

C1                C2              N

One measure of agreement is kappa κ, the amount of agreement beyond what
would be expected by chance:

observed agreement PO = (A + D) N

agreement expected by chance PC = (R 1C 1/N + R 2C 2/N ) N

P O − PC
kappa is deﬁned κ =
1 − PC

27

Proc Freq data=baker.qft;
tables QuantiferonTB * TST_Read / nopercent;
exact agree ;

QuantiferonTB(QuantiferonTB)

Frequency|
Row Pct |
Col Pct |         0|      1|       Total
---------+--------+--------+
0 |     67 |     23 |              90
| 74.44 | 25.56 |
| 77.01 | 21.30 |
---------+--------+--------+
1 |     20 |     85 |         105
| 19.05 | 80.95 |
| 22.99 | 78.70 |
---------+--------+--------+
Total          87      108           195

28
Simple Kappa Coefficient
--------------------------------
Kappa (K)                 0.5553
ASE                       0.0598         ASE = asymptotic standard error
95% Lower Conf Limit      0.4381
95% Upper Conf Limit      0.6725

Test of H0: Kappa = 0

ASE under H0                 0.0716
Z                            7.7579
One-sided Pr > Z             <.0001
Two-sided Pr > |Z|           <.0001

Exact Test
One-sided Pr >= K        2.905E-15
Two-sided Pr >= |K|      3.660E-15

Instead of relying on a large-sample-theory calculation for the conﬁdence interval,
we can use the bootstrap.

29

Writing %analyze for kappa

Proc Freq noprint data=baker.qft;
exact agree;
output out=B kappa ;

Obs    _KAPPA_      E_KAPPA      L_KAPPA           U_KAPPA    E0_KAPPA         Z_KAPPA    PL_KAPPA

1     0.55529      0.059784        0.43812         0.67246    0.071577         7.75794        .

Obs     PR_KAPPA        P2_KAPPA        XPL_KAPP         XPR_KAPP            XP2_KAPP

1     4.3299E-15      8.6597E-15             .         2.9045E-15       3.6599E-15

30
%macro analyze(data=,out=);
proc freq noprint data=&data;
exact agree;
output out=&out kappa ;
%bystmt;
run;
data &out;
set &out ;
drop E_KAPPA   L_KAPPA   U_KAPPA    E0_KAPPA   Z_KAPPA   PL_KAPPA
PR_KAPPA      P2_KAPPA     XPL_KAPP      XPR_KAPP      XP2_KAPP;
run;
%mend analyze;

%boot(data=baker.qft, samples = 2000, biascorr=0, chart=0);
%bootci(bca, alpha=.05);

31

32
Approximate   Approximate
Approximate         Lower         Upper
Observed   Bootstrap      Standard        Confidence    Confidence   Confidence
Name  Statistic      Mean         Error            Limit         Limit       Level (%)
_KAPPA_ 0.55529     0.55394       0.059870           0.43795       0.67263         95

Method for    Minimum   Maximum
Confidence   Resampled Resampled Number of
Name       Interval     Estimate Estimate Resamples LABEL OF FORMER VARIABLE
_KAPPA_ Bootstrap Normal 0.34386   0.75264     2000   Simple Kappa Coefficient

Approximate Approximate
Lower       Upper                         Method for
Observed Confidence Confidence Confidence             Confidence   Number of
Name  Statistic    Limit       Limit    Level (%)             Interval    Resamples
_KAPPA_ 0.55529    0.43093     0.66538        95               Bootstrap bca    2000

ASE was 0.0598, and 95% CI was 0.4381—0.6725. This is a large sample (n = 195) so
large-sample estimates should do well here.

33

TST
Positive          Negative

Quantiferon    Positive           A                 B           R1

Negative            C                 D           R2

C1                C2          N

Two other common measures are proportion of positive agreement and
proportion of negative agreement.

2A
proportion of positive agreement =
R 1 +C 1

2D
proportion of negative agreement =
R 2 +C 2

34
Bootstrap CIs for positive and negative agreement

SAS does not compute the proportions of positive and negative agreement. In our
%analyze macro, we use Proc Freq to get the counts in the 2 × 2 table, and then
compute positive and negative agreement proportions.

Proc Freq data=baker.qft;
tables QuantiferonTB * TST_Read/ out=A1              ;   writes table to dataset
exact agree;

35

QuantiferonTB(QuantiferonTB)
Col Pct |        0|       1|
---------+--------+--------+
0 |     67 |     23 |
---------+--------+--------+
1 |     20 |     85 |
---------+--------+--------+

Quantiferon

1          .              .          4       .
2          0              .          3       .
3          0              0         67     34.3590
4          0              1         23     11.7949
5          1              0         20     10.2564
6          1              1         85     43.5897

36
%macro analyze(data=,out=);
proc freq noprint data=&data;
tables QuantiferonTB * TST_Read /out=&out ;
%bystmt;
run;
data &out;
set &out;
retain a b c d p_pos p_neg;
if QuantiferonTB=0 and TST_Read=0 then a=count;
if QuantiferonTB=0 and TST_Read=1 then b=count;
if QuantiferonTB=1 and TST_Read=0 then c=count;
if (QuantiferonTB=1 and TST_Read=1) then do;
d=count;
p_pos=2*d/(2*d+b+c);
p_neg=2*a/(2*a+b+c);
output;
end;
drop a b c d QuantiferonTB TST_Read COUNT PERCENT;   2 numeric variables
run;
%mend analyze;

%boot(data=baker.qft,samples = 2000, biascorr=0,chart=0);
%bootci(bca, alpha=.05);

37

38
39

Approximate     Approximate
Lower           Upper                       Method for
Observed     Confidence      Confidence     Confidence     Confidence
Name     Statistic       Limit           Limit         Level (%)      Interval
p_neg     0.75706       0.67901         0.82353           95        Bootstrap bca
p_pos     0.79812       0.73118         0.85470             95        Bootstrap bca

Report the observed statistic(s) from the sample with the conﬁdence interval
(or SE) from the bootstrap:

The proportion of positive agreement between Quantiferon measurement and the
TST induration diameter was 80% (bootstrap 95% conﬁdence interval: 73–85%)
and the proportion of negative agreement was 76% (CI: 68–82%).

40

```
DOCUMENT INFO
Shared By:
Categories:
Stats:
 views: 33 posted: 3/12/2010 language: English pages: 20
How are you planning on using Docstoc?