Dirichlet回归优势分析：公式语法引发的报错问题排查

目标 #

对Dirichlet回归执行优势分析，估算一组预测变量（标准化连续预测变量、含样条的连续预测变量及因子）的相对重要性。Dirichlet回归是Beta回归的扩展，用于建模非计数来源且分为两类以上的比例（参考Douma&Weedon,2019）。

建模方法 #

使用DirichletReg包拟合Dirichlet回归，采用"alternative"参数化方式：该方式可同时估算参数与精度，语法为response ~ parameters | precision。参数估算的预测变量可与精度估算的不同，例如response ~ predictor1 + predictor2 | predictor3。若未声明精度部分，模型默认固定精度，即response ~ predictors，也可显式写为response ~ predictors | 1。

报错推测 #

报错可能与公式中分隔参数和精度预测变量的竖线有关。使用performance::r2()计算Nagelkerke伪R²作为模型质量指标，实际分析考虑使用McFadden或Estrella伪R²（参考Luchman,2014），因二者适用于多分类响应的优势分析。

遇到的问题 #

执行优势分析时收到报错："fitstat requires at least two elements"。

可复现示例 #

使用DirichletReg包内置数据，响应变量仅两类，仍会出现相同报错：

library(DirichletReg)
#> Warning: package 'DirichletReg' was built under R version 4.1.3
#> Loading required package: Formula
#> Warning: package 'Formula' was built under R version 4.1.1
library(domir)
library(performance)
#> Warning: package 'performance' was built under R version 4.1.3

# 整理数据
RS <- ReadingSkills
RS$acc <- DR_data(RS$accuracy)
#> only one variable in [0, 1] supplied - beta-distribution assumed.
#> check this assumption.
RS$dyslexia <- C(RS$dyslexia, treatment)

# 拟合Dirichlet回归
rs2 <- DirichReg(acc ~ dyslexia + iq | dyslexia + iq, data = RS, model = "alternative")

summary(rs2)
#> Call:
#> DirichReg(formula = acc ~ dyslexia + iq | dyslexia + iq, data = RS, model =
#> "alternative")
#> 
#> Standardized Residuals:
#>                   Min       1Q  Median      3Q     Max
#> 1 - accuracy  -1.5279  -0.7798  -0.343  0.6992  2.4213
#> accuracy      -2.4213  -0.6992   0.343  0.7798  1.5279
#> 
#> MEAN MODELS:
#> ------------------------------------------------------------------
#> Coefficients for variable no. 1: 1 - accuracy
#> - variable omitted (reference category) -
#> ------------------------------------------------------------------
#> Coefficients for variable no. 2: accuracy
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  2.22386    0.28087   7.918 2.42e-15 ***
#> dyslexiayes -1.81261    0.29696  -6.104 1.04e-09 ***
#> iq          -0.02676    0.06900  -0.388    0.698    
#> ------------------------------------------------------------------
#> 
#> PRECISION MODEL:
#> ------------------------------------------------------------------
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  1.71017    0.32697   5.230 1.69e-07 ***
#> dyslexiayes  2.47521    0.55055   4.496 6.93e-06 ***
#> iq           0.04097    0.27537   0.149    0.882    
#> ------------------------------------------------------------------
#> Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Log-likelihood: 61.26 on 6 df (33 BFGS + 1 NR Iterations)
#> AIC: -110.5, BIC: -99.81
#> Number of Observations: 44
#> Links: Logit (Means) and Log (Precision)
#> Parametrization: alternative
as.numeric(performance::r2(rs2))
#> [1] 0.4590758

# 执行优势分析：报错
# 未声明精度部分，模型默认固定精度：parameters | 1
domir::domin(acc ~ dyslexia + iq,
             reg =  function(y)  DirichletReg::DirichReg(y, data = RS, model = "alternative"),
             fitstat = list(\(x) list(r2.nagelkerke = as.numeric(performance::r2(x)), "r2.nagelkerke"))
)
#> Error in domir::domin(acc ~ dyslexia + iq, reg = function(y) DirichletReg::DirichReg(y, : fitstat requires at least two elements.

domir::domin(acc ~ dyslexia + iq | 1,
             reg =  function(y)  DirichletReg::DirichReg(y, data = RS, model = "alternative"),
             fitstat = list(\(x) list(r2.nagelkerke = as.numeric(performance::r2(x)), "r2.nagelkerke"))
             )
