/* ++++ Begin UPw01_Intro segment v. 2002.08.17a ++++ */ /* ---------------------------------------------------------------------- FOR INFORMATION AND LATEST VERSION: http://www.bio.ri.ccf.org/UnifyPow ---------------------------------------------------------------------- UnifyPow.sas, module (non-macro) version, 2002.08.17a 2002 Copyright (c) by Ralph G. O'Brien, PhD Department of Biostatistics and Epidemiology Cleveland Clinic Foundation Cleveland, OH 44195 Voice: 216-445-9451 robrien@bio.ri.ccf.org This single file contains all components of UnifyPow. You simply put it in an appropriate directory on your system and %include it in your SAS run. It is distributed so that it runs in the way I have been demonstrating in my workshops. ------------------------------------------- RUNNING UnifyPow AS DISTRIBUTED (NON-MACRO) ------------------------------------------- Regardless of the platform, to run this module version of UnifyPow and get the standard tables, simply follow this template: ---------------------------------------------------------------------- options ls=78 nosource2; %let UnifyPow = your_file_specification; %include "&UnifyPow"; title1 "A Title for Problem 1."; datalines4; {set of UnifyPow statements here} ;;;; %tables %include "&UnifyPow"; title1 "A Title for Problem 2."; datalines4; {second set of UnifyPow statements here} ;;;; %tables ---------------------------------------------------------------------- -------------------------------- THE FILE SPECIFICATION STATEMENT -------------------------------- For UNIX, use something like: ---------------------------------------------------------------------- %let UnifyPow = /home/robrien/SASmacros/UnifyPow.sas; ---------------------------------------------------------------------- With Windows, you could put UnifyPow in the C:\Sas\ folder and use: ---------------------------------------------------------------------- %let UnifyPow = C:\Sas\UnifyPow.sas; ---------------------------------------------------------------------- --------------------------------------------------- EXPERIENCED SAS USERS MAY CUSTOMIZE UnifyPow OUTPUT --------------------------------------------------- All results from each UnifyPow problem are stored in a temporary SAS data set called PowData. Knowing this, experienced SAS users may easily customize their output by merging results from two or more problems and by using their own PROC TABULATE or SAS/GRAPH code. The examples.sas file I distribute should have one or two examples of this. I recommend that you first examine the structure of PowData by just seeing what it holds. Non-macro version: ---------------------------------------------------------------------- options ls=78 nosource2; %let UnifyPow = file_specification; %include "&UnifyPow"; datalines4; {set of UnifyPow statements here} ;;;; proc print data=PowData; run; ---------------------------------------------------------------------- ---------------------- UnifyPow LEGAL NOTICES ---------------------- THIS SOFTWARE IS MADE AVAILABLE "AS IS". UnifyPow is a trademark of Ralph G. O'Brien. No commercial use of this trademark may be made without prior written permission of Ralph O'Brien. All UnifyPow software and its included text and accompanying documentation are Copyright 1997 by Ralph G. O'Brien. Permission to use, copy, modify, and distribute this software and its documentation for any purpose and without fee to Ralph O'Brien is hereby granted, provided that these legal notices appear in all copies and supporting documentation, that the name "UnifyPow" is retained, and that the names of Ralph O'Brien and the Cleveland Clinic Foundation are not used in advertising or publicity pertaining to distribution of the software without the specific, written prior permission of Ralph O'Brien. Although the above trademark and copyright restrictions do not convey the right to redistribute derivative works, Ralph O'Brien encourages unrestricted distribution of patches or ancillary code which can be applied to or used in conjunction with Ralph O'Brien's distribution. If this software is modified for local use, please denote this on all modified versions of the software by appending the letter "L" to the current version number and by noting the changes in the code itself and in the associated documentation. RALPH G. O'BRIEN AND THE CLEVELAND CLINIC FOUNDATION DISCLAIM ALL WARRANTIES, EXPRESS OR IMPLIED, WITH REGARD TO THIS SOFTWARE, INCLUDING WITHOUT LIMITATION ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, AND IN NO EVENT SHALL RALPH G. O'BRIEN OR THE CLEVELAND CLINIC FOUNDATION BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA, OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, TORT (INCLUDING NEGLIGENCE) OR STRICT LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. ----------------------------------------- REPORTING PROBLEMS AND MAKING SUGGESTIONS ----------------------------------------- UnifyPow's computations are checked thoroughly using multiple methods: comparing its results to those obtained by using other software, to those in published tables and examples, and to results obtained from Monte Carlo simulation. No true bugs have yet been reported to me about UnifyPow's numerical accuracy, but no software can be said to be totally free of problems. I really do appreciate knowing if you encounter difficulties or have suggestions for improvements. Most of the additions to UnifyPow come about this way. On the other hand, I cannot promise to address everyone's queries, especially those that are mostly related to consulting advice. I correspond best via email. Do not modify this code without obtaining written permission from me. If changes are needed for some purpose, have the wisdom and courtesy to consult me. --------------- CITING THE WORK --------------- If you use this work, please cite it--something like ... using UnifyPow, a module for the SAS System (O'Brien, 1998). This reference is O'Brien, RG (1998), "A Tour of UnifyPow: A SAS Macro for Sample-Size Analysis," Proceedings of the 23rd SAS Users Group International Conference, Cary NC: SAS Institute, 1346-1355. Another key reference for this work is: O'Brien, RG and Muller, KE (1993), "Unified Power Analysis for t Tests through Multivariate Hypotheses," in Edwards, EK. (Ed.), Applied Analysis of Variance in the Behavioral Sciences, New York: Marcel Dekker, pp. 297-344. Actually, that chapter covers freeware called OneWyPow.sas and other %include modules that are now unified in UnifyPow.sas, which offers a simpler interface and a far larger set of methods. --------------------------------------- OBTAINING FREEWARE SAS CODE AND UPDATES --------------------------------------- UnifyPow is freeware distributed via http://www.bio.ri.ccf.org/UnifyPow a website at the Department of Biostatistics and Epidemiology at the Cleveland Clinic Foundation. Register! Downloading new versions periodically ensures that you are getting the latest version I feel is safe for public distribution. At the website, you may register your name and email address, so that I can notify you of new releases. */ /* ++++ End UPw01_Intro segment ++++ */ /* ++++ Begin UPw02_Tables segment v. 2002.07.01 ++++ */ /* ---------------------------- SOURCE CODE for tables macro ---------------------------- */ %macro tables; /* All kinds of tables for displaying results. Easy to customize. */ %if &DoTables = no %then %do; %let DoTables = yes; %goto donetabl; %end; /* The tables are constructed so that they will handle either "Power" or "NTotal" (TotalN or TotalPairs or PairsTotal) statements. */ %let CmpSDFt1 = Note: values here are composites formed by taking a weighted; %let CmpSDFt2 = average over the results obtained across the SD specifications.; %if &ResltVar = power %then %do; %let SpecVar = NTotal; %let result = %str(Power="Power"*mean=" "); %let spec = %str(Ntotal="Total N"); %if &ProbType = McNemar %then %let spec = %str(Ntotal="Pairs"); %if (&ProbType = PairedMu) or (&ProbType = PairedMu:LogNormal) %then %let spec = %str(Ntotal="Total Pairs"); %let format = 4.3; %end; %if &ResltVar = NTotal %then %do; %let SpecVar = NomPower; %let result = %str(Ntotal="Total N"*mean=" "); %if &ProbType = McNemar %then %let result = %str(Ntotal="Pairs"*mean=" "); %if &ProbType = PairedMu %then %let result = %str(Ntotal="Total Pairs"*mean=" "); %let spec = %str(NomPower="Minimum Power"); %let format = 6.0; %end; %if %index(&TablType, WlcxnPow) > 0 %then %let AddRow2 = %str(* parent = "Parent"); %if &TablType = GnrlPow %then %do; /* all alphas in one table, no SD or parent */ proc tabulate format = &format order=data; class alpha EffcTitl testtype &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:", &MethLabl * testtype="Statistic", alpha * &spec * &result /rtspace=28; %if &fnote=2piChi2 or &fnote=2piBoth %then %do; %FtntDago %end; %if &fnote=2piXctUn or &fnote=2piBoth %then %do; %FtntSuis %end; run; footnote1; footnote2; footnote3; footnote4; footnote5; footnote6; footnote7; %end; /* GnrlPow */ %if &TablTyp3 = Pi1Specl %then %do; /* tables of TruAlpha and critical values for 1-group binomial problem */ proc tabulate format=6.0 order=data; class SpcTitl3 EffcTitl alpha testtype &SpecVar; format alpha &AlphaFmt; format NomPower &PowerFmt; var TruAlpha LoCrit HiCrit; table SpcTitl3 = " ", EffcTitl="Method" * &spec * testtype="Type", alpha*(TruAlpha = "Actual Alpha"*mean=" "*F=&AlphaFmt LoCrit="Lower Crit Value"*mean=" " HiCrit="Upper Crit Value"*mean=" ") /rtspace = 28; footnote1 "These critical values are part of the rejection region."; footnote2 "&Pi1Foot2"; footnote3 "&Pi1Foot3"; run; footnote1; footnote2; footnote3; %end; /* pi1Specl */ %if &TablType = Specl2pi %then %do; /* special 2 pi problem, including odds ratio and prob ratio spec */ proc tabulate format = &format order=data; class alpha EffcTitl testtype pi_1 &pi2SpcVr &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:" * EffcTitl = "Method:", alpha="Alpha" * testtype="Test Type", pi_1="Group 1 Probability" * &pi2SpcVr="&pi2SpcLb" * &spec * &result /rtspace=28; %if &fnote=2piChi2 or &fnote=2piBoth %then %do; %FtntDago %end; %if &fnote=2piXctUn or &fnote=2piBoth %then %do; %FtntSuis %end; run; footnote1; footnote2; footnote3; footnote4; footnote5; footnote6; footnote7; %end; %if &TablType = tPow %then %do; %if &TablSDwt = no %then %do; /* all alphas in one table, does not give parent */ proc tabulate format = &format order=data; class alpha EffcTitl testtype &SDorCV &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:"* EffcTitl="Effect:", alpha= "Alpha" * testtype="Type", &SDorCV="&SDType" * &spec * &result /rtspace=28; %end; %if &TablSDwt = yes %then %do; /* weight using SDwgt */ proc tabulate format = &format order=data; class alpha EffcTitl testtype &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; weight SDwgt; table scenario="Scenario:"* EffcTitl="Effect:", alpha= "Alpha" * testtype="Type", &spec * &result /rtspace=28; footnote1 "&CmpSDFt1"; footnote2 "&CmpSDFt2"; run; footnote1; footnote2; %end; %end; /* tpow */ %if &TablType = WlcxnPow %then %do; %if &TblMmnts = yes and &TablSDwt = no %then %do; data PowDatB; set PowData; if Wp1 ne .; proc tabulate data=PowDatB format = 4.3 order=data; class parent &SDorCV; var Wp1 Wp2 W&p3or4; table parent="Parent", mean="Nonparametric Moments"*(Wp1="p1" Wp2="p2" W&p3or4="&p3or4")* &SDorCV="&SDType" / rtspace=11 misstext='N/A'; %end; /* TblMmnts = yes (1) */ %if &TablSDwt = no %then %do; /* separate tables for each alpha, gives parent */ proc tabulate data=PowData format= &format order=data; class alpha EffcTitl testtype parent &SDorCV &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:"*alpha="Alpha:", EffcTitl="Method"* testtype="Type" *parent="Parent", &SDorCV="&SDType" * &spec * &result /rtspace=28; %end; %if &TablSDwt = yes %then %do; /* no SD, so no need for separate alpha tables */ proc tabulate data=PowData format= &format order=data; class alpha EffcTitl testtype parent &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; weight SDwgt; table scenario="Scenario:", EffcTitl="Method"* testtype="Type" *parent="Parent", alpha="Alpha"* &spec * &result /rtspace=28; footnote1 "&CmpFoot1"; footnote2 "&CmpFoot2"; run; footnote1; footnote2; %end; %end; /* TablType = WlcxnPow */ %if &TablType = 1or2WlcxnPow %then %do; %if &TblMmnts = yes %then %do; data PowDatB; set PowData; if Wp1 ne .; %if &SDtype ne none %then %do; proc tabulate data=PowDatB format=4.3 order=data; class parent &SDorCV; var Wp1 Wp2 W&p3or4; table parent="Parent", mean="Nonparametric Moments"*(Wp1="p1" Wp2="p2" W&p3or4="&p3or4")* &SDorCV="&SDType" / rtspace = 11; %end; /* SDtype ne none (1) */ %if &SDtype = none %then %do; proc tabulate data=PowDatB format=4.3 order=data; class parent; var Wp1 Wp2 W&p3or4; table parent = "Parent", mean="Nonparametric Moments" * (Wp1="p1" Wp2="p2" W&p3or4="&p3or4") /rtspace=11; %end; /* SDType = none (1) */ %end; /* TblMmnts = yes (2) */ %if &SDtype ne none %then %do; %if &TablSDwt = no %then %do; /* separate tables for each alpha, gives parent */ proc tabulate data=PowData format= &format order=data; class alpha EffcTitl testtype parent &SDorCV &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:"*alpha="Alpha:", EffcTitl="Method"* testtype="Type" *parent="Parent", &SDorCV="&SDtype" * &spec * &result /rtspace=28; %end; %if &TablSDwt = yes %then %do; /* table with composite over SD */ proc tabulate data=PowData format= &format order=data; class alpha EffcTitl testtype parent &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:", EffcTitl="Method"* testtype="Type" *parent="Parent", alpha="Alpha" * &spec * &result /rtspace=28; footnote1 "&CmpFoot1"; footnote2 "&CmpFoot2"; run; footnote1; footnote2; %end; %end; /* SDtype ne none (2) */ %if &SDtype = none %then %do; proc tabulate data=PowData format=&format order=data; class alpha EffcTitl &DistList testtype parent &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:" &OthScen, EffcTitl=" Method"* testtype="Type" *parent="Parent", alpha="Alpha" * &spec * &result /rtspace=28; %end; /* SDType = none (2) */ %end; /* TablType = 1or2WlcxnPow */ %if &TablType = WlcxMuPpow %then %do; %if &TblMmnts = yes %then %do; data PowDatB; set PowData; if Wp1 ne .; proc tabulate data=PowDatB format=6.3 order=data; class parent SDP Trials; var Wp1 Wp2 W&p3or4; table parent="Parent", mean="Nonparametric Moments"*(Wp1="p1" Wp2="p2" W&p3or4="&p3or4")* SDP="SD(P)" * Trials="&TrlName" / rtspace = 11; %end; /* TblMmnts = yes (3) */ /* separate tables for each alpha, gives parent */ proc tabulate data=PowData format=&format order=data; class alpha EffcTitl testtype parent SDP Trials &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:"*alpha="Alpha:", EffcTitl="Method"* testtype="Type" *parent="Parent", SDP="SD(P)" * Trials="&TrlName" * &spec * &result /rtspace=28; %end; /* TablType = WlcxMuPpow */ %if &TablType = WlcxPaMuPow %then %do; %if &TblMmnts = yes %then %do; data PowDatB; set PowData; if Wp1 ne .; proc tabulate data=PowDatB format=4.3 order=data; class parent SDMult Corr; var Wp1 Wp2 W&p3or4; table parent="Parent", SDMult="&MultLabl" * Corr="&CorrLabl" * mean="Nonparametric Moments"*(Wp1="p1" Wp2="p2" W&p3or4="&p3or4") / rtspace = 11; %end; /* TblMmnts = yes (4) */ proc tabulate data=PowData format=&format order=data; class alpha EffcTitl testtype parent SDMult Corr &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:" * alpha= "Alpha", EffcTitl="Method" * testtype="Type" * parent="Parent", SDMult="&MultLabl" * Corr="&CorrLabl" * &spec * &result /rtspace=28; %end; /* TablType = WlcxPaMuPow */ %if &TablType = PaMuPow %then %do; proc tabulate format=&format order=data; class alpha EffcTitl testtype SDMult Corr &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:"* EffcTitl="Effect:", alpha= "Alpha" * testtype="Type", SDMult="&MultLabl" * Corr="&CorrLabl" * &spec * &result /rtspace=28; %end; /* TablType = PaMuPow */ %if &TablType = muANCOVA %then %do; proc tabulate format=&format order=data; class alpha EffcTitl testtype &SDorCV &SDadjVar &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:", EffcTitl="Effect" * alpha= "Alpha" * testtype="Type", &SDorCV="&SDType" * &SDadjVar="&SDadLabl" * &spec * &result /rtspace=28; %end; /* TablType = muANCOVA */ %if &TablType = FPaMuPow %then %do; proc tabulate format=&format order=data; class alpha EffcTitl testtype SDMult Corr &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:" * EffcTitl="Effect:", alpha= "Alpha" * testtype="Type", SDMult="&MultLabl" * Corr="&CorrLabl" * &spec * &result /rtspace=28; %end; /* TablType = FPaMuPow */ %if &TablType = tmuP_Pow %then %do; proc tabulate format=&format order=data; class alpha EffcTitl testtype SDP Trials &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:"* EffcTitl="Effect:", alpha= "Alpha" * testtype="Test Type", SDP="SD(P)" * Trials="&TrlName" * &spec * &result /rtspace=28; %end; /* TablType = tmuP_Pow */ %if &TablType = FMuP_Pow %then %do; proc tabulate format=&format order=data; class alpha EffcTitl testtype SDP Trials &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:", EffcTitl=" " * alpha= "Alpha" * testtype="Test Type", SDP="SD(P)" * Trials="&TrlName" * &spec * &result /rtspace=28; %end; /* TablType = FMuP_Pow */ %if &TablType = FPow %then %do; %if &TablSDwt = no %then %do; /* condensed tables for F tests */ proc tabulate format = &format order=data; class alpha EffcTitl testtype &SDorCV &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table scenario="Scenario:", EffcTitl="Test" * alpha="Alpha" * testtype="Type", &SDorCV="&SDType" * &spec * &result /rtspace=28; %end; %if &TablSDwt = yes %then %do; /* condensed tables for F tests, composite over SD */ proc tabulate format = &format order=data; class alpha EffcTitl testtype &SpecVar scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; weight SDwgt; table scenario="Scenario:", EffcTitl="Test" * testtype="Type", &spec * &result /rtspace=28; footnote1 "&CmpSDFt1"; footnote2 "&CmpSDFt2"; run; footnote1; footnote2; %end; %end; /* TablType = FPow */ %if &TablType = 1betaOLS %then %do; proc tabulate format=&format order=data; class alpha BetaWt tolernce SDx testtype SD &SpecVar; format alpha &AlphaFmt; format NomPower &PowerFmt; var &ResltVar; table alpha="Alpha:"*BetaWt="Beta Coefficient:", tolernce="Tol(X)" *SDx="SD(X)" * testtype="Type", SD="SD(Resid)" * &spec * &result /rtspace=28; %end; /* TablType = 1betaOLS */ %if ((&TablTyp2 = FindNPow) and (&TablType ne none)) %then %do; /* Tables the powers after doing FindN problem */ proc tabulate format = 4.3 order=data; class TitlDlta alpha EffcTitl testtype NomPower parent scenario; format alpha &AlphaFmt; format NomPower &PowerFmt; var power; table TitlDlta="", EffcTitl= "Test" * testtype="Type" &AddRow2, alpha * NomPower="Minimum Power" * Power="Power"*mean=" " /rtspace=28; %end; /* FindNPow */ %donetabl: run; %mend tables; %macro FtntDago; footnote1 '*The Approximate Unconditional "chi**2" test corresponds to the'; footnote2 "common (Pearson) chi-square test for a 2 x 2 table. However,"; footnote3 "the method employed here uses a regular t test with Y = 0 (no)"; footnote4 "or 1 (yes), which has been shown to give more accurate p-levels"; footnote5 "than the common test and can be computed with any t test routine."; footnote6 "See D'Agostino, et. al. (1988), Am. Statistician, 42:198-202."; footnote7 " "; %mend FtntDago; %macro FtntSuis; footnote8 "**The Exact Unconditional corresponds to the test proposed by "; footnote9 "Suissa and Shuster (1985), J Royal Stat Soc A, 148:317-327). "; %mend FtntSuis; /* ++++ End UPw02_Tables segment ++++ */ /* ++++ Begin UPw03_MiscMacros segment v. 2002.07.01 ++++ */ %macro ChkNull (ChkVar, Endxx); /* Macro to check for null scenario (effect size = 0). If so, sets power = alpha and skips regular computation. Required to trap cases with no power, so that FindN will work OK. */ if &ChkVar = 0 then do; if not(FindingN) then power = alpha; if FindingN then do; if NulMsgOn then put // "============================================================" / "NOTE: There appears to be a situation in which the effect is" / " a null one, thus no Ntotal will achieve desired power." / " Ntotal has been set to missing (.)." / "============================================================" //; Ntotal = .; if piProb = -1 then do; TruAlpha = .; HiCrit = .; LoCrit = .; end; NulMsgOn = 0; end; /* FindingN */ go to &Endxx; end; /* ChkVar */ %mend ChkNull; %macro CheckWp (Wp1, Wp2, Wp3or4); /* Checks the Wilcoxon parameters and takes appropriate actions. */ if ((NumGrps = 1) and ((&Wp2 < &Wp1) or (&Wp2 < &Wp3or4))) or ((NumGrps = 2) and ((&Wp2 ge &Wp1) or (&Wp3or4 ge &Wp1))) then do; DoLH = 0; DoNoe = 1; if NumGrps = 1 then put // "WARNING: A situation makes p2 < p1 or p2 < p4."; if NumGrps = 2 then put // "WARNING: A situation fails to make p2 < p1 and p3 < p1."; put " Only Noether's method will be used."; end; %mend CheckWp; %macro CkCntrst (Cmat, Nrows, Ncols); do irowCC = 1 to &Nrows; sumCC = 0; do icolCC = 1 to &Ncols; sumCC = sumCC + &Cmat{irowCC, icolCC}; end; if abs(sumCC) gt .01 then put // "WARNING: Contrast row " irowCC " does not sum to 0.00."; end; %mend CkCntrst; %macro instruct (InfileNm=statmnts); /* Simple utility for input used when UnifyPow is a macro. */ data _null_; file &InfileNm; input; put _infile_; %mend instruct; %macro TroubMsg; put // "WARNING: Despite extensive testing, the iterative method" / " to find NTtotal for a specified power has failed" / " to converge. NTotal was set to missing." // " Please report this to" / " Ralph O'Brien" / " robrien@bio.ri.ccf.org"; Ntotal = .; power = .; go to &Cntinu; %mend TroubMsg; /* ++++ End UPw03_MiscMacros segment ++++ */ /* ++++ Begin UPw04_FindNMacros segment v. 2002.07.01 ++++ */ %macro FindN (GetPow, Cntinu); /* Secant method to find Ntotal such that power(Ntotal) = NomPower. Starting values are carefully found to ensure convergence. GetPow is the statement label to begin power computation. Cntinu is the statement label to go to after convergence. */ if Step = 1 then do; FindingN = 1; MinPossN = NumGrps + 1; if rhoProbT = 1 then MinPossN = 3; if rhoProb and (rhoProbT = 0) then MinPossN = NumGrps*4; if RSqProb or betLSPrb or SSH_Prob or ANCVProb then MinPossN = rXfull + 1; Ntotal = MinPossN; Step = 2; go to &GetPow; end; /* Step 1 */ if Step = 2 then do; if (Ntotal = MinPossN) and (power > NomPower) then do; put // "WARNING: Power exceeds " NomPower / " using smallest possible NTotal: " MinPossN / " Are your specifications correct?"; Ntotal = .; go to &Cntinu; end; Diff = power - NomPower; if Diff < 0 then do; Ntotal = 2*Ntotal; OldDifLo = Diff; go to &GetPow; end; NtotHi = Ntotal; DiffHi = Diff; NtotLo = Ntotal/2; DiffLo = OldDifLo; Ntotal = NtotLo + (NtotHi - NtotLo)/(1 - DiffHi/DiffLo); if not(Ntotal > 0) then do; %TroubMsg; end; Step = 3; go to &GetPow; end; /* Step 2 */ if Step = 3 then do; Diff = power - NomPower; if abs(Diff) < .00049 then do; FindingN = 0; go to &Cntinu; end; if Diff < 0 then do; NtotLo = Ntotal; DiffLo = Diff; end; else do; NtotHi = Ntotal; DiffHi = Diff; end; Ntotal = NtotLo + (NtotHi - NtotLo)/(1 - DiffHi/DiffLo); if not(Ntotal > 0) then do; %TroubMsg; end; go to &GetPow; end; /* Step 3 */ %mend FindN; %macro FindNBnl (GetPow, Cntinu); WndwSize = 10; /* half-size of window for point-by-point checking. */ /* Finds minimum N = Ntotal for binomial such that power(N) > NomPower. GetPow is the statement label to begin power computation. Cntinu is the statement label to go to after convergence. */ if Step = 1 then do; FindingN = 1; N_try = 1; Step = 2; go to &GetPow; end; /* Step 1 */ if Step = 2 then do; if (N_try = 1) and (power > NomPower) then do; put // "WARNING: Power of " NomPower "is reached with NTotal = 1." / " Are your specifications correct?"; FindingN = 0; go to &Cntinu; end; if power < NomPower then do; N_try = N_try*2; go to &GetPow; end; N_tryHi = N_try; N_tryLo = N_try/2; if (N_tryHi - N_tryLo) < 5 then Step = 4; else do; Step = 3; N_try = N_tryHi - .5*(N_tryHi - N_tryLo); go to &GetPow; end; end; /* Step 2 */ if Step = 3 then do; if power > NomPower then N_tryHi = N_try; else N_tryLo = N_try; if (N_tryHi - N_tryLo) < WndwSize then do; /* At this point, we know that N_tryLo gives insufficient power, and it is within WndwSize of a value (N_tryHi) that is sufficient. If the N vs. power function were monotonically increasing, we could just increase N_tryLo one at a time to find the minimally sufficient N. But the discreteness of the binomial makes the N vs. power function not perfectly monotone so we need to also check if values a little less than N_tryLo might satisfy power > NomPower. */ Step = 4; N_try = N_tryLo - WndwSize; go to &GetPow; end; else do; N_try = N_tryHi - .5*(N_tryHi - N_tryLo); go to &GetPow; end; end; /* Step 3 */ if Step = 4 then do; if power > NomPower then do; FindingN = 0; go to &Cntinu; end; N_try = N_try + 1; go to &GetPow; end; /* Step 4 */ %mend FindNBnl; /* ++++ End UPw04_FindNMacros segment ++++ *//* ++++ Begin UPw05_MtrxMacros segment v. 2002.07.01 ++++ */ /* Set of matrix routines for the data step: matDup ....... C = A %matAmnB(A, NrowA, NcolA, B, C) matAmnB ...... C = A - B %matAmnB(A, NrowA, NcolA, B, C) C may replace A or B matDiag ...... C = diag(vecA). Makes diagonal matrix from vector %matDiag(vecA, dim_vecA, C) matDiagI ...... C = diag(vecA). Makes diagonal matrix from reciprocals of elements of a vector %matDiagI(vecA, dim_vecA, C) mat_kA ....... C = k*A, where k is a scalar %mat_kA(A, NrowA, NcolA, k, C) C may replace A matAB ........ C = A*B %matAB(A, NrowA, NcolA, B, NcolB, C) C may replace A, not B matABt ....... C = A*B` %matABt(A, NrowA, NcolA, B, NrowB, C) C may replace A, not B. matAtB ....... C = A`*B %matAtB(A, NrowA, NcolA, B, NcolB, C) C may NOT replace A or B. matAtA ....... C = A`*A %matAtA(A, NrowA, NcolA, C) C may NOT replace A. matAv ....... C = A*vecB, where v is a 1-dimensional array %matAv(A, NrowA, NcolA, vecB, C) C may replace A. matAinv ...... C = inv(A) %matAinv(A, dimA, C, detA) C may replace A matHalf ...... finds L s.t. LL` = A, where A is PDS, L is lower triangular %matHalf(A, dimA, L) L may replace A matInput ...... inputs matrix using UnifyPow syntax %matInput (NewMatrx, Nrows, Ncols) CkCntrst checks whether rows of sum to 0.00 %CkCntrst (C, Nrows, Ncols) matEigVl: .... finds eigenvalues of symmetric matrix, EHE, of dimension dimEHE result is vector, phi. (difficult to make into macro) */ %macro matAmnB(A, NrowA, NcolA, B, C); /* Computes C = A - B, where all are NrowA x NcolA matrices. C may replace either A or B. */ do i_rowA = 1 to &NrowA; do i_colA = 1 to &NcolA; &C{i_rowA, i_colA} = &A{i_rowA, i_colA} - &B{i_rowA, i_colA}; end; end; %mend matAmnB; %macro matDiag(d, dim_d, A); /* Creates an dim_d x dim_d diagonal matrix, A, based on the vector d. */ do i_rowA = 1 to &dim_d; do i_colA = 1 to &dim_d; if i_rowA = i_colA then &A{i_rowA, i_colA} = &d{i_rowA}; else &A{i_rowA, i_colA} = 0; end; end; %mend matDiag; %macro matDiagI(d, dim_d, A); /* Creates an dim_d x dim_d diagonal matrix, A, based on the reciprocals of the elements in the vector d. */ do i_rowA = 1 to &dim_d; do i_colA = 1 to &dim_d; if i_rowA = i_colA then &A{i_rowA, i_colA} = 1/&d{i_rowA}; else &A{i_rowA, i_colA} = 0; end; end; %mend matDiagI; %macro mat_kA(A, NrowA, NcolA, k, kA); /* Multiplies an NrowA x NcolA matrix by a scalar k. Returns kA, an NrowA x NcolA matrix. kA may replace A. */ do i_rowA = 1 to &NrowA; do i_colA = 1 to &NcolA; &kA{i_rowA, i_colA} = &k*&A{i_rowA, i_colA}; end; end; %mend mat_kA; %macro matAB(A, NrowA, NcolA, B, NcolB, AB); /* Multiplies an NrowA x NcolA matrix by an NcolA x NcolB matrix Returns AB, an NrowA x NcolB matrix. AB may replace A. Must declare TempVec{max NcolB} in the main section. */ do i_rowA = 1 to &NrowA; do i_colB = 1 to &NcolB; TempVec{i_colB} = 0; do i_colA = 1 to &NcolA; TempVec{i_colB} = TempVec{i_colB} + &A{i_rowA, i_colA}*&B{i_colA, i_colB}; end; end; do i_colB = 1 to &NcolB; &AB{i_rowA, i_colB} = TempVec{i_colB}; end; end; %mend matAB; %macro matAv(A, NrowA, NcolA, v, Av); /* Multiplies NrowA x NcolA matrix by NcolA-element vector. Returns Av, an NrowA x 1 matrix. Av may replace A. */ do i_rowA = 1 to &NrowA; temp = 0; do i_colA = 1 to &NcolA; temp = temp + &A{i_rowA, i_colA}*&v{i_colA}; end; &Av{i_rowA, 1} = temp; end; %mend matAv; %macro matABt(A, NrowA, NcolA, B, NrowB, ABt); /* Multiplies an NrowA x NcolA matrix by the transpose of an NrowB x NcolA matrix. Returns ABt, an NrowA x NrowB matrix. ABt may replace A. */ do i_rowA = 1 to &NrowA; do i_rowB = 1 to &NrowB; tempVec{i_rowB} = 0; do i_colA = 1 to &NcolA; tempVec{i_rowB} = tempVec{i_rowB} + &A{i_rowA, i_colA}*&B{i_rowB, i_colA}; end; end; do i_colA = 1 to &NrowB; &ABt{i_rowA, i_colA} = tempVec{i_colA}; end; end; %mend matABt; %macro matAtB(A, NrowA, NcolA, B, NcolB, AtB); /* Multiplies the transpose of an NrowA x NcolA matrix by an NrowA x NcolB matrix Returns AtB, an NcolA x NcolB matrix. AtB may NOT replace A. */ do i_colA = 1 to &NcolA; do i_colB = 1 to &NcolB; &AtB{i_colA, i_colB} = 0; do i_rowA = 1 to &NrowA; &AtB{i_colA, i_colB} = &AtB{i_colA, i_colB} + &A{i_rowA, i_colA}*&B{i_rowA, i_colB}; end; end; end; %mend matAtB; %macro matAtA(A, NrowA, NcolA, AtA); /* A is NrowA x NcolA matrix. Returns A'A, an NcolA x NcolA matrix. AtA may NOT replace A. Must declare tempVec{max NcolB} in the main section. */ do i_rowAtA = 1 to &NcolA; do i_colAtA = i_rowAtA to &NcolA; &AtA{i_rowAtA, i_colAtA} = 0; do i_rowA = 1 to &NrowA; &AtA{i_rowAtA, i_colAtA} = &AtA{i_rowAtA, i_colAtA} + &A{i_rowA, i_rowAtA}*&A{i_rowA, i_colAtA}; end; if i_rowAtA < i_colAtA then &AtA{i_colAtA, i_rowAtA} = &AtA{i_rowAtA, i_colAtA}; end; end; %mend matAtA; %macro matAinv(Z, dimZ, Zinv, detZ); /* Z is square matrix of dimension dimZ. Returns inv(Z) and det(Z). Z is not altered, unless Zinv replaces Z, which is OK. */ %matDup(&Z, &dimZ, &dimZ, TempZ); %matAinv2(TempZ, &dimZ, &Zinv, &detZ); %mend matAinv; %macro matAinv2(A, dimA, Ainv, detA); /* A is square matrix of dimension dimA. Returns inv(A) and det(A). Ainv may NOT replace A. A is altered. Need to declare vv{max dimA} and indx{max dimA}. Algorithm taken from Press WH, et. al. Numerical recipes in FORTRAN: the art of scientific computing, 2nd Ed., Cambridge University Press, 1986, pp 34-41. */ tiny = 1.0E-12; do i = 1 to &dimA; do j = 1 to &dimA; &Ainv{i,j} = 0; end; &Ainv{i,i} = 1; end; /* LU decomposition: adapted from ludcmp subroutine in Numerical Recipes */ do i = 1 TO &dimA; maxA = 0; do j = 1 TO &dimA; if abs(&A{i,j}) > maxA then maxA = abs(&A{i,j}); end; if maxA = 0 then do; put "ERROR: A matrix singularity has suspended computations."; stop; end; vv{i} = 1.0/maxA; end; d = 1; do j = 1 to &dimA; do i = 1 to j-1; sum = &A{i,j}; do k = 1 to i-1; sum = sum - &A{i,k}*&A{k,j}; end; &A{i,j} = sum; end; maxA = 0; do i = j to &dimA; sum = &A{i,j}; do k = 1 to j-1; sum = sum - &A{i,k}*&A{k,j}; end; &A{i,j} = sum; dum = vv{i}*abs(sum); if dum >= maxA then do; max_i = i; maxA = dum; end; end; if j ne max_i then do; do k = 1 to &dimA; dum = &A{max_i, k}; &A{max_i,k} = &A{j,k}; &A{j,k} = dum; end; d = -d; vv{max_i} = vv{j}; end; indx{j} = max_i; if &A{j,j} = 0 then &A{j,j} = tiny; if j ne &dimA then do; dum = 1.0/&A{j,j}; do i = j+1 to &dimA; &A{i,j} = &A{i,j}*dum; end; end; end; /* j = 1 to dimA */ /* Compute determinant to see A if is singular. To understand this, read in Numerical Recipes. */ &detA = d; do i = 1 to &dimA; &detA = &detA*&A{i,i}; end; if abs(&detA) < tiny then do; put // "ERROR: A singular or near-singular matrix was encountered" / " when computing an inverse."; &detA = 0; do irowAinv = 1 to &dimA; do icolAinv = 1 to &dimA; &Ainv{irowAinv, icolAinv} = .; end; end; EndAinv = 1; end; if not(EndAinv) then do; /* functions as go to */ do m = 1 to &dimA; /* LU backsubstitution: adapted from lubksb subroutine in Numerical Recipes */ ii = 0; do i = 1 to &dimA; LL = indx{i}; sum = &Ainv{LL,m}; &Ainv{LL,m} = &Ainv{i,m}; if ii ne 0 then do j = ii TO (I-1); sum = sum - &A{i,j}*&Ainv{j,m}; end; else if sum ne 0 then ii = I; &Ainv{i,m} = sum; end; do i = &dimA to 1 by -1; sum = &Ainv{i,m}; do j = i+1 to &dimA; sum = sum - &A{I,j}*&Ainv{j,m}; end; &Ainv{i,m} = sum/&A{i,i}; end; end; end; /* endAinv */ %mend matAinv2; %macro matHalf(A, dimA, L); /* Choleski (square root) decomposition of dimA x dimA PDF symmetric matrix: A = LL`. Returns L, a dimA x dimA lower triangular matrix. You may set A and L to be the same, thus replacing A with L. Method adapted from choldc subroutine in Press WH, et. al. Numerical recipes in FORTRAN: the art of scientific computing, 2nd Ed., Cambridge University Press, 1986, p 90. Modified slightly to put results in L matrix and not use 'p' vector. */ tiny = 1.0E-12; do i_rowA = 1 to &dimA; do i_colA = i_rowA to &dimA; sum = &A{i_rowA, i_colA}; do k = i_rowA-1 to 1 by -1; sum = sum - &A{i_rowA, k}*&A{i_colA, k}; end; if (i_rowA = i_colA) then do; if (abs(sum) < tiny) then sum = 0; if (sum le 0) then do; put // "ERROR: Macro matHalf failed because 'A' is not positive definite."; stop; end; &L{i_rowA, i_rowA} = sqrt(sum); end; else &L{i_colA, i_rowA} = sum/&L{i_rowA, i_rowA}; end; end; do i_rowA = 1 to &dimA-1; do i_colA = i_rowA+1 to &dimA; &L{i_rowA, i_colA} = 0; end; end; %mend matHalf; %macro matDup(OldMatrx, NumRows, NumCols, NewMatrx); /* OldMatrix (NumRows x NumCols) is duplicated as NewMatrx. */ do irowSM = 1 to &NumRows; do icolSM = 1 to &NumCols; &NewMatrx{irowSM, icolSM} = &OldMatrx{irowSM, icolSM}; end; end; %mend matDup; %macro matPrint(A, NrowA, NcolA, label, format=best10. +2); /* Prints a matrix A with a given format for each value. Matrix is preceeded by a label. */ put // "&label"; do irow = 1 to &NrowA; put / +5 @; do icol = 1 to &NcolA; put &A{irow, icol} &format @; end; end; %mend matPrint; %macro vecPrint(A, dimA, label, format= best10. +2); /* Prints a vector A with a given format for each value. Matrix is preceeded by a label. */ put // "&label" / +5 @; do irow = 1 to &dimA; put &A{irow} &format @; end; %mend vecPrint; %macro matInput (NewMatrx, Nrows, Ncols); link GetMatrx; &Nrows = Nrows; &Ncols = Ncols; %matDup(Mvalues, &Nrows, &Ncols, &NewMatrx); %mend matInput; /* ++++ End UPw05_MtrxMacros segment ++++ */ /* ++++ Begin UPw06_BetaBnmlMacros segment v. 2002.07.01 ++++ */ %macro SDBetBnl (muP, SDP, Trials, p_beta, q_beta, SDBB); /* SDBB is std dev of Beta-Binomial outcome, Y/Trials, defined as follows. Subject's "true score" success rate, P, is distributed as standard beta random variable with mean MuP and std dev SDP. A check is made that this beta density is unimodal. Given P, Y is dist'd as binomial(Trials, P). Y is called a beta-binomial random variable, also known as the negative hypergeometric. The subject's success probability, Y/Trials, is the random variable of interest. */ /* Theory in Johnson, Kotz, Balakrishnan, [Continuous Univariate Distributions, Vol. 2, 2nd Ed., 1995], Equations 25.28 & 25.29. Note: JKB's "p" = p_beta and "q" = q_beta. */ pPLUSq = &muP*(1 - &muP)/&SDP**2 - 1; &p_beta = (&muP**2)*(1 - &muP)/(&SDP**2) - &muP; &q_beta = pPLUSq - &p_beta; if (&p_beta lt 1) or (&q_beta lt 1) then do; put // "ERROR: With mu(P) = " &muP "and SD(P) = " &SDP; put " the resulting beta distribution for P is not unimodal."; stop; end; /* Theory in Johnson, Kotz, and Kemp [Univariate Discrete Distributions, 2nd Ed., 1992], Equation 6.48. Note: JKK's "alpha" = p_beta and "beta" = q_beta. */ &SDBB = &Trials*&p_beta*&q_beta*(&p_beta+&q_beta+&Trials); &SDBB = &SDBB/(((&p_beta+&q_beta)**2)*(&p_beta+&q_beta+1)); &SDBB = sqrt(&SDBB)/&Trials; %mend SDBetBnl; %macro ProbBBnl(X, p_beta, q_beta, N, PrX); /* Computes probabilities for beta-binomial random variable X: PrX = Prob[X = x | P, N] where N fixed, P distd as beta(p_beta, q_beta). Uses (6.18) in Johnson, Kotz, and Kemp [Univariate Discrete Distributions, 2nd Ed., 1992]. */ %NBnlCoef(&p_beta, &X, result1); %NBnlCoef(&q_beta, &N-&X, result2); %NBnlCoef(&p_beta+&q_beta, &N, result3); &PrX = (result1/result3)*result2; %mend ProbBBnl; %macro NBnlCoef (NBn, NBr, NBresult); /* Finds /-n\ \ r/, using Eq. 1.10 in Johnson, Kotz, Kemp, 1992. */ %BnmlCoef(&NBn+&NBr-1, &NBr, &NBresult); &NBresult = ((-1)**(&NBr))*&NBresult; %mend NBnlCoef; %macro BnmlCoef(BCn, BCr, n_over_r); /* Finds /n\ \r/, using gamma function. */ &n_over_r = log(gamma(&BCn+1)) - log(gamma(&BCn-(&BCr)+1)) - log(gamma(&BCr+1)); &n_over_r = exp(&n_over_r); %mend BnmlCoef; %macro SetCumBB(p, q, Trials, CumDist); /* Sets the cumulative distribution function of a beta-binomial with parameters P ~ beta(p, q) and Y ~ binomial(Trials, P). */ %let i = i_CumDst; do &i = 0 to &Trials; %ProbBBnl(&i, &p, &q, &Trials, prob_i); if &i = 0 then &CumDist{0} = prob_i; else &CumDist{&i} = &CumDist{&i-1} + prob_i; end; %mend SetCumBB; %macro StWlBB2G(p_betaX, q_betaX, p_betaY, q_betaY, Trials, Wp1, Wp2, Wp3); /* Set p1, p2, p3 parameters for 2-group Wilcoxon problem for beta-binomial parent. Lots of ties require adjustments: p1 = Pr[Y > X] + .5*P[Y = X], p2 = Pr[{Yi > Xk} and {Yi' > Xk}] + 0.50*Pr[{Yi = Xk} and {Yi' > Xk}] + 0.25*Pr[{Yi = Xk} and {Yi' = Xk}], p3 = Pr[{Yi > Xk} and {Yi > Xk'}] + 0.50*Pr[{Yi = Xk} and {Yi > Xk'}] + 0.25*Pr[{Yi = Xk} and {Yi = Xk'}]. */ %let pX = pX_BB2G; %let qX = qX_BB2G; %let pY = pY_BB2G; %let qY = qY_BB2G; %let muX = muX_BB2G; %let muY = muY_BB2G; &pX = &p_betaX; &qX = &q_betaX; &pY = &p_betaY; &qY = &q_betaY; redoBB2G: &muX = &pX/(&pX+&qX); &muY = &pY/(&pY+&qY); %SetCumBB(&pX, &qX, &Trials, cumX); %SetCumBB(&pY, &qY, &Trials, cumY); %let X = X_SetCBB; %let Y = Y_SetCBB; do &Y = 0 to &Trials; if &Y = 0 then &Wp1 = .5*cumY{0}*cumX{0}; else &Wp1 = &Wp1 + (cumY{&Y} - cumY{&Y-1})*(0.5*(cumX{&Y} - cumX{&Y-1}) + cumX{&Y-1}); end; /* check if Wp1 makes sense */ &Wp1 = round(&Wp1, .001); if (&muX le &muY) and (&Wp1 < 0.500) or (&muX ge &muY) and (&Wp1 > 0.500) then do; put // "WARNING: The specified beta-binomial distributions"; put " are such that power computations for the"; put " Wilcoxon test would be misleading. Try"; put " increasing the effect size."; &Wp1 = .; &Wp2 = .; &Wp3 = .; go to outBB2G; end; /* UnifyPow is set only for Wp1 > 0.50 */ else if (&muX ge &muY) and (&Wp1 < 0.500) then do; /*"Reflect" problem by switching X and Y */ temp = &pX; &pX = &pY; &pY = temp; temp = &qX; &qX = &qY; &qY = temp; go to redoBB2G; end; do &X = 0 to &Trials; if &X = 0 then &Wp2 = cumX{0}*(.25*cumY{0}**2 + .5*2*cumY{0}*(1-cumY{0}) + (1-cumY{0})**2); else &Wp2 = &Wp2 + (cumX{&X} - cumX{&X-1})*(.25*(cumY{&X} - cumY{&X-1})**2 + .5*2*(cumY{&X} - cumY{&X-1})*(1-cumY{&X}) + (1-cumY{&X})**2); end; do &Y = 0 to &Trials; if &Y = 0 then &Wp3 = .25*cumY{0}*cumX{0}**2; else &Wp3 = &Wp3 + (cumY{&Y} - cumY{&Y-1})*(.25*(cumX{&Y} - cumX{&Y-1})**2 + .5*2*(cumX{&Y} - cumX{&Y-1})*cumX{&Y-1} + cumX{&Y-1}**2); end; outBB2G: %mend StWlBB2G; %macro StWlBB1G(p_beta, q_beta, Trials, NullVal, Wp1, Wp2, Wp4); /* Set p1, p2, p4 parameters for 1-group Wilcoxon problem for beta-binomial parent. Lots of ties require adjustments: p1 = Pr[Y > muYnull] + .5*P[Y = muYnull] > .50 p2 = Pr[(Yi + Yi')/2 > muYnull] + 0.50*Pr[(Yi + Yi')/2 = muYnull], p4 = Pr[{(Yi + Yi')/2 > muYnull)} and {(Yi + Yi'')/2 > muYnull)}] + 0.50*Pr[{(Yi + Yi')/2 > muYnull)} and {(Yi + Yi'')/2 = muYnull)}] + 0.25*Pr[{(Yi + Yi')/2 = muYnull)} and {(Yi + Yi'')/2 = muYnull)}]. */ %let pY = pY_BB1G; %let qY = qY_BB1G; %let NullValu = NV_BB1G; %let muY = muY_BB1G; &pY = &p_beta; &qY = &q_beta; &NullValu = &NullVal; if &NullValu = 0 then do; put // "ERROR: NULL of 0 is not allowed."; stop; end; if &NullValu = 1 then do; put // "ERROR: NULL of 1 is not allowed."; stop; end; SetCumY: %SetCumBB(&pY, &qY, &Trials, cumY); &muY = &pY/(&pY+&qY); %let muYnull = muYnull_; &muYnull = &Trials*&NullValu; if floor(&muYnull) = &muYnull then &Wp1 = .5*(cumY{&muYnull} - cumY{&muYnull - 1}) + (1 - cumY{&muYnull}); else &Wp1 = 1 - cumY{floor(&muYnull)}; /* check if Wp1 makes sense */ &Wp1 = round(&Wp1, .001); if ((&muY ge &NullValu) and (&Wp1 < 0.500)) or ((&muY le &NullValu) and (&Wp1 > 0.500)) then do; put // "WARNING: The specified beta-binomial distribution"; put " is such that power computations for the"; put " Wilcoxon test would be misleading. Try"; put " increasing the effect size."; &Wp1 = .; &Wp2 = .; &Wp4 = .; goto outBB1G; end; if ((&muY < &NullValu) and (&Wp1 < 0.500)) then do; /* UnifyPow is set only for Wp1 > 0.50, so "reflect" the problem. */ temp = &pY; &pY = &qY; &qY = temp; &NullValu = 1 - &NullValu; go to SetCumY; end; %let Y1 = Y1_SetBB; %let MGvY1 = MGvY1_; %let Wp2GvY1 = Wp2GvY1_; do &Y1 = 0 to &Trials; &MGvY1 = 2*&muYnull - &Y1; if &MGvY1 < 0 then &Wp2GvY1 = 1; else if &MGvY1 > &Trials then &Wp2GvY1 = 0; else if &MGvY1 = 0 then &Wp2GvY1 = .5*cumY{0} + (1 - cumY{&MGvY1}); else if floor(&MGvY1) = &MGvY1 then &Wp2GvY1 = .5*(cumY{&MGvY1} - cumY{&MGvY1 - 1}) + (1 - cumY{&MGvY1}); else &Wp2GvY1 = (1 - cumY{floor(&MGvY1)}); if &Y1 = 0 then &Wp2 = &Wp2GvY1*(cumY{0}); else &Wp2 = &Wp2 + &Wp2GvY1*(cumY{&Y1} - cumY{&Y1-1}); end; %let Wp4GvY1 = Wp4GvY1_; do &Y1 = 0 to &Trials; &MGvY1 = 2*&muYnull - &Y1; if &MGvY1 < 0 then &Wp4GvY1 = 1; else if &MGvY1 > &Trials then &Wp4GvY1 = 0; else if &MGvY1 = 0 then &Wp4GvY1 = .25*cumY{0}**2 + (1 - cumY{0})**2; else if floor(&MGvY1) = &MGvY1 then &Wp4GvY1 = .25*(cumY{&MGvY1} - cumY{&MGvY1 - 1})**2 + + 2*.5*(cumY{&MGvY1} - cumY{&MGvY1 - 1})*(1 - cumY{&MGvY1}) + (1 - cumY{&MGvY1})**2; else &Wp4GvY1 = (1 - cumY{floor(&MGvY1)})**2; if &Y1 = 0 then &Wp4 = &Wp4GvY1*(cumY{0}); else &Wp4 = &Wp4 + &Wp4GvY1*(cumY{&Y1} - cumY{&Y1-1}); end; outBB1G: %mend StWlBB1G; /* ++++ End UPw06_BetaBnmlMacros segment ++++ */ /* ++++ Begin UPw07_SetCatWlcxMacros segment v. 2002.07.01 ++++ */ %macro StWlCat1 (NullVal, Wp1, Wp2, Wp4); /* Set p1, p2, p3 parameters for 1-group Wilcoxon problem for general discrete parent specified by user. Lots of ties require adjustments: p1 = Pr[Y > NullVal] + .5*P[Y = NullVal] > .50 p2 = Pr[(Yi + Yi')/2 > NullVal] + 0.50*Pr[(Yi + Yi')/2 = NullVal], p4 = Pr[{(Yi + Yi')/2 > NullVal)} and {(Yi + Yi'')/2 > NullVal)}] + 0.50*Pr[{(Yi + Yi')/2 > NullVal)} and {(Yi + Yi'')/2 = NullVal)}] + 0.25*Pr[{(Yi + Yi')/2 = NullVal)} and {(Yi + Yi'')/2 = NullVal)}]. */ %let NullValu = NV_Cat1; &NullValu = &NullVal; if not(0 < &NullValu < NumCat) then do; put // "ERROR: NULL not within range of category values."; stop; end; do iCat = 1 to NumCat; if iCat = 1 then cumY{iCat} = ProbY{iCat}; else cumY{iCat} = cumY{iCat-1} + ProbY{iCat}; end; if abs(cumY{NumCat} - 1) > .001 then do; put // "ERROR: Distribution probabilities do not sum to 1.000"; stop; end; if floor(&NullValu) = &NullValu then &Wp1 = .5*(cumY{&NullValu} - cumY{&NullValu - 1}) + (1 - cumY{&NullValu}); else &Wp1 = 1 - cumY{floor(&NullValu)}; if (&Wp1 < 0.500) then do; &Wp1 = 1 - &Wp1; /* UnifyPow is set only for Wp1 > 0.50, so "reflect" the problem. */ do iCat = 1 to NumCat; if iCat = NumCat then cumX{iCat} = 1; /* using cumX as temp */ else cumX{iCat} = 1 - cumY{NumCat - iCat}; end; do iCat = 1 to NumCat; cumY{iCat} = cumX{iCat}; end; &NullValu = (1 + NumCat) - &NullValu; end; %let Y1 = Y1SetCat; %let MGvY1 = MGvY1_; %let Wp2GvY1 = Wp2GvY1_; do &Y1 = 1 to NumCat; &MGvY1 = 2*&NullValu - &Y1; if &MGvY1 < 1 then &Wp2GvY1 = 1; else if (&MGvY1 > NumCat) then &Wp2GvY1 = 0; else if &MGvY1 = 1 then &Wp2GvY1 = .5*cumY{1} + (1 - cumY{&MGvY1}); else if floor(&MGvY1) = &MGvY1 then &Wp2GvY1 = .5*(cumY{&MGvY1} - cumY{&MGvY1 - 1}) + (1 - cumY{&MGvY1}); else &Wp2GvY1 = (1 - cumY{floor(&MGvY1)}); if &Y1 = 1 then &Wp2 = &Wp2GvY1*(cumY{1}); else &Wp2 = &Wp2 + &Wp2GvY1*(cumY{&Y1} - cumY{&Y1-1}); end; %let Wp4GvY1 = Wp4GvY1_; do &Y1 = 1 to NumCat; &MGvY1 = 2*&NullValu - &Y1; if &MGvY1 < 1 then &Wp4GvY1 = 1; else if &MGvY1 > NumCat then &Wp4GvY1 = 0; else if &MGvY1 = 1 then &Wp4GvY1 = .25*cumY{1}**2 + (1 - cumY{1})**2; else if floor(&MGvY1) = &MGvY1 then &Wp4GvY1 = .25*(cumY{&MGvY1} - cumY{&MGvY1 - 1})**2 + + 2*.5*(cumY{&MGvY1} - cumY{&MGvY1 - 1})*(1 - cumY{&MGvY1}) + (1 - cumY{&MGvY1})**2; else &Wp4GvY1 = (1 - cumY{floor(&MGvY1)})**2; if &Y1 = 1 then &Wp4 = &Wp4GvY1*(cumY{1}); else &Wp4 = &Wp4 + &Wp4GvY1*(cumY{&Y1} - cumY{&Y1-1}); end; %mend StWlCat1; %macro StWlCat2(Wp1, Wp2, Wp3); /* Set p1, p2, p3 parameters for 2-group Wilcoxon problem for general discrete parents specified by user. Lots of ties require adjustments: p1 = Pr[Y > X] + .5*P[Y = X], p2 = Pr[{Yi > Xk} and {Yi' > Xk}] + 0.50*Pr[{Yi = Xk} and {Yi' > Xk}] + 0.25*Pr[{Yi = Xk} and {Yi' = Xk}], p3 = Pr[{Yi > Xk} and {Yi > Xk'}] + 0.50*Pr[{Yi = Xk} and {Yi > Xk'}] + 0.25*Pr[{Yi = Xk} and {Yi = Xk'}]. */ cumX{1} = ProbX{1}; cumY{1} = ProbY{1}; do iCat = 2 to NumCat; cumX{iCat} = cumX{iCat-1} + ProbX{iCat}; cumY{iCat} = cumY{iCat-1} + ProbY{iCat}; end; if abs(cumX{NumCat} - 1) > .001 then do; put // "ERROR: 1st distribution probabilities do not sum to 1.000"; stop; end; if abs(cumY{NumCat} - 1) > .001 then do; put // "ERROR: 2nd distribution probabilities do not sum to 1.000"; stop; end; %let X = XSetWCat; %let Y = YSetWCat; ReDoCat2: do &Y = 1 to NumCat; if &Y = 1 then &Wp1 = .5*cumY{1}*cumX{1}; else &Wp1 = &Wp1 + (cumY{&Y} - cumY{&Y-1})*(0.5*(cumX{&Y} - cumX{&Y-1}) + cumX{&Y-1}); end; /* UnifyPow is set only for Wp1 > 0.50 */ if (&Wp1 < 0.500) then do; /*"Reflect" problem by switching X and Y */ do iCat = 1 to NumCat; temp = CumX{iCat}; CumX{iCat} = CumY{iCat}; CumY{iCat} = temp; end; go to ReDoCat2; end; do &X = 1 to NumCat; if &X = 1 then &Wp2 = cumX{1}*(.25*cumY{1}**2 + .5*2*cumY{1}*(1-cumY{1}) + (1-cumY{1})**2); else &Wp2 = &Wp2 + (cumX{&X} - cumX{&X-1})*(.25*(cumY{&X} - cumY{&X-1})**2 + .5*2*(cumY{&X} - cumY{&X-1})*(1-cumY{&X}) + (1-cumY{&X})**2); end; do &Y = 1 to NumCat; if &Y = 1 then &Wp3 = .25*cumY{1}*cumX{1}**2; else &Wp3 = &Wp3 + (cumY{&Y} - cumY{&Y-1})*(.25*(cumX{&Y} - cumX{&Y-1})**2 + .5*2*(cumX{&Y} - cumX{&Y-1})*cumX{&Y-1} + cumX{&Y-1}**2); end; %mend StWlCat2; /* ++++ End UPw07_SetCatWlcxMacros segment ++++ */ /* ++++ Begin UPw08_Declares segment v. 2002.08.17a ++++ */ %global ProbType TablType TablTyp2 TablTyp3 SDType TblMmnts DistList OthScen; %global TrlName fnote WriteBox ResltVar format p3or4 DoTables MethLabl; %global SDorCV MultLabl CorrLabl InfileNm TablSDwt; %global AlphaFmt PowerFmt Pi1Foot2 Pi1Foot3; %let WriteBox = 0; %let InfileNm = cards; /* Activate for module version. */ /* Set overall maximum sizes for UnifyPow problems. */ %let MxGrpMes = 20; /* Max number of groups or measures */ %let MaxSDs = 20; /* Max number of specs for SDs or similar parameters */ %let MaxSDVV = 400; /* Make this the square of MaxSDs */ %let MaxAlphs = 10; /* Max number of alpha specs */ %let MaxNs = 100; /* Max number of Ntotals or powers */ %let ProbType = To be defined; %let OthScen = ; %let DistList = DistLst1 DistLst2; %let TablTyp2 = none; %let TablTyp3 = none; %let TablSDwt = no; %let fnote = none; %let SDorCV = SD; %let SDType = Std Dev; %let SDadLabl = Std Dev reduction; %let DltaPowr = .02; %let AddRow2 = ; %let WriteBox = %eval(&WriteBox+1); %let DoTables = yes; %let MethLabl = %str(EffcTitl = "Method"); %let SDname = Standard Deviation; %let MultLabl = x SD (SD Multiplier); %let CorrLabl = Corr(Y1, Y2); %let Pi1Foot2 = The note above describes how they are set for a non-directional; %let Pi1Foot3 = (two-tailed) test.; %let AlphaFmt = 6.3; %let PowerFmt = 6.3; data PowData; infile &InfileNm missover eof=EOFfound; file print; if _n_ = 1 then do; if &WriteBox = 1 then do; put /; put "--------------------------------------------------------------------"; put "| UnifyPow 2002.08.17a 2002 Copyright (c) by Ralph G. O'Brien |"; *put "| <<<< This version of UnifyPow is a work in progress. >>>> |"; put "| For information, see http://www.bio.ri.ccf.org/UnifyPow |"; put "| Report problems to robrien@bio.ri.ccf.org |"; put "--------------------------------------------------------------------"; end; put // "Verbatim specifications and processing notes" / "============================================" / ; end; keep alpha power NomPower SD CV NTotal EffcTitl dfH testtype tails; keep parent PrimSSHe scenario BetaWt SDx tolernce SDP Trials; keep Corr SDMult Wp1 Wp2 Wp3 Wp4 DistLst1 DistLst2 NullValu; keep TitlDlta TruAlpha LoCrit HiCrit SpcTitl3 rXfull SDwgt; keep CovPCorr SDredctn ConfLevl pi_1 pi_2 OddsRato ProbRato; array value {&MxGrpMes} _temporary_; /* temporary parameter values */ array mu {&MxGrpMes} _temporary_; array pi {&MxGrpMes} _temporary_; array pi1 {&MxGrpMes} _temporary_; array pi2 {&MxGrpMes} _temporary_; array pCatGGrp {&MxGrpMes,&MxGrpMes} _temporary_; /* Pr[Cat=col | Group=row] */ array jpCatGrp {&MxGrpMes,&MxGrpMes} _temporary_; /* Pr[Cat=col & Group=row] */ array pCat {&MxGrpMes} _temporary_; /* Pr[Cat=col] */ array jp_null {&MxGrpMes,&MxGrpMes} _temporary_; /* Pr[Cat=col & Group=row] under Ho */ array Cmat {&MxGrpMes,&MxGrpMes} _temporary_; array w {&MxGrpMes} _temporary_; array wDesign {&MxGrpMes} _temporary_; array alphaV {&MaxAlphs} _temporary_; array SDMultV {&MaxSDs} _temporary_; /* handles both SD and CV */ array CorrV {&MaxSDs} _temporary_; array SDadjV {&MaxSDs} _temporary_; array SDMultVV {&MaxSDVV} _temporary_; /* handles both SD and CV */ array CorrVV {&MaxSDVV} _temporary_; array SDV {&MaxSDs} _temporary_; array PriSDV {&MaxSDs} _temporary_; array CVV {&MaxSDs} _temporary_; array SDPV {&MaxSDs} _temporary_; array SDPVV {&MaxSDVV} _temporary_; /* do i = 1 to &MaxSDVV; SDPVV{i} = 0; end; Still need this?; */ array TrialsV {&MaxSDs} _temporary_; array TrialsVV {&MaxSDVV} _temporary_; /* do i = 1 to &MaxSDVV; TrialsVV{i} = 0; end; Still need this? */ array p_betaV{200} _temporary_; array q_betaV{200} _temporary_; array probX{1:25} _temporary_; array probY{1:25} _temporary_; array cumX {0:1000} _temporary_; array cumY {0:1000} _temporary_; array BetaWtV {&MaxSDs} _temporary_; array SDxV {&MaxSDs} _temporary_; array TolV {&MaxSDs} _temporary_; array NtotalV {&MaxNs} _temporary_; array NomPowrV {&MaxNs} _temporary_; array Cmu {&MxGrpMes,1} _temporary_; array Amat {&MxGrpMes,&MxGrpMes} _temporary_; array invCWCt {&MxGrpMes,&MxGrpMes} _temporary_; array indx{&MxGrpMes} _temporary_; array vv{&MxGrpMes} _temporary_; array DistList{2} $78 DistLst1 DistLst2 (" " " "); array Mvalues {&MxGrpMes, &MxGrpMes} _temporary_; array TempMat {&MxGrpMes,&MxGrpMes} _temporary_; array TempVec {&MxGrpMes} _temporary_; array TempZ {&MxGrpMes,&MxGrpMes} _temporary_; /* for %matAinv */ array MVmu {&MxGrpMes,&MxGrpMes} _temporary_; array CorrMat {&MxGrpMes,&MxGrpMes} _temporary_; array diagSD {&MxGrpMes, &MxGrpMes} _temporary_; array CovarMat {&MxGrpMes,&MxGrpMes} _temporary_; array Rmat {&MxGrpMes,&MxGrpMes} _temporary_; /* reduced model mat */ array invRtR {&MxGrpMes,&MxGrpMes} _temporary_; array diagW {&MxGrpMes,&MxGrpMes} _temporary_; array BetaMat {&MxGrpMes,&MxGrpMes} _temporary_; array ThetaMat {&MxGrpMes, &MxGrpMes} _temporary_; array Theta0 {&MxGrpMes, &MxGrpMes} _temporary_; array C_Ovrall {%eval(&MxGrpMes-1), &MxGrpMes} _temporary_; array A_Ovrall {%eval(&MxGrpMes-1), &MxGrpMes} _temporary_; array C_Ave {1, &MxGrpMes} _temporary_; array A_Ave {1, &MxGrpMes} _temporary_; array A_Identy {&MxGrpMes, &MxGrpMes} _temporary_; array invRtWR {&MxGrpMes, &MxGrpMes} _temporary_; array CBA {&MxGrpMes, &MxGrpMes} _temporary_; array invCRWRC {&MxGrpMes, &MxGrpMes} _temporary_; array PrimHmat {&MxGrpMes, &MxGrpMes} _temporary_; array HfInvASA {&MxGrpMes,&MxGrpMes} _temporary_; array EHE {&MxGrpMes,&MxGrpMes} _temporary_; array phi {&MxGrpMes} _temporary_; * exemplary eigenvalues; array e{&MxGrpMes} _temporary_; array ConfLevV {&MaxSDs} _temporary_; array ratio {&MxGrpMes} _temporary_; length statemnt $78; /* scratch variable to hold full statement */ length ProbStmt $78; /* holds entire problem statement */ length title $78; /* temporary variable for ReadTitl link */ length Scenario $78; length EffcTitl $40; /* effect title, up to 40 characters */ length vChar $20; length TitlDlta $78; format TitlDlta $78.; length SpcTitl3 $78; format SpcTitl3 $78.; length onechar $1; length keyword $25; length parent $8; length NomParnt $8; length TestType $11; format Wp1-Wp4 5.3; format NomPower 4.3; format alpha 4.3; length CurInput $15; length OldInput $15; length string1 $80; length AdString $40; EndJstIn = 0; /* END statement last thing read */ NewProb:; /* New Problem */ /***** set defaults on specifications *****/ /* indicator variables: 0=no, 1=yes */ _2piProb = 0; /* special syntax for 2 pi problem */ AlphaIn = 0; /* alpha levels input yet */ ANCVProb = 0; /* covariates added to means problem */ ARE = 1; /* asymp. rel. efficiency for given test vs. t */ BalTails = 0; /* Use optimal (not balanced) tails for binomial */ betLSPrb = 0; /* test single BetaWt in OLS reg model */ Chi2Prob = 0; /* Exemplary data method inputting Chi**2 values */ CnfLevIn = 0; /* Confidenece levels in for power CI problem */ CntrstIn = 0; /* Has CONTRASTS statement been read */ CorrIn = 0; /* correlations input yet */ CustXIn = 0; /* Custom X matrix for MVGLM module */ CVIn = 0; /* Coef of variations input yet(LogNormal problems) */ CVV{1} = 0; CvPCorIn = 0; /* Covariate parial correlations are in */ DoAdjstN = 1; /* Adjust found N to make it conform to cell weights */ DoARE = 0; /* Compute Wilcoxon power via asymp rel eff vs t */ DoChiSq = 0; /* do ChiSq instead of F in GetPower link */ DoHotell = 0; /* Compute power for Hotelling statistic */ DoLH = 1; /* Do Lehmann-Hettsmanperger for Wilcoxon power */ doLR = 1; /* Do LR test */ DoNoe = 0; /* Compute Wilcoxon power using Noether's crude approx */ DoPillai = 1; /* Compute power for Pillai statistic */ DoTails = 3; /* Do both 1- and 2-tailed tests (1 + 2 = 3) */ DoWilks = 0; /* Compute power for Wilks statistic */ EffectIn = 0; /* Has EFFECTS statement been read */ EndOFile = 0; EOFlink = ""; /* link, if any, to go to after hitting EOF */ ExemProb = 0; /* Exemplary data problem */ FindingN = 0; /* currently interating to find Ntotal */ ForceN = 0; /* Force Ntotal to be a multiple of CoreNtot */ GdFtProb = 0; /* test goodness of fit of prob distn */ GrpsCIn = 0; /* Has GROUPS contrast matrix been read */ LimitsIn = 0; /* Has LIMITS statement been read */ LogNProb = 0; /* means are from LogNormal distribution */ McNrProb = 0; /* McNemar test */ MeasCIn = 0; /* Has MEASURES contrast matrix been read */ MethodIn = 0; /* Has HETHODS statement been read */ muP2Spcl = 0; /* Special Wilcoxon for 2-grp MU(P) with NullValu ne 0 */ muProb = 0; /* test of means */ muP_Prob = 0; /* test diffs in mu(P)s, where Y bnml with random pi */ MVmuProb = 0; /* multivariate general linear model for cell means */ neighbor = 0; /* used in prBnml link function */ NewAlfTl = 0; /* new Alpha or Tails has just been set */ NexIn = 0; /* N for exemplary SSH or CHI**2 values in */ NmCovars = 0; /* no covariates added to means problem */ NmParmIn = 0; /* Number of parameters statement in */ NoEcho = 0; /* Temporarily stops echoing of input */ NoNotes = 0; /* Suppress printing of notes */ NoOveral = 0; /* skip overall test */ NtotalIn = 0; /* total N input yet */ NullEfct = 0; /* Is this effect a null one */ NullIn = 0; /* null value in yet */ NullValu = 0; /* null value for one- and two-group tests */ NulMsgOn = 1; /* If null is true, print message; only does one time */ NumCnfLv = 0; /* Number of confidence levels for power CI problem */ NumCV = 1; /* Number of coefficients of variation */ NumGrps = 1; /* Number of groups */ NumProp1 = 0; /* Number of possible proportions for group 1 */ NumProp2 = 0; /* Number of possible proportions for group 2 */ NumRSq = 0; /* Number of R**2 values entered. */ NumSDadj = 1; /* Number of SD adjustments for muANOVA problem */ NumSDmlt = 1; /* Number of SDMult factors; default is 1 */ NumTrls = 1; /* This allows looping on Trials even when irrelevant*/ OddsProb = 0; /* Odds ratio problem */ PaMuProb = 0; /* test comparing paired (2 correlated) means */ parent = "Normal"; /* Default that underlying distribution is Normal */ piProb = 0; /* test of proportions */ PowerIn = 0; /* power input yet */ PriSDin = 0; /* priors for SD in */ PrRtProb = 0; /* Probability ratio (relative risk, risk ratio) problem */ Prop1in = 0; /* proportions for group 1 have been read */ Prop2in = 0; /* proportions for group 2 have been read */ putSDwrn = 0; /* used SD or SIGMA, not SD(P); put note */ rhoProb = 0; /* test of simple correlations */ rhoProbT = 0; /* test of rho = 0 If so, uses t test */ RSqProb = 0; /* test comparing R**2 values */ SamFProb = 0; /* estimate power from sample F statistic */ SamTProb = 0; /* estimate power from sample t statistic */ SamCProb = 0; /* estimate power from sample Chi**2 statistic */ SampProb = 0; /* any kind of "Sam_Prob" problem */ ScenarIn = 0; /* Scenario title has been input */ RvrsTail = 0; /* Reversed tail for binomial calculations */ SDadjV{1} = 1; /* SD adjustment is default of 1.0 */ SDIn = 0; /* standard deviations input yet */ SDMultV{1} = 1; /* Default is 1.0, that is, just use SD specified. */ SDPIn = 0; /* std deviations for pi [mu(P) problem] input yet */ SDredIn = 0; /* SD reduced in for mumANCOVA problem */ SDvalid = 1; /* Is the SD statement valid for this problem */ SDwgt = 1; /* Prior weight for SD value */ SDxIn = 0; /* SDx values in for 1betaOLS problem */ SpclCrit = 0; /* Adapts crit_r link fnctn for getMcNpw link fnctn */ SpcTitl3 = "Critical values and actual alpha levels using binomial distribution."; SSH_Prob = 0; /* Exemplary data method inputting SSH values */ ThetaIn = 0; /* ThetaMat (null matrix) in yet */ TitlDlta = "Some actual powers are &DltaPowr greater than specified minimum powers."; TolIn = 0; /* Tolerence values in for 1betaOLS problem */ TrialsIn = 0; /* N_items or N_trials for mu(P) problem input yet */ ValidPrb = 0; /* Valid problem statement input yet */ wIn = 0; /* weights input yet */ WlcxCstm = 0; /* User specified p1, p2, and p3 or p4 for Wilcoxon */ WlcxOrCt = 0; /* User specified Ordered Categ pr's for Wilcoxon */ WlcxProb = 0; /* Wilcoxon one or 2-group test of location */ _2WCTPrb = 0; /* test independence in 2-way table */ /* ++++ End UPw08_Declares segment ++++ */ /* ++++ Begin UPw09_ProbPars1 segment v. 2002.07.01 ++++ */ /* Comment lines may be used anyplace. There are three kinds: Block of comment lines. /# Block of comments lines, at most 78 chars per line. If any contain ";" then use DATALINES4 statement and follow all UnifyPow input with ";;;;". #/ Single comment lines. // Comments lines, at most 78 chars per line. // Another comment line. // See note just above about using ";". // Either block comments or // comments may be used anyplace. Follow end of UnifyPow statement with open period ( . ) and add comment. Examples: pi .80 . Mystic Michelle claims 80% success rate. alpha .05 .0167 . 0.05/3 = 0.0167 for Bonferroni protection. */ ReadProb:; /* Read in problem and scenario statements */ CurInput = "ProbStmt"; link EchoInpt; EndJstIn = 0; input ProbStmt $char78. @; if not(ScenarIn) then scenario = ProbStmt; input @1 keyword $ @; keyword = upcase(keyword); /* Parse the END statement */ if index(keyword,'END') then do; input; EndJstIn = 1; go to NewProb; end; /* Parse the SCENARIO statement. This can happen either before of after parsing of problem statement. */ link ParsScen; /* Parse the mu and the mu:logNormal statements */ /* The "mu" statement leads to regular ANOVA (Normal-theory) tests on means. Unless blocked by a "NoOverall" command, UnifyPow automatically tests the overall hypothesis that all G means are equal. "Contrast" statements let users test any contrast of the form C*mu = 0, where C is a rank(C) x G contrast matrix and mu is the G x 1 vector of means. The "null" statement applies to one-group and two-group tests: Ho: mu_1 = NullValu Ho: mu_1 - mu_2 = NullValu ---- The mu:logNormal statement assumes that the means supplied are for positively-skewed data that have LogNormal distributions with homogeneous coefficients of variation. Then ln(Y) is Normal with homogeneous variances. This leads to some elegantly simple sample-size analyses. Reference: van Belle G, Martin DC, "Sample size as a function of coefficient of variation and ratio of means," American Statistician, 1993, 47:165-167. van Belle and Martin only considered the balanced two-group case, but the results easily generalize to the G-group, unbalanced case, as implemented in UnifyPow. Suppose that Y_ik is log-normally distributed with arithmetic (not geometric) mean mu_i and with a coefficient of variation, CV = SD_i/mu_i, that is common to all groups (homogeneous). Then, Y' = ln(Y_ik) is normal with mean mu'_i = ln(mu_i) - .50*SD'**2 and homogeneous standard deviation SD' = sqrt(ln(CV**2 + 1)). This provides a very simple way to handle positively-skewed data in which there is a "reference" group mean, mu_R and all study means are mu_i = f_i*mu_R. Sigma_i,the standard deviation of Y_ik, is proportional to the mu_i, i.e., sigma_i = CV*mu_i. Using ln(Y) for analysis renders the particular value of mu_R irrelevant: One can merely set mu_R == 1 and proceed. The problem is sufficiently defined by just specifying the G f_i values and the common CV. For example, in a 3-group case, we might have Group 1 being "usual care" and Groups 2 and 3 being two variations on an experimental treatment. If mu_2 = .80*mu_1 and mu_3 = .70*mu_1 then one would use mu:logNormal 1.00 0.80 0.70 . Means 2 & 3 are 80% & 70% of Grp 1. If the reference group is not being studied (perhaps an historical control is being used), then you could have something like mu:logNormal .80 .70 . Means are 80% & 70% of historical control. The NULL statement for one- and two-group tests conforms to the following: Ho: mu_1 = NullValu Ho: mu_1/mu_2 = NullValu In log-transformed terms, these are Ho: mu'_1 = ln(NullValu) - .50*SD'**2 Ho: mu'_1 - mu'_2 = ln(NullValu) The coefficient of variation (CV) works quite easily. Suppose that you expect 95% of the cases to fall within 20% of their group's mean, whatever that mean happens to be. Then the CV is about .20/2 = 0.10. In UnifyPow, it is easy to bracket that value, say with CV .075 .100 .125 One may combine mu:logNormal with Wilcoxon. Also, the same theory has been applied to the PairedMu problem, giving rise to the PairedMu:logNormal problem. */ if (keyword='MU' or index(keyword,'MEAN') or index(keyword,':LOG')) and not(index(keyword,'PAIR')) then do; muProb = 1; ValidPrb = 1; call symput('TablType','FPow'); link GetValus; NumGrps = count; rXfull = NumGrps; if NumGrps le 2 then call symput('TablType','tPow'); if not(index(keyword,':LOG')) then do; call symput('ProbType','mu'); call symput('SDType','Standard Deviation'); do i = 1 to NumGrps; mu{i} = value{i}; end; end; else do; LogNProb = 1; call symput('ProbType', 'mu:logNormal'); call symput('SDorCV', 'CV'); call symput('SDType', 'Coef of Variation'); SDvalid = 0; if NumGrps = 2 then NullOrig = 1; else NullOrig = .; do i = 1 to NumGrps; mu{i} = value{i}; /* mu{} are transformed below */ end; end; go to NextSpec; end; /* parse the PairedMu and PairedMu:logNormal statements */ /* The PairedMu problem handles G groups of paired means, (mu_j1, mu_j2), which are input as PairedMu mu_11 mu_12 > mu_21 mu_22 ... > mu_G1 mu_G2 with each ">" in column 1. These are transformed to G means of the difference score, (Y1 - Y2), giving mu{j} = mu_j1 - mu_j2. Only one pair of SDs is allowed, SD SD_1 SD_2 corresponding to the common "base" SDs for Y1 and Y2. A set of possible correlations between Y1 and Y2 is input via the CORR statement, Corr # # # # The base SD of (Y1 - Y2) is SDbase(Y1 - Y2) = sqrt(SD_1**2 + SD_2**2 - 2*Corr*SD_1*SD_2). SDbase is then varied using the SDMULTiplier statement, SDMULT # # # # which supplies a set of values, m, to adjust SDbase up or down via SD(Y1 - Y2) = m*SD_base(Y1 - Y2). If no SDMULTiplier statement is given, then m = 1. The TotalN (NTotal) statement is used in the usual way to give the total number of PAIRS of cases. Accordingly, the statements TotalN # # # # and TotalPairs # # # # are identical. (I.e., "TOTAL" is the operative key string for this statement.) ------------ The PairedMu:logNormal statement works the same way, albeit the distribution is assumed to be logNormal. The PairedMu:logNormal problem handles G groups of paired geometric means, (GeoMu_j1, GeoMu_j2), which are input as PairedMu:logNormal GeoMu_11 GeoMu_12 > GeoMu_21 GeoMu_22 ... > GeoMu_G1 GeoMu_G2 The hypothesis tested for the one-group test conforms to the following: Ho: GeoMu[Y1/Y2] = GeoMu_11/GeoMu_12 = 1 (default) where GeoMu[X] designates the geometric mean: exp{Mu[ln(X)]}, where Mu is the usual expected value. For the logNormal problem, the NULL statement lets the user work in the original scale, not in the log-scale. Thus, for the one-group PairedMu:logNormal problem, UnifyPow uses Ho: GeoMu_11/GeoMu_12 = NullValu Thus, the NullValu is the ratio of the geometric means, with "time=1" on the top. In log-transformed terms, we have Ho: Mu[ln(Y1) - ln(Y2)] = Mu[ln(Y1)] - Mu[ln(Y2)] = 0 or, if NULL is used, Ho: Mu[ln(Y1) - ln(Y2)] = ln(NullValu). The two-group case is constructed as Ho: (GeoMu_11/GeoMu_12)/(GeoMu_21/GeoMu_22) = 1 (default) and more generally as Ho: (GeoMu_11/GeoMu_12)/(GeoMu_21/GeoMu_22) = Nullvalu Thus, the NullValu is the ratio of the ratio of geometric means. In log-transformed terms, we have a default of Ho: Mu[ln(Y11) - ln(Y12)] - Mu[ln(Y21) - ln(Y22)] = 0 or, if NULL is used, Ho: Mu[ln(Y11) - ln(Y12)] - Mu[ln(Y21) - ln(Y22)] = log(Nullvalu) The CV statement must be used to give exactly two values for the base coefficient of variation, CV CV_1 CV_2 The CORR statement gives scenarios for the correlations (Corr') between ln(Y1) and ln(Y2). The CVMULTiplier statement parallels the SDMULTiplier statement above: CVMULT # # # # and so supplies a set of values m to adjust CV_1 and CV_2 up or down via CV_j{m} = m*CV_j Thus the SD of ln(Y1/Y2) = ln(Y1) - ln(Y2) is SDdiff'{m} = sqrt(SD'_1{m}**2 + SD'_2{m}**2 - 2*Corr'*SD'_1{m}*SD'_2{m}), where SD'_j = sqrt(ln((m*CV_j)**2 + 1)). */ if index(keyword,'PAIR') then do; PaMuProb = 1; ValidPrb = 1; if index(keyword, ':LOG') then LogNProb = 1; if not(LogNProb) then call symput('ProbType', 'PairedMu'); else do; call symput('ProbType', 'PairedMu:logNormal'); call symput('MultLabl', 'x CV (CV Multipier)'); call symput('CorrLabl', 'Corr(logY1, logY2)'); end; call symput('TablType', 'PaMuPow'); link GetValus; if (floor(count/2) ne count/2) then do; put // "ERROR: The PairedMu problem requires even number of means."; stop; end; NumGrps = count/2; rXfull = NumGrps; do i = 1 to NumGrps; if not(LogNProb) then mu{i} = value{2*i-1} - value{2*i}; else mu{i} = value{2*i-1}/value{2*i}; end; if NumGrps = 1 then do; /* check for more groups */ MoreGrps: link EchoInpt; input keyword @; if keyword = '>' then do; NumGrps = NumGrps + 1; rXfull = NumGrps; link GetValus; if count ne 2 then do; put // "ERROR: Means are not in pairs."; stop; end; if not(ScenarIn) then do; scenario = trim(ProbStmt)||';'; do i = 1 to 2; scenario = trim(scenario)||' '||compress(value{i}); end; end; if not(LogNProb) then mu{NumGrps} = value{1} - value{2}; else mu{NumGrps} = value{1}/value{2}; go to MoreGrps; end; /* keyword = '>' */ if keyword ne '>' then do; input @1 @; NoEcho = 1; end; end; /* check for more groups */ if NumGrps le 2 then NullOrig = 1; else NullOrig = .; go to NextSpec; end; /* PAIR parsing */ /* Parse the 1WILCOXON AND 2WILCOXON/MANN-WHITNEY statements. This requires input of parameters defined by Noether (1987, Sample size determination for some common nonparametric tests, JASA 82:645-647) and Lehmann (1975, Nonparametrics: Statistical Methods Based on Ranks, New York: Wiley, pp. 65-81, 164-171, 196, 400-401). These comments follow Hettsmansperger (1984, Statistical Inference Based On Ranks, New York: Wiley, pp. 47-62, 157-159), who presents methods identical to Lehmann's, but in a notation more suitable to SAS coding. For the one-group case, we take Y1, Y2, ...Yn to be i.i.d. continuous with a symmetric parent distribution centered at delta. The null hypothesis is Ho: delta = delta0, which equates to Ho: Prob[Y > delta0] = .5, Thus, the first parameter is the effect size, p1 = Prob[Y > delta0]. We take p1 > .5. It is the only parameter used by Noether. The Lehmann-Hettsmansperger method uses two more parameters, p2 = Prob[(Y1 + Y2)/2 > delta0], and p4 = Prob[{(Y1 + Y2)/2 > delta0} and {(Y1 + Y3)/2 > delta0}]. Note: Hettsmansperger's p3 = (p2 + p1**2)/2 is superfluous and so it is not an input parameter in UnifyPow. If the variance of Y is finite, p2 >= p1 and p4 < p2. Note that in the SAS coding, p1 is Wp1 (Wilcoxon p1), etc. For the two-group case, we take X1, X2, ..., Xm and Y1, Y2, ...,Yn to be independent random variables from parent distributions having identical spreads and shapes, but not necessarily symmetric. Let delta be the difference in the two medians, thus the null hypothesis is Ho: delta = 0.0, which equates to Ho: Prob[Y > X] = 0.5. Thus, the first parameter is the effect size p1 = Prob[Y > X]. We take p1 > .5. It is the only parameter used by Noether. The Lehmann-Hettsmansperger method uses two more parameters, p2 = Prob[{Yi > Xk} and {Yi' > Xk}], and p3 = Prob[{Yi > Xk} and {Yi > Xk'}]. Both p2 and p3 are less than p1. If the parent distribution is symmetric, p2 = p3. Note: The Lehmann-Hettsmansperger method employed here is NOT the same as the "Lehmann" method investigated by Lesaffre E, Scheys I, Frohlich J, Bluhmki E (1993, Calculation of power and sample size with bounded outcome scores, Statistics in Medicine, 12:1063-1078). That method is a crude one based on a Normal parent and Lehmann gave it almost as an aside without recommending it over his more complex one. The LH method used here is a nonparametric one, like Noether's, but more refined in that it uses all three "moments," not just p1. Thus it conforms to the Lesaffre, et. al. recommendations given in their final paragraph. For either the one or two-group case, UnifyPow requires only p1, the effect size. If one of the other parameters is not specified, it will automatically determine all of them based on theory for Normal, Logistic, and Laplace (double exponential) distributions, which are all symmetric, and have kurtoses of 0.0, 1.2, and 3.0, respectively. */ if (index(keyword,'WILC') gt 0) or index(keyword,'MANN') then do; call symput ('TablType', '1or2WlcxnPow'); call symput ('SDType', 'none'); if NoNotes = 0 then call symput ('TblMmnts','yes'); call symput ("DistList", " "); call symput ("OthScen", " "); SDIn = 1; NumSD = 1; SDV{1} = 1; WlcxCstm = 1; NomParnt = "Custom"; ValidPrb = 1; link GetValus; if count > 2 then do; SDvalid = 0; call symput ('TblMmnts','no'); end; if count = 0 then go to OrderCat; if count = 1 then do; p1OnlyIn = 1; NomParnt = "NORMAL"; /* Default is Normal parent */ WlcxCstm = 0; end; else p1OnlyIn = 0; do i = 1 to count; if not(0 < value{i} < 1) then go to WlcxErr; end; if (value{1} < .50) and (count = 1) then do; temp = 1 - value{1}; put // "WARNING: UnifyPow requires p1 > 0.50."; put " It has been reset to p1 = 1 - " value{1} "= " temp; value{1} = temp; end; else if (value{1} < .50) and (count > 1) then do; put // "ERROR: UnifyPow requires p1 > 0.50."; stop; end; if index(keyword,'1WILC') then do; call symput('ProbType','1Wilcoxon'); WlcxProb = 1; NumGrps = 1; rXfull = 1; if not(0 < count < 4) then go to WlcxErr; Wp1 = value{1}; Wp2 = value{2}; Wp3 = (Wp2 + Wp1**2)/2; Wp4 = value{3}; if (Wp2 = .) or (Wp4 = .) then do; Wp2 = .; Wp3 = .; Wp4 = .; p1OnlyIn = 1; WlcxCstm = 0; NomParnt = "NORMAL"; end; if (Wp2 ne .) and (Wp2 < Wp1) then do; put // "ERROR: p2 must be greater than p1."; stop; end; if (Wp2 ne .) and (Wp4 ne .) and (Wp2 le Wp4) then do; put // "ERROR: p4 must be less than p2."; stop; end; Wp1Custm = Wp1; Wp2Custm = Wp2; Wp3Custm = Wp3; Wp4Custm = Wp4; go to NextSpec; end; /* 1WIL do */ if index(keyword,'2WILC') or index(keyword,'MANN') then do; call symput('ProbType','2Wilcoxon'); WlcxProb = 2; NumGrps = 2; rXfull = 2; if not(0 < count < 4) then go to WlcxErr; Wp1 = value{1}; Wp2 = value{2}; Wp3 = value{3}; if (Wp2 = .) or (Wp3 = .) then do; Wp2 = .; Wp3 = .; p1OnlyIn = 1; WlcxCstm = 0; NomParnt = "NORMAL"; end; if (Wp2 ne .) and (Wp2 ge Wp1) then do; put // "ERROR: p2 must be less than p1."; stop; end; if (Wp3 ne .) and (Wp3 ge Wp1) then do; put // "ERROR: p3 must be less than p1."; stop; end; Wp1Custm = Wp1; Wp2Custm = Wp2; Wp3Custm = Wp3; Wp4Custm = Wp4; go to NextSpec; end; /* 2WIL or MANN do */ /* ++++ End UPw09_ProbPars1 segment ++++ */ /* ++++ Begin UPw10_ProbPars2 segment v. 2002.07.01 ++++ */ OrderCat: /* Ordered categorical problems via Wilcoxon testing, as per Lesaffre, et. al (1993); see note above. For 2-group problem, user specifies 2Wilcoxon # # # # # # # # # # # # # # # # where each row defines a probability dist'n, in this case with NumCat = 8 categories. For the 1-group problem, the category values are assumed to be interval in scale and are taken to be Y = {1, 2, ..., NumCat}. The user, however, may specify a different LOWer and UPPer value, but the sequence is assumed to be an arithmetic progression. The user must specify a NULL statement, where 0 < NullValu < NumCat, if LOW and UPP are not specified or LOWer < NullValu < UPP, if they are specified. */ if index(keyword,'1WILC') then do; NumGrps = 1; rXfull = 1; call symput("ProbType","1WilcOrdCat"); end; else if index(keyword,"2WILC") then do; call symput("ProbType","2WilcOrdCat"); NumGrps = 2; rXfull = 2; end; else go to WlcxErr; WlcxProb = NumGrps; WlcxOrCt = 1; SDvalid = 0; if NumGrps = 1 then do; call symput ("DistList", "DistLst1 NullValu"); call symput ("OthScen", "* DistLst1='' * NullValu='NULL:'"); end; if NumGrps = 2 then do; call symput ("DistList", "DistLst1 DistLst2"); call symput ("OthScen", "* DistLst1=''*DistLst2=''"); end; do iGrp = 1 to NumGrps; link EchoInpt; input DistList{iGrp} $ 1-78 @1 @; distlist{iGrp} = "{" || trim(distlist{iGrp}) || "}"; WlcxOrCt = 2; link GetValus; WlcxOrCt = 1; if iGrp = 1 then NumCat = count; else if (NumCat ne count) then do; put // "ERROR: Number of categories is not identical."; stop; end; do iCat = 1 to NumCat; if not(0 le value{iCat} le 1) then do; put "ERROR: Improper probability value."; stop; end; if NumGrps = 1 then probY{iCat} = value{iCat}; else if (NumGrps = 2) and (iGrp = 1) then probX{iCat} = value{iCat}; else probY{iCat} = value{iCat}; end; end; go to NextSpec; WlcxErr: put // "ERROR: Wilcoxon specification is bad."; stop; end; /* either WIL or MAN do */ /* parse the 2PI or 2PROPORTIONS statement */ if keyword = '2PI' or index(keyword, '2PRO') then do; call symput('TablType','Specl2pi'); _2piProb = 1; piProb = 1; ValidPrb = 1; NumGrps = 2; rXfull = 2; SDIn = 1; /* SD is fixed by the pis (computed later) and weights */ NumSD = 1; SDV{1} = 1; go to NextSpec; end; /* 2PI */ /* parse the PI or PROPORTIONS statement */ if keyword = 'PI' or index(keyword, 'PRO') then do; call symput('ProbType','pi'); call symput('TablType','GnrlPow'); piProb = 1; parent = "Bernoulli"; ValidPrb = 1; link GetValus; NumGrps = count; rXfull = NumGrps; if NumGrps = 1 then do; call symput ('TablTyp3','Pi1Specl'); end; if NumGrps = 2 then do; /* do all 4 tests, unless METHOD statement later changes this */ DoChi2pi = 1; DoXctUnC = 1; DoFshXCt = 1; DoLR = 1; end; do i = 1 to NumGrps; pi{i} = value{i}; if (pi{i} lt 0) or (pi{i} gt 1) then do; put // "ERROR: Improper pi: " pi{i}; stop; end; if (pi{i} = 0) then do; pi{i} = .001; put // "WARNING: pi = 0 is not permissable."; put " It has been automatically reset to pi = 0.001"; end; if (pi{i} = 1) then do; pi{i} = .999; put // "WARNING: pi = 0 is not permissable."; put " It has been automatically reset to pi = 0.999"; end; mu{i} = log(pi{i}/(1 - pi{i})); /* create logits */ end; /* i = 1 to NumGrps; */ SDIn = 1; /* SD is fixed by the pis and weights */ NumSD = 1; SDV{1} = 1; go to NextSpec; end; /* PI */ /* parse the 2WayContTable and GoodnessOfFit statements */ /* Handles (1) test of homogeneity of discrete distributions ("2-way independence") for a NumGrps x NumCat contingency table and (2) goodness of fit test of whether one group's prob distribution for a discrete Y comes from a specified null distribution. Both Pearson and likelihood ratio tests are considered. Examples: (1) 3 groups (rows) by 4 categories (columns), giving chi-square tests with 6 degrees of freedom: 2WayContTable .2 .3 .4 .1 > .2 .3 .1 .4 . Values must sum to 1.0. > .6 .2 .1 .1 weight .2 .4 .4 (2) Goodness of fit test on a 4-category probability distribution, giving chi-square tests with 3 degrees of freedom: GoodnessOfFit .6 .2 .1 .1 Null .3 .2 .2 .2 */ if index(keyword, '2WAY') or index(keyword, 'GOOD') then do; ValidPrb = 1; if index(keyword, '2WAY') then do; _2WCTPrb = 1; * test independence; call symput('ProbType','2WayCT'); end; if index(keyword, 'GOOD') then do; GdFtProb = 1; * test goodness of fit; call symput('ProbType','GoodOFit'); end; call symput('TablType','GnrlPow'); DoChiSq = 1; SDIn = 1; NumSD = 1; SDV{1} = 1; link GetValus; if count = 0 and _2WCTPrb then do; put // "ERROR: First group's probabilities must be on 2WayContTable line."; stop; end; if count = 0 and GdFtProb then do; put // "ERROR: Probabilities must be on GoodnessOfFit line."; stop; end; NumGrps = 1; NumCat = count; iGrp = 1; link SetProbs; if GdFtProb then go to NextSpec; /* Get probabilities for more groups */ MoreRows: link EchoInpt; input keyword @; if keyword = ">" then do; NumGrps = NumGrps + 1; link GetValus; if count ne NumCat then do; put // "ERROR: Number of categories is not identical."; stop; end; iGrp = NumGrps; link SetProbs; go to MoreRows; end; if keyword ne ">" then do; if NumGrps = 1 then do; put // "ERROR: Only one group of probabilities given."; stop; end; input @1 @; NoEcho = 1; go to NextSpec; end; end; /* 2WayCT and GoodOFit */ /* Parse the MU(P) statement. Handles designs in which subjects' data, Y, are binomial based on TRIALS trials per subject (fixed) and success probability, P, fixed to each subject but random across subjects. P is taken to be a standard beta(p_,q_) random variable with mean mu(P) and standard deviation SD(P). Thus, the user first specifies a set of J mu(P) values and also sets one or more values for SD(P) and Trials. Input is checked to assure that mu(P) = muP, 0 < MuP < 1 SD(P) = SDP, 0 < SDP < sqrt(1/12), with sqrt(1/12) limit imposed because SD(P) > sqrt(1/12) gives beta distributions for P that are not unimodal, a condition that makes little sense in building scenarios for power analyses. Power is determined as if the analysis will use ordinary Normal-theory methods on P_hat = Y/Trials. E(P_hat) = mu(P), but SD(P_hat) must be computed. For each mu(P) and pairing of SD(P) and Trials specifications, SD(P_hat) is determined using macro SDBetBnl, which contains the relevant references to distribution theory. So that all beta densities are unimodal, UnifyPow will abort runs if it finds that [mu(P), SD(P)] give (p_ - 1)*(q_ - 1) le 0. The common SD(P_hat) comes from the weighted average of Var(p_hat) across the J groups. WILCOXON can also be specified for the MU(P) problem. In this case, the Lehmann-Hettsmansperger parameters (p1 = Wp1, p2, p3, p4) are determined directly from exact beta-binomial distributions and the NULL statement. */ if keyword = 'MU(P)' then do; ValidPrb = 1; call symput('ProbType','muP'); call symput('TablType','FmuP_Pow'); muP_Prob = 1; link GetValus; NumGrps = count; rXfull = NumGrps; if NumGrps = 1 then NullValu = .; /* Must specify NullValu for this */ if NumGrps le 2 then call symput('TablType','tmuP_Pow'); do i = 1 to NumGrps; mu{i} = value{i}; if not(0 < mu{i} < 1) then do; put // "ERROR: Improper mu(P): " mu{i}; stop; end; if (mu{i} = 0) then do; /* 200 */ mu{i} = .001; put // "WARNING: mu(P) = 0 is not permissable."; put " It has been automatically reset to mu(P) = 0.001"; end; if (mu{i} = 1) then do; mu{i} = .999; put // "WARNING: mu(P) = 0 is not permissable."; put " It has been automatically reset to mu(P) = 0.999"; end; end; go to NextSpec; end; /* parse the MCNEMAR statement */ if keyword = 'MCNEMAR' then do; ValidPrb = 1; call symput('ProbType','McNemar'); call symput('TablType','GnrlPow'); McNrProb = 1; parent = "CorrBern"; link GetValus; if count ne 2 then do; put // "ERROR: Need exactly 2 pi values for McNemar's test."; stop; end; pi12 = value{1}; pi21 = value{2}; if not(0 < pi12 < 1) or not(0 < pi21 < 1) then do; put // "ERROR: A pi value is not between 0 and 1."; stop; end; if pi12 + pi21 > 1 then do; put // "ERROR: The sum of the pi values exceeds 1.0."; stop; end; SDIn = 1; /* SD is fixed by the the pis and weights */ go to NextSpec; end; /* parse the RHO statement */ if index(keyword, 'RHO') or index(keyword, 'CORR')=1 then do; ValidPrb = 1; call symput('ProbType','rho'); call symput('TablType','GnrlPow'); rhoProb = 1; link GetValus; NumGrps = count; if NumGrps > 1 then EffcTitl = "Comparing correlations (r-to-Z)"; do i = 1 to NumGrps; if i = 1 then rho1 = value{i}; if (value{i}**2 ge 1) then do; put "Improper rho: " value{i}; stop; end; else mu{i} = .5*log((1 + value{i})/(1 - value{i})); end; /* For RHO problem, the effective SD is fixed at SD = 1.0 */ SDIn = 1; NumSD = 1; SDV{1} = 1; go to NextSpec; end; /* RHO parser */ /* parse the EXEMPLARY statement */ if index(keyword, 'EXEMPLARY') then do; ValidPrb = 1; ExemProb = 1; call symput('ProbType','exemplary'); input keyword $; keyword = upcase(keyword); if index(keyword, 'SSH') then do; SSH_Prob = 1; call symput('TablType','FPow'); call symput('SDType','Standard Deviation'); end; else if index(keyword, 'CHI') then do; Chi2Prob = 1; call symput('TablType','GnrlPow'); SDIn = 1; NumSD = 1; SDV{1} = 1; end; else do; put // "ERROR: EXEMPLARY type must be SSH or CHI**2."; stop; end; go to NextSpec; end; /* EXEMPLARY parser */ /* parse the R**2 or PartialCorr statement */ if keyword = 'R**2' or index(keyword, 'PARTIAL') then do; ValidPrb = 1; call symput('ProbType','R2'); call symput('TablType','GnrlPow'); RSqProb = 1; link GetValus; /* PARTIALCORR problem is just transformed to R**2 problem, unless NullValu ne 0. In that case, it is switched to an r-to-Z problem with N-3 replaced by N-rXfull-1. Justification can be found in Section 27.22 of Kendall and Stuart, The Advanced Theory of Statistics, Vol. 2, 4th Ed, 1979. The main result is due to RA Fisher in "The distribution of the partial correlation coefficient," Metron, 3: 329, 1924. */ if index(keyword, 'PARTIAL') then do; PCorrPrb = 1; if count ne 1 then do; put // "ERROR: You must give only 1 partial correlation value."; stop; end; if (value{1}**2 ge 1) then do; put // "ERROR: Partial correlation value must be between -1 and +1."; stop; end; R2redc = 0; R2full = value{1}**2; mu{1} = .5*log((1 + value{1})/(1 - value{1})); go to FixSD1; end; /* keyword = PARTIAL */ NumRSq = count; if (NumRSq ne 1) and (NumRSq ne 2) then do; put // "ERROR: You must give 1 or 2 R**2 values."; stop; end; if NumRSq = 1 then do; R2redc = 0; R2full = value{1}; end; if NumRSq = 2 then do; R2full = max(value{1}, value{2}); R2redc = min(value{1}, value{2}); end; if not(0 < R2full < 1) or not(0 le R2redc < 1) then do; put // "ERROR: R**2 values must be non-negative and less than 1."; stop; end; /* For R**2, the effective SD is fixed at SD = 1.0 */ FixSD1: SDIn = 1; NumSD = 1; SDV{1} = 1; go to NextSpec; end; /* R**2 parsing */ /* parse the 1betaOLS statement */ if index(keyword,'1BETA') then do; ValidPrb = 1; call symput('TablType','1betaOLS'); betLSPrb = 1; link GetValus; NumBetas = count; do i = 1 to NumBetas; BetaWtV{i} = value{i}; end; go to NextSpec; end; /* parse the Sample statement */ if index(keyword, 'OBSERV') or index(keyword, 'SAMPLE') then do; call symput('TablType','PowCI'); ValidPrb = 1; SampProb = 1; SDIn = 1; NumSD = 1; SDV{1} = 1; call symput('ProbType','sample'); input keyword @; keyword = upcase(keyword); if keyword = 'F' then do; SamFProb = 1; input; end; else if index(keyword, 'CHI') then do; SamCProb = 1; input; end; else if index(keyword,'T') then do; SamTProb = 1; input; end; else if index(keyword, 'MEAN') then do; put "SORRY: This module not built yet."; stop; end; /* Note to myself. Here is where I might build the entry to let the user input sample statistics from a G group problem. The syntax might look something this: sample means # # # # sample SDs # # # # . or "sample pooledSD #" sample weights # # # # . or "sample n" sample Ntotal # . not needed if using "sample n" This would allow UnifyPow to estimate the noncentrality for any effect of interest. Most code to do this should already be in here, because the computations mirror those to compute the noncentrality. A similar thing could be done for other problems, e.g. proportions. */ else if index(keyword, 'PROP') then do; put "SORRY: This module not built yet."; stop; end; else do; put // "ERROR: SAMPLE statistic must be T, F, or CHI**2."; stop; end; /* set default confidence levels */ NumCnfLv = 2; ConfLevV{1} = .50; ConfLevV{2} = .75; go to NextSpec; end; /* SAMPLE parser */ /* ++++ End UPw10_ProbPars2 segment ++++ */ /* ++++ Begin UPw11_MVmuPars segment v. 2002.07.01 ++++ */ /* parse MVmu or MVmean statement */ /* Handles the multivariate general linear model quite generally. Theory is given in O'Brien and Shieh (manuscript available on request). This section inputs a NumGrps x NumMeas matrix of means. Example: /# 5 groups: 1 Control + 4 Treatments (A01, A02, B01, B02) 4 measurements (times): Hr00, Hr01, Hr02, Hr03 The complete MVmu specification defines B, the matrix of the group centroids. #/ MVmeans 40 43 38 33 . control > 40 40 35 30 . Treatment A, subtype 1 > 40 40 33 30 . Treatment A, subtype 2 > 40 38 35 30 . Treatment B, subtype 1 > 40 37 33 28 . Treatment B, subtype 2 /# The weights, Ntotal, and alpha statements work as usual. #/ weights 2 1 1 1 1 . balanced w.r.t. to treatments NTotal 36 72 144 alpha .01 .05 /# You must specify a basis set of SDs, exactly one SD for each measure. And you must specify the upper triangle of the within-groups correlation matrix. UnifyPow will compute the basis covariance matrix. #/ sd 12 11 10 8 . SDs of Hr00, Hr01, Hr02, Hr03 corr .7 .4 .2 . correlations of Hr00 with Hrs 01, 02, 03 > .7 .3 . correlations of Hr01 with Hrs 02, 03 > .8 . correlation of Hr02 with Hr03 #/ You can do sentivity analyses on the covariance matrix using the SDmult statement. Here, we will use the basis SD's (SDmult = 1.0) and SDs that are all 20% larger (SDmult = 1.20). #/ SDmult 1.0 1.2 . Assess effect of 20% larger SDs /# By default, UnifyPow computes the power for Pillai's statistic. The "method" statement overrides this and can specify method Wilks method Hotelling method Pillai method Wilks Hotelling . (or any other pair) method all #/ methods Pillai . not necessary, since this is default /# By default, UnifyPow assumes that the measures are distinct variables, *NOT* repeated measures on the same variable. This corresponds to the default specification measures distinct UnifyPow will automatically do a one-way MANOVA, unless a "NoOverall" statement is given. If the measures are repeated measures on a single factor, use measures repeated UnifyPow will then automatically perform power analyses for the overall Group and Repeated (here, Time) main effects and the overall Group x Repeated interaction. Here again, "NoOverall" will cause this step to be skipped. #/ measures repeated . Time factor /# UnifyPow's "contrasts" capability gives it great flexibility in handling MANOVA and repeated measures designs. You can obtain power for any hypothesis of the form Ho: CBA = Theta, where C is a full-rank contrast matrix that compares groups. B is the "MVmu" matrix of cell means; see above. A is a full-rank contrast matrix that compares measures. Theta is usually a zero matrix (default). Each special Ho: CBA = Theta contrast begins with a "title of contrast" statement and is specified using the Group (C), Measures (A), and Theta statements. Both Group and Measures contrast matrices must be given for the first contrast, but Theta is assumed to be a zero matrix. Thereafter, if any of these is absent, the one used previously is re-used. #/ contrasts /# The next contrast sharpens the overall interaction test by ignoring (averaging over) subtype, and only looks at linear and quadratic trends over time. #/ "Trtmnt x Time(L+Q)" Group 2 -1 -1 0 0 . control vs. A > 2 0 0 -1 -1 . control vs. B Measures -3 -1 1 3 . time(linear) > 1 -1 -1 1 . time(quadratic) /# The next contrast uses the previous group contrast, but now time only looks at the average of Hours 02 and 03 versus Hour 00. #/ "Trtmnt x Time(00 v 02+03)" Measures 2 0 -1 -1 . Hour00 vs ave(Hour02 & Hour03) /# The next contrasts test how the time(L+Q) trend varies over subtypes, separately for treatments A and B. The alpha level is increased substantially to better balance Type I and Type II error rates. #/ alpha .10 .20 . This replaces previous alpha statement. "Subtype x time(L+Q)" Group 0 1 -1 0 0 . subtype within A Measures -3 -1 1 3 . time(linear) > 1 -1 -1 1 . time(quadratic) "Subtype x time(L+Q)" Group 0 0 0 1 -1 . subtype within B */ /* begin parse MVmu */ if index(keyword, 'MVM') then do; ValidPrb = 1; MVmuProb = 1; call symput('ProbType','MVmu'); call symput('TablType','MVmu'); %matInput (MVmu, NumGrps, NumMeas); go to NextSpec; end; /* end parse MVmu */ if not(ValidPrb) then do; put // "ERROR: Invalid problem statement."; stop; end; /* ++++ End UPw11_MVmuPars segment ++++ */ /* ++++ Begin UPw12_OthrPars segment v. 2002.07.01 +++++ */ NextSpec: CurInput = "OthrStmt"; link EchoInpt; input keyword $ @; keyword = upcase(keyword); NxtSpec2:; /* Parse the SCENARIO statement. This can happen either before of after parsing of problem statement. */ link ParsScen; /* parse the WILCOXON option for mu and mu(P) problems */ /* For the one or two group MU problem, this initiates a Normal approximation to set p1 = Wp1 (defined above). For MU(P), this initiates exact theory to set p1, p2, p3, p4 unless NumGrps = 2 and NullValu NE 0. In that case, only the Normal approximation is used. */ if index(keyword,'WILCOX') or index(keyword,'MANN') then do; if WlcxProb = 0 then input; if ((muProb=0) and (muP_Prob=0) and (PaMuProb=0)) or (NumGrps > 2) then do; put // "ERROR: Rank test options not allowed here."; stop; end; if ANCVprob then do; put // "ERROR: Rank test options not allowed with covariates."; put " Only parametric tests will be evaluated."; go to NextSpec; end; WlcxProb = NumGrps; if MuProb then call symput('TablType','WlcxnPow'); if muP_Prob then call symput('TablType','WlcxMuPpow'); if PaMuProb then call symput('TablType','WlcxPaMuPow'); if NoNotes = 0 then call symput ('TblMmnts','yes'); go to NextSpec; end; /* parse the PROPORTION_1 statement */ if (index(keyword, 'PROP') or index(keyword, 'PI')) and (index(keyword, '1') > 3) then do; link GetValus; NumProp1 = count; Prop1in = 1; do i = 1 to NumProp1; pi1{i} = value{i}; if not (0 < pi1{i} < 1) then do; put // "ERROR: Improper proportion: " pi1{i}; stop; end; end; /* i = 1 to NumProp1; */ go to NextSpec; end; /* proportion_1 */ /* parse the PROPORTION_2 statement */ if (index(keyword, 'PROP') or index(keyword, 'PI')) and (index(keyword, '2') > 3) then do; if OddsProb + PrRtProb ne 0 then do; put // "ERROR: Either an OddsRatio or ProbRatio statement already used."; stop; end; call symput('pi2SpcVr','pi_2'); call symput('pi2SpcLb','Group 2 Prob'); Prop2in = 1; link GetValus; NumProp2 = count; do i = 1 to NumProp2; pi2{i} = value{i}; if not (0 < pi2{i} < 1) then do; put // "ERROR: Improper proportion: " pi2{i}; stop; end; end; /* i = 1 to NumProp2; */ go to NextSpec; end; /* proportion_2 */ /* parse the ODDSRATIO, PROBRATIO, RELRISK, RISKRATIO statements */ if index(keyword, 'RATIO') or index(keyword,'RISK') then do; if Prop2in ne 0 then do; put // "ERROR: PROPORTION_2 statement already used."; stop; end; parent = "Bernoulli"; NullOrig = 1; /* Do only LR test, unless METHOD statement later changes this or NullOrig becomes different from 1.0. */ DoChi2pi = 0; DoXctUnC = 0; DoFshXCt = 0; DoLR = 1; if index(keyword, 'ODDS') then do; call symput('pi2SpcVr','OddsRato'); call symput('pi2SpcLb','Odds Ratio'); OddsProb = 1; end; else if index(keyword, 'PRO') or index(keyword, 'RISK') then do; call symput('pi2SpcVr','ProbRato'); if index(keyword, 'PRO') then call symput('pi2SpcLb','Prob Ratio'); if index(keyword, 'RIS') then call symput('pi2SpcLb','Relative Risk'); PrRtProb = 1; end; else put "ERROR: Keyword must be OddsRatio or ProbRatio|RelRisk|RiskRatio."; link GetValus; NumProp2 = count; do i = 1 to NumProp2; ratio{i} = value{i}; if ratio{i} le 0 then do; if OddsProb then put // "ERROR: Improper odds ratio: " ratio{i}; if PrRtProb then put // "ERROR: Improper ratio of probabilities: " ratio{i}; stop; end; end; /* i = 1 to NumProp2; */ go to NextSpec; end; /* OddsRatio or ProbRatio */ /* parse the METHOD option */ /* Erase all defaults. */ if index(keyword,'METHOD') then do; MethodIn = 1; DoLH = 0; DoChi2pi = 0; DoXctUnC = 0; DoFshXCt = 0; DoLR = 0; NextMeth: input keyword $ @; if keyword = "" then do; input; go to NextSpec; end; keyword = upcase(keyword); if WlcxProb ne 0 then do; if index(keyword,'NOETHER') then DoNoe = 1; else if index(keyword,'ARE') then DoARE = 1; else if index(keyword,'LEHMANN') then DoLH = 1; else if index(keyword,'ALL') then do; DoLH = 1; DoNoe = 1; DoARE = 1; end; else do; put // "ERROR: METHOD statement has wrong specification."; stop; end; if WlcxCstm and DoARE then put // "NOTE: The ARE method is not appropriate for this problem."; go to NextMeth; end; if (piProb ne 0) or _2piProb then do; if index(keyword,'CHI') then DoChi2pi = 1; else if index(keyword,'EXACTUN') then DoXctUnc = 1; else if index(keyword,'FISHER') then DoFshXCt = 1; else if index(keyword,'LR') then DoLR = 1; else if index(keyword,'ALL') then do; DoChi2pi = 1; DoXctUnC = 1; DoFshXCt = 1; DoLR = 1; end; else do; put // "ERROR: METHOD statement has wrong specification."; stop; end; go to NextMeth; end; do; put // "ERROR: METHOD statement is not appropriate here." / " Check its placement in specicifications list."; stop; end; end; /* METHOD parser */ /* parse the PARENT option for 1Wilcoxon and 2Wilcoxon/Mann-Whitney problem */ /* Default is to assume NomParnt = "NORMAL" */ if index(keyword,'PARENT') then do; input keyword $; keyword = upcase(keyword); if index(keyword,'NORMAL') then NomParnt = "NORMAL"; else if index(keyword,'LOGISTIC') then NomParnt = "LOGISTIC"; else if index(keyword,'LAPLACE') then NomParnt = "LAPLACE"; else do; put // "ERROR: Parent must be Normal (default), Logistic, or Laplace."; stop; end; if WlcxProb = 0 then do; WlcxProb = NumGrps; keyword = 'WILCOXON'; go to NxtSpec2; end; go to NextSpec; end; /* parse the LIMITS statement */ if keyword = 'LIMITS' then do; LimitsIn = 1; link GetValus; if count ne 2 then do; put // "ERROR: Not exactly 2 LIMITS."; stop; end; UppLimit = max(value{1}, value{2}); LowLimit = min(value{1}, value{2}); go to NextSpec; end; /* parse the Nexemplary statement for SSH_Prob and Chi2Prob */ if index(keyword,'NEXE') then do; NexIn = 1; if SSH_Prob = 0 and Chi2Prob = 0 then do; put // "ERROR: EXEMPLARY statement is missing."; stop; end; link GetValus; if count ne 1 then do; put // "ERROR: There should only be one Nexemplary value."; stop; end; Nxmplr = value{1}; go to NextSpec; end; /* Nexemplary parsing */ /* parse the SampleN statement for the OBSERVED problem */ if index(keyword,'SAMPLEN') then do; SamNIn = 1; if SampProb = 0 then do; put // "ERROR: OBSERVED statement is missing."; stop; end; link GetValus; if count ne 1 then do; put // "ERROR: There should only be one SampleN value."; stop; end; SampleN = value{1}; go to NextSpec; end; /* SampleN parsing */ /* parse NumParms statement for R**2, PartialCorr, 1betaOLS, SSH SamF, and SamT problems */ if index(keyword,'NUMPARM') then do; NmParmIn = 1; if (RSqProb+PCorrPrb+betLSPrb+SSH_Prob+SamFProb+SamTProb = 0) then do; put // "ERROR: NumParm statement not defined for this problem."; stop; end; link GetValus; if PCorrPrb then do; if count ne 1 then do; put // "ERROR: Only one NumParm value allowed."; stop; end; rXfull = value{1}; rXredc = rXfull - 1; if rXfull < 2 then do; put // "ERROR: NumParm value must be 2 or more."; stop; end; if rXfull = 2 then put // "WARNING: NumParm value is normally 3 or more."; go to NextSpec; end; if RSqProb then do; if NumRSq ne count then do; put // "ERROR: Wrong number of NumParm values."; stop; end; if NumRSq = 1 then do; rXfull = value{1}; rXredc = 1; end; if NumRSq = 2 then do; rXfull = max(value{1}, value{2}); rXredc = min(value{1}, value{2}); end; if (rXfull = rXredc) or (rXredc < 1) or (rXfull ne round(rXfull)) or (rXredc ne round(rXredc)) then do; put // "ERROR: Improper NumParm value."; stop; end; go to NextSpec; end; if betLSPrb or SSH_Prob or SamFProb or SamTProb then do; if count ne 1 then do; put // "ERROR: Specify exactly 1 NumParm value."; stop; end; if (value{1} < 1) or (value{1} ne round(value{1})) then do; put // "ERROR: Improper NumParm value"; stop; end; rXfull = value{1}; go to NextSpec; end; end; /* NumParms parsing */ /* parse the SDX statement (SD of X variable in 1betaOLS Problem) */ if keyword = 'SDX' or keyword = 'SD_X' then do; if betLSPrb ne 1 then do; put // "SD_X statement must follow 1betaOLS statement."; stop; end; link GetValus; SDxIn = 1; NumSDx = count; do i = 1 to NumSDx; SDxV{i} = value{i}; if not(SDxV{i} > 0) then do; put // "ERROR: SDX value is not positive."; stop; end; end; go to NextSpec; end; /* parse the Tolerance statement (1 - R^2 of Xj predicted from other Xs) */ if index(keyword,'TOLER') then do; if betLSPrb ne 1 then do; put // "Tolerance statement must follow 1betaOLS statement."; stop; end; link GetValus; TolIn = 1; NumTol = count; do i = 1 to NumTol; TolV{i} = value{i}; end; go to NextSpec; end; /* parse the SDMULT and CVMULT statements */ if index(keyword, 'MULT') then do; if index(keyword, 'SD') then do; if not(muProb or MVmuProb or PaMuProb or SamFProb or SamTProb) or LogNProb then do; put // "ERROR: SDMULT statement does not fit this problem."; stop; end; go to SetSDmlt; end; /* keyword = SDMULT */ if index(keyword, 'CV') then do; if not(LogNprob) then do; put // "ERROR: CVMULT statement only fits logNormal problems."; stop; end; go to SetSDmlt; end; /* keyword = CVMULT */ SetSDmlt: link GetValus; NumSDmlt = count; do i = 1 to NumSDmlt; SDMultV{i} = value{i}; if SDMultV{i} le 0 then do; put // "ERROR: Multiplier must be positive."; stop; end; end; go to NextSpec; end; /* keyword = SDMULT or CVMULT */ /* Parse the STATISTIC statement for MVmu problem. Default is to only do for the Pillai. */ if index(keyword,'STATISTIC') and MVmuProb then do; NextSTAT: input keyword $ @; if keyword = "" then go to NextSpec; keyword = upcase(keyword); if keyword = 'NOWILKS' then DoWilks = 0; else if keyword = 'WILKS' then DoWilks = 1; else if index(keyword,'NOHOTEL') then DoHotel = 0; else if index(keyword,'HOTEL') then DoHotel = 1; else if index(keyword,'NOPILL') then DoPillai = 0; else if index(keyword,'PILL') then DoPillai = 1; else if index(keyword,'ALL') then do; DoWilks = 1; DoHotel = 1; DoPillai = 1; end; else do; put // "ERROR: STATISTIC statement has wrong specification."; stop; end; go to NextSTAT; end; /* STATISTIC parser */ /* parse the WEIGHT statement */ /* Note: Weights are relative. They are summed and converted to fractions. Thus use WEIGHT .333 .666 or, equivalently, WEIGHT 1 2 This is not the same as WEIGHT .333 .667 which may cause problems, especially when finding NTotal given a Power specification. */ if index(keyword,'WEIGHT') then do; link GetValus; wIn = 1; NumWght = count; if NumGrps ne NumWght then do; put // "ERROR: Improper number of WEIGHTS specified."; stop; end; SumWts = 0; do i = 1 to NumWght; w{i} = value{i}; SumWts = SumWts + w{i}; end; do i = 1 to NumWght; w{i} = w{i}/SumWts; wDesign{i} = w{i}; end; go to NextSpec; end; /* parse the TRIALS (or ITEMS) statement for the mu(P) problem */ if index(keyword,'TRIALS') or index(keyword,'ITEMS') then do; if index(keyword,'TRIALS') then call symput('TrlName','Trials'); if index(keyword,'ITEMS') then call symput('TrlName','Items'); if muP_Prob = 0 then do; put // "ERROR: " keyword $ "statement must be part of mu(P) problem."; stop; end; link GetValus; TrialsIn = 1; NumTrls = count; do i = 1 to NumTrls; TrialsV{i} = value{i}; if not(TrialsV{i} > 0) or (floor(TrialsV{i}) ne TrialsV{i}) then do; put // "ERROR: Improper " keyword $ "specification: " TrialsV{i}; stop; end; end; go to NextSpec; end; /* parse the ALPHA and TAILS statements */ link CkAlfTls; if NewAlfTl then go to NextSpec; /* parse BALANCE statement for binomial problem */ if index(keyword,'BALANCE') then do; BalTails = 1; go to NextSpec; end; /* parse the standard deviations statements: SD or SIGMA or RelativeSD*/ if (keyword = 'SD') or (keyword = 'SIGMA') or index(keyword,'RELATIVE') then do; if SDvalid = 0 then do; put // "WARNING: SD (SIGMA) statement is invalid for this problem."; put " It is being ignored."; go to NextSpec; end; if p1OnlyIn then call symput ('SDType','Relative Std Dev'); if (muP_Prob = 1) then do; keyword = 'SD(P)'; putSDwrn = 1; go to getSDP; end; link GetValus; SDIn = 1; NumSD = count; if PriSDin then do; if NumPriSD ne NumSD then do; put // "ERROR: Number of SDs does not equal number of priors for SD"; stop; end; end; if WlcxCstm then do; NumSD = 1; SDV{1} = 1; put // "WARNING: " keyword "is not appropriate for this problem."; put " Statement is ignored."; end; if PaMuProb and (NumSD ne 2) then do; put // "ERROR: The PairedMu problem requires exactly"; put " 2 standard deviations."; stop; end; if MVmuProb and (NumSD ne NumMeas) then do; put // "ERROR: The MVmu problem requires the same number of base SDs"; put " as there are measures, in this case: " NumMeas 2.0; stop; end; do i = 1 to NumSD; SDV{i} = value{i}; if not(SDV{i} > 0) then do; put // "ERROR: A standard deviation is not positive."; stop; end; end; if PaMuProb and not(ScenarIn) then do; scenario = trim(ProbStmt)||' & '||compress(keyword); do i = 1 to 2; scenario = trim(scenario)||' '||compress(SDV{i}); end; end; go to NextSpec; end; /* SD parser */ /* parse the priors_SD or priors_CV statements */ if index(keyword,'PRIOR') and (index(keyword,'SD') or index(keyword,'CV')) then do; call symput ('TablSDwt','yes'); link GetValus; PriSDIn = 1; NumPriSD = count; if SDin and (NumPriSD ne max(NumSD,NumCV)) then do; put // "ERROR: Number of priors does not equal number of SDs."; stop; end; link PriAdjst; do i = 1 to NumPriSD; PriSDV{i} = value{i}; end; go to NextSpec; end; /* priors_SD/cv parser */ /* parse the Coefficient of Variation statement: CV */ if keyword = 'CV' then do; CVIn = 1; SDIn = 1; if not(LogNProb) then do; put // "ERROR: CV statement is not valid for this problem."; stop; end; link GetValus; NumCV = count; if PriSDin then do; if NumPriSD ne NumCV then do; put // "ERROR: Number of CVs does not equal number of priors for SD"; stop; end; end; if PaMuProb and (NumCV ne 2) then do; put // "ERROR: The PairedMu problem requires exactly"; put " 2 coefficients of variation (CV)."; stop; end; do iCV = 1 to NumCV; CVV{iCV} = value{iCV}; if (CVV{iCV} le 0) then do; put // "ERROR: A coefficient of variation is not positive."; stop; end; end; if PaMuProb and not(ScenarIn) then do; scenario = trim(ProbStmt)||' & CV:' ; do i = 1 to 2; scenario = trim(scenario)||' '||compress(CVV{i}); end; end; go to NextSpec; end; /* CV parser */ /* parse the SD(P) statement */ if keyword = 'SD(P)' then do; getSDP: link GetValus; if putSDwrn = 1 then put // "NOTE: SD (or SIGMA) values assumed to be for pi" / " as part of mu(P) problem." //; SDPIn = 1; NumSDP = count; do i = 1 to NumSDP; SDPV{i} = value{i}; if (SDPV{i} > sqrt(1/12)) then do; put // "ERROR: SD(P) = " SDPV{i} "exceeds sqrt(1/12), resulting in "; put " beta distribution for P that is not unimodal."; stop; end; end; go to NextSpec; end; /* SD(P) parser */ /* Define "clean" ANCOVA problems as those in which the covariates are unrelated to the "treatments." This will be the case, for example, when the covariates are taken at baseline and treatment is randomized. means mu1 ... SD SD1 ... covariates #C . give one value for number of covariates added to cell means model PartialMultR #R1 #R2 #R3 The PartialMultR statement sets forth one or more possible partial multiple correlations (#R) of the #C covariates with the outcome after adjusting for the dummy predictors that model the adjusted group means. For simple ANCOVA, #C = 1, and #R is the pooled within-group correlation of the outcome measure with the covariate. One can use PartialCorr #R1 #R2 #R3 This will reduce the DFe term by #c and reduce the SD by 100*[1 - sqrt(1 - R**2)]%. Under the "clean" ANCOVA model, there is no change in the primary noncentrality for the mean differences. Alternatively, one can just use means mu1 ... SD SD1 ... covariates #C SDreduction #D1 #D2 #D3 Here, #D is the reduction in the SD achieved by adding the #C covariates to the model. */ /* parse the Covariates statement */ if index(keyword,'COVARIATE') then do; if WlcxProb then do; put // "NOTE: Covariate statement not associated with Wilcoxon-Mann-Whitney problem."; put " Evaluation will proceed without covariate adjustments." /; input; ANCVprob = -1; go to NextSpec; end; if not(muProb) then do; put // "NOTE: Covariate statement not associated with this problem."; put " Evaluation will proceed without covariate adjustments." /; input; ANCVprob = -1; go to NextSpec; end; if muProb then do; call symput('TablType', 'muANCOVA'); ANCVProb = 1; input NmCovars; if (round(NmCovars) ne NmCovars) or NmCovars < 1 then do; put // "ERROR: incorrect specification for number of covariates."; stop; end; rXfull = rXfull + NmCovars; go to NextSpec; end; end; /* Covariates */ /* parse the PartialMultR/PartialCorr/Corr statement for COVARIATES problem */ if index(keyword,'CORR') and PaMuProb then go to PaMuCorr; if index(keyword,'CORR') and MVmuProb then go to MVmuCorr; if (index(keyword,'PARTIAL') or index(keyword,'CORR') or index(keyword,'MULTR')) and ANCVprob=0 then do; put // "ERROR: " keyword "statement requires previous COVARIATES statement"; put " It is being ignored." /; input; go to NextSpec; end; if (index(keyword,'PARTIAL') or index(keyword,'CORR') or index(keyword,'MULTR')) and abs(ANCVprob)=1 then do; if ANCVProb = -1 then do; put / "NOTE: " keyword "statement is being ignored." /; input; go to NextSpec; end; if SDredIn then do; put // "ERROR: SDreduction statement already used," / " so correlation statement not accepted." /; go to NextSpec; end; CvPCorIn = 1; call symput ('SDadjVar', 'CovPCorr'); if NmCovars = 1 then call symput ('SDadLabl', 'Covariate Partial Correlation'); if NmCovars > 1 then call symput ('SDadLabl', 'Covariates Partial Multiple R'); link GetValus; NumSDadj = count; do i = 1 to NumSDadj; if not(-1 < value{i} < 1) then do; put // "ERROR: Partial correlation is not between -1 and +1"; stop; end; SDadjV{i} = sqrt(1 - value{i}**2); end; go to NextSpec; end; /* Partial corrs for ANCOVA */ /* parse the SDreduction statement for COVARIATES problem */ if index(keyword,'REDUC') and ANCVprob=0 then do; put // "ERROR: " keyword "statement requires previous COVARIATES statement"; put " It is being ignored."; input; go to NextSpec; end; if index(keyword,'REDUC') and abs(ANCVprob)=1 then do; if ANCVProb = -1 then do; put / "NOTE: " keyword "statement is being ignored."; input; go to NextSpec; end; if CvPCorIn then do; put // "ERROR: partial correlation statement already used," / " so SDreduction statement not accepted."; go to NextSpec; stop; end; SDredIn = 1; call symput ('SDadjVar', 'SDredctn'); call symput ('SDadLabl', 'Proportion of SD Reduced by Covariates'); link GetValus; NumSDadj = count; do i = 1 to NumSDadj; SDadjV{i} = 1 - value{i}; if not(0 < SDadjV{i} =< 1) then do; put // "ERROR: SDreduction is not between 0 and 1"; stop; end; end; go to NextSpec; end; /* SDreduction for ANCOVA */ /* parse the CORR statement for PairedMU problem */ PaMuCorr: if index(keyword,'CORR') and PaMuProb then do; link GetValus; CorrIn = 1; NumCorr = count; do i = 1 to NumCorr; CorrV{i} = value{i}; if not(-1 < CorrV{i} < 1) then do; put // "ERROR: Correlation is not between -1 and 1."; stop; end; end; go to NextSpec; end; /* Parse the CORR statement for MVmu problem. This is the upper triangle only. */ MVmuCorr: if index(keyword,'CORR') and MVmuProb then do; CorrIn = 1; link getCorrM; go to NextSpec; end; /* parse the TotalN or NTotal (total sample size) statement */ if (index(keyword,'TOTAL') gt 0) then do; call symput('ResltVar','power'); link GetValus; NtotalIn = 1; NumNtotl = count; do i = 1 to NumNtotl; NtotalV{i} = value{i}; if not(NtotalV{i} > 0) or not(floor(NtotalV{i}) = NtotalV{i}) then do; put // "ERROR: Total sample size is not correct."; stop; end; end; NumPower = 1; NomPowrV{1} = 0; go to NextSpec; end; /* parse the Power statement */ if keyword = 'POWER' then do; call symput('ResltVar','NTotal'); link GetValus; PowerIn = 1; NumPower = count; MinPower = 1; do i = 1 to NumPower; NomPowrV{i} = value{i}; if NomPowrV{i} > .999 then call symput('PowerFmt','6.4'); MinPower = min(of MinPower NomPowrV{i}); if AlphaIn and not(MaxAlpha < NomPowrV{i} < 1) then do; put // "ERROR: Power must be between alpha and 1."; stop; end; end; NumNtotl = 1; NTotalV{1} = 1; go to NextSpec; end; /* parse the NULL statement */ if keyword = 'NULL' then do; NullIn = 1; if GdFtProb then do; /* Set of null probs for goodness of fit problem. */ link GetValus; if (count ne NumCat) then do; put // "ERROR: NULL has incorrect number of categories."; stop; end; sum = 0; do iCat = 1 to NumCat; sum = sum + value{iCat}; jp_null{1, iCat} = value{iCat}; end; if abs(sum-1.0) gt .001 then do; put // "ERROR: NULL probabilities do not sum to 1.0."; stop; end; go to NextSpec; end; /* GdFtProb */ input NullValu; if LogNProb or OddsProb or PrRtProb then do; if NullValu le 0 then do; put // "ERROR: NULL value must exceed 0.00 for this problem."; stop; end; NullOrig = NullValu; NullValu = log(NullValu); if OddsProb or PrRtProb then NullValu = -NullValu; end; /* reset default method for special 2pi problems */ if (MethodIn=0) and (NullValu ne 0) and _2piProb then do; DoChi2pi = 1; DoXctUnC = 0; DoFshXCt = 0; DoLR = 0; end; if (RSqProb and not(PCorrPrb)) or (rhoProb and (NumGrps > 1)) or (NumGrps > 2) then do; put // "ERROR: NULL is not defined for this problem."; stop; end; if ((rhoProb or PCorrPrb) and (NullValu**2 ge 1)) then do; put // "ERROR: Improper NULL for rho or PartialCorr problem."; stop; end; /* Handle PCorrPrb with nonzero NullValu by switching to Fisher's r-to-Z */ if (PCorrPrb and (NullValu ne 0)) then do; RSqProb = 0; rhoProb = 1; end; go to NextSpec; end; /* NULL parser */ /* parse the NoOverall statement */ if index(keyword,'NOOVERALL') then do; NoOveral = 1; input; go to NextSpec; end; /* parse the NoAdjustN statement */ if index(keyword,'NOADJUSTN') then do; DoAdjstN = 0; input; go to NextSpec; end; /* parse the ForceN statement for EXEMPLARY problems*/ if index(keyword,'FORCEN') then do; if not(ExemProb) then do; put // "ERROR: ForceN statement invalid for this problem"; stop; end; ForceN = 1; link getValus; if count ne 1 then do; put // "ERROR: Need exactly one value for ForceN"; stop; end; CoreNtot = value{1}; go to NextSpec; end; /* parse the NoTables statement */ if index(keyword,'NOTABLE') then do; call symput('DoTables','no'); input; go to NextSpec; end; /* parse the NoNotes statement */ if index(keyword,'NONOTE') then do; NoNotes = 1; call symput ('TblMmnts','no'); input; go to NextSpec; end; /* Parse the (optional) END statement. This will cause UnifyPow to run what is in and recycle back to NewProb to start new problem. */ if index(keyword,'END') then do; input; EndJstIn = 1; go to CheckIt; end; /* parse the ConfidenceLevel statement for Sample and CI problems */ if index(keyword,'CONF') or index(keyword,'LEVEL') then do; link GetValus; CnfLevIn = 1; NumCnfLv = count; do i = 1 to NumCnfLv; ConfLevV{i} = value{i}; if not(0 < ConfLevV{i} < 1) then do; put // "ERROR: Confidence level is not between 0 and 1."; stop; end; if (ConfLevV{i} < 0.50) then put // "WARNING: Do you really want a confidence level < 0.50?"; end; go to NextSpec; end; /* parse the CONTRAST statement */ if index(keyword,'CONTRAST') then do; input; CntrstIn = 1; call symput ('MethLabl', 'EffcTitl = "Effect"'); go to CheckIt; end; if index(keyword,'EFFECT') then do; input; EffectIn = 1; call symput ('MethLabl', 'EffcTitl = "Effect"'); go to CheckIt; end; else do; put // "ERROR: " keyword $ "statement not understood."; stop; end; EOFfound: if EndJstIn then go to AllDone; EndOFile = 1; if CurInput = "ProbStmt" or CurInput = "GtXmpyVl_1" or CurInput = "GtSamVl_1" then do; put // "ERROR: End of UnifyPow statements was reached too soon."; stop; end; /* ++++ End UPw12_OthrPars segment +++++ */ /* ++++ Begin UPw13_ChkParms segment v. 2002.07.01 ++++ */ if CurInput = "OthrStmt" then go to CheckIt; if CurInput = "GetMatrx" then go to EndGtMat; if CurInput = "GetCorrM" then go to ChkGtCor; if CurInput = "GtXmpyVl_2" then go to PowXmpry; if CurInput = "muContrast" then go to ContPowr; if CurInput = "MVmuCBA" then go to EOFinCBA; if OvAlDone or EffectIn then go to AllDone; CheckIt: /* Check to see if parameters are in */ if NtotalIn and PowerIn then do; put // "ERROR: Cannot use both Ntotal and Power statements."; stop; end; if not(NtotalIn or PowerIn) then do; put // "ERROR: Missing Ntotal or Power statement for this problem."; stop; end; if ANCVprob=1 and not(CvPCorIN or SDredIn) then do; put // "ERROR: ANCOVA problem not fully specified."; stop; end; if muProb and not(LogNProb) and not(SDIn) then do; put // "ERROR: MU problem requires SD (SIGMA)."; stop; end; if LogNProb and not(CVIn) then do; put // "ERROR: LogNormal problem requires CV."; stop; end; if muProb and LogNProb and (NumGrps = 1) and not(NullIn) then do; put // "ERROR: This LogNormal problem requires NULL mean."; stop; end; if GdFtProb and not(NullIn) then do; put // "ERROR: Goodness of fit problem requires NULL."; stop; end; if PaMuProb and not(CorrIn) then do; put // "ERROR: PairedMu problem requires CORR statement."; stop; end; if WlcxOrCt and not(NullIn) and NumGrps = 1 then do; put // "ERROR: One-group ordered categorial problem requires NULL."; stop; end; if (betLSPrb or RSqProb or SSH_Prob or PCorrPrb) and (NmParmIn = 0) then do; put // "ERROR: Missing NumParms statement for this problem."; stop; end; if (SSH_Prob or Chi2Prob) and (NexIn = 0) then do; put // "ERROR: Missing Nexemplary statement for this problem."; stop; end; if betLSPrb and SDxIn=0 then do; put // "ERROR: Missing SD_X statement for this problem."; stop; end; if betLSPrb and TolIn=0 then do; put // "ERROR: Missing TOLERANCE statement for this problem."; stop; end; if (muP_Prob ne 1) and (SDIn ne 1) then do; put // "ERROR: Missing SD (or SIGMA) statement for this problem."; stop; end; if muP_Prob and (SDPIn ne 1) then do; put // "ERROR: Missing SD(P) statement for this problem."; stop; end; if muP_Prob and (TrialsIn ne 1) then do; put // "ERROR: Missing Trials or Items statement for this problem."; stop; end; if (muP_Prob = 1) and (NumGrps = 1) and (NullIn = 0) then do; put // "ERROR: Must have NULL statement for this problem."; stop; end; if MVmuProb and not(SDIn) then do; put // "ERROR: MVmu problem requires SD (or SIGMA) statement."; stop; end; if MVmuProb and not(CorrIn) then do; put // "ERROR: MVmu problem requires CORR statement."; stop; end; if SampProb and not(CnfLevIn) then put // "NOTE: No confidence level has been specified." / " Only confidence levels of 0.50 (the median estimator)" / " and 0.75 will be computed."; if (SamFProb or SamTProb) and not(NmParmIn) then do; put // "ERROR: No NumParms has been specified." / " How many parameters are in the statistical model?"; stop; end; if AlphaIn = 0 then do; AlphaIn = 1; NumAlpha = 1; alphaV{1} = .05; put /