1

我有兴趣从 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 函数的产品......

拟合值的示例输出。

4

1 回答 1

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 和统计信息。我怀疑有一些我不太了解的统计问题或解释问题。

于 2016-02-12T00:33:16.357 回答