在研究问题时,当变量属于正态分布的随机变量,服从样本独立假定时,我们会选择LM-线性模型;当变量之间不独立时,此时引入随机效应,考虑响应变量之间的相关性,选择LMM-线性混合模型;当响应变量不服从正态分布,但是样本之间独立时,使用GLM-广义线性模型;响应变量不服从正态分布,样本之间有相关性时,使用GLMM-广义线性混合模型。
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引入了连接函数来解决变量非正态分布的问题,同时可以引入随机效应来解决数据间的相关、过度离散和异质性问题。
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)
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)
很多时候,在文章中我们需要考虑交叉项,来体现文章的思考深度。
在上面的例子中,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
在R中也可以实现交互项的输出,在R中进行交互项的输出,首先是建立在已经输出模型结果的前提之下,具体操作的语句如下:
这个例子中,institute是一个定序变量,取值为1,2,3
summary(margins(fit1, variables = "inter", at = list(institute = 1:3),type="link"))
Stata可以直接给出回归结果中各个自变量系数的置信区间;
R中很多时候不能显示出置信区间,但是R的输出结果中会给出标准差的值,在z检验的情况下,这个值通常是1.96*SE得到的,通过系数的值加减这个得出来的结果,我们可以得到最终的置信区间。
至于为什么要用1.96*SE,解释应该是当样本容量很大的时候,置信水平95%时,α=0.05,Z(0.025)=1.96。
样本终究是样本,是随机抽取的结果,不能完全代表整体结果,所以有的时候会继续在我们选择的样本中进行Bootstrap抽样来得出对统计量的估计,Bootstrap抽样-是一种又放回的均匀抽样。这种抽样方法,加上我们的置信区间,就可以得到对样本统计量估计的准确范围。
均值的置信区间-Stata
bootstrap: mean salary
变量的置信区间-Stata
bootstrap: mean No_of_papers if age==1&inter==2
边际效应分析参考-区分连续性和虚拟变量
但在中午的时候,有的人还没吃饭, 有的人吃过饭了,有的人喝了酒,结果酒精和药物起了反应,有的人喝了醋,醋又和药物起了另一种反应。而对于随机效应而言,参数是服从正态分布的一个随机变量,也就是说对于两个不同的自变量的值,对应变量的影响不一定是相同的。rnd 一个描述模型随机效应部分的双边线性公式对象,在算符左侧有分组因子,右侧有 + 算符分隔的随机项,直接给出随机效应设计矩阵(列名合适)。GLM的缺点 GLM模型可以支持数据的多种分布,但是要求数据是独立不相关的,在遗传评估时,很显然数据不能满足要求。
GLM 一般是指 generalized linear model ,也就是广义线性模型;而非 general linear model,也就是一般线性模型;而GLMM (generalized linear mixed model)是广义线性混合模型。
广义线性模型GLM很简单,举个例子,药物的疗效和服用药物的剂量有关。这个相关性可能是多种多样的,可能是简单线性关系(发烧时吃一片...
广义线性混合模型GLMM(Generalized Linear Mixed Model),是广义线性模型GLM 和线性混淆模型LMM 的扩展形式,于二十世纪九十年代被提出。GLMM因其借鉴了混合模型的思想,其在处理纵向数据(重复测量资料)时,被认为具有独特的优势。GLMM不仅擅长处理重复测量资料,还可以用于任何层次结构的数据(因为本质上又是多水平模型)。
提到GLMM,有必要先介绍几个容易混淆的...
全文链接:http://tecdat.cn/?p=30914我们正和一位朋友讨论如何在R软件中用GLM模型处理全国的气候数据。本文获取了全国的2021年全国的气候数据(点击文末“阅读原文”获取完整代码数据)。采样时间:2021年1月1号~2021年12月31号采样地点:全国各地。本次调查搜集了2021年全国不同地区的风向、降雨量、风速、风速变化、最大风速、最大降雨量、闪电概率等数据。并对不同变量...
`glmmTMB`是一个建立在 Template Model Builder(https://github.com/kaskr/adcomp)自动分化引擎上的R包,用于拟合广义线性混合模型及其延伸。
三部分:固定项age,随机项(1|subject)和误差项ε。
为什么要加上一个随机项这部分呢?
在线性模型中我们将所有的不感兴趣的因素,非系统性的因素,不可预测的因素造成的误差统统由一个ε来代替。这样我们求出的...
在上面的示例中,我们首先创建一个包含固定效应和随机效应的数据框data,其中y表示因变量,x表示自变量,group表示分组变量。接着,我们使用glmer()函数拟合一个包含随机截距和斜率的混合效应模型,其中(1 + x | group)表示在分组变量group内拟合截距和斜率随自变量x变化的模型。需要注意的是,glmer()函数的模型公式与glm()函数的模型公式类似,但是在处理具有随机效应的数据时,需要考虑数据的层次结构和相关的随机效应。
glm price weight length, family(gaussian) link(identity)
2.Logistic回归
glm low age i.smoke i.race, family(binomial) link(logit)
glm low age i.smoke i.race, family(binomial) link(logit) eform
3.泊松回归
glm low age i.smoke i.race, family(poisson) link(log
0. 飞哥感言
这篇文章,主要是介绍了抗性数据,如何利用GLMM模型进行的分析,文中,他将9级分类性状变为了二分类性状,进行分析。
分析中用到了加性效应(A矩阵),空间分析(行列信息)。
对比了SAS和ASReml,结果基本一致。
其实,9分级性状,可以直接使用ASReml进行有序多分类性状分析,用累计Logistic模型分析,也可以考虑系谱数据和空间位置信息。这样效果应该更好。
回头找下数据,测试一下。
1. 文献
Genetic analysis of resistance to Pseudomona
重复测量数据有几个明显的特征,一是个体内数据是反复收集的,同一对象的多次观测结果往往不独立(存在相关性),二是变异来源上看有个体内变异和个体间变异,三是数据可能存在缺失值。有多个统计模型可以实现重复测量数据的分析:【1】一般线性模型中的重复测量方差分析,可以采用一元方差分析和多元方差分析。重复测量方差分析要求还是比较苛刻的,要求多元正态性、组间方差-协方差矩阵相等(Box’M检验),数据...
在第一个选项卡上,该函数显示用户选择的数据的预测区间。该函数通过从固定效应和随机效应项的模拟分布中抽样并组合这些模拟估计来快速计算预测区间,以产生每个观察的预测分布。对于每种情况,最多12个,在所选数据类型中,用户可以查看更改固定效应的影响。这允许用户比较变量之间的效果大小,以及相同数据之间的模型之间的效果大小。在下一个选项卡上,固定效应和组级效果的分布在置信区间图上显示。解释LMM和GLMM模型的结果很困难,尤其是不同参数对预测结果的相对影响。最简单的是得到固定和随机效应参数的后验分布。
进行数据分析时,会发现有时候一个模型中的变量之间可能具有相关性(correlation),比如面积和长度就具有高度的相关性,如果同时对这些参数建模,就存在共线性问题,所以一般是只针对其中一个参数建模。而这种相关性,其实还存在于数据之中,比如时间序列数据,在不同的时间,同一个对象的数据之间就是相互有联系的,那么我们应该怎么对这些具有相关性的数据进行建模分析呢。
在进一步分析之前,再次强调一下,这里...
R语言中,构建GLMM 模型时,一个好的选择是使用“lme4” 包中的“lmer()”函数。
前些天一个小伙伴问我 :"下面这个模型中,(1|car_type) 是啥意思啊?"
lmer(wear~wheel+(1|car_type))
我其实一直知道这是一个 "固定部分(wheel) + 随机部分 (car_type)" , 而且它仅仅是 随机截距而 没有随机斜...
快一个月没更新文章啦,今天收到好几个粉丝的催更私信,好的吧,实在对不住大家期待的眼神,看样子不能再拖啦,想想写啥好呢,大家咨询比较多的,混合模型算一个,今天就继续给大家写写混合模型如何做吧。
混合模型一般都可以用lme4这个包解决,lme4既可以做线性混合模型,也可以做广义线性混合模型还可以做非线性混合模型,大家有需要可以只研究这一个包就行。
所谓混合模型就是既有固定效应又有随机效应的模型:
“mixedeffects”, denotes a model that incorporates both