R语言分层线性模型案例

有许多分层数据的例子。例如，地理数据通常按层次分组，可能是顶级的全球数据，然后按国家和地区分组。一个生物学的例子是按物种分组的动物或植物的属性，或者属于一个级别的属性，然后是家族（例如，参见蓝山雀的Jarrod Hadfield的数据集）。一个金融的例子可能是按行业分组的特定国家的公司（例如，见Pak-Wing Fok等人使用分层多因素模型分析信贷组合风险的论文））。天体物理学的例子可能是由星系类型（椭圆形，螺旋形）形成恒星的速率。一个商业例子可能是业务部门和细分的员工满意度（参见Whitener Do的文章）。每个学科都有许多例子，其中观察以某种形式的层次结构进行分组。

在这里，我想解释使用一个简单的例子，如何使用R来构建分层线性模型。我在整个三组中使用简单的一维数据集。在每个组内，自变量x和因变量y之间存在强正关系。

library(ggplot2)

df$group <- as.factor(df$group)

tail(df)

g <- ggplot(df) + geom_point(aes(x=x,y=y,colour=group)) +

coord_cartesian(xlim=c(-3,1),ylim=c(-4,4)) +

geom_smooth(aes(x=x,y=y,group=group),method=lm,se=FALSE) +

theme_bw() + theme(legend.position="null")

g + geom_smooth(aes(x=x,y=y),method=lm,se=TRUE)

这就是数据的样子。这些组有不同的颜色。在本文的其余部分，我将展示如何使用层次模型来模拟这种情况，该模型确实考虑了组信息。

Gelman建议的分层线性模型的一个包是arm，它具有与lm（）函数非常相似的函数lmer（）。

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

lmer.both <- lmer(y~1+x+(1+x|group),data=df) summary(lmer.alpha) summary(lmer.beta) summary(lmer.both) # fixed effect is top-level intercept and slope # (Intercept) x # 1.978652 1.144952 # group random effects on intercept # > ranef(lmer.alpha)

# $group

# (Intercept)

# 1 3.4386106

# 2 -0.8360106

# 3 -2.6026000

# > group.alpha

# [1] 4.2883814 1.2134493 -0.5410049

# > ranef(lmer.alpha)$group[,1] + fixef(lmer.alpha)[1]

# [1] 5.4172624 1.1426413 -0.6239482

fixef(lmer.alpha)

ranef(lmer.alpha)

group.alpha

# fixed effect is top-level intercept

# (Intercept)

# 5.788223

# group random effects on intercept and slope

# (Intercept) x

# 1 -1.740225 0.518047

# 2 -4.564296 1.415710

# 3 -6.354477 1.231584

# > group.alpha

# [1] 4.2883814 1.2134493 -0.5410049

# > ranef(lmer.beta)$group + fixef(lmer.beta)[2]

# [1] 4.0479981 1.2239268 -0.5662542

fixef(lmer.beta)

ranef(lmer.beta)

group.beta

# > fixef(lmer.both)

# (Intercept) x

# 1.578741 1.059370

# > ranef(lmer.both)

# $group

# (Intercept) x

# 1 2.500014 -0.5272426

# 2 -0.355365 0.3545068

# 3 -2.144649 0.1727358

fixef(lmer.both)

ranef(lmer.both)

# if we have no pooling we would simply run 3 regressions, one per group

coef(lm(y~x,data=df[group==1,]))

coef(lm(y~x,data=df[group==2,]))

coef(lm(y~x,data=df[group==3,]))

# (Intercept) x

# 4.0653645 0.5259707

# 1.227969 1.428500

# -0.570280 1.225905

# true values for group.alpha are

# 4.2883814 1.2134493 -0.5410049

(ranef(lmer.alpha)$group[,1]) + fixef(lmer.alpha)[1]

(ranef(lmer.beta)$group[,1]) + fixef(lmer.beta)[1]

# alpha model plot

fit.lines <- data.frame(cbind(intercept=(ranef(lmer.alpha)$group[,1])+fixef(lmer.alpha)[[1]]

,slope=fixef(lmer.alpha)[[2]])

,group=factor(1:3))

g.alpha <- g + geom_abline(aes(intercept=intercept,slope=slope,hierarchicalcolour=group),data=fit.lines) +

ggtitle("Alpha Random Effects")

g.alpha

# beta model plot

fit.lines$group <- factor(fit.lines$group)

g.beta <- g + geom_abline(aes(intercept=intercept,slope=slope,colour=group),data=fit.lines) +

ggtitle("Beta Random Effects")

# get this handy function from http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/

multiplot(g.alpha,g.beta,g.both,cols = 3)

结果显示在这里。我有三个图，第一个是截距（alpha）依赖于组，第二个是斜率（β）依赖于组，第三个是截距和斜率都取决于组。你可能在想为什么不只是做三个单独的线性回归，因为第三个例子产生的系数非常接近于此。原因是基于这样的假设：alphas和beta是从顶级分布中提取的，因此是相关的。这意味着我们可以在组之间汇集信息，如果我们为其中一个组提供的数据非常少。

术语是回归系数是“固定效应”，组级别称为“随机效应”。

lmer.both.sim <- sim(lmer.both)

fixef.lmer.both <- fixef(lmer.both.sim)

nsamples <- 100

fit.lines.both$group <- factor(rep(1:3,each=nsamples))

# now do 3 separate linear regressions (one for each group)

lm.mcmc.1 <- MCMCglmm(y~1+x,data=df2[df2$group=="1",],thin = 100)

fit.lines.mcmc <- data.frame(rbind(lm.mcmc.1$Sol,lm.mcmc.2$Sol,lm.mcmc.3$Sol))

g.sim.mcmc <- g2 + geom_abline(aes(intercept=intercept,slope=slope,colour=group),data=fit.lines.mcmc,alpha=.2) +

结果如下所示。每组只有一个单独的线性回归。对于蓝色和红色组，线条在大多数情况下非常适合数据，但对于只有三个数据点的绿色组，线条遍布整个地方，因为没有任何先验信息，估计数据的斜率和偏移量非常不确定。右侧的图表显示因为该模型假设所有三组的斜率和偏移都是从一个公共分布中得出的，所以可以合理地假设斜率是正的。我们知道这适用于这个例子，因为我们设计了数据生成过程。

非常感谢您阅读本文，有任何问题请在下面留言！

非常感谢您阅读本文，有任何问题请在下面留言！

相关文章