运行 将 R 中的方差分析与两个因子内和一个因子间混合时出错

Error when running mixed anova in R with two within factor and one between factor

我想要一些帮助来使用 R 来分析一个实验,在这个实验中,三组参与者分别被展示了两种类型的刺激物三次。因变量是连续测量。这是数据集的外观示例。

 SubjectID  Group        Trial StimType    Measure
1          1 group1 trial3_stimA        A 0.55908866
2          2 group1 trial3_stimA        A 0.98884446
3          3 group2 trial3_stimA        A 0.00000000
4          4 group2 trial3_stimA        A 0.27067991
5          5 group3 trial3_stimA        A 0.37169285
6          6 group3 trial3_stimA        A 0.42113984
7          1 group1 trial3_stimB        B 0.00000000
8          2 group1 trial3_stimB        B 0.49892807
9          3 group2 trial3_stimB        B 0.14602589
10         4 group2 trial3_stimB        B 0.50946555
11         5 group3 trial3_stimB        B 0.25572820
12         6 group3 trial3_stimB        B 0.22932966
13         1 group1 trial1_stimA        A 0.42207604
14         2 group1 trial1_stimA        A 0.85599588
15         3 group2 trial1_stimA        A 0.36428381
16         4 group2 trial1_stimA        A 0.46679336
17         5 group3 trial1_stimA        A 0.69379734
18         6 group3 trial1_stimA        A 0.55607716
19         1 group1 trial1_stimB        B 0.24261465
20         2 group1 trial1_stimB        B 0.35176384
21         3 group2 trial1_stimB        B 0.21116215
22         4 group2 trial1_stimB        B 0.33112544
23         5 group3 trial1_stimB        B 0.00000000
24         6 group3 trial1_stimB        B 0.00000000
25         1 group1 trial2_stimA        A 0.05506943
26         2 group1 trial2_stimA        A 0.22537470
27         3 group2 trial2_stimA        A 0.00000000
28         4 group2 trial2_stimA        A 0.18511144
29         5 group3 trial2_stimA        A 0.15586156
30         6 group3 trial2_stimA        A 0.04467100
31         1 group1 trial2_stimB        B 0.03890585
32         2 group1 trial2_stimB        B 0.29787709
33         3 group2 trial2_stimB        B 0.00000000
34         4 group2 trial2_stimB        B 0.28971992
35         5 group3 trial2_stimB        B 0.12993238
36         6 group3 trial2_stimB        B 0.05066011

这是我的数据结构

'data.frame':   36 obs. of  5 variables:
 $ SubjectID: Factor w/ 6 levels "1","2","3","4",..: 1 2 3 4 5 6 1 2 3 4 ...
 $ Group    : Factor w/ 3 levels "group1","group2",..: 1 1 2 2 3 3 1 1 2 2 ...
 $ Trial    : Factor w/ 6 levels "trial1_stimA",..: 5 5 5 5 5 5 6 6 6 6 ...
 $ StimType : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 2 2 2 2 ...
 $ Measure  : num  0.559 0.989 0 0.271 0.372 ...

我需要 运行 混合方差分析,其中组作为受试者因素和试验之间的因素,刺激类型作为受试者因素内的因素。我已经尝试使用三个不同的 R 包和不同的语法,但 R 要么 returns 一条错误消息,要么输出缺少交互组 x 试验 x stim 类型。

例如,当我使用 rstatix 包中的 anova_test() 时

#Trying mixed ANOVA with rstatix
> 
> mixed.anova <- anova_test(
+   data = prepared_data, dv = Measure, wid = SubjectID,
+   between = Group, within = c(Trial,StimType)
+ )
Error in check.imatrix(X.design) : 
  Terms in the intra-subject model matrix are not orthogonal.
> get_anova_table(all_subjects)
Error in is.data.frame(x) : object 'all_subjects' not found

当我使用 afex 包中的 aov 时

> #using afex
> 
> 
> mixed.anova2 <- aov_car(Measure ~ Group*Trial*StimType + Error(1|SubjectID/(Trial*StimType)), prepared_data)
Error: Empty cells in within-subjects design  (i.e., bad data structure).
table(data[c("Trial", "StimType")])
#               StimType
# Trial          A B
#   trial1_stimA 6 0
#   trial1_stimB 0 6
#   trial2_stimA 6 0
#   trial2_stimB 0 6
#   trial3_stimA 6 0
#   trial3_stimB 0 6
> 
> aov.bww

