測(cè)量方差分析：兩組縱向數(shù)據(jù)實(shí)戰(zhàn)指南)
重復(fù)測(cè)量數(shù)據(jù)在臨床研究和實(shí)驗(yàn)科學(xué)里出現(xiàn)的頻率極高——同一個(gè)受試者在不同時(shí)間點(diǎn)被反復(fù)觀測(cè)比如術(shù)后第1天、第3天、第7天的某項(xiàng)指標(biāo)或者用藥前后多個(gè)隨訪節(jié)點(diǎn)的量表得分。這類數(shù)據(jù)的特點(diǎn)是觀測(cè)值之間不獨(dú)立同一個(gè)體內(nèi)部的多次測(cè)量存在相關(guān)性如果直接拿去做普通的獨(dú)立樣本t檢驗(yàn)或者單因素方差分析p值大概率是錯(cuò)的結(jié)論也就站不住腳。我見過不少做臨床數(shù)據(jù)分析的朋友數(shù)據(jù)明明是重復(fù)測(cè)量的結(jié)構(gòu)卻硬生生把每個(gè)時(shí)間點(diǎn)拆開做多次t檢驗(yàn)結(jié)果要么假陽(yáng)性泛濫要么把本來顯著的處理效應(yīng)給淹沒了。這篇內(nèi)容就圍繞“基于R進(jìn)行兩組患者重復(fù)測(cè)量數(shù)據(jù)的分析”這個(gè)主題把重復(fù)測(cè)量方差分析的完整思路、R語言落地細(xì)節(jié)、以及我踩過的坑一次性講透。適合有一定R基礎(chǔ)、正在處理縱向隨訪數(shù)據(jù)或者交叉實(shí)驗(yàn)數(shù)據(jù)的從業(yè)者參考也適合剛接觸重復(fù)測(cè)量設(shè)計(jì)、想搞清楚“為什么不能拆開做t檢驗(yàn)”的研究人員從頭看起。1. 重復(fù)測(cè)量方差分析到底在解決什么問題1.1 從數(shù)據(jù)結(jié)構(gòu)說起為什么普通方差分析不夠用先把數(shù)據(jù)結(jié)構(gòu)擺清楚。假設(shè)我們做了一個(gè)兩組對(duì)照的臨床觀察A組用新方案B組用常規(guī)方案每組各若干名患者在基線、第4周、第8周、第12周四個(gè)時(shí)間點(diǎn)測(cè)量同一個(gè)連續(xù)型指標(biāo)。數(shù)據(jù)排成“長(zhǎng)格式”之后每一行是“某患者在某時(shí)間點(diǎn)的一次測(cè)量”包含患者ID、分組、時(shí)間點(diǎn)、測(cè)量值四個(gè)核心字段。這種結(jié)構(gòu)下同一個(gè)患者貢獻(xiàn)了4行數(shù)據(jù)這4個(gè)觀測(cè)值天然相關(guān)——一個(gè)基線值偏高的患者后續(xù)幾個(gè)時(shí)間點(diǎn)大概率也偏高。普通單因素方差分析的前提是觀測(cè)獨(dú)立這個(gè)前提在這里直接被打破了。如果無視相關(guān)性把4個(gè)時(shí)間點(diǎn)當(dāng)成4個(gè)獨(dú)立組來做方差分析誤差項(xiàng)里混入了“患者個(gè)體差異”這部分變異導(dǎo)致誤差均方被高估檢驗(yàn)效能下降。更麻煩的是時(shí)間點(diǎn)之間的比較會(huì)失真因?yàn)槟銢]法區(qū)分“時(shí)間效應(yīng)”和“個(gè)體差異”。重復(fù)測(cè)量方差分析的核心思路就是把總變異拆成三塊組間變異處理效應(yīng)、組內(nèi)變異時(shí)間效應(yīng)、以及交互效應(yīng)處理×?xí)r間同時(shí)把“個(gè)體差異”從誤差項(xiàng)里單獨(dú)剝離出來讓誤差估計(jì)更干凈。這就是它比普通方差分析更適合縱向數(shù)據(jù)的根本原因。用生活化的類比假設(shè)你想比較兩個(gè)班級(jí)學(xué)生的考試成績(jī)變化如果每個(gè)學(xué)生只考一次那是獨(dú)立樣本問題但如果是同一批學(xué)生期中、期末各考一次學(xué)生自身的水平差異就是一個(gè)巨大的干擾項(xiàng)。重復(fù)測(cè)量方差分析相當(dāng)于先把“這個(gè)學(xué)生本來就行/本來就不行”這部分變異扣掉再去看“時(shí)間帶來的變化”和“班級(jí)之間的差異”結(jié)論自然更靠譜。1.2 兩組設(shè)計(jì)的特殊性組間、組內(nèi)、交互三重效應(yīng)兩組重復(fù)測(cè)量設(shè)計(jì)里我們關(guān)心的效應(yīng)有三個(gè)層面。第一是組間效應(yīng)也就是A組和B組整體上有沒有差異這個(gè)效應(yīng)在統(tǒng)計(jì)上等價(jià)于“把每個(gè)患者的多次測(cè)量取平均再做兩組比較”但它用的是合并后的誤差項(xiàng)比直接做t檢驗(yàn)更穩(wěn)健。第二是組內(nèi)效應(yīng)也就是時(shí)間主效應(yīng)所有患者合在一起看指標(biāo)隨時(shí)間有沒有顯著變化。第三是交互效應(yīng)這是最有臨床價(jià)值的一項(xiàng)——兩組隨時(shí)間變化的趨勢(shì)是否不同。舉個(gè)例子如果新方案組在第8周之后指標(biāo)快速下降而常規(guī)組下降平緩交互效應(yīng)就會(huì)顯著說明處理方案的效果依賴于時(shí)間這比單純看“兩組有沒有差異”信息量大得多。很多初學(xué)者只盯著組間p值看忽略了交互效應(yīng)這是很可惜的。臨床上真正有意義的往往是“處理組起效更快”或者“處理組在后期維持得更好”這類趨勢(shì)差異而這些恰恰體現(xiàn)在交互項(xiàng)上。所以在R里跑重復(fù)測(cè)量方差分析時(shí)交互項(xiàng)的顯著性一定要單獨(dú)拎出來解讀不能只看主效應(yīng)。1.3 球形假設(shè)這個(gè)前提不滿足會(huì)怎樣重復(fù)測(cè)量方差分析有一個(gè)繞不開的假設(shè)——球形假設(shè)sphericity。簡(jiǎn)單說它要求任意兩個(gè)時(shí)間點(diǎn)之間差值的方差相等。四個(gè)時(shí)間點(diǎn)的話就有6對(duì)差值這6個(gè)差值的方差理論上應(yīng)該差不多。如果球形假設(shè)被違反組內(nèi)效應(yīng)和交互效應(yīng)的F檢驗(yàn)自由度需要校正否則p值偏小容易假陽(yáng)性。球形假設(shè)在時(shí)間點(diǎn)較多、或者時(shí)間間隔不均勻時(shí)特別容易出問題。我的經(jīng)驗(yàn)是三個(gè)時(shí)間點(diǎn)以內(nèi)球形假設(shè)通常還過得去四個(gè)及以上時(shí)間點(diǎn)最好默認(rèn)做校正。R里car包的Anova()函數(shù)配合idata參數(shù)可以自動(dòng)輸出GG校正Greenhouse-Geisser和HF校正Huynh-Feldt的結(jié)果直接看校正后的p值就行。這里有個(gè)細(xì)節(jié)球形假設(shè)只影響組內(nèi)效應(yīng)和交互效應(yīng)組間效應(yīng)不受影響因?yàn)榻M間比較用的是合并的個(gè)體均值不涉及時(shí)間點(diǎn)之間的差值結(jié)構(gòu)。這一點(diǎn)很多人搞混以為球形假設(shè)不滿足整個(gè)分析就廢了其實(shí)不是。2. R語言落地的整體方案與工具選型2.1 寬格式還是長(zhǎng)格式數(shù)據(jù)整形的關(guān)鍵決策在R里做重復(fù)測(cè)量方差分析數(shù)據(jù)格式的選擇直接影響后續(xù)代碼怎么寫。寬格式是每個(gè)患者一行不同時(shí)間點(diǎn)的測(cè)量值各占一列比如id, group, t1, t2, t3, t4。長(zhǎng)格式是每個(gè)患者每個(gè)時(shí)間點(diǎn)一行比如id, group, time, value。兩種格式各有用途寬格式適合用aov()配合Error()項(xiàng)做經(jīng)典的重復(fù)測(cè)量方差分析長(zhǎng)格式適合用混合效應(yīng)模型或者car::Anova()配合idata參數(shù)。我個(gè)人的建議是如果時(shí)間點(diǎn)固定且不多比如四個(gè)隨訪節(jié)點(diǎn)用寬格式配合aov()最直接代碼短、輸出清晰如果時(shí)間點(diǎn)較多、或者有缺失值、或者想順便做趨勢(shì)分析長(zhǎng)格式更靈活。實(shí)際項(xiàng)目里我通常兩份都準(zhǔn)備——用tidyr::pivot_longer()和pivot_wider()互相轉(zhuǎn)換幾行代碼的事。需要注意的是寬格式轉(zhuǎn)長(zhǎng)格式時(shí)列名要規(guī)范比如week0, week4, week8, week12轉(zhuǎn)完之后時(shí)間因子要顯式設(shè)定水平順序否則R默認(rèn)按字母排序week12會(huì)排在week4前面圖和分析全亂套。2.2 核心函數(shù)對(duì)比aov、Anova與混合模型R里能做重復(fù)測(cè)量方差分析的函數(shù)不止一個(gè)選哪個(gè)取決于你的設(shè)計(jì)復(fù)雜度和對(duì)缺失值的容忍度。stats::aov()配合Error()項(xiàng)是最經(jīng)典的做法語法是aov(value ~ group * time Error(id/time))其中Error(id/time)告訴R“個(gè)體是隨機(jī)效應(yīng)時(shí)間是個(gè)體內(nèi)因子”。這個(gè)函數(shù)的好處是輸出直接給出組間、組內(nèi)、交互三張方差分析表球形假設(shè)的校正需要額外用summary()配合Mauchly檢驗(yàn)來看。car::Anova()是另一個(gè)常用選擇它需要把寬格式數(shù)據(jù)轉(zhuǎn)成mlm對(duì)象多元線性模型然后用idata指定個(gè)體內(nèi)因子結(jié)構(gòu)。它的優(yōu)勢(shì)是自動(dòng)輸出球形校正結(jié)果而且對(duì)不平衡設(shè)計(jì)更穩(wěn)健。缺點(diǎn)是數(shù)據(jù)準(zhǔn)備稍麻煩需要構(gòu)造idata數(shù)據(jù)框。如果數(shù)據(jù)有缺失值臨床隨訪里太常見了上面兩個(gè)函數(shù)都會(huì)直接丟掉含缺失的整行樣本量縮水嚴(yán)重。這時(shí)候混合效應(yīng)模型lme4::lmer()或nlme::lme()是更優(yōu)解它用限制性極大似然估計(jì)能利用部分信息不需要完整案例。不過混合模型的方差結(jié)構(gòu)設(shè)定更復(fù)雜隨機(jī)斜率、隨機(jī)截距的選擇需要結(jié)合設(shè)計(jì)來定不是一句兩句能說完的。這篇內(nèi)容聚焦經(jīng)典重復(fù)測(cè)量方差分析混合模型留到后續(xù)單獨(dú)展開。2.3 必備包清單與版本注意事項(xiàng)跑完整套分析我常用的包有這么幾個(gè)car負(fù)責(zé)Anova()和球形校正emmeans負(fù)責(zé)事后兩兩比較和簡(jiǎn)單效應(yīng)分析tidyr和dplyr負(fù)責(zé)數(shù)據(jù)整形ggplot2負(fù)責(zé)可視化rstatix提供了一些便捷的管道式統(tǒng)計(jì)函數(shù)。版本上要注意car包3.0之后的Anova()對(duì)idata參數(shù)的校驗(yàn)更嚴(yán)格idata的行數(shù)必須和個(gè)體內(nèi)因子的水平數(shù)完全一致順序也要對(duì)上否則報(bào)錯(cuò)信息很隱晦。emmeans在1.7版本之后對(duì)重復(fù)測(cè)量模型的mode參數(shù)有調(diào)整做簡(jiǎn)單效應(yīng)時(shí)記得顯式指定mode df.error或者根據(jù)模型類型選對(duì)自由度近似方式不然置信區(qū)間會(huì)偏窄。提示裝包時(shí)建議鎖定版本尤其是car和emmeans這兩個(gè)包更新頻繁不同版本之間輸出格式和默認(rèn)行為有差異團(tuán)隊(duì)協(xié)作時(shí)最好在項(xiàng)目里記錄sessionInfo()。3. 完整實(shí)操流程從數(shù)據(jù)準(zhǔn)備到結(jié)果解讀3.1 模擬一份兩組四時(shí)點(diǎn)的臨床數(shù)據(jù)為了把流程跑通我先構(gòu)造一份模擬數(shù)據(jù)。兩組各20名患者四個(gè)時(shí)間點(diǎn)指標(biāo)值設(shè)定為A組基線均值50隨時(shí)間逐步下降到第12周的38B組基線均值50下降到第12周的45。個(gè)體隨機(jī)效應(yīng)標(biāo)準(zhǔn)差設(shè)為5殘差標(biāo)準(zhǔn)差設(shè)為3。這樣設(shè)定下交互效應(yīng)應(yīng)該是顯著的因?yàn)閮山M下降斜率不同。set.seed(2024) n_per_group - 20 times - c(0, 4, 8, 12) groups - rep(c(A, B), each n_per_group) id - 1:(2 * n_per_group) # 個(gè)體隨機(jī)截距 rand_intercept - rnorm(2 * n_per_group, 0, 5) # 構(gòu)造長(zhǎng)格式 dat_long - expand.grid(id id, time times) dat_long$group - groups[match(dat_long$id, id)] dat_long$value - 50 rand_intercept[dat_long$id] ifelse(dat_long$group A, -1.0 * dat_long$time, -0.4 * dat_long$time) rnorm(nrow(dat_long), 0, 3) dat_long$time - factor(dat_long$time, levels times) dat_long$group - factor(dat_long$group, levels c(A, B))這段代碼里有兩個(gè)細(xì)節(jié)值得說。第一time轉(zhuǎn)因子時(shí)顯式指定了levels times保證順序是0、4、8、12而不是默認(rèn)的字符串排序。第二個(gè)體隨機(jī)截距用rand_intercept[dat_long$id]按id索引確保同一個(gè)體的四次測(cè)量共享同一個(gè)隨機(jī)截距這是模擬重復(fù)測(cè)量數(shù)據(jù)的關(guān)鍵漏了這一步數(shù)據(jù)就退化成獨(dú)立觀測(cè)了。轉(zhuǎn)成寬格式用于aov()分析library(tidyr) dat_wide - pivot_wider(dat_long, names_from time, values_from value, names_prefix t)轉(zhuǎn)完之后列名是t0, t4, t8, t12每個(gè)患者一行配合group列正好是aov()需要的結(jié)構(gòu)。3.2 用aov加Error項(xiàng)跑經(jīng)典重復(fù)測(cè)量方差分析寬格式數(shù)據(jù)準(zhǔn)備好之后核心代碼就一行l(wèi)ibrary(car) mod_aov - aov(cbind(t0, t4, t8, t12) ~ group, data dat_wide) summary(mod_aov)注意這里用的是cbind()把四個(gè)時(shí)間點(diǎn)綁成一個(gè)矩陣作為因變量group是組間因子。這個(gè)寫法下summary()輸出的是組間效應(yīng)的檢驗(yàn)結(jié)果。要拿到組內(nèi)效應(yīng)和交互效應(yīng)需要把模型轉(zhuǎn)成mlm對(duì)象再用Anova()mod_mlm - lm(cbind(t0, t4, t8, t12) ~ group, data dat_wide) idata - data.frame(time factor(c(t0, t4, t8, t12), levels c(t0, t4, t8, t12))) Anova(mod_mlm, idata idata, idesign ~time, type III)idesign ~time指定個(gè)體內(nèi)因子是timeidata的行數(shù)和順序必須和cbind()里列的排列完全一致這是最容易出錯(cuò)的地方。輸出里會(huì)包含time的球形檢驗(yàn)結(jié)果Mauchlys test以及未校正和GG/HF校正后的F值和p值。組間效應(yīng)看group那一行交互效應(yīng)看group:time那一行。實(shí)測(cè)這份模擬數(shù)據(jù)交互項(xiàng)的GG校正p值應(yīng)該在0.05以下說明兩組的時(shí)間趨勢(shì)確實(shí)不同。組間主效應(yīng)可能不顯著因?yàn)閮山M在第12周的均值差距被基線到中期的重疊稀釋了——這正好說明為什么不能只看組間p值。3.3 球形檢驗(yàn)與校正結(jié)果的取舍Anova()輸出里會(huì)有一行Mauchlys test for sphericity給出W統(tǒng)計(jì)量和p值。如果p 0.05說明球形假設(shè)被違反需要看校正后的結(jié)果。GG校正Greenhouse-Geisser的epsilon值通常偏保守HF校正Huynh-Feldt偏寬松。我的取舍原則是epsilon在0.75以上用HF校正0.75以下用GG校正。這個(gè)閾值不是絕對(duì)的但實(shí)踐中比較穩(wěn)妥。還有一個(gè)細(xì)節(jié)球形檢驗(yàn)本身對(duì)樣本量敏感小樣本時(shí)檢驗(yàn)效能低可能漏檢大樣本時(shí)又過于敏感輕微偏離就報(bào)顯著。所以不能完全依賴p值判斷結(jié)合epsilon值一起看更靠譜。如果epsilon接近1比如0.95以上即使Mauchly檢驗(yàn)顯著未校正結(jié)果也基本可用。注意球形假設(shè)只約束組內(nèi)效應(yīng)和交互效應(yīng)組間效應(yīng)的F檢驗(yàn)不需要校正。解讀結(jié)果時(shí)別把組間p值也拿去做校正那是畫蛇添足。3.4 事后比較簡(jiǎn)單效應(yīng)分析與兩兩對(duì)比交互效應(yīng)顯著之后下一步是拆解——在每個(gè)時(shí)間點(diǎn)上比較兩組或者在每組內(nèi)部比較各時(shí)間點(diǎn)。這就是簡(jiǎn)單效應(yīng)分析。用emmeans包可以很干凈地完成library(emmeans) emm - emmeans(mod_mlm, ~ group | time, mode df.error) pairs(emm)~ group | time的意思是“在每個(gè)time水平上比較group”輸出的是四個(gè)時(shí)間點(diǎn)各自的兩組對(duì)比。反過來如果想在每組內(nèi)部比較時(shí)間點(diǎn)寫成~ time | group。mode df.error指定用誤差項(xiàng)自由度做近似對(duì)于經(jīng)典的重復(fù)測(cè)量設(shè)計(jì)這是合適的如果用的是混合模型要改成mode satterthwaite或kenward-roger。多重比較校正別忘了。四個(gè)時(shí)間點(diǎn)各做一次兩組比較一共4次檢驗(yàn)用Bonferroni校正的話alpha除以4。emmeans里加adjust bonferroni或者h(yuǎn)olm都行我一般用Holm因?yàn)樗鼘?duì)檢驗(yàn)之間的相關(guān)性不敏感功效比Bonferroni略高。如果交互不顯著但時(shí)間主效應(yīng)顯著那就退一步只看所有患者合并后的時(shí)間點(diǎn)兩兩比較用emmeans(mod_mlm, ~ time)配合pairs()即可。這時(shí)候不需要按組分拆因?yàn)榻换ゲ伙@著意味著兩組的趨勢(shì)可以認(rèn)為是一致的。4. 常見問題與排查技巧實(shí)錄4.1 報(bào)錯(cuò)“number of items to replace is not a multiple of replacement length”這個(gè)報(bào)錯(cuò)在構(gòu)造idata或者模擬數(shù)據(jù)時(shí)特別常見根源是向量長(zhǎng)度不匹配。比如idata寫了4行但cbind()里只綁了3列或者模擬數(shù)據(jù)時(shí)rand_intercept的長(zhǎng)度和id對(duì)不上。排查方法很簡(jiǎn)單在報(bào)錯(cuò)的那一行之前用length()和nrow()把涉及的每個(gè)對(duì)象的維度打印出來逐個(gè)核對(duì)。我習(xí)慣在關(guān)鍵步驟后加一句stopifnot(length(x) expected)做斷言早發(fā)現(xiàn)早解決比事后debug省時(shí)間。另一個(gè)隱蔽的變體是因子水平順序不一致。idata里time的水平是t0, t4, t8, t12但cbind()里列的順序是t0, t8, t4, t12R不會(huì)報(bào)錯(cuò)但結(jié)果全錯(cuò)——因?yàn)樗次恢闷ヅ洳话疵制ヅ?。這種錯(cuò)誤最難查因?yàn)榇a能跑通輸出看起來也正常。我的做法是cbind()里的列順序和idata的行順序都用同一個(gè)變量控制比如先定義time_levels - c(t0,t4,t8,t12)兩邊都引用它從源頭上杜絕錯(cuò)位。4.2 缺失值導(dǎo)致樣本量驟減怎么辦臨床隨訪數(shù)據(jù)幾乎不可能沒有缺失患者失訪、漏檢、樣本溶血各種原因都會(huì)造成某個(gè)時(shí)間點(diǎn)缺值。aov()和Anova()默認(rèn)做完整案例分析只要一個(gè)患者四個(gè)時(shí)間點(diǎn)里缺了任何一個(gè)整行就被丟掉。20人的組缺了5個(gè)樣本量直接掉四分之一檢驗(yàn)效能大打折扣。應(yīng)對(duì)策略分兩層。第一層是分析前的缺失機(jī)制判斷如果缺失完全隨機(jī)MCAR完整案例分析雖然浪費(fèi)信息但結(jié)論無偏如果缺失與觀測(cè)值有關(guān)MAR完整案例分析就有偏了這時(shí)候應(yīng)該上混合效應(yīng)模型或者多重插補(bǔ)。第二層是實(shí)操層面的補(bǔ)救用mice包做多重插補(bǔ)插補(bǔ)后對(duì)每個(gè)插補(bǔ)數(shù)據(jù)集跑一遍重復(fù)測(cè)量方差分析再用Rubin規(guī)則合并結(jié)果。這套流程代碼量不小但mice和mitml配合起來還算順手。如果缺失比例很低比如5%以內(nèi)完整案例分析的影響可以接受但要在方法部分如實(shí)報(bào)告缺失情況別藏著。4.3 交互效應(yīng)顯著但事后比較全不顯著這種情況我遇到過好幾次第一次碰到時(shí)一頭霧水——交互p值0.03但四個(gè)時(shí)間點(diǎn)分別做兩組比較校正后一個(gè)都不顯著。原因通常是交互效應(yīng)體現(xiàn)在“趨勢(shì)”上而非“某個(gè)具體時(shí)間點(diǎn)的截距差異”上。比如兩組都是線性下降但斜率不同在中間某個(gè)時(shí)間點(diǎn)可能剛好交叉單點(diǎn)比較自然不顯著但整體趨勢(shì)差異是真實(shí)的。解決辦法有兩個(gè)方向。一是做趨勢(shì)分析把時(shí)間因子拆成線性、二次、三次正交多項(xiàng)式對(duì)比看交互效應(yīng)主要落在哪個(gè)趨勢(shì)成分上。emmeans里用contrast(emm, poly)或者手動(dòng)指定正交對(duì)比系數(shù)都行。二是做區(qū)域檢驗(yàn)或者模型-based inference比如擬合一個(gè)時(shí)間×組的線性回歸直接檢驗(yàn)斜率差異。這兩種方法都比逐點(diǎn)比較更貼合“趨勢(shì)差異”的本質(zhì)。我的經(jīng)驗(yàn)是如果研究假設(shè)本來就是“起效速度不同”那從一開始就該規(guī)劃趨勢(shì)對(duì)比而不是等到事后比較不顯著了再補(bǔ)救。4.4 問題速查表問題現(xiàn)象可能原因排查與解決Anova()報(bào)錯(cuò)維度不匹配idata行數(shù)與cbind()列數(shù)不一致用length()核對(duì)統(tǒng)一用變量控制順序結(jié)果與預(yù)期完全相反因子水平順序錯(cuò)位檢查factor(levels)設(shè)定idata與列順序?qū)R組內(nèi)效應(yīng)p值異常小球形假設(shè)違反未校正看Mauchly檢驗(yàn)和epsilon選GG或HF校正事后比較全不顯著交互體現(xiàn)在趨勢(shì)而非單點(diǎn)改用正交多項(xiàng)式對(duì)比或斜率檢驗(yàn)樣本量比預(yù)期少很多缺失值完整案例刪除評(píng)估缺失比例考慮多重插補(bǔ)或混合模型emmeans置信區(qū)間過窄自由度近似方式不對(duì)顯式指定mode參數(shù)匹配模型類型4.5 幾個(gè)容易被忽略的實(shí)操心得第一個(gè)心得先畫圖再跑模型。把每個(gè)患者的軌跡畫出來ggplot2里用geom_line(aes(group id))再疊加組均值軌跡數(shù)據(jù)的形態(tài)一目了然。有沒有離群個(gè)體、缺失集中在哪個(gè)時(shí)間點(diǎn)、兩組軌跡是平行還是交叉圖上一看便知。很多統(tǒng)計(jì)上的困惑畫個(gè)圖就解開了。我現(xiàn)在的習(xí)慣是任何重復(fù)測(cè)量分析圖一定跑在模型前面。第二個(gè)心得樣本量估算要在設(shè)計(jì)階段做。重復(fù)測(cè)量設(shè)計(jì)的樣本量估算和普通兩組比較不一樣要考慮時(shí)間點(diǎn)之間的相關(guān)性用預(yù)實(shí)驗(yàn)的相關(guān)系數(shù)估計(jì)和球形假設(shè)的epsilon。pwr包和WebPower包都能做但參數(shù)設(shè)定需要預(yù)實(shí)驗(yàn)數(shù)據(jù)支撐。事后補(bǔ)算意義不大但方法部分最好報(bào)告一下估算依據(jù)。第三個(gè)心得報(bào)告結(jié)果時(shí)把效應(yīng)量和置信區(qū)間帶上。光報(bào)p值信息量不夠尤其是交互效應(yīng)報(bào)個(gè)偏eta方或者兩組均值差的置信區(qū)間讀者才能判斷臨床意義。emmeans的confint()直接給置信區(qū)間effectsize包可以算偏eta方幾行代碼的事別省。第四個(gè)心得時(shí)間因子的編碼方式影響解讀。如果時(shí)間點(diǎn)間隔不均勻比如0、4、8、12周是均勻的但0、2、6、12就不均勻把時(shí)間當(dāng)分類因子做方差分析會(huì)損失“間隔信息”。這時(shí)候可以考慮把時(shí)間當(dāng)連續(xù)變量做趨勢(shì)分析或者用正交多項(xiàng)式對(duì)比。分類因子適合“每個(gè)時(shí)間點(diǎn)都關(guān)心”的場(chǎng)景連續(xù)編碼適合“關(guān)心變化速率”的場(chǎng)景選哪個(gè)取決于研究問題。5. 結(jié)果呈現(xiàn)與可視化建議5.1 均值軌跡圖讓交互效應(yīng)一眼可見重復(fù)測(cè)量分析的結(jié)果圖核心是展示“組均值隨時(shí)間的變化軌跡”同時(shí)把個(gè)體變異或置信區(qū)間疊上去。我常用的畫法是x軸時(shí)間點(diǎn)y軸指標(biāo)均值兩條線分別代表兩組每個(gè)點(diǎn)加誤差棒標(biāo)準(zhǔn)誤或95%置信區(qū)間底層用半透明細(xì)線畫出每個(gè)個(gè)體的軌跡。這樣既有群體趨勢(shì)又有個(gè)體異質(zhì)性的信息。library(ggplot2) dat_summary - dat_long %% group_by(group, time) %% summarise(mean mean(value), se sd(value)/sqrt(n()), .groups drop) ggplot() geom_line(data dat_long, aes(x time, y value, group id), color grey80, alpha 0.4) geom_line(data dat_summary, aes(x time, y mean, group group, color group), size 1.2) geom_point(data dat_summary, aes(x time, y mean, color group), size 3) geom_errorbar(data dat_summary, aes(x time, ymin mean - se, ymax mean se, color group), width 0.2) theme_bw()這張圖的信息密度比單純的兩條均值線高得多。審稿人或者讀者一眼就能看出兩組基線是否可比、軌跡是否交叉、個(gè)體變異有多大、哪個(gè)時(shí)間點(diǎn)開始分化。如果交互顯著圖上兩條線的分離趨勢(shì)會(huì)很明顯如果交互不顯著但時(shí)間主效應(yīng)顯著兩條線基本平行下降。5.2 球形校正結(jié)果在表格里怎么擺結(jié)果表格的呈現(xiàn)有個(gè)講究未校正和校正后的p值都要報(bào)但主結(jié)論以校正后為準(zhǔn)。我通常的表格結(jié)構(gòu)是效應(yīng)項(xiàng)組間、時(shí)間、交互、F值、未校正自由度、未校正p值、epsilon值、校正方法、校正后p值。這樣讀者既能看原始結(jié)果也能判斷校正的影響程度。如果epsilon接近1可以在表注里說明“球形假設(shè)滿足報(bào)告未校正結(jié)果”避免表格過于臃腫。效應(yīng)量方面偏eta方partial eta squared是重復(fù)測(cè)量方差分析里最常報(bào)的計(jì)算方式是效應(yīng)平方和除以效應(yīng)平方和誤差平方和。組間、組內(nèi)、交互各報(bào)一個(gè)。解釋閾值參考0.01小、0.06中、0.14大但這是通用參考具體領(lǐng)域可能有自己的慣例。5.3 簡(jiǎn)單效應(yīng)結(jié)果的文字化表達(dá)事后比較的結(jié)果別只甩一張表要用文字把關(guān)鍵對(duì)比說清楚。比如“在第8周A組均值比B組低6.2個(gè)單位95% CI: -9.8 至 -2.6校正后p 0.002而基線時(shí)兩組差異無統(tǒng)計(jì)學(xué)意義p 0.71”。這種寫法把時(shí)間點(diǎn)、組別、效應(yīng)量、置信區(qū)間、校正后p值全帶上了信息完整且可讀。如果對(duì)比很多挑研究假設(shè)相關(guān)的幾個(gè)重點(diǎn)寫其余的放補(bǔ)充材料。提示報(bào)告簡(jiǎn)單效應(yīng)時(shí)務(wù)必注明多重比較的校正方法。不同校正方法下p值差異可能很大讀者需要知道你的結(jié)論建立在哪種校正之上。6. 寫在最后的一點(diǎn)個(gè)人體會(huì)這套流程我在好幾個(gè)項(xiàng)目里反復(fù)用過從最初的照葫蘆畫瓢到后來能快速判斷“這個(gè)數(shù)據(jù)該用哪種方法”中間踩的坑基本都寫在上面了。如果只讓我留一句話給剛接觸重復(fù)測(cè)量分析的朋友那就是先把數(shù)據(jù)結(jié)構(gòu)想清楚再動(dòng)手寫代碼。寬格式長(zhǎng)格式怎么轉(zhuǎn)、因子水平怎么設(shè)、缺失值怎么處理這些“臟活”決定了后面模型跑出來可不可信。統(tǒng)計(jì)方法本身反而不復(fù)雜aov()加Error()項(xiàng)或者Anova()加idata核心代碼就那么幾行難的是對(duì)設(shè)計(jì)邏輯的理解和對(duì)結(jié)果的誠(chéng)實(shí)解讀。另外提一句重復(fù)測(cè)量方差分析不是萬能的。時(shí)間點(diǎn)很多、缺失較多、或者個(gè)體軌跡形態(tài)差異很大的時(shí)候混合效應(yīng)模型或者廣義加性模型往往更合適。經(jīng)典方法適合設(shè)計(jì)規(guī)整、時(shí)間點(diǎn)少、數(shù)據(jù)相對(duì)完整的情形超出這個(gè)范圍就該考慮升級(jí)工具了。工具沒有高低之分匹配研究問題才是關(guān)鍵。 R語言重復(fù)測(cè)量方差分析：兩組縱向數(shù)據(jù)實(shí)戰(zhàn)指南
重復(fù)測(cè)量數(shù)據(jù)在臨床研究和實(shí)驗(yàn)科學(xué)里出現(xiàn)的頻率極高——同一個(gè)受試者在不同時(shí)間點(diǎn)被反復(fù)觀測(cè)比如術(shù)后第1天、第3天、第7天的某項(xiàng)指標(biāo)或者用藥前后多個(gè)隨訪節(jié)點(diǎn)的量表得分。這類數(shù)據(jù)的特點(diǎn)是觀測(cè)值之間不獨(dú)立同一個(gè)體內(nèi)部的多次測(cè)量存在相關(guān)性如果直接拿去做普通的獨(dú)立樣本t檢驗(yàn)或者單因素方差分析p值大概率是錯(cuò)的結(jié)論也就站不住腳。我見過不少做臨床數(shù)據(jù)分析的朋友數(shù)據(jù)明明是重復(fù)測(cè)量的結(jié)構(gòu)卻硬生生把每個(gè)時(shí)間點(diǎn)拆開做多次t檢驗(yàn)結(jié)果要么假陽(yáng)性泛濫要么把本來顯著的處理效應(yīng)給淹沒了。這篇內(nèi)容就圍繞“基于R進(jìn)行兩組患者重復(fù)測(cè)量數(shù)據(jù)的分析”這個(gè)主題把重復(fù)測(cè)量方差分析的完整思路、R語言落地細(xì)節(jié)、以及我踩過的坑一次性講透。適合有一定R基礎(chǔ)、正在處理縱向隨訪數(shù)據(jù)或者交叉實(shí)驗(yàn)數(shù)據(jù)的從業(yè)者參考也適合剛接觸重復(fù)測(cè)量設(shè)計(jì)、想搞清楚“為什么不能拆開做t檢驗(yàn)”的研究人員從頭看起。1. 重復(fù)測(cè)量方差分析到底在解決什么問題1.1 從數(shù)據(jù)結(jié)構(gòu)說起為什么普通方差分析不夠用先把數(shù)據(jù)結(jié)構(gòu)擺清楚。假設(shè)我們做了一個(gè)兩組對(duì)照的臨床觀察A組用新方案B組用常規(guī)方案每組各若干名患者在基線、第4周、第8周、第12周四個(gè)時(shí)間點(diǎn)測(cè)量同一個(gè)連續(xù)型指標(biāo)。數(shù)據(jù)排成“長(zhǎng)格式”之后每一行是“某患者在某時(shí)間點(diǎn)的一次測(cè)量”包含患者ID、分組、時(shí)間點(diǎn)、測(cè)量值四個(gè)核心字段。這種結(jié)構(gòu)下同一個(gè)患者貢獻(xiàn)了4行數(shù)據(jù)這4個(gè)觀測(cè)值天然相關(guān)——一個(gè)基線值偏高的患者后續(xù)幾個(gè)時(shí)間點(diǎn)大概率也偏高。普通單因素方差分析的前提是觀測(cè)獨(dú)立這個(gè)前提在這里直接被打破了。如果無視相關(guān)性把4個(gè)時(shí)間點(diǎn)當(dāng)成4個(gè)獨(dú)立組來做方差分析誤差項(xiàng)里混入了“患者個(gè)體差異”這部分變異導(dǎo)致誤差均方被高估檢驗(yàn)效能下降。更麻煩的是時(shí)間點(diǎn)之間的比較會(huì)失真因?yàn)槟銢]法區(qū)分“時(shí)間效應(yīng)”和“個(gè)體差異”。重復(fù)測(cè)量方差分析的核心思路就是把總變異拆成三塊組間變異處理效應(yīng)、組內(nèi)變異時(shí)間效應(yīng)、以及交互效應(yīng)處理×?xí)r間同時(shí)把“個(gè)體差異”從誤差項(xiàng)里單獨(dú)剝離出來讓誤差估計(jì)更干凈。這就是它比普通方差分析更適合縱向數(shù)據(jù)的根本原因。用生活化的類比假設(shè)你想比較兩個(gè)班級(jí)學(xué)生的考試成績(jī)變化如果每個(gè)學(xué)生只考一次那是獨(dú)立樣本問題但如果是同一批學(xué)生期中、期末各考一次學(xué)生自身的水平差異就是一個(gè)巨大的干擾項(xiàng)。重復(fù)測(cè)量方差分析相當(dāng)于先把“這個(gè)學(xué)生本來就行/本來就不行”這部分變異扣掉再去看“時(shí)間帶來的變化”和“班級(jí)之間的差異”結(jié)論自然更靠譜。1.2 兩組設(shè)計(jì)的特殊性組間、組內(nèi)、交互三重效應(yīng)兩組重復(fù)測(cè)量設(shè)計(jì)里我們關(guān)心的效應(yīng)有三個(gè)層面。第一是組間效應(yīng)也就是A組和B組整體上有沒有差異這個(gè)效應(yīng)在統(tǒng)計(jì)上等價(jià)于“把每個(gè)患者的多次測(cè)量取平均再做兩組比較”但它用的是合并后的誤差項(xiàng)比直接做t檢驗(yàn)更穩(wěn)健。第二是組內(nèi)效應(yīng)也就是時(shí)間主效應(yīng)所有患者合在一起看指標(biāo)隨時(shí)間有沒有顯著變化。第三是交互效應(yīng)這是最有臨床價(jià)值的一項(xiàng)——兩組隨時(shí)間變化的趨勢(shì)是否不同。舉個(gè)例子如果新方案組在第8周之后指標(biāo)快速下降而常規(guī)組下降平緩交互效應(yīng)就會(huì)顯著說明處理方案的效果依賴于時(shí)間這比單純看“兩組有沒有差異”信息量大得多。很多初學(xué)者只盯著組間p值看忽略了交互效應(yīng)這是很可惜的。臨床上真正有意義的往往是“處理組起效更快”或者“處理組在后期維持得更好”這類趨勢(shì)差異而這些恰恰體現(xiàn)在交互項(xiàng)上。所以在R里跑重復(fù)測(cè)量方差分析時(shí)交互項(xiàng)的顯著性一定要單獨(dú)拎出來解讀不能只看主效應(yīng)。1.3 球形假設(shè)這個(gè)前提不滿足會(huì)怎樣重復(fù)測(cè)量方差分析有一個(gè)繞不開的假設(shè)——球形假設(shè)sphericity。簡(jiǎn)單說它要求任意兩個(gè)時(shí)間點(diǎn)之間差值的方差相等。四個(gè)時(shí)間點(diǎn)的話就有6對(duì)差值這6個(gè)差值的方差理論上應(yīng)該差不多。如果球形假設(shè)被違反組內(nèi)效應(yīng)和交互效應(yīng)的F檢驗(yàn)自由度需要校正否則p值偏小容易假陽(yáng)性。球形假設(shè)在時(shí)間點(diǎn)較多、或者時(shí)間間隔不均勻時(shí)特別容易出問題。我的經(jīng)驗(yàn)是三個(gè)時(shí)間點(diǎn)以內(nèi)球形假設(shè)通常還過得去四個(gè)及以上時(shí)間點(diǎn)最好默認(rèn)做校正。R里car包的Anova()函數(shù)配合idata參數(shù)可以自動(dòng)輸出GG校正Greenhouse-Geisser和HF校正Huynh-Feldt的結(jié)果直接看校正后的p值就行。這里有個(gè)細(xì)節(jié)球形假設(shè)只影響組內(nèi)效應(yīng)和交互效應(yīng)組間效應(yīng)不受影響因?yàn)榻M間比較用的是合并的個(gè)體均值不涉及時(shí)間點(diǎn)之間的差值結(jié)構(gòu)。這一點(diǎn)很多人搞混以為球形假設(shè)不滿足整個(gè)分析就廢了其實(shí)不是。2. R語言落地的整體方案與工具選型2.1 寬格式還是長(zhǎng)格式數(shù)據(jù)整形的關(guān)鍵決策在R里做重復(fù)測(cè)量方差分析數(shù)據(jù)格式的選擇直接影響后續(xù)代碼怎么寫。寬格式是每個(gè)患者一行不同時(shí)間點(diǎn)的測(cè)量值各占一列比如id, group, t1, t2, t3, t4。長(zhǎng)格式是每個(gè)患者每個(gè)時(shí)間點(diǎn)一行比如id, group, time, value。兩種格式各有用途寬格式適合用aov()配合Error()項(xiàng)做經(jīng)典的重復(fù)測(cè)量方差分析長(zhǎng)格式適合用混合效應(yīng)模型或者car::Anova()配合idata參數(shù)。我個(gè)人的建議是如果時(shí)間點(diǎn)固定且不多比如四個(gè)隨訪節(jié)點(diǎn)用寬格式配合aov()最直接代碼短、輸出清晰如果時(shí)間點(diǎn)較多、或者有缺失值、或者想順便做趨勢(shì)分析長(zhǎng)格式更靈活。實(shí)際項(xiàng)目里我通常兩份都準(zhǔn)備——用tidyr::pivot_longer()和pivot_wider()互相轉(zhuǎn)換幾行代碼的事。需要注意的是寬格式轉(zhuǎn)長(zhǎng)格式時(shí)列名要規(guī)范比如week0, week4, week8, week12轉(zhuǎn)完之后時(shí)間因子要顯式設(shè)定水平順序否則R默認(rèn)按字母排序week12會(huì)排在week4前面圖和分析全亂套。2.2 核心函數(shù)對(duì)比aov、Anova與混合模型R里能做重復(fù)測(cè)量方差分析的函數(shù)不止一個(gè)選哪個(gè)取決于你的設(shè)計(jì)復(fù)雜度和對(duì)缺失值的容忍度。stats::aov()配合Error()項(xiàng)是最經(jīng)典的做法語法是aov(value ~ group * time Error(id/time))其中Error(id/time)告訴R“個(gè)體是隨機(jī)效應(yīng)時(shí)間是個(gè)體內(nèi)因子”。這個(gè)函數(shù)的好處是輸出直接給出組間、組內(nèi)、交互三張方差分析表球形假設(shè)的校正需要額外用summary()配合Mauchly檢驗(yàn)來看。car::Anova()是另一個(gè)常用選擇它需要把寬格式數(shù)據(jù)轉(zhuǎn)成mlm對(duì)象多元線性模型然后用idata指定個(gè)體內(nèi)因子結(jié)構(gòu)。它的優(yōu)勢(shì)是自動(dòng)輸出球形校正結(jié)果而且對(duì)不平衡設(shè)計(jì)更穩(wěn)健。缺點(diǎn)是數(shù)據(jù)準(zhǔn)備稍麻煩需要構(gòu)造idata數(shù)據(jù)框。如果數(shù)據(jù)有缺失值臨床隨訪里太常見了上面兩個(gè)函數(shù)都會(huì)直接丟掉含缺失的整行樣本量縮水嚴(yán)重。這時(shí)候混合效應(yīng)模型lme4::lmer()或nlme::lme()是更優(yōu)解它用限制性極大似然估計(jì)能利用部分信息不需要完整案例。不過混合模型的方差結(jié)構(gòu)設(shè)定更復(fù)雜隨機(jī)斜率、隨機(jī)截距的選擇需要結(jié)合設(shè)計(jì)來定不是一句兩句能說完的。這篇內(nèi)容聚焦經(jīng)典重復(fù)測(cè)量方差分析混合模型留到后續(xù)單獨(dú)展開。2.3 必備包清單與版本注意事項(xiàng)跑完整套分析我常用的包有這么幾個(gè)car負(fù)責(zé)Anova()和球形校正emmeans負(fù)責(zé)事后兩兩比較和簡(jiǎn)單效應(yīng)分析tidyr和dplyr負(fù)責(zé)數(shù)據(jù)整形ggplot2負(fù)責(zé)可視化rstatix提供了一些便捷的管道式統(tǒng)計(jì)函數(shù)。版本上要注意car包3.0之后的Anova()對(duì)idata參數(shù)的校驗(yàn)更嚴(yán)格idata的行數(shù)必須和個(gè)體內(nèi)因子的水平數(shù)完全一致順序也要對(duì)上否則報(bào)錯(cuò)信息很隱晦。emmeans在1.7版本之后對(duì)重復(fù)測(cè)量模型的mode參數(shù)有調(diào)整做簡(jiǎn)單效應(yīng)時(shí)記得顯式指定mode df.error或者根據(jù)模型類型選對(duì)自由度近似方式不然置信區(qū)間會(huì)偏窄。提示裝包時(shí)建議鎖定版本尤其是car和emmeans這兩個(gè)包更新頻繁不同版本之間輸出格式和默認(rèn)行為有差異團(tuán)隊(duì)協(xié)作時(shí)最好在項(xiàng)目里記錄sessionInfo()。3. 完整實(shí)操流程從數(shù)據(jù)準(zhǔn)備到結(jié)果解讀3.1 模擬一份兩組四時(shí)點(diǎn)的臨床數(shù)據(jù)為了把流程跑通我先構(gòu)造一份模擬數(shù)據(jù)。兩組各20名患者四個(gè)時(shí)間點(diǎn)指標(biāo)值設(shè)定為A組基線均值50隨時(shí)間逐步下降到第12周的38B組基線均值50下降到第12周的45。個(gè)體隨機(jī)效應(yīng)標(biāo)準(zhǔn)差設(shè)為5殘差標(biāo)準(zhǔn)差設(shè)為3。這樣設(shè)定下交互效應(yīng)應(yīng)該是顯著的因?yàn)閮山M下降斜率不同。set.seed(2024) n_per_group - 20 times - c(0, 4, 8, 12) groups - rep(c(A, B), each n_per_group) id - 1:(2 * n_per_group) # 個(gè)體隨機(jī)截距 rand_intercept - rnorm(2 * n_per_group, 0, 5) # 構(gòu)造長(zhǎng)格式 dat_long - expand.grid(id id, time times) dat_long$group - groups[match(dat_long$id, id)] dat_long$value - 50 rand_intercept[dat_long$id] ifelse(dat_long$group A, -1.0 * dat_long$time, -0.4 * dat_long$time) rnorm(nrow(dat_long), 0, 3) dat_long$time - factor(dat_long$time, levels times) dat_long$group - factor(dat_long$group, levels c(A, B))這段代碼里有兩個(gè)細(xì)節(jié)值得說。第一time轉(zhuǎn)因子時(shí)顯式指定了levels times保證順序是0、4、8、12而不是默認(rèn)的字符串排序。第二個(gè)體隨機(jī)截距用rand_intercept[dat_long$id]按id索引確保同一個(gè)體的四次測(cè)量共享同一個(gè)隨機(jī)截距這是模擬重復(fù)測(cè)量數(shù)據(jù)的關(guān)鍵漏了這一步數(shù)據(jù)就退化成獨(dú)立觀測(cè)了。轉(zhuǎn)成寬格式用于aov()分析library(tidyr) dat_wide - pivot_wider(dat_long, names_from time, values_from value, names_prefix t)轉(zhuǎn)完之后列名是t0, t4, t8, t12每個(gè)患者一行配合group列正好是aov()需要的結(jié)構(gòu)。3.2 用aov加Error項(xiàng)跑經(jīng)典重復(fù)測(cè)量方差分析寬格式數(shù)據(jù)準(zhǔn)備好之后核心代碼就一行l(wèi)ibrary(car) mod_aov - aov(cbind(t0, t4, t8, t12) ~ group, data dat_wide) summary(mod_aov)注意這里用的是cbind()把四個(gè)時(shí)間點(diǎn)綁成一個(gè)矩陣作為因變量group是組間因子。這個(gè)寫法下summary()輸出的是組間效應(yīng)的檢驗(yàn)結(jié)果。要拿到組內(nèi)效應(yīng)和交互效應(yīng)需要把模型轉(zhuǎn)成mlm對(duì)象再用Anova()mod_mlm - lm(cbind(t0, t4, t8, t12) ~ group, data dat_wide) idata - data.frame(time factor(c(t0, t4, t8, t12), levels c(t0, t4, t8, t12))) Anova(mod_mlm, idata idata, idesign ~time, type III)idesign ~time指定個(gè)體內(nèi)因子是timeidata的行數(shù)和順序必須和cbind()里列的排列完全一致這是最容易出錯(cuò)的地方。輸出里會(huì)包含time的球形檢驗(yàn)結(jié)果Mauchlys test以及未校正和GG/HF校正后的F值和p值。組間效應(yīng)看group那一行交互效應(yīng)看group:time那一行。實(shí)測(cè)這份模擬數(shù)據(jù)交互項(xiàng)的GG校正p值應(yīng)該在0.05以下說明兩組的時(shí)間趨勢(shì)確實(shí)不同。組間主效應(yīng)可能不顯著因?yàn)閮山M在第12周的均值差距被基線到中期的重疊稀釋了——這正好說明為什么不能只看組間p值。3.3 球形檢驗(yàn)與校正結(jié)果的取舍Anova()輸出里會(huì)有一行Mauchlys test for sphericity給出W統(tǒng)計(jì)量和p值。如果p 0.05說明球形假設(shè)被違反需要看校正后的結(jié)果。GG校正Greenhouse-Geisser的epsilon值通常偏保守HF校正Huynh-Feldt偏寬松。我的取舍原則是epsilon在0.75以上用HF校正0.75以下用GG校正。這個(gè)閾值不是絕對(duì)的但實(shí)踐中比較穩(wěn)妥。還有一個(gè)細(xì)節(jié)球形檢驗(yàn)本身對(duì)樣本量敏感小樣本時(shí)檢驗(yàn)效能低可能漏檢大樣本時(shí)又過于敏感輕微偏離就報(bào)顯著。所以不能完全依賴p值判斷結(jié)合epsilon值一起看更靠譜。如果epsilon接近1比如0.95以上即使Mauchly檢驗(yàn)顯著未校正結(jié)果也基本可用。注意球形假設(shè)只約束組內(nèi)效應(yīng)和交互效應(yīng)組間效應(yīng)的F檢驗(yàn)不需要校正。解讀結(jié)果時(shí)別把組間p值也拿去做校正那是畫蛇添足。3.4 事后比較簡(jiǎn)單效應(yīng)分析與兩兩對(duì)比交互效應(yīng)顯著之后下一步是拆解——在每個(gè)時(shí)間點(diǎn)上比較兩組或者在每組內(nèi)部比較各時(shí)間點(diǎn)。這就是簡(jiǎn)單效應(yīng)分析。用emmeans包可以很干凈地完成library(emmeans) emm - emmeans(mod_mlm, ~ group | time, mode df.error) pairs(emm)~ group | time的意思是“在每個(gè)time水平上比較group”輸出的是四個(gè)時(shí)間點(diǎn)各自的兩組對(duì)比。反過來如果想在每組內(nèi)部比較時(shí)間點(diǎn)寫成~ time | group。mode df.error指定用誤差項(xiàng)自由度做近似對(duì)于經(jīng)典的重復(fù)測(cè)量設(shè)計(jì)這是合適的如果用的是混合模型要改成mode satterthwaite或kenward-roger。多重比較校正別忘了。四個(gè)時(shí)間點(diǎn)各做一次兩組比較一共4次檢驗(yàn)用Bonferroni校正的話alpha除以4。emmeans里加adjust bonferroni或者h(yuǎn)olm都行我一般用Holm因?yàn)樗鼘?duì)檢驗(yàn)之間的相關(guān)性不敏感功效比Bonferroni略高。如果交互不顯著但時(shí)間主效應(yīng)顯著那就退一步只看所有患者合并后的時(shí)間點(diǎn)兩兩比較用emmeans(mod_mlm, ~ time)配合pairs()即可。這時(shí)候不需要按組分拆因?yàn)榻换ゲ伙@著意味著兩組的趨勢(shì)可以認(rèn)為是一致的。4. 常見問題與排查技巧實(shí)錄4.1 報(bào)錯(cuò)“number of items to replace is not a multiple of replacement length”這個(gè)報(bào)錯(cuò)在構(gòu)造idata或者模擬數(shù)據(jù)時(shí)特別常見根源是向量長(zhǎng)度不匹配。比如idata寫了4行但cbind()里只綁了3列或者模擬數(shù)據(jù)時(shí)rand_intercept的長(zhǎng)度和id對(duì)不上。排查方法很簡(jiǎn)單在報(bào)錯(cuò)的那一行之前用length()和nrow()把涉及的每個(gè)對(duì)象的維度打印出來逐個(gè)核對(duì)。我習(xí)慣在關(guān)鍵步驟后加一句stopifnot(length(x) expected)做斷言早發(fā)現(xiàn)早解決比事后debug省時(shí)間。另一個(gè)隱蔽的變體是因子水平順序不一致。idata里time的水平是t0, t4, t8, t12但cbind()里列的順序是t0, t8, t4, t12R不會(huì)報(bào)錯(cuò)但結(jié)果全錯(cuò)——因?yàn)樗次恢闷ヅ洳话疵制ヅ?。這種錯(cuò)誤最難查因?yàn)榇a能跑通輸出看起來也正常。我的做法是cbind()里的列順序和idata的行順序都用同一個(gè)變量控制比如先定義time_levels - c(t0,t4,t8,t12)兩邊都引用它從源頭上杜絕錯(cuò)位。4.2 缺失值導(dǎo)致樣本量驟減怎么辦臨床隨訪數(shù)據(jù)幾乎不可能沒有缺失患者失訪、漏檢、樣本溶血各種原因都會(huì)造成某個(gè)時(shí)間點(diǎn)缺值。aov()和Anova()默認(rèn)做完整案例分析只要一個(gè)患者四個(gè)時(shí)間點(diǎn)里缺了任何一個(gè)整行就被丟掉。20人的組缺了5個(gè)樣本量直接掉四分之一檢驗(yàn)效能大打折扣。應(yīng)對(duì)策略分兩層。第一層是分析前的缺失機(jī)制判斷如果缺失完全隨機(jī)MCAR完整案例分析雖然浪費(fèi)信息但結(jié)論無偏如果缺失與觀測(cè)值有關(guān)MAR完整案例分析就有偏了這時(shí)候應(yīng)該上混合效應(yīng)模型或者多重插補(bǔ)。第二層是實(shí)操層面的補(bǔ)救用mice包做多重插補(bǔ)插補(bǔ)后對(duì)每個(gè)插補(bǔ)數(shù)據(jù)集跑一遍重復(fù)測(cè)量方差分析再用Rubin規(guī)則合并結(jié)果。這套流程代碼量不小但mice和mitml配合起來還算順手。如果缺失比例很低比如5%以內(nèi)完整案例分析的影響可以接受但要在方法部分如實(shí)報(bào)告缺失情況別藏著。4.3 交互效應(yīng)顯著但事后比較全不顯著這種情況我遇到過好幾次第一次碰到時(shí)一頭霧水——交互p值0.03但四個(gè)時(shí)間點(diǎn)分別做兩組比較校正后一個(gè)都不顯著。原因通常是交互效應(yīng)體現(xiàn)在“趨勢(shì)”上而非“某個(gè)具體時(shí)間點(diǎn)的截距差異”上。比如兩組都是線性下降但斜率不同在中間某個(gè)時(shí)間點(diǎn)可能剛好交叉單點(diǎn)比較自然不顯著但整體趨勢(shì)差異是真實(shí)的。解決辦法有兩個(gè)方向。一是做趨勢(shì)分析把時(shí)間因子拆成線性、二次、三次正交多項(xiàng)式對(duì)比看交互效應(yīng)主要落在哪個(gè)趨勢(shì)成分上。emmeans里用contrast(emm, poly)或者手動(dòng)指定正交對(duì)比系數(shù)都行。二是做區(qū)域檢驗(yàn)或者模型-based inference比如擬合一個(gè)時(shí)間×組的線性回歸直接檢驗(yàn)斜率差異。這兩種方法都比逐點(diǎn)比較更貼合“趨勢(shì)差異”的本質(zhì)。我的經(jīng)驗(yàn)是如果研究假設(shè)本來就是“起效速度不同”那從一開始就該規(guī)劃趨勢(shì)對(duì)比而不是等到事后比較不顯著了再補(bǔ)救。4.4 問題速查表問題現(xiàn)象可能原因排查與解決Anova()報(bào)錯(cuò)維度不匹配idata行數(shù)與cbind()列數(shù)不一致用length()核對(duì)統(tǒng)一用變量控制順序結(jié)果與預(yù)期完全相反因子水平順序錯(cuò)位檢查factor(levels)設(shè)定idata與列順序?qū)R組內(nèi)效應(yīng)p值異常小球形假設(shè)違反未校正看Mauchly檢驗(yàn)和epsilon選GG或HF校正事后比較全不顯著交互體現(xiàn)在趨勢(shì)而非單點(diǎn)改用正交多項(xiàng)式對(duì)比或斜率檢驗(yàn)樣本量比預(yù)期少很多缺失值完整案例刪除評(píng)估缺失比例考慮多重插補(bǔ)或混合模型emmeans置信區(qū)間過窄自由度近似方式不對(duì)顯式指定mode參數(shù)匹配模型類型4.5 幾個(gè)容易被忽略的實(shí)操心得第一個(gè)心得先畫圖再跑模型。把每個(gè)患者的軌跡畫出來ggplot2里用geom_line(aes(group id))再疊加組均值軌跡數(shù)據(jù)的形態(tài)一目了然。有沒有離群個(gè)體、缺失集中在哪個(gè)時(shí)間點(diǎn)、兩組軌跡是平行還是交叉圖上一看便知。很多統(tǒng)計(jì)上的困惑畫個(gè)圖就解開了。我現(xiàn)在的習(xí)慣是任何重復(fù)測(cè)量分析圖一定跑在模型前面。第二個(gè)心得樣本量估算要在設(shè)計(jì)階段做。重復(fù)測(cè)量設(shè)計(jì)的樣本量估算和普通兩組比較不一樣要考慮時(shí)間點(diǎn)之間的相關(guān)性用預(yù)實(shí)驗(yàn)的相關(guān)系數(shù)估計(jì)和球形假設(shè)的epsilon。pwr包和WebPower包都能做但參數(shù)設(shè)定需要預(yù)實(shí)驗(yàn)數(shù)據(jù)支撐。事后補(bǔ)算意義不大但方法部分最好報(bào)告一下估算依據(jù)。第三個(gè)心得報(bào)告結(jié)果時(shí)把效應(yīng)量和置信區(qū)間帶上。光報(bào)p值信息量不夠尤其是交互效應(yīng)報(bào)個(gè)偏eta方或者兩組均值差的置信區(qū)間讀者才能判斷臨床意義。emmeans的confint()直接給置信區(qū)間effectsize包可以算偏eta方幾行代碼的事別省。第四個(gè)心得時(shí)間因子的編碼方式影響解讀。如果時(shí)間點(diǎn)間隔不均勻比如0、4、8、12周是均勻的但0、2、6、12就不均勻把時(shí)間當(dāng)分類因子做方差分析會(huì)損失“間隔信息”。這時(shí)候可以考慮把時(shí)間當(dāng)連續(xù)變量做趨勢(shì)分析或者用正交多項(xiàng)式對(duì)比。分類因子適合“每個(gè)時(shí)間點(diǎn)都關(guān)心”的場(chǎng)景連續(xù)編碼適合“關(guān)心變化速率”的場(chǎng)景選哪個(gè)取決于研究問題。5. 結(jié)果呈現(xiàn)與可視化建議5.1 均值軌跡圖讓交互效應(yīng)一眼可見重復(fù)測(cè)量分析的結(jié)果圖核心是展示“組均值隨時(shí)間的變化軌跡”同時(shí)把個(gè)體變異或置信區(qū)間疊上去。我常用的畫法是x軸時(shí)間點(diǎn)y軸指標(biāo)均值兩條線分別代表兩組每個(gè)點(diǎn)加誤差棒標(biāo)準(zhǔn)誤或95%置信區(qū)間底層用半透明細(xì)線畫出每個(gè)個(gè)體的軌跡。這樣既有群體趨勢(shì)又有個(gè)體異質(zhì)性的信息。library(ggplot2) dat_summary - dat_long %% group_by(group, time) %% summarise(mean mean(value), se sd(value)/sqrt(n()), .groups drop) ggplot() geom_line(data dat_long, aes(x time, y value, group id), color grey80, alpha 0.4) geom_line(data dat_summary, aes(x time, y mean, group group, color group), size 1.2) geom_point(data dat_summary, aes(x time, y mean, color group), size 3) geom_errorbar(data dat_summary, aes(x time, ymin mean - se, ymax mean se, color group), width 0.2) theme_bw()這張圖的信息密度比單純的兩條均值線高得多。審稿人或者讀者一眼就能看出兩組基線是否可比、軌跡是否交叉、個(gè)體變異有多大、哪個(gè)時(shí)間點(diǎn)開始分化。如果交互顯著圖上兩條線的分離趨勢(shì)會(huì)很明顯如果交互不顯著但時(shí)間主效應(yīng)顯著兩條線基本平行下降。5.2 球形校正結(jié)果在表格里怎么擺結(jié)果表格的呈現(xiàn)有個(gè)講究未校正和校正后的p值都要報(bào)但主結(jié)論以校正后為準(zhǔn)。我通常的表格結(jié)構(gòu)是效應(yīng)項(xiàng)組間、時(shí)間、交互、F值、未校正自由度、未校正p值、epsilon值、校正方法、校正后p值。這樣讀者既能看原始結(jié)果也能判斷校正的影響程度。如果epsilon接近1可以在表注里說明“球形假設(shè)滿足報(bào)告未校正結(jié)果”避免表格過于臃腫。效應(yīng)量方面偏eta方partial eta squared是重復(fù)測(cè)量方差分析里最常報(bào)的計(jì)算方式是效應(yīng)平方和除以效應(yīng)平方和誤差平方和。組間、組內(nèi)、交互各報(bào)一個(gè)。解釋閾值參考0.01小、0.06中、0.14大但這是通用參考具體領(lǐng)域可能有自己的慣例。5.3 簡(jiǎn)單效應(yīng)結(jié)果的文字化表達(dá)事后比較的結(jié)果別只甩一張表要用文字把關(guān)鍵對(duì)比說清楚。比如“在第8周A組均值比B組低6.2個(gè)單位95% CI: -9.8 至 -2.6校正后p 0.002而基線時(shí)兩組差異無統(tǒng)計(jì)學(xué)意義p 0.71”。這種寫法把時(shí)間點(diǎn)、組別、效應(yīng)量、置信區(qū)間、校正后p值全帶上了信息完整且可讀。如果對(duì)比很多挑研究假設(shè)相關(guān)的幾個(gè)重點(diǎn)寫其余的放補(bǔ)充材料。提示報(bào)告簡(jiǎn)單效應(yīng)時(shí)務(wù)必注明多重比較的校正方法。不同校正方法下p值差異可能很大讀者需要知道你的結(jié)論建立在哪種校正之上。6. 寫在最后的一點(diǎn)個(gè)人體會(huì)這套流程我在好幾個(gè)項(xiàng)目里反復(fù)用過從最初的照葫蘆畫瓢到后來能快速判斷“這個(gè)數(shù)據(jù)該用哪種方法”中間踩的坑基本都寫在上面了。如果只讓我留一句話給剛接觸重復(fù)測(cè)量分析的朋友那就是先把數(shù)據(jù)結(jié)構(gòu)想清楚再動(dòng)手寫代碼。寬格式長(zhǎng)格式怎么轉(zhuǎn)、因子水平怎么設(shè)、缺失值怎么處理這些“臟活”決定了后面模型跑出來可不可信。統(tǒng)計(jì)方法本身反而不復(fù)雜aov()加Error()項(xiàng)或者Anova()加idata核心代碼就那么幾行難的是對(duì)設(shè)計(jì)邏輯的理解和對(duì)結(jié)果的誠(chéng)實(shí)解讀。另外提一句重復(fù)測(cè)量方差分析不是萬能的。時(shí)間點(diǎn)很多、缺失較多、或者個(gè)體軌跡形態(tài)差異很大的時(shí)候混合效應(yīng)模型或者廣義加性模型往往更合適。經(jīng)典方法適合設(shè)計(jì)規(guī)整、時(shí)間點(diǎn)少、數(shù)據(jù)相對(duì)完整的情形超出這個(gè)范圍就該考慮升級(jí)工具了。工具沒有高低之分匹配研究問題才是關(guān)鍵。