线性混合效应模型原理及应用-CSDN博客

本文链接：https://blog.csdn.net/weixin_45803246/article/details/134702363

文章目录

一般线性模型
广义线性模型
线性混合效应模型
广义线性混合模型
参考

一般线性模型

一般线性模型（general linear model，GLM）是一种常见的统计学习模型，基于一组自变量（预测变量）和一个连续型响应变量（因变量），利用最小二乘法进行回归分析从而得到各个自变量与因变量之间的关系。

#设置工作目录并加载包
setwd("C:\\Users\\Dell\\Desktop\\zp\\学习\\统计\\mixed effect model workshop-2019-12-21-22\\Mixed effect model-ZP") ##依据自己数据的实际存放位置设置工作目录

library(tidyverse)
library(car)
library(mgcv)
library(ciTools)
library(nlme)
library(AER)
library(MASS)
library(pscl)
library(lme4)
library(lmerTest)
library(MuMIn)
library(RLRsim)
library(sjPlot)
library(GLMMadaptive)

一般线性回归

RIKZ <- read.delim("RIKZ.txt")

plot(Richness~NAP,data=RIKZ, pch=19)
RIKZ<-RIKZ[order(RIKZ$NAP),]
lm0<-lm(Richness~NAP,data=RIKZ)
summary(lm0)
abline(lm0,col="red")
anova(lm0)

一般线性模型与t检验和方差分析的联系

### 结果解读
female<-rnorm(1000,100, 10)
male<-rnorm(1000,125,10)
weight<-(c(female,male))
d1<-as.data.frame(weight)
d1$sex[c(0:1000,0)]<-"female"
d1$sex[c(1001:2000,0)]<-"male"
View(d1)

lm1 <- lm(weight~sex-1, data =d1)
summary(lm1)  ###  单因素分类变量的线性模型结果等价于单因素的t检验
t.test(weight~1, data =d1 %>% filter(sex=='female')) ## 截距
#t.test(weight~1, data =d1 %>% filter(sex=='male')) ## 截距
t.test(weight~sex, data =d1)  ### 斜率
aov(weight~sex, data =d1) %>% summary()

模型评估

正态性检验

hist(lm0$residuals)
##对lm0模型的残差的正态性检验
shapiro.test(lm0$residuals)
#定量检验：P > 0.05， 符合正态

qqPlot(lm0)
#定性检验，大部分点落在置信区间，符合正态性

方差齐性检验


##方差齐次性检验--方法一
boxplot(weight~sex, data =d1, ylab="Weight")
bartlett.test(weight~sex, data =d1)

####方差齐次性检验--方法二
plot(lm1$fitted.values,lm1$residuals,main="Residuals VS fitted values")
abline(lm(lm1$residuals~lm1$fitted.values),col="red")

## 连续变量方差齐性检验
ncvTest(lm0)  ## P > 0.05, 方差齐性
spreadLevelPlot(lm0)  ## 拟合线没有明显趋势表明方差齐性

独立性检验

durbinWatsonTest(lm0)  ## P > 0.05, 误差独立

异常观测值或离群点

outlierTest(lm0) ## P > 0.05, 没有异常值

高杠杆值

hat.plot <- function(fit){
  p <- length(coefficients(fit))
  n <- length(fitted(fit))
  plot(hatvalues(fit),main = "Index Plot of Hat Value")
  abline(h=c(2,3)*p/n,col="red",lty=2)
  identify(1:n,hatvalues(fit),names(hatvalues(fit)))
}
hat.plot(lm0)

强影响点

cutoff <- 4/(nrow(RIKZ)-length(lm0$coefficients)-2)
plot(lm0,which = 4,cook.levels = cutoff)
abline(h=cutoff,lty=2,col="red")  ## 红线以上的是强影响点

influencePlot(lm0,id.method="identity",main="Influence Plot",
              sub="Circle size is proportional to Cook's distance")

模型优化

正态性解决–数据转换


#lm2 <- lm(log(Richness)~NAP,data=RIKZ)
lm2 <- lm((Richness)^(1/3)~NAP,data=RIKZ)
shapiro.test(lm2$residuals)
qqPlot(lm2)

方差齐性问题解决 – 指定方差变异趋势

##1) 允许方差随着解释变量的增大而增大(适用于连续变量)
library(nlme)
Squid <- read.delim("Squid.txt")
Squid$fMONTH=factor(Squid$MONTH)

M.gls1<-gls(Testisweight~DML*fMONTH,
            weights=varFixed(~DML),data=Squid)
