多元logit回归参数估计（多分类logit回归预测）

最新推荐文章于 2025-07-29 15:18:08 发布

原创最新推荐文章于 2025-07-29 15:18:08 发布 · 1.3w 阅读

44 ·

CC 4.0 BY-SA版权

文章标签：

#预测 #class

数据挖掘算法专栏收录该内容

68 篇文章

订阅专栏

本文介绍了多元离散选择模型的基础知识，特别是多元Logit模型，并通过R语言提供了一个详尽的样例程序来演示如何实现该模型。此外，还讨论了在处理大规模数据时可能遇到的问题，并给出了解决方案。

多元离散选择模型简介
样例程序
样例程序的问题
问题解决

多元离散选择模型简介

常用的离散选择模型有logit模型及probit模型，其区别就是假设不可观测量的分布不同。logit模型假设不可观测项服从Gumbel (Extreme Value Type I) Distribution 。多元logit模型是，多元离散选择模型的一种，适合于效用最大化时的分布选择。
如果决策者i在J项可供选择方案中选择了第j项，那么其效用模型为
这里写图片描述

样例程序

以下是基于R语言的一个程序

#引入所需要的包
require(foreign)
require(nnet)
require(ggplot2)
require(reshape2)
##读取数据集，这里数据集可以来找csv，mysql等
ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")
head(ml)  ## 查看数据的前六行
with(ml, table(ses, prog))  ## 以表格的形式统计数目
##函数tapply（进行分组统计）
with(ml, do.call(rbind, tapply(write, prog,function(x) c(M = mean(x), SD = sd(x)))))
## 重新排列因子水平，其中以academic作为基础标准
ml$prog2 <- relevel(ml$prog, ref = "academic")
ml$prog2
##训练模型
test <- multinom(prog2 ~ ses + write, data = ml)
summary(test)
#获取z值
z <- summary(test)$coefficients/summary(test)$standard.errors
z
#获取p值
p <- (1 - pnorm(abs(z), 0, 1))*2
p
## extract the coefficients from the model and exponentiate
exp(coef(test))
head(pp <- fitted(test))
##测试集上检测
dses <- data.frame(ses = c("low", "middle", "high"),write = mean(ml$write))
#获取预测的概率
predict(test, newdata = dses, "probs")
dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41),write = rep(c(30:70), 3))
## store the predicted probabilities for each value of ses and write
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))
## calculate the mean probabilities within each level of ses
by(pp.write[, 3:5], pp.write$ses, colMeans)
## melt data set to long for ggplot2
lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")
head(lpp) # view first few rows
## plot predicted probabilities across write values for
## each level of ses facetted by program type
ggplot(lpp, aes(x = write, y = probability, colour = ses)) +
  geom_line() +
  facet_grid(variable ~ ., scales="free")

样例程序的问题

但是，在使用上面的样例，如下，当X的维数较大时，类别较多时，就会出现问题，其主要问题是迭代时权重的问题，为此，我们需要修改源码。

test <- multinom(trainset$car_id2 ~ ., data = trainset[,-c(1,13,14)])

以下是容易出现的问题：

