0
$\begingroup$

Can someone tell me what is the best way to simulate a dataset with a binary target?

I understand the way in which a dataset can be simulated but what I'm looking for is to determine 'a-priori' the proportion of each class. What I thought was to change the intercept to achieve it but I couldn't do it and I don't know why. I guess because the average is playing a trick on me.

    set.seed(666)
    x1 = rnorm(1000)       
    x2 = rnorm(1000)

    p=0.25 # <<< I'm looking for a 25%/75%
    mean_z=log(p/(1-p))

    b0 = mean( mean_z - (4*x1 + 3*x2)) # = mean_z - mean( 2*x1 + 3*x2)
    z = b0 + 4*x1 + 3*x2  # = mean_z - (4*x1 + 3*x2)  + (4*x1 + 3*x2) = rep(mean_z,1000)
    mean( b0 + 4*x1 + 3*x2 ) == mean_z # TRUE!!

    pr = 1/(1+exp(-z))     
    y = rbinom(1000,1,pr)  
    mean(pr)  # ~ 40% << not achieved
    table(y) # 0:759 - 1:421

What I'm looking for is to simulate the typical "logistic" problem in which the binary target can be modeled as a linear combination of features.

These 'logistic' models assume that the log-odd ratio of the binary variable behaves linearly. That means:

log (p / (1-p)) = z = b0 + b1 * x1 + b2 * x2 where p = prob (y = 1)

Going back to my sample code, we could do, for example: z = 1.3 + 4 * x1 + 2 * x2 , but the probability of the class would be a result. Or instead we could choose coefficient b0 such that the probability is (statistically) similar to the one sought :

log (0.25 / 0.75) = b0 + 4 * x1 + 2 * x2

This is my approach, but there may be betters

$\endgroup$
5
  • 1
    $\begingroup$ Oh no ! The devil's number ! $\endgroup$
    – MysteryGuy
    Jul 30, 2018 at 13:40
  • $\begingroup$ maybe that's why! $\endgroup$
    – c0chi
    Jul 30, 2018 at 13:43
  • $\begingroup$ This seems very complex. Are you just trying to generate a data set with a given proportion of 0 & 1? That can be achieve much more simply with S = sample(0:1, 1000, replace=TRUE, prob=c(0.25, 0.75)) If that is not what you are looking for, please make it clear what your desired output would be. $\endgroup$
    – G5W
    Jul 30, 2018 at 23:57
  • $\begingroup$ Thanks @G5W. What I'm looking for is to simulate the typical "logistic" problem in which the binary target can be modeled as a linear combination of features. $\endgroup$
    – c0chi
    Jul 31, 2018 at 10:10
  • $\begingroup$ Actually, these 'logistic' models assume that the log-odd ratio of the binary variable behaves linearly. That means: log (p / (1-p)) = z = b0 + b1 * x1 + b2 * x2 where p = prob (y = 1). Going back to my sample code, we could do, for example: z = 1.3 + 4 * x1 + 2 * x2 , but the probability of the class would be a result. Or instead we could choose coefficient b0 such that the probability is (statistically) similar to the one sought :log (0.25 / 0.75) = b0 + 4 * x1 + 2 * x2. This is my approach, but there may be betters $\endgroup$
    – c0chi
    Jul 31, 2018 at 10:20

1 Answer 1

1
$\begingroup$

First suggestion is the coefficients you have used are very large. They correspond to odds ratios of 55 and 20. Unless you have specific reasons for doing this, $\ln(4)$ and $\ln(3)$ are probably more appropriate arbitrary choices.

The problem with your data simulation is you are only specifying the predicted probability when $x_1$ and $x_2$ are zero. This is easy to verify. Run your logistic regression model, and apply plogis() (inverse logit function) to the intercept coefficient. You will arrive at .25.

The issue is the transformation from logits to probabilities is non-linear. In your example, $$.25 = (1+e^{-\bar{z}})^{-1} \neq \mathrm{E}\big[(1+e^{-z})^{-1}\big]$$ as you observed.

I might approach your example problem using the latent variable approach. In this model, we assume:

\begin{equation} y^* = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \epsilon \quad\mathrm{and}\quad y = \begin{cases} 1, & \mathrm{if}\ y^* > 0\\ 0, & \mathrm{otherwise} \end{cases} \end{equation}