Call:
aov(formula = SCR ~ Group * Trial * CSType + Error(SubjectID) + 
    Group, data = sixPhasesAbs2)

Grand Mean: 397.1325

Stratum 1: SubjectID

Terms:
                     Group  Residuals
Sum of Squares   187283464 6399838881
Deg. of Freedom          2         81

Residual standard error: 8888.777
18 out of 20 effects not estimable
Estimated effects may be unbalanced

Stratum 2: Within

Terms:
                      Trial Group:Trial   Residuals
Sum of Squares    158766098   374583941 12799639740
Deg. of Freedom           5          10         405

Residual standard error: 5621.748
12 out of 27 effects not estimable
Estimated effects may be unbalanced

最后,我尝试使用 lmer

> #using lmer
> model = lmer(Measure ~Group*Trial*StimType +(1|SubjectID), data = prepared_data)
fixed-effect model matrix is rank deficient so dropping 18 columns / coefficients
> 
> summary(model)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Measure ~ Group * Trial * StimType + (1 | SubjectID)
   Data: prepared_data

REML criterion at convergence: -16.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.1854 -0.5424  0.0000  0.5424  1.1854 

Random effects:
 Groups    Name        Variance Std.Dev.
 SubjectID (Intercept) 0.024226 0.15565 
 Residual              0.006861 0.08283 
Number of obs: 36, groups:  SubjectID, 6

Fixed effects:
                               Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                    0.639036   0.124672  4.459241   5.126 0.005095 ** 
Groupgroup2                   -0.223497   0.176313  4.459241  -1.268 0.267088    
Groupgroup3                   -0.014099   0.176313  4.459241  -0.080 0.939728    
Trialtrial1_stimB             -0.341847   0.082829 15.000000  -4.127 0.000896 ***
Trialtrial2_stimA             -0.498814   0.082829 15.000000  -6.022 2.34e-05 ***
Trialtrial2_stimB             -0.470644   0.082829 15.000000  -5.682 4.35e-05 ***
Trialtrial3_stimA              0.134931   0.082829 15.000000   1.629 0.124124    
Trialtrial3_stimB             -0.389572   0.082829 15.000000  -4.703 0.000283 ***
Groupgroup2:Trialtrial1_stimB  0.197452   0.117138 15.000000   1.686 0.112550    
Groupgroup3:Trialtrial1_stimB -0.283091   0.117138 15.000000  -2.417 0.028865 *  
Groupgroup2:Trialtrial2_stimA  0.175831   0.117138 15.000000   1.501 0.154095    
Groupgroup3:Trialtrial2_stimA -0.025857   0.117138 15.000000  -0.221 0.828271    
Groupgroup2:Trialtrial2_stimB  0.199966   0.117138 15.000000   1.707 0.108413    
Groupgroup3:Trialtrial2_stimB -0.063997   0.117138 15.000000  -0.546 0.592870    
Groupgroup2:Trialtrial3_stimA -0.415129   0.117138 15.000000  -3.544 0.002946 ** 
Groupgroup3:Trialtrial3_stimA -0.363452   0.117138 15.000000  -3.103 0.007276 ** 
Groupgroup2:Trialtrial3_stimB  0.301779   0.117138 15.000000   2.576 0.021070 *  
Groupgroup3:Trialtrial3_stimB  0.007164   0.117138 15.000000   0.061 0.952043    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 18 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

