💂 个人信息：酷在前行👍 版权: 博文由【酷在前行】原创、需要转载请联系博主👀 如果博文对您有帮助，欢迎点赞、关注、收藏 + 订阅专栏🔖 本文收录于【R模型】，该专栏主要介绍R语言各类型机器学习，如线性回归模型、广义线性模型、混合线性模型、随机森林模型、支持向量机模型、决策树模型等。请大家多多关注点赞和支持，共同进步~ 欢迎大家订阅！
📋 文章目录
🐣 一、理解比较两个模型斜率差异的意义🐤 二、混合线性模型的R语言过程 🍎 1、方法1：利用simba包中的diffslope函数 🍌 2、方法2：在构建模型时添加固定变量与分类变量的交互项，判断交互项是否显著 🥝 3、方法3：定义一个新的函数来实现多类型模型及解决n阶循环的问题
📇 三、附上所有代码
想要判断两个模型是否存在显著性差异，即可以通过检验两者模型的斜率是否存在显著性差异来判断。这里，推荐了三种方式：1.通过R中
“simba包” 检验，但是存在其一弊端（似乎只能针对线性回归模型，并且不便于循环）；2.通过构建的分类模型中，
添加固定变量与分类变量的交互项 来判断，其不足为不能够获得置信区间；3.我重新定义一个能够识别
多类模型（包括线性回归模型、广义线性模型、混合线性模型和广义混合线性模型），并且能够很好地解决n阶循环的情况 ，函数名：
cal_linear_boot，其原理是利用bootstrap重采样对模型构建的斜率进行多次重采样的方式或者两个分类模型下的N次（通常N=999），得到两个分类的数据组，最后通过选择合适的统计方法来检验是否存在显著差异。但是再完美的方法也有不足的地方，可能需要时间较长，请您的耐心等待。
🐣 一、理解比较两个模型斜率差异的意义
首先，本篇推文沿用探究混合线性模型不同类型博文中的数据集来讲解R过程。这里先给大家解读一个案例来理解比较两个模型斜率显著性的意义：假如探究男生和女生的考试成绩受作业完成度的影响是否存在差异。毫无疑问，大家肯定能够想到通过判断模型中的固定变量对响应变量的系数及其p值来解读其中的差异；但是，这只能通过比较两者系数的大小来说明相对差异。如果要判断两个模型中响应变量随固定变量变化趋势的差异，那我们应该要检验两个模型的曲线斜率是否存在显著性，从而通过统计学角度量化和确定显著性差异。
针对不同的数据集及变量类型，采用不同的模型。比如针对上述问题我们可以考虑使用线性回归模型或者是广义线性模型（lm和glm）；如果数据集包括全国各个地区的情况，那这时候可能就存在地区差异了，比如在某一个地区男生和女生的差异与另外一个地区是相反的，所以可以将模型替换成混合线性模型或广义混合线性模型（lmer和glmer），来避免这种数据集的混乱。
🐤 二、混合线性模型的R语言过程
🍎 1、方法1：利用simba包中的diffslope函数
先加载R包和查看对应函数
# 加载R包
library(simba)
library(tidyverse)
# 查看diffslope函数介绍
?diffslope
主要参数的介绍：
x1 : 自变量x1的向量Y1 : 因变量Y2的向量，包含依赖于x1的变量(例如图之间的相似性。必须和x1的长度相同）x1 : 包含第二个自变量的X2向量(例如图之间的距离)。可以和x1一样。Y2 : 因变量Y2向量，包含依赖于x2的变量(例如图之间的相似性。必须和x2的长度相同）。
# 读取数据
mixdata <- read.csv("mixeddata.csv", header = T)
# 添加两列数据
mixdata$group1 <- sample(1:5, size = nrow(mixdata), replace = T)
mixdata$group2 <- sample(0:1, size = nrow(mixdata), replace = T)
# 查看数据结构
str(mixdata)
# 计算location中a与b分类下，广告停留时间与年龄的关系的显著性差异检验
diffslope((mixdata %>% filter(location == "a"))$age,
(mixdata %>% filter(location == "a"))$bounce_time,
(mixdata %>% filter(location == "b"))$age,
(mixdata %>% filter(location == "b"))$bounce_time)
从检验的结果来看，两个模型的斜率不存在显著性差异。其中结果显示了斜率差异的结果，还有对应的统计检验p值大小（分别为0.01165和0.423），和对应置信区间的斜率值。表明，如果要达到显著，应该斜率的差值需要大于0.1076。所以，从这个方法来看，很好理解，也很快的得到对应的结果。但似乎不能用于其他模型（仅限线性回归模型）。
🍌 2、方法2：在构建模型时添加固定变量与分类变量的交互项，判断交互项是否显著
# 构建线性回归模型及查看摘要
lm <- lm(bounce_time ~ age*location, data = mixdata)
summary(lm)
其他类型的模型：
# 以混合效应模型为例
library(lme4)
library(lmerTest)
lmer <- lmer(bounce_time ~ age*location + (1|county ), data = mixdata)
summary(lmer)
🥝 3、方法3：定义一个新的函数来实现多类型模型及解决n阶循环的问题
该函数中n阶循环是指可能你的数据集存在多种处理，也就是说存在不同处理变量间组合成新的子数据集。该函数的输出结果返回一个数据框，包含不同处理及bootstrap重采样的模型斜率的结果，最后可以利用统计学检验不同处理变量间的显著性差异。 第一 ，可以实现不同类型模型的斜率比较； 第二 ，可以解决不同变量处理间的多阶循环。
cal_linear_boot <- function(data,
R = 999,
sub = NULL,
methods = "lmer",
familys = NULL,
equa = bounce_time ~ age + ( 1 | county)) {
required_packages <- c("foreach", "doParallel", "parallel", "boot", "lme4", "lmerTest", "tidyverse")
install.lib <- required_packages[!required_packages %in% installed.packages()]
for(libs in install.lib) install.packages(libs, dependences = TRUE)
sapply(required_packages, require, character = TRUE) ###import the packages which are need af
stopifnot(is.data.frame(data))
stopifnot(is.numeric(R))
stopifnot(methods %in% c("lmer", "lm", "glm", "glmer"))
if (!is.null(sub)) {
data <- unite(data, sub, unique(sub))
}
cl <- makeCluster(detectCores()-2)
registerDoParallel(cl)
beta <- function(formula, data, indices){
d <- data[indices,]
fit <- tryCatch({
switch(methods,
"lmer" = lmer(equa, data = d),
"lm" = lm(equa, data = d),
"glm" = glm(equa, family = familys, data = d),
"glmer" = glmer(equa, family = familys, data = d))
},
error = function(e) {
message("Error: ", e$message)
NA
})
if (!is.na(fit)) {
return(summary(fit)$coefficients[2])
} else {
return(NA)
}
}
parallel_fun <- function(data, formula) {
r <- boot(data = data, statistic = beta, R = R, formula = equa)$t
treat <- rep(unique(data$sub), length(r))
data.frame(treat = treat, slope = r)
}
data_list <- split(data, data$sub)
result_list <- foreach(i = 1:length(data_list), .combine = rbind,
.packages = c("boot", "lme4", "lmerTest", "doParallel", "parallel") ) %dopar% {
parallel_fun(data_list[[i]])
}
stopCluster(cl)
final_result <- data.frame(result_list)
if (!is.null(sub)) {
final_result[, sub] <- t(as.data.frame(strsplit(as.character(final_result$treat), "_")))
}
return(final_result)
}
主要参数的介绍：
data : 数据集R : bootstrap的次数sub : 需要对哪些变量进行划分成子数据集的循环equa : 对应构建模型的fomulamethods : 使用哪个模型，包括lm、glm、lmer、glmerfamilys : 当使用glm和glmer时，需要定义响应变量的分布类型，可以参考glm中family参数的介绍
# source("cal_linear_boot.R")
result <- cal_linear_boot(data = mixdata, R = 999, sub = c("location", "group1", "group2"),
equa = bounce_time ~ age + (1|county), methods = "lmer")
head(result)
对结果进行可视化：
library(ggplot2)
library(ggsignif)
# 设置对比顺序
compaired <- list(c("1","2"),c("1","3"),c("1","4"),c("1","5"),
c("2","3"),c("2","4"),c("2","5"),c("3","4"),
c("4","5"))
# 结果可视化并进行统计学检验
ggplot(result, aes(group1, slope, fill = group2))+
geom_boxplot()+
facet_wrap(~ location, scale = "free")+
theme_bw()+
geom_signif(y_position = c(1.1, 1.4, 1.7, 2.0, 2.3, 2.6, 2.9, 3.2, 3.5),
comparisons = compaired,
map_signif_level = T,
test = t.test,
tip_length = 0,
vjust = 0)
可视化结果：
📇 三、附上所有代码
library(simba)
library(tidyverse)
?diffslope
mixdata <- read.csv("mixeddata.csv", header = T)
mixdata$group1 <- sample(1:5, size = nrow(mixdata), replace = T)
mixdata$group2 <- sample(0:1, size = nrow(mixdata), replace = T)
str(mixdata)
diffslope((mixdata %>% filter(location == "a"))$age,
(mixdata %>% filter(location == "a"))$bounce_time,
(mixdata %>% filter(location == "b"))$age,
(mixdata %>% filter(location == "b"))$bounce_time)
lm <- lm(bounce_time ~ age*location, data = mixdata)
summary(lm)
library(lme4)
library(lmerTest)
lmer <- lmer(bounce_time ~ age*location + (1|county ), data = mixdata)
summary(lmer)
cal_linear_boot <- function(data,
R = 999,
sub = NULL,
methods = "lmer",
familys = NULL,
equa = bounce_time ~ age + ( 1 | county)) {
required_packages <- c("foreach", "doParallel", "parallel", "boot", "lme4", "lmerTest", "tidyverse")
install.lib <- required_packages[!required_packages %in% installed.packages()]
for(libs in install.lib) install.packages(libs, dependences = TRUE)
sapply(required_packages, require, character = TRUE) ###import the packages which are need af
stopifnot(is.data.frame(data))
stopifnot(is.numeric(R))
stopifnot(methods %in% c("lmer", "lm", "glm", "glmer"))
if (!is.null(sub)) {
data <- unite(data, sub, unique(sub))
}
cl <- makeCluster(detectCores()-2)
registerDoParallel(cl)
beta <- function(formula, data, indices){
d <- data[indices,]
fit <- tryCatch({
switch(methods,
"lmer" = lmer(equa, data = d),
"lm" = lm(equa, data = d),
"glm" = glm(equa, family = familys, data = d),
"glmer" = glmer(equa, family = familys, data = d))
},
error = function(e) {
message("Error: ", e$message)
NA
})
if (!is.na(fit)) {
return(summary(fit)$coefficients[2])
} else {
return(NA)
}
}
parallel_fun <- function(data, formula) {
r <- boot(data = data, statistic = beta, R = R, formula = equa)$t
treat <- rep(unique(data$sub), length(r))
data.frame(treat = treat, slope = r)
}
data_list <- split(data, data$sub)
result_list <- foreach(i = 1:length(data_list), .combine = rbind,
.packages = c("boot", "lme4", "lmerTest", "doParallel", "parallel") ) %dopar% {
parallel_fun(data_list[[i]])
}
stopCluster(cl)
final_result <- data.frame(result_list)
if (!is.null(sub)) {
final_result[, sub] <- t(as.data.frame(strsplit(as.character(final_result$treat), "_")))
}
return(final_result)
}
source("cal_linear_boot.R")
result <- cal_linear_boot(data = mixdata, R = 999, sub = c("location", "group1", "group2"),
equa = bounce_time ~ age + (1|county), methods = "lmer")
head(result)
library(ggplot2)
library(ggsignif)
compaired <- list(c("0","1"))
ggplot(result, aes(group2, slope, fill = group1))+
geom_boxplot()+
facet_wrap(~ location, scale = "free")+
theme_bw()+
geom_signif(
comparisons = compaired,
map_signif_level = T,
test = t.test,
tip_length = 0,
vjust = 0)