I am trying to fit a conditional logit using mlogit::mlogit() at the end leaves of the tree generated by the MOB algorithm partykit::mob(). Apparently, it cannot be made directly using the partykit::mob() function (below my attempts). However, I found the LORET algorithm, but I couldn't find any documentation with examples, so I tried guessing which function I need from the source code, but unfortunately, I couldn't make it work.
Do you know how (1) where I could find examples for the LORET library and (2) if it is possible to use the partykit:mob() function to work together with mlogit::mlogit? Thanks in advance.
For illustration, please gently consider the following data. It represents data from 5 individuals (id_ind ) that choose among 3 alternatives (altern). Each of the five individuals chose three times; hence we have 15 choice situations (id_choice). Each alternative is represented by two generic attributes (x1 and x2), and the choices are registered in y (1 if selected, 0 otherwise). Finally, z1 is a candidate partition variable.
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()This function runs the tree, but it is not what I want because there is a missing constant for alternative 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
I don't have a complete solution but hopefully enough feedback to get you started:
It is important to distinguish so-called "conditional logit" models with alternative-specific variables (which you are interested in) and classic "multinomial logit models" with only subject-specific (or individual-specific) variables. mlogit::mlogit() can fit both (and also mixed versions) while nnet::multinom() only supports the latter.
For fitting conditional logit models you could have the data in long vs. wide form. You specify your data in long form and also have a splitting variable z1 that is alternative-specific. This means that data from the same individual could end up in different nodes of the tree which would be rather awkward.
Instead it would be better to have the data in wide form so that each row corresponds to an individual and then you could only consider individual-specific variables for splitting. This would also match the view of the $gradient element of a fittedmlogit object which provides the individual-specific gradient contributions. (This is what sandwich::estfun() extracts which in turns is the essential information for partykit::mob().)
It might also be possible to do sensible recursive partitioning based on alternative-specific variables but I find it hard to see what kind of models this would yield and what these would mean. In any case, you then would have to write your own code to extract the estfun from the fitted-model object that provides the fully disaggregated alternative-specific gradient contributions.
The loret package is somewhat unfinished and hasn't been updated in quite a while. Hence, I would not recommend using it "in production" at the moment. Also nnet::multinom() (underlying loret::multinomtree()) does not fit the model you need (as mentioned above) and ordinal::clm() (underlying loret::clmtree()) is for a completely different model.
One specific aspect that we wanted to build into loret but did not finish yet, is automatic detection of (quasi-)complete separation in logistic models and handling it appropriately in the tree.
Your mlogit + partykit::mob() approach does not work because the fitting function does not have the right interface (as you are correctly informed). See vignette("mob", package = "partykit") for the two supported interfaces.
To write an appropriate interface you need to make sure that you have all the variables you need in each subset. Note that the response variable plus regressor matrix is not enough for this but you need the index variables as well! I would recommend to include these variables either through the y variables or the x variables of the formula specified for partykit::mob(). In mob_control() you can then set ytype = "data.frame" and xtype = "data.frame". Then both y and x are provided as data.frame objects and can be combined again prior to calling mlogit::mlogit(). The formula and idx arguments for mlogit() have to be provided in some way then. In the example below I have hard-coded them.
Illustration based on your example: You can set up a model-fitting function my_fit() that expects y and x to be data frames and then uses the formula = y ~ x1 + x2 and idx = c("id_choice", "altern") to fit the mlogit() model.
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"))
}
To specifify the tree you then need to pass all required variables for the mlogit() model through the x variables of the mob() formula:
tr <- mob(y ~ x1 + x2 + id_choice + altern | z1, data = df, fit = my_fit,
control = mob_control(ytype = "data.frame", xtype = "data.frame"))
Superficially, this works. At least it fits the desired model in the root node as you can see by checking coef(tree). However, the only reason it works is that it does not even try to do any partitioning as the sample size is judged to be too small compared to the number of parameters.
But when you additionally set minsize = 10 in mob_control() then the partitioning process will fail. The reason for this is that the splitting variable is of length 45 but the gradient/estfun from mlogit() is only of length 15. This is the long alternative-specific format vs. the short individual-specific format.
Thus, to make things work you have to use the data in wide form and then adapt the mob() formula and the mlogit() call inside my_fit() accordingly.