rm(list = ls())
#Load in libraries
library(readr)
library(lme4)
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 3.5.3
library(car)
## Warning: package 'car' was built under R version 3.5.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.5.2
library(emmeans)
## Warning: package 'emmeans' was built under R version 3.5.3
setwd("c:/users/Paul/Documents/Rwork")
Deer_data<- read.csv(file="Deer_data.csv")
Deer_data$SamplingperiodNum = as.numeric(Deer_data$Samplingperiod)
str(Deer_data)
## 'data.frame': 80 obs. of 21 variables:
## $ Sampling.season : Factor w/ 4 levels "fall2018","Spring2018",..: 2 2 2 2 2 4 4 4 4 4 ...
## $ Site : Factor w/ 5 levels "Bishop","College",..: 3 4 1 5 2 3 4 1 5 2 ...
## $ Samplingperiod : Factor w/ 4 levels "S1","S2","S3",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ Deer : Factor w/ 2 levels "Con","Excl": 1 1 1 1 1 1 1 1 1 1 ...
## $ Honeysuckle : Factor w/ 2 levels "H ","NH": 1 1 1 1 1 1 1 1 1 1 ...
## $ InorgN : num 21.1 16.3 26.7 19.7 47.4 ...
## $ NitrateN : num 22.4 16.8 27.9 20.7 34 ...
## $ Nmin : num 1.24 0.96 1.57 1.16 2.79 1.2 0.44 0.39 0.87 1.08 ...
## $ NetNitr : num 1.32 0.99 1.64 1.22 2 1.17 0.54 0.56 0.99 0.75 ...
## $ SMR : num 1.6 1.85 1.47 1.63 2.18 0.93 1.26 0.91 1.05 1.4 ...
## $ Nitr_index : num 1.06 1.03 1.04 1.05 0.72 0.98 1.23 1.44 1.14 0.69 ...
## $ sqrNetNitr : num 1.15 0.99 1.28 1.1 1.41 1.08 0.73 0.75 0.99 0.87 ...
## $ lgNmin : num 0.09 -0.02 0.2 0.06 0.45 0.08 -0.36 -0.41 -0.06 0.03 ...
## $ lgSMR : num 0.2 0.27 0.17 0.21 0.34 -0.03 0.1 -0.04 0.02 0.15 ...
## $ sqrSMR : num 1.26 1.36 1.21 1.28 1.48 0.96 1.12 0.95 1.02 1.18 ...
## $ Nmin_OM : num 29.9 22.9 36.6 27.7 66.1 ...
## $ NetNitr_OM : num 32 23.6 38.4 29.2 47.4 ...
## $ NetNitr_OM_index : num 1.07 1.03 1.05 1.06 0.72 0.98 1.24 1.43 1.14 0.69 ...
## $ smr_OM : num 39 43.5 34.2 39 51.7 ...
## $ lgSMR_OM : num 1.59 1.64 1.53 1.59 1.71 1.66 1.8 1.62 1.69 1.87 ...
## $ SamplingperiodNum: num 1 1 1 1 1 2 2 2 2 2 ...
##Col names:
# Site - site
#Season of sampling - sampling time of the year
#Samplingperiod - Number of sampling
# Deer - Acc (Control) and Ex (Exclosure)
# Honeysuckle - H (control) and NH (Removal)
# # The rest are response variables
#using the time variable as continuous v. categorical/discrete
##########################
#InorgN model
#Samplingperiod as a continuous variable
###This will give you a trend line through time (e.g. classic regression)
repeated.split.plot.model1 <- lmer(InorgN ~ Deer * Honeysuckle * SamplingperiodNum +
(SamplingperiodNum | Site/Deer/Honeysuckle),
data = Deer_data,
contrasts = list(Deer = contr.sum,
Honeysuckle = contr.sum),
REML = TRUE)
# This gives you the correct type 3 Anova; interpret P-values from here!
# Must use the capital "A" Anova function from the "car" package
Anova(repeated.split.plot.model1, type = 3, test.statistic = "F")
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: InorgN
## F Df Df.res Pr(>F)
## (Intercept) 41.9123 1 4 0.002933 **
## Deer 1.1583 1 4 0.342408
## Honeysuckle 0.0192 1 8 0.893128
## SamplingperiodNum 6.5561 1 4 0.062606 .
## Deer:Honeysuckle 0.0797 1 8 0.784900
## Deer:SamplingperiodNum 0.2140 1 4 0.667661
## Honeysuckle:SamplingperiodNum 0.2082 1 8 0.660322
## Deer:Honeysuckle:SamplingperiodNum 0.2634 1 8 0.621680
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#pairwise comparisons of deer*honeysuckle interaction
emmeans(repeated.split.plot.model1, pairwise ~ Deer * Honeysuckle)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
## Deer Honeysuckle emmean SE df lower.CL upper.CL
## Con H 17.7 1.7 12.9 14.0 21.4
## Excl H 14.6 1.7 12.9 11.0 18.3
## Con NH 15.0 1.7 12.9 11.4 18.7
## Excl NH 13.3 1.7 12.9 9.6 16.9
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Con,H - Excl,H 3.080 2.03 10.7 1.517 0.4614
## Con,H - Con,NH 2.667 2.03 8.0 1.313 0.5804
## Con,H - Excl,NH 4.438 2.03 10.7 2.186 0.1887
## Excl,H - Con,NH -0.414 2.03 10.7 -0.204 0.9968
## Excl,H - Excl,NH 1.357 2.03 8.0 0.669 0.9060
## Con,NH - Excl,NH 1.772 2.03 10.7 0.872 0.8188
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
##NitrateN model
#Samplingperiod as a continuous variable
###This will give you a trend line through time (e.g. classic regression)
repeated.split.plot.model2 <- lmer(NitrateN ~ Deer * Honeysuckle * SamplingperiodNum +
(SamplingperiodNum | Site/Deer/Honeysuckle),
data = Deer_data,
contrasts = list(Deer = contr.sum,
Honeysuckle = contr.sum),
REML = TRUE)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues
# This gives you the correct type 3 Anova; interpret P-values from here!
# Must use the capital "A" Anova function from the "car" package
Anova(repeated.split.plot.model2, type = 3, test.statistic = "F")
## Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
##
## Response: NitrateN
## F Df Df.res Pr(>F)
## (Intercept) 70.1243 1 4 0.001112 **
## Deer 2.0248 1 4 0.227835
## Honeysuckle 0.2740 1 8 0.614878
## SamplingperiodNum 8.4543 1 4 0.043782 *
## Deer:Honeysuckle 1.1789 1 8 0.309203
## Deer:SamplingperiodNum 0.2978 1 4 0.614253
## Honeysuckle:SamplingperiodNum 2.2715 1 8 0.170203
## Deer:Honeysuckle:SamplingperiodNum 1.1976 1 8 0.305658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#pairwise comparisons of deer*honeysuckle interaction
emmeans(repeated.split.plot.model2, pairwise ~ Deer * Honeysuckle)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
## Deer Honeysuckle emmean SE df lower.CL upper.CL
## Con H 18.2 1.67 12.5 14.62 21.9
## Excl H 16.0 1.67 12.5 12.40 19.7
## Con NH 16.0 1.67 12.5 12.41 19.7
## Excl NH 13.3 1.67 12.5 9.69 16.9
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Con,H - Excl,H 2.2120 2.22 7.21 0.998 0.7553
## Con,H - Con,NH 2.2095 1.67 8.00 1.326 0.5733
## Con,H - Excl,NH 4.9225 2.22 7.21 2.221 0.2046
## Excl,H - Con,NH -0.0025 2.22 7.21 -0.001 1.0000
## Excl,H - Excl,NH 2.7105 1.67 8.00 1.627 0.4170
## Con,NH - Excl,NH 2.7130 2.22 7.21 1.224 0.6319
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates