partykit::mob() + mlogit：` `Error没有指定的合适的拟合函数‘

Question

我试图使用mlogit::mlogit()在MOB算法 partykit::mob()生成的树的尾叶上拟合一个有条件的日志。显然，它不能直接使用partykit::mob()函数(低于我的尝试)。然而，我找到了洛雷特算法，但是我找不到任何带有示例的文档，所以我尝试从源代码中猜测我需要哪个函数，但不幸的是，我无法使它工作。

您知道(1)在哪里可以找到LORET库的示例，以及(2)如果可以使用partykit:mob()函数与mlogit::mlogit一起工作的话？提前谢谢。

示例数据

为便于说明，请仔细考虑以下数据。它代表了来自5个个人(id_ind )的数据，他们选择了3种备选方案(altern)。五个人中的每一个选择了三次，因此我们有15个选择情境(id_choice)。每个选项都由两个泛型属性(x1和x2)表示，选择在y中注册(如果被选中，则在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

Answer 1

我没有一个完整的解决方案，但希望有足够的反馈让你开始：

• 区分所谓的“条件logit”模型和经典的“多项logit模型”(只有特定主题(或个人特定)变量)是很重要的。mlogit::mlogit()可以同时适用于这两种版本(也可以是混合版本)，而nnet::multinom()只支持后者。
• 为了拟合有条件的logit模型，您可以有长的和宽的数据。您可以以长形式指定数据，并且有一个拆分变量z1，该变量是可选的特定变量。这意味着来自同一个体的数据最终可能会出现在树的不同节点上，这将是相当尴尬的。
• 相反，最好让数据以宽的形式出现，这样每一行都对应于一个个体，然后您只能考虑用于拆分的特定于个人的变量。这也将匹配一个合适的$gradient对象的mlogit元素的视图，该对象提供了特定于个人的梯度贡献。(这就是sandwich::estfun()提取的内容，它依次是partykit::mob()的基本信息。)
• 基于可选的特定变量进行合理的递归分区也是可能的，但我发现很难看出这会产生什么样的模型，以及这些模型意味着什么。在任何情况下，您都必须编写自己的代码从贴合模型对象中提取estfun，该对象提供完全分类的可选的梯度贡献。
• loret包有些未完成，很长一段时间没有更新。因此，我不建议在“生产中”使用它。此外，nnet::multinom() (基础loret::multinomtree())不适合您需要的模型(如上所述)，ordinal::clm() (底层loret::clmtree())适用于完全不同的模型。
• 我们希望构建到loret中但尚未完成的一个具体方面是在逻辑模型中自动检测(准)完全分离，并在树中适当地处理它。
• 您的mlogit + partykit::mob()方法不起作用，因为拟合函数没有正确的接口(正如您正确地了解到的那样)。有关两个受支持的接口，请参见vignette("mob", package = "partykit")。
• 要编写适当的接口，需要确保每个子集中都有所需的所有变量。请注意，响应变量加上回归矩阵是不够的，但是您也需要索引变量！我建议通过y变量或为partykit::mob()指定的公式的x变量来包含这些变量。在mob_control()中，您可以设置ytype = "data.frame"和xtype = "data.frame"。然后，y和x都作为data.frame对象提供，并且可以在调用mlogit::mlogit()之前再次组合。因此，必须以某种方式提供formula和idx的mlogit()参数。在下面的例子中，我对它们进行了硬编码。

基于您的示例:您可以设置一个模型拟合函数my_fit()，它期望y和x是数据帧，然后使用formula = y ~ x1 + x2和idx = 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"))
}

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

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)所看到的那样。但是，它工作的唯一原因是它甚至不尝试进行任何分区，因为与参数数相比，样本大小被判断为太小。

但是，当您在minsize = 10中另外设置mob_control()时，分区过程将失败。其原因是分裂变量的长度为45，但mlogit()的梯度/自定义长度仅为15，这是一种特殊的长格式，而非独立的短格式。

因此，为了使事情正常工作，您必须使用宽格式的数据，然后相应地调整mob()公式和my_fit()中的mlogit()调用。

Answer 2

数据准备

这里有一个使work mlogit::mlogit与基于Achim的mlogit::mlogit的partykit::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算法。

第一阶段: fit功能my_fit_for_mlogit()。

调整提供给partykit::mob的fit函数，我们以宽格式显示数据，但是我们通过使用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()函数

下面我们在my_fit_for_mlogit()函数中调用mob()函数，使用一个非常大的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