如何从 clogit 模型中获取拟合值

Question

我有兴趣从 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 模型中，activity 在一个层内没有变化，因此没有 activity 主效应。

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，因此，预测函数试图预测同一层内成对受试者之间的相对风险（在此case= 使用过）。

我的问题是，是否有办法解决这个问题。这很容易在 Stata 中完成（使用边距命令），并在 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))

我的最终产品理想情况下看起来像这样，但它是预测函数的产物....

Example output of fitted values.

Answer 1

显然 Terry Therneau 在 clogit 模型的预测问题上不那么纯粹：http://markmail.org/search/?q=list%3Aorg.r-project.r-help+predict+clogit#query:list%3Aorg.r-project.r-help%20predict%20clogit%20from%3A%22Therneau%2C%20Terry%20M.%2C%20Ph.D.%22+page:1+mid:tsbl3cbnxywkafv6+state:results

这里是对您的代码的修改，它确实生成了 51 个预测。确实需要放入虚拟 Strata 列。

 newdata <- data.frame(Open = seq(0,50,1),
                   Activity = rep(max(data$Activity),51), Strata=1)

 risk <- predict(mod,newdata=newdata,type = "risk")

> risk/(risk+1)
        1         2         3         4         5         6         7 
0.5194350 0.5190029 0.5185707 0.5181385 0.5177063 0.5172741 0.5168418 
        8         9        10        11        12        13        14 
0.5164096 0.5159773 0.5155449 0.5151126 0.5146802 0.5142478 0.5138154 
       15        16        17        18        19        20        21 
0.5133829 0.5129505 0.5125180 0.5120855 0.5116530 0.5112205 0.5107879 
       22        23        24        25        26        27        28 
0.5103553 0.5099228 0.5094902 0.5090575 0.5086249 0.5081923 0.5077596 
       29        30        31        32        33        34        35 
0.5073270 0.5068943 0.5064616 0.5060289 0.5055962 0.5051635 0.5047308 
       36        37        38        39        40        41        42 
0.5042981 0.5038653 0.5034326 0.5029999 0.5025671 0.5021344 0.5017016 
       43        44        45        46        47        48        49 
0.5012689 0.5008361 0.5004033 0.4999706 0.4995378 0.4991051 0.4986723 
       50        51 
0.4982396 0.4978068

{警告}：凡人实际上很难确定哪个 R 神相信这个。我从这些专家那里学到了很多 R 和统计数据。我怀疑有一些我不太了解的统计问题或解释。

如何从 clogit 模型中获取拟合值

How to get fitted values from clogit model

predict

cox-regression

survival-analysis