rm(list = ls())
setwd("c:/users/Paul/Documents/Rwork")
getwd()
## [1] "c:/users/Paul/Documents/Rwork"
library(readr)
library(qqplotr)
## Warning: package 'qqplotr' was built under R version 3.5.3
## Loading required package: ggplot2
##
## Attaching package: 'qqplotr'
## The following objects are masked from 'package:ggplot2':
##
## stat_qq_line, StatQqLine
library(lavaan)
## This is lavaan 0.6-3
## lavaan is BETA software! Please report any bugs.
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(carData)
#Load the Deer_projectdata
Deer_projectdata<- read.csv(file="Deer_project.csv")
str(Deer_projectdata)
## 'data.frame': 80 obs. of 23 variables:
## $ Time : Factor w/ 4 levels "fall2018","Spring2018",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Site : Factor w/ 5 levels "Bishop","College",..: 3 3 3 3 4 4 4 4 1 1 ...
## $ Deer : Factor w/ 2 levels "Con","Excl": 1 1 2 2 1 1 2 2 1 1 ...
## $ Honeysuckle : Factor w/ 2 levels "H ","NH": 2 1 2 1 2 1 2 1 2 1 ...
## $ Nmin : num 1.39 1.24 1.15 0.85 0.84 0.96 0.67 0.88 1.3 1.57 ...
## $ Deer.effect : num 0.23 0.39 NA NA 0.18 0.09 NA NA 0.25 -0.29 ...
## $ Honeysuckle.effect : num -0.14 NA -0.3 NA 0.12 NA 0.21 NA 0.27 NA ...
## $ NetNitr : num 1.32 1.32 1.17 1.16 0.88 0.99 0.76 0.98 1.34 1.64 ...
## $ Deer.effect.1 : num 0.15 0.16 NA NA 0.12 0.01 NA NA 0.23 0.01 ...
## $ Honeysuckle.effect.1: num 0 NA -0.01 NA 0.11 NA 0.23 NA 0.31 NA ...
## $ SMR : num 1.63 1.6 1.55 1.27 1.61 1.85 1.57 1.78 1.44 1.47 ...
## $ Deer.effect.2 : num 0.07 0.33 NA NA 0.04 0.07 NA NA -0.47 -0.02 ...
## $ Honeysuckle.effect.2: num -0.03 NA -0.28 NA 0.24 NA 0.21 NA 0.04 NA ...
## $ Nitr_index : num 0.95 1.06 1.02 1.36 1.05 1.03 1.13 1.11 1.03 1.04 ...
## $ sqrNetNitr : num 1.15 1.15 1.08 1.08 0.94 0.99 0.87 0.99 1.16 1.28 ...
## $ lgNmin : num 0.14 0.09 0.06 -0.07 -0.08 -0.02 -0.17 -0.06 0.11 0.2 ...
## $ lgSMR : num 0.21 0.2 0.19 0.1 0.21 0.27 0.2 0.25 0.16 0.17 ...
## $ sqrSMR : num 1.28 1.26 1.24 1.13 1.27 1.36 1.25 1.33 1.2 1.21 ...
## $ Nmin_OM : num 34 29.9 27.2 17.6 19.1 ...
## $ NetNitr_OM : num 32.4 32 27.8 23.4 20.2 ...
## $ NetNitr_OM_index : num 0.95 1.07 1.02 1.33 1.06 1.03 1.14 1.12 1.02 1.05 ...
## $ SMR_OM : num 40 39 37 25.9 37.1 ...
## $ lgSMR_OM : num 1.6 1.59 1.57 1.41 1.57 1.64 1.57 1.61 1.52 1.53 ...
head(Deer_projectdata)
## Time Site Deer Honeysuckle Nmin Deer.effect
## 1 Spring2018 Kramer Con NH 1.39 0.23
## 2 Spring2018 Kramer Con H 1.24 0.39
## 3 Spring2018 Kramer Excl NH 1.15 NA
## 4 Spring2018 Kramer Excl H 0.85 NA
## 5 Spring2018 Rhinehart Con NH 0.84 0.18
## 6 Spring2018 Rhinehart Con H 0.96 0.09
## Honeysuckle.effect NetNitr Deer.effect.1 Honeysuckle.effect.1 SMR
## 1 -0.14 1.32 0.15 0.00 1.63
## 2 NA 1.32 0.16 NA 1.60
## 3 -0.30 1.17 NA -0.01 1.55
## 4 NA 1.16 NA NA 1.27
## 5 0.12 0.88 0.12 0.11 1.61
## 6 NA 0.99 0.01 NA 1.85
## Deer.effect.2 Honeysuckle.effect.2 Nitr_index sqrNetNitr lgNmin lgSMR
## 1 0.07 -0.03 0.95 1.15 0.14 0.21
## 2 0.33 NA 1.06 1.15 0.09 0.20
## 3 NA -0.28 1.02 1.08 0.06 0.19
## 4 NA NA 1.36 1.08 -0.07 0.10
## 5 0.04 0.24 1.05 0.94 -0.08 0.21
## 6 0.07 NA 1.03 0.99 -0.02 0.27
## sqrSMR Nmin_OM NetNitr_OM NetNitr_OM_index SMR_OM lgSMR_OM
## 1 1.28 34.03 32.38 0.95 40.02 1.60
## 2 1.26 29.87 32.02 1.07 39.01 1.59
## 3 1.24 27.21 27.75 1.02 36.96 1.57
## 4 1.13 17.60 23.44 1.33 25.88 1.41
## 5 1.27 19.06 20.22 1.06 37.09 1.57
## 6 1.36 22.86 23.60 1.03 43.45 1.64
#Checking for correlation between SMR_OM and Nmin_OM
cor<-cor.test(Deer_projectdata$SMR_OM,Deer_projectdata$Nmin_OM)
cor
##
## Pearson's product-moment correlation
##
## data: Deer_projectdata$SMR_OM and Deer_projectdata$Nmin_OM
## t = 4.4198, df = 78, p-value = 3.153e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2526577 0.6075110
## sample estimates:
## cor
## 0.4475308
#Checking to make sure that variances across samples are equal (Nmin_OM, lgSMR_OM, NetNitr_OM)
bartlett.test(Deer_projectdata$NetNitr_OM,Deer_projectdata$Deer)
##
## Bartlett test of homogeneity of variances
##
## data: Deer_projectdata$NetNitr_OM and Deer_projectdata$Deer
## Bartlett's K-squared = 0.73827, df = 1, p-value = 0.3902
bartlett.test(Deer_projectdata$Nmin_OM,Deer_projectdata$Deer)
##
## Bartlett test of homogeneity of variances
##
## data: Deer_projectdata$Nmin_OM and Deer_projectdata$Deer
## Bartlett's K-squared = 0.58108, df = 1, p-value = 0.4459
bartlett.test(Deer_projectdata$lgSMR_OM,Deer_projectdata$Deer)
##
## Bartlett test of homogeneity of variances
##
## data: Deer_projectdata$lgSMR_OM and Deer_projectdata$Deer
## Bartlett's K-squared = 2.4134, df = 1, p-value = 0.1203
#The outputs ofn the above codes show that variances across samples are equal (because all p-values>0.05)
#Normality check for NetNitr_OM
hist(Deer_projectdata$NetNitr_OM)
#Run a lm on NetNitr_OM
model1<-lm(NetNitr_OM~Deer*Honeysuckle*Time, data=Deer_projectdata,na.action=na.omit)
summary(model1)
##
## Call:
## lm(formula = NetNitr_OM ~ Deer * Honeysuckle * Time, data = Deer_projectdata,
## na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.272 -4.274 0.052 4.610 23.818
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.400 3.957 6.166 5.54e-08
## DeerExcl 20.472 5.597 3.658 0.000521
## HoneysuckleNH 1.292 5.597 0.231 0.818175
## TimeSpring2018 9.718 5.597 1.736 0.087378
## Timespring2019 -10.360 5.597 -1.851 0.068839
## Timesummer2018 14.978 5.597 2.676 0.009477
## DeerExcl:HoneysuckleNH -22.304 7.915 -2.818 0.006448
## DeerExcl:TimeSpring2018 -28.176 7.915 -3.560 0.000712
## DeerExcl:Timespring2019 -27.122 7.915 -3.427 0.001081
## DeerExcl:Timesummer2018 -23.144 7.915 -2.924 0.004794
## HoneysuckleNH:TimeSpring2018 -7.024 7.915 -0.887 0.378209
## HoneysuckleNH:Timespring2019 -3.587 8.158 -0.440 0.661679
## HoneysuckleNH:Timesummer2018 3.178 7.915 0.402 0.689390
## DeerExcl:HoneysuckleNH:TimeSpring2018 25.324 11.193 2.262 0.027127
## DeerExcl:HoneysuckleNH:Timespring2019 22.313 11.367 1.963 0.054063
## DeerExcl:HoneysuckleNH:Timesummer2018 14.636 11.193 1.308 0.195766
##
## (Intercept) ***
## DeerExcl ***
## HoneysuckleNH
## TimeSpring2018 .
## Timespring2019 .
## Timesummer2018 **
## DeerExcl:HoneysuckleNH **
## DeerExcl:TimeSpring2018 ***
## DeerExcl:Timespring2019 **
## DeerExcl:Timesummer2018 **
## HoneysuckleNH:TimeSpring2018
## HoneysuckleNH:Timespring2019
## HoneysuckleNH:Timesummer2018
## DeerExcl:HoneysuckleNH:TimeSpring2018 *
## DeerExcl:HoneysuckleNH:Timespring2019 .
## DeerExcl:HoneysuckleNH:Timesummer2018
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.849 on 63 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.6895, Adjusted R-squared: 0.6156
## F-statistic: 9.328 on 15 and 63 DF, p-value: 5.891e-11
plot(model1)
anova(model1)
## Analysis of Variance Table
##
## Response: NetNitr_OM
## Df Sum Sq Mean Sq F value Pr(>F)
## Deer 1 167.9 167.85 2.1436 0.14814
## Honeysuckle 1 251.7 251.74 3.2149 0.07777 .
## Time 3 8611.7 2870.57 36.6592 7.847e-14 ***
## Deer:Honeysuckle 1 219.4 219.42 2.8021 0.09910 .
## Deer:Time 3 932.2 310.74 3.9684 0.01178 *
## Honeysuckle:Time 3 296.5 98.82 1.2621 0.29506
## Deer:Honeysuckle:Time 3 476.5 158.83 2.0284 0.11892
## Residuals 63 4933.2 78.30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
resid1<-resid(model1)
hist(resid1)
#Run linear mixed effects model setting "Site and Time" as random effects
#Also compare the AICs
library(lme4)
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 3.5.3
model1b<-lmer(NetNitr_OM~Deer +Honeysuckle + (1|Site),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model1b)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: NetNitr_OM ~ Deer + Honeysuckle + (1 | Site)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 651.1 662.9 -320.5 641.1 74
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.80298 -0.73568 -0.03644 0.44561 2.98069
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.0 0.00
## Residual 195.8 13.99
## Number of obs: 79, groups: Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 29.850 2.716 10.992
## DeerExcl -2.870 3.149 -0.911
## HoneysuckleNH -3.571 3.149 -1.134
##
## Correlation of Fixed Effects:
## (Intr) DrExcl
## DeerExcl -0.580
## HoneyscklNH -0.565 -0.013
resid1b<-resid(model1b)
hist(resid1b)
model1ba<-lmer(NetNitr_OM~Deer +Honeysuckle + (1|Site/Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model1ba)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: NetNitr_OM ~ Deer + Honeysuckle + (1 | Site/Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 626.1 640.3 -307.1 614.1 73
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7379 -0.5527 -0.1016 0.4270 3.8280
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time:Site (Intercept) 109.8 10.477
## Site (Intercept) 0.0 0.000
## Residual 88.9 9.429
## Number of obs: 79, groups: Time:Site, 20; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 29.624 2.973 9.963
## DeerExcl -2.419 2.126 -1.138
## HoneysuckleNH -4.021 2.126 -1.892
##
## Correlation of Fixed Effects:
## (Intr) DrExcl
## DeerExcl -0.357
## HoneyscklNH -0.346 -0.016
resid1ba<-resid(model1ba)
hist(resid1ba)
model1bb<-lmer(NetNitr_OM~Deer +Honeysuckle + (1|Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model1bb)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: NetNitr_OM ~ Deer + Honeysuckle + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 603.6 615.5 -296.8 593.6 74
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5905 -0.5861 -0.1212 0.3883 4.0101
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 106.56 10.323
## Residual 91.44 9.562
## Number of obs: 79, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 29.638 5.485 5.403
## DeerExcl -2.447 2.153 -1.137
## HoneysuckleNH -3.993 2.153 -1.855
##
## Correlation of Fixed Effects:
## (Intr) DrExcl
## DeerExcl -0.196
## HoneyscklNH -0.191 -0.013
resid1bb<-resid(model1bb)
hist(resid1bb)
model1bc<-lmer(NetNitr_OM~Deer +Honeysuckle + Time + (1|Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model1bc)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: NetNitr_OM ~ Deer + Honeysuckle + Time + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 592.8 611.8 -288.4 576.8 71
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6851 -0.5937 -0.1539 0.3874 4.0992
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 0.00 0.000
## Residual 86.81 9.317
## Number of obs: 79, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 32.926 2.552 12.904
## DeerExcl -2.428 2.097 -1.158
## HoneysuckleNH -4.012 2.097 -1.913
## TimeSpring2018 -1.551 2.946 -0.526
## Timespring2019 -20.292 2.986 -6.796
## Timesummer2018 8.654 2.946 2.937
##
## Correlation of Fixed Effects:
## (Intr) DrExcl HnysNH TS2018 Tm2019
## DeerExcl -0.406
## HoneyscklNH -0.406 -0.013
## TmSprng2018 -0.577 0.000 0.000
## Tmsprng2019 -0.570 -0.019 0.019 0.493
## Timsmmr2018 -0.577 0.000 0.000 0.500 0.493
resid1bc<-resid(model1bc)
hist(resid1bc)
#split-plot model
library(car)
library(emmeans)
## Warning: package 'emmeans' was built under R version 3.5.3
## Welcome to emmeans.
## NOTE -- Important change from versions <= 1.41:
## Indicator predictors are now treated as 2-level factors by default.
## To revert to old behavior, use emm_options(cov.keep = character(0))
# This sets up your split-plot model
# Deer represents Exclosure or Control/Access
# Honeysuckle represents Removal or Control/Retention
# Site is the site (Bachelor, Reinhart, Western, College, Kramer)
# The random effect is needed to correctly account for degrees of freedom
split.plot.model1c <- lmer(NetNitr_OM ~ Deer * Honeysuckle +
(1 | Site/Time),
data = Deer_projectdata,
contrasts = list(Deer = contr.sum,
Honeysuckle = contr.sum))
summary(split.plot.model1c)
## Linear mixed model fit by REML ['lmerMod']
## Formula: NetNitr_OM ~ Deer * Honeysuckle + (1 | Site/Time)
## Data: Deer_projectdata
##
## REML criterion at convergence: 602
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5555 -0.5803 -0.1834 0.4112 3.6301
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time:Site (Intercept) 115.34 10.740
## Site (Intercept) 0.00 0.000
## Residual 90.01 9.487
## Number of obs: 79, groups: Time:Site, 20; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.430 2.629 10.054
## Deer1 1.235 1.070 1.155
## Honeysuckle1 1.985 1.070 1.856
## Deer1:Honeysuckle1 -1.666 1.070 -1.558
##
## Correlation of Fixed Effects:
## (Intr) Deer1 Hnysc1
## Deer1 0.007
## Honeysuckl1 -0.007 -0.016
## Dr1:Hnysck1 -0.007 -0.016 0.016
# This gives you the correct type 3 Anova; interpret P-values from here!
anova(split.plot.model1c, type = 3, test.statistic = "F")
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 121.11 121.11 1.3455
## Honeysuckle 1 318.75 318.75 3.5412
## Deer:Honeysuckle 1 218.38 218.38 2.4261
# This allows you to do post-hoc comparisons;
# specifically this line give pairwise comparisons of deer*honeysuckle interaction
emmeans(split.plot.model1c, pairwise ~ Deer * Honeysuckle)
## $emmeans
## Deer Honeysuckle emmean SE df lower.CL upper.CL
## Con H 28.0 3.20 8.68 20.7 35.3
## Excl H 28.8 3.20 8.68 21.6 36.1
## Con NH 27.3 3.25 9.16 20.0 34.7
## Excl NH 21.5 3.20 8.68 14.3 28.8
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Con,H - Excl,H -0.862 3.00 56.0 -0.287 0.9917
## Con,H - Con,NH 0.638 3.05 56.2 0.209 0.9967
## Con,H - Excl,NH 6.441 3.00 56.0 2.147 0.1512
## Excl,H - Con,NH 1.500 3.05 56.2 0.492 0.9606
## Excl,H - Excl,NH 7.302 3.00 56.0 2.434 0.0823
## Con,NH - Excl,NH 5.802 3.05 56.2 1.902 0.2390
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
###########################################
##############################
#Run a lm on Nmin_OM
model2<-lm(Nmin_OM~Deer*Honeysuckle*Time, data=Deer_projectdata,na.action=na.omit)
summary(model2)
##
## Call:
## lm(formula = Nmin_OM ~ Deer * Honeysuckle * Time, data = Deer_projectdata,
## na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.282 -5.568 -0.324 5.838 29.470
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.314 4.420 4.143 0.000102
## DeerExcl 13.302 6.251 2.128 0.037203
## HoneysuckleNH 2.686 6.251 0.430 0.668876
## TimeSpring2018 18.316 6.251 2.930 0.004694
## Timespring2019 -5.720 6.251 -0.915 0.363618
## Timesummer2018 21.128 6.251 3.380 0.001240
## DeerExcl:HoneysuckleNH -12.848 8.841 -1.453 0.151028
## DeerExcl:TimeSpring2018 -22.200 8.841 -2.511 0.014570
## DeerExcl:Timespring2019 -18.982 8.841 -2.147 0.035577
## DeerExcl:Timesummer2018 -18.372 8.841 -2.078 0.041712
## HoneysuckleNH:TimeSpring2018 -11.552 8.841 -1.307 0.195992
## HoneysuckleNH:Timespring2019 -7.914 8.841 -0.895 0.374041
## HoneysuckleNH:Timesummer2018 -5.278 8.841 -0.597 0.552602
## DeerExcl:HoneysuckleNH:TimeSpring2018 18.978 12.503 1.518 0.133955
## DeerExcl:HoneysuckleNH:Timespring2019 15.672 12.503 1.254 0.214581
## DeerExcl:HoneysuckleNH:Timesummer2018 13.662 12.503 1.093 0.278604
##
## (Intercept) ***
## DeerExcl *
## HoneysuckleNH
## TimeSpring2018 **
## Timespring2019
## Timesummer2018 **
## DeerExcl:HoneysuckleNH
## DeerExcl:TimeSpring2018 *
## DeerExcl:Timespring2019 *
## DeerExcl:Timesummer2018 *
## HoneysuckleNH:TimeSpring2018
## HoneysuckleNH:Timespring2019
## HoneysuckleNH:Timesummer2018
## DeerExcl:HoneysuckleNH:TimeSpring2018
## DeerExcl:HoneysuckleNH:Timespring2019
## DeerExcl:HoneysuckleNH:Timesummer2018
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.884 on 64 degrees of freedom
## Multiple R-squared: 0.6103, Adjusted R-squared: 0.5189
## F-statistic: 6.681 on 15 and 64 DF, p-value: 2.444e-08
resid2<-resid(model2)
hist(resid2)
#Run linear mixed effects model setting "Site and Time" as random effects
#Also compare the AICs
library(lme4)
model2b<-lmer(Nmin_OM~Deer +Honeysuckle + (1|Site),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model2b)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Nmin_OM ~ Deer + Honeysuckle + (1 | Site)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 659.2 671.1 -324.6 649.2 75
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9419 -0.6379 -0.1030 0.5431 2.8089
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.9964 0.9982
## Residual 194.7906 13.9567
## Number of obs: 80, groups: Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.938 2.739 9.834
## DeerExcl -1.971 3.121 -0.632
## HoneysuckleNH -3.885 3.121 -1.245
##
## Correlation of Fixed Effects:
## (Intr) DrExcl
## DeerExcl -0.570
## HoneyscklNH -0.570 0.000
resid2b<-resid(model2b)
hist(resid2b)
model2ba<-lmer(Nmin_OM~Deer +Honeysuckle + (1|Site/Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model2ba)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Nmin_OM ~ Deer + Honeysuckle + (1 | Site/Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 626.5 640.8 -307.2 614.5 74
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.13227 -0.60642 -0.09609 0.36941 2.92938
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time:Site (Intercept) 118.0 10.86
## Site (Intercept) 0.0 0.00
## Residual 77.8 8.82
## Number of obs: 80, groups: Time:Site, 20; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.938 2.969 9.072
## DeerExcl -1.971 1.972 -1.000
## HoneysuckleNH -3.885 1.972 -1.970
##
## Correlation of Fixed Effects:
## (Intr) DrExcl
## DeerExcl -0.332
## HoneyscklNH -0.332 0.000
resid2ba<-resid(model2ba)
hist(resid2ba)
model2bb<-lmer(Nmin_OM~Deer +Honeysuckle + (1|Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model2bb)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Nmin_OM ~ Deer + Honeysuckle + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 612.3 624.2 -301.2 602.3 75
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0797 -0.5790 -0.2108 0.5946 3.5355
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 102.65 10.132
## Residual 93.14 9.651
## Number of obs: 80, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.937 5.400 4.989
## DeerExcl -1.971 2.158 -0.914
## HoneysuckleNH -3.885 2.158 -1.800
##
## Correlation of Fixed Effects:
## (Intr) DrExcl
## DeerExcl -0.200
## HoneyscklNH -0.200 0.000
resid2bb<-resid(model2bb)
hist(resid2bb)
model2bc<-lmer(Nmin_OM~Deer +Honeysuckle + Time + (1|Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model2bc)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Nmin_OM ~ Deer + Honeysuckle + Time + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 601.7 620.7 -292.8 585.7 72
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1882 -0.5969 -0.1507 0.6068 3.6030
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 0.00 0.000
## Residual 88.48 9.406
## Number of obs: 80, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.024 2.576 10.103
## DeerExcl -1.971 2.103 -0.937
## HoneysuckleNH -3.885 2.103 -1.847
## TimeSpring2018 6.184 2.975 2.079
## Timespring2019 -15.250 2.975 -5.127
## Timesummer2018 12.718 2.975 4.276
##
## Correlation of Fixed Effects:
## (Intr) DrExcl HnysNH TS2018 Tm2019
## DeerExcl -0.408
## HoneyscklNH -0.408 0.000
## TmSprng2018 -0.577 0.000 0.000
## Tmsprng2019 -0.577 0.000 0.000 0.500
## Timsmmr2018 -0.577 0.000 0.000 0.500 0.500
resid2bc<-resid(model2bc)
hist(resid2bc)
#split-plot model
library(car)
library(emmeans)
# This sets up your split-plot model
# Deer represents Exclosure or Control/Access
# Honeysuckle represents Removal or Control/Retention
# Site is the site (Bachelor, Reinhart, Western, College, Kramer)
# The random effect is needed to correctly account for degrees of freedom
split.plot.model2c <- lmer(Nmin_OM ~ Deer * Honeysuckle +
(1 | Site/Time),
data = Deer_projectdata,
contrasts = list(Deer = contr.sum,
Honeysuckle = contr.sum))
summary(split.plot.model2c)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Nmin_OM ~ Deer * Honeysuckle + (1 | Site/Time)
## Data: Deer_projectdata
##
## REML criterion at convergence: 605.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0578 -0.5706 -0.1150 0.3815 2.8772
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time:Site (Intercept) 124.21 11.145
## Site (Intercept) 0.00 0.000
## Residual 81.84 9.047
## Number of obs: 80, groups: Time:Site, 20; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 24.0093 2.6895 8.927
## Deer1 0.9857 1.0114 0.975
## Honeysuckle1 1.9425 1.0114 1.921
## Deer1:Honeysuckle1 -0.1925 1.0114 -0.190
##
## Correlation of Fixed Effects:
## (Intr) Deer1 Hnysc1
## Deer1 0.000
## Honeysuckl1 0.000 0.000
## Dr1:Hnysck1 0.000 0.000 0.000
# This gives you the correct type 3 Anova; interpret P-values from here!
anova(split.plot.model2c, type = 3, test.statistic = "F")
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 77.736 77.736 0.9498
## Honeysuckle 1 301.864 301.864 3.6884
## Deer:Honeysuckle 1 2.965 2.965 0.0362
# This allows you to do post-hoc comparisons;
# specifically this line give pairwise comparisons of deer*honeysuckle interaction
emmeans(split.plot.model2c, pairwise ~ Deer * Honeysuckle)
## $emmeans
## Deer Honeysuckle emmean SE df lower.CL upper.CL
## Con H 26.7 3.21 8.01 19.3 34.1
## Excl H 25.2 3.21 8.01 17.8 32.6
## Con NH 23.2 3.21 8.01 15.8 30.6
## Excl NH 20.9 3.21 8.01 13.5 28.3
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Con,H - Excl,H 1.59 2.86 57 0.555 0.9449
## Con,H - Con,NH 3.50 2.86 57 1.223 0.6147
## Con,H - Excl,NH 5.86 2.86 57 2.047 0.1831
## Excl,H - Con,NH 1.91 2.86 57 0.669 0.9084
## Excl,H - Excl,NH 4.27 2.86 57 1.493 0.4486
## Con,NH - Excl,NH 2.36 2.86 57 0.824 0.8430
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
#####################
#Run a lm on lgSMR_OM
model3<-lm(lgSMR_OM~Deer*Honeysuckle*Time, data=Deer_projectdata,na.action=na.omit)
summary(model3)
##
## Call:
## lm(formula = lgSMR_OM ~ Deer * Honeysuckle * Time, data = Deer_projectdata,
## na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4260 -0.0690 0.0040 0.0845 0.5080
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.70400 0.06656 25.599 < 2e-16
## DeerExcl 0.09600 0.09414 1.020 0.31166
## HoneysuckleNH -0.03400 0.09414 -0.361 0.71915
## TimeSpring2018 -0.09200 0.09414 -0.977 0.33209
## Timespring2019 -0.25000 0.09414 -2.656 0.00998
## Timesummer2018 0.02400 0.09414 0.255 0.79958
## DeerExcl:HoneysuckleNH -0.02400 0.13313 -0.180 0.85750
## DeerExcl:TimeSpring2018 -0.15400 0.13313 -1.157 0.25166
## DeerExcl:Timespring2019 -0.20800 0.13313 -1.562 0.12312
## DeerExcl:Timesummer2018 -0.14600 0.13313 -1.097 0.27689
## HoneysuckleNH:TimeSpring2018 -0.01200 0.13313 -0.090 0.92846
## HoneysuckleNH:Timespring2019 0.04600 0.13313 0.346 0.73083
## HoneysuckleNH:Timesummer2018 0.10800 0.13313 0.811 0.42023
## DeerExcl:HoneysuckleNH:TimeSpring2018 0.12200 0.18827 0.648 0.51930
## DeerExcl:HoneysuckleNH:Timespring2019 -0.08400 0.18827 -0.446 0.65698
## DeerExcl:HoneysuckleNH:Timesummer2018 -0.04400 0.18827 -0.234 0.81596
##
## (Intercept) ***
## DeerExcl
## HoneysuckleNH
## TimeSpring2018
## Timespring2019 **
## Timesummer2018
## DeerExcl:HoneysuckleNH
## DeerExcl:TimeSpring2018
## DeerExcl:Timespring2019
## DeerExcl:Timesummer2018
## HoneysuckleNH:TimeSpring2018
## HoneysuckleNH:Timespring2019
## HoneysuckleNH:Timesummer2018
## DeerExcl:HoneysuckleNH:TimeSpring2018
## DeerExcl:HoneysuckleNH:Timespring2019
## DeerExcl:HoneysuckleNH:Timesummer2018
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1488 on 64 degrees of freedom
## Multiple R-squared: 0.5731, Adjusted R-squared: 0.473
## F-statistic: 5.727 on 15 and 64 DF, p-value: 3.069e-07
plot(model3)
anova(model3)
## Analysis of Variance Table
##
## Response: lgSMR_OM
## Df Sum Sq Mean Sq F value Pr(>F)
## Deer 1 0.03828 0.03828 1.7280 0.19336
## Honeysuckle 1 0.00253 0.00253 0.1143 0.73645
## Time 3 1.63387 0.54462 24.5839 1.063e-10 ***
## Deer:Honeysuckle 1 0.00325 0.00325 0.1468 0.70292
## Deer:Time 3 0.17046 0.05682 2.5649 0.06232 .
## Honeysuckle:Time 3 0.02491 0.00830 0.3749 0.77141
## Deer:Honeysuckle:Time 3 0.02983 0.00994 0.4489 0.71894
## Residuals 64 1.41784 0.02215
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
resid3<-resid(model3)
hist(resid3)
#Run linear mixed effects model setting "Site and Time" as random effects
#Also compare the AICs
library(lme4)
model3b<-lmer(SMR_OM~Deer +Honeysuckle + (1|Site),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model3b)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: SMR_OM ~ Deer + Honeysuckle + (1 | Site)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 747.3 759.2 -368.7 737.3 75
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3899 -0.5629 -0.1529 0.3624 6.4047
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 24.09 4.908
## Residual 570.55 23.886
## Number of obs: 80, groups: Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 46.167 5.120 9.017
## DeerExcl -4.060 5.341 -0.760
## HoneysuckleNH 1.320 5.341 0.247
##
## Correlation of Fixed Effects:
## (Intr) DrExcl
## DeerExcl -0.522
## HoneyscklNH -0.522 0.000
resid3b<-resid(model3b)
hist(resid3b)
model3ba<-lmer(SMR_OM~Deer +Honeysuckle + (1|Site/Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model3ba)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: SMR_OM ~ Deer + Honeysuckle + (1 | Site/Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 729.6 743.9 -358.8 717.6 74
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5111 -0.3455 -0.0609 0.1492 6.5740
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time:Site (Intercept) 279.0 16.70
## Site (Intercept) 0.0 0.00
## Residual 315.6 17.77
## Number of obs: 80, groups: Time:Site, 20; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 46.167 5.078 9.091
## DeerExcl -4.060 3.973 -1.022
## HoneysuckleNH 1.320 3.973 0.332
##
## Correlation of Fixed Effects:
## (Intr) DrExcl
## DeerExcl -0.391
## HoneyscklNH -0.391 0.000
resid3ba<-resid(model3ba)
hist(resid3ba)
model3bb<-lmer(SMR_OM~Deer +Honeysuckle + (1|Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model3bb)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: SMR_OM ~ Deer + Honeysuckle + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 731.4 743.3 -360.7 721.4 75
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7182 -0.4435 -0.1158 0.2511 7.0120
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 160.9 12.68
## Residual 433.8 20.83
## Number of obs: 80, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 46.167 7.516 6.143
## DeerExcl -4.060 4.657 -0.872
## HoneysuckleNH 1.320 4.657 0.283
##
## Correlation of Fixed Effects:
## (Intr) DrExcl
## DeerExcl -0.310
## HoneyscklNH -0.310 0.000
resid3bb<-resid(model3bb)
hist(resid3bb)
model3bc<-lmer(SMR_OM~Deer +Honeysuckle + Time + (1|Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model3bc)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: SMR_OM ~ Deer + Honeysuckle + Time + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 724.7 743.8 -354.4 708.7 72
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8325 -0.5136 -0.0615 0.2617 7.1159
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 0.0 0.0
## Residual 412.1 20.3
## Number of obs: 80, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 58.079 5.559 10.447
## DeerExcl -4.060 4.539 -0.894
## HoneysuckleNH 1.320 4.539 0.291
## TimeSpring2018 -17.826 6.419 -2.777
## Timespring2019 -31.274 6.419 -4.872
## Timesummer2018 1.449 6.419 0.226
##
## Correlation of Fixed Effects:
## (Intr) DrExcl HnysNH TS2018 Tm2019
## DeerExcl -0.408
## HoneyscklNH -0.408 0.000
## TmSprng2018 -0.577 0.000 0.000
## Tmsprng2019 -0.577 0.000 0.000 0.500
## Timsmmr2018 -0.577 0.000 0.000 0.500 0.500
resid3bc<-resid(model3bc)
hist(resid3bc)
#split-plot model
library(car)
library(emmeans)
# This sets up your split-plot model
# Deer represents Exclosure or Control/Access
# Honeysuckle represents Removal or Control/Retention
# Site is the site (Bachelor, Reinhart, Western, College, Kramer)
# The random effect is needed to correctly account for degrees of freedom
split.plot.model3c <- lmer(SMR_OM ~ Deer * Honeysuckle +
(1 | Site/Time),
data = Deer_projectdata,
contrasts = list(Deer = contr.sum,
Honeysuckle = contr.sum))
summary(split.plot.model3c)
## Linear mixed model fit by REML ['lmerMod']
## Formula: SMR_OM ~ Deer * Honeysuckle + (1 | Site/Time)
## Data: Deer_projectdata
##
## REML criterion at convergence: 702.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5759 -0.3355 -0.0637 0.1156 6.3473
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time:Site (Intercept) 294.6 17.16
## Site (Intercept) 0.0 0.00
## Residual 328.7 18.13
## Number of obs: 80, groups: Time:Site, 20; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 44.7971 4.3403 10.321
## Deer1 2.0301 2.0269 1.002
## Honeysuckle1 -0.6601 2.0269 -0.326
## Deer1:Honeysuckle1 -1.5976 2.0269 -0.788
##
## Correlation of Fixed Effects:
## (Intr) Deer1 Hnysc1
## Deer1 0.000
## Honeysuckl1 0.000 0.000
## Dr1:Hnysck1 0.000 0.000 0.000
# This gives you the correct type 3 Anova; interpret P-values from here!
anova(split.plot.model3c, type = 3, test.statistic = "F")
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 329.71 329.71 1.0032
## Honeysuckle 1 34.86 34.86 0.1061
## Deer:Honeysuckle 1 204.19 204.19 0.6213
# This allows you to do post-hoc comparisons;
# specifically this line give pairwise comparisons of deer*honeysuckle interaction
emmeans(split.plot.model3c, pairwise ~ Deer * Honeysuckle)
## $emmeans
## Deer Honeysuckle emmean SE df lower.CL upper.CL
## Con H 44.6 5.58 10.6 32.2 56.9
## Excl H 43.7 5.58 10.6 31.4 56.0
## Con NH 49.1 5.58 10.6 36.7 61.4
## Excl NH 41.8 5.58 10.6 29.5 54.2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Con,H - Excl,H 0.865 5.73 57 0.151 0.9988
## Con,H - Con,NH -4.516 5.73 57 -0.788 0.8598
## Con,H - Excl,NH 2.740 5.73 57 0.478 0.9637
## Excl,H - Con,NH -5.380 5.73 57 -0.939 0.7843
## Excl,H - Excl,NH 1.875 5.73 57 0.327 0.9878
## Con,NH - Excl,NH 7.255 5.73 57 1.266 0.5883
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates