关于GLMM（generalized linear mixed model）广义线性混合模型

在研究问题时，当变量属于正态分布的随机变量，服从样本独立假定时，我们会选择LM-线性模型；当变量之间不独立时，此时引入随机效应，考虑响应变量之间的相关性，选择LMM-线性混合模型；当响应变量不服从正态分布，但是样本之间独立时，使用GLM-广义线性模型；响应变量不服从正态分布，样本之间有相关性时，使用GLMM-广义线性混合模型。

1、广义线性混合模型

GLMM（generalized linear mixed model）广义线性混合模型中的关键是“mixed”，“mixed”是区别于一般的GLM（generalized linear model）的显著体现。

一般的GLM指的就是要求因变量符合“指数分布族”即可。关于GLM的详细解释可以在stata的help文档中看到，GLM的两个核心是 Family 和 Link。其中Family指的就是因变量的分布函数，常见的几种因变量的分布如下：
连续变量——Gaussian分布/正态分布
binary变量（0，1）——二项分布（一个变量0，1分布是属于伯努利分布，多个伯努利就是二项分布）
count/rate变量（1，2，3…计数分布）——possion分布

其中link指的是将因变量进行转化使模型呈线性分布的方法，我们经常使用的logit回归的link就是logit link，它的方法是对因变量进行如下变换：ln(P(Y=1)/P(Y=0))=aX+b。除此之外，还有其他多种不同的回归就是不同的Family和link的组合，logistic回归只是GLM回归的一种特殊情况，这就是为什么GLM被称为广义的线性模型，是所有线性模型的最初模样。

GLMM模型，它的特殊之处在于它的“mixed”，这个”mixed"本质上说的是效应的混合，因为回归模型中同时有固定效应和随机效应的存在，固定效应说的是自变量每次变化一个单位对因变量的影响是一样的，但是随机效应说的是自变量变化一个单位对因变量的影响是不固定的。而这种情况时有发生，如果在一个模型里面，多个不同的自变量对因变量的影响同时存在这两种效应，我们使用到的模型就是所谓的GLMM模型。GLMM引入了连接函数来解决变量非正态分布的问题，同时可以引入随机效应来解决数据间的相关、过度离散和异质性问题。

2、GLMM的基本程序实现

2.1 Stata基本实现

Stata提供了可视化操作界面，这有助于我们找到自己想要的模型。操作如下：
“ 导入数据 > 统计 > 多层混合效应模型 > 广义线性模型（GLM）”
以工资为例，这里显示的是工资高低，工资高为1，工资低记为0，所以回归在这里是一个logit回归，由上面的解释，可以知道分布为二项分布，因此family和link设置如下。
GLMM回归的基本常规语法：，可以从菜单界面选择，也可以自己输入下列语法到工作台输入栏

meglm salary hours education gender age region || year:, family(binomial) link(logit)

当工资表示为连续变量时，根据工资的实际分布情况，我们这里暂时认为工资是正态分布的，也就是高斯分布，此时可以使用以下回归方程：

meglm salary hours education gender age region || year:, family(gaussian) link(log)

2.2 R基本实现方法

R的绘图效果没有Stata好，实现的结果和Stata基本一致。
给出R中实现GLMM回归的代码如下：
R中关于广义线性混合模型的主要包库是lme4和lmerTest

setwd("D://research")
library(readxl)
library(lme4)
library(car)
library(lmerTest)
data1 = read_excel("Data.xlsx")
data1$gender <- factor(data1$gender)    #将变量设置为分类变量/虚拟变量
fit1=glmer(salary~hour+education+gender+age+region(1|year),
data=data1,family=binomial("logit"),nAGQ=1)    #定义模型方程
summary(fit1)    #输出模型结果

3、GLMM的交叉项分析

3.1 Stata实现

很多时候，在文章中我们需要考虑交叉项，来体现文章的思考深度。
在上面的例子中，hour和education都是连续变量，希望在回归中体现他们的交叉项，因此通过c.hour##c.education这样设计，或者是生成交乘变量再加入模型，后面还考虑了可能影响工资高低的性别，年龄，区域和年份因素，年份作为随机条件。（这里假设工资随年份变化不大）。i.gender是将性别设为虚拟变量。

meglm salary c.hour##c.education hour education i.gender age region || year:, family(binomial) link(logit)

得到结果之后，反应的主要结论是：不同工作时间的人的工资高低。此外，交互项我们希望得出不同教育水平的工人，工作时间对他们工资的影响高低。因此，我们希望得出边际效应。这里的education=(1(1)4) 表示我希望结果输出时，横坐标开始是1，步长为1，终点为1。

margins, eydx(hour) at(education=(1(1)4))

这里需要注意的是，如果是logit模型，需要对边际效应进行一个转换。使用eydx的指令，如果是一般回归，我们使用dydx。
参考文章：二值模型的Stata命令
最后输出带有置信区间的，边际效应折线图，类似下图

marginsplot

3.2 R实现

在R中也可以实现交互项的输出，在R中进行交互项的输出，首先是建立在已经输出模型结果的前提之下，具体操作的语句如下：
这个例子中，institute是一个定序变量，取值为1，2，3

summary(margins(fit1, variables = "inter", at = list(institute = 1:3),type="link"))

4、置信区间的计算方法

4.1 回归系数的置信区间

Stata可以直接给出回归结果中各个自变量系数的置信区间；
R中很多时候不能显示出置信区间，但是R的输出结果中会给出标准差的值，在z检验的情况下，这个值通常是1.96*SE得到的，通过系数的值加减这个得出来的结果，我们可以得到最终的置信区间。
至于为什么要用1.96*SE，解释应该是当样本容量很大的时候，置信水平95%时，α=0.05，Z(0.025)=1.96。

4.2 变量/统计量的置信区间

样本终究是样本，是随机抽取的结果，不能完全代表整体结果，所以有的时候会继续在我们选择的样本中进行Bootstrap抽样来得出对统计量的估计，Bootstrap抽样-是一种又放回的均匀抽样。这种抽样方法，加上我们的置信区间，就可以得到对样本统计量估计的准确范围。

均值的置信区间-Stata

bootstrap: mean salary

变量的置信区间-Stata

bootstrap: mean No_of_papers if age==1&inter==2    #对其他变量进行限制

5、参考文献

边际效应分析参考-区分连续性和虚拟变量