Can someone help convert this SAS simulation code to Stata?

[Suhail Doi]

Hi everyone,

I have a SAS code for simulating the effect size and standard error for a meta-analytic data-set of K studies over L iterations. Would be very grateful if someone could help convert this to Stata code or point me to some automated way of doing so.
Many thanks
Suhail Doi

%let OR = 2;
%let phi2 = 0;
%let K = 10;
%let L=1000;
*    Create a set of data with K points;
data output;
/*Start loop*/
do r=1 to &L;
     do idx = 1 to &K. ;
 
*    Step a;
           theta = log(&OR.);
*    Step b;
/*No Delaporte distribution implemented in SAS, using Native Binomial distribution instead*/
           ub = rand("Uniform");
           bl=100;
           bh=1000;
           n = bl + (bh-bl)*ub;
 
*    Step c;
 
           uc = rand("Uniform");
           cl=0.82;
           ch=0.993;
           p_nc = cl + (ch-cl)*uc;
 
 
    
*    Step d;
 
           ud = RAND('UNIFORM');
           dl=0.454;
           dh=0.545;
           tmpd=dl + (dh-dl)*ud;
           b = tmpd*n* p_nc;
           d = n* p_nc-b;
 
*    Step e;
 
           ue = RAND('UNIFORM');
           el=0.1;
           eh=0.55;
 
           tmpe=el + (eh-el)*ue;
           a=tmpe*n*(1 - p_nc);
           c = n * (1 - p_nc) - a;
 
*    Step f;
           if mod(idx,2)=1 then do;  
                p_2 = RAND('BETA', b, d);
                p_1 = (p_2 * exp(theta)) / (1 - p_2 + (p_2 * exp(theta)));
    
           end;
           if mod(idx,2)=0 then do;
                p_1 = RAND('BETA', a, c);
                p_2 = (p_1 / exp(theta)) / (1 - p_1 + (p_1 / exp(theta)));
           end;
 
*    Step g;
           a_s = p_1 * n * (1 - p_nc);
           c_s = (n * (1 - p_nc)) - a_s;
           b_s = p_2 * n *p_nc;
           d_s = (n * p_nc) - b_s;
 
*    Step h;
           p_1s = RAND('BETA', a_s, c_s);
           p_2s = RAND('BETA', b_s, d_s);
*    Step i;
           theta_d = log((p_1s / (1 - p_1s))/(p_2s / (1 - p_2s)));
 
*    Step j;
 
          uj = RAND('UNIFORM');
          jl=1;
           jh=9;
           q=(jl + (jh-jl)*uj)/10;
*    Step k;
           phi2 = 0.000000000000001 + (&phi2. * (1 - q) / q);
 
*    Step l;
           theta_hat = RAND('NORMAL', theta_d, phi2);
 
*    Step m;
           p_1ss = (p_2 * exp(theta_hat)) / (1 - p_2 + (p_2 * exp(theta_hat)));
*    Step n - repeat step g with p_1ss and p_2;
           a_ss = p_1ss * n * (1 - p_nc);
           c_ss = n * (1 - p_nc) - a_ss;
 
*    Step o;
           if a_ss gt 0 and c_ss gt 0 and b_s gt 0 and d_s gt 0 then do;
           var_theta_hat = 1 / a_ss + 1 / c_ss + 1 / b_s + 1 / d_s;
               se=sqrt(var_theta_hat);
           end;
           output;
     end;
end;
run;

[Christophe Kolodziejczyk]

See help simulate and help runiform

[Suhail Doi]

Btw, the code above was incorrect .....the correct code is:

%let OR = 2;
%let phi2 = 2;
%let K = 10;
%let L=1000;

