使用PROC NLMIXED时如何获得带有置信带的拟合值

问题描述 投票：0回答：1

我正在分析一个数据集，其中包含树木测量值，例如直径（厘米）以及它是死是活（0/1）。测量收集是不定期的，即从 1960 年开始，至今仍在测量中。迄今为止，已在广阔的地理区域内总共建立了 10,000 个样地，所有高度超过 1.3 m 的树木都被永久标记和监控。这也意味着，如果该地块中长出了新树并且它们超过了 1.3m，它们也会包含在集合中。每个地块平均包含约 100 棵树。大多数地块每 5-7 年左右测量一次，有些地块只测量两次，有些地块超过 10 次，这意味着我对独特的树木进行了重复测量。一旦一棵树死亡，它的最终直径就会被记录下来，然后排除在进一步的测量之外。

这是数据子集的一些结构：

> glimpse(d)
Rows: 15,472
Columns: 12
$ plot_number          <dbl> 272.1, 272.1, 272.1, 272.1, 272.1, 272.1, 272.1, 27…
$ establishment_year   <dbl> 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 196…
$ measurement_year     <dbl> 1977, 1977, 1977, 1977, 1977, 1977, 1977, 1977, 197…
$ measurement_interval <dbl> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,…
$ tree_number          <dbl> 1, 5, 6, 7, 12, 15, 18, 19, 20, 22, 24, 26, 27, 29,…
$ species              <chr> "Pl", "Pl", "Pl", "Pl", "Pl", "Pl", "Pl", "Pl", "Pl…
$ dbh                  <dbl> 27.9, 30.2, 13.0, 14.2, 10.9, 19.8, 9.1, 36.8, 24.4…
$ survival             <dbl> 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ tree_id              <chr> "272.1-1", "272.1-5", "272.1-6", "272.1-7", "272.1-…

> d %>% 
+   group_by(plot_number, measurement_year) %>% 
+   summarise(number_of_trees = n())
`summarise()` has grouped output by 'plot_number'. You can override using
the `.groups` argument.
# A tibble: 210 × 3
# Groups:   plot_number [72]
           plot_number measurement_year number_of_trees
                 <dbl>            <dbl>           <int>
 1                 6.1             1965              33
 2                 6.1             1968              37
 3                 6.1             1973              34
 4                 6.1             1982              23
 5                 6.1             1992              11
 6                 6.2             1965              35
 7                 6.2             1968              34
 8                 6.2             1973              34
 9                 6.2             1982              28
10                 6.2             1992              20
# ℹ 200 more rows
# ℹ Use `print(n = ...)` to see more rows

> d %>% 
+   group_by(plot_number, tree_number) %>% 
+   summarise(tree_measurements = n()) %>% 
+   arrange(plot_number, tree_number) 
`summarise()` has grouped output by 'plot_number'. You can override using
the `.groups` argument.
# A tibble: 6,725 × 3
# Groups:   plot_number [72]
           plot_number tree_number tree_measurements
                 <dbl>       <dbl>             <int>
 1                 6.1         531                 4
 2                 6.1         532                 3
 3                 6.1         533                 5
 4                 6.1         534                 3
 5                 6.1         536                 3
 6                 6.1         537                 4
 7                 6.1         540                 2
 8                 6.1         542                 4
 9                 6.1         543                 2
10                 6.1         546                 3
# ℹ 6,715 more rows
# ℹ Use `print(n = ...)` to see more rows

我想使用这些数据来生成按物种的生存模型，即生存作为直径的函数（带有置信带的拟合值）。

在R中，我是这样做的：

require(glmmTMB)
require(ggeffects)
m <- glmmTMB(survival ~ dbh + (1|tree_id) , 
             family = binomial(link = "logit"), 
             data = d)

m_pred <- ggpredict(m, terms = "dbh [all]")

plot(m_pred)

但是，需要根据同一棵树的后续测量之间的时间来调整预测值。在文献中，对于 5-7 年的间隔长度，建议使用以下广义逻辑模型（A Generalized Mixed Logistic Model for Predicting individual Tree Survival Probability with UnequalMeasurement Intervals）：

$$P = \left( rac{e^ ta}{1 + e^ ta} 右)^{\Delta t}$$

这是我必须转向 SAS proc nlmixed 的地方，因为我无法找到在 R 中编码的方法（尽管有一个 @BenBolker hack here 似乎适用于我的数据）。不管怎样，我也想在 SAS 中运行它，然后与 Ben Bolker 的代码进行比较。这是我使用的 SAS 代码：

