[学习笔记] Mplus实现(多分类)Logistic回归
废话少说版
-
Logistic回归是适用于用连续变量或类别变量作为预测变量,类别变量作为结果变量的回归模型。
-
对结果变量采取logit变换,若结果变量为二分变量,变换形式为\(ln\frac{P}{1-P}\),若结果变量为多分类变量,变换形式为\(ln\frac{P(A)}{P(B)}\),\(P(A)\)和\(P(B)\)表示归属于A类别和B类别的概率。
-
预测变量中的类别变量需要转换为哑变量,再计算每个哑变量对应的回归系数。
-
回归系数\(\beta\),但通常看的是\(OR\)(称作“优势比”或“比值比”),也就是\(e^{\beta}\),表示预测变量变化一个单位,结果变量不同事件的发生比\(\frac{P}{1-P}\)或者\(\frac{P(A)}{P(B)}\)变化多少倍。
-
Mplus中,需要用
NOMINAL
语句声明结果变量为类别变量,否则会默认将其视为连续变量处理。 -
Mplus中,预测变量中的类别变量不是用
NOMINAL
语句声明,否则在后续建立模型时会报错这不是因变量。此处需要用DEFINE语句将类别变量转化为哑变量。 -
Mplus中,
MODEL
部分的定义与线性回归类似,用ON
表示预测变量与结果变量的关系。
以下是废话多说版。
Logistic回归简介
第一次听说Logistic回归,我的反应是“啊是不是跟线性回归差不多的啦,看看几个变量对另一个变量的预测作用,每个变量都有一个回归系数,随着自变量的增减,因变量也相应地增减嘛。”可以这么理解Logistic回归吗?对,也不对。
Logistic回归确实和线性回归有相似之处,都是建立一个回归模型,有预测变量,有相应的回归系数。但是,Logistic回归与线性回归的区别在于它们的因变量的类型。Logistic回归的因变量不是连续数据,而是类别数据——可能是二分变量,对应的就是二分类Logistic回归(Binary Logistic Regression);也可能不止两类,是包含了多个分类的变量,对应的就是多分类Logistic回归(Multinomial Logistic Regression)。
结果变量的logit转换
那么,如何建立起这样的回归模型,来描述预测变量(可以是类别变量,也可以是连续变量)对结果变量(类别变量)的影响呢?Logistic回归采用的办法是logit转换:
\[logit(P) = ln\frac{P}{1-P} \]如果结果变量是二分变量,\(P\)表示的是其中一种类别出现的概率,而\(1-P\)则表示另一种类别出现的概率,\(\frac{P}{1 - P}\)被称作事件的发生比(odds)。例如,某种疾病的发病率还是未发病率;又例如,人为将某一个指标得分划分为高、低两个水平,\(P\)和\(1-P\)就分别表示属于高水平和低水平的概率。
那如果结果变量是多分类变量呢?几个分类两两之间就不是对立事件,不能直接用\(P\)和\(1-P\)去表示它们的概率。于是,多分类Logistic回归的logit转换就写成了:\(ln\frac{P(B)}{P(A)}\), \(ln\frac{P(C)}{P(A)}\), \(ln\frac{P(D)}{P(A)} \dots\)其中,\(P(A)\)、\(P(B)\)、\(P(C)\)、\(P(D)\)分别代表归属于类别A、B、C、D的概率(或者说,事件A、B、C、D发生的概率)。在这里我们可以看到\(P(A)\)总是出现在分母的位置,为什么呢?因为多分类Logistic回归模型中通常需要指定一个参照组别(reference group),拿其他类别与这个参照组一一作比较,在上面的例子中就是以类别A为参照,用类别B、C、D依次与类别A对比。
经过logit转换后,\(logit(P)\)就拥有了连续的范围。这时候,我们可以参照多元线性回归模型:
\[y = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + \dots +\beta_{k}x_{k} + \mu \]建立一个Logistic回归模型:
\[logit(P) = \alpha + \beta_{1}x_{1} + \beta_{2}x_{2} + \dots + \beta_{k}x_{k} \]可以看出,Logistic回归是假设了各个预测变量与logit转换后的结果变量存在线性关系。
预测变量中类别变量的处理
了解了Logistic回归中对结果变量的处理方式,接下来让我们将目光转向Logistic回归模型中等号的右边,预测变量。在前面已经说过,预测变量可以是连续变量,也可以是类别变量,如果是连续变量,用\(\beta x\)来表示回归系数和预测变量就类似于多元线性回归,很好理解,变量\(x\)的值增加一个单位,结果变量的值(这里是\(ln\frac{P}{1-P}\)改变\(\beta\)个单位。那如果是类别变量呢?类别变量的值显然不会这样连续变化。因此,在Logistic回归中,对于类别变量,需要将其转化为哑变量(dummy variable,也叫虚拟变量)。
具体而言,如果一个类别变量\(X\)共有\(N\)种类别,那么需要建立\(N - 1\)个哑变量\(A_{1}, A_{2}, A_{3}, \dots, A_{N - 1}\),根据\(X\)的值对这\(N\)个哑变量赋值,对于前\(N - 1\)个类别,只有对应的那个哑变量赋值为\(1\),其他哑变量都赋值为\(0\);而当\(X\)值为第\(N\)个类别时,所有\(N - 1\)个哑变量都赋值为\(0\),以此表示\(X\)的所有\(n\)种可能取值:
\[X = 1: A_{1} = 1, A_{2} = A_{3} = \dots = A_{N - 1} = 0 \]\[X = k, k \neq N: A_{k} = 1, A_{1} = A_{2} = \dots = A_{k - 1} = A_{k + 1} = \dots = A_{N - 1} = 0 \]\[X = N: A_{1} = A_{2} = A_{3} = \dots = A_{N - 1} = 0 \]这样转化后,原先的Logistic回归模型\(logit(P) = \alpha + \beta x\)就转化为了如下形式,每个哑变量都有一个单独对应的回归系数。
\[logit(P) = \alpha + \beta_{A_{1}} A_{1} + \beta_{A_{2}} A_{2} + \dots + \beta_{A_{N - 1}} A_{N - 1} \]模型结果的判读
从上面的Logistic回归模型可以看出,每个变量对应的回归系数\(\beta\)表示的是这个变量增加一个单位,等号左边的\(logit(P)\)的值随之增加\(\beta\)个单位(\(\beta\)为正,\(logit(P)\)增加;\(\beta\)为负,\(logit(P)\)降低)。但是,这样的结果仅说明了预测变量与\(logit(P)\)的关系,然而我们分析需要了解的是预测变量与结果变量的关系,这样的结果还不够直观。因此,通常需要将等号两边同时作为\(e\)的指数,化为如下形式。(在二分类Logistic回归中,下式的\(\frac{P(B)}{P(A)}\)即为\(\frac{P}{1 - P}\)。)
\[\frac{P(B)}{P(A)} = e^{\alpha} \cdot e^{\beta_{1}x_{1}} \cdot e^{\beta_{2}x_{2}} \cdot \dots \cdot e^{\beta_{k}x_{k}} \]由此可以看到,变量\(X_{k}\)增加一个单位,\(\frac{P(B)}{P(A)}\)变为原来的\(e^{\beta_{k}}\)倍。\(e^{\beta}\)即为Logistic回归分析中需要报告的\(OR\)值(odds ratio),称为“比值比”或者“优势比”。显然,\(OR\)比\(\beta\)更能反映预测变量对结果变量的作用。
预测变量为连续变量
预测变量为连续变量时,\(OR\)大于1,说明预测变量\(X_{k}\)增大,\(\frac{P}{1 - P}\)增大,以\(P\)表示发病率为例,即为“发病率”:“不发病率”的比值增大,也就是\(X_{k}\)增大,发病的概率增大。在多分类Logistic回归中,则是\(X_{k}\)增大,\(\frac{P(B)}{P(A)}\)增大,说明随着\(X_{k}\)增大,事件B比事件A更有可能发生。
同理,\(OR\)小于1,说明预测变量\(X_{k}\)增大,\(\frac{P}{1 - P}\)或\(\frac{P(B)}{P(A)}\)减小。
预测变量为类别变量
预测变量为类别变量时,对\(OR\)的解读会比连续变量更复杂一点点。还记得在之前定义的哑变量吗?对于有\(N\)种可能类别的变量\(X\),我们建立了\(N - 1\)个哑变量\(A_{1}, A_{2}, \dots, A_{N - 1}\),每个哑变量都有各自的\(\beta\)和\(OR\)值,于是这个变量\(X\)也就有了\(N - 1\)个\(OR\)值。要如何解释这些\(OR\)值的现实含义呢?
首先我们把具体的变量\(X\)和哑变量的取值代入Logistic回归模型来看,当\(X = k, k \neq N\)时,由于只有\(A_{k}\)值为\(1\),其他哑变量取值皆为\(0\),因此有:
\[logit(P) = \alpha + \beta_{k} \]此处\(logit(P)\)表示的是,在\(X = k\)的条件下,\(ln\frac{P}{1 - P}\)的值,我们把它记作\(ln\frac{P_{k}}{1 - P_{k}}\)好了。在多分类Logistic回归中,则是在\(X = k\)的条件下,\(ln\frac{P(B)}{P(A)}\)的值,也就是\(ln\frac{P(B | X = k)}{P(A | X = k)}\)。
那当\(X = N\)时呢?所有哑变量取值皆为\(0\),有:
\[logit(P) = \alpha \]同理,我们把\(logit(P)\)写成\(ln\frac{P_{N}}{1 - P_{N}}\),在多分类Logistic回归中可以写成\(ln\frac{P(B | X = N)}{P(A | X = N)}\)。
两式相减,得到了:
\[\beta_{k} = ln\Big\{\frac{P_{k}}{1 - P_{k}} \Big/ \frac{P_{N}}{1 - P_{N}}\Big\} \]取\(e\)的指数,化为\(OR\)值的形式,得到:
\[OR_{k} = \frac{P_{k}}{1 - P_{k}} \Big/ \frac{P_{N}}{1 - P_{N}} \]在多分类Logistic回归中,\(\beta\)和\(OR\)则为:
\[\beta_{k} = ln\Big\{\frac{P(B | X = k)}{P(A | X = k)} \Big/ \frac{P(B | X = N)}{P(A | X = N)}\Big\}, OR_{k} = \frac{P(B | X = k)}{P(A | X = k)} \Big/ \frac{P(B | X = N)}{P(A | X = N)} \]可见,当预测变量为类别变量时,\(OR\)表示的是类别k和类别N的发生比之比。\(OR\)大于1,说明在类别k的条件下,事件的发生比高于类别N条件下的发生比。例如,用变量\(X\)表示性别,\(X = 1\)表示男性,\(X = 0\)表示女性,\(P\)表示发病率,\(1 - P\)表示未发病率,\(OR_{1} = \frac{P_{1}}{1 - P_{1}} \Big/ \frac{P_{0}}{1 - P_{0}}\)表示的是男性和女性发病:不发病的比数比,\(OR\)大于1,说明男性比女性更有可能发病。
这里还有一点需要注意,所有哑变量都赋值为\(0\)的类别N其实就成为了其他类别的参照,所有的\(OR\)值解释都是基于“某某类别相对于类别N”的比值得到的,因此类别N内包含的样本量也有一点讲究,如果样本量太小,可能导致\(OR\)值有较大的偏差。其实不止作为参照的类别N,所有类别都可能存在样本量太小导致\(OR\)不可靠这样的问题,因此在定义哑变量之前可以考虑先把样本量较小的几个类别合并一下,以减小类别数量、提高类别内的样本量。
OR的显著性
对Logistic回归系数的显著性判断主要有两种依据,一是\(p\)值,二是置信区间。在Mplus中输出的结果给出了\(\beta\)(在MODEL RESULTS部分)和\(OR\)(在LOGISTIC REGRESSION ODDS RATIO RESULTS部分)的估计值、标准误、\(p\)值,Mplus对\(p\)值的计算是根据估计值和标准误的比值\(\frac{Est.}{S.E.}\)(计算\(OR\)的\(p\)值用的是\(\frac{Est. - 1}{S.E.}\))近似正态分布。在R语言中使用glm()
函数也是用这种方法计算\(p\)值。而在SPSS中,则是用Wald统计量(\(W = \Big(\frac{\beta}{S.E.}\Big)^{2}\),符合卡方分布)计算\(\beta\)的\(p\)值。
结合\(OR\)的置信区间判断显著性主要是看置信区间整体是否在\(1\)的一侧,一般显著的\(OR\)值对应的置信区间都会出现在\(1\)的同侧,比如\([0.7, 0.9]\)或者\([1.23, 1.56]\)(数字我随便打的哈)。而不显著的置信区间是什么样子呢?左端点小于\(1\)而右端点大于\(1\),就像是一个人“骑墙”,例如\([0.8, 1.2]\),从这样的置信区间看不出随着预测变量变化,发生比到底是增大还是减小。要注意,Mplus中默认不输出置信区间,需要在OUTPUT
中添加一句CINTERVAL
以获得置信区间。
Mplus代码实现
使用Mplus只需要几行代码便能够简单便捷地完成Logistic回归分析,具体操作如下。
在inp文件开头是DATA
、VARIABLE
、USEVARIABLES
分别声明数据文件、文件中每列对应的变量名、分析需要用到的变量,都是Mplus的常规操作,与其他分析一样,此处略去不表。可以用MISSING
声明哪些数据被标为缺失值,在Mplus的Logistic回归分析中,缺失值都会被直接以列删法(listwise)排除。
首先要注意的是,由于Logistic回归的结果变量为类别变量,需要用NOMINAL
或CATEGORICAL
语句加以说明,如果是有序类别变量或者二分类别变量,用CATEGORICAL
;如果是无序类别变量,用NOMINAL
。否则,Mplus会将结果变量视作连续变量处理,进行线性回归分析而非Logistic回归分析。
随后,如果预测变量中有类别变量,需要我们自己用DEFINE
语句定义哑变量,格式为哑变量 = 原有变量 == 类别编号
,也就是说,当原有变量值为指定类别时(“原有变量 == 类别编号”判断结果为True),这个哑变量赋值为1,否则为0。记得将这里新定义的哑变量也加入前面的USEVARIABLES
里面。
接下来是指定分析类型,ANALYSIS
语句下需要声明的是ESTIMATOR = ML
,指定采用极大似然估计值。加上这一句,Mplus会进行Logistic回归分析,如果不加这一句则是进行Probit回归分析。
最后是模型定义,与线性回归一样,用ON
命令来表示回归模型,结果变量放在前面,一系列预测变量(连续变量和哑变量)放在ON
后面。没有特殊说明的情况下,Mplus会默认将结果变量的最后一个类别作为参照组,但是在实际应用中,我们可能会选取不同类别作为参照组,多做几次Logistic回归分析,于是我们可以用#
指定参与分析的类别,例如,结果变量\(u1\)有编号为1、2、3、4的4个类别,首先我们选取第4类作为参照组,那么可以写成u1 ON 某某预测变量
让Mplus默认定义最后一类第4类为参照,也可以写成u1#1 u1#2 u1#3 ON 某某预测变量
指定计算前3类(相对第4类)的回归系数。那假如是选取第3类作为参照呢?同理,用u1#1 u1#2 ON 某某预测变量
即可,注意,这里不用把u1#4
放进模型,会报错,因为第3类和第4类之间的对比在之前设定第4类为参照就可以实现了。
到这里,Mplus就已经可以顺利跑起来,在MODEL RESULTS部分和LOGISTIC REGRESSION ODDS RATIO RESULTS部分输出\(\beta\)和\(OR\)的估计值、\(p\)值了。但是此时是不输出置信区间的,如果想要再看看\(\beta\)和\(OR\)的置信区间,可以再加上OUTPUT: CINTERVAL
语句,就能得到99%、95%和90%置信区间了。
下面是具体代码。这个例子是为了分析青少年的年龄、性别、家庭收入、父母教育程度对其心理行为问题的预测作用,变量含义如下:
- index:样本编号。
- group:心理问题类型,结果变量,无序多分类变量,包含了4类心理症状,编号为1、2、3、4。
- age:年龄,预测变量,作为连续变量处理。
- gender:性别,预测变量,1为男性,0为女性。因为本身只有2类,不用再定义哑变量,最后得到的\(OR\)即为男性:女性的比值比。
- income:家庭收入,预测变量,作为连续变量处理。
- edu:教育程度,预测变量,1为初中及以下,2为高中大专,3为本科及以上。
TITLE: Multinomial Logistic Regression
! DATA 说明你的数据文件是哪一个。
DATA: FILE IS BlahBlahBlah.dat;
! VARIABLE 依次列出数据文件中的所有变量名。
VARIABLE: NAMES ARE index group age gender income edu
! USEVARIABLES 列出这次分析中要用到的变量名,包括原先数据文件中的变量,也包括后面自己新定义的变量。
USEVARIABLES ARE group age gender income middle high;
! MISSING 说明数据缺失值都已被赋值为1000。
MISSING = ALL(1000);
! 结果变量要用NOMINAL语句说明,否则会按线性回归分析。
NOMINAL = group;
! 此处定义2个哑变量:1初中及以下,2高中大专。第3类“本科及以上”作为教育程度的参照。
DEFINE:
middle = edu == 1;
high = edu == 2;
ANALYSIS:
! ML, logistic regression
ESTIMATOR = ML;
MODEL:
! 用ON语句定义回归模型。
group ON age gender income middle high;
OUTPUT:
! 输出置信区间。
CINTERVAL;
参考文献
刘红云. (2019). 高级心理统计. 中国人民大学出版社.
汪海波, 罗莉, & 汪海玲. (2018). R语言统计分析与应用. 人民邮电出版社.
王孟成. (2014). 潜变量建模与Mplus应用 基础篇. 重庆大学出版社.
Muthén, L.K., & Muthén, B.O. (2017). Mplus User's Guide (8th ed.). Los Angeles, CA: Muthén & Muthén.
UCLA: Statistical Consulting Group. Multinomial logistic regression Mplus data Analysis. https://stats.oarc.ucla.edu/mplus/dae/multinomiallogistic-regression/
最后闲扯几句
之所以想要写这一篇,是因为之前写论文分析数据需要做Logistic回归,然而在中文网上找了很久都没找到Mplus如何处理预测变量中的类别变量,(可能大家都用SPSS省事?)最后是在UCLA: Statistical Consulting Group的网页中找到了方法。因此想在这里发一篇中文版的笔记,或许未来有人也有同样的疑问呢?
这一篇拖了好久了,写开头的时候我还在分析数据,现在写结尾的时候我论文都写完提交盲审了。早知道写一篇博文这么累我就不写了hhhhhh。不过,用自己的语言把学到的内容重新讲述一遍,这也是梳理知识点的过程吧。
关于Logistic回归的内容还有很多很多,像回归分析中的同时或逐步进入方法、模型拟合评价涉及到的似然值(likelihood value)、负2倍对数似然值(\(-2LL\))、伪测定系数(\(pseudo-R^{2}\))、\(AIC\)等等都没有在这篇文中介绍。日后我有时间琢磨琢磨再来补充吧。
这篇文主要是基于我自学之后对Logistic回归的理解,难免存在纰漏谬误或是含糊不清之处,欢迎批评讨论!
标签:frac,变量,回归,笔记,beta,Logistic,类别,Mplus From: https://www.cnblogs.com/GhouGHpsy/p/17404334.html