SAS syntax to fit the MREM model by writing the likelihood in the SAS procedure NLMIXED

###

*PROC NLMIXED = call SAS procedure NLMIXED, and the options included represents:

df = degree of freedom, we specified 1000 to get Wald test in stead of t-test

qpoints = number of quadrature points, if not specified SAS automatically does

miniter = is the minimum number of iteration;

PROC NLMIXED DATA = cage DF = 1000 MINITER = 30 QPOINTS = 5;

*Specifies initial values for the maximum likelihood estimates

if not specified, all the initial values automatically assigned to be one;

PARMS ma = 2.6 b = 0.8 b2 = 0 mxi1 = -5.2 mxi2 = -3.7 mxi3 = -2.2 mxi4 = -1.0

va = 0.4 cavdi = -0.2 vdi = 0.5 vdij = 0.02;

*Constraint to make sure that xi's are in the right order;

BOUNDS mxi1-mxi2< = 0, mxi2-mxi3< = 0, mxi3-mxi4< = 0;

*Between-studies model (equation (5));

eta1 = a + b*xi1 + b2*z1;

eta2 = a + b*xi2 + b2*z1;

eta3 = a + b*xi3 + b2*z1;

eta4 = a + b*xi4 + b2*z1;

*Within-study model (equation (8));

p01 = 1/(1+exp(-(xi1)));

p02 = 1/(1+exp(-(xi2))) - 1/(1+exp(-(xi1)));

p03 = 1/(1+exp(-(xi3))) - 1/(1+exp(-(xi2)));

p04 = 1/(1+exp(-(xi4))) - 1/(1+exp(-(xi3)));

p05 = 1 - 1/(1+exp(-(xi4)));

*Within-study model (equation (9));

p11 = 1/(1+exp(-(eta1)));

p12 = 1/(1+exp(-(eta2))) - 1/(1+exp(-(eta1)));

p13 = 1/(1+exp(-(eta3))) - 1/(1+exp(-(eta2)));

p14 = 1/(1+exp(-(eta4))) - 1/(1+exp(-(eta3)));

p15 = 1 - 1/(1+exp(-(eta4)));

*Log-likelihood (equation (10));

if (p01^ = 0 and p02^ = 0 and p03^ = 0 and p04^ = 0 and p05^ = 0 and

p11^ = 0 and p12^ = 0 and p13^ = 0 and p14^ = 0 and p15^ = 0) then

ll = n01*log(p01)+n02*log(p02)+n03*log(p03)+n04*log(p04)+n05*log(p05)+

n11*log(p11)+n12*log(p12)+n13*log(p13)+n14*log(p14)+n15*log(p15);

else ll = -1**100;

*SAS maximized the likelihood using the general distribution;

MODEL n11 ~general(ll);

*Random effects parameters alpha_i and xi_ij's are normally distributed around their mean and between-studies variance covariance matrix (equation 7);

RANDOM a xi1 xi2 xi3 xi4 ~normal([ma, mxi1, mxi2, mxi3, mxi4],

[va,

cavdi, vdi+vdij,

cavdi, vdi, vdi+vdij,

cavdi, vdi, vdi, vdi+vdij,

cavdi, vdi, vdi, vdi, vdi+vdij])

SUBJECT = study OUT = study_specific; *Empirical Bayes estimate;

*Estimate sensitivities and specificities at each of known thresholds (1,2,3,4);

ESTIMATE 'sensitivity_1' exp(ma + mxi1)/(1 + exp(ma + mxi1));

ESTIMATE 'sensitivity_2' exp(ma + mxi2)/(1 + exp(ma + mxi2));

ESTIMATE 'sensitivity_3' exp(ma + mxi3)/(1 + exp(ma + mxi3));

ESTIMATE 'sensitivity_4' exp(ma + mxi4)/(1 + exp(ma + mxi4));

ESTIMATE 'specificity_1' exp(1)/(1 + exp(mxi1));

ESTIMATE 'specificity_2' exp(1)/(1 + exp(mxi2));

ESTIMATE 'specificity_3' exp(1)/(1 + exp(mxi3));

ESTIMATE 'specificity_4' exp(1)/(1 + exp(mxi4));

run;