fit warnings:
fixed-effect model matrix is rank deficient so dropping 18 columns / coefficients
> 
> anova(model)
Missing cells for: Trialtrial1_stimB:StimTypeA, Trialtrial2_stimB:StimTypeA, Trialtrial3_stimB:StimTypeA, Trialtrial1_stimA:StimTypeB, Trialtrial2_stimA:StimTypeB, Trialtrial3_stimA:StimTypeB, Groupgroup1:Trialtrial1_stimB:StimTypeA, Groupgroup2:Trialtrial1_stimB:StimTypeA, Groupgroup3:Trialtrial1_stimB:StimTypeA, Groupgroup1:Trialtrial2_stimB:StimTypeA, Groupgroup2:Trialtrial2_stimB:StimTypeA, Groupgroup3:Trialtrial2_stimB:StimTypeA, Groupgroup1:Trialtrial3_stimB:StimTypeA, Groupgroup2:Trialtrial3_stimB:StimTypeA, Groupgroup3:Trialtrial3_stimB:StimTypeA, Groupgroup1:Trialtrial1_stimA:StimTypeB, Groupgroup2:Trialtrial1_stimA:StimTypeB, Groupgroup3:Trialtrial1_stimA:StimTypeB, Groupgroup1:Trialtrial2_stimA:StimTypeB, Groupgroup2:Trialtrial2_stimA:StimTypeB, Groupgroup3:Trialtrial2_stimA:StimTypeB, Groupgroup1:Trialtrial3_stimA:StimTypeB, Groupgroup2:Trialtrial3_stimA:StimTypeB, Groupgroup3:Trialtrial3_stimA:StimTypeB.  
Interpret type III hypotheses with care.
Type III Analysis of Variance Table with Satterthwaite's method
                      Sum Sq  Mean Sq NumDF DenDF F value    Pr(>F)    
Group                0.00723 0.003614     2     3  0.5267 0.6367151    
Trial                0.96171 0.192343     5    15 28.0356 4.172e-07 ***
Group:Trial          0.44102 0.044102    10    15  6.4283 0.0007411 ***
StimType                                                               
Group:StimType                                                         
Trial:StimType                                                         
Group:Trial:StimType                                                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我搜索过类似的问题,但没有找到任何可以帮助我解决这个问题的答案。我怀疑(鉴于错误消息)我可能必须更改数据集的结构,但我不知道该怎么做,因为我是 R 的初学者。我如何 运行 混合方差分析数据集?

由于 Trial 包含有关 TrialStimType 的信息,因此固定效应模型矩阵秩亏。从原始 post 中的 afex::aov_car() 输出可以看出这一点,它说明了空单元格是 Trial 的效果,其中包含来自两个不同变量的信息。

> mixed.anova2 <- aov_car(Measure ~ Group*Trial*StimType + Error(1|SubjectID/(Trial*StimType)), prepared_data)
Error: Empty cells in within-subjects design  (i.e., bad data structure).
table(data[c("Trial", "StimType")])
#               StimType
# Trial          A B
#   trial1_stimA 6 0
#   trial1_stimB 0 6
#   trial2_stimA 6 0
#   trial2_stimB 0 6
#   trial3_stimA 6 0
#   trial3_stimB 0 6
> 
> aov.bww

我们可以通过将 Trial 变量拆分为两个变量来纠正排名不足。

使用 lme4 包和 tidyr::separate() 将试验值与刺激值分开,分析如下所示:

textFile <- "rowId SubjectID  Group        Trial StimType    Measure
1          1 group1 trial3_stimA        A 0.55908866
2          2 group1 trial3_stimA        A 0.98884446
3          3 group2 trial3_stimA        A 0.00000000
4          4 group2 trial3_stimA        A 0.27067991
5          5 group3 trial3_stimA        A 0.37169285
6          6 group3 trial3_stimA        A 0.42113984
7          1 group1 trial3_stimB        B 0.00000000
8          2 group1 trial3_stimB        B 0.49892807
9          3 group2 trial3_stimB        B 0.14602589
10         4 group2 trial3_stimB        B 0.50946555
11         5 group3 trial3_stimB        B 0.25572820
12         6 group3 trial3_stimB        B 0.22932966
13         1 group1 trial1_stimA        A 0.42207604
14         2 group1 trial1_stimA        A 0.85599588
15         3 group2 trial1_stimA        A 0.36428381
16         4 group2 trial1_stimA        A 0.46679336
17         5 group3 trial1_stimA        A 0.69379734
18         6 group3 trial1_stimA        A 0.55607716
19         1 group1 trial1_stimB        B 0.24261465
20         2 group1 trial1_stimB        B 0.35176384
21         3 group2 trial1_stimB        B 0.21116215
22         4 group2 trial1_stimB        B 0.33112544
23         5 group3 trial1_stimB        B 0.00000000
24         6 group3 trial1_stimB        B 0.00000000
25         1 group1 trial2_stimA        A 0.05506943
26         2 group1 trial2_stimA        A 0.22537470
27         3 group2 trial2_stimA        A 0.00000000
28         4 group2 trial2_stimA        A 0.18511144
29         5 group3 trial2_stimA        A 0.15586156
30         6 group3 trial2_stimA        A 0.04467100
31         1 group1 trial2_stimB        B 0.03890585
32         2 group1 trial2_stimB        B 0.29787709
33         3 group2 trial2_stimB        B 0.00000000
34         4 group2 trial2_stimB        B 0.28971992
35         5 group3 trial2_stimB        B 0.12993238
36         6 group3 trial2_stimB        B 0.05066011
"
data <- read.table(text = textFile,header = TRUE)

