R模拟场景下分组向量构建与处理效应模型拟合问题咨询

问题1：生成分组标识向量的实现方法 #

你可以使用rep()函数生成长度和vect完全匹配的分组向量，两种常用编码方式可选：

• 哑变量编码（1=处理组，0=对照组）：

# 因为treat和cont长度均为3，each参数指定每组重复次数
group <- rep(c(1, 0), each = length(treat))

• 因子编码（更推荐，回归时自动设置参照组）：

# 设置对照组为参照水平，方便直接解读处理效应系数
group <- factor(rep(c("处理组", "对照组"), each = length(treat)), levels = c("对照组", "处理组"))

生成的group向量顺序和vect完全对应，前3位对应处理组观测，后3位对应对照组观测。

问题2：处理效应检验的代码修正与实现 #

原代码存在的问题 #

现有代码的模型设定完全不符合处理效应检验的逻辑：

1. 用glm(treat ~ cont, family = poisson)的设定完全错误：你生成的是正态分布连续值，不应该用泊松回归；且将处理组向量作为因变量、对照组向量作为自变量，样本量仅为3，完全无法体现两组的差异对比。
2. 没有用到合并后的vect向量和分组变量，自然无法检验处理效应。

修正后的实现方案 #

拟合线性模型时，以合并后的观测值vect为因变量，分组变量group为自变量，group的回归系数就是你要检验的处理效应（即两组均值的差值）。修正后的完整模拟代码如下：

library(car)

nsims = 1000
# 预定义存储向量，指定长度提升运行效率
p.value.saved = vector(length = nsims)
treat_effect.saved = vector(length = nsims)

for (i in 1:nsims) {
  treat=rnorm(3, mean = 460, sd = 110)
  cont=rnorm(3, mean = 415, sd = 110)
  vect=c(treat, cont)
  # 生成分组变量
  group = rep(c(1, 0), each = 3)
  # 拟合线性模型检验处理效应
  model = lm(vect ~ group)
  # 保存group系数的p值和效应值
  p.value.saved[i] = Anova(model)$"Pr(>F)"[1]
  treat_effect.saved[i] =  coef(model)[["group"]]
}

# 可以查看模拟的统计功效：效应显著的模拟次数占比
power = mean(p.value.saved < 0.05)
power

如果你需要用泊松回归的场景（比如你的观测值实际是计数数据，生成逻辑要对应调整为泊松分布生成），可以将lm(vect ~ group)替换为glm(vect ~ group, family = poisson)即可。

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