/* INPUT PARAMETERS nreps = number of simulations nclst = number of clusters numobs = number of observations per cluster intvar1 = variance of random cluster effects (intercepts) sigeps1 = residual variance */ %macro varcompsim(nreps, nclst, numobs, intvar1, sigeps1); %do reps = 1 %to &nreps; /* Simulate data based on the specified model */ data simsample (drop = j k u); do j = 1 to &nclst; clstid = j; u = sqrt(&intvar1)*rannor(0); /* u ~ N(0,intvar1) */ do k = 1 to &numobs; y = 45 + u + sqrt(&sigeps1)*rannor(0); /* e ~ N(0,sigeps1) */ output; end; end; run; /* Analyze the simulated data set, and output the p-value from the LRT */ proc glimmix data = simsample; class clstid; model y = / solution; random int / subject = clstid; covtest glm; ods output covtests = lrtresult; run; /* Determine whether or not the LRT was significant */ %if &reps = 1 %then %do; data lrtresults (keep = ProbChiSq power_cont); set lrtresult; if ProbChiSq < 0.05 then power_cont = 1; else power_cont = 0; run; %end; %else %do; data lrtresult (keep = ProbChiSq power_cont); set lrtresult; if ProbChiSq < 0.05 then power_cont = 1; else power_cont = 0; run; data lrtresults; set lrtresults lrtresult; run; %end; %end; /* End of replicated simulations */ /* Compute the proportion of tests that were significant */ proc means data = lrtresults mean; var power_cont; run; %mend varcompsim; %varcompsim(nreps=100, nclst=30, numobs=40, intvar1=2, sigeps1=64);