9
$\begingroup$

I was trying to create some test data for logistic regression and I found this post How to simulate artificial data for logistic regression?

It is a nice answer but it creates only continuous variables. What about a categorical variable x3 with 5 levels (A B C D E) associated with y for the same example as in the link?

$\endgroup$
  • $\begingroup$ sample(x=c(1, 2, 3), size=1, prob=rep(1/3, 3)) generates one of "1", "2", or "3" with equal probability. $\endgroup$ – ocram Feb 13 '13 at 17:42
  • $\begingroup$ thanks for your comment, but how do I associate the probabilities here with the y of the post I mentioned? I copy some code from that post 'code' > set.seed(666) > x1 = rnorm(1000) # some continuous variables > x2 = rnorm(1000) > z = 1 + 2*x1 + 3*x2 # linear combination with a bias > pr = 1/(1+exp(-z)) # pass through an inv-logit function > y = rbinom(1000,1,pr) # bernoulli response variable 'code' $\endgroup$ – user1301295 Feb 13 '13 at 21:47
13
$\begingroup$

The model

Let $x_B = 1$ if one has category "B", and $x_B = 0$ otherwise. Define $x_C$, $x_D$, and $x_E$ similary. If $x_B = x_C = x_D = x_E = 0$, then we have category "A" (i.e., "A" is the reference level). Your model can then be written as

$$ \textrm{logit}(\pi) = \beta_0 + \beta_B x_B + \beta_C x_C + \beta_D x_D + \beta_E x_E $$ with $\beta_0$ an intercept.


Data generation in R

(a)

x <- sample(x=c("A","B", "C", "D", "E"), 
              size=n, replace=TRUE, prob=rep(1/5, 5))

The x vector has n components (one for each individual). Each component is either "A", "B", "C", "D", or "E". Each of "A", "B", "C", "D", and "E" is equally likely.

(b)

library(dummies)
dummy(x)

dummy(x) is a matrix with n rows (one for each individual) and 5 columns corresponding to $x_A$, $x_B$, $x_C$, $x_D$, and $x_E$. The linear predictors (one for each individual) can then be written as

linpred <- cbind(1, dummy(x)[, -1]) %*% c(beta0, betaB, betaC, betaD, betaE)

(c)

The probabilities of success follows from the logistic model:

pi <- exp(linpred) / (1 + exp(linpred))

(d)

Now we can generate the binary response variable. The $i$th response comes from a binomial random variable $\textrm{Bin}(n, p)$ with $n = 1$ and $p =$ pi[i]:

y <- rbinom(n=n, size=1, prob=pi)

Some quick simulations to check this is OK

> #------ parameters ------
> n <- 1000 
> beta0 <- 0.07
> betaB <- 0.1
> betaC <- -0.15
> betaD <- -0.03
> betaE <- 0.9
> #------------------------
> 
> #------ initialisation ------
> beta0Hat <- rep(NA, 1000)
> betaBHat <- rep(NA, 1000)
> betaCHat <- rep(NA, 1000)
> betaDHat <- rep(NA, 1000)
> betaEHat <- rep(NA, 1000)
> #----------------------------
> 
> #------ simulations ------
> for(i in 1:1000)
+ {
+   #data generation
+   x <- sample(x=c("A","B", "C", "D", "E"), 
+               size=n, replace=TRUE, prob=rep(1/5, 5))  #(a)
+   linpred <- cbind(1, dummy(x)[, -1]) %*% c(beta0, betaB, betaC, betaD, betaE)  #(b)
+   pi <- exp(linpred) / (1 + exp(linpred))  #(c)
+   y <- rbinom(n=n, size=1, prob=pi)  #(d)
+   data <- data.frame(x=x, y=y)
+   
+   #fit the logistic model
+   mod <- glm(y ~ x, family="binomial", data=data)
+   
+   #save the estimates
+   beta0Hat[i] <- mod$coef[1]
+   betaBHat[i] <- mod$coef[2]
+   betaCHat[i] <- mod$coef[3]
+   betaDHat[i] <- mod$coef[4]
+   betaEHat[i] <- mod$coef[5]
+ }
> #-------------------------
> 
> #------ results ------
> round(c(beta0=mean(beta0Hat), 
+         betaB=mean(betaBHat), 
+         betaC=mean(betaCHat), 
+         betaD=mean(betaDHat), 
+         betaE=mean(betaEHat)), 3)
 beta0  betaB  betaC  betaD  betaE 
 0.066  0.100 -0.152 -0.026  0.908 
> #---------------------
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ @ocram -- could give some intuition for good choices of parameters and choice of probabilities of components (part a)? How would changes to these affect the exercise? $\endgroup$ – d_a_c321 Aug 12 '14 at 18:10
  • $\begingroup$ @dchandler: The parameters and probabilities were chosen arbitrarily, for the sake of illustration. $\endgroup$ – ocram Aug 12 '14 at 18:16
  • 2
    $\begingroup$ @ocram -- understood. However, I'm looking for intuition on what would be good coefficients so that I could run more extensive simulations. For instance, if I wanted to simulate lasso regressions, I may be interested in adding meaningless variables (w/ zero coefficients) and seeing how the # of meaningless variables and the magnitude of the non-zero coefficients on meaningful variables affect the simulation. $\endgroup$ – d_a_c321 Aug 13 '14 at 16:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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