知識背景
廣義線性混合模型可以看做是廣義線性模型（GLM）以及線性混合模型(LMM)的擴展琴庵，為了更好地理解GLMM迷殿，肯定要對普通線性模型、廣義線性模型以及線性混合模型有個理解蚊夫。
普通線性模型長大概這個樣子（以只有一個因變量舉例）：
或
也就是固定效應+隨機誤差知纷，且因變量必須要滿足正態(tài)性导披、獨立性以及方差齊性撩匕。廣義線性模型（GLM）可以長這個樣子：
也可以長這個樣子：...
也就是對因變量y進行了變換，使得變換后的值適用普通線性回歸模蜡。而GLM之所以稱為“廣義”忍疾，關鍵就在于這種變換——連接函數(shù)（就是上邊兩個等式左邊那一堆）谨朝，使得模型的應用范圍變得更寬，只要因變量服從指數(shù)分布就行则披；而線性混合模型（LMM）則形如（以矩陣形式表示）：
相較普通線性模型士复，LMM其實就是多了隨機效應而已。
如何劃分固定效應和隨機效應便贵？
（1）固定效應的因子水平明確承璃，在抽樣數(shù)據(jù)集中包含其所有類別俏竞，通常是在實驗中加以控制的因素魂毁，例如組別。
（2）隨機效應的因子水平不清晰咬崔，在抽樣數(shù)據(jù)集中并不包含該自變量的所有情況烦秩，通常是實驗中無法控制的因素只祠，例如個體。
換句話說熊杨，以固定方式對因變量產(chǎn)生影響的因素為固定效應盗舰，相反钻趋，則為隨機效應。
- 廣義線性混合模型
至此较沪，GLMM應該就好理解了尸曼，其實就是在等式左側(cè)對因變量進行變換，在等式右側(cè)將固定效應與隨機效應進行結合。但公式寫起來蠻復雜解幽，我就不列出來了烘苹。
示例
光說不練假把式镣衡，下面我們以縱向研究的宏基因組數(shù)據(jù)（分實驗組和對照組，且每一個個體在3個時間點采樣）為例來說明GLMM的具體應用：
假設你已經(jīng)獲得了物種豐度結果taxa.raw望浩，形如：
以及樣本的各種信息sample.info磨德，形如：
- 步驟1：多濾掉低豐度數(shù)據(jù)
#過濾掉物種豐度為0的樣本數(shù)過多的行
> filter.index <- apply(taxa.raw, 1, function(X){sum(X>0)>0.1*length(X)})
> taxa.filter <- taxa.raw[filter.index, ]
#確保樣本豐度和相加為100
> taxa.filter <- 100*sweep(taxa.filter, 2, colSums(taxa.filter), FUN="/")
> taxa.data <- as.data.frame(t(taxa.filter))
- 步驟2：創(chuàng)建協(xié)變量矩陣
我們可以用如下代碼為group典挑、time以及group和time的交互項創(chuàng)建協(xié)變量矩陣：
> reg.cov <- tibble::rownames_to_column(sample.info, var = 'Sample') %>%
dplyr::select(Sample, group, indID, time) %>%
dplyr::mutate(time = ifelse(time == 'w0', 0, ifelse(time == 'w2', 1, ifelse(time == 'w4', 2,NA)))) %>% #創(chuàng)建時間變量代碼：w0是0您觉，w2是1琳水，w4是2褒墨；
dplyr::mutate(group = ifelse(group == 'ck', 0, 1)) %>% #創(chuàng)建分組變量代碼：ck組為0郁妈，A組為1
dplyr::mutate(timeXgroup = time*group) %>% #創(chuàng)建時間與分組交互代碼
dplyr::mutate(indID = as.character(indID)) %>%
as.data.frame
#過濾掉采樣不全的個體數(shù)據(jù)
> ind.count <- table(reg.cov$indID)
> reg.cov <- subset(reg.cov, indID %in% names(ind.count)[ind.count == 3])
> head(reg.cov)
Sample group indID time timeXgroup
1 FN3T 1 A2-1 1 1
2 FN40 0 D2-4 0 0
3 FN42 1 A1-4 1 1
4 FN44 1 A1-5 0 0
5 FN45 1 A2-3 1 1
7 FN49 0 D2-2 1 0
- 步驟3：擬合模型
先將w0和w1+w2的數(shù)據(jù)分別取出備用:
> reg.cov.w0 <- subset(reg.cov, time == 0)
> rownames(reg.cov.w0) <- reg.cov.w0$indID
> head(reg.cov.w0)
Sample group indID time timeXgroup
D2-4 FN40 0 D2-4 0 0
A1-5 FN44 1 A1-5 0 0
D2-5 FN4D 0 D2-5 0 0
D1-2 FN54 0 D1-2 0 0
A1-4 FN5J 1 A1-4 0 0
A2-4 FN60 1 A2-4 0 0
> reg.cov.w12 <- subset(reg.cov, time != 0)
> reg.cov.w12 <- na.omit(data.frame(baseline.sample = reg.cov.w0[reg.cov.w12$indID, 'Sample'],
baseline.subject = reg.cov.w0[reg.cov.w12$indID, 'indID'],
reg.cov.w12,
stringsAsFactors = F))
> head(reg.cov.w12)
baseline.sample baseline.subject Sample group indID time timeXgroup
1 FN70 A2-1 FN3T 1 A2-1 1 1
3 FN5J A1-4 FN42 1 A1-4 1 1
5 FN68 A2-3 FN45 1 A2-3 1 1
7 FN71 D2-2 FN49 0 D2-2 1 0
8 FN54 D1-2 FN4A 0 D1-2 2 0
9 FN6T D1-5 FN4B 0 D1-5 1 0
然后單獨為每一個物種進行擬合（為了便于說明顾彰，此處我們僅以其中一個物種豐度舉例：
library(dplyr)
library(nlme)
> taxa.all <- colnames(taxa.data)
> p.taxa.list.lme <- list()
> spe <- taxa.all[9]
> spe
[1] "s__Bacteroides_vulgatus"
> X <- data.frame(Baseline = taxa.data[reg.cov.w12$baseline.sample, spe]/100,
reg.cov.w12[, c('time', 'group')])
> rownames(X) <- reg.cov.w12$Sample
> head(X)
Baseline time group
FN3T 0.0004122127 1 1
FN42 0.0001939168 1 1
FN45 0.0026208886 1 1
FN49 0.0003091166 1 0
FN4A 0.0007188949 2 0
FN4B 0.0005986076 1 0
> Z <- X
> subject.ind <- reg.cov.w12$indID
> time.ind <- reg.cov.w12$time
> Y <- taxa.data[reg.cov.w12$Sample, spe]/100
> cat('Zeros/All',sum(Y==0),'/',length(Y),'\n')
Zeros/All 0 / 30
#對Y進行變換（先開方在反正弦涨享？）
> tdata <- data.frame(Y.tran = asin(sqrt(Y)), X, SID = subject.ind)
> head(tdata)
Y.tran Baseline time group SID
FN3T 0.04531814 0.0004122127 1 1 A2-1
FN42 0.05269246 0.0001939168 1 1 A1-4
FN45 0.01527299 0.0026208886 1 1 A2-3
FN49 0.02599619 0.0003091166 1 0 D2-2
FN4A 0.01817293 0.0007188949 2 0 D1-2
FN4B 0.03147346 0.0005986076 1 0 D1-5
#進行擬合（以Baseline厕隧、time以及group為固定效應，以SID為隨機效應）
> lme.fit <- try(lme(Y.tran ~ Baseline + time + group, random = ~1|SID, data = tdata))
隨機效應的表達方式：
（1）1|SID：SID是隨機截距髓迎，Baseline建丧、time翎朱、group是固定斜率，也就是他們前面的參數(shù)是固定的拴曲；
（2）0 + time|SID：SID是固定截距争舞，time是隨機斜率，Baseline和time是固定斜率;
（3）1 + time|SID：SID是隨機截距疗韵，time是隨機斜率兑障，Baseline和time是固定斜率;
> summary(lme.fit)
...
Random effects:
Formula: ~1 | SID
(Intercept) Residual
StdDev: 0.007468231 0.009009141
Fixed effects: Y.tran ~ Baseline + time + group
Value Std.Error DF t-value p-value
(Intercept) 0.0236327 0.0060813 14 3.886108 0.0016
Baseline -1.2738533 1.3537572 12 -0.940976 0.3653
time 0.0023881 0.0032897 14 0.725932 0.4798
group 0.0016992 0.0053765 12 0.316037 0.7574
...
這樣我們就得到了每個自變量的參數(shù)，以及P值蕉汪，后續(xù)我們再使用多重檢驗矯正就可以得到矯正后的P值了流译。
