线性混合模型R实现的更多实例

最新推荐文章于 2024-04-03 16:41:45 发布
最新推荐文章于 2024-04-03 16:41:45 发布
阅读量2.8k 收藏 21
点赞数 1
分类专栏： 统计分析 R语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_42458954/article/details/115910269
版权 
library(nlme)
## Warning: package 'nlme' was built under R version 3.6.3
library(lme4)
## Warning: package 'lme4' was built under R version 3.6.3
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 3.6.3
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(AED)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.6.3

使用示例数据集，比较不同协变量调整下线性回归模型的效能

data("RIKZ")

mod_lm1<-lm(Richness~NAP,data=RIKZ)
mod_lm2<-lm(Richness~NAP+Exposure,data=RIKZ)
mod_lm3<-lm(Richness~NAP+Exposure+Beach,data=RIKZ)

summary(mod_lm1)
## 
## Call:
## lm(formula = Richness ~ NAP, data = RIKZ)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0675 -2.7607 -0.8029  1.3534 13.8723 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.6857     0.6578  10.164 5.25e-13 ***
## NAP          -2.8669     0.6307  -4.545 4.42e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.16 on 43 degrees of freedom
## Multiple R-squared:  0.3245, Adjusted R-squared:  0.3088 
## F-statistic: 20.66 on 1 and 43 DF,  p-value: 4.418e-05
summary(mod_lm2)
## 
## Call:
## lm(formula = Richness ~ NAP + Exposure, data = RIKZ)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3083 -1.7107 -0.8489  0.7674 13.3264 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2909     5.1878   7.188 7.83e-09 ***
## NAP          -2.7252     0.4716  -5.779 8.26e-07 ***
## Exposure     -2.9988     0.5060  -5.926 5.07e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.106 on 42 degrees of freedom
## Multiple R-squared:  0.6321, Adjusted R-squared:  0.6146 
## F-statistic: 36.09 on 2 and 42 DF,  p-value: 7.577e-10
summary(mod_lm3)
## 
## Call:
## lm(formula = Richness ~ NAP + Exposure + Beach, data = RIKZ)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0393 -1.6114 -0.5544  0.7196 13.5599 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.3706     5.1212   7.102 1.18e-08 ***
## NAP          -2.5119     0.4810  -5.223 5.46e-06 ***
## Exposure     -2.7647     0.5170  -5.348 3.64e-06 ***
## Beach        -0.3095     0.1906  -1.623    0.112    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.048 on 41 degrees of freedom
## Multiple R-squared:  0.6543, Adjusted R-squared:  0.6291 
## F-statistic: 25.87 on 3 and 41 DF,  p-value: 1.484e-09

同时调整NAP,Exposure及Beach时，Beach与Richness的关联无统计学意义

绘制model1的诊断性曲线图

par(mfrow=c(2,2))
plot(mod_lm1)

model1的诊断性曲线提示model1中自变量不符合正态性及方差齐性，比较各个model的AIC值

AIC(mod_lm1,mod_lm2,mod_lm3)
##         df      AIC
## mod_lm1  3 259.9535
## mod_lm2  4 234.6083
## mod_lm3  5 233.8053

mod_lm3的AIC值最小，其中Beach无统计学意义，按Beach进行分组，绘制回归曲线

ggplot(RIKZ,aes(NAP,Richness,group=Beach,color=Beach))+
  geom_point(size=3)+
  geom_smooth(formula = "y~x",method = "lm",se=F,size=0.5)

发现不同Beach分组的回归曲线斜率与截距不一致，提示Beach仍然是Richness的影响因素,单独调整Beach时，结果如下

mod_lm4<-lm(Richness~NAP+Beach,data = RIKZ)
summary(mod_lm4)
## 
## Call:
## lm(formula = Richness ~ NAP + Beach, data = RIKZ)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8980 -2.4944 -0.4541  1.2401 14.2386 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.5053     1.2793   7.430 3.55e-09 ***
## NAP          -2.4363     0.6188  -3.937 0.000305 ***
## Beach        -0.5939     0.2357  -2.520 0.015622 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.923 on 42 degrees of freedom
## Multiple R-squared:  0.4133, Adjusted R-squared:  0.3853 
## F-statistic: 14.79 on 2 and 42 DF,  p-value: 1.372e-05