library(dplyr)
library(tidyr)
# split trial from Stim
data %>% separate(Trial,into = c("trialName","stimName")) -> data
library(lme4)
model = lmer(Measure ~Group*trialName*StimType +(1|SubjectID), data = data)

summary(model)

...输出:

> summary(model)
Linear mixed model fit by REML ['lmerMod']
Formula: Measure ~ Group * trialName * StimType + (1 | SubjectID)
   Data: data

REML criterion at convergence: -16.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.1854 -0.5424  0.0000  0.5424  1.1854 

Random effects:
 Groups    Name        Variance Std.Dev.
 SubjectID (Intercept) 0.024226 0.15565 
 Residual              0.006861 0.08283 
Number of obs: 36, groups:  SubjectID, 6

Fixed effects:
                                      Estimate Std. Error t value
(Intercept)                            0.63904    0.12467   5.126
Groupgroup2                           -0.22350    0.17631  -1.268
Groupgroup3                           -0.01410    0.17631  -0.080
trialNametrial2                       -0.49881    0.08283  -6.022
trialNametrial3                        0.13493    0.08283   1.629
StimTypeB                             -0.34185    0.08283  -4.127
Groupgroup2:trialNametrial2            0.17583    0.11714   1.501
Groupgroup3:trialNametrial2           -0.02586    0.11714  -0.221
Groupgroup2:trialNametrial3           -0.41513    0.11714  -3.544
Groupgroup3:trialNametrial3           -0.36345    0.11714  -3.103
Groupgroup2:StimTypeB                  0.19745    0.11714   1.686
Groupgroup3:StimTypeB                 -0.28309    0.11714  -2.417
trialNametrial2:StimTypeB              0.37002    0.11714   3.159
trialNametrial3:StimTypeB             -0.18266    0.11714  -1.559
Groupgroup2:trialNametrial2:StimTypeB -0.17332    0.16566  -1.046
Groupgroup3:trialNametrial2:StimTypeB  0.24495    0.16566   1.479
Groupgroup2:trialNametrial3:StimTypeB  0.51946    0.16566   3.136
Groupgroup3:trialNametrial3:StimTypeB  0.65371    0.16566   3.946

Correlation matrix not shown by default, as p = 18 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

> 

更新:解析效果名称

在对我的回答的评论中,发问者解释说,用于描述试验和刺激类型的数据与 post 随问题一起提供的数据中表示的数据不同。鉴于评论的内容,这里有一个 dplyr 的方法来解析试验条件和刺激类型。

# solving the data cleaning problem in comments

df <- data.frame(treatmentName = c("MeanCondPlusLate",
                   "MeanCondMinusLate", 
                   "MeanExtPlusLate", 
                   "MeanExtMinusLate", 
                   "TestPlus1",
                   "TestMinus1"))

library(dplyr)
df %>% mutate(trial = case_when(grepl("Cond",treatmentName) == TRUE ~ 1,
                                grepl("Ext",treatmentName) == TRUE ~ 2,
                                grepl("Test",treatmentName) == TRUE ~ 3), 
              stimLevel = case_when(grepl("Plus",treatmentName) == TRUE ~ "A",
                                    grepl("Minus",treatmentName) == TRUE ~ "B"))

...输出:

      treatmentName trial stimLevel
1  MeanCondPlusLate     1         A
2 MeanCondMinusLate     1         B
3   MeanExtPlusLate     2         A
4  MeanExtMinusLate     2         B
5         TestPlus1     3         A
6        TestMinus1     3         B
>