/******************************************************************************************* Created By: Brandy Ringham Date: 6/22/13 Description: Create power curve for hypothetical oral cancer biomarkers study. Use inputs from Elashoff et al., 2012. Fixed predictor = disease status (normal, case) Dependent variables = mRNA (molecule1, molecule2, molecule3); Hypothesis test = group difference for any mRNA *******************************************************************************************/ title1 h = 12pt "Oral Cancer Biomarkers Demonstration"; footnote1 "&sysdate"; footnote2 "OCBiomarkersDemonstration.sas"; /*initiate SAS/IML*/ proc iml; /*define planned per group sample size*/ unadjustRepNt = 75; /*vector of anticipated proportion of missing data*/ piVec = { 0, .05, .1 }; /*correlation between repeated measures*/ rho = .4; /*define 7 inputs of power analysis*/ /*1. alpha*/ alpha = .05; /*2. between subject contrast matrix*/ /*dimensions 1x2*/ c = { 1 -1 }; /*3. within subject contrast matrix*/ /*dimensions 3x2*/ u = { 1 0 0, 0 1 0, 0 0 1 }; /*4. sigma*/ /*see below--values will vary so it resides in the do loop*/ /*5. beta */ /*see below--values will vary so it resides in the do loop*/ /*6. null hypothesis*/ /*set by POWERLIB as a 1x2 matrix of 0's by default*/ /*7. essenceX (Muller and Stewart, 2006)*/ /*dimensions 2x2*/ essenceX = { 1 0, 0 1 }; /*calculate total planned sample size - unadjusted*/ unadjustNt = unadjustRepNt * nrow( essenceX ); /*include NMK code*/ /*NMK module calculate adjusted sample size, or E(Nmk)*/ %include "NMK.SAS" / nosource2; /*include POWERLIB code*/ /*POWERLIB module calculates power*/ %include "POWERLIB22.IML" / nosource2; /*set options for POWERLIB*/ opt_off = { HF UN BOX WARN }; opt_on = { UNIFORCE DS FRACREPN NOPRINT }; /*loop through rows of proportion missing vector*/ do piID = 1 to nrow( piVec ); /*define proportion missing as the value in the current row of the vector*/ pi = piVec[ piID, ]; /*loop through different scale factors for the beta matrix*/ do k = 0 to 2 by .005; /*variances are from Table 3, Elashoff et al., 2012*/ var = 2.9**2; /*define sigma matrix*/ oneMat = j( 3, 1, 1 ); sigma = var * ( oneMat * oneMat` * rho + I( 3 ) * ( 1 - rho ) ); /*cell means are from for 1) IL-1B, 2) IL-8, and 3) SAT, Cohort 4 reported in Table 3, Elashoff et al., 2012*/ /*first set means for normals*/ mu1n = 20.1; mu2n = 19.8; mu3n = 21.3; /*define difference between normals and cases*/ delta1 = -1.3; delta2 = -2.1; delta3 = -1.4; /*calculate mean for cases based on scale factor*/ mu1c = mu1n + k * delta1; mu2c = mu2n + k * delta2; mu3c = mu3n + k * delta3; /*define beta matrix from means*/ beta = ( mu1c || mu2c || mu3c ) // ( mu1n || mu2n || mu3n ); /*calculate sample sizes adjusted for anticipated amount of missing data*/ /*use Nm1*/ call NMK( unadjustNt, 2, pi, 1, 0, ENmk, error ); /*calculated adjusted number of participants per group*/ /*set repn equal to adjusted repn to pass into POWERLIB*/ repn = ENmk / nrow( essenceX ); /*save power results in a dataset*/ dsname = { work adjustPower }; /*calculate power adjusted for missing data using POWERLIB*/ run power; use work.adjustPower; /*create an IML matrix from the POWERLIB output dataset*/ read all var{ ALPHA TOTAL_N POWER_MULT EPSILON EXEPS_GG POWER_GG } into adjustPower[ colname = name ]; /*delete dataset after use*/ call delete( "work", "adjustPower" ); /*row of output containing experimental conditions and power results*/ powerRow = rho || k || pi || repn || 2 || adjustPower[ , 1:3 ]; /*create matrix of results over all proportion missing*/ power = power // powerRow; /*free variables to be reused*/ free powerRow ENmk repn dsname adjustPower; end; /*end k loop*/ end; /*end pi loop*/ /*list of column names for output dataset*/ expNames = { "rho" "k" "PropMiss" "unadjustRepNt" "NumDepVars" }; allNames = expNames || name[ , 1:3 ]; /*output matrix as SAS dataset*/ create out01.appPow from power[ colname = allNames ]; append from power; quit; /*set size of graph*/ goptions hsize = 4in vsize = 4in; /*plot power curve*/ proc gplot data = out01.appPow; title2 h = 12pt "Power Curve"; /*plot power by scale factor*/ plot power_mult * k = propMiss / noframe vaxis = axis1 haxis = axis2 legend = legend1 href = 1 vref = .9 cvref = red chref = red; /*define linetypes for different proportion missing*/ symbol1 i = j l = 1 c = black w = 1; symbol2 i = j l = 41 c = black w = 1; symbol3 i = j l = 35 c = black w = 1; /*define axes options*/ axis1 order = ( 0 to 1 by .2 ) minor = none major = none label = ( angle = 90 font = times h = 12pt "Power" ) value = ( h = 12pt font = times ); axis2 order = ( 0 to 2 by .5 ) minor = none major = none label = ( font = times h = 12pt "k" ) value = ( h = 12pt font = times ); /*define legend*/ legend1 value = ( height = 12pt font = "times" "0%" "5%" "10%" ) position = ( inside bottom right ) label = ( font = times h = 12pt position = top justify = center "Percent Missing" ) mode = protect across = 1 down = 3 frame shape = line( .3in ); run; quit; title; footnote;