Beach仍有统计学意义，提示Exposure与Beach可能存在交互作用

par(mfrow=c(2,2))
plot(mod_lm4)

使用lme函数，添加Beach的随机截距项

mod_lmm1<-lme(Richness ~NAP,random = ~1|Beach,data=RIKZ)
summary(mod_lmm1)
## Linear mixed-effects model fit by REML
##  Data: RIKZ 
##        AIC     BIC    logLik
##   247.4802 254.525 -119.7401
## 
## Random effects:
##  Formula: ~1 | Beach
##         (Intercept) Residual
## StdDev:    2.944065  3.05977
## 
## Fixed effects: Richness ~ NAP 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  6.581893 1.0957618 35  6.006682       0
## NAP         -2.568400 0.4947246 35 -5.191574       0
##  Correlation: 
##     (Intr)
## NAP -0.157
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.4227495 -0.4848006 -0.1576462  0.2518966  3.9793918 
## 
## Number of Observations: 45
## Number of Groups: 9

预测Richness的值，为fit1

RIKZ$fit1<-predict(mod_lmm1)
head(RIKZ)
##   Sample Richness Exposure    NAP Beach      fit1
## 1      1       11       10  0.045     1  9.087834
## 2      2       10       10 -1.036     1 11.864274
## 3      3       13       10 -1.336     1 12.634793
## 4      4       11       10  0.616     1  7.621278
## 5      5       10       10 -0.684     1 10.960197
## 6      6        8        8  1.190     2  8.725105

绘制线性回归图

ggplot(RIKZ,aes(NAP,Richness,color=factor(Beach)))+
  geom_point(size=3,alpha=0.25)+
  geom_abline(aes(intercept=6.5819),slope=-2.5684)+
  geom_point(aes(y=fit1),size=5,alpha=0.5)+
  geom_line(aes(y=fit1),size=1)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_color_brewer(palette = "Set1")
## Warning: geom_abline(): Ignoring `mapping` because `slope` and/or `intercept`
## were provided.

斜率相同，但截距不同，与线性回归模型进行比较

RIKZ$fit0<-predict(mod_lm1)
head(RIKZ)
##   Sample Richness Exposure    NAP Beach      fit1      fit0
## 1      1       11       10  0.045     1  9.087834  6.556653
## 2      2       10       10 -1.036     1 11.864274  9.655722
## 3      3       13       10 -1.336     1 12.634793 10.515778
## 4      4       11       10  0.616     1  7.621278  4.919680
## 5      5       10       10 -0.684     1 10.960197  8.646589
## 6      6        8        8  1.190     2  8.725105  3.274107

一般线性回归图

ggplot(RIKZ,aes(NAP,Richness,color=factor(Beach)))+
  geom_point(size=3,alpha=0.25)+
  geom_abline(aes(intercept=6.5819),slope=-2.5684)+
  geom_point(aes(y=fit0),size=5,alpha=0.5)+
  geom_line(aes(y=fit0),size=1)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_color_brewer(palette = "Set1")
## Warning: geom_abline(): Ignoring `mapping` because `slope` and/or `intercept`
## were provided.

线性混合模型更能体现数据分布特征，添加随机斜率项目

head(RIKZ,10)
##    Sample Richness Exposure    NAP Beach      fit1      fit0
## 1       1       11       10  0.045     1  9.087834  6.556653
## 2       2       10       10 -1.036     1 11.864274  9.655722
## 3       3       13       10 -1.336     1 12.634793 10.515778
## 4       4       11       10  0.616     1  7.621278  4.919680
## 5       5       10       10 -0.684     1 10.960197  8.646589
## 6       6        8        8  1.190     2  8.725105  3.274107
## 7       7        9        8  0.820     2  9.675413  4.334842
## 8       8        8        8  0.635     2 10.150567  4.865210
## 9       9       19        8  0.061     2 11.624829  6.510784
## 10     10       17        8 -1.334     2 15.207746 10.510044
mod_lmm2<-lme(Richness~NAP,random=~1+NAP|Beach,data=RIKZ)

