option ps=500 ls=90 nonumber nodate; /* This program is designed for between-clinical centre heterogeniety in the probability of a binary outcome in the Standard Arm which follows a beta distribution with a given mean and coeffecient of variation. Simulation results are analysed using the risk difference. */ %let reps1 = 10000; * reps1 = number of repetitions in simulation of the trial for suroggate outcome; %let reps2 = 10000; * reps2 = number repetitions in sampling from posterior distributions; %let a1 = 167; %let b1 = 365; * for beta distribution for Pr(PPM|ptb); %let a2 = 6; %let b2 = 929; * for beta distribution for Pr(PPM|tb); %let piSmean = 0.3; * piSmean = mean prob. of BAD outcome on Standard; %let coefvarS = 0.3; * the between-center distribution of the prob. of outcome on Standard is beta with mean=piSmean and CV=coefvarS if coefvarS = 0 then there is no between-center variation with respect to the prob. of outcome on Standard; %let piTmean = 0.225; * piTmean = mean prob. of BAD outcome on Treatment; %let coefvarT = 0.0; * the between-center distribution of the prob. of outcome on Treatment is beta with mean=piTmean and CV=coefvarT if coefvarT = 0 then there is no between-center variation with respect to the prob. of outcome on Treatment; %let minSS = 6; * minimum sample size per clinical center; %let all_to_CS = 'random'; * if equal 'even' then patients evenly allocated across clinical centers if equal 'random' then patients randomly allocated across clinical center; %let strat_by_CS = 'yes'; * if equal 'yes' then randomization is stratified by clinical center if equal 'no' then randomization is NOT stratified by clinical center; %let k = 80; * k = number of clinical centers; %let NperArm = 700; * NperArm = number of patients per arm Must be at least 6*k Should be multiples of k if all_to_CS = 'even'; %macro RD; proc iml; print "Number of repetitions in simulation of the trial for suroggate outcome:" &reps1; print "Number repetitions in sampling from posterior distributions:" &reps2; print "Mean prob. of BAD suroggate outcome in the Standard arm:" &piSmean; print "CV in the Standard arm:" &coefvarS; print "Mean prob. of BAD suroggate outcome in the Treatment arm:" &piTmean; print "CV in the Treatment arm:" &coefvarT; print "Number of clinical centers:" &k; print "Number of patients per arm:" &NperArm; n = j(&k,1,0); bb = j(&k,1,0); piS = j(&k,1,0); piT = j(&k,1,0); nS = j(&k,1,0); nS1 = j(&k,1,0); nS0 = j(&k,1,0); nT = j(&k,1,0); nT1 = j(&k,1,0); nT0 = j(&k,1,0); pT = j(&k,1,0); pS = j(&k,1,0); d = j(&k,1,0); v = j(&k,1,0); w = j(&k,1,0); wAA = j(&k,1,0); vStar = j(&k,1,0); wStar = j(&k,1,0); dwStar = j(&k,1,0); vStarF = j(&k,1,0); wStarF = j(&k,1,0); dwStarF = j(&k,1,0); outout = j(&reps1,1,0); delta = j(&reps2,1,0); Ntotal = 2*&NperArm; *******************************************************************; * To determine aS, bS, modeS and percentiles of the Beta distribution and to verify mean and coefficient of variation for Standard arm; if &coefvarS > 0 then do; * Loop 1A; oddsS = &piSmean/(1 - &piSmean); bS = ( 1 - (&coefvarS**2)*oddsS ) / ( (&coefvarS**2)*oddsS*(oddsS + 1) ); aS = oddsS*bS; betaABs = beta(aS,bS); meanS = aS/(aS+bS); varS = (aS*bS)/((aS+bS)*(aS+bS)*(aS+bS+1)); medianS = (aS - 1/3)/(aS + bS - 2/3); modeS = (aS - 1)/(aS + bS - 2); maxS = (modeS**(aS-1))*((1 - modeS)**(bS-1))/betaABs; cvS = sqrt(varS)/meanS; percentile_05S = betainv(0.05,aS,bS); percentile_95S = betainv(0.95,aS,bS); probGpiTmean = 1 - probbeta(&piTmean,aS,bS); print "5th percential for the prob. of a BAD outcome in the Standard arm:" percentile_05S; print "95th percential for the prob. of a BAD outcome in the Standard arm:" percentile_95S; end; * Loop 1A *******************************************************************; * To determine aT, bT, modeT and percentilesT of the Beta distribution and to verify mean and coefficient of variation for Treatment arm; if &coefvarT > 0 then do; * Loop 1B; oddsT = &piTmean/(1 - &piTmean); bT = ( 1 - (&coefvarT**2)*oddsT ) / ( (&coefvarT**2)*oddsT*(oddsT + 1) ); aT = oddsT*bT; betaABt = beta(aT,bT); meanT = aT/(aT+bT); varT = (aT*bT)/((aT+bT)*(aT+bT)*(aT+bT+1)); medianT = (aT - 1/3)/(aT + bT - 2/3); modeT = (aT - 1)/(aT + bT - 2); maxT = (modeT**(aT-1))*((1 - modeT)**(bT-1))/betaABt; cvT = sqrt(varT)/meanT; percentile_05T = betainv(0.05,aT,bT); percentile_95T = betainv(0.95,aT,bT); print "5th percential for the prob. of a BAD outcome in the Treatment arm:" percentile_05T; print "95th percential for the prob. of a BAD outcome in the Treatment arm:" percentile_95T; end; * Loop 1B *******************************************************************; do ii = 1 to &reps1; * Loop 2; * To allocate patients between the k clinical centers; n = j(&k,1,&minSS); if &all_to_CS = 'even' then do; * Loop 3; nSum = 0; do i=1 to &k; * Loop 4; n[i] = round(Ntotal/&k,1); nSum = nSum + n[i]; end; * Loop 4; end; * Loop 3; if &all_to_CS = 'random' then do; * Loop 5; bbSum = 0; do i=1 to &k; * Loop 6; bb[i] = ranuni(0); bbSum = bbSum + bb[i]; end; * Loop 6; nSum = 0; do i=1 to &k; * Loop 7; n[i] = n[i] + round((Ntotal - &k*&minSS)*bb[i]/bbSum,1); nSum = nSum + n[i]; end; * Loop 7; end; * Loop 5; dwSum = 0; wSum = 0; wSqdSum = 0; piSsum = 0; *******************************************************************; do i=1 to &k; * Loop 8; * To randomize patients between arms, within clinical centre; if &strat_by_CS = 'yes' then nS[i] = round(n[i]/2,1); if &strat_by_CS = 'no' then nS[i] = round(ranbin(0,n[i],0.5),1); if nS[i] = 0 then nS[i] = round(n[i]/2,1); nT[i] = n[i] - nS[i]; *******************************************************************; * To generate a Beta(a,b) random variable for the ith clinical centre for Standard arm; if &coefvarS > 0 then do; * Loop 9A; U = 1; YY = 0; do while (U > YY); * Loop 10; Y = ranuni(0); U = ranuni(0); YY = (Y**(aS-1))*((1 - Y)**(bS-1))/(maxS*betaABs); end; * Loop 10; piS[i] = Y; end; * Loop 9A; *******************************************************************; * To generate a Beta(a,b) random variable for the ith clinical centre for Treatment arm; if &coefvarT > 0 then do; * Loop 9B; U = 1; YY = 0; do while (U > YY); * Loop 10; Y = ranuni(0); U = ranuni(0); YY = (Y**(aT-1))*((1 - Y)**(bT-1))/(maxT*betaABt); end; * Loop 10; piT[i] = Y; end; * Loop 9B; *******************************************************************; if &coefvarS = 0 then piS[i] = &piSmean; if &coefvarT = 0 then piT[i] = &piTmean; * To generate those patients who had preterm birth; if nS[i] > 0 then nS1[i] = ranbin(0,nS[i],piS[i]); nS0[i] = nS[i] - nS1[i]; if nT[i] > 0 then nT1[i] = ranbin(0,nT[i],piT[i]);; nT0[i] = nT[i] - nT1[i]; end; * Loop 8; *******************************************************************; do i=1 to &k; * Loop 11; * To calculate fixed effects solution; pT[i] = nT1[i]/nT[i]; pS[i] = nS1[i]/nS[i]; d[i] = pT[i] - pS[i]; v[i] = max(0.01,pT[i])*(1 - min(0.99,pT[i]))/nT[i] + max(0.01,pS[i])*(1 - min(0.99,pS[i]))/nS[i]; w[i] = 1/v[i]; dwSum = dwSum + d[i]*w[i]; wSum = wSum + w[i]; wSqdSum = wSqdSum + w[i]**2; piSsum = piSsum + piS[i]; end; * Loop 11; piSaverage = piSsum/&k; RDfixed = dwSum/wSum; VarRDfixed = 1/wSum; *******************************************************************; * To calculate random effects solution; piSvar = 0; AAsum = 0; do i=1 to &k; * Loop 12; piSvar = piSvar + (piS[i] - piSaverage)**2; AAsum = AAsum + w[i]*((d[i] - RDfixed)**2); end; * Loop 12; piScoefVar = sqrt(piSvar/&k)/piSaverage; sigmaSQ = max(0, (AAsum - (&k - 1))/(WSum - WSqdSum/Wsum)); dwStarSum = 0; wStarSum = 0; do i=1 to &k; * Loop 13; vStar[i] = v[i] + sigmaSQ; wStar[i] = 1/vStar[i]; dwStar[i] = d[i]*wStar[i]; wStarSum = wStarSum + wStar[i]; dwStarSum = dwStarSum + dwStar[i]; end; * Loop 13; RDrandom = dwStarSum/wStarSum; VarRDrandom = 1/wStarSum; do iii = 1 to &reps2; * Loop 14; *******************************************************************; * To generate a Beta(a1,b1) random variable for prob(PMM|ptb); beta1 = beta(&a1,&b1); mode1 = (&a1 - 1)/(&a1 + &b1 - 2); max1 = (mode1**(&a1-1))*((1 - mode1)**(&b1-1))/beta1; U1 = 1; YY1 = 0; do while (U1 > YY1); * Loop 10A; Y1 = ranuni(0); U1 = ranuni(0); YY1 = (Y1**(&a1-1))*((1 - Y1)**(&b1-1))/(max1*beta1); end; * Loop 10A; prppm1 = Y1; *******************************************************************; * To generate a Beta(a2,b2) random variable for prob(PMM|tb); beta2 = beta(&a2,&b2); mode2 = (&a2 - 1)/(&a2 + &b2 - 2); max2 = (mode2**(&a2-1))*((1 - mode2)**(&b2-1))/beta2; U2 = 1; YY2 = 0; do while (U2 > YY2); * Loop 10B; Y2 = ranuni(0); U2 = ranuni(0); YY2 = (Y2**(&a2-1))*((1 - Y2)**(&b2-1))/(max2*beta2); end; * Loop 10B; prppm2 = Y2; *******************************************************************; RDsim = RDrandom + sqrt(VarRDrandom)*rannor(0); delta[iii] = RDsim*(prppm1 - prppm2); * Output to SAS database; end; * Loop 14; propNeg = 0; deltaRanks = rank(delta); do iFind = 1 to &reps2; if delta[iFind] < 0 then propNeg = propNeg + 1; end; outout[ii,1] = propNeg/&reps2; end; * Loop 2; vars = "propNeg"; create delta from outout [colname=vars]; append from outout; *******************************************************************; quit; data delta; set delta; label propNeg = "(Proportion of Negative Delta's)"; run; proc means data=delta n p5 mean median; var propNeg; title "Statistics for the Proportion of Negative Delta's"; run; %mend; %RD;