如何基于逻辑回归结果在R中模拟多分类变量？

要模拟包含4个水平的Class变量，核心是先确定该变量的概率分布（边际或联合分布），再基于分布抽样生成变量，最后转换为模型所需的哑变量代入逻辑回归预测式。以下是具体实现步骤：

方法1：基于边际概率模拟Class #

先从原数据中提取Class各水平的出现概率，再用sample()函数按概率抽样，之后转换为模型所需的哑变量：

library(dplyr)
library(tidyr)
library(sjPlot)
# Step 1 - Load the dataset
data(Titanic)
# Step 2 - Transform the Titanic dataset
mydata <- reshape2::melt(Titanic) %>%
  uncount(value) %>%
  as_tibble()
# Step 3 - Create the dummy variables
mydata$Crew <- ifelse(mydata$Class == "Crew", 1, 0)
mydata$First <- ifelse(mydata$Class == "1st", 1, 0)
mydata$Second <- ifelse(mydata$Class == "2nd", 1, 0)
mydata$Third <- ifelse(mydata$Class == "3rd", 1, 0)
mydata$Male <- ifelse(mydata$Sex == "Male", 1, 0)
mydata$Adult <- ifelse(mydata$Age == "Adult", 1, 0)
# Step 4 - Fit the logistic regression model
glm.1 <- glm(Survived ~ Second + Third + Crew + Adult + Male, family = binomial("logit"), data = mydata )
# Step 5 – View the model summary（确认系数准确）
summary(glm.1)

# -------------------------- 模拟多分类变量Class --------------------------
set.seed(016752277)
# 1. 计算原数据中Class各水平的边际概率
class_probs <- table(mydata$Class) / nrow(mydata)
# 2. 按概率模拟2000个Class样本
sim_Class <- sample(names(class_probs), size = 2000, replace = TRUE, prob = class_probs)
# 3. 转换为模型所需的哑变量（参考水平为1st）
sim_Second <- as.integer(sim_Class == "2nd")
sim_Third <- as.integer(sim_Class == "3rd")
sim_Crew <- as.integer(sim_Class == "Crew")
# 4. 模拟其他二分类变量
sim_Male <- sample(c(0,1), size = 2000, replace = TRUE)
sim_Adult <- sample(c(0,1), size = 2000, replace = TRUE)
# 5. 代入原模型系数计算线性预测值（系数需与summary(glm.1)输出一致）
xb <- 3.1054 + (-1.1902)*sim_Second + (-1.5000)*sim_Third + (-0.8379)*sim_Crew + (-1.0615)*sim_Adult + (-2.4201)*sim_Male
p <- 1/(1 + exp(-xb))
sim_Survived <- rbinom(n = 2000, size = 1, prob = p)
# 6. 组合模拟数据集并拟合模型
sim_data <- tibble(Class = sim_Class, Male = sim_Male, Adult = sim_Adult, Survived = sim_Survived,
                   Second = sim_Second, Third = sim_Third, Crew = sim_Crew)
glm.1simulated <- glm(Survived ~ Second + Third + Crew + Adult + Male, family = binomial(), data = sim_data)
# 7. 对比原模型与模拟模型
tab_model(glm.1, glm.1simulated)

方法2：基于联合概率模拟（更贴合原数据关联） #

如果要保留Class与Sex、Age的关联，可以基于三者的联合概率抽样：

# 替换上述模拟部分的代码
set.seed(016752277)
# 1. 计算Class、Sex、Age的联合概率
joint_probs <- prop.table(table(mydata$Class, mydata$Sex, mydata$Age))
joint_df <- as.data.frame(joint_probs) %>% filter(Freq > 0)
# 2. 按联合概率抽样
sim_joint_idx <- sample(1:nrow(joint_df), size = 2000, replace = TRUE, prob = joint_df$Freq)
# 3. 提取抽样后的变量
sim_Class <- joint_df$Var1[sim_joint_idx]
sim_Sex <- joint_df$Var2[sim_joint_idx]
sim_Age <- joint_df$Var3[sim_joint_idx]
# 4. 转换为模型所需的哑变量
sim_Male <- as.integer(sim_Sex == "Male")
sim_Adult <- as.integer(sim_Age == "Adult")
sim_Second <- as.integer(sim_Class == "2nd")
sim_Third <- as.integer(sim_Class == "3rd")
sim_Crew <- as.integer(sim_Class == "Crew")
# 5. 后续计算xb、Survived及拟合模型的步骤与方法1一致
xb <- 3.1054 + (-1.1902)*sim_Second + (-1.5000)*sim_Third + (-0.8379)*sim_Crew + (-1.0615)*sim_Adult + (-2.4201)*sim_Male
p <- 1/(1 + exp(-xb))
sim_Survived <- rbinom(n = 2000, size = 1, prob = p)
sim_data <- tibble(Class = sim_Class, Sex = sim_Sex, Age = sim_Age, Survived = sim_Survived,
                   Second = sim_Second, Third = sim_Third, Crew = sim_Crew, Male = sim_Male, Adult = sim_Adult)
glm.1simulated <- glm(Survived ~ Second + Third + Crew + Adult + Male, family = binomial(), data = sim_data)
tab_model(glm.1, glm.1simulated)

注意事项 #

• 务必使用summary(glm.1)输出的准确系数值代入线性预测式，确保模拟逻辑与原模型一致。
• 方法2生成的模拟数据会保留原数据中变量间的关联（比如船员多为男性成人），模拟结果更贴近真实数据分布。

内容的提问来源于stack exchange，提问作者Paul