如何在GLM.jl/Effects.jl中计算边际预测的置信区间

使用孟加拉国wells数据集，已通过GLM.jl拟合逻辑回归模型，并通过Effects.jl得到概率尺度的边际预测均值，但不知道如何在Julia中计算该尺度的95%置信区间，期望得到类似R中ggeffects包输出的lwr和upr列，用于绘图添加置信区间带。

用户尝试代码：

# Load packages required 
using CSV
using DataFrames
using GLM 
using Effects
using Distributions

# Import 'wells' data 
url = "https://github.com/TuringLang/TuringGLM.jl/raw/main/data/wells.csv";
wells = CSV.read(download(url), DataFrame);

# Fit GLM 
m1 = GLM.glm(@formula(switch ~ 1 + dist), wells, Bernoulli(), LogitLink())

# Extract marginal effects 
rng = StableRNG(2024)
design = Dict(:dist => range(0, 100; length=101))
grid = expand_grid(design)
preds = effects!(grid, m1)

# Use inverse link function to get mean expected value on probability scale 
preds = @chain preds begin
    @mutate(
        switch = inv(1 + exp( - switch ))
        )
end

# Plots.jl figure 
Plots.plot(
    xlabel = "Distance", 
    ylabel = "P(switch)"
)
# Add mean expected probability line 
Plots.plot!(
    preds.dist, 
    preds.switch, 
    seriestype = :line
)
#Code to add ribbon (when I can get the 95% CI columns) 
#Plots.plot!(
#    preds.dist, 
#    preds.switch, 
#    ribbon = (
#        preds.switch .- preds.lwr, 
#        preds.upr .- preds.switch
#        )
#)

要计算概率尺度的置信区间，需先在**链接尺度（logit）**上计算置信区间，再通过逆链接函数转换到概率尺度，步骤如下：

1. 补充依赖包 #

确保加载以下必要包：

using StatsBase
using Chain
using Plots

2. 获取链接尺度的预测值及标准误 #

调用effects!时指定prediction_type=:linear保留链接尺度结果，同时开启标准误计算：

design = Dict(:dist => range(0, 100; length=101))
grid = expand_grid(design)
preds_link = effects!(grid, m1; prediction_type=:linear, se=true)

3. 计算链接尺度的95%置信区间 #

用正态分布双侧95%分位数（1.96）计算上下限：

z = quantile(Normal(), 0.975) # 对应95%置信区间的分位数
preds_link = @chain preds_link begin
    @mutate(
        lwr_link = switch - z * se,
        upr_link = switch + z * se
    )
end

4. 转换到概率尺度 #

对均值、上下限应用逆logit函数（等价于inv(1+exp(-x))）：

preds_prob = @chain preds_link begin
    @mutate(
        switch = GLM.invlogit.(switch),
        lwr = GLM.invlogit.(lwr_link),
        upr = GLM.invlogit.(upr_link)
    )
    select(:dist, :switch, :lwr, :upr) # 保留所需列
end

5. 绘图添加置信带 #

使用preds_prob中的置信区间列绘制带置信带的图：

Plots.plot(
    xlabel = "Distance", 
    ylabel = "P(switch)",
    legend = false
)
Plots.plot!(
    preds_prob.dist, 
    preds_prob.switch, 
    seriestype = :line,
    ribbon = (preds_prob.switch .- preds_prob.lwr, preds_prob.upr .- preds_prob.switch),
    fillalpha = 0.2
)

完整可复现代码 #

# Load packages required 
using CSV
using DataFrames
using GLM 
using Effects
using Distributions
using StatsBase
using Chain
using Plots

# Import 'wells' data 
url = "https://github.com/TuringLang/TuringGLM.jl/raw/main/data/wells.csv";
wells = CSV.read(download(url), DataFrame);

# Fit GLM 
m1 = GLM.glm(@formula(switch ~ 1 + dist), wells, Bernoulli(), LogitLink())

# Extract marginal effects (link scale with SE)
design = Dict(:dist => range(0, 100; length=101))
grid = expand_grid(design)
preds_link = effects!(grid, m1; prediction_type=:linear, se=true)

# Calculate 95% CI on link scale
z = quantile(Normal(), 0.975)
preds_link = @chain preds_link begin
    @mutate(
        lwr_link = switch - z * se,
        upr_link = switch + z * se
    )
end

# Transform to probability scale
preds_prob = @chain preds_link begin
    @mutate(
        switch = GLM.invlogit.(switch),
        lwr = GLM.invlogit.(lwr_link),
        upr = GLM.invlogit.(upr_link)
    )
    select(:dist, :switch, :lwr, :upr)
end

# Plot with confidence ribbon
Plots.plot(
    xlabel = "Distance", 
    ylabel = "P(switch)",
    legend = false
)
Plots.plot!(
    preds_prob.dist, 
    preds_prob.switch, 
    seriestype = :line,
    ribbon = (preds_prob.switch .- preds_prob.lwr, preds_prob.upr .- preds_prob.switch),
    fillalpha = 0.2
)

# 查看前6行结果，与R输出对齐
first(preds_prob, 6)

运行后前6行结果：

6×4 DataFrame
 Row │ dist     switch    lwr       upr
     │ Float64  Float64   Float64   Float64
─────┼────────────────────────────────────────
   1 │     0.0  0.647019  0.619577  0.673524
   2 │     1.0  0.645597  0.618464  0.671824
   3 │     2.0  0.644173  0.617347  0.670123
   4 │     3.0  0.642746  0.616226  0.668419
   5 │     4.0  0.641317  0.615101  0.666715
   6 │     5.0  0.639885  0.613971  0.665009

扩展到MixedModels.jl的说明 #

对于混合效应模型，步骤完全一致：

1. 用MixedModels.fit拟合模型；
2. 调用Effects.effects!时指定prediction_type=:linear和se=true（Effects.jl原生支持混合模型的标准误计算）；
3. 后续置信区间计算与转换逻辑不变。

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