7
$\begingroup$

Suppose I want to simulate a survey variable in R which values derive from four binary cases $c_i$ with $i=\{1,2,3, 4\}$, each with it's probability $\Pr(c_i=1)=p_i$. Let's assume that the cases should correlate pairwise with each other in this (idealized) structure:

$\begin{matrix}& c1 & c2 & c3 & c4\\\ c1 & 1 & -1 & 0 & 0\\\ c2 & -1 & 1 & 0 & 0\\\ c3 & 0 & 0 & 1 & -1\\\ c4 & 0 & 0 & -1 & 1 \end{matrix}$

For the simulation I want to take $n$ draws (i. e. a sample with size $n$) from the correlation structure above, where in each draw the number of $c_i=1$ are added. This should result in something like $X=(0, 2, 3, 0, 2, 0, 4, 1, 2 ,0 , \dots)$.

My attempt in R coding so far looks like the following. (Thereby I orientated myself on this tutorial, unfortunately it doesn't fit to my needs until the end.)

set.seed(961)
mu <- rep(0, 4)
Sigma <- matrix(c(1, -1, 0, 0, 
                  -1, 1, 0, 0, 
                  0, 0, 1, -1, 
                  0, 0, -1, 1), nrow = 4, ncol = 4)
Sigma
#      [,1] [,2] [,3] [,4]
# [1,]    1   -1    0    0
# [2,]   -1    1    0    0
# [3,]    0    0    1   -1
# [4,]    0    0   -1    1

library(MASS)
rawvars <- mvrnorm(n = 1e4, mu = mu, Sigma = Sigma)

# cov(rawvars)

cor(rawvars)
#             [,1]         [,2]         [,3]         [,4]
# [1,]  1.000000000 -1.000000000 -0.006839596  0.006839597
# [2,] -1.000000000  1.000000000  0.006839597 -0.006839598
# [3,] -0.006839596  0.006839597  1.000000000 -1.000000000
# [4,]  0.006839597 -0.006839598 -1.000000000  1.000000000

pvars <- pnorm(rawvars)

Until here I think it looks good and it seems I am on the right way. In the following the author draws Poisson, exponential, and other data and of all things his code seems to be flawed in the binary example.

I made an attempt myself but I could not specify the probabilities $p_i$ so as a workaround I chose a value $p=.3$ which seems OK by rule of thumb, and it looks like this:

binvars <- qbinom(pvars, 4, .3)

I then get following distribution:

summary(as.factor(binvars))
#    0     1     2     3     4 
# 9704 16283 10676  2994   343

hist(binvars)

enter image description here

Now I'm facing three major problems. First, how should the resulting vector resp. the distribution look like in general (e. g. with respect to the zeroes)? Second, does the attempt so far make sense at all? Third, how could I solve this problem to the end?

Any help is very appreciated.

Improve this question
$\endgroup$
4
  • $\begingroup$ I understand what you are doing, but I don't understand your goal. Do you want to generate 4 correlated bernulli or 1 binomial with 4 trials? because it seems to me that you hae done the latter $\endgroup$
    – gioxc88
    Jan 28, 2018 at 8:59
  • $\begingroup$ I want to generate one resulting variable $X=(0,2,3,0,2,0,4,1,2,0,…)$ which is derived from four binaries $(0, 1)$, from which two pairs are negatively correlated somehow (sure not by exactly -1!) and together they add up to $X$. E.g. $X11=1, X12=0, X13=0, X14=1, X21=0, X22=0, X23=1, X24=0$; result $X = \{X11, X21\} = \{2, 1\}$. $\endgroup$
    – jay.sf
    Jan 28, 2018 at 9:10
  • $\begingroup$ @HardCore Actually there will be a second step: a second variable $Z$ a fifth binary $\pi$–four binaries of $Z$ are pairwise positively correlated to $X$, $\pi$ can be independent. Finally by subtraction of the means $\pi$ is to be estimated (because it is a simulation of a case where $\pi$ is unknown). $\endgroup$
    – jay.sf
    Jan 28, 2018 at 9:29
  • 1
    $\begingroup$ I think you are on the right track. check my answer to see if it makes sense! $\endgroup$
    – gioxc88
    Jan 28, 2018 at 10:15

4 Answers 4

Reset to default
6
+50
$\begingroup$

I can see the problem in you experiment.

You did everything right until:

binvars <- qbinom(pvars, 4, .3)

Why? Let's try to understand what pvars represents first.

pvars is a matrix that contains Uniforms between $0$ and $1$ that have the correlation structure you specified earlier in the variable sigma.

If you feed those Uniforms to any desired inverse CDF, say Bernulli, you get 4 vectors of correlated Bernulli.

This is an application of a GAUSSIAN COPULA.

The problem with your code is that you feed the pvars to the inverse CDF of a Binomial distribution with 4 trials, $Bin(n,p) = Bin(4, 0.3)$ in your case.

You simulated 4 correlated Binomial(4, 0.3) not 4 correlated Bernulli(0.3)

Think about it. Given the correlation matrix you chose how is it possible that the sum of your 4 Bernulli gives 4? It is not possible, nevertheless you obtained some 4 in your summary.

Change the code like this:

binvars <- qbinom(pvars, 1, .3)
X = apply(binvars,1,sum)

X is the vector you are looking for.

EDIT:

Following @Xi'an answer it's much simpler:

Basically you simulate only two independent bernoulli say b1 and b3 with respective probability of success p1 and p3.

Finally you set the remaining two bernoulli b2 = 1 - b1 and b4 = 1 - b3, ans @Xi'an suggested.

This way b2 and b4 end up being perfectly negatively correlated to b1 and b3.

set.seed(961)    

p1 = 0.3
p3 = 0.7 

b1 =rbinom(1e4, 1, p1)
b3 =rbinom(1e4, 1, p3)
b2 = 1 - b1
b4 = 1 - b3

binvars = cbind(b1,b2,b3,b4)

X = apply(binvars,1,sum)

cor(binvars)

          b1            b2            b3            b4
b1  1.0000000000 -1.0000000000  0.0007217015 -0.0007217015
b2 -1.0000000000  1.0000000000 -0.0007217015  0.0007217015
b3  0.0007217015 -0.0007217015  1.0000000000 -1.0000000000
b4 -0.0007217015  0.0007217015 -1.0000000000  1.0000000000
Improve this answer
$\endgroup$
6
  • $\begingroup$ This perfectly makes sense! I don't understand what probability $p$ I define in qbinom(pvars, 1, p) and which value would be appropriate. $\endgroup$
    – jay.sf
    Jan 28, 2018 at 11:53
  • $\begingroup$ you define the probability of success that is the probability of obtaining 1 in a bernulli trial $\endgroup$
    – gioxc88
    Jan 28, 2018 at 12:11
  • $\begingroup$ For sure. Since there are at least two probabilities for each bernoulli pairs, and the value of $p$ significantly shifts the vertex of the PDF, what is the meaning of this probability in the simulation, and what is the right value I want? $\endgroup$
    – jay.sf
    Jan 28, 2018 at 12:20
  • $\begingroup$ I edited my answer. Check it out $\endgroup$
    – gioxc88
    Jan 28, 2018 at 21:14
  • $\begingroup$ I like the first part of your answer much more because I talked about survey data ro be simulated, and these are never perfectly negatively correlated. This is what I meant by "idealized structure". Sorry I wasn't more clear in this point. The second part of your answer is a very special case, ok. Now we don't have yet a solution for the case where $\mathrm{Cor}(b_1,\ b_2)=[.6; .9]$ something ($\mathrm{Cor}(b_2,\ b_3)$. accordingly). $\endgroup$
    – jay.sf
    Jan 29, 2018 at 7:39
10
$\begingroup$

With the values of the correlation matrix that are proposed in the question, no simulation is needed to solve the question, as the answer is deterministic and available. Here is the reason why:

To generate two Bernoulli variates that are perfectly correlated, i.e., when$$\mathrm{corr}(C_1,C_2)=-1$$ one needs to find the conditional distribution that fits this constraint: \begin{align*} \mathrm{cov}(C_1,C_2) &=\mathbb{E}[C_1C_2]-\mathbb{E}[C_1]\mathbb{E}[C_2]\\ &=\mathbb{P}(C_1=C_2=1)-\mathbb{P}(C_1=1)\mathbb{P}(C_2=1)\\ &=\mathbb{P}(C_1=1|C_2=1)\mathbb{P}(C_2=1)-\mathbb{P}(C_1=1)\mathbb{P}(C_2=1)\\ &=q_{11}p_2-p_1p_2\\ &=p_2(q_{11}-p_1) \end{align*} which should satisfy$$p_2(q_{11}-p_1)=-\sqrt{p_1p_2(1-p_1)(1-p_2)}$$or$$q_{11}=p_1-\sqrt{p_1(1-p_1)(1-p_2)/p_2}=p_1\left\{1- \sqrt{(1-p_1)(1-p_2)/p_1p_2}\right\}$$which allows for a solution in $(0,1)$ if and only if$$(1-p_1)(1-p_2)\le p_1p_2$$To generate $C_1$ when $C_2=0$, one needs the conditional probability $q_{10}=\mathbb{P}(C_1=1|C_2=0)$ which is given by$$q_{10}(1-p_2)+q_{11}p_2=p_1$$or $$q_{10}=(p_1-p_2q_{11})/(1-p_2)=p_1+\sqrt{p_1p_2(1-p_1)/(1-p_2)}$$which allows for a solution in $(0,1)$ only if $$\sqrt{\frac{p_2(1-p_1)}{(1-p_2)p_1}}\le\frac{1}{p_1}-1=\frac{1-p_1}{p_1}$$ i.e.,$$\sqrt{\frac{p_2}{(1-p_2)}}\le\sqrt{\frac{1-p_1}{p_1}}$$which amounts to$$p_1p_2\le (1-p_1)(1-p_2)$$Therefore, the only case when a negative correlation of $-1$ is feasible is when$$p_1p_2 = (1-p_1)(1-p_2)$$or$$\frac{p_1}{1-p_1}=\frac{1-p_2}{p_2}$$i.e.,$$p_2=1-p_1$$ This leads to $$q_{11}=p_1-\sqrt{p_1(1-p_1)(1-p_2)/p_2}=p_1-\sqrt{p_1(1-p_1)p_1/(1-p_1)}=0$$and$$q_{10}=p_1+\sqrt{p_1p_2(1-p_1)/(1-p_2)}=p_1+\sqrt{p_1(1-p_1)(1-p_1)/p_1}=1$$meaning that $C_1$ is equal to zero when $C_2$ is equal to one and vice-versa. Therefore$$C_1+C_2=1$$with probability $1$.

QED: no simulation is needed!

Note that this property would extend to the Binomial case in that $C_2=n-C_1$ is the only Binomial $\mathcal{B}(n,1-p)$ perfectly and negatively correlated with $C_2\sim\mathcal{B}(n,p)$

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ great answer. I would only suggest to correct $q_{10}$ which accidentally you wrote as $q_{10}=\mathbb{P}(C_1=1|C_2=1)$ $\endgroup$
    – gioxc88
    Jan 28, 2018 at 17:23
4
$\begingroup$

Using a mathematical method described by whuber in this related question, I have programmed a function that generates pairs of correlated binomial random variables using the standard syntax for distributions in R. You can call this function to generate any desired number of correlated Bernoulli random variables, with specified probabilities prob1 and prob1 and specified corelation corr. Note that the correlation coefficient is the correlation of the individual Bernoulli values that sum to the binomial, not the correlation between the binomial values themselves.

rcorrbinom <- function(n, size = 1, prob1, prob2, corr = 0) {
  
  #Check inputs
  if (!is.numeric(n))             { stop('Error: n must be numeric') }
  if (length(n) != 1)             { stop('Error: n must be a single number') }
  if (as.integer(n) != n)         { stop('Error: n must be a positive integer') }
  if (n < 1)                      { stop('Error: n must be a positive integer') }
  if (!is.numeric(size))          { stop('Error: n must be numeric') }
  if (length(size) != 1)          { stop('Error: n must be a single number') }
  if (as.integer(size) != size)   { stop('Error: n must be a positive integer') }
  if (size < 1)                   { stop('Error: n must be a positive integer') }
  if (!is.numeric(prob1))         { stop('Error: prob1 must be numeric') }
  if (length(prob1) != 1)         { stop('Error: prob1 must be a single number') }
  if (prob1 < 0)                  { stop('Error: prob1 must be between 0 and 1') }
  if (prob1 > 1)                  { stop('Error: prob1 must be between 0 and 1') }
  if (!is.numeric(prob2))         { stop('Error: prob2 must be numeric') }
  if (length(prob2) != 1)         { stop('Error: prob2 must be a single number') }
  if (prob2 < 0)                  { stop('Error: prob2 must be between 0 and 1') }
  if (prob2 > 1)                  { stop('Error: prob2 must be between 0 and 1') }
  if (!is.numeric(corr))          { stop('Error: corr must be numeric') }
  if (length(corr) != 1)          { stop('Error: corr must be a single number') }
  if (corr < -1)                  { stop('Error: corr must be between -1 and 1') }
  if (corr > 1)                   { stop('Error: corr must be between -1 and 1') }
  
  #Compute probabilities
  P00   <- (1-prob1)*(1-prob2) + corr*sqrt(prob1*prob2*(1-prob1)*(1-prob2))
  P01   <- 1 - prob1 - P00
  P10   <- 1 - prob2 - P00
  P11   <- P00 + prob1 + prob2 - 1
  PROBS <- c(P00, P01, P10, P11)
  if (min(PROBS) < 0)       { stop('Error: corr is not in the allowable range') }
  
  #Generate the output
  RAND <- array(sample.int(4, size = n*size, replace = TRUE, prob = PROBS),
                dim = c(n, size))
  VALS <- array(0, dim = c(2, n, size))
  OUT  <- array(0, dim = c(2, n))
  for (i in 1:n)    { 
  for (j in 1:size) { 
    VALS[1,i,j] <- (RAND[i,j] %in% c(3, 4))
    VALS[2,i,j] <- (RAND[i,j] %in% c(2, 4)) } 
    OUT[1, i]   <- sum(VALS[1,i,])
    OUT[2, i]   <- sum(VALS[2,i,]) }
  
  #Give output
  OUT }

Here is an example of using this function to produce a sample array containing a large number of correlated Bernoulli random variables. We can confirm that, for a large sample, the sampled values have sample means and sample correlation that is close to the specified parameters.

#Set parameters
n     <- 10^6
PROB1 <- 0.3
PROB2 <- 0.8
CORR  <- 0.2

#Generate sample of correlated Bernoulli random variables
set.seed(1)
SAMPLE <- rcorrbinom(n = n, prob1 = PROB1, prob2 = PROB2, corr = CORR)

#Check the properties of the sample
str(SAMPLE)
 num [1:2, 1:10000] 0 1 0 1 1 1 0 0 0 1 ...

mean(SAMPLE[1,])
[1] 0.300122

mean(SAMPLE[2,])
[1] 0.800145

cor(SAMPLE[1,], SAMPLE[2,])
[1] 0.20018
Improve this answer
$\endgroup$
0
$\begingroup$

Here is another code, specifically for generating 2 correlated Bernoulli variables following the method by @whuber:

  1. Given p, q and rho, determine the probabilities of (0,0), (0,1), (1,0) and (1,2)
  2. Then sample the four combinations with probabilities determined in 1
sim_bernoulli <- function(rho=-4/5, p=1/3, q=3/4, N=10){
  
  ## get prob for each combination x-y
  p_00 <- rho * sqrt(p*q*(1-p)*(1-q)) + (1-p)*(1-q)
  prob <- c(`(0,0)`=p_00, `(1,0)`=1-q-p_00, `(0,1)`=1-p-p_00, `(1,1)`=p_00+p+q-1)
  if (min(prob) < 0) {
    print(prob)
    stop("Error: a probability is negative.")
  }

  ## now simulate each combination x-y from prob above
  u <- sample.int(4, N, replace=TRUE, prob=prob)
  
  ## now get coresponding x-y combos
  XY <- cbind(c(0,1,0,1), c(0,0,1,1))[u,]
  XY
}

## example
set.seed(123)
out <- sim_bernoulli(N=10000)
c(colMeans(out), p=1/3, q=3/4)
#>                             p         q 
#> 0.3287000 0.7501000 0.3333333 0.7500000
c(cor(out)[2,1], rho = -4/5)
#>                   rho 
#> -0.8101124 -0.8000000

Created on 2023-09-29 with reprex v2.0.2

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.