devbox
运维做开发
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

线性混合模型R实现的更多实例_boundary (singular) fit: see help('issingular')

作者：运维做开发 | 2024-08-01 04:56:48

赞

踩

boundary (singular) fit: see help('issingular')
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
  1. ## 
  2. ## Attaching package: 'lme4'
  1. ## The following object is masked from 'package:nlme':
  2. ## 
  3. ##     lmList
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 3.6.3
  1. ## 
  2. ## Attaching package: 'lmerTest'
  1. ## The following object is masked from 'package:lme4':
  2. ## 
  3. ##     lmer
  1. ## The following object is masked from 'package:stats':
  2. ## 
  3. ##     step
  1. library(AED)
  2. library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.6.3

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

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

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

绘制model1的诊断性曲线图

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

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

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

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

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

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

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

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

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

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

  1. mod_lmm1<-lme(Richness ~NAP,random = ~1|Beach,data=RIKZ)
  2. summary(mod_lmm1)
  1. ## Linear mixed-effects model fit by REML
  2. ##  Data: RIKZ 
  3. ##        AIC     BIC    logLik
  4. ##   247.4802 254.525 -119.7401
  5. ## 
  6. ## Random effects:
  7. ##  Formula: ~1 | Beach
  8. ##         (Intercept) Residual
  9. ## StdDev:    2.944065  3.05977
  10. ## 
  11. ## Fixed effects: Richness ~ NAP 
  12. ##                 Value Std.Error DF   t-value p-value
  13. ## (Intercept)  6.581893 1.0957618 35  6.006682       0
  14. ## NAP         -2.568400 0.4947246 35 -5.191574       0
  15. ##  Correlation: 
  16. ##     (Intr)
  17. ## NAP -0.157
  18. ## 
  19. ## Standardized Within-Group Residuals:
  20. ##        Min         Q1        Med         Q3        Max 
  21. ## -1.4227495 -0.4848006 -0.1576462  0.2518966  3.9793918 
  22. ## 
  23. ## Number of Observations: 45
  24. ## Number of Groups: 9

预测Richness的值，为fit1

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

绘制线性回归图

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

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

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

一般线性回归图

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

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

head(RIKZ,10)
  1. ##    Sample Richness Exposure    NAP Beach      fit1      fit0
  2. ## 1       1       11       10  0.045     1  9.087834  6.556653
  3. ## 2       2       10       10 -1.036     1 11.864274  9.655722
  4. ## 3       3       13       10 -1.336     1 12.634793 10.515778
  5. ## 4       4       11       10  0.616     1  7.621278  4.919680
  6. ## 5       5       10       10 -0.684     1 10.960197  8.646589
  7. ## 6       6        8        8  1.190     2  8.725105  3.274107
  8. ## 7       7        9        8  0.820     2  9.675413  4.334842
  9. ## 8       8        8        8  0.635     2 10.150567  4.865210
  10. ## 9       9       19        8  0.061     2 11.624829  6.510784
  11. ## 10     10       17        8 -1.334     2 15.207746 10.510044
  1. mod_lmm2<-lme(Richness~NAP,random=~1+NAP|Beach,data=RIKZ)
  2.  
  3. summary(mod_lmm2)
  1. ## Linear mixed-effects model fit by REML
  2. ##  Data: RIKZ 
  3. ##        AIC      BIC    logLik
  4. ##   244.3839 254.9511 -116.1919
  5. ## 
  6. ## Random effects:
  7. ##  Formula: ~1 + NAP | Beach
  8. ##  Structure: General positive-definite, Log-Cholesky parametrization
  9. ##             StdDev   Corr  
  10. ## (Intercept) 3.549064 (Intr)
  11. ## NAP         1.714963 -0.99 
  12. ## Residual    2.702820       
  13. ## 
  14. ## Fixed effects: Richness ~ NAP 
  15. ##                 Value Std.Error DF   t-value p-value
  16. ## (Intercept)  6.588706  1.264761 35  5.209448   0e+00
  17. ## NAP         -2.830028  0.722940 35 -3.914610   4e-04
  18. ##  Correlation: 
  19. ##     (Intr)
  20. ## NAP -0.819
  21. ## 
  22. ## Standardized Within-Group Residuals:
  23. ##        Min         Q1        Med         Q3        Max 
  24. ## -1.8213326 -0.3411043 -0.1674617  0.1921216  3.0397049 
  25. ## 
  26. ## Number of Observations: 45
  27. ## Number of Groups: 9

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

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

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

补充 lme4包实现

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

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

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/运维做开发/article/detail/912378
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号