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")
seasonal_data<- read.csv(file="seasonal_data.csv")
seasonal_data$SeasonNum = as.numeric(seasonal_data$Season)
str(seasonal_data)
## 'data.frame': 80 obs. of 9 variables:
## $ Season : int 64 64 64 64 64 183 183 183 183 183 ...
## $ Time : Factor w/ 4 levels "S1","S2","S3",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ Site : Factor w/ 5 levels "Bachelor","College",..: 3 4 1 5 2 3 4 1 5 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 ...
## $ Treatment : Factor w/ 4 levels "ConH","ConNH",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ min_OM : num 29.9 22.9 36.6 27.7 66.1 ...
## $ NetNitr_OM : num 32 23.6 38.4 29.2 47.4 ...
## $ SeasonNum : num 64 64 64 64 64 183 183 183 183 183 ...
##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
##########################
#Nmin model
#Season as a continuous variable
###This will give you a trend line through time (e.g. classic regression)
repeated.split.plot.model1 <- lmer(min_OM ~ Treatment* SeasonNum + (SeasonNum | Site/Deer/Honeysuckle),
data = seasonal_data,
contrasts = list(Deer = contr.sum,
Honeysuckle = contr.sum), na.action=na.omit,
REML = TRUE)
## Warning in model.matrix.default(fixedform, fr, contrasts): variable 'Deer'
## is absent, its contrast will be ignored
## Warning in model.matrix.default(fixedform, fr, contrasts): variable
## 'Honeysuckle' is absent, its contrast will be ignored
## Warning in optwrap(optimizer, devfun, getStart(start, rho$lower, rho$pp), :
## convergence code 1 from bobyqa: bobyqa -- maximum number of function
## evaluations exceeded
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 2.24112 (tol =
## 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
summary(repeated.split.plot.model1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## min_OM ~ Treatment * SeasonNum + (SeasonNum | Site/Deer/Honeysuckle)
## Data: seasonal_data
##
## REML criterion at convergence: 560.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4301 -0.3754 0.1121 0.4706 3.6772
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Honeysuckle:(Deer:Site) (Intercept) 2.434e+01 4.93328
## SeasonNum 4.258e-04 0.02064 -1.00
## Deer:Site (Intercept) 1.492e+01 3.86206
## SeasonNum 2.562e-04 0.01601 -1.00
## Site (Intercept) 2.294e+01 4.78991
## SeasonNum 3.261e-04 0.01806 -1.00
## Residual 6.013e+01 7.75450
## Number of obs: 79, groups:
## Honeysuckle:(Deer:Site), 20; Deer:Site, 10; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 37.78549 5.28379 7.151
## TreatmentConNH -8.66716 6.42252 -1.349
## TreatmentExcH -12.08342 6.83079 -1.769
## TreatmentExcNH -14.59860 6.83263 -2.137
## SeasonNum -0.12626 0.02524 -5.002
## TreatmentConNH:SeasonNum 0.03656 0.03233 1.131
## TreatmentExcH:SeasonNum 0.05233 0.03382 1.547
## TreatmentExcNH:SeasonNum 0.05630 0.03385 1.663
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH TrtENH SesnNm TCNH:S TEH:SN
## TretmntCnNH -0.600
## TretmntExcH -0.646 0.464
## TrtmntExcNH -0.646 0.464 0.564
## SeasonNum -0.928 0.580 0.617 0.617
## TrtmnCNH:SN 0.551 -0.913 -0.426 -0.426 -0.638
## TrtmntEH:SN 0.596 -0.433 -0.921 -0.514 -0.670 0.476
## TrtmnENH:SN 0.595 -0.433 -0.514 -0.921 -0.669 0.476 0.544
## convergence code: 1
## Model failed to converge with max|grad| = 2.24112 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
# 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: min_OM
## F Df Df.res Pr(>F)
## (Intercept) 51.1397 1 14.0206 4.885e-06 ***
## Treatment 1.5161 3 7.5089 0.2872357
## SeasonNum 25.0242 1 15.0615 0.0001557 ***
## Treatment:SeasonNum 0.9940 3 7.3281 0.4473099
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#pairwise comparisons of deer*honeysuckle interaction
emmeans(repeated.split.plot.model1, pairwise ~ Treatment* SeasonNum)
## $emmeans
## Treatment SeasonNum emmean SE df lower.CL upper.CL
## ConH 170 16.4 2.07 14.5 11.94 20.8
## ConNH 170 13.9 2.12 15.5 9.40 18.4
## ExcH 170 13.2 2.07 14.5 8.74 17.6
## ExcNH 170 11.3 2.14 14.2 6.74 15.9
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## ConH,169.658227848101 - ConNH,169.658227848101 2.464 2.65 8.10
## ConH,169.658227848101 - ExcH,169.658227848101 3.205 2.71 9.85
## ConH,169.658227848101 - ExcNH,169.658227848101 5.047 2.76 9.75
## ConNH,169.658227848101 - ExcH,169.658227848101 0.741 2.75 10.26
## ConNH,169.658227848101 - ExcNH,169.658227848101 2.583 2.80 10.16
## ExcH,169.658227848101 - ExcNH,169.658227848101 1.842 2.67 7.68
## t.ratio p.value
## 0.929 0.7910
## 1.181 0.6513
## 1.825 0.3186
## 0.269 0.9927
## 0.922 0.7942
## 0.690 0.8979
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
##NetNitr model
#Season as a continuous variable
###This will give you a trend line through time (e.g. classic regression)
repeated.split.plot.model2 <- lmer(NetNitr_OM ~ Treatment* SeasonNum * SeasonNum + (SeasonNum | Site/Deer/Honeysuckle),
data = seasonal_data,
contrasts = list(Deer = contr.sum,
Honeysuckle = contr.sum), na.action=na.omit,
REML = TRUE)
## Warning in model.matrix.default(fixedform, fr, contrasts): variable 'Deer'
## is absent, its contrast will be ignored
## Warning in model.matrix.default(fixedform, fr, contrasts): variable
## 'Honeysuckle' is absent, its contrast will be ignored
## Warning in optwrap(optimizer, devfun, getStart(start, rho$lower, rho$pp), :
## convergence code 1 from bobyqa: bobyqa -- maximum number of function
## evaluations exceeded
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 2.00316 (tol =
## 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
# 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: NetNitr_OM
## F Df Df.res Pr(>F)
## (Intercept) 87.0619 1 14.9904 1.240e-07 ***
## Treatment 2.5397 3 7.4794 0.1347
## SeasonNum 34.1113 1 15.8049 2.635e-05 ***
## Treatment:SeasonNum 1.6636 3 7.2425 0.2579
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#pairwise comparisons of deer*honeysuckle interaction
emmeans(repeated.split.plot.model2, pairwise ~ Treatment* SeasonNum)
## $emmeans
## Treatment SeasonNum emmean SE df lower.CL upper.CL
## ConH 170 16.5 1.67 14.8 12.95 20.1
## ConNH 170 15.1 1.71 15.9 11.52 18.8
## ExcH 170 14.0 1.67 14.8 10.45 17.6
## ExcNH 170 11.4 1.73 14.2 7.65 15.1
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df
## ConH,169.658227848101 - ConNH,169.658227848101 1.36 2.18 8.11
## ConH,169.658227848101 - ExcH,169.658227848101 2.51 2.22 9.89
## ConH,169.658227848101 - ExcNH,169.658227848101 5.15 2.27 9.75
## ConNH,169.658227848101 - ExcH,169.658227848101 1.15 2.25 10.34
## ConNH,169.658227848101 - ExcNH,169.658227848101 3.79 2.30 10.19
## ExcH,169.658227848101 - ExcNH,169.658227848101 2.65 2.20 7.63
## t.ratio p.value
## 0.624 0.9216
## 1.131 0.6799
## 2.275 0.1703
## 0.509 0.9551
## 1.649 0.3961
## 1.205 0.6419
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates