在R中的多项逻辑多层次模型

Question

在R中的多项逻辑多层次模型

49

问题：我需要估计一组多项Logistic混合效应模型，但找不到合适的R包。哪个R包是用于估算这种模型的最佳选择？STATA 13最近添加了这个功能到他们的混合效应模型中——所以估算这种模型的技术似乎已经可用。

细节：许多研究问题需要估算多项Logistic回归模型，其中结果变量是分类的。例如，生物学家可能有兴趣调查哪种类型的树（例如松树、枫树、橡树）最受酸雨的影响。市场研究人员可能会想知道客户年龄与在Target、Safeway或Walmart的购物频率之间是否存在关系。这些情况的共同之处在于结果变量是分类的（无序的），而多项Logistic回归是首选的估算方法。在我的情况下，我正在研究人类迁移类型的差异，结果变量（mig）编码为0=未迁移，1=内部迁移，2=国际迁移。这是我的数据集的简化版本：

migDat=data.frame(hhID=1:21,mig=rep(0:2,times=7),age=ceiling(runif(21,15,90)),stateID=rep(letters[1:3],each=7),pollution=rep(c("high","low","moderate"),each=7),stringsAsFactors=F)

   hhID mig age stateID pollution
1     1   0  47       a      high
2     2   1  53       a      high
3     3   2  17       a      high
4     4   0  73       a      high
5     5   1  24       a      high
6     6   2  80       a      high
7     7   0  18       a      high
8     8   1  33       b       low
9     9   2  90       b       low
10   10   0  49       b       low
11   11   1  42       b       low
12   12   2  44       b       low
13   13   0  82       b       low
14   14   1  70       b       low
15   15   2  71       c  moderate
16   16   0  18       c  moderate
17   17   1  18       c  moderate
18   18   2  39       c  moderate
19   19   0  35       c  moderate
20   20   1  74       c  moderate
21   21   2  86       c  moderate

我的目标是估计年龄（自变量）对于以下情况的概率影响：（1）内部迁移与不迁移，（2）国际迁移与不迁移，（3）内部迁移与国际迁移。另一个复杂因素在于我的数据在不同的聚合级别上运作（例如，污染在州一级运作），我也有兴趣预测空气污染（pollution）对于采取特定类型运动的概率影响。

笨拙的解决方案：一个方法是通过将每个模型的数据集减少到仅包含两种迁移类型来估计一组不同的逻辑回归模型（例如，模型1：只有编码为mig=0和mig=1的情况；模型2：只有编码为mig=0和mig=2的情况；模型3：只有编码为mig=1和mig=2的情况）。虽然可以使用lme4估计这样一个简单的多层逻辑回归模型，但这种方法并不完全理想，因为它没有适当地考虑省略的情况的影响。第二种解决方法是使用R中的R2MLwiN包通过MLWiN运行多项逻辑多层模型。但由于MLWiN不是开源的，且生成的对象难以使用，我更喜欢避免这个选项。根据全面的互联网搜索，似乎有一些对这样的模型的需求，但我不知道是否有一个好的R包。因此，如果已经运行过这样的模型的专家可以提供建议，如果有多个包，也许可以指出一些优缺点。我相信这样的信息会成为多个R用户非常有用的资源。谢谢！

最好的祝福，拉斐尔

- Raphael

13

两个建议：(1) 研究一下“MCMCglmm”包；(2) 你所说的“笨拙方法”其实是标准方法（例如，参见Dobson和Barnett的《广义线性模型导论》第3版），即将多项式模型参数化为一系列二项式对比（水平1与水平2，水平1与水平3），并拟合一系列模型。这实际上是一个完整的模型，因为多项式模型的任何两类子集都是有条件的二项式（即如果您知道它是A或B，则A是从(A+B)中抽取的二项式样本）；任何完整的一组对都是有效的参数化。 - Ben Bolker

2

