请在阅读之前阅读此内容:
对 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))
这似乎可行,但是 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