参考的文章为：《The host phylogeny determines viral infectivity and replication across Staphylococcus host species》

首先简单介绍一下 MCMCglmm 这个包的用法：

1.基本用法

library(dplyr)
data(Traffic, package = 'MASS')
Traffic$year <- as.factor(Traffic$year)
head(Traffic)

library(MCMCglmm)

# 定义先验分布
prior <- list(R=list(V=1, nu=0.002))

# 本例只有固定效应
m2a.5 <- MCMCglmm(y ~ limit + year + day, 
                  family = 'poisson',
                  data = Traffic, 
                  prior = prior, 
                  nitt=5010, thin=100, burnin=100,
                  pl = TRUE)

# 查看区间采样的参数值
m2a.5$Sol

其中

1. nitt 代表 MCMC 的迭代次数，采样次数
2. thin 代表用于计算的数据区间，thin = 100，代表使用100个数据计算一次待估计参数，thin = 1000，代表使用1000个数据计算一次待估计参数
3. burnin 代表丢弃前，后 burnin 次的采样数据（据说这样做收敛更准确？）；e.g：nitt = 1100，burnin = 100，则丢弃前100次和后100次取样，实际采样为101-1000次
   因此当 nitt = 10100，thin = 100，burnin = 1000 时有效采样数据为1001-99100次，并且计算的参数（截距，固定效应斜率和随机效应斜率的参数为一个组合）组合有 981 种

prior0 <- list(R = list(V = 1, nu = 0.002))
prior1 <- list(R = list(V = 1, nu = 0.002), B = list(mu = 0, V = 1e+08))
prior2 <- list(R = list(V = 1, nu = 0.002), 
               G = list(G1 = list(V = 1, nu = 0.002)))
prior3 <- list(R = list(V = 1, nu = 0.002), 
               G = list(G1 = list(V = 1,nu = 0.002), 
                           G2 = list(V = 1, nu = 0.002)))

1. prior 0 只包含了残差的先验概率分布（默认固定效应服从正态分布）
2. prior 1 只包含残差和固定效应的先验概率分布
3. prior 2 包含随机效应的先验概率分布（只有一个随机效应）
4. prior 3 包含随机效应的先验概率分布（两个随机效应）

一般认为随机效应的先验分布服从 inverse-wishart 分布，有两个参数 v 和 nu；固定效应一般服从正态分布

2. 固定效应和随机效应的写法

data
##         tarsus       back  animal     dam fosternest  hatchdate  sex
## 1: -1.89229718  1.1464212 R187142 R187557      F2102 -0.6874021  Fem
## 2:  1.13610981 -0.7596521 R187154 R187559      F1902 -0.6874021 Male
## 3:  0.98468946  0.1449373 R187341 R187568       A602 -0.4279814 Male
## 4:  0.37900806  0.2555847 R046169 R187518      A1302 -1.4656641 Male
## 5: -0.07525299 -0.3006992 R046161 R187528      A2602 -1.4656641  Fem
## 6: -1.13519543  1.5577219 R187409 R187945      C2302  0.3502805  Fem

prior = list(R = list(V = 1, nu = 0.002), 
                 G = list(G1 = list(V = 1,nu = 0.002), 
                          G2 = list(V = 1, nu = 0.002))
             )

MCMCglmm(tarsus ~ sex, 
         random = ~ dam + fosternest,
         prior = prior,
         pr=TRUE,
         data = data, verbose = FALSE)

1. random=∼dam+fosternest，可以改写了 random = ~us(1):dam + us(1):fosternest，random = ~idh(1):dam + idh(1):fosternest，1 代表有独立的截距项，即 dam 和 fosternest 单独为分类变量，dam 单独的分类变量为(R187142，R187341，R187409 ......)，fosternest (R187557，R187559，R187528 ......)，随机变量的分组就分别按照 dam 和 fosternest 进行广义线性模型分析
2. 如果是 us(sex):dam 或者 idh(sex):dam 代表 sex 与 dam 组合的分类变量，sex 相当于随机因子，sex 和 dam (Fem-R187142，Male-R187154，Male-R187518，Fem-R046161)，随机变量的分组就分别按照组合变量（sex 和 dam，Fem-R187142，Male-R187154，Male-R187518，Fem-R046161这样的分类） 进行广义线性模型分析
3. idh() 与 us() 的区别在于， idh() 只需要对角线值，即各个水平之间不存在交互作用（即协方差矩阵仅有对角线元素，其余元素为 0 ，而 us() ，则允许有对角线值和其余元素的值，即各个水平之间存在交互作用（协方差）
4. 例如 condition = A，B，C；则 idh(condition) 代表（condition 为随机因子）：
   idh(condition)
   σ²_A，σ²_B，σ²_C，代表自身与自身的协方差（某种角度就是相关性）
5. 例如 condition = A，B，C；则 us(condition) 代表（condition 为随机因子）：
   us(condition)
   σ²_A，σ²_B，σ²_C，代表自身与自身的协方差（某种角度就是相关性）；而 σ²_A,B，σ²_A,C，σ²_B,C 分别代表 A 和 B，A 和 C，B 和 C 的协方差（某种角度就是相关性），也表征 A 和 B，A 和 C，B 和 C 这些因素互作的强度

再举一个例子，

