Generate simulated correlated binary variables

[Jen Ward]

Hi all,

I am interesting in creating a simulated dataset; I would like to generate two binary variables which are correlated (rho=0.4) and I thought the -corr2data- command could achieve this; however, this does not seem to work well as the generated variables are not correlated. Any advice on what's going wrong?

clear

set seed 6710789
local n = 1000
local intercept6 = 0.20
local intercept12 = 0.35
local or_sex = 1.05
local or_score = 1.10

local b_sex = log(`or_sex')
local b_score = log(`or_score')
matrix input C = (1 0.20 -0.10 \ 0.20 1 -0.10 \ -0.10 -0.10 1)
matlist C
drawnorm age lat sexr, double corr(C) n(`n')
corr age lat sexr

generate byte score = 0
forvalues cut = 1/10 {
    quietly replace score = score + 1 if invnormal(`cut' / 11) < lat 
}

g sex = runiform() < logistic(sexr)

* Adults aged 18 to 65
quietly replace age = floor((65 - 18 + 1) * normal(age) + 18) 

*** Generate binary variables
corr2data r1 r2, corr(1 0.4 \ 0.4 1) means(`intercept6', `intercept12')

* Binary variables
sum r1 
local r6mean = r(mean)
sum r2
local r12mean = r(mean)

g temp = (logit(`r12mean') + ///
sex* `b_sex' + ///
score*`b_score')

generate recovery12= runiform() < logistic(temp)

capture drop temp
g temp = (logit(`r6mean') + ///
sex* `b_sex' + ///
score*`b_score')

generate recovery6= runiform() < logistic(temp)
drop temp

tab1 recovery6 recovery12
polychoric recovery6 recovery12

[Leonardo Guizzetti]

Didn’t you post this question on Reddit?

[Jen Ward]

Leonardo Guizzetti , I didn't. It is a follow-up from an earlier query on here re. simulating data, but I figured it deserved a separate thread. Does the Reddit post have a solution?

[Leonardo Guizzetti]

I don’t know now, the person deleted their thread.

[Bert Lloyd]

As a rough approximation, could you do something like

loc rho = 0.4
mat myCorr = [1, `rho' \ `rho' , 1]
drawnorm x1_star x2_star, corr(myCorr)
gen byte x1 = x1_star>0
gen byte x2 = x2_star>0

and then experiment with different values of rho until the correlation of x1 and x2 is close to what you want?

[Jen Ward]

Thanks Bert Lloyd , apologies if my question wasn't clear. The two binary variables are meant to have 35% and 20% probability; I also know that they are influenced by age and a Likert scale score - this is why I generate the temp variables in my code to be able to specify this. I am just not sure my approach is working well - any suggestion?

[Bert Lloyd]

I would try user-written rbinary, SJ article here:

Chen, M. (2015). Generating Nonnegatively Correlated Binary Random Variates. The Stata Journal, 15(1), 301-308. https://doi.org/10.1177/1536867X1501500118

[Joseph Coveney]

You can also use the user-written ovbd from SSC. It comes with an ancillary do-file that illustrates various applications in its use, including use cases of negatively correlated binomial data.

I illustrate it below for the particulars that you mention, namely, correlation coefficient of 0.4 for proportions of 0.35 and 0.20.

.ÿ
.ÿversionÿ18.0

.ÿ
.ÿclearÿ*

.ÿ
.ÿ//ÿseedem
.ÿsetÿseedÿ884193559

.ÿ
.ÿtempnameÿMeansÿCorr

.ÿmatrixÿinputÿ`Means'ÿ=ÿ(0.35ÿ0.20)

.ÿmatrixÿdefineÿ`Corr'ÿ=ÿJ(2,ÿ2,ÿ0.4)ÿ+ÿI(2)ÿ*ÿ0.6

.ÿovbdÿout1ÿout2,ÿmeans(`Means')ÿcorr(`Corr')ÿn(10000)ÿclear
(obsÿ10,000)

.ÿ
.ÿforeachÿvarÿofÿvarlistÿout?ÿ{
ÿÿ2.ÿÿÿÿÿsummarizeÿ`var',ÿmeanonly
ÿÿ3.ÿÿÿÿÿdisplayÿinÿsmclÿasÿtextÿ"`var':ÿ"ÿasÿresultÿ%04.2fÿr(mean)
ÿÿ4.ÿ}
out1:ÿ0.35
out2:ÿ0.20

.ÿ
.ÿquietlyÿcorrelateÿout?

.ÿdisplayÿinÿsmclÿasÿtextÿ"rhoÿ=ÿ"ÿasÿresultÿ%04.2fÿr(rho)
rhoÿ=ÿ0.39

.ÿ
.ÿgenerateÿintÿpidÿ=ÿ_n

.ÿquietlyÿreshapeÿlongÿout,ÿi(pid)ÿj(itm)

.ÿxtgeeÿoutÿi.itm,ÿi(pid)ÿt(itm)ÿfamily(binomial)ÿlink(logit)ÿcorr(unstructured)ÿnolog

GEEÿpopulation-averagedÿmodelÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿ=ÿ20,000
Groupÿandÿtimeÿvars:ÿpidÿitmÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿgroupsÿ=ÿ10,000
Family:ÿBinomialÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿObsÿperÿgroup:ÿÿ
Link:ÿÿÿLogitÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿminÿ=ÿÿÿÿÿÿ2
Correlation:ÿunstructuredÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿavgÿ=ÿÿÿÿ2.0
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿmaxÿ=ÿÿÿÿÿÿ2
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(1)ÿÿÿÿÿ=ÿ804.06
Scaleÿparameterÿ=ÿ1ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿ=ÿ0.0000

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿoutÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿ2.itmÿ|ÿÿ-.7263824ÿÿÿ.0256166ÿÿÿ-28.36ÿÿÿ0.000ÿÿÿÿÿÿ-.77659ÿÿÿ-.6761748
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-.6406483ÿÿÿ.0210349ÿÿÿ-30.46ÿÿÿ0.000ÿÿÿÿ-.6818759ÿÿÿ-.5994207
------------------------------------------------------------------------------

.ÿxtcorr

Estimatedÿwithin-pidÿcorrelationÿmatrixÿR:

ÿÿÿÿÿÿÿÿc1ÿÿÿÿÿÿc2
r1ÿÿ1.0000
r2ÿÿ0.3864ÿÿ1.0000

.ÿ
.ÿexit

endÿofÿdo-file


.

It does have a dependency on ridder available on StataCorp's website (search ridder).

You can render the proportions dependent upon predictor variables in the manner that you show above in #1.

[Jen Ward]

thanks Joseph Coveney - I had a go at editing your code so that the proportions can also be dependent on predictor variables; but I do get an error message "Root not found" as it seems -ovbd- does not seem to like the way I input the means as local variables. How can I overcome this?

I also noticed that ovbd creates a new dataset which only include out1 and out2; but I would like to keep the other variables age, sex and score as those are covariates in my model. I suppose if I use a seed that I can create the covariates and the binary variables in two separate files and then merge them. Is there an easier way?

clear *

// ssc install ovbd
// search ridder

set seed 6710789

local n = 1000
local intercept6 = 0.20
local intercept12 = 0.35
local or_sex = 1.05
local or_score = 1.10

local b_sex = log(`or_sex')
local b_score = log(`or_score')

matrix input C = (1 0.20 -0.10 \ 0.20 1 -0.10 \ -0.10 -0.10 1)
matlist C

drawnorm age lat sexr, double corr(C) n(`n')
corr age lat sexr

generate byte score = 0
forvalues cut = 1/10 {
    quietly replace score = score + 1 if invnormal(`cut' / 11) < lat
}

g sex = runiform() < logistic(sexr)

* Adults aged 18 to 65
quietly replace age = floor((65 - 18 + 1) * normal(age) + 18)



*** BINARY OUTCOMES

local r6mean = 0.35
local r12mean = 0.20

g temp1 = (logit(`r12mean') + ///
sex* `b_sex' + ///
score*`b_score')


g temp2 = (logit(`r6mean') + ///
sex* `b_sex' + ///
score*`b_score')

tempname Means Corr
sum temp1
local mean1 = r(mean) 

sum temp2
local mean2 = r(mean)


matrix input `Means' = (`mean1' `mean2')
matlist `Means'
matrix define `Corr' = J(2,2,0.4)+I(2)*0.6
matlist `Corr'
ovbd out1 out2, means(`Means') corr(`Corr') n(10000) clear

foreach var of varlist out? {
    sum `var', meanonly

    display in smcl as text "`var': " as result %04.2f r(mean)
}

corr out*

[Joseph Coveney]

. . . -ovbd- does not seem to like the way I input the means as local variables. How can I overcome this?

. . . ovbd creates a new dataset which only include out1 and out2 . . . Is there an easier way?

ovbd is intended for longitudinal data with semiarbitrary covariance structure. Because you have a simple covariance between pairs of binomial outcome values and you're using logistic regression, I now think it would be better in your case to use a more conventional random effects setup to model your correlated binary variables, sorry. You can specify the correlation between the six- and twelve-month binary outcomes in terms of an intraclass correlation coefficient. Threads about how to simulate such correlated binary outcomes have come up before, for example, here.

Using the approach shown there in that thread to first specify the correlation (see the first lines of code below that are shown in red), you can then specify the linear predictor for a logistic regression model that includes sex and score as predictors with given odds ratios (pertinent line of code also shown in red below). Because you're using simulated predictor values, it's simpler and more convenient to just leave the predictor values in their standard normal form that you get from drawnorm and that's what I've done below. The obtained odds ratios for the two predictors will vary somewhat but their confidence interval should usually include your population values of 1.05 and 1.10; likewise, the observed proportions will vary (the ones obtained below are spot on by chance), but should hover around 35% and 20% for six- and twelve-month time points. The intraclass correlation coefficient (rho in the regression table from xtlogit) should be close to the 0.40 specified—it's 0.39 for the seed below. Output from the regression model is shown below the code section, but I've also attached both the do-file and log-file if you're interested in looking at the output in detail or in pursuing this approach further.

version 18.0

clear *

// seedem
set seed 549607219

tempname Corr
matrix input `Corr' = (1 0.20 -0.10 \ 0.20 1 -0.10 \ -0.10 -0.10 1)
quietly drawnorm age sco sex, double corr(`Corr') n(1000)

generate int pid = _n

local icc 0.4
generate double pid_u = rnormal(0, sqrt(`icc' * _pi^2 / (3 - 3 *`icc')))

quietly expand 2
bysort pid: generate byte tim = _n * 6

generate double xbu = logit(cond(tim == 6, 0.35, 0.20)) + ///
    ln(1.05) * sex + ln(1.10) * sco + ///
        pid_u

generate byte out = rbinomial(1, invlogit(xbu))

xtlogit out c.(age sex sco) i.tim , i(pid) re or nolog
display in smcl as text %04.2f invlogit(_b[_cons]), ///
    %04.2f invlogit(_b[_cons] + _b[12.tim])

exit

.ÿxtlogitÿoutÿc.(ageÿsexÿsco)ÿi.timÿ,ÿi(pid)ÿreÿorÿnolog

Random-effectsÿlogisticÿregressionÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿ=ÿÿ2,000
Groupÿvariable:ÿpidÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿgroupsÿ=ÿÿ1,000

Randomÿeffectsÿu_iÿ~ÿGaussianÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿObsÿperÿgroup:
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿminÿ=ÿÿÿÿÿÿ2
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿavgÿ=ÿÿÿÿ2.0
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿmaxÿ=ÿÿÿÿÿÿ2

Integrationÿmethod:ÿmvaghermiteÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿIntegrationÿpts.ÿ=ÿÿÿÿÿ12

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(4)ÿÿÿÿÿ=ÿÿ46.89
Logÿlikelihoodÿ=ÿ-1208.8635ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿ=ÿ0.0000

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿoutÿ|ÿOddsÿratioÿÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿageÿ|ÿÿÿ1.158607ÿÿÿ.0899614ÿÿÿÿÿ1.90ÿÿÿ0.058ÿÿÿÿÿ.9950466ÿÿÿÿ1.349052
ÿÿÿÿÿÿÿÿÿsexÿ|ÿÿÿ1.132853ÿÿÿ.0851516ÿÿÿÿÿ1.66ÿÿÿ0.097ÿÿÿÿÿ.9776703ÿÿÿÿ1.312667
ÿÿÿÿÿÿÿÿÿscoÿ|ÿÿÿ1.157913ÿÿÿÿ.085396ÿÿÿÿÿ1.99ÿÿÿ0.047ÿÿÿÿÿ1.002074ÿÿÿÿ1.337987
ÿÿÿÿÿÿ12.timÿ|ÿÿÿ.4843162ÿÿÿ.0567754ÿÿÿÿ-6.18ÿÿÿ0.000ÿÿÿÿÿ.3848969ÿÿÿÿ.6094156
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿ.5274069ÿÿÿ.0499675ÿÿÿÿ-6.75ÿÿÿ0.000ÿÿÿÿÿ.4380275ÿÿÿÿ.6350242
-------------+----------------------------------------------------------------
ÿÿÿÿ/lnsig2uÿ|ÿÿÿ.7598952ÿÿÿ.1981203ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.3715866ÿÿÿÿ1.148204
-------------+----------------------------------------------------------------
ÿÿÿÿÿsigma_uÿ|ÿÿÿ1.462208ÿÿÿ.1448465ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ1.204173ÿÿÿÿ1.775535
ÿÿÿÿÿÿÿÿÿrhoÿ|ÿÿÿÿ.393899ÿÿÿ.0472997ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.3059206ÿÿÿÿ.4893407
------------------------------------------------------------------------------
Note:ÿEstimatesÿareÿtransformedÿonlyÿinÿtheÿfirstÿequationÿtoÿoddsÿratios.
Note:ÿ_consÿestimatesÿbaselineÿoddsÿ(conditionalÿonÿzeroÿrandomÿeffects).
LRÿtestÿofÿrho=0:ÿchibar2(01)ÿ=ÿ68.11ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>=ÿchibar2ÿ=ÿ0.000

.ÿdisplayÿinÿsmclÿasÿtextÿ%04.2fÿinvlogit(_b[_cons]),ÿ///
>ÿÿÿÿÿ%04.2fÿinvlogit(_b[_cons]ÿ+ÿ_b[12.tim])
0.35ÿ0.20

Because of its correlation with the two predictors that are, the odds ratio for age is also different from one (P = 0.06) even though it is not included in the linear prediction generation.

[Jen Ward]

Thanks Joseph Coveney , I will use this approach. I noticed you used

generate byte out = rbinomial(1, invlogit(xbu))

while I have been using

generate recovery12= runiform() < logistic(temp)

I am wondering what is the main difference and what's the advantage of using -rbinomial-; I looked at the help file but it didn't seem obvious to me how it works.

[Joseph Coveney]

I noticed you used

generate byte out = rbinomial(1, invlogit(xbu))

while I have been using

generate recovery12= runiform() < logistic(temp)

I am wondering what is the main difference and what's the advantage of using -rbinomial-;

In the Bernoulli case, they're essentially equivalent and nearly identical.

version 18.0

clear *

// seedem
set seed 1359979969

quietly set obs 10000

quietly {
    generate byte out1 = .
    generate byte out2 = .
    
}

foreach xb in -sqrt(2) 0 sqrt(2) {
    
    display in smcl as text _newline(2) "`xb' " ///
        as result invlogit(`xb') - logistic(`xb')

    local state = c(rngstate)
    quietly replace out1 = rbinomial(1, invlogit(`xb'))

    set rngstate `state'
    quietly replace out2 = runiform() < logistic(`xb')

    tabulate out?, matcell(M)
    display in smcl as text _newline(1) ///
        "Method 1: " as result %04.2f (M[2, 1] + M[2, 2]) / r(N)
    display in smcl as text ///
        "Method 2: " as result %04.2f (M[1, 2] + M[2, 2]) / r(N)
}

exit

-sqrt(2)ÿ-2.776e-17

ÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿout2
ÿÿÿÿÿÿout1ÿ|ÿÿÿÿÿÿÿÿÿ0ÿÿÿÿÿÿÿÿÿÿ1ÿ|ÿÿÿÿÿTotal
-----------+----------------------+----------
ÿÿÿÿÿÿÿÿÿ0ÿ|ÿÿÿÿÿ6,075ÿÿÿÿÿÿ1,964ÿ|ÿÿÿÿÿ8,039ÿ
ÿÿÿÿÿÿÿÿÿ1ÿ|ÿÿÿÿÿ1,961ÿÿÿÿÿÿÿÿÿÿ0ÿ|ÿÿÿÿÿ1,961ÿ
-----------+----------------------+----------
ÿÿÿÿÿTotalÿ|ÿÿÿÿÿ8,036ÿÿÿÿÿÿ1,964ÿ|ÿÿÿÿ10,000ÿ

Methodÿ1:ÿ0.20
Methodÿ2:ÿ0.20


0ÿ0

ÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿout2
ÿÿÿÿÿÿout1ÿ|ÿÿÿÿÿÿÿÿÿ0ÿÿÿÿÿÿÿÿÿÿ1ÿ|ÿÿÿÿÿTotal
-----------+----------------------+----------
ÿÿÿÿÿÿÿÿÿ0ÿ|ÿÿÿÿÿÿÿÿÿ0ÿÿÿÿÿÿ5,089ÿ|ÿÿÿÿÿ5,089ÿ
ÿÿÿÿÿÿÿÿÿ1ÿ|ÿÿÿÿÿ4,911ÿÿÿÿÿÿÿÿÿÿ0ÿ|ÿÿÿÿÿ4,911ÿ
-----------+----------------------+----------
ÿÿÿÿÿTotalÿ|ÿÿÿÿÿ4,911ÿÿÿÿÿÿ5,089ÿ|ÿÿÿÿ10,000ÿ

Methodÿ1:ÿ0.49
Methodÿ2:ÿ0.51


sqrt(2)ÿ0

ÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿout2
ÿÿÿÿÿÿout1ÿ|ÿÿÿÿÿÿÿÿÿ0ÿÿÿÿÿÿÿÿÿÿ1ÿ|ÿÿÿÿÿTotal
-----------+----------------------+----------
ÿÿÿÿÿÿÿÿÿ0ÿ|ÿÿÿÿÿ1,981ÿÿÿÿÿÿÿÿÿÿ0ÿ|ÿÿÿÿÿ1,981ÿ
ÿÿÿÿÿÿÿÿÿ1ÿ|ÿÿÿÿÿÿÿÿÿ0ÿÿÿÿÿÿ8,019ÿ|ÿÿÿÿÿ8,019ÿ
-----------+----------------------+----------
ÿÿÿÿÿTotalÿ|ÿÿÿÿÿ1,981ÿÿÿÿÿÿ8,019ÿ|ÿÿÿÿ10,000ÿ

Methodÿ1:ÿ0.80
Methodÿ2:ÿ0.80

.ÿ
.ÿexit

endÿofÿdo-file


.

The advantage of the former is when you want to simulate multiple-trial binomial outcomes.