Appendix: Sample SAS Codes for PMM Approximation and Bayesian Analysis

Appendix: Sample SAS codes for PMM Approximation and Bayesian Analysis

Assume the dataset “InDat” contain variables: patient (patient id), therapy (treatment code: DRUG or PLACEBO), basval (baseline value), y1, y2, y3, and y4 (repeated measures after baseline). Assume monotone missing and no missing for y1.

**** PMM Approximation using mixed model estimates ****;

datalda; *** longitudinal data for MMRM mixed model ***;

setInDat;

arrayyy{4} y1 y2 y3 y4;

do time=1to4;

y=yy{time};

output;

end;

drop y1 – y4;

run;

dataforMCMC; *** dataset for Bayesian MCMC analysis ***;

setInDat end=eof;

retainntrt0 np1 0 np2 0 np3 0;

*total subjects in Drug, and numbers of missing in patterns 1, 2, 3 *;

trt=(therapy='DRUG');

iftrt=1thendo;

ntrt=ntrt+1;

if y2=.thendo; np1=np1+1; end; ** missing pattern 1 **;

elseif y3=.thendo; np2=np2+1; end; ** pattern 2 **;

elseif y4=.thendo; np3=np3+1; end; ** pattern 3 **;

end;

ifeofthendo;

mp1=round(np1/ntrt,.00000001); ** missing proportion in pattern 1 **;

mp2=round(np2/ntrt,.00000001); ** missing proportion in pattern 2 **;

mp3=round(np3/ntrt,.00000001); ** missing proportion in pattern 3 **;

mp4=1-mp1-mp2-mp3; ** proportion of completers **;

callsymput('ntrt', ntrt); *** sample size in Drug ***;

callsymput('mp1', mp1);

callsymput('mp2', mp2);

callsymput('mp3', mp3);

callsymput('mp4', mp4);

end;

keep patient basval y1 y2 y3 y4trt;

run;

**** Primary analysis using MMRM mixed model ****;

**** Coding for Therapy: Drug should be ordered before Placebo ****;

procmixeddata=lda;

class therapy time;

model y=basval*time therapy*time / nointscovbddfm=kr;

estimate"MAR" therapy*time 0001 000 -1;

lsmeans therapy*time / ecov;

repeated time / subject=patient type=un r=1;

** assume first subject is a completer, otherwise change “r=” option

to a patient of completer so the covariance matrix can be outputed **;

odsoutputlsmeans=lsmR=covestimates=mix;

title"MMRM analysis with mixed model";

run;

**** PMM approximation for CIR, CR, and JR ****;

prociml;

c_mar = {0001 000 -1}; * MAR analysis at last time point *;

c_jr = {000 &mp4 000 -&mp4}; * JR analysis at last time point *;

c_cir = {&mp1 &mp2 &mp3 &mp4 -&mp1 -&mp2 -&mp3 -&mp4}; * for CIR *;

usecov; read all var {col1 col2 col3 col4} into v;

uselsm; read all var{estimate} into beta;

uselsm;read all var{cov1 cov2 cov3 cov4 cov5 cov6 cov7 cov8}into w;

v_mar = (c_mar)*w*T(c_mar);

out=(c_mar*beta)||sqrt(v_mar); *** Mixed model estimate and SE ***;

v_Jr = (c_jr)*w*T(c_jr) + &mp4*(1-&mp4)*(beta[8]-beta[4])**2/&ntrt;

*** combined variance for JR, ntrt=sample size in Drug group ***;

out=out//( (c_jr*beta)||sqrt(v_jr) ); ** JR estimate and SE **;

ww=(c_cir[1:4])*T(c_cir[1:4]);

ww=diag(c_cir[1:4]) - ww;

diff = beta[1:4] - beta[5:8];

v_cir = (c_cir)*w*T(c_cir) + T(diff)*ww*diff/&ntrt;

*** combined variance for CIR, ntrt= sample size in Drug group ***;

out=out//( (c_cir*beta)||sqrt(v_cir) ); ** CIR estimate and SE **;

mp1=&mp1*v[4,1]/v[1,1] + &mp2*(v[4,1:2]*inv(V[1:2,1:2]))[1] +

&mp3*(v[4,1:3]*inv(V[1:3,1:3]))[1];

mp2=&mp2*(v[4,1:2]*inv(V[1:2,1:2]))[2] +

&mp3*(v[4,1:3]*inv(V[1:3,1:3]))[2];

mp3=&mp3*(v[4,1:3]*inv(V[1:3,1:3]))[3];

c_cr=(mp1)||(mp2)||(mp3)||{&mp4};

c_cr = c_cr||(-c_cr);

