Simulation de données pour la régression logistique avec une variable catégorielle

J'essayais de créer des données de test pour la régression logistique et j'ai trouvé ce post Comment simuler des données artificielles pour la régression logistique?

C'est une bonne réponse mais elle ne crée que des variables continues. Qu'en est-il d'une variable catégorielle x3 à 5 niveaux (ABCDE) associée à y pour le même exemple que dans le lien?

Le modèle

$x_B = 1$$x_B = 0$$x_C$$x_D$$x_E$$x_B = x_C = x_D = x_E = 0$

Génération de données dans R

(une)

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

Le xvecteur a des ncomposants (un pour chaque individu). Chaque composant est soit "A", "B", "C", "D" ou "E". Chacun de "A", "B", "C", "D" et "E" est également probable.

b)

library(dummies)
dummy(x)

dummy(x)n$x_A$$x_B$$x_C$$x_D$$x_E$

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

(c)

Les probabilités de réussite découlent du modèle logistique:

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

(ré)

$i$$\textrm{Bin}(n, p)$$n = 1$$p =$ pi[i]

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

Quelques simulations rapides pour vérifier que c'est 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 
> #---------------------