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)
#load data for One way linear mixed models, repeated measures ##########################
glm_test <- read_excel("C:/Users/maber/Documents/Clover Valley/CLOVER VALLEY/veg data/LPI/UPL_cover.xlsx")
view(glm_test)
head(glm_test, 3)
## # A tibble: 3 x 4
## plot transect year cover
## <chr> <dbl> <dbl> <dbl>
## 1 D11_NORTH 1 1 24
## 2 D11_SOUTH 1 1 16
## 3 R27_NORTH 1 1 18
#factor data#############tell it you have factors and what they are called##############
glm_test$year <- factor(glm_test$year,
levels = c(1, 2, 3),
labels = c("2018", "2019", "2020"))
glm_test$transect <- factor(glm_test$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(glm_test)
## plot transect year cover
## Length:66 Lower Meadow :18 2018:22 Min. : 0.00
## Class :character Middle Meadow:24 2019:22 1st Qu.: 2.00
## Mode :character Upper Meadow :24 2020:22 Median :15.00
## Mean :21.24
## 3rd Qu.:35.50
## Max. :82.00
#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(cover ~ as.factor(year) * as.factor(transect) + (1|plot),
data = glm_test)
#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) 3526.5 1763.27 2 38 8.6301
## as.factor(transect) 1456.5 728.25 2 19 3.5643
## as.factor(year):as.factor(transect) 730.5 182.64 4 38 0.8939
## Pr(>F)
## as.factor(year) 0.0008129 ***
## as.factor(transect) 0.0484784 *
## as.factor(year):as.factor(transect) 0.4770714
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(cover ~ as.factor(year) * as.factor(transect) + Error (plot),
data = glm_test))
##
## Error: plot
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(transect) 2 5084 2542.1 3.564 0.0485 *
## Residuals 19 13551 713.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(year) 2 3964 1982.1 9.701 0.000395 ***
## as.factor(year):as.factor(transect) 4 731 182.6 0.894 0.477071
## Residuals 38 7764 204.3
## ---
## 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 31.8 4.16 40.4 23.40 40.2
## 2019 14.5 4.16 40.4 6.12 22.9
## 2020 18.6 4.16 40.4 10.18 27.0
##
## 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 17.28 4.35 38 3.972 0.0009
## 2018 - 2020 13.22 4.35 38 3.040 0.0116
## 2019 - 2020 -4.06 4.35 38 -0.932 0.6234
##
## 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
# 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 30.67 7.89 40.4 14.72 46.6
## 2019 23.33 7.89 40.4 7.38 39.3
## 2020 24.00 7.89 40.4 8.05 40.0
##
## transect = Middle Meadow:
## year emmean SE df lower.CL upper.CL
## 2018 19.75 6.84 40.4 5.94 33.6
## 2019 2.75 6.84 40.4 -11.06 16.6
## 2020 6.75 6.84 40.4 -7.06 20.6
##
## transect = Upper Meadow:
## year emmean SE df lower.CL upper.CL
## 2018 45.00 6.84 40.4 31.19 58.8
## 2019 17.50 6.84 40.4 3.69 31.3
## 2020 25.00 6.84 40.4 11.19 38.8
##
## 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 7.333 8.25 38 0.889 0.6507
## 2018 - 2020 6.667 8.25 38 0.808 0.7005
## 2019 - 2020 -0.667 8.25 38 -0.081 0.9964
##
## transect = Middle Meadow:
## contrast estimate SE df t.ratio p.value
## 2018 - 2019 17.000 7.15 38 2.379 0.0572
## 2018 - 2020 13.000 7.15 38 1.819 0.1771
## 2019 - 2020 -4.000 7.15 38 -0.560 0.8421
##
## transect = Upper Meadow:
## contrast estimate SE df t.ratio p.value
## 2018 - 2019 27.500 7.15 38 3.848 0.0013
## 2018 - 2020 20.000 7.15 38 2.798 0.0214
## 2019 - 2020 -7.500 7.15 38 -1.049 0.5509
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
#Visulize data #########################################################
#ggboxplot(glm_test, x = "year", y = "cover", title = "Total % Cover of Upland Species Across All Years",
# ylab = " % Cover UPL Species",
# color = "transect", ggtheme = theme_gray(base_size = 14))+
#stat_compare_means(comparisons = my_comparisons, label.y = c(85, 90, 95))+
#stat_compare_means(label.y = 110)