data
##    item1 item2
## 1:  a      A
## 2:  a      B
## 3:  b      A
## 4:  b      B
## 5:  b      A
## 6:  a      B

详细的随机效应的设计矩阵用法可见：如何定义两物种之间基因表达量的保守性

那么如果是单独的 item1 和 item2 水平用于随机效应，那么 random=∼dam+fosternest 的随机效应的设计矩阵为，假设 item1 对照为 a，item2 的对照为 A：

若写法为 us(sex):dam 或者 idh(sex):dam，将条件设为 item，item 对照为 a-A，sex 相当于随机因子，则它的设计矩阵为：

3. 固定效应加随机效应的矩阵表示方式

混合线性模型

X代表固定效应设计矩阵；Z 代表随机效应设计矩阵；b 为固定效应系数矩阵，u 代表随机效应系数矩阵

b 的求解方式为：

a 的求解方式：

G 代表随机因子协方差矩阵：例如 condition = A，B，C；利用标签为 A，B，C 的数据计算 A与A，A与B，A与C，B与B，B与C，C与C 之间的协方差矩阵

V 矩阵为：

推导为：

由于 a 和 e 的协方差为 0，因此后两项为 0；所以用户仅需要定义随机因子（condition，us(condition) 或 idh(condition)）的协方差矩阵 G 即可（us(condition) 考虑各个水平之间存在交互作用；idh(condition) 不考虑各个水平之间存在交互作用

矩阵表示：

文章例子

library(MCMCglmm)
library(ape)
library(ggtree)

# 读取分析数据
Data <- read.csv("PhageSusceptibilityData.csv", header = TRUE)

# 定义 trait 即 "OD", "qPCR", "ST_Cont", "ST_Bin"
Data$trait <- as.factor(Data$trait)
Data$trait <- factor(Data$trait, levels = c("OD", "qPCR", "ST_Cont", "ST_Bin"))

# 加载树文件
Tree <- read.tree("Staph_Phylogeny.nwk")

# 计算树的平均枝长
mean(diag(vcv(Tree)))

# 定义先验分布，考虑 Species
Prior.With.Species <- list(G = list(
  G1 = list(V = diag(4), nu = 4, alpha.mu = rep(0,4), alpha.V = diag(4) * 1000),
  G2 = list(V = diag(4), nu = 4, alpha.mu = rep(0,4), alpha.V = diag(4) * 1000)),
  R = list(V = diag(4), nu = 0.002, fix = 4))

# 定义先验分布，不考虑 Species
Prior.Without.Species <- list(G = list(
  G1 = list(V = diag(4), nu = 4, alpha.mu = rep(0,4), alpha.V = diag(4) * 1000)),
  R = list(V = diag(4), nu = 0.002, fix = 4))

# 用 inverseA 计算 phylogenetic variance-covariance 矩阵
# 基于树文件进行计算
Ainv<-inverseA(Tree, scale=FALSE)$Ainv

# MCMC
## 随机效应考虑了 trait 分别与 animal 和 Sample_Full 的组合
WithSpeciesModel <- 
  MCMCglmm(ISP ~ trait, random = ~us(trait) : animal + us(trait) : Sample_Full,
           rcov = ~idh(trait) : units, ginv = list(animal=Ainv), prior = Prior.With.Species, data = Data,
           family = NULL, nitt = 13000000, thin = 5000, burnin = 3000000, pr = TRUE)

## 随机效应考虑了 trait 和 animal 的组合
WithoutSpeciesModel <- 
  MCMCglmm(ISP ~ trait, random = ~us(trait) : animal,
           rcov = ~idh(trait) : units, ginv = list(animal=Ainv), prior = Prior.Without.Species, data = Data,
           family = NULL, nitt = 13000000, thin = 5000, burnin = 3000000, pr = TRUE)

1. 其中，~us(trait) : animal + us(trait) : Sample_Full 代表随机效应有两项，一项是 us(trait) : animal考虑了trait各个水平之间的互作（同时满足trait和animal两个标签的数据方可加入计算），另一项是 us(trait) : Sample_Full考虑了trait各个水平之间的互作（同时满足trait和Sample_Full两个标签的数据方可加入计算）。
2. rcov = ~idh(trait) : units，残差矩阵 e，来说明trait的水平可以有不同的残差大小，即残差矩阵与trait的水平有关。一般来说，idh(trait) 表征残差矩阵 e 只有主对角线有值，其余元素为 0（对角线数值表征trait各个水平之间的SE值）；us(trait) 表征残差矩阵 e 只有主对角线有值，其余元素也有数值
3. 补充 2，例如 trait = a，b，c；idh(trait) 所得的残差矩阵 e 对角线元素表示为标签分别为 a，b，c 的数据的SE值：
   其中，SE_a 代表标签为 a 的数据所计算的SE值
4. ginv = list(animal=Ainv) 推测是输入的随机效应 animal 的协方差矩阵，因为 Ainv 矩阵表示的是两个序列之间的遗传距离
   Ainv 矩阵
   
   Data 如下：
   Data
   因此矩阵 Ainv 的行列名需要在 Data animal 列中存在，那么随机效应计算协方差时满足：

参考：https://nhcooper123.github.io/pcm-primer-online/phylogenetic-generalised-linear-mixed-models-in-r.html