library(emmeans) ### estimated marginal means and p values
library(sjstats) ### partial eta squared and cohens f effect size
library(lme4) ### estimated the multi level model (random intercept for participants)
## Loading required package: Matrix
library(lmerTest) ### gives more comprehsive anova output with p values
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(MuMIn) ### R2 for the model
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.4 v dplyr 1.0.2
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x tidyr::pack() masks Matrix::pack()
## x tidyr::unpack() masks Matrix::unpack()
library(ggpubr)
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
library(dplyr)
library(readxl)
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
library(ggsignif)
library(sjPlot)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
#load data for One way linear mixed models, repeated measures ##########################
facw <- read_excel("C:/Users/maber/Documents/Clover Valley/CLOVER VALLEY/veg data/LPI/FACW_cover.xlsx")
view(facw)
head(facw, 3)
## # A tibble: 3 x 4
## plots transect Year Total
## <chr> <dbl> <dbl> <dbl>
## 1 D11_NORTH 1 1 20
## 2 D11_SOUTH 1 1 12
## 3 D13_NORTH 1 1 0
#factor data#############tell it you have factors and what they are called##############
facw$Year <- factor(facw$Year,
levels = c(1, 2, 3),
labels = c("2018", "2019", "2020"))
facw$transect <- factor(facw$transect,
levels = c(1, 2,3),
labels = c("Lower Meadow", "Middle Meadow", "Upper Meadow"))
###R knows your factors, the below command brakes down your variables
summary(facw)
## plots transect Year Total
## Length:74 Lower Meadow :18 2018:22 Min. : 0.00
## Class :character Middle Meadow:32 2019:26 1st Qu.: 4.00
## Mode :character Upper Meadow :24 2020:26 Median : 8.00
## Mean :12.89
## 3rd Qu.:20.00
## Max. :54.00
## NA's :1
#run linear model with cover as dependent varible, year as #############################
#independent var. for the plot by stating + (1| plot), This makes it a repeated measures
#we call our linear model lmer1
#set up the model use the following code
lmer2 <- lmer(Total ~ as.factor(Year) * as.factor(transect) + (1|plots),
data = facw)
#linear model dep variable predicted by the indep variable (year)#####################
#run the model with the following code
anova(lmer2)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## as.factor(Year) 3458.0 1729.01 2 43.895 24.6886
## as.factor(transect) 596.6 298.28 2 31.365 4.2592
## as.factor(Year):as.factor(transect) 505.1 126.29 4 44.514 1.8032
## Pr(>F)
## as.factor(Year) 6.542e-08 ***
## as.factor(transect) 0.0231 *
## as.factor(Year):as.factor(transect) 0.1451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(Total ~ as.factor(Year) * as.factor(transect) + Error (plots),
data = facw))
##
## Error: plots
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Year) 2 49 24.4 0.145 0.866
## as.factor(transect) 2 797 398.4 2.372 0.120
## as.factor(Year):as.factor(transect) 2 395 197.6 1.176 0.330
## Residuals 19 3192 168.0
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Year) 2 3321 1660.4 24.884 1.08e-07 ***
## as.factor(transect) 2 475 237.6 3.561 0.038 *
## as.factor(Year):as.factor(transect) 4 550 137.5 2.060 0.105
## Residuals 39 2602 66.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#post Hoc comparisons###################################################
#emmeans(lmer2, list(pairwise ~ Year), adjust = "bonferroni")
emmeans(lmer2, list(pairwise ~ Year), adjust = "tukey")
## Warning in model.frame.default(formula, data = data, ...): variable 'Year' is
## not a factor
## Warning in model.frame.default(formula, data = data, ...): variable 'transect'
## is not a factor
## NOTE: Results may be misleading due to involvement in interactions
## $`emmeans of Year`
## Year emmean SE df lower.CL upper.CL
## 2018 11.43 2.14 58.4 7.14 15.7
## 2019 23.01 2.06 55.7 18.87 27.1
## 2020 6.28 2.04 55.4 2.18 10.4
##
## Results are averaged over the levels of: transect
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $`pairwise differences of Year`
## 1 estimate SE df t.ratio p.value
## 2018 - 2019 -11.58 2.53 43.7 -4.570 0.0001
## 2018 - 2020 5.15 2.53 44.2 2.036 0.1156
## 2019 - 2020 16.73 2.43 41.1 6.877 <.0001
##
## Results are averaged over the levels of: transect
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
#emmeans(lmer2, list(pairwise ~ Year), adjust = "holm")
# if there was an interaction your code would be this###################
emmeans(lmer2, pairwise ~ Year | transect)
## Warning in model.frame.default(formula, data = data, ...): variable 'Year' is
## not a factor
## Warning in model.frame.default(formula, data = data, ...): variable 'transect'
## is not a factor
## $emmeans
## transect = Lower Meadow:
## Year emmean SE df lower.CL upper.CL
## 2018 18.41 4.03 61.8 10.349 26.5
## 2019 32.52 4.08 58.0 24.349 40.7
## 2020 8.85 4.08 58.0 0.683 17.0
##
## transect = Middle Meadow:
## Year emmean SE df lower.CL upper.CL
## 2018 4.22 3.49 62.5 -2.763 11.2
## 2019 20.01 3.00 60.3 14.010 26.0
## 2020 5.50 2.88 59.0 -0.267 11.3
##
## transect = Upper Meadow:
## Year emmean SE df lower.CL upper.CL
## 2018 11.66 3.47 63.2 4.721 18.6
## 2019 16.49 3.53 58.9 9.428 23.6
## 2020 4.49 3.53 58.9 -2.572 11.6
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## transect = Lower Meadow:
## contrast estimate SE df t.ratio p.value
## 2018 - 2019 -14.11 4.94 44.0 -2.854 0.0177
## 2018 - 2020 9.56 4.94 44.0 1.934 0.1412
## 2019 - 2020 23.67 4.83 40.8 4.898 <.0001
##
## transect = Middle Meadow:
## contrast estimate SE df t.ratio p.value
## 2018 - 2019 -15.80 3.97 43.2 -3.983 0.0007
## 2018 - 2020 -1.28 3.92 44.5 -0.327 0.9429
## 2019 - 2020 14.51 3.52 42.1 4.123 0.0005
##
## transect = Upper Meadow:
## contrast estimate SE df t.ratio p.value
## 2018 - 2019 -4.83 4.33 46.2 -1.115 0.5101
## 2018 - 2020 7.17 4.33 46.2 1.654 0.2340
## 2019 - 2020 12.00 4.18 40.8 2.868 0.0175
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
#Visulize data #########################################################
#ggboxplot(facw, x = "Year", y = "Total" , title = "Total % Cover of Facultative Wet Species Across All Years",
# ylab = " % Cover FACW Species",
# color = "transect", ggtheme = theme_gray(base_size = 14))+
#stat_compare_means(comparisons = my_comparisons, label.y = c(65, 75, 80))+
#stat_compare_means(label.y = 80)
ggboxplot(facw, x = "Year", y = "Total", add = "jitter", title = "Total % Cover of Facultative Wet Species Across All Years",
ylab = " % Cover FACW Species",
color = "transect", ggtheme = theme_gray(base_size = 14))+
geom_hline(yintercept = mean(facw$Total), linetype = 2)+ # Add horizontal line at base mean
stat_compare_means(method = "anova", label.y = 85)+ # Add global annova p-value
stat_compare_means(label = "p.signif", method = "t.test",
ref.group = ".all.") # Pairwise comparison against all
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1 rows containing non-finite values (stat_compare_means).
## Warning: Removed 1 rows containing non-finite values (stat_compare_means).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_hline).
#stat_compare_means(comparisons = my_comparisons, label.y = c(65, 75, 80))+
#stat_compare_means(label.y = 80)