Simulieren von Daten für die logistische Regression mit einer kategorialen Variablen

Ich habe versucht, einige Testdaten für die logistische Regression zu erstellen, und diesen Beitrag gefunden. Wie simuliere ich künstliche Daten für die logistische Regression?

Es ist eine schöne Antwort, aber es werden nur kontinuierliche Variablen erstellt. Was ist mit einer kategorialen Variablen x3 mit 5 Ebenen (ABCDE), die y für dasselbe Beispiel wie im Link zugeordnet sind?

Das Model

Sei wenn man die Kategorie "B" hat, und sonst x B = 0 . Definieren Sie x C , x D und x E ähnlich. Wenn x B = x C = x D = x E = 0 ist , haben wir die Kategorie "A" (dh "A" ist der Referenzpegel). Ihr Modell kann dann als geschrieben werden$x_B = 1$$x_B = 0$$x_C$$x_D$$x_E$$x_B = x_C = x_D = x_E = 0$

mit β 0 als Achsenabschnitt.

Datengenerierung in R.

(ein)

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

Der xVektor hat nKomponenten (eine für jedes Individuum). Jede Komponente ist entweder "A", "B", "C", "D" oder "E". Jedes von "A", "B", "C", "D" und "E" ist gleich wahrscheinlich.

(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)

Die Erfolgswahrscheinlichkeiten ergeben sich aus dem Logistikmodell:

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

(d)

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

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

Einige schnelle Simulationen, um dies zu überprüfen, sind in Ordnung

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