在您的情况下，由于您的类别有点有序，我可能会将它们参数化为（无迁移 vs. 内部或国际迁移），（内部 vs. 国际迁移）；这也为您提供了与序数模型进行比较的机会（请参见“ordinal”包）。 - Ben Bolker

非常感谢Ben Bolker！两个建议确实非常有帮助，我会进一步探索它们。 - Raphael

7个回答

22

虽然这是一个老问题，但我认为最近出现了一种可行的选择，即使用 Bayesian 的 Stan 程序来运行模型，brms 就是这样的工具。例如，如果您想在 iris 数据上运行多项逻辑回归：

b1 <- brm (Species ~ Petal.Length + Petal.Width + Sepal.Length + Sepal.Width,
           data=iris, family="categorical",
           prior=c(set_prior ("normal (0, 8)")))

为了进行序数回归——当然不适用于iris——您需要将family="categorical"更改为family="acat"（或者，根据所需的序数回归类型，更改为cratio或sratio），并确保因变量是ordered。

根据Raphael的评论澄清：此brm调用将您的公式和参数编译成Stan代码。 Stan将其编译为C ++并使用您系统的C ++编译器，这是必需的。例如，在Mac上，您可能需要安装免费的开发人员工具来获取C ++。不确定Windows如何操作。Linux应默认安装了C ++。）

根据Qaswed的评论澄清：brms也可以使用R公式(1 | groupvar)轻松处理多层模型，以添加组（随机）截距，(1 + foo | groupvar)来添加随机截距和斜率等。

- Wayne

“需要 C++” 是指我必须知道如何使用 C++ 编程，还是只需安装某个程序？此外，如果观测数据存在聚类现象，例如在 iris 数据集中按“花瓣宽度”进行聚类，那么在上述 brm 调用中如何指定随机效应/水平？ - Raphael

在Linux上安装BRM的快速说明：我遇到了依赖问题，需要PKI包。通过安装开发SSL包sudo apt-get install libssl-dev解决了这个问题。 - Simon

1

你能详细介绍一下将先验设定为normal(0,8)吗？我知道这意味着均值为0，标准差为8的正态分布，但我想知道在一个简单模型中，这是否是一个相对安全的选择，预计会表现良好。我试图查阅Gelman关于弱信息先验的指南，但这超出了我的理解范围。 - Szasulja

1

@Szasulja 我不是专家，也不能说 normal (0, 8) 是有原则的。看看 curve (dnorm (x, 0, 8), from=-20, to=20) 并注意到 qnorm (c(0.05, 0.95), 0, 8) 是 -13.15883 13.15883，这基本上给了斜率和截距足够的空间，但消除了高度不可能的值。 - Wayne

@Qaswed：我已经编辑了答案以考虑到这一点。谢谢！（brms轻松处理多层模型。） - Wayne

要注意因子标签。我遇到了一些奇异的标签“^”和“@”，出现了问题。我将它们重命名为更具通用性的名称（分别为“x”和“q”）。使用原始标签时会出现语法错误。 - tmn

2

我很困惑这种技术被描述为“标准”和“等效”，尽管它可能是一个好的实际解决方案。(我猜最好检查一下Allison和Dobson & Barnett的参考文献) 对于简单的多项式情况(没有聚类、重复测量等)，Begg and Gray (1984)建议使用k-1个二项logit相对于参考类别作为近似(虽然在许多情况下是很好的近似)到完整的多项logit。他们证明了在使用单个参考类别时会有一些效率损失，但在使用单个高频基线类别作为参考时，这种损失很小。Agresti (2002: p. 274)提供了一个例子，在一个五类别的例子中，即使基线类别占219个案例的70%以上，标准误差也会略微增加。也许这不是什么大不了的事，但我不知道在添加第二层随机性时，近似会变得更好。 参考文献 Agresti, A. (2002). Categorical data analysis. Hoboken NJ: Wiley. Begg, C. B., & Gray, R. (1984). Calculation of polychotomous logistic regression parameters using individualized regressions. Biometrika, 71(1), 11–18.

