Wrong effects computed for nestedLogit models about ggeffects

[friendly]

More straightforward comparisons:

library(nestedLogit)
library(ggeffects)
library(effects)
ggpredict(m, c("hincome", "children")) |> plot() # OK 
ggeffect(m, "hincome") |> plot() # bugged!

Effect("hincome", m) |> plot() # OK 

Also tried a multinomial model; that seems OK

[strengejacke]

I must confess that the two functions "ggeffect()" and "ggemmeans()" should deserve a little more love. I started working on this package when I was rather new to R and had no idea of method dispatch. Actually, like John pointed out, I should refactor large parts of the package and make the functions generics, so the code will become a bit clearer and easier to maintain.

But that's a long-term goal, I think related to this issue, I'll indeed "special-casing" nestedLogit.

[strengejacke]

btw, is this in the CRAN or GitHub version of nestedLogit?

[strengejacke]

Ok, got it, the CRAN version was just updated. Btw, can you explain the large differences between what predict() and Effect() return?

library(nestedLogit)
data("Womenlf", package = "carData")
m <- nestedLogit(partic ~ hincome + children,
  logits(
    work = dichotomy("not.work", c("parttime", "fulltime")),
    full = dichotomy("parttime", "fulltime")
  ),
  data = Womenlf
)

effects::Effect("hincome", m, xlevels = list(hincome = seq(10, 40, 10)))
#> 
#> hincome effect (probability) for not.work
#> hincome
#>        10        20        30        40 
#> 0.5472608 0.6485550 0.7380345 0.8113586 
#> 
#> hincome effect (probability) for parttime
#> hincome
#>        10        20        30        40 
#> 0.1656323 0.2206211 0.2177867 0.1764002 
#> 
#> hincome effect (probability) for fulltime
#> hincome
#>         10         20         30         40 
#> 0.28710686 0.13082390 0.04417872 0.01224116
predict(m, newdata = ggeffects::new_data(m, terms = c("hincome[10:40, by=10]")))
#> 
#> predicted response-category probabilties
#>    not.work   parttime  fulltime
#> 1 0.2864384 0.05907221 0.6544894
#> 2 0.3799734 0.12943673 0.4905899
#> 3 0.4833619 0.22495878 0.2916793
#> 4 0.5881947 0.28527228 0.1265330
#> 
#> predicted response-category logits
#>      not.work   parttime    fulltime
#> 1 -0.91274548 -2.7681058  0.63883223
#> 2 -0.48966118 -1.9059482 -0.03764493
#> 3 -0.06657687 -1.2369990 -0.88724209
#> 4  0.35650744 -0.9184576 -1.93196717
#> 
#> standard errors of predicted probabilities
#>     not.work   parttime   fulltime
#> 1 0.05355137 0.02685851 0.05580606
#> 2 0.06155651 0.04714845 0.06655909
#> 3 0.09546538 0.09563873 0.10160578
#> 4 0.13317026 0.12816405 0.09792827
#> 
#> standard errors of predicted logits
#>    not.work  parttime  fulltime
#> 1 0.2620040 0.4832173 0.2467842
#> 2 0.2612826 0.4184172 0.2663307
#> 3 0.3822848 0.5485371 0.4917936
#> 4 0.5497867 0.6285879 0.8860491

^{Created on 2023-06-23 with reprex v2.0.2}

[friendly]

I think John will have to look at your example where you use:

predict(m, newdata = ggeffects::new_data(m, terms = c("hincome[10:40, by=10]")))

predict simply calculates the predicted probababilities, SEs etc evaluated at the values given by newdata. (It can also do this for the individual nested dichotomies, but that's not the issue here.)

Effect, when not all predictors are specified, calculates the predicted probs. at the margins/average of the omitted variables.

[john-d-fox]

As I understand it, and without checking in detail, predict() sets children to "absent" (its first level) while Effect() computes effects marginalized over the two levels of children. Labour-force participation is higher when children are absent.

BTW, to see confidence limits computed by Effect() for comparison, use summary().

[john-d-fox]

With respect to making the functions in ggeffects generic: I'm not familiar enough with your package to be confident in offering advice, but from what I've seen, the main functions, like ggeffect(), could be simple S3 generics and the current versions of the functions, or something very close to them, could just become the default methods. That way, what works now should continue to work, and it would be possible to introduce specific methods for cases, like "nestedLogit" objects that currently don't work.

[friendly]

Thanks for following up!