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