summary(mod_lmm2)
## Linear mixed-effects model fit by REML
##  Data: RIKZ 
##        AIC      BIC    logLik
##   244.3839 254.9511 -116.1919
## 
## Random effects:
##  Formula: ~1 + NAP | Beach
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 3.549064 (Intr)
## NAP         1.714963 -0.99 
## Residual    2.702820       
## 
## Fixed effects: Richness ~ NAP 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  6.588706  1.264761 35  5.209448   0e+00
## NAP         -2.830028  0.722940 35 -3.914610   4e-04
##  Correlation: 
##     (Intr)
## NAP -0.819
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.8213326 -0.3411043 -0.1674617  0.1921216  3.0397049 
## 
## Number of Observations: 45
## Number of Groups: 9

绘制混合线性模型回归曲线

RIKZ$fit2<-predict(mod_lmm2)

ggplot(RIKZ,aes(NAP,Richness,color=factor(Beach)))+
  geom_point(size=3,alpha=0.25)+
  geom_abline(aes(intercept=6.5819),slope=-2.5684)+
  geom_point(aes(y=fit2),size=5,alpha=0.5)+
  geom_line(aes(y=fit2),size=1)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_color_brewer(palette = "Set1")
## Warning: geom_abline(): Ignoring `mapping` because `slope` and/or `intercept`
## were provided.

可以看到不同NAP组的斜率及结局都不尽相同，与其他模型相比，更能体现数据特征

补充 lme4包实现

library(lme4)
mod_lmm4<-lmer(Richness~NAP+(1+NAP|Beach),REML=F,data = RIKZ)
## boundary (singular) fit: see ?isSingular
summary(mod_lmm4)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Richness ~ NAP + (1 + NAP | Beach)
##    Data: RIKZ
## 
##      AIC      BIC   logLik deviance df.resid 
##    246.7    257.5   -117.3    234.7       39 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7985 -0.3418 -0.1827  0.1749  3.1389 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  Beach    (Intercept) 10.949   3.309         
##           NAP          2.502   1.582    -1.00
##  Residual              7.174   2.678         
## Number of obs: 45, groups:  Beach, 9
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   6.5818     1.1883  8.8936   5.539 0.000377 ***
## NAP          -2.8293     0.6849  7.9217  -4.131 0.003366 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## NAP -0.810
## convergence code: 0
## boundary (singular) fit: see ?isSingular
AIC(mod_lmm1,mod_lmm2,mod_lmm4)
##          df      AIC
## mod_lmm1  4 247.4802
## mod_lmm2  6 244.3839
## mod_lmm4  6 246.6561

mod_lmm2的AIC最小，模型最好 注意：lme和lmer函数添加随机项目的表达式有点差别，同时，lmer函数中默认REML=TRUE）

