R语言之mvpart多元回归分析

时间尺度下，比较不同处理下微生物群落抗干扰能力

示例数据
提取码：zyxp

工具包的安装

mvpart属于比较老的工具包了，2014年之后再无更新，CRAN上已不支持，费了点周折。下面的解决方法很赞。看这里

#这里会下载一个100M大小的软件，略耗时，稍安勿躁。
install.packages("devtools")
devtools::install_github("cran/mvpart")
devtools::install_github("cran/MVPARTwrap")

数据准备

library(mvpart)
otu = read.table('L6.txt', header = T, row.names = 1, sep = '\t')   #物种数据，行为物种，列为样点
map = read.table('map.txt', header = T, row.names = 1, sep = '\t')  #实验设计，行为样点，列为分组信息

library('dplyr')
otu = otu[,rownames(map)] %>% .[rowSums(.) != 0,] %>% t() # 数据过滤

多元回归分析。

详细参数介绍，查看帮助??mvpart()

fit <- mvpart(as.matrix(otu) ~ map$Week, map, #物种变化对时间进行回归
              xv="pick", xval=nrow(otu), xvmult=500, #交叉验证系列参数
              legend=TRUE, margin=0.01, cp=0，which=4, big.pts=T, bars=F) #定义出图效果

自定义分组数量
多元回归树

对多元回归结果进行主成分分析，呈现为2维点阵。

详细参数介绍，查看帮助??rpart.pca()

rpart.pca(fit)
rpart.pca(fit,speclabs = F) #不显示物种标签

默认出图效果

输出默认出图

# pdf
pdf(file="rpart.pca.pdf", width=12, height=8)
rpart.pca(fit,speclabs = F)
dev.off()

# jpg
jpeg(file="rpart.pca.jpg", width=4000, height=3000, pointsize=8,units = 'px',res = 600)
rpart.pca(fit,speclabs = F,cex.lab = 4)
dev.off()

---

默认出图模式，操作起来不太友好，很多属性无法自定义设置，故将其中的坐标信息提取出来，使用ggplot2进行绘图

获取pca结果

rpart.pca(fit,speclabs = F) -> pts

# 获取物种的坐标
spec.ord = pts$y %>% data.frame()
colnames(spec.ord) = paste0('PC',seq(4))
rownames(spec.ord) = rownames(otu)

# 获取样点坐标
samp.ord = pts$xx
colnames(samp.ord) = paste0('PC',seq(4))
samp.ord = samp.ord %>% 
  cbind(subset(map,select = c(Week)))

# 求组内样点均值坐标
samp.mean = aggregate(.~ Week, samp.ord, mean)

绘图

library('ggplot2')
library('ggalt') #用于将各组样点分别连线

p <- ggplot(samp.ord,aes(PC1,PC2,color = as.factor(Week))) + 
  theme_minimal() +
  geom_point(size = 3) + # 样点
  geom_point(data = samp.mean,size = 5) + # 组内样点均值
  geom_encircle(s_shape = 1, expand = 0, lwd =1, fill = NA, show.legend = FALSE) + # 分组轮廓线
  geom_segment(data = spec.ord, aes(x=0, y=0, xend=PC1, yend=PC2),color = 'grey') + # 物种
  labs(x = 'Dim1 58.29%', y = 'Dim2 23.93%', color = 'Week') + #坐标轴和图例标题
  theme(axis.title = element_text(size = 15,color = 'black'), 
        axis.text = element_blank(),
        panel.grid = element_blank())

ggplot2出图效果

注意

-1. 使用ggplot2和原图有些不同，各组均值之间的联系不太清楚怎么连接的，但貌似也没啥用。
-3. ggplot2可以比较自由的设置图的属性，原图就比较受限制
-3. 所以，用原图还是ggplot2,随意