where $\epsilon \sim L(0, 1)$, $x_1\sim N(0,1)$, and $x_2\sim N(0,1)$ keeping with your approach.

With this approach, we can readily manipulate $\mathrm{E}[y]$ by changing the intercept. So we can re-write the equation as:

\begin{equation} y = \begin{cases} 1, & \mathrm{if}\ \beta_1 x_1 + \beta_2 x_2 + \epsilon > -\beta_0\\ 0, & \mathrm{otherwise} \end{cases} \end{equation}

Since the equation is additive on the logit scale, we can assume the quantity, $\beta_1 x_1 + \beta_2 x_2 + \epsilon$, converges to normality. The variance of such a distribution would be: $\beta_1^2 \times 1^2 + \beta_2^2 \times 1^2 + \pi^2/3$, $\pi^2/3$ is the variance of $L(0, 1)$.

If we want the proportion of 1's to be 25%, we can set $-\beta_0$ to the 75th percentile of a normal distribution with the variance above. The mean of this distribution would be: $\beta_1\times \bar{x}_1 + \beta_2 \times \bar{x}_2$. In our example, that would be 0. So $-\beta_0$ would be the 75% percentile of $N(0, \beta_1^2 \times 1^2 + \beta_2^2 \times 1^2 + \pi^2/3)$. Following my own recommendation, I use $\ln(4)$ and $\ln(3)$ as $\beta_1$ and $\beta_2$.

In R syntax, this would be:

set.seed(666)
b0 <- qnorm(.25, sd = sqrt(log(4) ^ 2 + log(3) ^ 2 + pi ^ 2 / 3))
res <- t(replicate(2500, {
  x1 <- rnorm(1000)
  x2 <- rnorm(1000)
  z <- b0 + log(4) * x1 + log(3) * x2
  y <- rbinom(1000, 1, plogis(z))
  c(mean(y), unname(coef(glm(y ~ x1 + x2, binomial))))
}))
# The first element is the mean of y
# Others are regression coefficients
colMeans(res)
# [1]  0.244050 -1.718929  1.397497  1.109592

Here's a more complicated example to see if it still works:

\begin{equation} y = \begin{cases} 1, & \mathrm{if}\ x_1 + 1.5 \times x_2 + .25 \times x_3 + \epsilon > -\beta_0\\ 0, & \mathrm{otherwise} \end{cases} \end{equation}

where $x_1\sim N(1, 1)$, $x_2 \sim N(0, 1)$ and $x_3 \sim \mathrm{Bern}(.5)$.

In this example, I want the mean of $y$ to be 10%. So $-\beta_0$ would be the 90% percentile of $N(1 + .25 \times .5, 1 + 1.5^2 + .25 \times .5 ^ 2 + \pi^2/3)$:

b0 <- -qnorm(.9, 1 + .25 * .5, sqrt(1 + 1.5 ^ 2 + .25 ^ 2 * .5 * .5 + pi ^ 2 / 3))
res <- t(replicate(2500, {
  x1 <- rnorm(1000, 1)
  x2 <- rnorm(1000)
  x3 <- rbinom(1000, 1, .5)
  z <- b0 + x1 + 1.5 * x2 + .25 * x3
  y <- rbinom(1000, 1, plogis(z))
  c(mean(y), unname(coef(glm(y ~ x1 + x2 + x3, binomial))))
}))
# The first element is the mean of y
# Others are regression coefficients
colMeans(res)
# [1]  0.0964992 -4.4560205  1.0125437  1.5117373  0.2611867

I hope this helps. I have tested it a couple of times and it appears to do a decent job. This is one approach I attempted that worked not so badly.


EDIT: The normal approximation requires that there are many additive features with none having an inordinate influence on the latent variable relative to others. If you simulate a large number of variables, then it works relatively well.

set.seed(12345)
mean(replicate(10000, {
  X <- matrix(rnorm(10000), 100)
  target <- qnorm(.25, 0, sqrt(ncol(X) + pi ^ 2 / 3))
  mean(y <- rbinom(nrow(X), 1, plogis(target + X %*% rep(1, ncol(X)))))
}))
# [1] 0.249637
$\endgroup$
0

Your Answer

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

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