partykit::mob() + mlogit: `错误没有指定合适的拟合函数`

partykit::mob() + mlogit: `Error no suitable fitting function specified`

我试图在 MOB algorithm partykit::mob(). Apparently, it cannot be made directly using the partykit::mob() function (below my attempts). However, I found the LORET 算法生成的树的末端叶子处使用 mlogit::mlogit() 来拟合条件 logit,但我找不到任何包含示例的文档,所以我尝试了从源码中猜测我需要哪个功能,可惜没能实现

您是否知道 (1) 我在哪里可以找到 LORET 库的示例以及 (2) 是否可以使用 partykit:mob() 函数与 mlogit::mlogit 一起工作?提前致谢。


示例数据

为了便于说明,请仔细考虑以下数据。它表示来自 5 个人 (id_ind ) 的数据,他们在 3 个备选方案 (altern) 中做出选择。五个人每人选择三次;因此我们有 15 种选择情况 (id_choice)。每个备选方案由两个通用属性(x1x2)表示,并且选项在 y 中注册(如果选择 1,否则 0)。最后,z1是候选分区变量。

df <- read.table(header = TRUE, text = "
id_ind id_choice altern           x1          x2 y
1       1         1      1  1.586788801  0.11887832 1
2       1         1      2 -0.937965347  1.15742493 0
3       1         1      3 -0.511504401 -1.90667519 0
4       1         2      1  1.079365680 -0.37267925 0
5       1         2      2 -0.009203032  1.65150370 1
6       1         2      3  0.870474033 -0.82558651 0
7       1         3      1 -0.638604013 -0.09459502 0
8       1         3      2 -0.071679538  1.56879334 0
9       1         3      3  0.398263302  1.45735788 1
10      2         4      1  0.291413453 -0.09107974 0
11      2         4      2  1.632831160  0.92925495 0
12      2         4      3 -1.193272276  0.77092623 1
13      2         5      1  1.967624379 -0.16373709 1
14      2         5      2 -0.479859282 -0.67042130 0
15      2         5      3  1.109780885  0.60348187 0
16      2         6      1 -0.025834772 -0.44004183 0
17      2         6      2 -1.255129594  1.10928280 0
18      2         6      3  1.309493274  1.84247199 1
19      3         7      1  1.593558740 -0.08952151 0
20      3         7      2  1.778701074  1.44483791 1
21      3         7      3  0.643191170 -0.24761157 0
22      3         8      1  1.738820924 -0.96793288 0
23      3         8      2 -1.151429915 -0.08581901 0
24      3         8      3  0.606695064  1.06524268 1
25      3         9      1  0.673866953 -0.26136206 0
26      3         9      2  1.176959443  0.85005871 1
27      3         9      3 -1.568225496 -0.40002252 0
28      4        10      1  0.516456176 -1.02081089 1
29      4        10      2 -1.752854918 -1.71728381 0
30      4        10      3 -1.176101700 -1.60213536 0
31      4        11      1 -1.497779616 -1.66301234 0
32      4        11      2 -0.931117325  1.50128532 1
33      4        11      3 -0.455543630 -0.64370825 0
34      4        12      1  0.894843784 -0.69859139 0
35      4        12      2 -0.354902281  1.02834859 0
36      4        12      3  1.283785176 -1.18923098 1
37      5        13      1 -1.293772990 -0.73491317 0
38      5        13      2  0.748091387  0.07453705 1
39      5        13      3 -0.463585127  0.64802031 0
40      5        14      1 -1.946438667  1.35776140 0
41      5        14      2 -0.470448172 -0.61326604 1
42      5        14      3  1.478763383 -0.66490028 0
43      5        15      1  0.588240775  0.84448489 1
44      5        15      2  1.131731049 -1.51323232 0
45      5        15      3  0.212145247 -1.01804594 0
")
df$z1 <- rnorm(n= nrow(df),mean = 0,sd = 1)

mlogit + partykit::mob()

library(mlogit)
library(partykit)
mo <-  mlogit(formula =  y ~ x1 + x2 , 
              data =  df,
              idx  =  c("id_choice", "altern"))
# Coefficients:
#   (Intercept):2  (Intercept):3             x1             x2  
#        0.036497       0.293254       0.821173       1.062794

mlogit_function <-  function(y, x,
                             offset = NULL,
                             ...){ mlogit(y ~  x ,
                                          data =  df)}
formula <-  y ~ x1 + x2 | z1 
mob(formula = formula,
    data    = df,
    fit     = mlogit_function,
    control = mob_control(minsize = 10, 
                          alpha = 0.01))
# Error in mob(formula = formula, data = df, fit = mlogit_function, control = mob_control(minsize = 10,  
# no suitable fitting function specified

mlogit + loret::multinomtree()

这个函数运行树,但它不是我想要的,因为备选方案 2 缺少常量。

loret::multinomtree(formula = formula,
                    data = df)
# Model-based recursive partitioning (NULL)
# Model formula:
#   y ~ x1 + x2 | z1
# 
# Fitted party:
#   [1] root: n = 45
# 1:(Intercept)          1:x1          1:x2 
# -1.1046317     0.7663315     1.0418296  
# 
# Number of inner nodes:    0
# Number of terminal nodes: 1
# Number of parameters per node: 3
# Objective function: 22.62393

mlogit + loret::clmtree()

loret::clmtree(formula = formula,
               data = df)
# Error in clm.fit.default(y = c(1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L,  : 
#                                  is.list(y) is not TRUE

我没有完整的解决方案,但希望有足够的反馈帮助您入门:

  • 重要的是要区分具有替代特定变量(您感兴趣的)的所谓“条件 logit”模型和仅具有特定主题(或个体)的经典“多项 logit 模型” -特定的)变量。 mlogit::mlogit() 可以同时适用(以及混合版本),而 nnet::multinom() 只支持后者。

  • 为了拟合条件 Logit 模型,您可以使用长格式和宽格式的数据。您以长格式指定数据,并且还有一个特定于备选方案的拆分变量 z1。这意味着来自同一个人的数据可能会出现在树的不同节点中,这会很尴尬。

  • 取而代之的是让数据以宽形式出现会更好,这样每一行都对应一个个体,然后你就可以只考虑个体特定的变量来进行拆分。这也将与提供个体特定梯度贡献的拟合 mlogit 对象的 $gradient 元素的视图相匹配。 (这就是 sandwich::estfun() 提取的内容,而这又是 partykit::mob() 的基本信息。)

  • 也可以根据替代特定变量进行合理的递归划分,但我发现很难看出这会产生什么样的模型以及这些模型意味着什么。在任何情况下,您都必须编写自己的代码以从提供完全分解的替代特定梯度贡献的拟合模型对象中提取 estfun

  • loret包有点未完成,好久没更新了。因此,我目前不建议在“生产环境”中使用它。此外 nnet::multinom()(底层 loret::multinomtree())不适合您需要的模型(如上所述)并且 ordinal::clm()(底层 loret::clmtree())适用于完全不同的模型。

  • 我们想要构建到 loret 但尚未完成的一个具体方面是自动检测逻辑模型中的(准)完全分离并在树中适当地处理它.

  • 您的 mlogit + partykit::mob() 方法不起作用,因为拟合函数没有正确的接口(如您正确告知的那样)。有关两个受支持的接口,请参阅 vignette("mob", package = "partykit")

  • 要编写适当的接口,您需要确保在每个子集中拥有所需的所有变量。请注意,响应变量加上回归矩阵是不够的,但您还需要索引变量!我建议通过为 partykit::mob() 指定的公式的 y 变量或 x 变量来包含这些变量。在 mob_control() 中,您可以设置 ytype = "data.frame"xtype = "data.frame"。然后 yx 都作为 data.frame 对象提供,并且可以在调用 mlogit::mlogit() 之前再次组合。 mlogit()formulaidx 参数必须以某种方式提供。在下面的示例中,我对它们进行了硬编码。

