partykit::mob() + mlogit：`未指定合适的拟合函数时出错` 示例数据第一阶段：拟合函数my_fit_for_mlogit()第二阶段：使用partykit::mob()函数

如何解决partykit::mob() + mlogit：`未指定合适的拟合函数时出错` 示例数据第一阶段：拟合函数my_fit_for_mlogit()第二阶段：使用partykit::mob()函数

我试图在由 MOB algorithm mlogit::mlogit() 生成的树的末端叶子上使用 partykit::mob() 来拟合条件 logit。显然，它不能直接使用 partykit::mob() 函数（在我的尝试之下）。然而，我找到了LORET算法，但是我找不到任何带有示例的文档，所以我尝试从源代码中猜测我需要哪个函数，但不幸的是，我无法使其工作。

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

---

示例数据

为了说明，请仔细考虑以下数据。它表示来自 5 个个体 (id_ind ) 的数据，这些个体在 3 个备选方案 (altern) 中进行选择。五个人各选择三遍；因此我们有 15 种选择情况 (id_choice)。每个备选方案由两个通用属性（x1 和 x2）表示，并且选项注册在 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,# 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,1L,: 
#                                  is.list(y) is not TRUE

解决方法

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

• 区分具有替代特定变量（您感兴趣的）的所谓“条件 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 指定的公式的 x 变量或 partykit::mob() 变量包含这些变量。在 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() 的梯度/estfun 的长度仅为 15。这是长的替代特定格式与短的个人特定格式。

因此，要使事情顺利进行，您必须以宽格式使用数据，然后相应地调整 mob() 公式和 mlogit() 内的 my_fit() 调用。

,

数据准备

以下是基于 Achim 的 answer 使 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 算法。

第一阶段：拟合函数my_fit_for_mlogit()。

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

library(mlogit)
library(partykit)
#### First define the function sketched by Achim in his comment.
my_fit_for_mlogit <- function(y,...) {

  #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）和个体索引变量（{{ 1}}) 和分区变量 (id_ind) 放在 z1 符号之后。

|