proc nlmixed data=d;
    parms b0=1 b1=0.01 ss=1; 
    eta = (b0+u1)+b1*dbh ; 
    prob = (exp(eta)/(1+exp(eta)))**measurement_interval; /*ADJUSTED LOGISTIC FUNCTION*/
    model survival ~ binary(prob);  
     random u1 ~ normal(0,ss*ss) subject= tree_id;
      predict prob out=predv;
      predict u1 out=ran;
      
run;

我现在无法弄清楚的是如何获得与上面相同的图（带有置信带的拟合值）...我确实使用了SAS中的

predict
选项，但预测值似乎是个体的预测值观测值而不是拟合值。虽然，我可以使用估计的模型参数创建拟合线（如上所述），但我不知道如何为这些拟合值添加置信带。

在此网站上（https://stats.oarc.ucla.edu/sas/faq/how-can-i-run-simple-linear-and-nonlinear-models-using-nlmixed/）是一个链接PROC NLMIXED 的 SAS 玩具数据集。使用这个数据集，我运行了以下代码，生成了这个图表（我不想要）：

proc nlmixed data="C:\hsbdemo.sas7bdat";
  parms b0=5 b1=0 ss=1;
  xb=(b0+u1)+b1*read;
  prob = exp(xb)/(1+exp(xb));
  model honors ~ binary(prob);
  random u1 ~ normal(0,ss*ss) subject = ses;
  predict prob out=pred;
run;

proc sgplot data=pred;
    series x=read y=pred / lineattrs=(color=blue);
    band x=read lower=lower upper=upper / transparency=0.5;
    xaxis label="Read";
    yaxis label="Predicted Probability";
run;

所以我的问题是如何使用 SAS 和 PROC NLMIXED 复制 R 中生成的图（上面的第一个图）？

编辑

根据@Tom的建议，这是一种计算所需置信区间的方法。以下示例是在 R 中使用的，但我希望它可以在 SAS 中工作，以在 PROC NLMIXED 中运行模型。我添加了注释，希望能够足够清楚地解释这些步骤，以便那些不了解 R 但了解 SAS 的人可以复制它。然而，理想情况下，我想要一个 SAS 程序，通过编程来执行此操作，而不是手动执行此操作。

require(haven)
require(ggeffects)
require(glmmTMB)

# read data (https://stats.oarc.ucla.edu/sas/faq/how-can-i-run-simple-linear-and-nonlinear-models-using-nlmixed/)
d <- read_sas("hsbdemo.sas7bdat")

# run logistic regression
m <- glmmTMB(HONORS ~ READ + (1|SES), data = d, family = binomial(link = "logit"))

# extract fixed effects parameter estimates
coef <- fixef(m)$cond

# create data frame with sequence of predictions
new_preds <- data.frame(INTERCEPT = 1, READ = seq(min(d$READ),max(d$READ), by = 0.1))

# create design matrix
X <- as.matrix(new_preds)

# matrix multiply design matrix with parameter vector to get predictions
eta <- X %*% coef

# extract variance-covariance matrix
vcov_model <- vcov(m)$cond

# matrix multiply design matrix with variance covariance matrix with the transpose of the design matrix 
vcov_preds <- diag(X %*% vcov_model %*% t(X))

# calculate standard deviations (errors)
se_logit <- sqrt(vcov_preds)

# transform prediction from logit to probability scale
prob <- plogis(eta)

# calculate confidence intervals and transform from logit to probability scale
lower <- plogis(eta - 1.96 * se_logit)
upper <- plogis(eta + 1.96 * se_logit)

# plot results
plot(new_preds$READ, prob)
lines(new_preds$READ, lower, col = "red")
lines(new_preds$READ, upper, col = "red")

这就是它应该的样子（不像上面的那样）。

这个图可以重现为适合 PROC NLMIXED 的模型吗？

r sas logistic-regression
1个回答
0
投票

我在寻找一种使用 nlmixed 在回归线上添加置信区间带的方法时看到了你的帖子。如果通过 read(x 变量) 对预测输出进行排序，由于随机效应，您将得到合理但锯齿状的线条和带。

proc sort data=pred_n;
by read;
run;

[绘图前对 x 变量进行排序][1]

https://agstats.io/tutorials/sas-nonlinear的帖子（8混合模型估计下的示例）展示了一个调整随机效应预测函数的示例，无需预测语句中的随机参数即可得到像这样的：

predict  exp(b0+b1*(read))/(1+exp(b0+b1*(read))) out=pred;

[预测函数中没有随机参数][2]

我是 nlmixed 的新手，也许一些高级用户可以为此目的发布正确的模型设置？
[1]：https://i.stack.imgur.com/URZXO.png [2]：https://i.stack.imgur.com/Lb4uw.png

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.