Error in nnet.default(X, Y, w, mask = mask, size = 0, skip = TRUE, softmax = TRUE,  : 
  too many (11745) weights

问题解决

针对问题too many (11745) weights，这类问题，解决的办法是修改源码。即新写一个方法。如下，我修改的部分是这里面的

MaxNWts=12000 ##将程序中的MaxNWts变大就行了

##估计参数的函数
nnetfunction<-function (formula, data, weights, subset, na.action, contrasts = NULL, 
                        Hess = FALSE, summ = 0, censored = FALSE, model = FALSE, 
                        ...) 
{
  class.ind <- function(cl) {
    n <- length(cl)
    x <- matrix(0, n, length(levels(cl)))
    x[(1L:n) + n * (as.integer(cl) - 1L)] <- 1
    dimnames(x) <- list(names(cl), levels(cl))
    x
  }
  summ2 <- function(X, Y) {
    X <- as.matrix(X)
    Y <- as.matrix(Y)
    n <- nrow(X)
    p <- ncol(X)
    q <- ncol(Y)
    Z <- t(cbind(X, Y))
    storage.mode(Z) <- "double"
    z <- .C(VR_summ2, as.integer(n), as.integer(p), as.integer(q), 
            Z = Z, na = integer(1L))
    Za <- t(z$Z[, 1L:z$na, drop = FALSE])
    list(X = Za[, 1L:p, drop = FALSE], Y = Za[, p + 1L:q])
  }
  call <- match.call()
  m <- match.call(expand.dots = FALSE)
  m$summ <- m$Hess <- m$contrasts <- m$censored <- m$model <- m$... <- NULL
  m[[1L]] <- quote(stats::model.frame)
  m <- eval.parent(m)
  Terms <- attr(m, "terms")
  X <- model.matrix(Terms, m, contrasts)
  cons <- attr(X, "contrasts")
  Xr <- qr(X)$rank
  Y <- model.response(m)
  if (!is.matrix(Y)) 
    Y <- as.factor(Y)
  w <- model.weights(m)
  if (length(w) == 0L) 
    if (is.matrix(Y)) 
      w <- rep(1, dim(Y)[1L])
  else w <- rep(1, length(Y))
  lev <- levels(Y)
  if (is.factor(Y)) {
    counts <- table(Y)
    if (any(counts == 0L)) {
      empty <- lev[counts == 0L]
      warning(sprintf(ngettext(length(empty), "group %s is empty", 
                               "groups %s are empty"), paste(sQuote(empty), 
                                                             collapse = " ")), domain = NA)
      Y <- factor(Y, levels = lev[counts > 0L])
      lev <- lev[counts > 0L]
    }
    if (length(lev) < 2L) 
      stop("need two or more classes to fit a multinom model")
    if (length(lev) == 2L) 
      Y <- as.integer(Y) - 1
    else Y <- class.ind(Y)
  }
  if (summ == 1) {
    Z <- cbind(X, Y)
    z1 <- cumprod(apply(Z, 2L, max) + 1)
    Z1 <- apply(Z, 1L, function(x) sum(z1 * x))
    oZ <- order(Z1)
    Z2 <- !duplicated(Z1[oZ])
    oX <- (seq_along(Z1)[oZ])[Z2]
    X <- X[oX, , drop = FALSE]
    Y <- if (is.matrix(Y)) 
      Y[oX, , drop = FALSE]
    else Y[oX]
    w <- diff(c(0, cumsum(w))[c(Z2, TRUE)])
    print(dim(X))
  }
  if (summ == 2) {
    Z <- summ2(cbind(X, Y), w)
    X <- Z$X[, 1L:ncol(X)]
    Y <- Z$X[, ncol(X) + 1L:ncol(Y), drop = FALSE]
    w <- Z$Y
    print(dim(X))
  }
  if (summ == 3) {
    Z <- summ2(X, Y * w)
    X <- Z$X
    Y <- Z$Y[, 1L:ncol(Y), drop = FALSE]
    w <- rep(1, nrow(X))
    print(dim(X))
  }
  offset <- model.offset(m)
  r <- ncol(X)
  if (is.matrix(Y)) {
    p <- ncol(Y)
    sY <- Y %*% rep(1, p)
    if (any(sY == 0)) 
      stop("some case has no observations")
    if (!censored) {
      Y <- Y/matrix(sY, nrow(Y), p)
      w <- w * sY
    }
    if (length(offset) > 1L) {
      if (ncol(offset) != p) 
        stop("ncol(offset) is wrong")
      mask <- c(rep(FALSE, r + 1L + p), rep(c(FALSE, rep(TRUE, 
                                                         r), rep(FALSE, p)), p - 1L))
      X <- cbind(X, offset)
      Wts <- as.vector(rbind(matrix(0, r + 1L, p), diag(p)))
      fit <- nnet(X, Y, w, Wts = Wts, mask = mask, 
                  size = 0, skip = TRUE, softmax = TRUE, censored = censored, 
                  rang = 0,MaxNWts=12000, ...)
    }
    else {
      mask <- c(rep(FALSE, r + 1L), rep(c(FALSE, rep(TRUE, 
                                                     r)), p - 1L))
      fit <- nnet(X, Y, w, mask = mask, size = 0, 
                  skip = TRUE, softmax = TRUE, censored = censored, 
                  rang = 0,MaxNWts=12000, ...)
    }
  }
  else {
    if (length(offset) <= 1L) {
      mask <- c(FALSE, rep(TRUE, r))
      fit <- nnet(X, Y, w, mask = mask, size = 0, 
                  skip = TRUE, entropy = TRUE, rang = 0,MaxNWts=12000, ...)
    }
    else {
      mask <- c(FALSE, rep(TRUE, r), FALSE)
      Wts <- c(rep(0, r + 1L), 1)
      X <- cbind(X, offset)
      fit <- nnet(X, Y, w, Wts = Wts, mask = mask, 
                  size = 0, skip = TRUE, entropy = TRUE, rang = 0, MaxNWts=12000, 
                  ...)
    }
  }
  fit$formula <- attr(Terms, "formula")
  fit$terms <- Terms
  fit$call <- call
  fit$weights <- w
  fit$lev <- lev
  fit$deviance <- 2 * fit$value
  fit$rank <- Xr
  edf <- ifelse(length(lev) == 2L, 1, length(lev) - 1) * Xr
  if (is.matrix(Y)) {
    edf <- (ncol(Y) - 1) * Xr
    if (length(dn <- colnames(Y)) > 0) 
      fit$lab <- dn
    else fit$lab <- 1L:ncol(Y)
  }
  fit$coefnames <- colnames(X)
  fit$vcoefnames <- fit$coefnames[1L:r]
  fit$na.action <- attr(m, "na.action")
  fit$contrasts <- cons
  fit$xlevels <- .getXlevels(Terms, m)
  fit$edf <- edf
  fit$AIC <- fit$deviance + 2 * edf
  if (model) 
    fit$model <- m
  class(fit) <- c("multinom", "nnet")
  if (Hess) 
    fit$Hessian <- multinomHess(fit, X)
  fit
}