时间戳

首页 > 解决方案 > 使用 rstatix 运行单向重复测量 anova 但无法正常工作

r - 使用 rstatix 运行单向重复测量 anova 但无法正常工作

问题描述

请在阅读之前阅读此内容：

对 RMANOVA 进行组比较 Shapiro-Wilks 测试时遇到问题

如您所见，我已经解决了获取图和正态性检查的问题。现在我只是想让实际的方差分析本身起作用。据我所知，anova_test 的帮助页面说这是默认公式：

Within-Ss ANOVA (repeated measures ANOVA): y ~ w1*w2 + Error(id/(w1*w2))

所以我使用了这个命令：

anova_test(data=weight,
           formula = score ~ trial + Error(id/trial))

它似乎不起作用，只给我这个错误：

Error: Each row of output must be identified by a unique combination of keys.
Keys are shared for 144 rows:
* 1, 13, 25, 37
* 2, 14, 26, 38
* 3, 15, 27, 39
* 4, 16, 28, 40
* 5, 17, 29, 41
* 6, 18, 30, 42
* 7, 19, 31, 43
* 8, 20, 32, 44
* 9, 21, 33, 45
* 10, 22, 34, 46
* 11, 23, 35, 47
* 12, 24, 36, 48
* 49, 61, 73, 85
* 50, 62, 74, 86
* 51, 63, 75, 87
* 52, 64, 76, 88
* 53, 65, 77, 89
* 54, 66, 78, 90
* 55, 67, 79, 91
* 56, 68, 80, 92
* 57, 69, 81, 93
* 58, 70, 82, 94
* 59, 71, 83, 95
* 60, 72, 84, 96
* 97, 109, 121, 133
* 98, 110, 122, 134
* 99, 111, 123, 135
* 100, 112, 124, 136
* 101, 113, 125, 137
* 102, 114, 126, 138
* 103, 115, 127, 139
* 104, 116, 128, 140
* 105, 117, 129, 141
* 106, 118, 130, 142
* 107, 119, 131, 143
* 108, 120, 132, 144

或者，我试过这个：

anova_test(data=weight,
           formula = score ~ t1*t2*t3 + Error(id/t1*t2*t3))

但是，它只会给我另一个错误：

Error: Can't subset columns that don't exist.
x Columns `t1`, `t2`, and `t3` don't exist.

现在这很明显，因为 t1:t3 不再是列，但我不确定我的原始代码可能存在哪些其他问题。最后，我尝试通过 aov 执行此操作：

aov(data=weight, score~trial+Error(id/trial))

aov

这似乎可行，但是 1）我不能像 rstatix 那样添加重要性和 2）我觉得它比 rstatix 更难阅读。

如果有帮助，这是整个 R 脚本：

##################################################
#                   Libraries                      
##################################################

library(rstatix) # summary stats functions
library(dplyr) # piping
library(ggpubr) # dunno what this actually does
library(datarium) # loads weightloss dataset
library(ggplot2) # creates plots
library(tidyverse) # loads pivot longer

##################################################
#         Creating Summaries/Dataframes                       
##################################################

# Create Data Frame for Dataset:


weight <- structure(list(id = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L,  9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L,  12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L,  3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L), .Label = c("1", "2",  "3", "4", "5", "6", "7", "8", "9", "10", "11", "12"), class = "factor"), 
    diet = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("no", "yes"), class = "factor"), 
    exercises = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("no", "yes"
    ), class = "factor"), t1 = c(10.43, 11.59, 11.35, 11.12, 
    9.5, 9.5, 11.12, 12.51, 11.35, 11.12, 11.12, 10.2, 11.12, 
    9.96, 12.05, 8.11, 12.05, 12.05, 12.28, 10.66, 11.35, 10.2, 
    10.2, 9.5, 10.2, 12.98, 13.21, 10.2, 11.59, 12.05, 11.59, 
    12.05, 11.82, 11.12, 12.51, 11.59, 10.43, 11.35, 11.82, 10.2, 
    13.67, 10.66, 10.2, 12.05, 11.82, 10.43, 12.74, 11.35), t2 = c(13.21, 
    10.66, 11.12, 9.5, 9.73, 12.74, 12.51, 12.28, 11.59, 10.66, 
    13.44, 11.35, 12.51, 12.74, 13.67, 14.37, 14.6, 12.98, 12.05, 
    14.37, 14.6, 11.82, 14.13, 13.21, 12.51, 12.98, 11.12, 9.73, 
    13.44, 13.67, 12.98, 11.35, 12.05, 15.29, 11.82, 12.05, 12.51, 
    14.83, 13.9, 13.21, 14.13, 15.06, 12.98, 11.35, 12.51, 14.13, 
    12.74, 11.35), t3 = c(11.59, 13.21, 11.35, 11.12, 12.28, 
    10.43, 11.59, 12.74, 9.96, 11.35, 10.66, 11.12, 15.76, 16.68, 
    17.84, 14.6, 17.84, 17.61, 18.54, 16.91, 15.52, 17.38, 19, 
    14.13, 14.6, 14.6, 12.05, 15.52, 13.9, 12.74, 13.21, 14.83, 
    14.6, 10.89, 15.52, 12.98, 14.37, 15.06, 13.44, 14.13, 15.29, 
    14.6, 15.06, 15.52, 13.9, 14.37, 15.06, 15.06)), class = c("tbl_df",  "tbl", "data.frame"), row.names = c(NA, -48L)) 
    weight


# Create Data Frame for Summary Stats:
summary <- weight %>% 
  get_summary_stats(type = "mean_sd")
summary_stats <- data.frame(summary)
summary_stats

# Pivot Longer Data to Create Factors and Scores：
weight <- weight %>% 
  pivot_longer(names_to = 'trial', # creates factor (x)
               values_to = 'value', # creates value (y)
               cols = t1:t3) # finds which cols to factor
weight

##################################################
#     Initial Boxplot and Normality Checks
##################################################
  
# Plot Means in Boxplot:
ggplot(weight,
       aes(x=trial,y=value))+
  geom_boxplot()+
  labs(title = "Trial Means",
       subtitle = "Boxplot of Each Trial's Scores",
       caption = "Data obtained from 'weightloss' dataset in R.")+
  theme(plot.title = element_text(face="bold"),
        plot.caption = element_text(face = "italic"),
        panel.background = element_rect(fill = "khaki"),
        plot.background = element_rect(fill = "bisque3"))
# As can be predicted, inc w/time

# Identify Outliers (Should be None Given Boxplot):
outlier <- weight %>% 
  group_by(trial) %>% 
  identify_outliers(value)
outlier_frame <- data.frame(outlier) 
outlier_frame # none found :)

# Normality (General Shapiro):
model <- lm(value~trial,
            data = weight) # creates model
shapiro_test_weight <- shapiro_test(residuals(model)) %>% 
  add_significance() # measures Shapiro
shapiro_test_weight # non sig

# Normality (General QQ):
ggqqplot(residuals(model))+
  labs(title = "QQ Plot of Residuals") 

# Normality (Group Shapiro):
weight %>%
  group_by(trial) %>%
  shapiro_test(value) %>% 
  add_significance() #####SOMETHING WRONG HERE


# Normality (Group QQ):
ggqqplot(weight, "value", ggtheme = theme_bw())+
  facet_wrap(~trial)+
labs(title = "QQPlot of Each Trial",
     subtitle = "Normality Check of Each Testing Period",
     caption = "*Data obtained from 'weightloss' dataset in R.")+
  theme(plot.title = element_text(face="bold"),
        plot.caption = element_text(face = "italic"),
        panel.background = element_rect(fill = "khaki"),
        plot.background = element_rect(fill = "bisque3"))
#looks normal


##################################################
#         One-way repeated measures ANOVA
##################################################

### FIND ANOVA AFTER POSTING ####

anova_test(data= weight,
           formula = value ~ trial + Error(id/trial))

# Pairwise Comparison Post-Test
pairwise <- weight %>% 
  pairwise_t_test(value~trial,
                  paired=TRUE, 
                  p.adjust.method = "bonferroni") 
data.frame(pairwise) # all pairs significant

标签： rtidyverseanovarstatix

解决方案

推荐阅读