0

我无法使用后验生存预测（）生成后验预测。我正在尝试使用新的数据框，但它没有使用新的数据框，而是使用我用来拟合模型的数据集中的值。模型中的拟合变量是 New.Treatment（6 个处理 = 分类）、Openness（连续光指数 min=2.22、mean=6.903221 和 max=10.54）、subplot_by_site（categorical-720 个站点）、New.Species.name（分类 - 165 种）。我的新数据框有 94 行，而后验生存 () 给了我 3017800 行。请帮忙！

head(nd)
      New.Treatment Openness
1          BE                  5
2          BE                  6
3          BE                  7
4          BE                  8
5          BE                  9
6          BE                 10


fit= stan_surv(formula = Surv(days, Status_surv) ~ New.Treatment*Openness + (1 |subplot_by_site)+(1|New.Species.name),
  data = dataset,
  basehaz = "weibull",
  chains=4,
  iter = 2000,
  cores =4 )

Post=posterior_survfit(fit, type="surv",
                        newdata=nd5)

head(Post)
  id cond_time    time median  ci_lb  ci_ub
1  1        NA 62.0000 0.9626 0.9623 1.0000
2  1        NA 69.1313 0.9603 0.9600 0.9997
3  1        NA 76.2626 0.9581 0.9579 0.9696
4  1        NA 83.3939 0.9561 0.9557 0.9665
5  1        NA 90.5253 0.9541 0.9537 0.9545
6  1        NA 97.6566 0.9522 0.9517 0.9526

##Here some reproducible code to explain my problem:

library(rstanarm)

data_NHN<- expand.grid(New.Treatment = c("A","B","C"), Openness = c(seq(2, 11, by=0.15)))
data_NHN$subplot_by_site=c(rep("P1",63),rep("P2",60),rep("P3",60))
data_NHN$Status_surv=sample(0:1,183, replace=TRUE) 
data_NHN$New.Species.name=c(rep("sp1",10),rep("sp2",40),rep("sp1",80),rep("sp2",20),rep("sp1",33))
data_NHN$days=sample(10, size = nrow(data_NHN), replace = TRUE)

nd_t<- expand.grid(New.Treatment = c("A","B","C"), Openness = c(seq(2, 11, by=1)))


mod= stan_surv(formula = Surv(days, Status_surv) ~ New.Treatment+Openness + (1 |subplot_by_site)+(1|New.Species.name),
                  data =data_NHN,
                  basehaz = "weibull",
                  chains=4,
                  iter = 30,
                  cores =4)

summary(mod)
pos=posterior_survfit(mod, type="surv",
                        newdataEvent=nd_t,
                      times = 0)
head(pos)

 #I am interested in predicting values for specific Openess values  
 #(nd_t=20 rows)but I am getting instead values for each point in time 
 #(pos=18300rows)

操作系统：Mac OS Catalina 10.15.6 R 版本：4.0 rstan 版本：2.21.2 rstanarm 版本：rstanarm_2.21.2 关于为什么它不起作用的任何建议。目前尚不清楚如何给出一个变量在交互中的影响的某种图，因为其他变量和相关的不确定性（即边际效应图）。在我的示例中，我有兴趣获取特定“开放性”值的值，而不是在后验结果中出现的每个特定时间。TIA。

4

1 回答 1

0

我只能给你部分答案。您应该知道这是非常前沿的。尽管 ArXiv 论文（日期为 2020 年 2 月）说

希望在您阅读本文时，该功能将在综合 R 存档网络 (CRAN) 上的稳定版本中提供

但到目前为止，这甚至都不是真的；它甚至不在主 GitHub 分支中，所以我曾经remotes::install_github("stan-dev/rstanarm@feature/survival")从源代码安装它。

最接近的问题是您应该将新数据框指定为newdata，而不是newdataEvent。主分支和这个分支之间以及文档和代码之间似乎有很多不匹配......newdataEvent在旧方法中用于stanjm模型，但不适用于stansurv模型。你可以看看这里的代码，或者使用formals(rstanarm:::posterior_survfit.stansurv). 不幸的是，因为这个方法有一个（未使用的，未经检查的）...参数，这意味着任何错误命名的参数都将被忽略。

下一个问题是，如果您在示例中指定新数据，newdata您将得到

错误：数据中缺少以下变量：subplot_by_site、New.Species.name

也就是说，似乎没有一种明显的方法来生成人口水平的后验预测。NA（不允许将随机效应分组变量设置为。）如果您想这样做，您可以：

  • 扩展您newdata以包含数据集中分组变量的所有组合，并自行对各个级别的结果进行平均；
  • 在 GitHub 上发布问题或联系维护人员...
于 2020-08-04T20:04:50.527 回答