如何在 Julia 中创建连续变量的范围?
How to create a range out of a continuous variable in Julia?
我正在 运行使用 Effects 包计算 Julia 中的平均边际效应。我的目标是了解不同年龄段男性和女性的体重变化情况。正如您在下面的输出中看到的,它 运行 是男性和女性每个年龄段的平均边际效应。但是,我想采用一系列年龄变量,而不是单独采用每年。例如,我想要 0:5、5:10、10:15 等年龄范围。这必须在回归模型 运行 之后而不是事先完成。我尝试自己完成它,但我对 Julia 不够流利。
因此唯一需要纠正的行如下:
d1 = Dict(:sex => ["male","female"],:age => [0:5; 6:20])
代码如下:
using DataFrames, Effects, GLM, StatsModels, StableRNGs
rng = StableRNG(42)
growthdata = DataFrame(; age=[13:20; 13:20],
sex=repeat(["male", "female"], inner=8),
weight=[range(100, 155; length=8); range(100, 125; length=8)] .+ randn(rng, 16))
mod_uncentered = lm(@formula(weight ~ 1 + sex * age), growthdata)
d1 = Dict(:sex => ["male","female"],:age => [0:5; 6:20])
ave = effects(d1, mod_uncentered)
输出
sex age weight err lower upper
String Int64 Float64 Float64 Float64 Float64
1 male 0 0.287822 2.88762 -2.5998 3.17545
2 female 0 56.4387 2.88762 53.5511 59.3263
3 male 1 8.00869 2.71603 5.29266 10.7247
4 female 1 59.8481 2.71603 57.1321 62.5641
5 male 2 15.7296 2.54468 13.1849 18.2742
6 female 2 63.2575 2.54468 60.7128 65.8022
7 male 3 23.4504 2.37361 21.0768 25.824
8 female 3 66.6669 2.37361 64.2933 69.0405
9 male 4 31.1713 2.2029 28.9684 33.3742
10 female 4 70.0763 2.2029 67.8734 72.2792
11 male 5 38.8922 2.03264 36.8595 40.9248
12 female 5 73.4857 2.03264 71.4531 75.5184
13 male 6 46.613 1.86295 44.7501 48.476
14 female 6 76.8951 1.86295 75.0322 78.7581
15 male 7 54.3339 1.69399 52.6399 56.0279
16 female 7 80.3046 1.69399 78.6106 81.9985
17 male 8 62.0548 1.52602 60.5288 63.5808
18 female 8 83.714 1.52602 82.1879 85.24
19 male 9 69.7756 1.3594 68.4162 71.135
20 female 9 87.1234 1.3594 85.764 88.4828
21 male 10 77.4965 1.19469 76.3018 78.6912
22 female 10 90.5328 1.19469 89.3381 91.7275
23 male 11 85.2174 1.03282 84.1846 86.2502
24 female 11 93.9422 1.03282 92.9094 94.975
25 male 12 92.9383 0.875345 92.0629 93.8136
26 female 12 97.3516 0.875345 96.4762 98.2269
27 male 13 100.659 0.72515 99.934 101.384
28 female 13 100.761 0.72515 100.036 101.486
29 male 14 108.38 0.587838 107.792 108.968
30 female 14 104.17 0.587838 103.583 1
我不知道 Effects 包,但在 [0:5; 6:20]
中,范围由 julia 自动扩展。你也试过 [0:5, 6:20]
?
对于那些熟悉 R 的人来说,Effects.jl 相当于 effects
包,而不是 emmeans
包。虽然 effects
和 emmeans
存在一定程度的重叠,但 effects
“仅”对预测变量的特定值进行预测,而 emmeans
能够计算边际平均值预测变量的多个值(例如,范围)。
Effects.jl 本质上是做一些事情的包装器:
- 计算一小组预测变量的完全交叉的“参考网格”
- 找到该模型中所有其他预测变量的典型值。 (通常是均值,但您可以使用不同的汇总函数,请注意,您需要考虑您的汇总函数对与分类预测变量相关的对比意味着什么有解释)
- 将这些典型值添加到参考网格中以获得完全指定的数据集,以对
进行预测
- 根据模型参数估计值 (
vcov
) 的 variance-covariance 矩阵计算预测和相关的标准误差。请注意,对于混合模型,这意味着只有固定效应起作用。 (这同样适用于在适用于 lme4 的模型中使用 R 中的效果包。)
换句话说,Effects.jl 不理解范围,它只理解一组值。它不知道如何对 0:5
进行预测,但它确实知道如何对 0
、1
等进行预测
由于您对某个范围内的平均预测感兴趣,因此您可以只计算预测的平均值:
julia> using Statistics
julia> transform!(ave, :age => ByRow(x -> x <= 5 ? "0:5" : "6:20") => :age_bin)
42×7 DataFrame
Row │ sex age weight err lower upper age_bin
│ String Int64 Float64 Float64 Float64 Float64 String
─────┼────────────────────────────────────────────────────────────────────
1 │ male 0 0.287822 2.88762 -2.5998 3.17545 0:5
2 │ female 0 56.4387 2.88762 53.5511 59.3263 0:5
3 │ male 1 8.00869 2.71603 5.29266 10.7247 0:5
4 │ female 1 59.8481 2.71603 57.1321 62.5641 0:5
5 │ male 2 15.7296 2.54468 13.1849 18.2742 0:5
6 │ female 2 63.2575 2.54468 60.7128 65.8022 0:5
7 │ male 3 23.4504 2.37361 21.0768 25.824 0:5
8 │ female 3 66.6669 2.37361 64.2933 69.0405 0:5
9 │ male 4 31.1713 2.2029 28.9684 33.3742 0:5
10 │ female 4 70.0763 2.2029 67.8734 72.2792 0:5
11 │ male 5 38.8922 2.03264 36.8595 40.9248 0:5
12 │ female 5 73.4857 2.03264 71.4531 75.5184 0:5
13 │ male 6 46.613 1.86295 44.7501 48.476 6:20
14 │ female 6 76.8951 1.86295 75.0322 78.7581 6:20
....
julia> rms(x) = sqrt(mean(abs2, x))
rms (generic function with 1 method)
julia> combine(groupby(ave, [:sex, :age_bin]), :weight => mean, :err => rms; renamecols=false)
4×4 DataFrame
Row │ sex age_bin weight err
│ String String Float64 Float64
─────┼────────────────────────────────────
1 │ male 0:5 19.59 2.47686
2 │ female 0:5 64.9622 2.47686
3 │ male 6:20 100.659 1.04247
4 │ female 6:20 100.761 1.04247
对于错误,我使用了 root-mean-square(RMS):换句话说,取相关方差的平均值,然后转换回标准偏差尺度。 (标准误差是检验统计量的抽样分布的标准偏差。)
对于这个特定的模型(well-balanced 数据,没有讨厌的协变量,没有响应的非线性变换),这与您通过取预测变量的平均值然后计算得到的预测相同单个预测:
julia> d2 = Dict(:sex => ["male","female"],:age => [ mean(0:5); mean(6:20)])
Dict{Symbol, Vector} with 2 entries:
:sex => ["male", "female"]
:age => [2.5, 13.0]
julia> effects(d2, mod_uncentered)
4×6 DataFrame
Row │ sex age weight err lower upper
│ String Float64 Float64 Float64 Float64 Float64
─────┼────────────────────────────────────────────────────────
1 │ male 2.5 19.59 2.4591 17.1309 22.0491
2 │ female 2.5 64.9622 2.4591 62.5031 67.4213
3 │ male 13.0 100.659 0.72515 99.934 101.384
4 │ female 13.0 100.761 0.72515 100.036 101.486
误差小一些,因为这里的误差反映的是单个预测的不确定性,而上面的误差反映的是多个预测的不确定性。
我正在 运行使用 Effects 包计算 Julia 中的平均边际效应。我的目标是了解不同年龄段男性和女性的体重变化情况。正如您在下面的输出中看到的,它 运行 是男性和女性每个年龄段的平均边际效应。但是,我想采用一系列年龄变量,而不是单独采用每年。例如,我想要 0:5、5:10、10:15 等年龄范围。这必须在回归模型 运行 之后而不是事先完成。我尝试自己完成它,但我对 Julia 不够流利。
因此唯一需要纠正的行如下:
d1 = Dict(:sex => ["male","female"],:age => [0:5; 6:20])
代码如下:
using DataFrames, Effects, GLM, StatsModels, StableRNGs
rng = StableRNG(42)
growthdata = DataFrame(; age=[13:20; 13:20],
sex=repeat(["male", "female"], inner=8),
weight=[range(100, 155; length=8); range(100, 125; length=8)] .+ randn(rng, 16))
mod_uncentered = lm(@formula(weight ~ 1 + sex * age), growthdata)
d1 = Dict(:sex => ["male","female"],:age => [0:5; 6:20])
ave = effects(d1, mod_uncentered)
输出
sex age weight err lower upper
String Int64 Float64 Float64 Float64 Float64
1 male 0 0.287822 2.88762 -2.5998 3.17545
2 female 0 56.4387 2.88762 53.5511 59.3263
3 male 1 8.00869 2.71603 5.29266 10.7247
4 female 1 59.8481 2.71603 57.1321 62.5641
5 male 2 15.7296 2.54468 13.1849 18.2742
6 female 2 63.2575 2.54468 60.7128 65.8022
7 male 3 23.4504 2.37361 21.0768 25.824
8 female 3 66.6669 2.37361 64.2933 69.0405
9 male 4 31.1713 2.2029 28.9684 33.3742
10 female 4 70.0763 2.2029 67.8734 72.2792
11 male 5 38.8922 2.03264 36.8595 40.9248
12 female 5 73.4857 2.03264 71.4531 75.5184
13 male 6 46.613 1.86295 44.7501 48.476
14 female 6 76.8951 1.86295 75.0322 78.7581
15 male 7 54.3339 1.69399 52.6399 56.0279
16 female 7 80.3046 1.69399 78.6106 81.9985
17 male 8 62.0548 1.52602 60.5288 63.5808
18 female 8 83.714 1.52602 82.1879 85.24
19 male 9 69.7756 1.3594 68.4162 71.135
20 female 9 87.1234 1.3594 85.764 88.4828
21 male 10 77.4965 1.19469 76.3018 78.6912
22 female 10 90.5328 1.19469 89.3381 91.7275
23 male 11 85.2174 1.03282 84.1846 86.2502
24 female 11 93.9422 1.03282 92.9094 94.975
25 male 12 92.9383 0.875345 92.0629 93.8136
26 female 12 97.3516 0.875345 96.4762 98.2269
27 male 13 100.659 0.72515 99.934 101.384
28 female 13 100.761 0.72515 100.036 101.486
29 male 14 108.38 0.587838 107.792 108.968
30 female 14 104.17 0.587838 103.583 1
我不知道 Effects 包,但在 [0:5; 6:20]
中,范围由 julia 自动扩展。你也试过 [0:5, 6:20]
?
对于那些熟悉 R 的人来说,Effects.jl 相当于 effects
包,而不是 emmeans
包。虽然 effects
和 emmeans
存在一定程度的重叠,但 effects
“仅”对预测变量的特定值进行预测,而 emmeans
能够计算边际平均值预测变量的多个值(例如,范围)。
Effects.jl 本质上是做一些事情的包装器:
- 计算一小组预测变量的完全交叉的“参考网格”
- 找到该模型中所有其他预测变量的典型值。 (通常是均值,但您可以使用不同的汇总函数,请注意,您需要考虑您的汇总函数对与分类预测变量相关的对比意味着什么有解释)
- 将这些典型值添加到参考网格中以获得完全指定的数据集,以对 进行预测
- 根据模型参数估计值 (
vcov
) 的 variance-covariance 矩阵计算预测和相关的标准误差。请注意,对于混合模型,这意味着只有固定效应起作用。 (这同样适用于在适用于 lme4 的模型中使用 R 中的效果包。)
换句话说,Effects.jl 不理解范围,它只理解一组值。它不知道如何对 0:5
进行预测,但它确实知道如何对 0
、1
等进行预测
由于您对某个范围内的平均预测感兴趣,因此您可以只计算预测的平均值:
julia> using Statistics
julia> transform!(ave, :age => ByRow(x -> x <= 5 ? "0:5" : "6:20") => :age_bin)
42×7 DataFrame
Row │ sex age weight err lower upper age_bin
│ String Int64 Float64 Float64 Float64 Float64 String
─────┼────────────────────────────────────────────────────────────────────
1 │ male 0 0.287822 2.88762 -2.5998 3.17545 0:5
2 │ female 0 56.4387 2.88762 53.5511 59.3263 0:5
3 │ male 1 8.00869 2.71603 5.29266 10.7247 0:5
4 │ female 1 59.8481 2.71603 57.1321 62.5641 0:5
5 │ male 2 15.7296 2.54468 13.1849 18.2742 0:5
6 │ female 2 63.2575 2.54468 60.7128 65.8022 0:5
7 │ male 3 23.4504 2.37361 21.0768 25.824 0:5
8 │ female 3 66.6669 2.37361 64.2933 69.0405 0:5
9 │ male 4 31.1713 2.2029 28.9684 33.3742 0:5
10 │ female 4 70.0763 2.2029 67.8734 72.2792 0:5
11 │ male 5 38.8922 2.03264 36.8595 40.9248 0:5
12 │ female 5 73.4857 2.03264 71.4531 75.5184 0:5
13 │ male 6 46.613 1.86295 44.7501 48.476 6:20
14 │ female 6 76.8951 1.86295 75.0322 78.7581 6:20
....
julia> rms(x) = sqrt(mean(abs2, x))
rms (generic function with 1 method)
julia> combine(groupby(ave, [:sex, :age_bin]), :weight => mean, :err => rms; renamecols=false)
4×4 DataFrame
Row │ sex age_bin weight err
│ String String Float64 Float64
─────┼────────────────────────────────────
1 │ male 0:5 19.59 2.47686
2 │ female 0:5 64.9622 2.47686
3 │ male 6:20 100.659 1.04247
4 │ female 6:20 100.761 1.04247
对于错误,我使用了 root-mean-square(RMS):换句话说,取相关方差的平均值,然后转换回标准偏差尺度。 (标准误差是检验统计量的抽样分布的标准偏差。)
对于这个特定的模型(well-balanced 数据,没有讨厌的协变量,没有响应的非线性变换),这与您通过取预测变量的平均值然后计算得到的预测相同单个预测:
julia> d2 = Dict(:sex => ["male","female"],:age => [ mean(0:5); mean(6:20)])
Dict{Symbol, Vector} with 2 entries:
:sex => ["male", "female"]
:age => [2.5, 13.0]
julia> effects(d2, mod_uncentered)
4×6 DataFrame
Row │ sex age weight err lower upper
│ String Float64 Float64 Float64 Float64 Float64
─────┼────────────────────────────────────────────────────────
1 │ male 2.5 19.59 2.4591 17.1309 22.0491
2 │ female 2.5 64.9622 2.4591 62.5031 67.4213
3 │ male 13.0 100.659 0.72515 99.934 101.384
4 │ female 13.0 100.761 0.72515 100.036 101.486
误差小一些,因为这里的误差反映的是单个预测的不确定性,而上面的误差反映的是多个预测的不确定性。