Generate random data for logistic regression with a categorical independent variable

Question

I am trying to generate a data frame of fake data for exploratory purposes. Specifically, I am trying to produce data with a binary dependent variable (say, failure/success), and a categorical independent variable called 'picture' with 5 levels (pict1, pict2, etc.). I am following the answer provided here, which allows me to successfully generate the data. However, I need each level of 'picture' to occur the same number of times (i.e. 11 repetitions of each level = 55 total observations per subject).

Here is a reproducible example of what has worked so far (code from user: ocram):

library(dummies)

#------ 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("pict1","pict2", "pict3", "pict4", "pict5"), 
              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(picture=x, choice=y)

  #fit the logistic model
  mod <- glm(choice ~ picture, 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]
}

However, as you can see from the output, each level of the factor 'picture' does not occur the same number of times (i.e. 200 times each).

> summary(data)
picture     choice     
pict1:200   Min.   :0.000  
pict2:207   1st Qu.:0.000  
pict3:217   Median :1.000  
pict4:163   Mean   :0.559  
pict5:213   3rd Qu.:1.000  
            Max.   :1.000 

Moreover, it is not entirely clear to me how to manipulate the initial beta values as to determine the probability of success/failure for each level of 'picture'. I cannot comment the original question because I do not yet have the necessary reputation points.

Answer 1

1. If you want 200 copies of each of 5 levels in a polytomous variable in random order, then do this instead:

   x <- sample(rep(paste0('pict', 1:5), 200))
   

2. If you want to control for overall prevalence of a specific outcome, then you must choose which beta you will fudge. I usually do beta0.

   MM         <- model.matrix(~x)
   betas      <- rnorm(4)
   prevTarget <- 0.3
   prevDiff   <- function(beta0)  prevTarget - 
                                  mean(binomial()$linkinv(MM%*%c(beta0, betas)))
   beta0      <- uniroot(prevDiff, c(-100, 100))$root
   mean(binomial()$linkinv(MM%*%c(beta0, betas)))
   

Answer 2

I eventually found an acceptable answer to my question. This may or may not be the best/most sophisticated way to handle this, but it enabled me to get what I was after.

set.seed(1)

intercept1 = -6.0
beta1      = 2.75
xtest1     = rnorm(1300, 1, 3)
linpred1   = intercept1 + xtest1 * beta1
prob1      = exp(linpred1)/(1 + exp(linpred1))
runis1     = runif(1300, 0, 1)
ytest1     = ifelse(runis1 < prob1, 1, 0)

intercept2 = -7.0
beta2      = 2.75
xtest2     = rnorm(1300, 1, 3)
linpred2   = intercept2 + xtest2 * beta2
prob2      = exp(linpred2)/(1 + exp(linpred2))
runis2     = runif(1300, 0, 1)
ytest2     = ifelse(runis2 < prob2, 1, 0)

intercept3 = -7.0
beta3      = 2.75
xtest3     = rnorm(1300, 1, 3)
linpred3   = intercept3 + xtest3 * beta3
prob3      = exp(linpred3)/(1 + exp(linpred3))
runis3     = runif(1300, 0, 1)
ytest3     = ifelse(runis3 < prob3, 1, 0)


newdf1 <- data.frame(prop = ytest1, x = xtest1, week = "w1")
newdf1 <- newdf1[with(newdf1, order(week, x)), ]
newdf1$stim <- rep(seq(-60, 60, by = 10), each = 100)

newdf2 <- data.frame(prop = ytest2, x = xtest2, week = "w3")
newdf2 <- newdf2[with(newdf2, order(week, x)), ]
newdf2$stim <- rep(seq(-60, 60, by = 10), each = 100)

newdf3 <- data.frame(prop = ytest3, x = xtest3, week = "w6")
newdf3 <- newdf3[with(newdf3, order(week, x)), ]
newdf3$stim <- rep(seq(-60, 60, by = 10), each = 100)

temp <- rbind(newdf1, newdf2)
df <- rbind(temp, newdf3)

I am now able to manipulate the data frame, perform logistic regression, and produce plots. Thank you to everybody who provided input.