一个人旅行*-*
关注
专栏目录
r语言 最优线性无偏预测_R语言如何解决线性混合模型中畸形拟合(Singular fit)的问题...
weixin_39700625的博客
11-23 576
原文链接：如何解决线性混合模型中畸形拟合(Singular fit)的问题​tecdat.cn假设我们有一个模型mod <- Y ~ X*Condition + (X*Condition|subject)# Y = logit variable # X = continuous variable # Condition = values A and B, dummy coded; the...
混合线性效应模型
04-23
混合线性模型的应用 介绍混合线性模型的结构，固定效应项和随机效应的含义。对具有内部相关性的资料，宜选用混合线性模型进行配合。方法：用一个具有聚集性结构的例子和一个重复测量的例子说明混合线性模型的方法和步骤。
基于R语言混合效应模型（mixed model）案例研究
拓端研究室TRL
11-17 1911
全文链接：http://tecdat.cn/?p=2596在本文中，我们描述了灵活的竞争风险回归模型。回归模型被指定为转移概率，也就是竞争性风险设置中的累积发生率（点击文末“阅读原文”获取完整代码数据）。相关视频1.混合模型是否适合您的需求？混合模型在很多方面与线性模型相似。它估计一个或多个解释变量对因变量的影响。混合模型的输出将为解释值列表，它们的效果大小的估计值和置信区间，每种效果的p值以及...
关于R语言多水平线性回归分析
qq_46191800的博客
09-11 5362
我再使用R语言分析嵌套模型的时候，遇到了一下问题，希望可以有大神指导一下。 数据是学校，班级，学生组成的嵌套模型。在模型构建的开始，我先确定模型中应该使用几水平。首先我先把三个水平均作为随机项放进模型中 model2.1<-lmer(disease~ age+gender+height+weight+ (1|school)+(1|class)+(1|stu), data=mydata,REML = FALSE, control = lmerControl( optimizer = “optimx”,
R语言数据分析报错解决办法和有用命令（8.3-8.10）
meruru的博客
08-10 4770
报错 pls 报错： 是否混淆pls，plsda， gplsda等，每个model所需输入不同 lmer： output 和input可以不在一个data，但行数必须相同 boundary (singular) fit: see ?isSingular fixed-effect model matrix is rank deficient so dropping 1 column / coefficient： check every adjust variable anova.mer Mod(mod_nul
全同胞家系如何计算遗传力及育种值
邓飞----育种数据分析之放飞自我
08-04 1895
全同胞家系，由同父同母所生子女的集合体称为全同胞家系。比如父本是A1，A2，A3，母本是B1，B2，B3，如果A1B1，A2A1*B1，有10个后代，分别是：A1B1_1, A1B1_2……，A1B1_10，那么这10个后代为一个全同胞家系上面共有30个个体，属于3个全同胞家系。...
广义线性混合模型pdf
11-09
### 广义线性混合模型概述 广义线性混合模型（Generalized Linear Mixed Models, GLMMs）是一种统计分析工具，它结合了广义线性模型（Generalized Linear Model, GLM）与混合效应模型的优点，适用于处理复杂的、非...
R语言 线性混合效应模型实战案例
拓端研究室TRL
12-19 2448
最近我们被客户要求撰写关于性混合效应模型的研究报告，包括一些图形和统计输出。首先，请注意，围绕多层次模型的术语有很大的不一致性。例如，多层次模型本身可能被称为分层线性模型、随机效应模型、多层次模型、随机截距模型、随机斜率模型或集合模型。根据不同的学科、使用的软件和学术文献，许多这些术语可能指的是相同的一般建模策略。线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例时长12:13。
R语言广义线性混合模型GLMMs在生态学中应用可视化2实例合集|附数据代码
最新发布
拓端研究室TRL
04-03 800
全文链接 :https://tecdat.cn/?p=35607分析师：Kaizong Ye，Colin Ge在生态学研究领域，广义线性混合模型（Generalized Linear Mixed Models，简称GLMMs）是一种强大的统计工具，能够同时处理固定效应和随机效应，从而更准确地揭示生态系统中复杂关系的本质（点击文末“阅读原文”获取完整代码数据）。随着数据分析技术的不断发展，R语言已成...
R语言广义线性混合模型（GLMM）bootstrap预测置信区间可视化
拓端研究室TRL
03-29 1229
全文链接：https://tecdat.cn/?p=35552通过线性模型和广义线性模型（GLM），预测函数可以返回在观测数据或新数据上预测值的标准误差（点击文末“阅读原文”获取完整代码数据）。相关视频然后，利用这些标准误差绘制出拟合回归线周围的置信区间或预测区间。置信区间（CI）的重点在于回归线，其可以解释为（假设我们绘制的是95%的置信区间）：“如果我们重复抽样X次，那么回归线将有95%的概率...
【视频】线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例
大数据部落
03-03 5097
原文链接：http://tecdat.cn/?p=23050 原文出处：拓端数据部落公众号 如果您熟悉线性模型，意识到它们的局限，那么您应该学习线性混合模型mixed-model。本视频中，我们讨论了线性混合模型并在R软件中进行应用。 视频：线性混合效应模型(LMM,Linear Mixed Models)和R语言实现 线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例 时长12:13 什么是混合效应建模，为什么要使用？ 统计分析中许多问题的传统方法.
蛛网模型收敛条件_混合线性模型的实现（更新20190607）
weixin_30609253的博客
12-28 1076
本文最早发布在本人的GitHub上，后来在R语言中文社区的公共号上发布过。在之后对其内容进行过几次更新，这一版为最新版，修改了一些错误的地方（如调整比较方式部分），增添了新的内容（随机斜率取舍部分）。R 中混合线性模型可依靠lme4或者lmerTest数据包（强烈推荐后者，因为会输出显著性）library(lmerTest)基本表达式fit = lmer(data = , formula = DV...
R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例
拓端研究室TRL
02-15 597
原文链接：http://tecdat.cn/?p=23050在本文中，我们将用R语言对数据进行线性混合效应模型的拟合，然后可视化你的结果。线性混合效应模型是在有随机效应时使用的，随机效应发生在对随机抽样的单位进行多次测量时。来自同一自然组的测量结果本身并不是独立的随机样本。因此，这些单位或群体被假定为从一个群体的 "人口 "中随机抽取的。示例情况包括当你划分并对各部分进行单独实验时（随机组）。当你...
R统计绘图-线性混合效应模型详解(理论、模型构建、检验、选择、方差分解及结果可视化)
qq_39859424的博客
09-24 2万+
​此文主要涉及线性混合效应模型分析，主要包含以下几部分内容： 1）混合模型基础知识； 2）线性混合效应模型构建(lme4)、检验、解读及可视化； 随机截距与随机斜率模型 3）线性混合效应模型方差分解及可视化； 4）线性混合效应模型比较及选择；
R语言︱线性混合模型理论与案例探究（固定效应&随机效应）
热门推荐
素质云笔记
06-11 11万+
线性混合模型与普通的线性模型不同的地方是除了有固定效应外还有随机效应。 ___________________________________________________________________________________ 一、线性混合模型理论 由两个部分来决定，固定效应部分+随机效应部分。 与普通线性回归模型以及广义线性模型
Advanced ShaderLab topics
wolf96的专栏
09-10 2531
Implementing Fixed Function TexGen in Shaders Before Unity 5, texture properties could have options inside thecurly brace block, e.g. TexGen CubeReflect. These were controlling fixed functionte
R 语言【ashr】warning: solve(): system is singular (rcond:****)； attempting approx solution
鬼话
06-26 675
最近更新了R，升级到了4.2，在使用 ashr 对 DESeq2 结果进行 LFC shrinkage 的时候，突然出现了warning……warning: solve(): system is singular (rcond:****); attempting approx solution
R 多变量数据预处理_用R进行多层中介的完整流程：从数据预处理到置信区间的生成...
weixin_39670849的博客
11-22 1647
1. 中介与多层中介1.1 中介效应中介效应(mediation effect)指的是变量X对Y的影响，是否经过另一个中介变量M。在上图中，a表示的是自变量X到中介变量M的回归系数，b表示的是从中介变量到因变量Y的回归系数，而c’表示的是去除中介效应后，自变量X对因变量Y剩余的直接效应。在中介分析中：传统的中介通常分析遵循Baron和Kenny(1986)的四步法：1) 忽略中介变量M...
线性混合模型
XIAOFEI@IDO
02-27 8871
背景 线性模型需要满足正态性、独立性和同方差性等假设，其中独立性是线性模型最重要的假设之一，独立性要求每一个数据点必须来自于不同的总体。但由于重复测量数据、区组数据以及空间相关数据不能满足独立性假设，因此常常利用线性混合效应模型对上述数据进行分析（如相关性分析）。 以一个探索基因表达量与采样距离相关性为例。 该数据来自于“ntegrated transcriptomic analysis of distance-related field cancerization in rectal cancer pat
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值