Cmodelingproject_Paul
PAUL OJO
March 14, 2020
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, StatQqLinelibrary(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.2library(carData)
#Load the Deer_projectdata
Deer_projectdata<- read.csv(file="Deer_project2.csv")
str(Deer_projectdata)## 'data.frame': 80 obs. of 23 variables:
## $ Time : 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 ...
## $ 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 ...
## $ Nmin : num 1.24 0.96 1.57 1.16 2.79 1.2 0.44 0.39 0.87 1.08 ...
## $ Deer.effect : num 0.39 0.09 -0.29 0.13 1.34 0.67 0.03 -0.5 0.04 0.26 ...
## $ Honeysuckle.effect : num NA NA NA -0.21 1.75 NA NA NA NA NA ...
## $ NetNitr : num 1.32 0.99 1.64 1.22 2 1.17 0.54 0.56 0.99 0.75 ...
## $ Deer.effect.1 : num 0.16 0.01 0.01 0.52 0.64 0.48 0 -0.38 0.25 -0.02 ...
## $ Honeysuckle.effect.1: num NA NA NA NA NA NA NA NA NA NA ...
## $ SMR : num 1.6 1.85 1.47 1.63 2.18 0.93 1.26 0.91 1.05 1.4 ...
## $ Deer.effect.2 : num 0.33 0.07 -0.02 0.03 0.34 0.12 0.1 0.06 -0.14 0.5 ...
## $ Honeysuckle.effect.2: num NA NA NA NA NA NA NA NA NA NA ...
## $ 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 ...head(Deer_projectdata)## Time Site Deer Honeysuckle Nmin Deer.effect
## 1 Spring2018 Kramer Con H 1.24 0.39
## 2 Spring2018 Rhinehart Con H 0.96 0.09
## 3 Spring2018 Bishop Con H 1.57 -0.29
## 4 Spring2018 Western Con H 1.16 0.13
## 5 Spring2018 College Con H 2.79 1.34
## 6 summer2018 Kramer Con H 1.20 0.67
## Honeysuckle.effect NetNitr Deer.effect.1 Honeysuckle.effect.1 SMR
## 1 NA 1.32 0.16 NA 1.60
## 2 NA 0.99 0.01 NA 1.85
## 3 NA 1.64 0.01 NA 1.47
## 4 -0.21 1.22 0.52 NA 1.63
## 5 1.75 2.00 0.64 NA 2.18
## 6 NA 1.17 0.48 NA 0.93
## Deer.effect.2 Honeysuckle.effect.2 Nitr_index sqrNetNitr lgNmin lgSMR
## 1 0.33 NA 1.06 1.15 0.09 0.20
## 2 0.07 NA 1.03 0.99 -0.02 0.27
## 3 -0.02 NA 1.04 1.28 0.20 0.17
## 4 0.03 NA 1.05 1.10 0.06 0.21
## 5 0.34 NA 0.72 1.41 0.45 0.34
## 6 0.12 NA 0.98 1.08 0.08 -0.03
## sqrSMR Nmin_OM NetNitr_OM NetNitr_OM_index SMR_OM lgSMR_OM
## 1 1.26 29.87 32.02 1.07 39.01 1.59
## 2 1.36 22.86 23.60 1.03 43.45 1.64
## 3 1.21 36.63 38.36 1.05 34.16 1.53
## 4 1.28 27.69 29.22 1.06 39.00 1.59
## 5 1.48 66.10 47.39 0.72 51.66 1.71
## 6 0.96 59.11 57.80 0.98 45.89 1.66#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.3902bartlett.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.4459bartlett.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-11plot(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 ' ' 1model1aov<-aov(NetNitr_OM~Deer*Honeysuckle*Time, data=Deer_projectdata)
posthoc1 <- TukeyHSD(x=model1aov, 'Time', conf.level=0.95)
posthoc1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = NetNitr_OM ~ Deer * Honeysuckle * Time, data = Deer_projectdata)
##
## $Time
## diff lwr upr p adj
## Spring2018-fall2018 -1.55100 -8.935552 5.833552 0.9450517
## spring2019-fall2018 -20.26919 -27.750272 -12.788099 0.0000000
## summer2018-fall2018 8.65400 1.269448 16.038552 0.0152687
## spring2019-Spring2018 -18.71819 -26.199272 -11.237099 0.0000001
## summer2018-Spring2018 10.20500 2.820448 17.589552 0.0029623
## summer2018-spring2019 28.92319 21.442099 36.404272 0.0000000plot(posthoc1)
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.3model1b<-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.013anova(model1b)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 167.85 167.85 0.8572
## Honeysuckle 1 251.74 251.74 1.2856resid1b<-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.016anova(model1ba)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 121.44 121.44 1.3660
## Honeysuckle 1 318.21 318.21 3.5793resid1ba<-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.013anova(model1bb)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 123.39 123.39 1.3494
## Honeysuckle 1 314.69 314.69 3.4415resid1bb<-resid(model1bb)
hist(resid1bb)
model1bc<-lmer(NetNitr_OM~Deer +Honeysuckle + Time + (1|Site),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 | Site)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 591.5 610.5 -287.8 575.5 71
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8950 -0.5074 -0.1525 0.4032 4.2367
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 4.958 2.227
## Residual 81.819 9.045
## Number of obs: 79, groups: Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 32.926 2.670 12.333
## DeerExcl -2.425 2.037 -1.191
## HoneysuckleNH -4.016 2.037 -1.972
## TimeSpring2018 -1.551 2.860 -0.542
## Timespring2019 -20.300 2.900 -7.000
## Timesummer2018 8.654 2.860 3.025
##
## Correlation of Fixed Effects:
## (Intr) DrExcl HnysNH TS2018 Tm2019
## DeerExcl -0.376
## HoneyscklNH -0.376 -0.014
## TmSprng2018 -0.536 0.000 0.000
## Tmsprng2019 -0.528 -0.019 0.019 0.493
## Timsmmr2018 -0.536 0.000 0.000 0.500 0.493anova(model1bc)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 168.9 168.85 2.0637
## Honeysuckle 1 250.5 250.52 3.0618
## Time 3 8608.6 2869.52 35.0716anova(model1b,model1ba,model1bb,model1bc)## Data: Deer_projectdata
## Models:
## model1b: NetNitr_OM ~ Deer + Honeysuckle + (1 | Site)
## model1bb: NetNitr_OM ~ Deer + Honeysuckle + (1 | Time)
## model1ba: NetNitr_OM ~ Deer + Honeysuckle + (1 | Site/Time)
## model1bc: NetNitr_OM ~ Deer + Honeysuckle + Time + (1 | Site)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model1b 5 651.09 662.94 -320.54 641.09
## model1bb 5 603.65 615.49 -296.82 593.65 47.443 0 < 2.2e-16 ***
## model1ba 6 626.11 640.32 -307.05 614.11 0.000 1 1
## model1bc 8 591.51 610.46 -287.75 575.51 38.600 2 4.151e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resid1bc<-resid(model1bc)
hist(resid1bc)
#split-plot model
library(car)
library(emmeans)## Warning: package 'emmeans' was built under R version 3.5.3# 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 + Deer * Honeysuckle +(1 | Time/Site),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 + Deer * Honeysuckle + (1 | Time/Site)
## Data: Deer_projectdata
##
## REML criterion at convergence: 579.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4608 -0.4972 -0.0766 0.3693 3.8333
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site:Time (Intercept) 2.723 1.650
## Time (Intercept) 142.446 11.935
## Residual 89.808 9.477
## Number of obs: 79, groups: Site:Time, 20; Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.433 6.073 4.352
## Deer1 1.239 1.067 1.161
## Honeysuckle1 1.982 1.067 1.857
## Deer1:Honeysuckle1 -1.669 1.067 -1.565
##
## Correlation of Fixed Effects:
## (Intr) Deer1 Hnysc1
## Deer1 0.002
## Honeysuckl1 -0.002 -0.014
## Dr1:Hnysck1 -0.002 -0.014 0.014# 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.94 121.94 1.3578
## Honeysuckle 1 317.09 317.09 3.5308
## Deer:Honeysuckle 1 219.89 219.89 2.4485# 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 6.34 3.57 9.50 46.5
## Excl H 28.8 6.34 3.57 10.36 47.3
## Con NH 27.4 6.36 3.61 8.92 45.8
## Excl NH 21.5 6.34 3.57 3.06 40.0
##
## 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.9916
## Con,H - Con,NH 0.624 3.04 56.7 0.205 0.9969
## Con,H - Excl,NH 6.441 3.00 56.0 2.149 0.1504
## Excl,H - Con,NH 1.486 3.04 56.7 0.489 0.9613
## Excl,H - Excl,NH 7.302 3.00 56.0 2.437 0.0818
## Con,NH - Excl,NH 5.816 3.04 56.7 1.913 0.2343
##
## 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-08anova(model2)## Analysis of Variance Table
##
## Response: Nmin_OM
## Df Sum Sq Mean Sq F value Pr(>F)
## Deer 1 77.7 77.74 0.7957 0.37572
## Honeysuckle 1 301.9 301.86 3.0899 0.08356 .
## Time 3 8584.7 2861.56 29.2907 4.858e-12 ***
## Deer:Honeysuckle 1 3.0 2.96 0.0303 0.86226
## Deer:Time 3 528.7 176.24 1.8040 0.15536
## Honeysuckle:Time 3 32.9 10.98 0.1124 0.95256
## Deer:Honeysuckle:Time 3 261.1 87.05 0.8910 0.45070
## Residuals 64 6252.5 97.70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resid2<-resid(model2)
hist(resid2)
model2aov<-aov(Nmin_OM~Deer*Honeysuckle*Time, data=Deer_projectdata)
posthoc2 <- TukeyHSD(x=model2aov, 'Time', conf.level=0.95)
posthoc2## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Nmin_OM ~ Deer * Honeysuckle * Time, data = Deer_projectdata)
##
## $Time
## diff lwr upr p adj
## Spring2018-fall2018 6.1845 -2.060399 14.429399 0.2067939
## spring2019-fall2018 -15.2500 -23.494899 -7.005101 0.0000432
## summer2018-fall2018 12.7185 4.473601 20.963399 0.0007458
## spring2019-Spring2018 -21.4345 -29.679399 -13.189601 0.0000000
## summer2018-Spring2018 6.5340 -1.710899 14.778899 0.1672289
## summer2018-spring2019 27.9685 19.723601 36.213399 0.0000000plot(posthoc2)
#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.000anova(model2b)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 77.736 77.736 0.3991
## Honeysuckle 1 301.865 301.865 1.5497resid2b<-resid(model2b)
plot(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.000anova(model2ba)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 77.736 77.736 0.9992
## Honeysuckle 1 301.865 301.865 3.8801resid2ba<-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.000anova(model2bb)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 77.736 77.736 0.8347
## Honeysuckle 1 301.865 301.865 3.2411resid2bb<-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.500anova(model2bc)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 77.7 77.74 0.8786
## Honeysuckle 1 301.9 301.86 3.4117
## Time 3 8584.7 2861.56 32.3420anova(model2b,model2ba,model2bb,model2bc)## Data: Deer_projectdata
## Models:
## model2b: Nmin_OM ~ Deer + Honeysuckle + (1 | Site)
## model2bb: Nmin_OM ~ Deer + Honeysuckle + (1 | Time)
## model2ba: Nmin_OM ~ Deer + Honeysuckle + (1 | Site/Time)
## model2bc: Nmin_OM ~ Deer + Honeysuckle + Time + (1 | Time)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2b 5 659.18 671.09 -324.59 649.18
## model2bb 5 612.30 624.21 -301.15 602.30 46.874 0 < 2.2e-16 ***
## model2ba 6 626.47 640.76 -307.23 614.47 0.000 1 1
## model2bc 8 601.65 620.71 -292.83 585.65 28.816 2 5.529e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resid2bc<-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-07plot(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 ' ' 1resid3<-resid(model3)
hist(resid3)
model3aov <- aov(lgSMR_OM~Deer*Honeysuckle*Time, data=Deer_projectdata)
posthoc3 <- TukeyHSD(x=model3aov, 'Time', conf.level=0.95)
posthoc3## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = lgSMR_OM ~ Deer * Honeysuckle * Time, data = Deer_projectdata)
##
## $Time
## diff lwr upr p adj
## Spring2018-fall2018 -0.1445 -0.26865718 -0.02034282 0.0161751
## spring2019-fall2018 -0.3520 -0.47615718 -0.22784282 0.0000000
## summer2018-fall2018 -0.0060 -0.13015718 0.11815718 0.9992521
## spring2019-Spring2018 -0.2075 -0.33165718 -0.08334282 0.0002332
## summer2018-Spring2018 0.1385 0.01434282 0.26265718 0.0228944
## summer2018-spring2019 0.3460 0.22184282 0.47015718 0.0000000plot(posthoc3)
#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.000resid3b<-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.000resid3ba<-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.000resid3bb<-resid(model3bb)
hist(resid3bb)
anova(model3bb)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 329.71 329.71 0.7601
## Honeysuckle 1 34.86 34.86 0.0804model3bc<-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.500anova(model3bc)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Deer 1 329.7 329.7 0.8001
## Honeysuckle 1 34.9 34.9 0.0846
## Time 3 14605.3 4868.4 11.8144anova(model3b,model3ba,model3bb,model3bc)## Data: Deer_projectdata
## Models:
## model3b: SMR_OM ~ Deer + Honeysuckle + (1 | Site)
## model3bb: SMR_OM ~ Deer + Honeysuckle + (1 | Time)
## model3ba: SMR_OM ~ Deer + Honeysuckle + (1 | Site/Time)
## model3bc: SMR_OM ~ Deer + Honeysuckle + Time + (1 | Time)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model3b 5 747.34 759.25 -368.67 737.34
## model3bb 5 731.35 743.26 -360.68 721.35 15.9877 0 < 2e-16 ***
## model3ba 6 729.64 743.93 -358.82 717.64 3.7153 1 0.05391 .
## model3bc 8 724.73 743.78 -354.36 708.73 8.9095 2 0.01162 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resid3bc<-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```