1

我正在尝试使用多项逻辑回归模型，其中公式或线性预测变量对于三个结果之一不同。

这是一个示例数据集。抱歉，创建数据集的代码有点长：

my.data <- read.table(text = '
  obs  cov  cov2 n.a  n.b  n.c
    1   -7   49   40   60    0
    2   -6   36   40   60    0
    3   -5   25   40   60    0
    4   -4   16   40   60    0
    5   -3    9   40   59    1
    6   -2    4   40   57    3
    7   -1    1   40   47   13
    8    0    0   40   27   33
    9    1    1   40    9   51
   10    2    4   40    2   58
   11    3    9   40    1   59
   12    4   16   40    0   60
   13    5   25   40    0   60
   14    6   36   40    0   60
   15    7   49   40    0   60
', header = TRUE, stringsAsFactors = FALSE)

# duplicate rows
n.times <- my.data$n.a
data.a  <- my.data[rep(seq_len(nrow(my.data)), n.times),]
data.a$stage <- 'a'
n.times <- my.data$n.b
data.b  <- my.data[rep(seq_len(nrow(my.data)), n.times),]
data.b$stage <- 'b'
n.times <- my.data$n.c
data.c  <- my.data[rep(seq_len(nrow(my.data)), n.times),]
data.c$stage <- 'c'

# combine data sets
my.data <- rbind(data.a, data.b)
my.data <- rbind(my.data, data.c)

my.data <- my.data[order(my.data$cov, my.data$stage),]
head(my.data)
dim(my.data)

nnet下面是用包和包创建模型的代码mlogit： 在这个模型阶段b和c用相同的公式（截距cov和cov2）建模。阶段a是参考。这两个包返回的估计值非常相似。

# first with package nnet

library(nnet)

my.data$stage  <- as.factor(my.data$stage)
my.data$stage2 <- relevel(my.data$stage, ref = "a")
model1 <- multinom(stage2 ~ cov + cov2, data = my.data)
summary(model1)

#
# Call:
# multinom(formula = stage2 ~ cov + cov2, data = my.data)
#
# Coefficients:
#   (Intercept)        cov       cov2
# b  -0.7180498 -0.6390276 -0.0735323
# c  -0.5639989  0.5903990 -0.0701099
#
# Std. Errors:
#   (Intercept)        cov        cov2
# b   0.1191425 0.06643554 0.010191801
# c   0.1109950 0.05976451 0.009468451
#
# Residual Deviance: 2301.073 
# AIC: 2313.073 
#

fitted(model1)[1:10,]

# now with package mlogit 

library(mlogit)

my.datad <- my.data

my.datad <- my.data[,c('stage', 'cov', 'cov2')]
rownames(my.datad) <- NULL
head(my.datad)

my.datae <- mlogit.data(my.datad, shape = "wide", choice = "stage")
head(my.datae)

summary(mlogit(stage ~ 0 | cov + cov2, data = my.datae))

#
# Call:
# mlogit(formula = stage ~ 0 | cov + cov2, data = my.datae, method = "nr", 
#     print.level = 0)
#
# Frequencies of alternatives:
#       a       b       c 
# 0.40000 0.29467 0.30533 
#
# nr method
# 8 iterations, 0h:0m:0s 
# g'(-H)^-1g = 8.63E-06 
# successive function values within tolerance limits 
#
# Coefficients :
#                 Estimate Std. Error t-value  Pr(>|t|)    
# b:(intercept) -0.7189757  0.1192246 -6.0304 1.635e-09 ***
# c:(intercept) -0.5634641  0.1109489 -5.0786 3.802e-07 ***
# b:cov         -0.6398978  0.0665175 -9.6200 < 2.2e-16 ***
# c:cov          0.5898187  0.0597128  9.8776 < 2.2e-16 ***
# b:cov2        -0.0736489  0.0102012 -7.2197 5.211e-13 ***
# c:cov2        -0.0700294  0.0094624 -7.4008 1.352e-13 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Log-Likelihood: -1150.5
# McFadden R^2:  0.29554 
# Likelihood ratio test : chisq = 965.34 (p.value = < 2.22e-16)
#

但是，我想做的是使用 stageb作为参考，将 model stagec作为一个 intercept 的函数，cov如上所述cov2，但 model stagea只是作为一个 intercept 的函数。请注意，在数据集中，协变量不会影响在阶段结束的试验次数：无论协变量的值如何，a40 次试验都在阶段结束。a

这样的模型可能吗？我相信是这样，但我无法弄清楚如何使用这些软件包中的任何一个来做到这一点。我尝试使用各种指标变量从阶段公式中删除协变量，a但无论如何总是估计系数并且标准误差变得巨大。有时点估计也变得非常大。

我在问一个相关的问题Cross Validated，但我认为这个问题主要是关于编程的。如果有兴趣，这是我关于交叉验证的相关问题的链接：

https://stats.stackexchange.com/questions/183427/modeling-probability-with-the-multinomial-logit-link/184004#184004

谢谢你的任何建议。

编辑 2015 年 11 月 30 日

我现在已经从另外两个软件程序中获得了估计。这些估计值是我希望从中看到的可能目标值R。虽然，我怀疑最终可能会有更好的估计。

来自一个应用程序的估计：

Parameter       Beta           SE    Lower 95%CI  Upper 95%CI
state a: B0   0.305620     0.062682     0.182764     0.428476
state c: B0  -0.094760     0.113606    -0.317428     0.127908
state c: B1   0.750266     0.038993     0.673841     0.826692
state d: B2  -0.085494     0.012216    -0.109437    -0.061551

来自第二个应用程序的估计：

Parameter       Beta           SE      Lower 95%CI   Upper 95%CI
state a: B0   0.3056197    0.0626826    0.1827618    0.4284777
state c: B0  -0.0947603    0.1124746   -0.3152105    0.1256900
state c: B1   0.7502663    0.0601626    0.6323476    0.8681850     
state c: B2  -0.0854941    0.0095836   -0.1042780   -0.0667102

编辑二 2015 年 11 月 30 日

如果我对两个状态a和c两个协变量进行建模，我会从两个R包和其他两个软件应用程序中获得以下信息：

#
# model data with stage 'b' as reference
#
# model stage 'a' as function of intercept, cov and cov2
# model stage 'c' as function of intercept, cov and cov2
#
# model: a(cov, cov2) c(cov1, cov2)
#
# Parameter        Beta       SE      95%CI Lower     95%CI Upper
#
#    1:        0.1555116   0.1390947  -0.1171141      0.4281373     
#    2:        0.7189753   0.1192245   0.4852953      0.9526554     
#    3:        1.2297161   0.0853667   1.0623974      1.3970347     
#    4:        0.0036194   0.0147607  -0.0253116      0.0325505     
#    5:        0.6398974   0.0665175   0.5095231      0.7702717     
#    6:        0.0736488   0.0102012   0.0536545      0.0936431     
#

library(nnet)
my.data2 <- my.data

my.data2$stage  <- as.factor(my.data2$stage)
my.data2$stage2 <- relevel(my.data2$stage, ref = "b")

model1.nnet <- multinom(stage2 ~ cov + cov2, data = my.data2)
summary(model1.nnet)

# Call:
# multinom(formula = stage2 ~ cov + cov2, data = my.data2)
#
# Coefficients:
#   (Intercept)       cov        cov2
# a   0.7189754 0.6398974 0.073648810
# c   0.1555120 1.2297159 0.003619449
#
# Std. Errors:
#   (Intercept)        cov       cov2
# a   0.1192246 0.06651748 0.01020116
# c   0.1390947 0.08536677 0.01476072
#
# Residual Deviance: 2301.073 
# AIC: 2313.073 

library(mlogit)
my.data2b <- my.data2[,c('stage', 'cov', 'cov2')]
rownames(my.data2b) <- NULL
head(my.data2b)

my.data2.mlogit <- mlogit.data(my.data2b, shape = "wide", choice = "stage")
head(my.data2.mlogit)

summary(mlogit(stage ~ 0 | cov + cov2, data = my.data2.mlogit, reflevel = "b"))

# Coefficients :
#               Estimate Std. Error t-value  Pr(>|t|)    
# a:(intercept) 0.7189757  0.1192246  6.0304 1.635e-09 ***
# c:(intercept) 0.1555116  0.1390948  1.1180    0.2636    
# a:cov         0.6398978  0.0665175  9.6200 < 2.2e-16 ***
# c:cov         1.2297166  0.0853668 14.4051 < 2.2e-16 ***
# a:cov2        0.0736489  0.0102012  7.2197 5.211e-13 ***
# c:cov2        0.0036195  0.0147607  0.2452    0.8063    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#

但是，如果我尝试a仅使用截距对状态进行建模，我仍然没有得到与R其他两个应用程序中的任何一个包类似的估计值：

#
# model data with stage 'b' as reference
#
# model stage 'a' as function of intercept only
# model stage 'c' as function of intercept, cov and cov2
#
#  Parameter        Beta           SE     95%CI Lower  95%CI Upper
#
#  stage a: B0    0.305620     0.062682     0.182764     0.428476
#  state c: B0   -0.094760     0.113606    -0.317428     0.127908
#  state c: B1    0.750266     0.038993     0.673841     0.826692
#  state c: B2   -0.085494     0.012216    -0.109437    -0.061551
#

library(nnet)

my.data3 <- my.data

my.data3$stage  <- as.factor(my.data3$stage)
my.data3$stage2 <- relevel(my.data3$stage, ref = "b")

my.data3$cov  <- ifelse(my.data3$stage == 'a', 0, my.data3$cov )
my.data3$cov2 <- ifelse(my.data3$stage == 'a', 0, my.data3$cov2)

model2.nnet <- multinom(stage2 ~ cov + cov2, data = my.data3)
summary(model2.nnet)

#     Call:
# multinom(formula = stage2 ~ cov + cov2, data = my.data3)
#
# Coefficients:
#   (Intercept)       cov         cov2
# a   3.1129805 0.5936333 -13.85909619
# c   0.2221975 1.5220859  -0.01343098
#
# Std. Errors:
#   (Intercept)        cov        cov2
# a   0.1694357 33.9858262 33.98601992
# c   0.1834233  0.1339483  0.06296883
#
# Residual Deviance: 661.0351 
# AIC: 673.0351 

library(mlogit)

my.data3b <- my.data3[,c('stage', 'cov', 'cov2')]
rownames(my.data3b) <- NULL
head(my.data3b)

my.data3.mlogit <- mlogit.data(my.data3b, shape = "wide", choice = "stage")
head(my.data3.mlogit)

summary(mlogit(stage ~ 0 | cov + cov2, data = my.data3.mlogit, reflevel = "b"))

# Coefficients :
#                 Estimate  Std. Error t-value Pr(>|t|)    
# a:(intercept)    3.112970    0.169436 18.3726   <2e-16 ***
# c:(intercept)    0.222162    0.183426  1.2112   0.2258    
# a:cov            0.829259 2276.499314  0.0004   0.9997    
# c:cov            1.522129    0.133954 11.3631   <2e-16 ***
# a:cov2         -22.295201 2276.499317 -0.0098   0.9922    
# c:cov2          -0.013431    0.062973 -0.2133   0.8311    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
4

2 回答 2

0

下面是使用when 阶段R估计参数的代码，每个阶段都与协变量相关，并且当阶段仅使用截距建模时。optimaccovcov2a

鉴于我现在已经能够a用三种不同的截距来建模阶段，我不清楚为什么我不能用mlogitornnet R包获得相同的估计。

首先像以前一样创建数据集：

my.data <- read.table(text = '
  obs  cov  cov2 n.a  n.b  n.c
    1   -7   49   40   60    0
    2   -6   36   40   60    0
    3   -5   25   40   60    0
    4   -4   16   40   60    0
    5   -3    9   40   59    1
    6   -2    4   40   57    3
    7   -1    1   40   47   13
    8    0    0   40   27   33
    9    1    1   40    9   51
   10    2    4   40    2   58
   11    3    9   40    1   59
   12    4   16   40    0   60
   13    5   25   40    0   60
   14    6   36   40    0   60
   15    7   49   40    0   60
', header = TRUE, stringsAsFactors = FALSE)

# duplicate rows
n.times.a <- my.data$n.a
data.a  <- my.data[rep(seq_len(nrow(my.data)), n.times.a),]
data.a$stage <- 'a'
n.times.b <- my.data$n.b
data.b  <- my.data[rep(seq_len(nrow(my.data)), n.times.b),]
data.b$stage <- 'b'
n.times.c <- my.data$n.c
data.c  <- my.data[rep(seq_len(nrow(my.data)), n.times.c),]
data.c$stage <- 'c'

# combine data sets
my.data <- rbind(data.a, data.b)
my.data <- rbind(my.data, data.c)

my.data <- my.data[order(my.data$cov, my.data$stage),]

# Here are a few additional lines to prepare the data set for my `optim` functions.

cov  <- my.data$cov
cov2 <- my.data$cov2

n.a  <- ifelse(my.data$stage == 'a', 1, 0)
n.b  <- ifelse(my.data$stage == 'b', 1, 0)
n.c  <- ifelse(my.data$stage == 'c', 1, 0)

这是多项逻辑回归的代码，它返回与and包和其他两个软件应用程序optim相同的估计值（即，stage和每个都使用截距和and效果建模）：mlogitnnetaccov1cov2

my.function <- function(betas, cov, cov2, n.a, n.b, n.c){

     b0a = betas[1]
     b1a = betas[2]
     b2a = betas[3]
     b0c = betas[4]
     b1c = betas[5]
     b2c = betas[6]

     n = nrow(my.data)

     llh = 0

     for(i in 1:n){

          y <- (

                (n.b[i] * (1 - exp(b0a + b1a * cov[i] + b2a * cov2[i]) / 
                          (1 + exp(b0a + b1a * cov[i] + b2a * cov2[i]) + exp(b0c + b1c * cov[i] + b2c * cov2[i]))    - 
                               exp(b0c + b1c * cov[i] + b2c * cov2[i]) / 
                          (1 + exp(b0a + b1a * cov[i] + b2a * cov2[i]) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) )) +

                (n.c[i] * (    exp(b0c + b1c * cov[i] + b2c * cov2[i]) / 
                          (1 + exp(b0a + b1a * cov[i] + b2a * cov2[i]) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) )) +

                (n.a[i] * (    exp(b0a + b1a * cov[i] + b2a * cov2[i]) / 
                          (1 + exp(b0a + b1a * cov[i] + b2a * cov2[i]) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) ))

               )

          y <- log(y)

          y <- ifelse(is.na(y), 0.0000000001, y)

          llh = llh + y

     }

     -1 * llh

}

Nstar <- optim(c(0,0,0,0,0,0), my.function, cov = cov, cov2 = cov2, n.a = n.a, n.b = n.b, n.c = n.c, method = "BFGS", hessian = TRUE)
Nstar$par

# [1] 0.718951850 0.639832930 0.073637858 0.155471765 1.229635652 0.003612455

这是当阶段使用截距和效果建模optim时的多项逻辑回归代码，但阶段仅使用截距建模。返回的估计值与我使用其他两个软件应用程序获得的估计值相匹配，但与使用以下或中的包获得的估计值不匹配：ccov1cov2amlogitnnetR

my.other.function <- function(betas, cov, cov2, n.a, n.b, n.c){

     b0a = betas[1]
     b0c = betas[2]
     b1c = betas[3]
     b2c = betas[4]

     n = nrow(my.data)

     llh = 0

     for(i in 1:n){

          y <- (

                (n.b[i] * (1 - exp(b0a                               ) / (1 + exp(b0a) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) - 
                               exp(b0c + b1c * cov[i] + b2c * cov2[i]) / (1 + exp(b0a) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) )) +

                (n.c[i] * (    exp(b0c + b1c * cov[i] + b2c * cov2[i]) / (1 + exp(b0a) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) )) +

                (n.a[i] * (    exp(b0a                               ) / (1 + exp(b0a) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) ))


               )

          y <- log(y)

          y <- ifelse(is.na(y), 0.0000000001, y)

          llh = llh + y

     }

     -1 * llh

}

Nstar <- optim(c(0,0,0,0), my.other.function, cov = cov, cov2 = cov2, n.a = n.a, n.b = n.b, n.c = n.c, method = "BFGS", hessian = TRUE)
Nstar$par

# [1]  0.30561794 -0.09473753  0.75021769 -0.08548674

也许我采用的方法存在根本性的错误，optim这解释了为什么mlogitandnnet包不允许创建这种模型结构？或者也许我只是还没有弄清楚与mlogitandnnet包一起使用的正确语法？

我可能需要提取和研究mlogitandnnet包正在使用的源代码，以查看是否可以修改它，或者至少弄清楚当我尝试a仅使用截取来建模阶段时它在做什么。

如果我弄清楚如何使用or （或）包a仅使用截距来建模阶段，那么我将发布更新。mlogitnnetmnlogitR

编辑：2015 年 12 月 7 日

我现在已经能够用来optim重现由mlogit. R代码如下。结论是，到目前为止，我使用的三种方法都涉及我修改设计矩阵以从 stage 中删除协变量a。简单地将协变量的数据设置为0不会从设计矩阵中删除这些协变量。

cov  <- ifelse(my.data$stage == 'a', 0, cov )
cov2 <- ifelse(my.data$stage == 'a', 0, cov2)

my.third.function <- function(betas, cov, cov2, n.a, n.b, n.c){

     b0a = betas[1]
     b1a = betas[2]
     b2a = betas[3]
     b0c = betas[4]
     b1c = betas[5]
     b2c = betas[6]

     n = nrow(my.data)

     llh = 0

     for(i in 1:n){

          y <- (

                (n.b[i] * (1 - exp(b0a + b1a * cov[i] + b2a * cov2[i]) / 
                          (1 + exp(b0a + b1a * cov[i] + b2a * cov2[i]) + exp(b0c + b1c * cov[i] + b2c * cov2[i]))    - 
                               exp(b0c + b1c * cov[i] + b2c * cov2[i]) / 
                          (1 + exp(b0a + b1a * cov[i] + b2a * cov2[i]) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) )) +

                (n.c[i] * (    exp(b0c + b1c * cov[i] + b2c * cov2[i]) / 
                          (1 + exp(b0a + b1a * cov[i] + b2a * cov2[i]) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) )) +

                (n.a[i] * (    exp(b0a + b1a * cov[i] + b2a * cov2[i]) / 
                          (1 + exp(b0a + b1a * cov[i] + b2a * cov2[i]) + exp(b0c + b1c * cov[i] + b2c * cov2[i])) ))

               )

#         y <- ifelse(is.na(y) | y <= 0, 0.0000000001, y)

          y <- log(y)

          llh = llh + y

     }

     -1 * llh

}

model3 <- optim(c(0,0,0,0,0,0), my.third.function, cov = cov, cov2 = cov2, n.a = n.a, n.b = n.b, n.c = n.c, method = "BFGS", hessian = TRUE)
model3$par

#
# [1]   3.11296505   0.61033815 -13.89223292   0.22214130   1.52209746  -0.01344045
#

我还剖析了mlogit包的源代码。到目前为止，我已经能够从该源代码中的设计矩阵中删除协变量，但只是这样做并不会返回正确的估计值。我对设计矩阵的更改一定会导致源代码稍后出现错误。

如果我能够修改源代码的其余部分以返回正确的估计值，或者如果我能够找出删除mlogit原始帖子中声明中协变量的正确语法，我将发布更新。

于 2015-12-01T10:33:26.883 回答
0

在我看来，解决这个问题的一个好方法是将其分解为两个模型。您需要 stage = a 独立于协变量的概率。然后你想知道，给定阶段！= a，概率阶段 = b 或 c 取决于协变量。

#pr(stage=a)
my.data$stageA.BC = my.data$stage=="a"
glm(my.data$stageA.BC ~ 1,family=binomial)



#pr(stage=c|cov,cov2,stage!= a)
my.data.BC = my.data[my.data$stageA.BC==0,]
my.data.BC = relevel(my.data.BC$stage,ref="b")
glm(stage ~cov + cov2, data=my.data.BC,family=binomial)

由于 pr(stage = b OR c) = 1 - pr(stage=a)，那么您将拥有：

pr(stage = a) 

pr(stage = b) = (1 - pr(stage = a)) * pr(stage=b|cov,cov2,stage!= a)

pr(stage = c) = (1 - pr(stage = a)) * pr(stage=c|cov,cov2,stage!= a)
于 2015-11-30T15:41:19.523 回答