根据您的示例进行说明:您可以设置一个模型拟合函数 my_fit(),期望 yx 成为数据框,然后使用 formula = y ~ x1 + x2idx = c("id_choice", "altern") 以适应 mlogit() 模型。

my_fit <- function(y, x = NULL, start = NULL, weights = NULL, offset = NULL, ...) {
  mlogit::mlogit(formula = y ~ x1 + x2, data = cbind(y, x), idx = c("id_choice", "altern"))
}

要指定树,您需要通过 mob() 公式的 x 变量传递 mlogit() 模型所需的所有变量:

tr <- mob(y ~ x1 + x2 + id_choice + altern | z1, data = df, fit = my_fit,
  control = mob_control(ytype = "data.frame", xtype = "data.frame"))

从表面上看,这是有效的。至少它符合根节点中所需的模型,正如您通过检查 coef(tree) 可以看到的那样。然而,它起作用的唯一原因是它甚至不尝试进行任何分区,因为样本量被判断为与参数数量相比太小。

但是当你在 mob_control() 中另外设置 minsize = 10 时,分区过程将失败。这样做的原因是拆分变量的长度为 45,但 mlogit() 中的 gradient/estfun 的长度仅为 15。这是长的特定于替代的格式与短的特定于个人的格式。

因此,要使事情正常进行,您必须使用宽格式的数据,然后相应地调整 mob() 公式和 my_fit() 中的 mlogit() 调用。

资料准备

这是一个基于 Achim 的 的解决方法,可以使 mlogit::mlogitpartykit::mob() 函数一起工作。首先,当 partykit::mob() 函数使用 estfun 提取数据时,我们需要使用宽格式的数据来正确解析评分函数。此外,我在个人级别生成了一个分区变量(z1)(这是我在原始问题中犯的错误(!))。