*    Create a set of data with K points;
data output;
    /*Start loop*/
    do r=1 to &L;
        do idx = 1 to &K.;
            *    Step a;
            theta = log(&OR.);

            *    Step b;
            /*No Delaporte (0.1,8000,160) distribution implemented in SAS, using Native Poisson and gamma distributions instead*/
            gamma=8000*rangam(100,0.1);
            n = RAND('POISSON',gamma)+RAND('POISSON',160);

            *    Step c;
            uc = rand("Uniform");
            cl=0.82;
            ch=0.993;
            p_nc = cl + (ch-cl)*uc;

            *    Step d;
            ud = RAND('UNIFORM');
            dl=0.454;
            dh=0.545;
            tmpd=dl + (dh-dl)*ud;
            b = tmpd*n* p_nc;
            d = n* p_nc-b;

            *    Step e;
            ue = RAND('UNIFORM');
            el=0.1;
            eh=0.55;
            tmpe=el + (eh-el)*ue;
            a=tmpe*n*(1 - p_nc);
            c = n * (1 - p_nc) - a;

            *    Step f;
            if mod(idx,2) =0 then
                do;
                    p_2 = RAND('BETA', b, d);
                    p_1 = (p_2 * exp(theta)) / (1 - p_2 + (p_2 * exp(theta)));
                end;

            if mod(idx,2) ne 0 then
                do;
                    p_1 = RAND('BETA', a, c);
                    p_2 = (p_1 / exp(theta)) / (1 - p_1 + (p_1 / exp(theta)));
                end;

            *    Step g;
            a_s = (int(p_1 * n * (1 - p_nc)));
            c_s = (int((n * (1 - p_nc)) - a_s));
            b_s = (int(p_2 * n *p_nc));
            d_s = (int((n * p_nc) - b_s));

            if (a_s < 1 or c_s<1 or b_s<1 or d_s<1) then
                do;
                    a_s=a_s+0.5;
                    c_s=c_s+0.5;
                    b_s=b_s+0.5;
                    d_s=d_s+0.5;
                end;

            varx = (1/a_s) + (1/b_s) +(1/c_s) +(1/d_s);

            *    Step i;
            *    Step j;
            q=idx/11;

            *    Step k;
            phi2 = 0.000000000000001 + (&phi2. * (1 - q) / q);

            *    Step l;
            theta_hat = RAND('NORMAL', theta, sqrt(phi2+varx));

            *    Step m;
            p_1ss = (p_2 * exp(theta_hat)) / (1 - p_2 + (p_2 * exp(theta_hat)));

            *    Step n - repeat step g with p_1ss and p_2;
            a_ss = (int(p_1ss * n * (1 - p_nc)));
            c_ss = (int(n * (1 - p_nc) - a_ss));

            if (a_ss <1 or c_ss <1) then
                do;
                    a_ss=a_ss+0.5;
                    c_ss=c_ss+0.5;
                end;

            *    Step o;
            var_theta_hat = (1 / a_ss) + (1 / c_ss) + (1 / b_s) + (1 / d_s);
            se=sqrt(var_theta_hat);
            output;
        end;
    end;
run;

[Clyde Schechter]

I haven't used SAS in a very long time, but this SAS code is fairly straightforward. It builds up a data set with 1000 replications of 10 random observations generated according to the formulas shown there. The Stata version below does that, and saves the result in a data set called simulated_data.dta. In addition to the variables generated by the SAS code, it includes a variable called iteration, which indexes the 1,000 replications.

version 14.1

set more off

clear*

local OR 2
local phi2 2
local K 10
local L 1000

set seed 1234  // OR YOUR FAVORITE RANDOM NUMBER SITE

