使用 Julia GLM 的分类变量如何 select 参考水平?

With Julia GLM for categorical variable how to select the reference level?

参考水平似乎被 selected 作为分类值的第一个唯一元素。但是,就我而言,引用是 P(而不是 A)。

ols_all= lm( @formula( value ~ Treatment ), s_treat)

给予

value ~ 1 + Treatment

Coefficients:
────────────────────────────────────────────────────────────────────────
                Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
────────────────────────────────────────────────────────────────────────
(Intercept)   19.0845    0.531803  35.89    <1e-99   18.039     20.13
Treatment: P   5.5775    0.752082   7.42    <1e-12    4.09895    7.05605

我真正想要的是治疗:A(P是安慰剂或对照组)。 当然,我可以重命名变量的值。但是在 SAS 和 R 中,可以 select 引用,因此我希望也有一种方法可以用 Julia GLM 来实现。

GLM.jl不取CategoricalVector的第一个唯一元素,而是以该列的第一个level作为reference。因此,如果您重新排序级别,您可以更改参考以及级别在输出中的出现顺序。这是一个例子:

julia> using CategoricalArrays

julia> using DataFrames

julia> using GLM

julia> y = rand(10)
10-element Vector{Float64}:
 0.6680787249599323
 0.4405942175942186
 0.012595806803754828
 0.21742822324104805
 0.4588945549282415
 0.05463125900208077
 0.5889309471773907
 0.014531957298235865
 0.8444514228200215
 0.13148010370633267

julia> x = categorical(rand(["a", "b", "c"], 10))
10-element CategoricalArray{String,1,UInt32}:
 "b"
 "b"
 "a"
 "a"
 "c"
 "a"
 "c"
 "c"
 "a"
 "b"

julia> df = DataFrame(x=x, y=y)
10×2 DataFrame
 Row │ x     y
     │ Cat…  Float64
─────┼─────────────────
   1 │ b     0.668079
   2 │ b     0.440594
   3 │ a     0.0125958
   4 │ a     0.217428
   5 │ c     0.458895
   6 │ a     0.0546313
   7 │ c     0.588931
   8 │ c     0.014532
   9 │ a     0.844451
  10 │ b     0.13148

julia> lm(@formula(y~x), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}}}}, Matrix{Float64}}

y ~ 1 + x

Coefficients:
────────────────────────────────────────────────────────────────────────
                 Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
────────────────────────────────────────────────────────────────────────
(Intercept)  0.282277     0.165972  1.70    0.1328  -0.110185   0.674739
x: b         0.131108     0.253527  0.52    0.6210  -0.468388   0.730603
x: c         0.0718425    0.253527  0.28    0.7851  -0.527653   0.671338
────────────────────────────────────────────────────────────────────────

julia> levels(df.x)
3-element Vector{String}:
 "a"
 "b"
 "c"

julia> levels!(df.x, ["c", "b", "a"])
10-element CategoricalArray{String,1,UInt32}:
 "b"
 "b"
 "a"
 "a"
 "c"
 "a"
 "c"
 "c"
 "a"
 "b"

julia> lm(@formula(y~x), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}}}}, Matrix{Float64}}

y ~ 1 + x

Coefficients:
───────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)   Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept)   0.354119     0.191648   1.85    0.1071  -0.0990568   0.807295
x: b          0.0592652    0.271031   0.22    0.8331  -0.581622    0.700153
x: a         -0.0718425    0.253527  -0.28    0.7851  -0.671338    0.527653
───────────────────────────────────────────────────────────────────────────

这里描述了更高级的策略:https://juliastats.org/StatsModels.jl/stable/contrasts/

使用对比文档https://juliastats.org/StatsModels.jl/stable/contrasts/ .

还有另一种方法可以设置分类变量的引用:

ols_all= lm(@formula(value ~ Treatment), s_treat, contrasts= Dict(:Treatment => DummyCoding(base="P")))

优点是当级别列表相对较长时重新排序可能很乏味。所以我认为在答案中包含这两个选项会有所帮助。