使用R survey包计算连续协变量风险比的变量类型异常问题

我想在调整awards、comp.imp协变量后，计算连续协变量avg.ed每增加2个单位的患病率比（PR）。尽管avg.ed本身是数值型，代码里也用as.numeric()指定了类型，但不同区间的2单位增量算出的比值存在明显差异，看起来R似乎把它当成了因子处理。

可复现代码 #

library(survey)
data(api)

# 应用调查设计
dstrat <- svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
# 拟合不含目标连续预测变量的模型
model <- svyglm(I(sch.wide=="Yes") ~ awards+comp.imp , design=dstrat, family=quasibinomial()) 
# 获取预测边际值
predmarg<-svypredmeans(model, ~as.numeric(avg.ed), predictat=c(0, 2, 4,6, 8,10 ))
# 尽管指定为连续预测变量，不同2单位增量的比值仍不相同
# 2 vs. 0的PR: 
svycontrast(predmarg,quote(`2`/`0`))
# 4 vs. 2的PR: 
svycontrast(predmarg,quote(`4`/`2`))
# 6 vs. 4的PR: 
svycontrast(predmarg,quote(`6`/`4`))
# 8 vs. 6的PR: 
svycontrast(predmarg,quote(`8`/`6`))
# 10 vs. 8的PR: 
svycontrast(predmarg,quote(`10`/`8`))

输出结果 #

> svycontrast(predmarg,quote(`2`/`0`))
          nlcon    SE
contrast 1.1232 0.055
> svycontrast(predmarg,quote(`4`/`2`))
          nlcon     SE
contrast 1.1282 0.0811
> svycontrast(predmarg,quote(`6`/`4`))
          nlcon     SE
contrast 1.0689 0.0144
> svycontrast(predmarg,quote(`8`/`6`))
         nlcon     SE
contrast 1.025 0.0207
> svycontrast(predmarg,quote(`10`/`8`))
          nlcon     SE
contrast 1.0078 0.0131

核心原因 #

你当前的模型完全没有将avg.ed纳入预测变量，svypredmeans()只是在已拟合的无avg.ed模型基础上，对avg.ed的不同取值做预测边际值计算——本质是把avg.ed当作分类变量拆分处理，自然不同区间的2单位增量比值不会一致。

正确实现方式 #

要得到连续变量avg.ed每增加2单位的恒定患病率比，必须先将avg.ed作为连续预测变量纳入模型，再通过参数转换计算比值：

library(survey)
data(api)

# 应用调查设计
dstrat <- svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
# 将avg.ed作为连续变量纳入模型
model <- svyglm(I(sch.wide=="Yes") ~ awards + comp.imp + avg.ed, design=dstrat, family=quasibinomial()) 

# 计算avg.ed增加2单位的患病率比：利用指数转换线性预测值
pr_ratio <- exp(2 * coef(model)["avg.ed"])
# 用delta方法计算标准误
pr_ratio_se <- 2 * exp(2 * coef(model)["avg.ed"]) * sqrt(vcov(model)["avg.ed", "avg.ed"])

# 输出结果
cat("avg.ed每增加2单位的患病率比：", round(pr_ratio, 4), "\n")
cat("标准误：", round(pr_ratio_se, 4), "\n")

补充验证 #

如果需要确认avg.ed与结局的线性关系是否成立，可以加入多项式项（如poly(avg.ed, 2)）检验非线性显著性：

# 拟合含二次项的模型
model_nonlinear <- svyglm(I(sch.wide=="Yes") ~ awards + comp.imp + poly(avg.ed, 2), design=dstrat, family=quasibinomial())
# 对比线性与非线性模型
anova(model, model_nonlinear)

若非线性显著，再考虑分段回归或样条模型处理；若线性假设成立，上述连续变量模型得到的就是恒定的2单位增量比值。

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