# First generate the choice variable in the wide format
df$y_wide<- with(df, 
                 ave(x = altern * y , 
                     by = id_choice, 
                     FUN = max))
# drop the long formate choice variable
df <- subset( df, select = -y )
# reshape the data.
df_wide <- reshape(df, 
                   idvar = "id_choice", 
                   timevar = "altern", 
                   direction = "wide",
                   v.names = c("x1", "x2"))
# Add the partition variable
set.seed(777)
# Here I am generating a partition variable at the individual level.
# I made a mistake in my original question generating the partition variable 
# at the choice situation level.
df_wide$z1 <- rep(rnorm(max(df_wide$id_ind)),each = 3)


head(df_wide)
#    id_ind id_choice y_wide        x1.1        x2.1         x1.2       x2.2       x1.3       x2.3         z1
# 1       1         1      1  1.58678880  0.11887832 -0.937965347  1.1574249 -0.5115044 -1.9066752  0.4897862
# 4       1         2      2  1.07936568 -0.37267925 -0.009203032  1.6515037  0.8704740 -0.8255865  0.4897862
# 7       1         3      3 -0.63860401 -0.09459502 -0.071679538  1.5687933  0.3982633  1.4573579  0.4897862
# 10      2         4      3  0.29141345 -0.09107974  1.632831160  0.9292550 -1.1932723  0.7709262 -0.3985414
# 13      2         5      1  1.96762438 -0.16373709 -0.479859282 -0.6704213  1.1097809  0.6034819 -0.3985414
# 16      2         6      3 -0.02583477 -0.44004183 -1.255129594  1.1092828  1.3094933  1.8424720 -0.3985414

mlogit一起实现MOB算法。

第一阶段:拟合函数my_fit_for_mlogit().

调整提供给 partykit::mob 的拟合函数,我们以宽格式显示数据,但我们通过使用 dfidx() 函数“格式化”以供 mlogit::mlogit() 使用。这种处理方式意味着有一个隐藏的 stats::reshape() 进程在幕后工作,因此如果数据集很大,这可能会减慢整个进程。

library(mlogit)
library(partykit)
#### First define the function sketched by Achim in his comment.
my_fit_for_mlogit <- function(y, 
                              x = NULL, 
                              start = NULL, 
                              weights = NULL, 
                              offset = NULL, ...) {

  #First declare dfidx data 
  xy <- cbind(y,x)
  # Properly declare the data to be used by mlogit.
  d <- dfidx(data    = xy, 
             shape   = "wide", 
             choice  = "y_wide",
             varying = 2:7,
             sep     = ".",
             idx     = list(c("id_choice", "id_ind")), 
             idnames = c(NA, "altern"))
  
  # fit mlogit using data equal "d"
  model <- mlogit::mlogit(formula = y_wide ~ x1 + x2|1 ,
                         data    =  d)

    return(model)
}

第二阶段:使用partykit::mob()函数

下面我们调用 mob() 函数内部的 my_fit_for_mlogit(),使用非常大的 alpha (alpha = 0.99) 只是为了说明。这里需要注意的是,我们必须包括宽格式的所有解释变量(例如 x1.1 + x2.1 + x1.2)以及索引变量、选择情况变量(id_choice)和个人索引变量(id_ind) 和分区变量 (z1) 放在 | 符号之后。

# the trick is to include all the variables in wide format and then
# let dfidx() do the magic INSIDE the my_fit_for_mlogit() function. 


tr <- mob(formula = y_wide ~  
                    x1.1 + x2.1 + x1.2 + 
                    x2.2 + x1.3 + x2.3 + 
                    id_choice + id_ind + 0 | z1, 
          data    = df_wide, 
          fit     = my_fit_for_mlogit,
          cluster = id_ind, 
          control = mob_control(ytype = "data.frame", 
                                xtype = "data.frame",
                                alpha = 0.99, # This is just for illustration
                                minsize = 6))

# Model-based recursive partitioning (my_fit_for_mlogit)
# 
# Model formula:
#   y_wide ~ x1.1 + x2.1 + x1.2 + x2.2 + x1.3 + x2.3 + id_choice + 
#   id_ind + 0 | z1
# 
# Fitted party:
#   [1] root
# |   [2] z1 <= 0.48979: n = 9
# |       (Intercept):2 (Intercept):3            x1            x2 
# |          -3.3989117     0.3696964     0.9647511     2.5129869 
# |   [3] z1 > 0.48979: n = 6
# |       (Intercept):2 (Intercept):3            x1            x2 
# |            45.77874     -11.51773      21.72540      31.29221 
# 
# Number of inner nodes:    1
# Number of terminal nodes: 2
# Number of parameters per node: 4
# Objective function: 4.605233