- tmn

1

我正在处理相同的问题，我发现一个可能的解决方案似乎是采用泊松（对数线性/计数）等效于多项式逻辑模型，如邮件列表中所述，这些漂亮的幻灯片或Agresti（2013年：353-356）。因此，应该可以使用lme4包中的glmer（... family = poisson）函数对数据进行一些聚合。 参考资料：
Agresti, A.（2013）分类数据分析。霍博肯，NJ：Wiley。

- non-numeric_argument

你能否简要澄清一下，为什么认为可以将泊松模型用于多项结果？（我浏览了邮件列表和幻灯片，但不完全理解这种方法）。名义结果变量上的不同值意味着完全不同的含义（例如，0=无迁移，1=国内迁移；2=国际迁移）。泊松模型将分配的值视为计数，因此我认为这种模型的系数估计是相当无意义的。 - Raphael

1

以下是关于所谓的“二项式-多项式变换”的一些参考资料。1）Baker，S.G.（1994）。多项式泊松变换。《统计学家》，43（4），495-504。2）Guimaraes，P.（2004）。理解多项式泊松变换。《Stata杂志》，4，265-273。该方法要求在预测变量中为每个唯一的解释变量组合添加一个指示变量。它对于小模型足够好，但否则会变得非常笨拙。 - tmn

1

我建议您使用“mlogit”软件包

- Anom

1

也许可以给出一个链接（http://cran.r-project.org/web/packages/mlogit/）以及该软件包使用的算法的一些描述？ - Ben Bolker

3

据我所知，mlogit软件包不允许包含随机效应或指定多层结构。因此，我不能在我的研究中使用这个软件包。 - Raphael

2

mlogit包清晰地处理（至少是两级，例如重复测量或面板数据）混合效应多项模型。其算法基于模拟最大似然方法，用于随机效用模型，并在Train（2009）中得到了详细记录。这些方法本质上与Demidenko（2005）所称的“固定样本模拟似然”相同。参考文献 Demidenko，E.（2005）。混合模型：理论和应用。Wiley-Interscience.Train，K。（2009）。使用模拟的离散选择方法（第二版）。剑桥：剑桥大学出版社。 - tmn

尝试使用bayesm::rhierMnlRwMixture。它可以处理具有（混合）高斯随机斜率和截距的简单纵向/重复测量模型。有一本书包含许多示例：Rossi，P.E.和Allenby，G.M.（2003）。贝叶斯统计和营销。Marketing Science，22（3），304-328。 - tmn

对于贝叶斯方法：请尝试使用bayesm包。[链接]（https://cran.r-project.org/web/packages/bayesm/index.html）。函数`bayesm :: rhierMnlRwMixture`将处理具有（混合）高斯随机斜率和截距的简单纵向/重复测量模型。包作者有一本带有许多示例的书籍：Rossi，P.E.，＆Allenby，G.M。（2003）。贝叶斯统计和营销。Marketing Science，22（3），304-328。 - tmn

1

自从我遇到同样的问题后，最近发现了这个问题。我发现有一个叫做ordinal的软件包，其中有一个cumulative link mixed model函数(clmm2)，看起来类似于提出的brm函数，但使用频率主义方法。

基本上，您需要设置链接函数(例如作为logit)，可以选择具有名义变量(即不符合比例几率假设的那些变量)，如果要允许具有非结构化切点，则将阈值设置为“flexible”，最后添加参数“random”以指定任何您希望具有随机效应的变量。

我还发现了一本书Multilevel Modeling Using R, W. Holmes Finch Jocelyn E. Bolin, Ken Kelley，它详细说明了如何使用第151页的函数，并提供了很好的例子。

- s.cerioli

0

这里有一个实现（不是我自己的）。我会在这份代码基础上进行工作。而且，这样你就真正知道引擎内部正在发生什么了。