diff=( (v[4,1]/v[1,1])#beta[1]//

((v[4,1:2]*inv(V[1:2,1:2]))*beta[1:2])//

((v[4,1:3]*inv(V[1:3,1:3]))*beta[1:3])//beta[4] ) -

( (v[4,1]/v[1,1])#beta[5]//

((v[4,1:2]*inv(V[1:2,1:2]))*beta[5:6])//

((v[4,1:3]*inv(V[1:3,1:3]))*beta[5:7])//beta[8] );

v_cr = c_cr*w*T(c_cr) + T(diff)*ww*diff/&ntrt;

out=out//( (c_cr*beta)||sqrt(v_cr) );

out = {1,2,3,4}||out;

create sum1 from out[colname={'Code''estimate''Se'}]; appendfrom out;

quit;

dataoutall; set sum1;

if code=1then Method='MAR_MMRM ';

elseif code=2then Method='JR_Apprx ';

elseif code=3then Method='CIR_Apprx ';

elseif code=4then Method='CR_Apprx ';

drop code;

run;

**** Bayesian Analysis with MCMC ****;

data_null_; set mix;

callsymput('mixDF', df);

*** get DF from MMRM mixed model analysis ***;

*** get Dirichlet posterior parameters using Geofrey's prior ***;

nn1=&ntrt*&mp1+0.5;

nn2=&ntrt*&mp2+0.5;

nn3=&ntrt*&mp3+0.5;

nn4=&ntrt*&mp4+0.5;

callsymput('nn1', nn1);

callsymput('nn2', nn2);

callsymput('nn3', nn3);

callsymput('nn4', nn4);

run;

*** Bayesian approach ***;

odshtml;

odsgraphics;

ProcMCMCdata=forMCMCnbi=2000ntu=2000nmc=200000seed=7892354thin=10propcov=CONGRA diag=allplots(smooth)=allmonitor=(MAR_MCMC CR__MCMC JR__MCMC CIR_MCMC);

arrayyy{4} y1 y2 y3 y4;

array mu{4} mu1 - mu4;

arraymp[4] mp1 - mp4;

arraynn[4] (&nn1 &nn2 &nn3 &nn4);

arraymup{4} mup1-mup4; *** mean for placebo ***;

arraymut{4} mut1-mut4; *** mean for test drug ***;

array beta[4] beta1-beta4; *** slope for baseline over time ***;

arraycov[4,4];

arraycov_a[3,3]; *** sub matrix for lower dimension ***;

arraycov_b[2,2];

arraybcov_a[1,3]; *** sub vector for cov of y4 with y1-y3 ***;

arraybcov_b[1,2];

arrayiv_a[3,3]; *** inverse matrix for cov_a ***;

arrayiv_b[2,2];

arraybiv_a[1,3]; *** bcov_a*iv_a ***;

arraybiv_b[1,2];

arraymtp_a[3]; *** subvector for mut - mup ***;

arraymtp_b[2];

array bmu3[1]; *** biv_a*mtp_a = bcov_a*inv(cov_a)*(mut-mup) ***;

array bmu2[1];

array s0[4,4] (1000 0100 0010 0001);

BEGINNODATA;

MAR_MCMC = mut4-mup4; *** mean difference at time 4 under MAR ***;

doi=1to3; mtp_a[i]=mut[i]-mup[i]; bcov_a[1,i]=cov[4,i];

do j=1to3; cov_a[i,j]=cov[i,j]; end;

end;

doi=1to2; mtp_b[i]=mut[i]-mup[i]; bcov_b[1,i]=cov[4,i];

do j=1to2; cov_b[i,j]=cov[i,j]; end;

end;

callinv(cov_a,iv_a);

callmult(bcov_a, iv_a, biv_a);

callmult(biv_a, mtp_a, bmu3);

callinv(cov_b,iv_b);

callmult(bcov_b, iv_b, biv_b);

callmult(biv_b, mtp_b, bmu2);

bmu1 = cov[1,4]*(mut1-mup1)/cov[1,1];

CR__MCMC = (mp4*(mut4-mup4) + mp3* bmu3[1]

+ mp2* bmu2[1] + mp1* bmu1); *** copy reference ***;

JR__MCMC = mp4*(mut4-mup4); *** jump to reference ***;

CIR_MCMC = (mp4*(mut4-mup4) + mp3* (mut3-mup3)

+ mp2*(mut2-mup2) + mp1*(mut1-mup1)); *** CIR ***;

ENDNODATA;

*** initial values ***;

parms mup1 -1.6 mup2 -4.5 mup3 -6.1 mup4 -7.1

mut1 -2.0 mut2 -5.5 mut3 -8.0 mut4 -9.6;

parms beta1 -1. beta2 -1. beta3 -1.0 beta4 -1.4;

parmscov;

parmsmp {&mp1 &mp2 &mp3 &mp4};

priormp ~ dirich(nn);

priormup: mut: beta: ~ normal(0, var=10000);

priorcov ~ iwish(4,S0);

do i = 1to4;

mu[i] = mup[i]*(1-trt)+mut[i]*trt + beta[i]*basval;

end;

modelyy ~ mvn(mu, cov);

title"Bayesian MCMC analysis";

odsoutputPostSumInt=sum2;

run;

**** Combine results and use DF from mixed model ****;

dataoutall;

setoutall

sum2(rename=(mean=Estimate stddev=se parameter=Method)

drop=N HPDLowerHPDUpper);

DF = &mixDF;

probt=2-2*probt(abs(estimate/se), df);

LCL = estimate - tinv(0.975, DF)*se;

UCL = estimate + tinv(0.975, DF)*se;

run;

odslisting;

procprintdata=outall;

var Method estimate se probt LCL UCL;

title"Summary of Analysis Results";

run;

1