Incorrect result by predict.mmblogit() in case of multiple random formulas

[fweber144]

I'm not 100% sure, but I think line

mclogit/pkg/R/mblogit.R

Line 872 in dadb418

ZD <- blockMatrix(ZD)

needs to be removed. Consider the following reprex:

data(VA, package = "MASS")
ngrp <- 8L
VA$grp <- gl(n = ngrp,
             k = floor(nrow(VA) / ngrp),
             length = nrow(VA),
             labels = paste0("gr", seq_len(ngrp)))
set.seed(2643)
VA$grp <- sample(VA$grp)

nsbj <- 13L
VA$sbj <- gl(n = nsbj,
             k = floor(nrow(VA) / nsbj),
             length = nrow(VA),
             labels = paste0("subj", seq_len(nsbj)))
VA$sbj <- sample(VA$sbj)

fit_multi <- mclogit::mblogit(cell ~ treat + age + Karn,
                              random = list(~ 1 | sbj,
                                            ~ 1 | grp),
                              data = VA)

VA_new <- head(VA, 1)
VA_new$sbj <- "subj1"
VA_new$grp <- "gr1"
lpreds_multi <- predict(fit_multi, newdata = VA_new)

# Manual check:
VA_new$treat2 <- as.numeric(VA_new$treat == "2")
fixnms_b <- c("treat2", "age", "Karn")
nlats <- nlevels(VA$cell) - 1L
mfixef <- fit_multi$coefmat
mranef <- setNames(fit_multi$random.effects, names(fit_multi$groups))
# Actual result:
lpreds_multi_man <- matrix(mfixef[, "(Intercept)"],
                           nrow = nrow(VA_new),
                           ncol = nlats, byrow = TRUE) +
  as.matrix(VA_new[, fixnms_b, drop = FALSE]) %*% t(mfixef[, fixnms_b]) +
  t(mranef$sbj[1:3, ])
all.equal(unname(lpreds_multi), unname(lpreds_multi_man), tolerance = 1e-15)
## --> Gives TRUE.
# Expected result:
lpreds_multi_man_expect <- lpreds_multi_man +
  t(mranef$grp[1:3, ])
all.equal(unname(lpreds_multi), unname(lpreds_multi_man_expect), tolerance = 1e-15)
## --> Gives "Mean relative difference: 0.2396385".

# After commenting line `ZD <- blockMatrix(ZD)` in mclogit:::predict.mmblogit() and
# running
lpreds_multi <- predict(fit_multi, newdata = VA_new)
# I get
all.equal(unname(lpreds_multi), unname(lpreds_multi_man), tolerance = 1e-15)
## --> Gives "Mean relative difference: 0.2759492".
# and
all.equal(unname(lpreds_multi), unname(lpreds_multi_man_expect), tolerance = 1e-15)
## --> Gives TRUE.
# as expected.

EDIT: I was using mclogit version 0.9.4.1 (at commit dadb418) for this.

[melff]

Thanks for the report and the vigorous testing. I'll see to it ASAP.

[melff]

This appears to be fixed in release 0.9.4.2

[fweber144]

Yes, I can confirm this. Thank you!