Mensimulasikan data untuk regresi logistik dengan variabel kategori

Saya mencoba membuat beberapa data uji untuk regresi logistik dan saya menemukan posting ini. Bagaimana mensimulasikan data buatan untuk regresi logistik?

Ini adalah jawaban yang bagus tetapi hanya menciptakan variabel kontinu. Bagaimana dengan variabel kategori x3 dengan 5 level (ABCDE) yang dikaitkan dengan y untuk contoh yang sama seperti pada tautan?

Model

Misalkan jika seseorang memiliki kategori "B", dan x B = 0 sebaliknya. Tentukan x C , x D , dan x E similary. Jika x B = x C = x D = x E = $x_B = 1$$x_B = 0$$x_C$$x_D$$x_E$$x_B = x_C = x_D = x_E = 0$ , maka kita memiliki kategori "A" (yaitu, "A" adalah level referensi). Model Anda kemudian dapat ditulis sebagai

dengan β 0 intersep.

Pembuatan data dalam R

(Sebuah)

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

The xvektor memiliki nkomponen (satu untuk setiap individu). Setiap komponen dapat berupa "A", "B", "C", "D", atau "E". Setiap "A", "B", "C", "D", dan "E" memiliki kemungkinan yang sama.

(b)

library(dummies)
dummy(x)

dummy(x)adalah matriks dengan nbaris (satu untuk setiap individu) dan 5 kolom sesuai dengan , x B , x C , x D , dan x E . Prediktor linier (satu untuk setiap individu) kemudian dapat ditulis sebagai$x_A$$x_B$$x_C$$x_D$$x_E$

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

(c)

Peluang keberhasilan mengikuti dari model logistik:

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)

Beberapa simulasi cepat untuk memeriksa ini 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 
> #---------------------