我有兴趣从 clogit 模型中获取设定位置的拟合值。这包括总体水平响应及其周围的置信区间。例如,我的数据大致如下所示:
set.seed(1)
data <- data.frame(Used = rep(c(1,0,0,0),1250),
Open = round(runif(5000,0,50),0),
Activity = rep(sample(runif(24,.5,1.75),1250, replace=T), each=4),
Strata = rep(1:1250,each=4))
在 Clogit 模型中,活动在层内没有变化,因此没有活动主效应。
mod <- clogit(Used ~ Open + I(Open*Activity) + strata(Strata),data=data)
我想要做的是建立一个新的数据框架,我最终可以在 Open 的指定位置绘制边缘拟合值,类似于传统 glm 模型中的新数据设计:例如,
newdata <- data.frame(Open = seq(0,50,1),
Activity = rep(max(data$Activity),51))
但是,当我尝试在 clogit 上运行预测函数时,出现以下错误:
fit<-predict(mod,newdata=newdata,type = "expected")
Surv (rep(1, 5000L), Used) 中的错误:找不到对象“Used”
我意识到这是因为 r 中的 clogit 正在通过 Cox.ph 运行,因此,预测功能试图预测同一层内的成对受试者之间的相对风险(在这种情况下 = 已使用)。
然而,我的问题是是否有办法解决这个问题。这很容易在 Stata 中完成(使用 Margins 命令),并在 Excel 中手动完成,但是我想在 R 中自动化,因为其他所有东西都是在那里编程的。我也在 R 中手动构建了这个(下面的示例代码),但是我一直在我的真实数据中发现似乎不正确的 CI,因此我希望尽可能依赖预测函数。我的手动预测代码是:
coef<-data.frame(coef = summary(mod)$coefficients[,1],
se= summary(mod)$coefficients[,3])
coef$se <-summary(mod)$coefficients[,4]
coef$UpCI <- coef[,1] + (coef[,2]*2) ### this could be *1.96 but using 2 for simplicity
coef$LowCI <-coef[,1] - (coef[,2]*2) ### this could be *1.96 but using 2 for simplicity
fitted<-data.frame(Open= seq(0,50,2),
Activity=rep(max(data$Activity),26))
fitted$Marginal <- exp(coef[1,1]*fitted$Open +
coef[2,1]*fitted$Open*fitted$Activity)/
(1+exp(coef[1,1]*fitted$Open +
coef[2,1]*fitted$Open*fitted$Activity))
fitted$UpCI <- exp(coef[1,3]*fitted$Open +
coef[2,3]*fitted$Open*fitted$Activity)/
(1+exp(coef[1,3]*fitted$Open +
coef[2,3]*fitted$Open*fitted$Activity))
fitted$LowCI <- exp(coef[1,4]*fitted$Open +
coef[2,4]*fitted$Open*fitted$Activity)/
(1+exp(coef[1,4]*fitted$Open +
coef[2,4]*fitted$Open*fitted$Activity))
理想情况下,我的最终产品看起来像这样,但它是 predict 函数的产品......