I am trying to run an anova test using the anova_test() function. I keep getting one of two errors.
It tells me "Error: Can't subset columns that don't exist. x Column trait doesn't exist" for columns that are clearly in my data sheet.
Or I get "Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) : contrasts can be applied only to factors with 2 or more levels" but I all of my factors have over 2 levels (I triple checked).
I am adding my code with sample data below. My real data is similar but with many more data points. I have tried the solutions people gave for similar issues but none of them work.
If anyone knows how to fix these errors please let me know!
---
title: "morphometric traits ltreb analysis-sample"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.
When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
```{r #load data}
library(dplyr)
library(tidyverse)
library(ggpubr)
library(rstatix)
sample_data = data.frame(place = c('TX','TX', 'TX', 'TX', 'TX', 'TX','TX','CA','CA','CA','CA','CA','CA','FL','FL','FL','FL','FL','FL','FL'),
date = c('1/5/15', '1/5/15', '1/5/15', '1/4/18','1/4/18','1/5/21','1/5/21','1/6/15','1/6/15','1/5/18', '1/5/18', '1/8/19', '1/8/19', '5/25/15', '5/25/15', '1/7/18', '1/7/18', '5/14/21','5/14/21','5/14/21'),
id = c(1,2,3,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,3),
trait = c(5.2,4.5,2.3,6.7,2.4,4.5,8.9, 6.7,3.5,6.7,0.7,4.5,1.2,4.5,6.7,0.5,6.7,4.5,2.3,4.5))
```
```{r #sample_data preparation}
sample_data %>% sample_n_by(place, size = 1)
```
```{r #convert id and date and trait into factor variables}
sample_data <- sample_data %>%
convert_as_factor(id, date, place)
# Inspect some random rows of the sample_data by groups
set.seed(123)
sample_data %>% sample_n_by(place, date, size = 1)
```
```{r #summary stat}
sample_data %>%
group_by(date, place) %>%
get_summary_stats(trait, type = "mean_sd")
```
```{r #visualization}
bxp <- ggboxplot(sample_data, x = "date", y = "trait", color = "place", palette = "jco")
bxp
```
```{r #check assumptions - outliers}
outliers = sample_data %>%
group_by(date, place) %>%
identify_outliers(trait)
outliers
```
```{r #normality test: If the sample_data is normally distributed, the p-value should be greater than 0.05.}
normality = sample_data %>%
group_by(date, place) %>%
shapiro_test(trait)
normality
```
```{r #Homogneity of variance assumption}
homogeneity = sample_data %>%
group_by(date) %>%
levene_test(trait ~ place)
homogeneity
```
#If this test is statistically significant (i.e., p < 0.001), you do not have equal covariances, but if the test is not statistically significant (i.e., p > 0.001), you have equal covariances and you have not violated the assumption of homogeneity of covariances.}
```{r #Homogeneity of covariances assumption}
box_m(sample_data[, "trait", drop = FALSE], sample_data$place)
```
```{r #Assumption of sphericity: Two-way mixed ANOVA test, p < 0.0001}
res.aov <- anova_test(data, dv = trait, wid = id, between = place, within = date)
get_anova_table(res.aov)
```
```{r #post hoc tests: significant two way interaction, Effect of place at each time point}
one.way <- data %>%
group_by(date) %>%
anova_test(dv = trait, wid = id, within = date) %>%
get_anova_table() %>%
adjust_pvalue(method = "bonferroni")
one.way
```
```{r # Pairwise comparisons between group levels}
pwc <- data %>%
group_by(date) %>%
pairwise_t_test(trait ~ place, p.adjust.method = "bonferroni")
pwc
```
```{r # Effect of time at each level of place}
one.way2 <- data %>%
mutate(place = 'all', date = 'all', id = 'all', trait = 'all') %>%
group_by(place) %>%
anova_test(dv = trait, wid = id, within = date) %>%
get_anova_table() %>%
adjust_pvalue(method = "bonferroni")
one.way2
```
```{r # Pairwise comparisons between time points at each place levels: Paired t-test is used because we have repeated measures by time}
pwc2 <- data %>%
group_by(place) %>%
pairwise_t_test(
trait ~ place, paired = TRUE,
p.adjust.method = "bonferroni"
) %>%
select(-df, -statistic, -p) # Remove details
pwc2
```
```{r # Visualization: boxplots with p-values}
pwc <- pwc %>% add_xy_position(x = "date")
pwc.filtered <- pwc %>% filter(time != "2015")
bxp +
stat_pvalue_manual(pwc.filtered, tip.length = 0, hide.ns = TRUE) +
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)
```