Closed fweber144 closed 2 years ago
I'm not 100% sure, but I think line https://github.com/melff/mclogit/blob/dadb418e3354824019ac72f175d82c18d1b92a59/pkg/R/mblogit.R#L872 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.
Thanks for the report and the vigorous testing. I'll see to it ASAP.
This appears to be fixed in release 0.9.4.2
Yes, I can confirm this. Thank you!
I'm not 100% sure, but I think line https://github.com/melff/mclogit/blob/dadb418e3354824019ac72f175d82c18d1b92a59/pkg/R/mblogit.R#L872 needs to be removed. Consider the following reprex:
EDIT: I was using mclogit version 0.9.4.1 (at commit dadb418) for this.