我正在尝试使用该gmnl包估计一个混合混合多项 logit 模型。它在不包含替代特定常量 (ASC) 时完美运行，但在合并它们时会产生奇怪的错误。下面的代码取自（并改编）自该软件包发布的原始文章。
数据准备
options(digits = 3)
library("gmnl")
library("mlogit")
data("Electricity", package = "mlogit")
Electr <- mlogit.data(Electricity,
id.var = "id",
choice = "choice",
varying = 3:26,
shape = "wide",
sep = "")
####Alternative Specific Constants
Electr$asc2 <- as.numeric(Electr$alt == 2)
Electr$asc3 <- as.numeric(Electr$alt == 3)
Electr$asc4 <- as.numeric(Electr$alt == 4)
潜在类模型（使用 ASC）
下面的代码可以完美运行，即使在公式的第二部分 ( LC_ASC_in_formula) 中包含 ASC 或显式包含回归量 ( LC_ASC_in_variables)。
LC_ASC_in_formula <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 1 | 0 | 0 | 1,
data = Electr,
subset = 1:3000,
model = "lc",
panel = TRUE,
Q = 2)
summary(LC_ASC_in_formula)
LC_ASC_in_variables <- gmnl( choice ~ pf + cl + loc + wk + tod + seas +asc2 +asc3 +asc4 | 0 | 0 | 0 | 1,
data = Electr,
subset = 1:3000,
model = "lc",
panel = TRUE,
Q = 2)
summary(LC_ASC_in_variables)
## Are they the same?
logLik(LC_ASC_in_variables) == logLik(LC_ASC_in_formula)
## [1] TRUE
混合混合MNL模型
该模型基本上是一个潜在类模型，但在每个类中，参数是随机的（遵循先前指定的参数分布）。
没有ASC的混合混合 MNL 模型
当省略 ASC 时，该模型工作得很好。
MM_no_ASC <- gmnl(choice ~ pf + cl + loc + wk + tod + seas | 0 | 0 | 0 | 1,
data = Electr,
subset = 1:3000,
model = "mm",
R = 5,
panel = TRUE,
ranp = c(pf = "n",cl = "n",loc = "n",wk = "n", tod = "n",seas= "n"),
Q = 2,
iterlim = 500)
但是，当包含 ASC 时，它无法估计模型：
- 作为模型中变量的一部分。
MM_ASC_in_variables <- gmnl( choice ~ pf + cl + loc + wk + tod + seas +
asc2 +asc3 +asc4 | 0 | 0 | 0 | 1 ,
data = Electr,
subset = 1:3000,
model = "mm",
R = 5,
panel = TRUE,
ranp = c(pf = "n",cl = "n",loc = "n",wk = "n", tod = "n",seas= "n"),
Q = 2,
iterlim = 500)
> Error in if (distr == "n") { : missing value where TRUE/FALSE needed
- 当将它们包含在
formula.
MM_ASC_in_formula <- gmnl( choice ~ pf + cl + loc + wk + tod + seas | 1 | 0 | 0 | 1 ,
data = Electr,
subset = 1:3000,
model = "mm",
R = 5,
panel = TRUE,
ranp = c(pf = "n",cl = "n",loc = "n",wk = "n", tod = "n",seas= "n"),
Q = 2,
iterlim = 500)
> Error in if (distr == "n") { : missing value where TRUE/FALSE needed
然而，两种包含 ASC 参数的方法都无法初始化模型估计。希望有人可以帮助我解决这个问题。先感谢您。
Bonus1：错误的追溯。我减少了估计 ( subset = 1:20) 中包含的观察数，以便更好地查看traceback()下面显示的误差。但我自己无法发现错误。
MM_ASC_in_formula <- gmnl( choice ~ pf + cl + loc + wk + tod + seas | 1 | 0 | 0 | 1 ,
data = Electr,
subset = 1:20,
model = "mm",
R = 5,
panel = TRUE,
ranp = c(pf = "n",cl = "n",loc = "n",wk = "n", tod = "n",seas= "n"),
Q = 2,
iterlim = 500)
# Error in if (distr == "n") { : missing value where TRUE/FALSE needed
traceback()
# Estimating MM-MNL model
# Error in if (distr == "n") { : missing value where TRUE/FALSE needed
# > traceback()
# 14: Makeh.rcoef(beta[, q], stds[, q], ranp, Omega[, ((i - 1) * R +
# 1):(i * R), drop = FALSE], correlation, Pi = NULL, Slist = NULL,
# mvar = NULL)
# 13: fnOrig(theta, ...)
# 12: logLikFunc(theta, fnOrig = function (theta, y, X, H, Q, id = NULL,
# ranp, R, correlation, weights = NULL, haltons = NULL, seed = 12345,
# gradient = TRUE, get.bi = FALSE)
# {
# K <- ncol(X[[1]])
# J <- length(X)
# N <- nrow(X[[1]])
# panel <- !is.null(id)
# if (panel) {
# n <- length(unique(id))
# if (length(weights) == 1)
# weights <- rep(weights, N)
# }
# beta <- matrix(theta[1L:(K * Q)], nrow = K, ncol = Q)
# nstds <- if (!correlation)
# K * Q
# else (0.5 * K * (K + 1)) * Q
# stds <- matrix(theta[(K * Q + 1):(K * Q + nstds)], ncol = Q)
# rownames(beta) <- colnames(X[[1]])
# colnames(beta) <- colnames(stds) <- paste("class", 1:Q, sep = ":")
# gamma <- theta[-c(1L:(K * Q + nstds))]
# ew <- lapply(H, function(x) exp(crossprod(t(x), gamma)))
# sew <- suml(ew)
# Wnq <- lapply(ew, function(x) {
# v <- x/sew
# v[is.na(v)] <- 0
# as.vector(v)
# })
# Wnq <- Reduce(cbind, Wnq)
# set.seed(seed)
# Omega <- make.draws(R * ifelse(panel, n, N), K, haltons)
# XBr <- vector(mode = "list", length = J)
# for (j in 1:J) XBr[[j]] <- array(NA, dim = c(N, R, Q))
# nind <- ifelse(panel, n, N)
# if (panel)
# theIds <- unique(id)
# if (get.bi)
# bi <- array(NA, dim = c(nind, R, Q, K), dimnames = list(NULL,
# NULL, NULL, colnames(X[[1]])))
# for (i in 1:nind) {
# if (panel) {
# anid <- theIds[i]
# theRows <- which(id == anid)
# }
# else theRows <- i
# for (q in 1:Q) {
# bq <- Makeh.rcoef(beta[, q], stds[, q], ranp, Omega[,
# ((i - 1) * R + 1):(i * R), drop = FALSE], correlation,
# Pi = NULL, Slist = NULL, mvar = NULL)
# for (j in 1:J) {
# XBr[[j]][theRows, , q] <- crossprod(t(X[[j]][theRows,
# , drop = FALSE]), bq$br)
# }
# if (get.bi)
# bi[i, , q, ] <- t(bq$br)
# }
# }
# EXB <- lapply(XBr, function(x) exp(x))
# SEXB <- suml.array(EXB)
# Pntirq <- lapply(EXB, function(x) x/SEXB)
# Pnrq <- suml.array(mapply("*", Pntirq, y, SIMPLIFY = FALSE))
# if (panel)
# Pnrq <- apply(Pnrq, c(2, 3), tapply, id, prod)
# Pnq <- apply(Pnrq, c(1, 3), mean)
# WPnq <- Wnq * Pnq
# Ln <- apply(WPnq, 1, sum)
# if (get.bi)
# Qir <- list(wnq = Wnq, Ln = Ln, Pnrq = Pnrq)
# lnL <- if (panel)
# sum(log(Ln) * weights[!duplicated(id)])
# else sum(log(Ln) * weights)
# if (gradient) {
# lambda <- mapply(function(y, p) y - p, y, Pntirq, SIMPLIFY = FALSE)
# Wnq.mod <- aperm(repmat(Wnq/Ln, dimen = c(1, 1, R)),
# c(1, 3, 2))
# Qnq.mod <- Wnq.mod * Pnrq
# if (panel)
# Qnq.mod <- Qnq.mod[id, , ]
# eta <- lapply(lambda, function(x) x * Qnq.mod)
# dUdb <- dUds <- vector(mode = "list", length = J)
# for (j in 1:J) {
# dUdb[[j]] <- array(NA, dim = c(N, K, Q))
# dUds[[j]] <- array(NA, dim = c(N, nrow(stds), Q))
# }
# for (i in 1:nind) {
# if (panel) {
# anid <- theIds[i]
# theRows <- which(id == anid)
# }
# else theRows <- i
# for (q in 1:Q) {
# bq <- Makeh.rcoef(beta[, q], stds[, q], ranp,
# Omega[, ((i - 1) * R + 1):(i * R), drop = FALSE],
# correlation, Pi = NULL, Slist = NULL, mvar = NULL)
# for (j in 1:J) {
# dUdb[[j]][theRows, , q] <- tcrossprod(eta[[j]][theRows,
# , q, drop = TRUE], bq$d.mu)
# dUds[[j]][theRows, , q] <- tcrossprod(eta[[j]][theRows,
# , q, drop = TRUE], bq$d.sigma)
# }
# }
# }
# if (correlation) {
# vecX <- c()
# for (i in 1:K) {
# vecX <- c(vecX, i:K)
# }
# Xac <- lapply(X, function(x) x[, vecX])
# }
# else {
# Xac <- X
# }
# Xr <- lapply(X, function(x) x[, rep(1:K, Q)])
# Xacr <- lapply(Xac, function(x) x[, rep(1:ncol(Xac[[1]]),
# Q)])
# dUdb <- lapply(dUdb, function(x) matrix(x, nrow = N))
# dUds <- lapply(dUds, function(x) matrix(x, nrow = N))
# grad.beta <- suml(mapply("*", Xr, dUdb, SIMPLIFY = FALSE))/R
# grad.stds <- suml(mapply("*", Xacr, dUds, SIMPLIFY = FALSE))/R
# Qnq <- WPnq/Ln
# if (panel) {
# Wnq <- Wnq[id, ]
# H <- lapply(H, function(x) x[id, ])
# Qnq <- Qnq[id, ]
# }
# Wg <- vector(mode = "list", length = Q)
# IQ <- diag(Q)
# for (q in 1:Q) Wg[[q]] <- rowSums(Qnq * (repRows(IQ[q,
# ], N) - repCols(Wnq[, q], Q)))
# grad.gamma <- suml(mapply("*", H, Wg, SIMPLIFY = FALSE))
# gari <- cbind(grad.beta, grad.stds, grad.gamma)
# colnames(gari) <- names(theta)
# attr(lnL, "gradient") <- gari * weights
# }
# if (get.bi) {
# Pnjq <- lapply(Pntirq, function(x) apply(x, c(1, 3),
# mean))
# if (panel)
# Wnq <- Wnq[id, ]
# Pw <- lapply(Pnjq, function(x) x * Wnq)
# attr(lnL, "prob.alt") <- sapply(Pw, function(x) apply(x,
# 1, sum))
# attr(lnL, "prob.ind") <- Ln
# attr(lnL, "bi") <- bi
# attr(lnL, "Qir") <- Qir
# attr(lnL, "Wnq") <- Wnq
# }
# lnL
# },# weights = 1, R = 5, seed = 12345, ranp = c(pf = "n", cl = "n",
# loc = "n", wk = "n", tod = "n", seas = "n"), id = structure(c(1L,
# 1L, 1L, 1L, 1L), .Label = "1", class = "factor"), H = list(
# `1` = structure(0, .Dim = c(1L, 1L), .Dimnames = list(
# "1", "(class)2")), `2` = structure(1, .Dim = c(1L,
# 1L), .Dimnames = list("2", "(class)2"))), correlation = FALSE,
# haltons = NA, Q = 2)
# 11: eval(f, sys.frame(sys.parent()))
# 10: eval(f, sys.frame(sys.parent()))
# 9: callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
# ...)
# 8: (function (theta, fName, ...)
#
# 7: do.call(callWithoutSumt, argList)
# 6: maxOptim(fn = fn, grad = grad, hess = hess, start = start, method = "BFGS",
# fixed = fixed, constraints = constraints, finalHessian = finalHessian,
# parscale = parscale, control = mControl, ...)
# 5: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
# constraints = constraints, ...)
# 4: maxLik(method = "bfgs", iterlim = 500, start = c(`class.1.2:(intercept)` = -4.85114128700713,
# `class.1.3:(intercept)` = -7.69322200825539, `class.1.4:(intercept)` = 5.01582959989182,
# class.1.pf = -1.60963678008691, class.1.cl = 0.109892050051351,
# class.1.loc = 18.3461318629584, class.1.wk = 5.01552145983325,
# class.1.tod = 6.12905713997904, class.1.seas = -4.37562129235275,
# `class.2.2:(intercept)` = -4.81114128700713, `class.2.3:(intercept)` = -7.6532220082554,
# `class.2.4:(intercept)` = 5.05582959989182, class.2.pf = -1.56963678008691,
# class.2.cl = 0.149892050051351, class.2.loc = 18.3861318629584,
# class.2.wk = 5.05552145983325, class.2.tod = 6.16905713997903,
# class.2.seas = -4.33562129235275, class.1.sd.pf = 0.08, class.1.sd.cl = 0.08,
# class.1.sd.loc = 0.08, class.1.sd.wk = 0.08, class.1.sd.tod = 0.08,
# class.1.sd.seas = 0.08, class.2.sd.pf = 0.12, class.2.sd.cl = 0.12,
# class.2.sd.loc = 0.12, class.2.sd.wk = 0.12, class.2.sd.tod = 0.12,
# class.2.sd.seas = 0.12, `(class)2` = 0), X = Xl, y = yl, gradient = gradient,
# weights = weights, logLik = ll.mnlogit, R = R, seed = seed,
# ranp = ranp, id = id, H = Hl, correlation = correlation,
# haltons = haltons, Q = Q)
# 3: eval(opt, sys.frame(which = nframe))
# 2: eval(opt, sys.frame(which = nframe))
# 1: gmnl(choice ~ pf + cl + loc + wk + tod + seas | 1 | 0 | 0 | 1,
# data = Electr, subset = 1:20, model = "mm", R = 5, panel = TRUE,
# ranp = c(pf = "n", cl = "n", loc = "n", wk = "n", tod = "n",
# seas = "n"), Q = 2, iterlim = 500)
奖金2：sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)
Matrix products: default
attached base packages:
[1] grid stats graphics grDevices utils datasets
[7] methods base
other attached packages:
[1] here_1.0.1 strucchange_1.5-2 sandwich_3.0-1
[4] zoo_1.8-9 partykit_1.2-15 mvtnorm_1.1-3
[7] libcoin_1.0-9 mlogit_1.1-1 dfidx_0.0-4
[10] gmnl_1.1-3.2 Formula_1.2-4 maxLik_1.5-2
[13] miscTools_0.6-26 dplyr_1.0.7 nnet_7.3-17
先感谢您。