##2)允许方差随着随着解释变量（指分类变量）水平的不同而不同呢
M.gls2 <- gls(Testisweight ~ DML*fMONTH,
              weights=varIdent(form= ~ 1|fMONTH), data =Squid)

##3) 组合型方差结构：允许同时方差随着连续变量和分类变量而变化
M.gls3<-gls(Testisweight ~ DML * fMONTH,
            weights =varComb(varIdent(form= ~ 1 | fMONTH) ,
                             varExp(form =~ DML) ), data = Squid)

summary(M.gls1)
summary(M.gls2)
summary(M.gls3)

lm3 <- gls(Richness~NAP,weights=varFixed(~NAP),data=RIKZ)
summary(lm3)

非线性关系

par(mfrow = c(1,1))
isit <- read.delim("ISIT.txt")
View(isit)
isit<-isit[order(isit$SampleDepth),]
plot(Sources~SampleDepth, data=isit)

可加模型

am <- gam(Sources~s(SampleDepth),data=isit) ## 可加模型
summary(am)
plot(Sources~SampleDepth,data=isit)
am_ci <- add_ci(isit, am, alpha = 0.05, names = c("lwr", "upr"))
lines(am_ci$SampleDepth, am_ci$pred, col = "blue", lwd = 2,lty="solid")
lines(am_ci$SampleDepth, am_ci$upr, col = "blue", lwd = 2,lty="dotted")
lines(am_ci$SampleDepth, am_ci$lwr, col = "blue", lwd = 2,lty="dotted")

AIC(lm0,lm2)

广义线性模型

广义线性模型（Generalized Linear Model，GLIM）是一种基于广义概率分布进行建模、以预测某个连续或离散响应变量的统计学习方法。与传统的线性回归模型不同，GLIM可以同时处理连续型和离散型数据，并且有更大的灵活性来拟合各种类型的响应变量。

GLIM模型结构包括三部分：随机部分、系统部分和链接函数。其中随机部分指的是响应变量所服从的概率分布，例如二项式分布、泊松分布、正态分布等；系统部分则表示关键自变量及其系数，在这里可使用线性或非线性函数进行自变量转换；而链接函数则是将随机部分和系统部分连接起来的桥梁，用来确保随机变量具有最佳拟合度并满足带约束的条件。

相对于传统的线性模型，GLM 主要具有以下三个特点：

可以处理不符合正态分布的数据，例如二项式、泊松、正态和伽马等分布。
可以实现非线性函数的变换，使得建模能够更好地拟合数据，并且可以探索因素的交互影响。
能够处理离散因素和定类数值因素，同时也支持多重共线性太高的情况。
可以这样理解，GLM 将响应变量看做是某个连续性随机变量的结果，这个随机变量服从某个概率分布，在此基础上，建立广义的多元线性回归模型，并通过最大似然函数等方法来对参数进行估计，进而分析不同变量对响应变量的影响程度和形态关系。

泊松分布

RIKZ<-RIKZ[order(RIKZ$NAP),]
plot(Richness~NAP,data=RIKZ, pch=19)

计数数据：泊松回归案例

hist(RIKZ$Richness)
pm1<- glm(Richness~NAP, family="poisson", data=RIKZ)
summary(pm1)

plot(Richness~NAP,data=RIKZ, pch=19)
ci_pm1<-add_ci(RIKZ, pm1, alpha = 0.05, names = c("lwr", "upr"))
lines(ci_pm1$NAP, ci_pm1$pred, col = "red", lwd = 1,lty="solid")
lines(ci_pm1$NAP, ci_pm1$upr, col = "red", lwd = 1,lty="dotted")
lines(ci_pm1$NAP, ci_pm1$lwr, col = "red", lwd = 1,lty="dotted")
title("glm-poisson regression")

模型诊断

plot(pm1)
plot(pm1$residuals~pm1$fitted.value,main="Model diagnosis for glm-poisson model")
abline(lm(pm1$residuals~pm1$fitted.values),col="red")
summary(lm(pm1$residuals~pm1$fitted.values))

查看线性回归与泊松回归的AIC值

AIC(lm0,pm1)

查看Richness的均值与方差

mean(RIKZ$Richness)
var(RIKZ$Richness)

泊松分布的均值为5.688889时的理论分布

poisson2<-rpois(44,5.688889)
hist(poisson2)
hist(RIKZ$Richness)  ## 实际分布

对过度离散的检验

pm1<- glm(Richness~NAP, family="poisson", data=RIKZ)
dispersiontest(pm1)