最近在做重復測量方差分析，真的是走了很多彎路，足足花費了我兩周的時間，因此在此寫一篇博文，希望能給其他人提供一些參考。
先說建議：
建議使用SPSS，不要使用R，會省很多精力，我用R做了3天，失敗了，然后改用SPSS，花了1天就搞定，一方面是因為SPSS確實對用戶很友好，而且很簡單，另一方面也是因為SPSS有很多的教程，照着用就行了，很方便。
接下來，我首先介紹我的項目背景吧，我是獲得了某月一個城市的車輛數隨時間的變化情況，然后我希望比較各天的車輛數有沒有顯著性差異，所以就做了方差分析。
數據類似於下面這樣（部分）：
其中‘周x’是一個組間因素，因為每個樣本只會做一個處理（這邊的樣本就是該城市運行的車輛），然后時間（7:00-7:15這些）就是組內因素，因為每個樣本都在兩個或兩個以上的處理下進行了測量，這就是一個十分典型的重復測量方差分析了，然后我這邊各個處理的樣例數目是不一樣的，因此說，我做的是非均衡重復測量方差分析。
前期介紹就到這里啦，如果對重復測量方差分析還有疑問的可以網上查一下，其中對於組間因素和組內因素我參考的是這篇文章：https://wenku.baidu.com/view/d17dc3130b4e767f5acfce33.html
然后結合R語言實戰中的方差分析一章節基本搞懂了。
先說明下重復測量方差分析要滿足的假設：
使用兩因素重復測量方差分析(Two-way Repeated Measures Anova)進行分析時，需要考慮5個假設。
假設1：因變量唯一，且為連續變量；
假設2：有兩個受試者內因素（Within-Subject Factor），每個受試者內因素有2個或以上的水平。
注：在重復測量的方差分析模型中，對同一個體相同變量的不同次觀測結果被視為一組，用於區分重復測量次數的變量被稱為受試者內因素，受試者內因素實際上是自變量。
對數據的假設：
假設3：受試者內因素的各個水平，因變量沒有極端異常值；
假設4：受試者內因素的各個水平，因變量需服從近似正態分布；
假設5：對於受試者內因素的各個水平組合而言，因變量的方差協方差矩陣相等，也稱為球形假設。
摘自——https://mp.weixin.qq.com/s?__biz=MzI2OTQyMzc5MA==&mid=2247488797&idx=1&sn=34f6d5c3dabf2988b7b85071ebd312f3&chksm=eae1d0dcdd9659ca501644823def4d9d40777939813a381f50f6c1a117863ca6f7eb89bae289&scene=21#wechat_redirect
然后說一下用R做的過程中踩過的一些吭吧：
1、R語言的資料真的很少，所以我基本上都是參考R語言實戰來做的，在重復測量方差分析一節直接講了怎么做重復測量方差分析，但是並沒有對數據進行檢驗，重復測量方差分析一下子就成功做出來了，很簡單，但是我看方差分析都是要做檢驗的，所以我就去網上搜檢驗的方法。（詳見前述的5個假設）
2、知道了要做這些假設后，基本也就明白了最關鍵的就是球形假設了，而球形假設最常見的就是使用mauchly test，但是網上搜啊搜，就是沒搜到這方面的資料，非常尷尬，沒辦法，只好直接上R的官方文檔了，網址如下：
https://www.rdocumentation.org/packages/stats/versions/3.4.1/topics/mauchly.test
然后要看懂這個還需要以下一些文檔：
https://www.rdocumentation.org/packages/stats/versions/3.4.1/topics/SSD
https://www.rdocumentation.org/packages/stats/versions/3.4.1/topics/mauchly.test
由於不是學統計的，對立面的一些術語也不是很理解，所以只有依靠官方給的例子了：
# Lifted from Baron+Li:
# "Notes on the use of R for psychology experiments and questionnaires"
# Maxwell and Delaney, p. 497
reacttime <- matrix(c(
420, 420, 480, 480, 600, 780,
420, 480, 480, 360, 480, 600,
480, 480, 540, 660, 780, 780,
420, 540, 540, 480, 780, 900,
540, 660, 540, 480, 660, 720,
360, 420, 360, 360, 480, 540,
480, 480, 600, 540, 720, 840,
480, 600, 660, 540, 720, 900,
540, 600, 540, 480, 720, 780,
480, 420, 540, 540, 660, 780),
ncol = 6, byrow = TRUE,
dimnames = list(subj = 1:10,
cond = c("deg0NA", "deg4NA", "deg8NA",
"deg0NP", "deg4NP", "deg8NP")))
mlmfit <- lm(reacttime ~ 1)
summary(mlmfit)
SSD(mlmfit)
$SSD
cond
cond deg0NA deg4NA deg8NA deg0NP deg4NP deg8NP
deg0NA 29160 30600 26640 23760 32400 25560
deg4NA 30600 66600 32400 7200 36000 30600
deg8NA 26640 32400 56160 41040 57600 69840
deg0NP 23760 7200 41040 70560 72000 63360
deg4NP 32400 36000 57600 72000 108000 100800
deg8NP 25560 30600 69840 63360 100800 122760
$call
lm(formula = reacttime ~ 1)
$df
[1] 9
attr(,"class")
[1] "SSD"
estVar(mlmfit)
cond
cond deg0NA deg4NA deg8NA deg0NP deg4NP deg8NP
deg0NA 3240 3400 2960 2640 3600 2840
deg4NA 3400 7400 3600 800 4000 3400
deg8NA 2960 3600 6240 4560 6400 7760
deg0NP 2640 800 4560 7840 8000 7040
deg4NP 3600 4000 6400 8000 12000 11200
deg8NP 2840 3400 7760 7040 11200 13640
### traditional test of intrasubj. contrasts
mauchly.test(mlmfit, X = ~1)
Mauchly's test of sphericity
Contrasts orthogonal to
~1'
data: SSD matrix from lm(formula = reacttime ~ 1)
W = 0.031084, p-value = 0.04765
### tests using intra-subject 3x2 design
idata <- data.frame(deg = gl(3, 1, 6, labels = c(0,4,8)),
noise = gl(2, 3, 6, labels = c("A","P")))
mauchly.test(mlmfit, X = ~ deg + noise, idata = idata)
Mauchly's test of sphericity
Contrasts orthogonal to
~deg + noise'
data: SSD matrix from lm(formula = reacttime ~ 1)
W = 0.89378, p-value = 0.6381
mauchly.test(mlmfit, M = ~ deg + noise, X = ~ noise, idata = idata)
Mauchly's test of sphericity
Contrasts orthogonal to
~noise
Contrasts spanned by
~deg + noise'
data: SSD matrix from lm(formula = reacttime ~ 1)
W = 0.96011, p-value = 0.8497
簡單說明下其中的一些東西吧，輸入一個矩陣，一列是一個sample，做一個線性回歸，但是這個線性回歸~1，其實就是求每一列的均值，形成一個object mlmfit，SSD（）函數就是求出一個矩陣，名字叫這玩意：The residual sums of squares and products matrix，我自己翻譯為“殘差平方和-積和矩陣”也不知對不對，不好意思，沒學過這玩意。
estVar(mlmfit)這個函數其實就是求出協方差矩陣啦！然后就可以做球形檢驗了。
我基本上就是做到這一步，之后就沒再做下去了，因為我發現我的例子中不滿足球形假設，因此實際上需要做出調整，但是R語言沒有給出調整后的值，所以如果堅持要使用R，還需要花費很多時間，所以就放棄了，改用SPSS，因為SPSS可以直接給出球形檢驗的結果和調整值，十分簡單方便。
SPSS的使用在此就不介紹了，因為網上介紹很多很多了，這里推薦兩個，寫的很好，很完整，大家自己去看大佬寫的東西吧：
https://www.sohu.com/a/123940962_489312
https://mp.weixin.qq.com/s?__biz=MzI2OTQyMzc5MA==&mid=2247488797&idx=1&sn=34f6d5c3dabf2988b7b85071ebd312f3&chksm=eae1d0dcdd9659ca501644823def4d9d40777939813a381f50f6c1a117863ca6f7eb89bae289&scene=21#wechat_redirect
不過，用spss我還是遇到了一個坑，在此提醒大家注意下：
如果重復檢測次數過多，也就是組內因素time設的過多時，球形檢驗那一項的sig會是空的，因為無法計算出來。原因是自由度太多了，但是給的樣例又太少了，這樣就無法計算出結果，遇到這個問題時，需要減少time的檢測次數，或是增加樣例數，但是一般增加樣例數比較難，所以就減少重復檢測次數就行啦，然后就可以成功計算出球形檢驗值啦！
其他應該沒什么需要注意的了，希望大家一切順利。