/* suggest the following for standard use by others */ options ls=78 formdlim='-' nocenter nodate nonumber nosource2; /* use following for running examples */ *options ls=78 formdlim='-' nocenter nodate nonumber symbolgen mlogic mprint nosource2; /* use following for debugging */ *options ls=78 formdlim='-' nocenter symbolgen mlogic mprint source2; %let UnifyPow = ~robrien/UnifyPow.new/UnifyPow980823.sas; *-----------------------------------------------------------------------; * 23 August 1998 ; * --Examples: MODULE Version-- ; * This input conforms to the "module" (non-macro) version of UnifyPow. ; * Another file available from the UnifyPow website has input conforming ; * to the macro version. ; *-----------------------------------------------------------------------; *-----------------------------------------------------------------------; * UnifyPow website: ; * http://www.bio.ri.ccf.org/Power ; *-----------------------------------------------------------------------; *-----------------------------------------------------------------------; * --Notes-- ; * The first example is used in most workshops and one-hour talks. ; * The rest of the workshop examples follows. The last few examples ; * are for the March 1998 presentation at SUGI 23 in Nashville and a ; * new paper on power for generalized linear models, co-authored ; * with Gwowen Shieh, presented in August 1998 at the Joint Statistical ; * Meetings in Dallas. ; * ; * The UnifyPow website has several other files, including Acrobat ; * (.pdf) files of the workshop notes, the latest update of the SUGI 23 ; * proceedings paper, the manuscript for the O'Brien-Muller (1993) ; * chapter, and the handout for the Shieh-O'Brien JSM '98 invited paper ; * on power for generalized linear models. ; *-----------------------------------------------------------------------; title1 "UnifyPow Workshop and SUGI 23 Presentation"; *-----------------------------------------------------------------------; * Testing single proportion (sign test) ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Mystic Michelle: Ordinary (pi = .50) or Gifted (pi = .80)?"; title3 "Testing a single proportion"; datalines; pi .80 . Mystic Michelle claims 80% accuracy. null .50 . (This is default--not really needed.) Ntotal 20 40 alpha .01 .05 %tables; %include "&UnifyPow"; datalines; /# Same problem, but now find minimum N to achieve specified power at given alphas. #/ pi .80 power .99 .995 alpha .005 tails 1 %tables; title1 "UnifyPow Workshop"; *-----------------------------------------------------------------------; * Testing single proportion (custom null) ; *-----------------------------------------------------------------------; title2 "DCA for Lactic Acidosis in Children with Severe Malaria"; title3 "Uncontrolled pilot study, comparing mortality to 28%"; %include "&UnifyPow"; datalines; pi .14 . DCA cuts mortality 50%. null .28 . estimated mortality rate from epi studies alpha .05 .20 . high alpha OK for pilot study Ntotal 40 60 100 %tables; *-----------------------------------------------------------------------; * Comparing 2 independent proportions ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Placebo vs. DCA for Lactic Acidosis in Children with Malaria"; title3 "28% die untreated. What if DCA cuts this by 25%?"; datalines; /# o Actual study of children with _severe_ malaria. o 28% base mortality rate supported by meta-analysis of published surveys. o 25% reduction in mortality supported by animal-model study. o Small pilot study had 2/10 deaths in each group. o 2:1 randomization to DCA pleases local health officials. o One interim look using alpha = .01 & final analysis at .045. o Able to randomize 1500. Should interim look come at 750 or 999? #/ pi .28 .21 weight 1 2 alpha .01 .045 Ntotal 750 999 1500 %tables; *-----------------------------------------------------; * Testing Goodness of Fit to Multinomial Distribution ; *-----------------------------------------------------; %include "&UnifyPow"; title2 "Do the Biostat Boys use fair dice?"; datalines; /# A trio of shady characters called the Biostat Boys run a craps game at lunch near the FDA. One of the regular patrons, Lucky Luke, has been losing lately and wonders if the dice are fair. The known distribution for a fair pair of dice is: X: 2 3 4 5 6 7 8 9 10 11 12 p(X): .028 .056 .083 .111 .139 .166 .139 .111 .083 .056 .028 Luke believes that one dice is weighted so that the "1" comes up 1/4 of the time, not 1/6. So he computes that the distribution may be: X: 2 3 4 5 6 7 8 9 10 11 12 p(X): .042 .067 .092 .116 .142 .166 .125 .100 .075 .050 .025 How many tosses must he observe in order to have 90% power (alpha = .05) to find such a discrepancy, based on the full distribution? #/ GoodnessOfFit .042 .067 .092 .116 .142 .166 .125 .100 .075 .050 .025 Null .028 .056 .083 .111 .139 .166 .139 .111 .083 .056 .028 alpha .05 .10 power .80 .90 %tables; %include "&UnifyPow"; datalines; /# What if Lucky Luke only counts how many "1"s appear on a throw? The distributions are X: 0 1 2 fair p(X): .694 .278 .028 biased p(X): .625 .333 .042 #/ GoodnessOfFit .625 .333 .042 Null .694 .278 .028 alpha .05 .10 power .80 .90 %tables; *----------------------------------------------------; * Testing association in an R x C contingency table ; *----------------------------------------------------; %include "&UnifyPow"; title2 "Variation in sarcoma type by region"; /* Stimulated by example found in Section 3.4 of Freeman DH (1987). Applied Categorical Data Analysis, New York: Marcel Dekker. */ datalines; /# There are 3 different types of soft-tissue sarcomas of the arms and legs: o Fibroid o Lipoid o Mixed (or other) Incidence data can be obtained from good cancer registries. The Question: Do the relative proportions of these 3 types differ among 4 geographic regions? Scenario- SARCOMA SOFT TISSUE TYPE % OF TOTAL ------------------------- SAMPLE REGION Fibroid Lipoid Mixed ---------- A .50 .20 .30 30% B .60 .25 .15 20% C .35 .35 .30 25% D .45 .20 .35 25% How large must Ntotal be to have 90% or 95% power? #/ 2WayContTable .50 .20 .30 > .60 .25 .15 > .35 .35 .30 > .45 .20 .35 weight .30 .20 .25 .25 alpha .01 .05 power .90 .95 %tables; *-----------------------------------------------------------------------; * Ordinary 2-group t test ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Dr. Alalgia's Two-Group Design"; title3 "Traditional t-Tests"; datalines; /# Scenario represents reductions of 45% vs. 25%. #/ mu -.86 -.42 SD .45 .57 .65 alpha .05 .01 weight 2 1 NTotal 15 21 33 %tables; %include "&UnifyPow"; datalines; mu -.86 -.42 SD .45 .57 .65 alpha .05 .01 weight 2 1 power .80 .90 %tables; *----------------------------------------------; * Matched-pairs t test ; *----------------------------------------------; %include "&UnifyPow"; title2 "Dr. Alalgia's Matched-Pairs Design"; title3 "Traditional t-Based Tests"; /* The datalines4 statement allows UnifyPow's comment blocks to contain semicolons. But this also requires that UnifyPow statements end with ;;;; */ datalines4; /# ---- NOTE ---- UnifyPow requires: mu1 mu2 ........ at least one pair of means; see 2-group example below: AB/BA cross-over design SD1 SD2 ........ exactly one pair of SDs Corr(Y1,Y2) .... at least one correlation For each value of Corr(Y1,Y2), UnifyPow computes SD of diff = Y1 - Y2. Call this: SDdiff[Corr(Y1,Y2)]. The SDMult variable (m) allows the user to further vary SDdiff[Corr(Y1,Y2)], i.e. SDdiff[Corr(Y1,Y2), m] = m*SDdiff[Corr(Y1,Y2)]. Default is m = 1.00 only. #/ PairedMu -.86 -.42 SD .60 .80 . give exactly 2 SDs corr .5 .6 SDMult 1 1.2 . may give several, default = 1.0 alpha .05 .01 TotalPairs 17 25 . Any word with 'total' is OK ;;;; %tables; *-----------------------------------------------------------------------; * 2-group, Wilcoxon-Mann-Whitney (using mu, SD) ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Dr. Alalgia's Two-Group Design "; title3 "Wilcoxon Tests Based on mu & SD, Assuming Logistic Parent"; datalines; mu -.86 -.42 SD .45 .57 .65 alpha .05 .01 weight 2 1 NTotal 15 21 33 Wilcoxon parent logistic methods all %tables; *-----------------------------------------------------------------------; * 2-group, Wilcoxon-Mann-Whitney (using p1) ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Dr. Alalgia's Two-Group Design "; title3 "Wilcoxon Test, Specifying p1 Directly, Assuming Logistic Parent"; datalines; 2Wilcoxon .75 SD 1.00 1.25 1.40 weight 2 1 NTotal 15 21 33 methods Noether parent logistic %tables; *-----------------------------------------------------------------------; * Matched-pairs design, Wilcoxon signed-rank ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Dr. Alalgia's Matched-Pairs Design"; title3 "One-Group Wilcoxon Tests, Assuming Normal Parent"; datalines; PairedMu -.86 -.42 SD .60 .80 corr .5 .6 SDMult 1 1.2 TotalPairs 17 25 Wilcoxon %tables; *-----------------------------------------------------------------------; * Wilcoxon-Mann-Whitney ; * 2 groups, ordered categorical outcome ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Dr. Alalgia's 2-Group Design"; title3 "7-pt Likert Scale for Improvement in Headache"; /* Much Somewhat A Little No A Little Somewhat Much Worse Worse Worse Change Better Better Better -3 -2 -1 0 1 2 3 */ datalines; 2Wilcoxon .03 .04 .08 .10 .27 .28 .20 .10 .10 .15 .25 .20 .15 .05 weight 2 1 power .90 .95 alpha .05 .01 %tables; *-----------------------------------------------------------------------; * Wilcoxon signed-rank ; * 1 group, interval-level categorical outcome ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Dr. Alalgia's Matched-Pairs Design/ 7-pt Likert scale"; title3 "Ordered Categorical Outcome, 1-Group Wilcoxon Test"; /* D = Difference of pair of 7-pt Likert scores, scored Y = {-3, -2, -1, 0, 1, 2, 3}. Thus, D can range from -3 - 3 = -6 (worst outcome favoring ET therapy) to 3 - (-3) = +6 (best outcome favoring ET therapy). */ datalines; 1Wilcoxon .02 .03 .05 .05 .05 .06 .09 .14 .14 .14 .09 .08 .06 Limits -6 6 Null 0 Ntotal 50 75 100 alpha .05 method Noether %tables; *-----------------------------------------------------------------------; * One-way ANOVA with complex contrasts ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Dr. Mindy Bowdy: 4-Group Design"; title3 "Example 3 in O'Brien and Muller (1993, Section 8.2.3)"; title4 "(1) Dominators, (2) Ordinaries, (3) Loners, (4) Friendlies"; datalines; mu .35 .50 .52 .60 weight .20 .50 .10 .20 SD .16 .19 alpha .05 Ntotal 60 80 100 contrasts "Ord vs Loners" 0 1 -1 0 "Almost Overall (2 DF)" 1 -.833 -.167 0 > 0 -.833 -.167 1 ; /* Temporarily store .05 results */ data PowData1; set PowData; %include "&UnifyPow"; title2 "Dr. Mindy Bowdy: 4-Group Design"; title3 "(1) Dominators, (2) Ordinaries, (3) Loners, (4) Friendlies"; datalines; mu .35 .50 .52 .60 weight .20 .50 .10 .20 SD .16 .19 alpha .0167 .05 Ntotal 60 80 100 NoOverall contrasts "{Ord & Loners} vs Fr" 0 -.833 -.167 1 "Dom vs {Ord & Loners}" 1 -.833 -.167 0 "Fr vs Dom" -1 0 0 1 ; /* Concatenate the two sets of results and table them altogether. */ data PowData; set PowData1 PowData; %tables; *-----------------------------------------------------------------------; * 2-group (AB/BA) cross-over design via t tests on differences ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Asthma (Active) Treaments A vs. B, Outcome is FEV1"; title4 "Two-Group (AB/BA) Cross-Over Design with Continuous Outcome"; /* This follows Patel (1983) example (2.1) in Jones B, Kenward MG (1989, Design and Analysis of Cross-Over Trials, Chapman and Hall). Two active drugs, patients have acute bronchial asthma, design is ordinary AB/BA crossover, outcome measure is FEV1. Scenario: Mean FEV1 Drug A Drug B | ------------------------ AB | mu11 = 1.6 mu12 = 1.7 Order | BA | mu21 = 1.9 mu22 = 2.3 SD(A) = .60, SD(B) = .75, SDMult = {1.0, 1.2}, corr(A,B) = {.60, .70, 80} */ datalines; PairedMu 1.6 1.7 > 1.9 2.3 SD .60 .75 Corr .6 .7 .8 . conjectures for corr(A, B) SDMult 1.0 1.2 . gives 100% and 120% of base SDdiff sides 2 TotalPairs 50 NoOverall contrast "Drug x Order" 1 -1 "Drug A vs. Drug B" 1 1 %tables; *-----------------------------------------------------------------------; * Blackwelder's one-sided "equivalency" testing strategy ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Blackwelder's one-sided 'equivalency' testing strategy"; title3; datalines4; /# Abciximab plus a low-dose of heparin reduces complications and increases 30-day survival rates in patients undergoing high-risk coronary angioplasty or other kinds of revascularization (EPILOG Invetigators, 1997 [N Engl J Med, 336:1689]). But abciximab is expensive. Giving bivalirudin (Bittl et. al., 1995 [N Engl J Med, 333:764-9]) albeit with provisional ("bail-out") use of abciximab may give equivalent efficacy and safety at lower costs. Design. Bivalirudin plus provisional abciximab (B+provA) vs. abciximab plus a low-dose of heparin (A+lowH). Randomization ratio is 2:1, that is, 2/3 get B+provA. Primary end-point. As per EPILOG trial: "death by any cause, myocardial infarction, or urgent revascularization within 30 days." Statistical considerations. This equivalency trial will be handled using Blackwelder's method (1982 [Controlled Clinical Trials, 3:345-53]), which tests Ho: pi(B+provA) - pi(A+lowH) GE delta [B+provA appreciably worse] Ha: pi(B+provA) - pi(A+lowH) < delta [B+provA not appreciably worse; could be better] Delta sets the region of equivalency, which we take to be within 25% of the 30-day event rate of 0.052 found in the EPILOG trial. Thus, take delta = 0.012 < 0.013 = 0.25*0.052. Our sample-size analysis will assume that pi(B+provA) = pi(A+lowH) = 0.052. We seek a power of 0.90 for an alpha = .05 test of this one-tailed hypothesis. We shall use either the traditional Pearson ("approximate") unconditional test or its exact counterpart (Suissa & Shuster, 1985 [J Royal Stat Soc A, 48:317-27]), whichever is more powerful. [Even though their true alpha-levels are practically identical in large studies, the two methods will have different powers when the group pi conjectures are unequal and the sample-sizes are not balanced. See O'Brien and Muller (1993).] #/ pi .052 .052 weight 2 1 null .012 Ntotal 8001 9999 12000 14001 tails 1 ;;;; /* store results to later build custom table */ data StorePow; set PowData; %include "&UnifyPow"; datalines4; /# Sponsor actually believes that bivaliridin may cut primary events by at least 5%. That is, they think that pi(B+provA) = .95*.052 = .0494. #/ pi .0494 .052 weight 2 1 null .012 Ntotal 8001 9999 12000 14001 tails 1 ;;;; /* merge results with first set */ data PowData; set StorePow PowData; /* build custom table */ proc tabulate format = 6.3 order=data; class alpha effctitl NullValu NTotal ProbStmt; var power; table ProbStmt="True State" * NullValu="Delta" * effctitl = "Method", alpha = "Alpha" * NTotal ="Total N" * power="Power"*mean=" " /rtspace=45; *-----------------------------------------------------------------------; * McNemar's test of 2 correlated proportions ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Thromboembolism (case/control) & The Pill ('outcome')"; title3 "Matched Pairs, Binary Outcome"; /* Example follows study by Sartwell PE, et. al. (1969, Thromoembolism and oral contraceptives: an epidemiologic case-control study. Am J Epi, 90: 365-380.) Cases: Women of childbearing age who have some kind of thromoembolism. Matched Controls: One per case--treated in the same hosptal, matched on age, number of prior pregnancies, income level. Scenario: Matched Control Using The Pill? No Yes | ------------------------ Thromboembolism Case No | p00 = .54 p01 = .08 Using The Pill? | Yes | p10 = .32 p11 = .06 */ datalines; McNemar .32 .08 Ntotal 50 to 150 by 25 %tables; *-----------------------------------------------------------------------; * Testing a single Pearson correlation ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Jean Netticks: ARX Enzyme Levels in Father and Son"; title3 "Long been assumed that rho[ARX(dad), ARX(son)] < .55."; title4 "New genetic theory predicts a correlation of rho > .70."; datalines; rho .70 null .55 Ntotal 50 75 100 tails 1 %tables; *-----------------------------------------------------------------------; * Comparing two Pearson correlations ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "rho[ARX(dad), ARX(son)]:"; title3 "Is there a difference between obese and normal fathers?"; title4 "Obese: rho = .55 Normal: rho = .70"; datalines; rho .55 .70 weight 1 3 Ntotal 400 to 1600 by 400 tails 2 %tables; *-----------------------------------------------------------------------; * Comparing nested (full vs. reduced) OLS regression models ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Corey Latour: OLS Modeling to Predict Job Satisfaction"; title3 "Adding a 4-level nominal predictor (3 dummy variables)."; title4; datalines; R**2 .45 .50 NumParms 5 8 Ntotal 100 125 150 alpha .05 .01 %tables; *-----------------------------------------------------------------------; * Testing a single beta in a k-variable OLS regression model ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Corey Latour: Predicting Salaries:"; title3 "Will the TeamValu scale signicantly improve the model?"; /* Predicting Y = annual salary (example: 32.4 = $32,400 US). X3 = TeamValu is a 0-4 visual analog (continuous) scale measuring the employee's value to his/her teammates as determined by an external evaluator. TeamValu's (unstandardized) beta coefficient, beta3, is the only one under examination here. It corresponds to the increase in annual salary (in $1000 US) associated with a 1-point increase in TeamValu given 2 other predictors, X1 and X2, whatever they are. The R**2 from predicting TeamValu from X1 and X2 is 0.30-0.40, e.g. TeamValuÕs tolerance is 0.60-0.70. */ datalines; 1betaOLS 1.5 2.0 . One or more conjectures for beta3 SDX .6 .8 . One or more conjectures for SD(TeamValu) Tolerance .6 .7 . One or more conjectures for Tol(TeamValu) SD 3 3.5 . One or more conjectures for SD(error) NumParms 4 . (= 3 for Xs + 1 for intercept) Ntotal 50 75 tails 1 %tables; ; /* I didn't like the default table for this problem, so I wrote a custom one. This also works nicely for using SAS/GRAPH. */ proc tabulate format=4.3 order=data; class alpha BetaWt tolernce SDx testtype SD Ntotal; var power; table alpha='Alpha:',BetaWt='Beta for TeamValu'* tolernce='Tol(X)' *SDx='SD(X)' * testtype='Type', SD='SD(Resid)' * Ntotal='Total N' * Power='Power'*mean=' ' /rtspace=35; run; *-----------------------------------------------------------------------; * General linear models (intro.): traditional t-test the hard way ; *-----------------------------------------------------------------------; /* A trick to handle any situation involving the normal-theory univariate general linear model. */ title2 "Dr. Alalgia's Two-Group Design (revisited)"; title3 "Simple Example to Introduce EXEMPLARY SSH problem,"; /* Set the groups, their conjectured means, and give sample sizes that reflect the cell weights. This creates an exemplary data set. For this problem, an exemplary data set is one having: sample means = "true" means, no error variance. The theory justifying this computational trick is described in O'Brien and Muller (1993, Section 8.3). The EXEMPLARY SSH problem type in UnifyPow replaces the PowSetUp module described in O'Brien and Muller. */ data; input feedback $ VHchange n_exemp ; datalines; ET -.86 20 SP -.42 10 proc print; var feedback VHchange n_exemp; /* "analyze" exemplary data to get SSHe values */ proc glm order=data; class feedback; freq n_exemp; model VHchange = feedback; /* Key part of output ------------------ Dependent Variable: VHCHANGE Frequency: N_EXEMP Source DF Type III SS Mean Square F Value Pr > F FEEDBACK 1 1.29066667 1.29066667 99999.99 0.0001 */ /* Use UnifyPow to compute and table powers or sample sizes */ %include "&UnifyPow"; datalines; /# Nobody should use UnifyPow this way to do a traditional t test! It is only done here to introduce the use of the EXEMPLARY SSH method to compute and table power or sample sizes for complex general linear models. #/ exemplary SSH NumParms 2 Nexemplary 30 sd .45 .57 .65 alpha .05 .01 NTotal 15 21 33 effects "ET vs SP on VHchange" 1 1.29066667 %tables; %include "&UnifyPow"; datalines; exemplary SSH NumParms 2 Nexemplary 30 sd .45 .57 .65 alpha .05 .01 power .80 .90 effects "ET vs SP on VHchange" 1 1.29066667 %tables; *-----------------------------------------------------------------------; * General linear models: one-way analysis of covariance, with contrasts ; *-----------------------------------------------------------------------; title2 "Dr. Mindy Bowdy: 4-Group Design + Covariate"; title3 "Example 4 in O'Brien and Muller (1993, Section 8.3)"; title4 "(D) Dominators, (R) Regulars {Ordinary or Loner}, (F) Friendlies"; /* Three groups, one covariate (LESI), unequal slopes. LESI is the Life Events Stress Index: -2 = low stress, -1 = below average, 0 = average, 1 = above average, 2 = high stress. Order of exemplary data: D = Dominator R = Regulars = {Ordinary or Loner} F = Friendly */ data; input lysis0 DRFgrp $ LESI n_exemp; /* create exemplary data */ if DRFgrp = 'D' then beta = -.03; if DRFgrp = 'R' then beta = -.01; if DRFgrp = 'F' then beta = .00; lysis = lysis0 + beta*LESI; * lysis at ; * LESI=0 DRF LESI n_exemp ; datalines; .3350 D -2 02 .3350 D -1 03 .3350 D 0 04 .3350 D 1 05 .3350 D 2 06 .5033 R -2 12 .5033 R -1 12 .5033 R 0 12 .5033 R 1 12 .5033 R 2 12 .6000 F -2 04 .6000 F -1 04 .6000 F 0 04 .6000 F 1 04 .6000 F 2 04 proc print; var lysis DRFgrp LESI n_exemp; /* "analyze" exemplary data to get SSHe values */ proc glm order=data; class DRFgrp; freq n_exemp; model lysis = DRFgrp DRFgrp*LESI/noint solution; contrast 'DvRvF main | LESI=0' DRFgrp 1 -1 0, DRFgrp 0 1 -1; contrast 'Ave LESI slope' DRFgrp*LESI .333 .333 .333; contrast 'DvRvF x LESI' DRFgrp*LESI 1 -1 0, DRFgrp*LESI 0 1 -1; contrast 'DvR | LESI=0' DRFgrp 1 -1 0; contrast 'RvF | LESI=0' DRFgrp 0 1 -1; contrast 'Slopes: DvR' DRFgrp*LESI 1 -1 0; contrast 'Slopes: RvF' DRFgrp*LESI 0 1 -1; /* Key parts of the output: -------------------------------------------------------------------------- Number of observations in data set = 100 Contrast DF Contrast SS Mean Square F Value Pr > F DvRvF main | LESI=0 2 0.6722149 0.3361074 99999.99 0.0001 Ave LESI slope 1 0.0258462 0.0258462 99999.99 0.0001 DvRvF x LESI 2 0.0175385 0.0087692 99999.99 0.0001 DvR | LESI=0 1 0.3837566 0.3837566 99999.99 0.0001 RvF | LESI=0 1 0.1402634 0.1402634 99999.99 0.0001 Slopes: DvR 1 0.0108387 0.0108387 99999.99 0.0001 Slopes: RvF 1 0.0030000 0.0030000 99999.99 0.0001 -------------------------------------------------------------------------- This information is then used in UnifyPow, as follows. */ %include "&UnifyPow"; datalines; Exemplary SSH Nexemplary 100 alpha .05 SD .12 .15 Ntotal 200 300 500 NumParms 6 tails 2 effects "DvRvF main | LESI=0" 2 0.6722149 "Ave LESI slope" 1 0.0258462 "DvRvF x LESI" 2 0.0175385 ; /* store these .05 results; */ data PowData1; set PowData; /* now run .025 tests (Bonferroni adjustment) */ %include "&UnifyPow"; datalines; Exemplary SSH Nexemplary 100 alpha .025 SD .12 .15 Ntotal 200 300 500 NumParms 6 tails 2 effects "DvR | LESI=0" 1 0.3837566 "RvF | LESI=0" 1 0.1402634 "Slopes: DvR" 1 0.0108387 "Slopes: RvF" 1 0.0030000 ; /* concantenate datasets and produce one table */ data PowData; set PowData1 PowData; %tables; *-----------------------------------------------------------------------; * Logit analysis (intro.): comparing 2 indep. proportions the hard way ; *-----------------------------------------------------------------------; /* A trick to handle any situation involving likelihood ratio tests in categorical modeling. */ title2 "Lactic Acidosis in Children with Malaria (revisited)"; title3 "28% die untreated. What if DCA cuts this by 25%?"; /* Input the groups (here, Trtment), their conjectured mortality rates (PrDeath), and give sample sizes that reflect the cell weights (n_exemp). Then compute exemplary frequencies for each outcome category (ExempFrq). This creates an exemplary data set (i.e., sample proportions = "true" proportions). The theory justifying this approach is given in O'Brien (1986, SUGI-11 Proceedings, 778-784) and summarized in Agresti (1990, Analysis of Categorical Data, Wiley, p 243). */ data; input Trtment $ PrDeath n_exemp; if Trtment = "Placebo" then DCA = 0; else DCA = 1; died = 1; ExempFrq = PrDeath*n_exemp; output; died = 0; ExempFrq = (1-PrDeath)*n_exemp; output; /* Trtment PrDeath n_exemp */ datalines; Placebo .28 100 DCA .21 200 proc print; var Trtment DCA died PrDeath n_exemp ExempFrq; /* "analyze" exemplary data to get -2lnL values */ proc logistic; weight ExempFrq; model died = DCA; /* Key part of output ------------------ Criterion Only Covariates Chi-Square for Covariates -2 LOG L 325.964 324.173 1.790 with 1 DF (p=0.1809) */ /* Note: The value of 1.790 will be used in the UnifyPow EFFECTS specification. */ /* Use UnifyPow to compute and table powers or sample sizes */ %include "&UnifyPow"; datalines; /# Nobody should use UnifyPow this way to test two independent proportions. It is only done here to introduce the use of the EXEMPLARY CHI**2 method to compute and table powers or sample sizes for complex logistic regression or log-linear models, or similar methods using -2lnL chi**2 test statistics. #/ exemplary chi**2 Nexemplary 300 alpha .01 .045 NTotal 750 999 1500 effects "DCA vs. Placebo" 1 1.790 %tables; *-----------------------------------------------------------------------; * Logit analysis ; *-----------------------------------------------------------------------; data babies; title2 "Get Exemplary D values for Study of Lyon's Nonstress Test"; title3 "See O'Brien (1986, SUGI-11 Proceedings, 778-784)"; Title4 "Summarized in Agresti (1990, Analysis of Categorical Data, Wiley, p 243)."; /* "PrNRcare" stands for Probability that "Non-Routine" cardiac care is given at birth. The question is whether a new diagnostic test ("Lyons") provides additional value in predicting NRcare beyond that povided by using only the Standard test. Outcome at birth "NonRoutn" 0 = no non-routine care was given 1 = non-routine care was given Predictors before birth "Standard" test 1 = worrisome result 2 = reassuring result "Lyons" test 1 = very worrisome result 2 = somewhat worrisome result 3 = reassuring result */ input Standard Lyons CellProb PrNRcare; * expected number of routine outcomes; NonRoutn=0; ExempFrq = (1-PrNRcare)*CellProb*1000; output; * expected number of nonroutine outcomes; NonRoutn=1; ExempFrq = PrNRcare*CellProb*1000; output; *Standard Lyons CellProb PrNRcare ; datalines; 1 1 .04 .40 1 2 .08 .32 1 3 .04 .27 2 1 .02 .30 2 2 .18 .22 2 3 .64 .15 data babies; set babies; *build dummy vars for PROC LOGIST; XStd = (Standard=2) - (Standard=1); XLyons1 = (Lyons=1) - (Lyons=3); XLyons2 = (Lyons=2) - (Lyons=3); XSL1 = XStd*XLyons1; XSL2 = XStd*XLyons2; proc print; var Standard Lyons CellProb PrNRCare NonRoutn ExempFrq; /* "analyze" exemplary data to get -2lnL values */ proc logistic; weight ExempFrq; model NonRoutn = Xstd; proc logistic; weight ExempFrq; model NonRoutn = Xstd XLyons1; proc logistic; weight ExempFrq; model NonRoutn = Xstd XLyons1 XLyons2; proc logistic; weight ExempFrq; model NonRoutn = Xstd XSL1 XSL2; proc logistic; weight ExempFrq; model NonRoutn = Xstd XLyons1 XLyons2 XSL1 XSL2; /* These analyses provide the -2LnL values for UnifyPow's EFFECTS specifications. */ %include "&UnifyPow"; datalines; Exemplary chi**2 Nexemplary 1000 alpha .01 .05 Ntotal 400 600 1000 effects "Standard Only [LR]" 1 983.943 964.437 "Lyons(lin) Given Standand [LR]" 1 964.437 956.319 "Lyons Given Standard [LR]" 2 964.437 956.277 "Lyons Main [LR, Sat Mdl]" 2 960.910 955.990 ; /* save results based on likelihood ratio tests */ data PowData2; set PowData; /* The following UnifyPow problem handles the Nonstress Test study using Wald-type contrasts on the 6 logits. This is regarded as a cruder (yet one-step) substitute for the two-step EXEMPLARY CHI**2 method just given. */ %include "&UnifyPow"; datalines; pi .40 .32 .27 .30 .22 .15 weight .04 .08 .04 .02 .18 .64 alpha .01 .05 Ntotal 400 600 1000 NoOverall contrasts "Lyons Main [Wald, Sat Mdl]" 1 -1 0 1 -1 0 > 0 1 -1 0 1 -1 "Lyons(lin) SubMain [Wald, Sat Mdl]" 1 0 -1 1 0 -1 ; /* Concatenate LRT and Wald-type results and make one table */ data PowData; set PowData2 PowData; ProbStmt = "pi = {.40 .32 .27 .30 .22 .15}, CellProbs = {.04 .08 .04 .02 .18 .64}"; %tables; title1 "SUGI 23 (Nashville, 1998) Presentation"; *-----------------------------------------------------------------------; * Comparing two independent proportions ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "BeeBop's XDM-X vs. XDM running shoe: Serious injury (Y/N)"; title3 "Comparing Two Independent Proportions"; datalines; pi .06 .03 weight 1 2 NTotal 201 270 alpha .05 .01 %tables; *-----------------------------------------------------------------------; * Testing two means (locations) ; *-----------------------------------------------------------------------; options pagesize = 100; %include "&UnifyPow"; title2 "XDM-X vs. XDM running shoe: arcsin(sqrt(%)) days injured"; title3 "Testing Two Means (Locations)"; datalines; /# median(%): 9% vs. 7% #/ mu .30 .27 weight 1 2 SD .08 .10 .125 NTotal 201 270 Wilcoxon methods all %tables; *-----------------------------------------------------------------------; * One-way ANOVA with complex contrasts ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "XDM-S vs XDM-X1, XDM-X2, XDM-X3 subtypes"; title3 "Testing on 4 means via ANOVA with complex contrasts"; datalines; mu .300 .275 .270 .265 /# 2/3 of the runners will get some version of the experimental shoe. #/ weight 3 2 2 2 SD .08 .10 .125 NTotal 207 270 NoOverall contrasts "XDM-S vs. XDM-X (all versions)" 3 -1 -1 -1 "Variation among XDM-X subtypes" 0 1 -1 0 > 0 0 1 -1 %tables; *-----------------------------------------------------------------------; * McNemar's test of 2 correlated proportions ; *-----------------------------------------------------------------------; %include "&UnifyPow"; title2 "Thromboembolism (case/control) & The Pill ('outcome')"; title3 "Matched Pairs, Binary Outcome"; /* Example follows study by Sartwell PE, et. al. (1969, Thromoembolism and oral contraceptives: an epidemiologic case-control study. Am J Epi, 90: 365-380.) Cases: Women of childbearing age who have some kind of thromoembolism. Matched Controls: One per case--treated in the same hosptal, matched on age, number of prior pregnancies, income level. Scenario: Matched Control Using The Pill? No Yes | ------------------------ Thromboembolism Case No | p00 = .54 p01 = .08 Using The Pill? | Yes | p10 = .32 p11 = .06 */ datalines; McNemar .32 .08 alpha .05 .01 power .90 .95 %tables; *-----------------------------------------------------------------------; * Examples for: ; * "Power for Categorical Data Analysis: A Simpler Method to Compute ; * Power for Likelihood Ratio Tests in Generalized Linear Models" ; * Gwowen Shieh and Ralph O'Brien ; * Invited paper for 1998 Joint Statistical Meetings (Dallas) ; *-----------------------------------------------------------------------; *-----------------------------------------------------------------------; * Logistic regression. ; *-----------------------------------------------------------------------; /* Suppose that MA710808 is a possible monoclonal antibody therapy for patients in critical condition with severe sepsis. The proposed design is a randomized, placebo controlled, blinded Phase III trial in which 2/3 will receive MA710808 and 1/3 will get placebo. Several baseline measures are collected and combined to create a baseline severity index, SevIndx0. SevIndx0 is a continuous measure developed in the previous Phase II trial and in this population ranges from 1.5 (least severe) to 9.5 with a tri-modal shape. For the sake of pragmatism, we shall simplify somewhat and pretend that its distribution is discrete: SevIndx0: 2 3 4 5 6 7 8 9 Prob: .10 .10 .05 .25 .25 .05 .10 .10 Code the dummy variable MA710808 to be 0 if placebo and 1 of MA710808. Thus, Pr[SevIndx0 = 2, MA710808 = 0] = 1/10 * 1/3 = 1/30 = 2/60. From the Phase II results we assume that the logit is logit = 2.5 - 0.5*SevIndx0 + ln(3)*MA710808. This sets an odds ratio of 3 in favor of MA710808, a large effect, but actually less than the odds ratio of almost 5.0 observed in the Phase II study. */ data sepsis; title1 "Power for Categorical Data Analysis, Shieh-O'Brien, JSM 1998"; title2 "Get Exemplary LR test values for trial of MA710808 for severe sepsis"; keep SevIndx0 ProbSevI MA710808 probX PrAliv25 Alive25 ExempCnt; input SevIndx0 ProbSevI; TotNExmp = 10002; *Total N of exemplary data; RanRateP = 1/3; *Randomization Rate for the Placebo group; betaInt = 2.5; betaSev = -0.5; betaMA = log(3); *Note: ln(3) = 1.0986; /* Each dataline gives one SevIndx0 group. From each we generate "exemplary" counts (ExempCnt) for 4 types of cases based on the total exemplary sample size, TotNExmp: (1) Cases on placebo who are not alive at Day 25. (2) Cases on placebo who are alive at Day 25. (3) Cases on MA710808 who are not alive at Day 25. (4) Cases on MA710808 who are alive at Day 25. */ do MA710808 = 0 to 1; * 1 = randomized to MA710808; logit = betaInt + betaSev*SevIndx0 + betaMA*MA710808; PrAliv25 = 1/(1 + exp(-logit)); *Prob that case is alive at Day 25; if MA710808 = 0 then RandRate = RanRateP; else RandRate = 1 - RanRateP; probX = RandRate*ProbSevI; do Alive25 = 0 to 1; * 1 = alive at Day 25; if Alive25 = 0 then probY = (1-PrAliv25); else probY = PrAliv25; ExempCnt = round(TotNExmp*probX*probY); *Exemplary count; output; end; end; *SevIndx0 ProbSevI ; datalines; 2 .10 3 .10 4 .05 5 .25 6 .25 7 .05 8 .10 9 .10 proc print; var SevIndx0 ProbSevI MA710808 probX PrAliv25 Alive25 ExempCnt; /* "analyze" exemplary data to get LR chi**2 values */ proc genmod; model Alive25 = SevIndx0 MA710808 / dist = binomial link = logit type3; freq ExempCnt; run; title2 "Sample-size analysis for MA710808 Phase III study"; %include "&UnifyPow"; datalines; Exemplary chi**2 Nexemplary 10005 . Not exactly 10000 due to rounding of exemplary events. alpha .01 .05 power .80 .90 .95 effects "MA710808 vs. Placebo" 1 528.9458 %tables; %include "&UnifyPow"; title3 "Get powers for NTotals that fit 1:2 randomization"; datalines; Exemplary chi**2 Nexemplary 10005 . Not exactly 10000 due to rounding of exemplary events. alpha .01 NTotal 222 282 339 effects "MA710808 vs. Placebo" 1 528.9458 %tables; %include "&UnifyPow"; title3 "Get powers for NTotals that fit 1:2 randomization"; datalines; Exemplary chi**2 Nexemplary 10005 . Not exactly 10000 due to rounding of exemplary events. alpha .05 NTotal 150 201 249 effects "MA710808 vs. Placebo" 1 528.9458 %tables; *-----------------------------------------------------------------------; * Poisson regression. ; *-----------------------------------------------------------------------; /* The Reliable Web Server Company is relying on 8 large identical computers to support all of its clients. Since Reliable upgraded its operating system 12 months ago to JupiterOS v3.0, it has been experiencing a number of unexplained system crashes, averaging over 1 crash per week per computer, but with some variation over computers: Computer: A B C D E F G H Crashes/week: 1.2 2.0 0.6 0.8 1.4 1.2 1.0 0.4 This is barely acceptable. Now there is has a new release, JupiterOS v3.1, which is supposed to greatly reduce the number of unexplained crashes. Reliable's systems managers are skeptical of such claims and worried that the costs involved in upgrading to 3.1 may not justify the increased reliability, or, worse, that the new version may be even more crash prone. They decide to upgrade 2 of the 8 machines so they can compare v3.0 and v3.1 over a number of weeks. Any other changes to the computers will be done identically to all 8. How long should they test v3.1? We take the number of crashes per week on Computer i running on v3.x to be distributed as a Poisson random variable with a mean of m(i,v3.x). For the sammple-size analysis we shall presume that v3.1 has 1/3 the number of crashes as v3.0, i.e. m(i, v3.1) = m(i, v3.0)/3. We shall also assume that the past crash rates (PastCrRt) given above are predictive of the future rates. Thus, the Poisson linear model for m(i, v3.x) takes the form ln{m(i, v3.x)} = betaInt + betaPCR*ln[PastCrRt(i)] + betaV31*V31(i), where V31(i) = 1 if the computer is running v3.1, and V31(i) = 0 otherwise. The true model we presume has betaInt = 0.0 No discrete jump in crash rates. betaPCR = 1.0 Past rates carry over to future rates. betaV31 = -ln(3) v3.1 crash rates are 1/3 those of v3.0. This true model can also be described as m(i, v3.0) = PastCrRt(i) m(i, v3.1) = PastCrRt(i)/3. Computers B and D are selected for v3.1 upgrades. */ data crashes; title2 "Get exemplary LR test value for computer crashing problem"; keep Computer PastCrRt V31 m Ncrashes lnPCR lnTotWks; input Computer $ PastCrRt V31; TotWeeks = 100; betaInt = 0.0; betaPCR = 1.0; betaV31 = -log(3); *Note: -ln(3) = = -1.0986; m = exp(betaInt + betaPCR*log(PastCrRt) + betaV31*V31); Ncrashes = TotWeeks*m; lnPCR = log(PastCrRt); lnTotWks = log(TotWeeks); * Computer PastCrRt V31; datalines; A 1.2 0 B 2.0 1 C 0.6 0 D 0.8 1 E 1.4 0 F 1.2 0 G 1.0 0 H 0.4 0 proc print; var Computer PastCrRt V31 m Ncrashes; /* "analyze" exemplary data to get LR chi**2 values */ proc genmod; model Ncrashes = lnPCR V31 / dist = Poisson link = log offset = lnTotWks type3; run; title2 "Sample-size analysis for JupiterOS v3.1 evaluation"; title3 "Units = weeks (all eight machines)"; %include "&UnifyPow"; datalines; Exemplary chi**2 Nexemplary 100 . weeks (all eight machines) alpha .05 .10 power .90 .95 tails 2 effects "v3.1 vs. v3.0" 1 104.5144 %tables; title3 "Units = weeks (running all eight machines)"; %include "&UnifyPow"; datalines; Exemplary chi**2 Nexemplary 100 . weeks (running all eight machines) alpha .05 .10 NTotal 9 11 13 tails 2 effects "v3.1 vs. v3.0" 1 104.5144 0 %tables;