此文介绍如何使用MuMIn包使用最优子集法进行多重线性回归的模型筛选以及模型平均。多重线性回归需要进行的数据检验过程都写在R统计绘图-多重线性回归(最优子集法特征筛选，leaps)中了。大家可以自行查看。

一、数据准备

# 设置工作路径
#knitr::opts_knit$set(root.dir="D:\\EnvStat\\EcoEvoPhylo\\MLR")# 使用Rmarkdown进行程序运行
Sys.setlocale('LC_ALL','C') # Rmarkdown全局设置
#options(stringsAsFactors=F)# R中环境变量设置，防止字符型变量转换为因子# 1.1 导入数据
library(tidyverse)
da <- read.csv("data.csv",row.names = 1)
da %>%head()da.scale <- da %>%mutate(across(.cols = !contains(c("condition","depth","grazing")),.fns = scale), # 对数值数据进行标准化grazing = factor(grazing, levels = c("CK","LG","MG","HG")) # 将因子变量设置为实验设计顺序。#若要将分类变量用作线性回归分析，一定要进行此设置，将对照组排在第一位，分析时其会作为哑变量，即，它在回归方程中的回归系数会为0。) %>% #summarise_if(is.numeric,scale) %>% # 对数值数据进行标准化,且不保留分类变量data.frame()
da.scale %>%head()

图1|环境因子及分组信息表，data.csv。行为样品名称，列为环境因子名称和分组信息，共有11个环境变量，3个分组信息。

二、 模型构建

最优子集回归：算法使用所有可能的特征组合来拟合模型。分析者需要应用自己的判断和统计分析来选择最优模型。当特征数多于观测数时，这个方法的效果就不会好。此处介绍如何使用MuMIn包进行最优子集回归，筛选特征，根据模型筛选准则筛选模型以及构建平均加权模型。

2.1 最优子集特征选择-MuMIn包

library(MuMIn)
# 2.1.1 全模型构建
fm1 <- lm(# 可在此处选择自变量pH ~ grazing + depth + TN + TP + TK + Nitrate + Ammonia + AHN + AP + AK + OM + OC,# 不能在此处筛选变量，否则后续构建模型子集会报错。# 即，不能使用da.scale[-1]等data = da.scale,# 不能设置为“na.omit"，会报错。na.action = "na.fail" )# 2.1.2 生成模型选择表
## 生成具有全局模型中固定效应项的组合(子集)的模型的模型选择表，具有可选的模型包含规则。
dd <- MuMIn::dredge(# 全模型fm1,# 模型子集的排序基准。rank = "BIC", # 回归系数标准化方法。beta = "partial.sd",# 设置输出的模型筛选指数。extra = alist(AICc,AIC, BIC, ICOMP, Cp,"R^2","adjR^2"),# 设置某变量包含在所有模型中。#fixed = "AP",# 设置模型筛选方式,只输出包含AP的模型子集。#subset = expression(AP))
dd %>%summary()# 2.1.3 模型筛选方法一：提取最佳拟合模型
fm_best <- get.models(dd, 1)[[1]]
fm_best %>%summary()

图2|基于BIC筛选的最优模型。输出结果中Estimate， Std. Error， t value, Pr(>|t|):依次表示对应自变量的回归系数及其标准误，回归系数检验的t检验统计量及其p值。F-statistic: 40.4 on 7 and 28 DF, p-value:5.502e-13 则是是对整个模型的显著性检验。分类因子作为自变量时，因子中的每个分类水平都会有一个回归系数，第一个分类水平作为哑变量，其回归系数为0。此模型中grazingCK和depth0-10cm的回归系数为0。

2.2 平均加权模型-MuMIn包

利用MuMIn包dredge()生成所有自变量可能组合出的模型，并根据BIC排序。然后使用MuMIn::model.avg()利用模型平均的方法获得所有满足筛选条件的模型加权后的最优模型参数估计。

# 2.2.1 生成所有模型统计量表-包含模型F检验结果
d2 <- MuMIn::dredge(fm1,# 设置模型包含的terms(回归系数)数目(截距除外)# (1,NA)表示模型中至少要有一个回归系数。m.lim = c(1, NA), rank = "BIC",extra = list("BIC","AIC","adjR^2","*" = function(x) {s <- summary(x)c(Rsq = s$r.squared, adjRsq = s$adj.r.squared,F = s$fstatistic[[1]])})
)
d2 %>%summary()# 2.2.2 模型筛选方法二：模型平均加权获得最优模型
##根据delta AICc筛选模型:ΔAICc<2的模型被认为具有相同的优良性。
## 提取delta AICc<2的所有模型
fm_d2 <- get.models(d2, delta < 2)
fm_d2 # 输出结果是一个模型列表数据。## 根据信息准则模型平均
##运行dredge()，未设置rank时，最好先使用get.models()提取模型，再传递给model.avg()，否则predict()等函数可能会对平均后模型不起作用。
##平均模型涉及因子变量时，需要注意，模型子集可能包含因子的不同分类，目前对此进行平均会产生毫无意义的结果。
fm_avg1 <- model.avg(fm_d2,beta = "partial.sd")fm_avg1 %>%summary()

图3|输出满足筛选条件的模型子集，fm_d2。

图4|输出满足筛选条件的模型子集，fm_avg1。平均加权模型中包含的模型子集以及"full"和““subset””平均模型的回归系数检验结果。

## 直接使用model.avg()进行模型筛选和平均，结果是一样的。
##dredge()未设置rank时，运行此函数设置fit=TRUE，否则平均加权后模型使用不了predict()等函数。
fm_avg2 <- model.avg(d2, subset = delta < 2,beta = "partial.sd",# 评估模型fit = TRUE)fm_avg2 %>%summary()## 其他筛选条件：根据95% confidence筛选
model.avg(dd, subset = cumsum(weight) <= .95) # get averaged coefficients# 2.2.3 将所有优良模型加权后的最优模型参数估计
## 平均加权模型汇总统计
fm_avg <- fm_avg2 %>%summary()fm_avg %>%str() # 查看模型中包含的信息## 自变量数据表
fm_avg$x## 每个变量被多少个模型包括及其权重
fm_avg$sw## 模型子集中包含的回归方程式，自由度和模型筛选信息准则。
fm_avg$msTable## 平均模型系数：包含“full"和”subset”两种形式。
##"full"平均表示计算平均回归系数时，假设变量在所有子模型中都存在，只是在部分模型(回归公式中不存在该变量的模型)中系数为0。
##"full"平均"与“subset”平均不同的是，它没有偏离零的倾向。“full”平均是一种收缩估计器，对于与响应变量关系较弱的子变量，它比“subset”平均得到的值更小。
##‘subset’ (‘conditional’) 平均表示只平均该变量回归系数存在的模型。
fm_avg$coefficients## "full"平均加权模型的变量回归系数、标准误，z检验值和p值。
fm_avg$coefmat.full## "subset"平均加权模型的变量回归系数、标准误，z检验值和p值。
fm_avg$coefmat.subset## 每个变量被多少个模型包括
fm_avg$coef.nmod## 运行公式
fm_avg$call## 反应式
fm_avg$formula ## 模型残差
fm_avg$residuals## 三维数组：包含每个模型参数的统计量，自由度和标准误差。
fm_avg$coefArray## “full"平均模型预测
predict(fm_avg2,full = TRUE)## 计算“full"平均模型置信区间
confint(fm_avg2,full = TRUE)## 提取“full"平均模型回归系数
predict(fm_avg2,full = TRUE)## 计算“full"平均模型的方差-协方差矩阵
vcov(fm_avg2,full = TRUE)

数据及代码可以QQ交流群文件夹中下载，或微信公众号后台发送“MLR2”获取。

原文链接：R统计绘图-多重线性回归(平均加权模型/最优子集筛选，MuMIn)

R绘图-物种、环境因子相关性网络图(简单图、提取子图、修改图布局参数、物种-环境因子分别成环径向网络图)

R统计绘图-分子生态相关性网络分析(拓扑属性计算，ggraph绘图)

R统计绘图-变量分组相关性网络图(igraph)

机器学习-分类随机森林分析(randomForest模型构建、参数调优、特征变量筛选、模型评估和基础理论等)

R统计绘图-多重线性回归(最优子集法特征筛选，leaps)

R统计绘图 | 物种组成堆叠柱形图(绝对/相对丰度)

R统计绘图 | 物种组成冲积图(绝对/相对丰度，ggalluvial)

R统计绘图 | 物种组成桑基图(不同分类单元微生物丰度及其从属关系，ggalluvial) -免费送书

R统计绘图 | 物种组成堆叠面积图(绝对/相对丰度，ggalluvial)

R统计绘图 | 物种组成气泡图和聚类热图(ggplot2/ComplexHeatmap)

R统计绘图-随机森林分类分析及物种丰度差异检验组合图

机器学习-多元分类/回归决策树模型(tree包)

R统计绘图-环境因子相关性+mantel检验组合图(linkET包介绍1)

R统计绘图-NMDS、环境因子拟合(线性和非线性)、多元统计(adonis2和ANOSIM)及绘图(双因素自定义图例)

R统计绘图-RDA分析、Mantel检验及绘图

R绘图-RDA排序分析R统计绘图-VPA(方差分解分析)

R统计绘图-PCA详解1(princomp/principal/rcomp/rda等)

R统计-PCA/PCoA/db-RDA/NMDS/CA/CCA/DCA等排序分析教程

R统计绘图-PCA分析绘图及结果解读(误差线，多边形，双Y轴图、球形检验、KMO和变量筛选等)

R统计-微生物群落结构差异分析及结果解读

R统计绘图-PCA分析及绘制双坐标轴双序图

R统计绘图-分子生态相关性网络分析

R中进行单因素方差分析并绘图

R统计-多变量单因素参数、非参数检验及多重比较

R绘图-相关性分析及绘图

R绘图-相关性系数图

R统计绘图-环境因子相关性热图

R统计绘图-corrplot绘制热图及颜色、字体等细节修改

R统计绘图-corrplot热图绘制细节调整2(更改变量可视化顺序、非相关性热图绘制、添加矩形框等)

R数据可视化之美-节点链接图

R统计绘图-rgbif包下载GBIF数据及绘制分布图

R统计-单因素ANOVA/Kruskal-Wallis置换检验

R统计-正态性分布检验[Translation]

R统计-数据正态分布转换[Translation]

R统计-方差齐性检验[Translation]

R统计-Mauchly球形检验[Translation]

R统计绘图-单、双、三因素重复测量方差分析[Translation]

R统计绘图-混合方差分析[Translation]

R统计绘图-协方差分析[Translation]

R统计绘图-One-Way MANOVA

文献管理 | Endnote使用小技巧(修改引文格式、合并文献库、删除重复文献和更新文献信息等)

Alignment--本地blast使用详解1-数据库序列检索下载及比对