Documentation Example 11 for PROC POWER
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: pwrex11 */
/* TITLE: Documentation Example 11 for PROC POWER */
/* (Logistic Regression with the CUSTOM Statement) */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: power */
/* sample size */
/* power analysis */
/* noncentrality */
/* logistic regression */
/* PROCS: POWER */
/* DATA: */
/* */
/* SUPPORT: John Castelloe */
/* REF: */
/* MISC: */
/****************************************************************/
* Compute P(Yj|Xi) values from conjectured betas;
data Exemplary;
_B0 = 0; _B1 = -1.5; _B2 = 0.5; _B3 = 2;
do X1 = 0 to 1;
do X2 = 0 to 1;
_pi = logistic(_B0 + _B1*X1 + _B2*X2 + _B3*X1*X2);
do Y = 0 to 1;
PY = _pi**Y * (1-_pi)**(1-Y);
output;
end;
end;
end;
keep X1 X2 Y PY;
run;
proc print data=Exemplary;
run;
* Expand according to relative allocations of design profiles;
data Exemplary;
set Exemplary;
if (X1 = 0 and X2 = 0) then _alloc=4;
else if (X1 = 0 and X2 = 1) then _alloc=5;
else if (X1 = 1 and X2 = 0) then _alloc=5;
else _alloc=6;
do _i = 1 to _alloc;
output;
end;
keep X1 X2 Y PY;
run;
proc print data=Exemplary;
run;
* Run full model;
proc logistic data=Exemplary;
ods output parameterestimates=fitFull fitstatistics=statsFull;
weight PY;
model Y = X1 X2 X1*X2;
run;
* Fetch Wald chi-square statistic for test of X1*X2;
data fitFull;
set fitFull;
where Variable = "X1*X2";
keep DF WaldChiSq;
run;
* Fetch -2LogL for full model;
data statsFull;
set statsFull;
where Criterion = "-2 Log L";
Neg2LogLFull = InterceptAndCovariates;
keep Neg2LogLFull;
run;
* Run reduced model;
proc logistic data=Exemplary;
ods output fitstatistics=statsRed;
weight PY;
model Y = X1 X2;
run;
* Fetch -2LogL for reduced model;
data statsRed;
set statsRed;
where Criterion = "-2 Log L";
Neg2LogLRed = InterceptAndCovariates;
keep Neg2LogLRed;
run;
* Compute primary noncentralities for Wald and LR tests;
data PrimNC;
merge fitFull statsFull statsRed;
WaldPrimNC = WaldChiSq / 20;
LRPrimNC = -(Neg2LogLFull-Neg2LogLRed) / 20;
keep DF WaldPrimNC LRPrimNC;
call symput ("DF", DF);
call symput ("WaldPrimNC", WaldPrimNC);
call symput ("LRPrimNC", LRPrimNC);
run;
proc print data=PrimNC;
run;
proc power;
custom
dist = chisquare
primnc = &WaldPrimNC &LRPrimNC
testdf = &DF
ntotal = 100
power = .;
run;
data Power;
set PrimNC;
Crit = quantile("chisq", 1-0.05, DF);
WaldPower = sdf("chisq", Crit, DF, 100 * WaldPrimNC);
LRPower = sdf("chisq", Crit, DF, 100 * LRPrimNC);
keep WaldPower LRPower;
run;
proc print data=Power;
run;