如何用SAS PROC NLMIXED生成带置信带的拟合值图（复刻R结果）

问题背景 #

我正在分析一套1960年起的树木不规则重复监测数据集：

• 包含10000个样地，仅对高度≥1.3m的树木做永久标记监测，新达标树木会纳入
• 每个样地平均约100棵树，测量间隔多为5-7年，部分仅测2次、部分超10次，树木死亡后不再测量
• 核心需求：按物种构建**存活概率（胸径dbh的函数）**模型，生成带置信带的拟合值图，且需根据测量间隔$\Delta t$调整预测值，采用文献推荐的广义logistic模型：
  $$P = \left( \frac{e^\eta}{1 + e^\eta} \right)^{\Delta t}$$

现有尝试 #

1. R中初始实现（未考虑测量间隔） #

用glmmTMB+ggeffects可直接生成带置信带的图，但未按测量间隔调整预测值：

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) 

2. SAS中改用广义logistic模型（考虑测量间隔） #

用PROC NLMIXED实现了调整后的模型，但仅能输出单个观测的预测值，无法生成带置信带的拟合曲线：

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; /* 调整后的logistic函数 */
    model survival ~ binary(prob);  
    random u1 ~ normal(0,ss*ss) subject= tree_id;
    predict prob out=predv;
    predict u1 out=ran;
run;

3. R中手动计算置信区间的逻辑 #

已在R中实现手动推导置信区间并绘图的流程，希望在SAS中复刻：

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

# 读取示例数据
d <- read_sas("hsbdemo.sas7bdat")

# 运行逻辑回归
m <- glmmTMB(HONORS ~ READ + (1|SES), data = d, family = binomial(link = "logit"))

# 提取固定效应参数
coef <- fixef(m)$cond

# 生成预测序列
new_preds <- data.frame(INTERCEPT = 1, READ = seq(min(d$READ),max(d$READ), by = 0.1))

# 计算logit尺度预测值
X <- as.matrix(new_preds)
eta <- X %*% coef

# 计算标准误
vcov_model <- vcov(m)$cond
vcov_preds <- diag(X %*% vcov_model %*% t(X))
se_logit <- sqrt(vcov_preds)

# 转换为概率尺度并计算置信区间
prob <- plogis(eta)
lower <- plogis(eta - 1.96 * se_logit)
upper <- plogis(eta + 1.96 * se_logit)

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

SAS复刻方案 #

以下步骤完全对应R中的手动计算逻辑，生成带置信带的拟合曲线：

步骤1：生成预测序列数据集 #

先提取原数据中dbh的范围，生成均匀间隔的预测序列，并指定目标测量间隔：

/* 获取dbh的取值范围 */
proc sql noprint;
    select min(dbh), max(dbh) into :min_dbh, :max_dbh from d;
quit;

/* 生成dbh预测序列，指定测量间隔（示例用5年） */
data pred_grid;
    do dbh = &min_dbh to &max_dbh by 0.1;
        measurement_interval = 5;
        output;
    end;
run;

步骤2：获取模型参数估计与方差-协方差矩阵 #

修改NLMIXED代码，输出参数估计和协方差矩阵：

proc nlmixed data=d covest outest=param_est;
    parms b0=1 b1=0.01 ss=1; 
    eta = (b0+u1)+b1*dbh ; 
    prob = (exp(eta)/(1+exp(eta)))**measurement_interval;
    model survival ~ binary(prob);  
    random u1 ~ normal(0,ss*ss) subject= tree_id;
run;

/* 提取固定效应参数b0、b1 */
proc sql noprint;
    select estimate into :b0 from param_est where parm='b0';
    select estimate into :b1 from param_est where parm='b1';
quit;

/* 将协方差矩阵转换为数据集格式 */
proc transpose data=param_est(where=(type='COV')) out=cov_matrix;
    id parm;
    var b0 b1;
run;

步骤3：计算logit尺度的预测值与标准误 #

在预测序列中计算logit尺度的线性预测值、方差和标准误：

data pred_calc;
    set pred_grid;
    /* 计算logit尺度的边际预测值（不含随机效应） */
    eta = &b0 + &b1 * dbh;
    
    /* 设计矩阵元素 */
    x0 = 1;
    x1 = dbh;
    
    /* 读取协方差矩阵元素 */
    if _n_=1 then set cov_matrix;
    cov_b0b0 = b0;
    cov_b0b1 = b1;
    cov_b1b1 = col2; /* 注意：需根据cov_matrix实际列名调整 */
    
    /* 计算eta的方差和标准误 */
    var_eta = x0*x0*cov_b0b0 + x1*x1*cov_b1b1 + 2*x0*x1*cov_b0b1;
    se_eta = sqrt(var_eta);
    
    /* 计算logit尺度的95%置信区间上下限 */
    eta_lower = eta - 1.96 * se_eta;
    eta_upper = eta + 1.96 * se_eta;
run;

步骤4：转换为概率尺度并绘图 #

将logit尺度的结果转换为存活概率，并用PROC SGPLOT绘制带置信带的曲线：

data pred_prob;
    set pred_calc;
    /* 应用测量间隔调整，转换为概率尺度 */
    prob = (exp(eta)/(1+exp(eta)))**measurement_interval;
    prob_lower = (exp(eta_lower)/(1+exp(eta_lower)))**measurement_interval;
    prob_upper = (exp(eta_upper)/(1+exp(eta_upper)))**measurement_interval;
run;

/* 绘制拟合曲线与置信带 */
proc sgplot data=pred_prob;
    series x=dbh y=prob / lineattrs=(color=black thickness=2);
    series x=dbh y=prob_lower / lineattrs=(color=red pattern=dash);
    series x=dbh y=prob_upper / lineattrs=(color=red pattern=dash);
    yaxis label="存活概率" min=0 max=1;
    xaxis label="胸径(cm)";
    title="胸径与树木存活概率关系（测量间隔=5年）";
run;

关键说明 #

• 上述代码计算的是边际预测值（整合随机效应后的平均预测），与R中ggpredict默认结果一致
• 若需针对不同测量间隔生成曲线，可在pred_grid中添加多个measurement_interval值，绘图时用group参数区分
• 协方差矩阵的列名可能因PROC TRANSPOSE结果略有差异，需根据实际输出调整引用
• 若需包含随机效应的条件预测，需额外处理随机效应分布，通常拟合曲线用边际预测更合适

内容的提问来源于stack exchange，提问作者Stefan