#> Error in domir::domin(acc ~ dyslexia + iq | 1, reg = function(y) DirichletReg::DirichReg(y, : fitstat requires at least two elements.

domir::domin(acc ~ dyslexia + iq | dyslexia + iq,
             reg =  function(y)  DirichletReg::DirichReg(y, data = RS, model = "alternative"),
             fitstat = list(\(x) list(r2.nagelkerke = as.numeric(performance::r2(x)), "r2.nagelkerke"))
             )
#> Error in domir::domin(acc ~ dyslexia + iq | dyslexia + iq, reg = function(y) DirichletReg::DirichReg(y, : fitstat requires at least two elements.

domir::domin(acc ~ dyslexia + iq,
             reg =  function(y)  DirichletReg::DirichReg(y, data = RS, model = "alternative"),
             fitstat = list(\(x) list(r2.nagelkerke = as.numeric(performance::r2(x)), "r2.nagelkerke")),
             consmodel = "| dyslexia + iq"
             )
#> Error in domir::domin(acc ~ dyslexia + iq, reg = function(y) DirichletReg::DirichReg(y, : fitstat requires at least two elements.

sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19045)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252   
#> [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                  
#> [5] LC_TIME=Spanish_Spain.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] performance_0.10.0 domir_1.0.1        DirichletReg_0.7-1 Formula_1.2-4     
#> 
#> loaded via a namespace (and not attached):
#>  [1] rstudioapi_0.13  knitr_1.38       magrittr_2.0.3   insight_0.19.1  
#>  [5] lattice_0.20-44  rlang_1.1.0      fastmap_1.1.0    stringr_1.5.0   
#>  [9] highr_0.9        tools_4.1.0      grid_4.1.0       xfun_0.30       
#> [13] cli_3.6.0        withr_2.5.0      htmltools_0.5.2  maxLik_1.5-2    
#> [17] miscTools_0.6-28 yaml_2.3.5       digest_0.6.29    lifecycle_1.0.3 
#> [21] vctrs_0.6.1      fs_1.5.2         glue_1.6.2       evaluate_0.15   
#> [25] rmarkdown_2.13   sandwich_3.0-1   reprex_2.0.1     stringi_1.7.6   
#> [29] compiler_4.1.0   generics_0.1.2   zoo_1.8-9

参考文献 #

• Luchman. Relative Importance Analysis With Multicategory Dependent Variables: An Extension and Review of Best Practices (2014) Organizational research methods
• Douma & Weedon. Analysing continuous proportions in ecology and evolution: A practical introduction to beta and Dirichlet regression (2019) Methods in Ecology and Evolution

报错原因 #

domir包的fitstat参数要求传入的列表必须包含至少两个独立元素：第一个是计算拟合统计量的函数，第二个是该统计量的名称（字符串）。你的写法错误地将函数返回值和名称嵌套在同一个列表里，不符合参数格式要求。

修正后的代码 #

调整fitstat的结构，将函数和名称作为列表的两个独立元素传入：

# 固定精度模型的优势分析
domir::domin(acc ~ dyslexia + iq,
             reg = function(y) DirichletReg::DirichReg(y, data = RS, model = "alternative"),
             fitstat = list(
               function(x) as.numeric(performance::r2(x)),  # 元素1：计算统计量的函数
               "r2.nagelkerke"  # 元素2：统计量名称
             )
)

带精度项的模型处理 #

如果要使用包含精度预测变量的公式（如acc ~ dyslexia + iq | dyslexia + iq），直接将完整公式传入即可，DirichReg支持这种双部分公式：

domir::domin(acc ~ dyslexia + iq | dyslexia + iq,
             reg = function(y) DirichletReg::DirichReg(y, data = RS, model = "alternative"),
             fitstat = list(
               function(x) as.numeric(performance::r2(x)),
               "r2.nagelkerke"
             )
)

更换伪R²类型 #

若要使用McFadden或Estrella伪R²，只需替换fitstat中的计算函数：

• McFadden伪R²：使用performance::r2_mcfadden()
• Estrella伪R²：使用performance::r2_estrella()

示例：

# 使用McFadden伪R²的优势分析
domir::domin(acc ~ dyslexia + iq,
             reg = function(y) DirichletReg::DirichReg(y, data = RS, model = "alternative"),
             fitstat = list(
               function(x) as.numeric(performance::r2_mcfadden(x)),
               "r2.mcfadden"
             )
)

内容的提问来源于stack exchange，提问作者M. Riera