tempfile building
save `building', emptyok

forvalues r = 1/`L' {
    clear
    set obs `K'
    
    // STEP a
    local theta = log(`OR')
    
    // STEP b
    gen gamma = 8000*rgamma(100, 0.1)
    gen n = rpoisson(gamma) + rpoisson(160)
    
    // STEP c
    gen p_nc = 0.82 + (0.993-0.82)*runiform()
    
    // STEP d
    gen tmpd = 0.454 + (0.545-0.454)*runiform()
    gen b = tmpd*n*p_nc
    gen d = n*p_nc-b
    
    // STEP e
    gen tmpe = 0.1 + (0.55-0.1)*runiform()
    gen a = tmpe*n*(1-p_nc)
    gen c = n* (1-p_nc) - a
    
    // STEP f
    gen p_2 = rbeta(b, d) if mod(_n, 2) == 0
    gen p_1 = p_2*exp(`theta')/(1-p_2 + (p_2/exp(`theta'))) if mod(_n, 2) == 0
    replace p_1 = rbeta(a, c) if mod(_n, 2) == 1
    replace p_2 = (p_1/exp(`theta')) / (1-p_1 + p_1/exp(`theta')) if mod(_n, 2) == 1
    
    // STEP g
    gen a_s = int(p_1 * n * (1-p_nc))
    gen c_s = int((n*(1-p_nc)) - a_s)
    gen b_s = int(p_2*n*p_nc)
    gen d_s = int(n*p_nc - b_s)
    
    foreach v of varlist a_s b_s c_s d_s {
        replace `v' = `v' + 0.5 if min(a_s, b_s, c_s, d_s) < 1
    }
    gen varx = 1/a_s + 1/b_s + 1/c_s + 1/d_s
    
    // NO STEP h ???
    // STEP i DOES NOTHING???
    
    // STEP j
    gen q = _n/11
    
    // STEP k
    gen phi2 = 0.000000000000001 + (`phi2' * (1 - q) / q)
    
    // STEP l
    gen theta_hat = rnormal(`theta', sqrt(phi2+varx))
    
    // STEP m
    gen p_lss = (p_2*exp(theta_hat))/(1-p_2 + p_2*exp(theta_hat))
    
    // STEP n
    gen a_ss = int(p_lss*n*(1-p_nc))
    gen c_ss = int(n*(1-p_nc)-a_ss)
    foreach v of varlist a_ss c_ss {
        replace `v' = `v' + 0.5 if min(a_ss, c_ss) < 1
    }
    
    // STEP o
    gen var_theta_hat = 1/a_ss + 1/c_ss + 1/b_s + 1/d_s
    gen se = sqrt(var_theta_hat)
    
    gen int iteration = `r'
    order iteration, first
    
    append using `building'
    save `"`building'"', replace
}

use `building'
sort iteration
quietly compress
save simulated_data, replace

exit

Notes: The -set seed 1234- command just guarantees that you will always get the same results each time you run this. You don't have to use 1234--substituting any positive integer will generate a different set of results, but then you will get that same set of results each time.

There is no Step h in the SAS code. Did you leave something out?

Step i in the SAS code does nothing at all. Again, was something left out?

The SAS code seems to have a lot of redundant sets of parentheses. I have omitted most of those.

The code could be made somewhat more efficient. For example, most occurrences of `theta' are as exp(`theta'). Since local macro theta itself is defined as log(`OR'), all occurrences of exp(`theta') could be replaced by just `OR'. But I didn't bother doing that--the code doesn't take very long to run in any case.

I did this on version 14.1, and put that version at the top of the code. But as far as I can tell, it will run on many older versions of Stata. I don't remember at which version the current nomenclature for the random number functions was introduced (?9, ?10), but any version from that point on should be able to run this code if you change the version statement accordingly.

I have run this on my own Stata installation, so I know that it is free of syntax errors. But I can't guarantee that the code is correct. There are a lot of variable names that resemble each other, and I won't promise that I didn't substitute one for another here or there. I did my best to be careful. But since I don't really get what these variables are supposed to be, I don't have intuition to guide me that something "doesn't look right." So do double-check this code against the original: the comments in the code identify where each step is carried out, and it should be fairly simple for you to verify if I got everything right. Also, I would recommend that once you run it, you explore the resulting data to make sure that it all makes sense. As I don't know what these variables are, I have no intuition about that.

[Suhail Doi]

Dear Clyde,
Thanks very much - really appreciated. Was away for a while so unable to respond in a timely manner. I will go through and fix any variable issues you have pointed out and post here again.
Once again, thanks a lot
Regards
Suhail Doi

[Suhail Doi]

Hi all,

Here it is at last and this was made possible by Clyde's help (above) and David Fisher's admetan program without which this would have taken eight days to run. Now it runs in under 30 minutes
For the record, this simulation is under review at an epidemiology journal but given that it has been turned down by both AJE and EJE it seems right to deliver this to Stata users without whom this would not have been possible (if accepted somewhere a link will replace the code). Despite the fact that this code may change our scientific perspective forever yet fails to sway journal editors seems indicative of the problems scientists have to grapple with (thankfully not on Stata list).
Look forward to your thoughts
Regards
Suhail

NB to make it run faster change local iterations to a smaller number

**************************
**************************
*** Effect size is standardised mean difference***
*** tes (true effect size)***
*** maxscore (maximum score on the RoB scale)***
*** Beta distribution on Qi is meant to reflect variability
*** in assesment of Qi in real life ***
*** n is assumed equal in both arms***

**************************
**************************

clear all



*************************
** Initial parameters *** note tes = true effect size
*************************
local numb_studies 10
local tes 1
local runs 10
local iterations 5000
local maxscore 20

set output error

***********************
*** Data generation ***
***********************
tempfile building
save `building', emptyok
forvalues y = 1/`runs' {
forvalues meta = 1/`iterations' {

clear
set obs `numb_studies'
gen tes = `tes'

gen study_id = _n

***Study quality information***
***************************
gen qscore = 0+int((`maxscore'-0)*runiform() )
gsort -qscore
gen qscore_max=qscore[1]
gen qp = qscore/`maxscore'
gen qi = rbeta((qscore+0.15),(qscore_max-qscore+0.15))

***Study sample size***
********************
gen sample_size=(1+int((8-1)*runiform() ))*25

***Adding error***
********************
gen phi = `y'-1
gen var_bias = phi*((1-qp)/(qp))
gen serr = sqrt((4/sample_size)+(tes^2/(2*sample_size)))
gen tot_err = sqrt((serr^2)+var_bias)

***Effect size for each study***
**********************************
gen es = rnormal(tes,tot_err)

***Standard error for each study***
***************************************
replace serr = sqrt((4/sample_size)+(es^2/(2*sample_size)))



***************************
*** IVhet meta-analysis ***
***************************
admetan es serr, re(ivhet) nograph
gen es_ivhet = r(eff)
gen se_ivhet = r(se_eff)
gen lci_ivhet = es_ivhet - (1.96*se_ivhet)
gen uci_ivhet = es_ivhet + (1.96*se_ivhet)
gen width_ivhet = uci_ivhet - lci_ivhet
gen coverage_ivhet = 0
replace coverage_ivhet = 1 if tes >= lci_ivhet & tes<= uci_ivhet


************************
*** RE meta-analysis ***
************************
admetan es serr, re(re) nograph
gen es_re = r(eff)
gen se_re = r(se_eff)
gen tsq = r(tausq)
gen lci_re = es_re - (1.96*se_re)
gen uci_re = es_re + (1.96*se_re)
gen width_re = uci_re - lci_re
gen coverage_re = 0
replace coverage_re = 1 if tes >= lci_re & tes<= uci_re


************************
*** QE meta-analysis ***
************************
admetan es serr, qe(qi) nograph
gen es_qe = r(eff)
gen se_qe = r(se_eff)
gen lci_qe = es_qe - (1.96*se_qe)
gen uci_qe = es_qe + (1.96*se_qe)
gen width_qe = uci_qe - lci_qe
gen coverage_qe = 0
replace coverage_qe = 1 if tes >= lci_qe & tes<= uci_qe


***************
*** Dataset ***
***************
drop _*
keep in 1
gen meta = `meta'
gen runs = `y'
keep runs meta tes es_* se_* lci_* uci_* coverage_* tsq width_*
append using `building'
save `building', replace
}
}


********************************
*** Create output variables***
********************************
gen c_ivhet=.
gen width_ci_ivhet=.
gen bias2_ivhet=.
gen var_ivhet=.
gen mse_ivhet=.

gen c_re=.
gen width_ci_re=.
gen bias2_re=.
gen var_re=.
gen mse_re=.

gen c_qe=.
gen width_ci_qe=.
gen bias2_qe=.
gen var_qe=.
gen mse_qe=.

gen med_tsq =.

****************************
***Do analysis by runs***
****************************
forvalues x = 1/`runs' {

*** Analyse IVhet output ***

summ coverage_ivhet if `x'==runs, det
replace c_ivhet = r(mean) in `x'
summ width_ivhet if `x'==runs, det
replace width_ci_ivhet = r(p50) in `x'
summ es_ivhet if `x'==runs, det
replace bias2_ivhet = (r(mean) - tes)^2 in `x'
replace var_ivhet = (r(sd))^2 in `x'
replace mse_ivhet = bias2_ivhet + var_ivhet in `x'

*** Analyse RE output ***

summ coverage_re if `x'==runs, det
replace c_re = r(mean) in `x'
summ width_re if `x'==runs, det
replace width_ci_re = r(p50) in `x'
summ es_re if `x'==runs, det
replace bias2_re = (r(mean) - tes)^2 in `x'
replace var_re = (r(sd))^2 in `x'
replace mse_re = bias2_re + var_re in `x'

*** Analyse QE output ***

summ coverage_qe if `x'==runs, det  
replace c_qe = r(mean) in  `x'
summ width_qe if `x'==runs, det
replace width_ci_qe = r(p50) in  `x'
summ es_qe if `x'==runs, det
replace bias2_qe = (r(mean) - tes)^2 in `x'
replace var_qe = (r(sd))^2 in `x'
replace mse_qe = bias2_qe + var_qe in `x'


***Compute median tau squared***

summ tsq if `x'==runs, det
replace med_tsq = r(p50) in `x'
}


*********************
*** drop dataset ***
*********************

**drop tes es_ivhet se_ivhet lci_ivhet uci_ivhet coverage_ivhet es_re se_re tsq lci_re uci_re coverage_re es_qe se_qe lci_qe uci_qe coverage_qe meta runs width_ivhet width_re width_qe

set output proc

****************
*** Plots ***
****************

sort med_tsq

twoway (conn bias2_ivhet med_tsq) (conn  bias2_re med_tsq) (conn  bias2_qe med_tsq), sch(s1mono) name(bias) ylabel(, angle(0))
twoway (conn mse_ivhet med_tsq) (conn  mse_re med_tsq) (conn  mse_qe med_tsq), sch(s1mono) name(mse) ylabel(, angle(0))
twoway (conn width_ci_ivhet med_tsq, ) (conn  width_ci_re med_tsq) (conn  width_ci_qe med_tsq), sch(s1mono) name(width) ylabel(, angle(0))
twoway (conn c_ivhet med_tsq, yline(0.95)) (conn  c_re med_tsq) (conn  c_qe med_tsq), sch(s1mono) name(cover) ylabel(, angle(0))
graph combine bias mse width cover,  xcommon
graph drop bias width cover mse

[Suhail Doi]

Hi everyone,

Code updated due to slight changes in admetan syntax.

**************************
**************************
*** Effect size is standardised mean difference***
*** tes (true effect size)***
*** maxscore (maximum score on the RoB scale)***
*** Beta distribution on Qi is meant to reflect variability
*** in assesment of Qi in real life ***
*** n is assumed equal in both arms***

**************************
**************************

clear all



*************************
** Initial parameters *** note tes = true effect size
*************************
local numb_studies 10
local tes 1
local runs 10
local iterations 5000
local maxscore 20

set output error

***********************
*** Data generation ***
***********************
tempfile building
save `building', emptyok
forvalues y = 1/`runs' {
forvalues meta = 1/`iterations' {

clear
set obs `numb_studies'
gen tes = `tes'

gen study_id = _n

***Study quality information***
***************************
gen qscore = 0+int((`maxscore'-0)*runiform() )
gsort -qscore
gen qscore_max=qscore[1]
gen qp = qscore/`maxscore'
gen qi = rbeta((qscore+0.15),(qscore_max-qscore+0.15))

***Study sample size***
********************
gen sample_size=(1+int((8-1)*runiform() ))*25

***Adding error***
********************
gen phi = `y'-1
gen var_bias = phi*((1-qp)/(qp))
gen serr = sqrt((4/sample_size)+(tes^2/(2*sample_size)))
gen tot_err = sqrt((serr^2)+var_bias)

***Effect size for each study***
**********************************
gen es = rnormal(tes,tot_err)

***Standard error for each study***
***************************************
replace serr = sqrt((4/sample_size)+(es^2/(2*sample_size)))



***************************
*** IVhet meta-analysis ***
***************************
admetan es serr, model(ivhet) nograph
gen es_ivhet = r(eff)
gen se_ivhet = r(se_eff)
gen lci_ivhet = es_ivhet - (1.96*se_ivhet)
gen uci_ivhet = es_ivhet + (1.96*se_ivhet)
gen width_ivhet = uci_ivhet - lci_ivhet
gen coverage_ivhet = 0
replace coverage_ivhet = 1 if tes >= lci_ivhet & tes<= uci_ivhet


************************
*** RE meta-analysis ***
************************
admetan es serr, re nograph
gen es_re = r(eff)
gen se_re = r(se_eff)
gen tsq = r(tausq)
gen lci_re = es_re - (1.96*se_re)
gen uci_re = es_re + (1.96*se_re)
gen width_re = uci_re - lci_re
gen coverage_re = 0
replace coverage_re = 1 if tes >= lci_re & tes<= uci_re


************************
*** QE meta-analysis ***
************************
admetan es serr, qe(qi) nograph
gen es_qe = r(eff)
gen se_qe = r(se_eff)
gen lci_qe = es_qe - (1.96*se_qe)
gen uci_qe = es_qe + (1.96*se_qe)
gen width_qe = uci_qe - lci_qe
gen coverage_qe = 0
replace coverage_qe = 1 if tes >= lci_qe & tes<= uci_qe


***************
*** Dataset ***
***************
drop _*
keep in 1
gen meta = `meta'
gen runs = `y'
keep runs meta tes es_* se_* lci_* uci_* coverage_* tsq width_*
append using `building'
save `building', replace
}
}


********************************
*** Create output variables***
********************************
gen c_ivhet=.
gen width_ci_ivhet=.
gen bias2_ivhet=.
gen var_ivhet=.
gen mse_ivhet=.

gen c_re=.
gen width_ci_re=.
gen bias2_re=.
gen var_re=.
gen mse_re=.

gen c_qe=.
gen width_ci_qe=.
gen bias2_qe=.
gen var_qe=.
gen mse_qe=.

gen med_tsq =.

****************************
***Do analysis by runs***
****************************
forvalues x = 1/`runs' {

*** Analyse IVhet output ***

summ coverage_ivhet if `x'==runs, det
replace c_ivhet = r(mean) in `x'
summ width_ivhet if `x'==runs, det
replace width_ci_ivhet = r(p50) in `x'
summ es_ivhet if `x'==runs, det
replace bias2_ivhet = (r(mean) - tes)^2 in `x'
replace var_ivhet = (r(sd))^2 in `x'
replace mse_ivhet = bias2_ivhet + var_ivhet in `x'

*** Analyse RE output ***

summ coverage_re if `x'==runs, det
replace c_re = r(mean) in `x'
summ width_re if `x'==runs, det
replace width_ci_re = r(p50) in `x'
summ es_re if `x'==runs, det
replace bias2_re = (r(mean) - tes)^2 in `x'
replace var_re = (r(sd))^2 in `x'
replace mse_re = bias2_re + var_re in `x'

*** Analyse QE output ***

summ coverage_qe if `x'==runs, det  
replace c_qe = r(mean) in  `x'
summ width_qe if `x'==runs, det
replace width_ci_qe = r(p50) in  `x'
summ es_qe if `x'==runs, det
replace bias2_qe = (r(mean) - tes)^2 in `x'
replace var_qe = (r(sd))^2 in `x'
replace mse_qe = bias2_qe + var_qe in `x'


***Compute median tau squared***

summ tsq if `x'==runs, det
replace med_tsq = r(p50) in `x'
}


*********************
*** drop dataset ***
*********************

**drop tes es_ivhet se_ivhet lci_ivhet uci_ivhet coverage_ivhet es_re se_re tsq lci_re uci_re coverage_re es_qe se_qe lci_qe uci_qe coverage_qe meta runs width_ivhet width_re width_qe

set output proc

****************
*** Plots ***
****************

sort med_tsq

twoway (conn bias2_ivhet med_tsq) (conn  bias2_re med_tsq) (conn  bias2_qe med_tsq), sch(s1mono) name(bias) ylabel(, angle(0))
twoway (conn mse_ivhet med_tsq) (conn  mse_re med_tsq) (conn  mse_qe med_tsq), sch(s1mono) name(mse) ylabel(, angle(0))
twoway (conn width_ci_ivhet med_tsq, ) (conn  width_ci_re med_tsq) (conn  width_ci_qe med_tsq), sch(s1mono) name(width) ylabel(, angle(0))
twoway (conn c_ivhet med_tsq, yline(0.95)) (conn  c_re med_tsq) (conn  c_qe med_tsq), sch(s1mono) name(cover) ylabel(, angle(0))
graph combine bias mse width cover,  xcommon
graph drop bias width cover mse

Regards
Suhail