http://www.nhsilbert.net/docs/rcode/multilevel_multinomial_logistic_regression.R

- Henry David Thorough

1

谢谢分享！但我想我的统计知识太有限了，无法处理像这样的原始代码。代码需要更多的注释，以便我知道他们为什么这样做，并确保没有问题/错误存在。看起来他们正在使用类似于Ben Bolker上面提到的MCMCglmm包的Monte Carlo方法，但我不是很确定... - Raphael

好的！你必须在R中进行分析吗？你也可以在Stata中进行分析。否则，我想不到其他的了。 - Henry David Thorough

是的，STATA将是我的最后选择。但由于我在R中进行所有变量构建和数据预处理，因此我想坚持使用一个软件。 - Raphael

1

我可能误解了问题，但为什么不将处理后的数据框以csv格式写出，并将其导入Stata？ - Henry David Thorough

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Raphael · Accepted Answer

一般有两种方法来拟合具有J个组的分类变量的多项式模型: (1) 同时估计J-1个对比度; (2) 为每个对比度估计单独的logit模型。

这两种方法会产生相同的结果吗？不，但结果通常相似。

哪种方法更好？同时拟合更精确(下面解释原因)。

为什么有人会使用单独的logit模型呢？(1) lme4包没有用于同时拟合多项式模型的例程，也没有其他可以在R中执行此操作的多级R包。因此，如果有人想在R中估计多级多项式模型，则单独的logit模型是目前唯一的实际解决方案。(2)正如一些有权威的统计学家所认为的(Begg and Gray, 1984; Allison, 1984, p. 46-47)，单独的logit模型更加灵活，因为它们允许对每个对比度的模型方程进行独立规定。

使用单独的logit模型是否合法？是的，但需要注意一些事项。这种方法被称为“Begg and Gray Approximation”。Begg和Gray(1984，p.16)表明，这种“个性化方法非常有效”。但是，存在一些效率损失，Begg和Gray Approximation会产生较大的标准误差(Agresti 2002，p.274)。因此，使用此方法更难获得显著结果，并且结果可以被认为是保守的。当参考类别很大时，这种效率损失最小(Begg and Gray, 1984; Agresti 2002)。采用Begg and Gray Approximation(非多级)的R包括mlogitBMA(Sevcikova and Raftery, 2012)。

为什么一系列单独的logit模型不够精确？在我的初始示例中，我们有一个变量(migration)，它可以有三个值A(没有迁移)、B(内部迁移)、C(国际迁移)。只有一个预测变量x(age)，则多项式模型被参数化为一系列二项对比度(如Long and Cheng, 2004 p.277所示)。

Eq. 1:  Ln(Pr(B|x)/Pr(A|x)) = b0,B|A + b1,B|A (x) 
Eq. 2:  Ln(Pr(C|x)/Pr(A|x)) = b0,C|A + b1,C|A (x)
Eq. 3:  Ln(Pr(B|x)/Pr(C|x)) = b0,B|C + b1,B|C (x)

为了这些对比，以下方程必须成立：

Eq. 4: Ln(Pr(B|x)/Pr(A|x)) + Ln(Pr(C|x)/Pr(A|x)) = Ln(Pr(B|x)/Pr(C|x))
Eq. 5: b0,B|A + b0,C|A = b0,B|C
Eq. 6: b1,B|A + b1,C|A = b1,B|C

这些方程式（方程式4-6）的问题在于，实际上它们不会完全成立，因为系数是基于略有不同的样本估计出来的，仅使用两个对比组的病例，第三组的病例被省略了。同时估计多项对比的程序确保方程式4-6成立（Long和Cheng，2004年第277页）。我不确定这种“同时”模型求解是如何工作的 - 或许有人能够提供一个解释？同时拟合多层次多项式模型的软件包括MLwiN（Steele 2013，第4页）和STATA（xlmlogit命令，Pope，2014年）。