# Simulating data for logistic regression with a categorical variable

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?

• sample(x=c(1, 2, 3), size=1, prob=rep(1/3, 3)) generates one of "1", "2", or "3" with equal probability. Commented Feb 13, 2013 at 17:42
• 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' Commented Feb 13, 2013 at 21:47

## 1 Answer

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
> #---------------------

• @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? Commented Aug 12, 2014 at 18:10
• @dchandler: The parameters and probabilities were chosen arbitrarily, for the sake of illustration. Commented Aug 12, 2014 at 18:16
• @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. Commented Aug 13, 2014 at 16:06