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(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_project2.csv")
str(Deer_projectdata)
## 'data.frame': 80 obs. of 24 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 ...
## $ Treatment : Factor w/ 4 levels "ConH","ConNH",..: 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 Treatment Nmin Deer.effect
## 1 Spring2018 Kramer Con H ConH 1.24 0.39
## 2 Spring2018 Rhinehart Con H ConH 0.96 0.09
## 3 Spring2018 Bishop Con H ConH 1.57 -0.29
## 4 Spring2018 Western Con H ConH 1.16 0.13
## 5 Spring2018 College Con H ConH 2.79 1.34
## 6 summer2018 Kramer Con H ConH 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.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~Treatment*Time, data=Deer_projectdata,na.action=na.omit)
summary(model1)
##
## Call:
## lm(formula = NetNitr_OM ~ Treatment * 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 ***
## TreatmentConNH 1.292 5.597 0.231 0.818175
## TreatmentExclH 20.472 5.597 3.658 0.000521 ***
## TreatmentExclNH -0.540 5.597 -0.096 0.923440
## 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 **
## TreatmentConNH:TimeSpring2018 -7.024 7.915 -0.887 0.378209
## TreatmentExclH:TimeSpring2018 -28.176 7.915 -3.560 0.000712 ***
## TreatmentExclNH:TimeSpring2018 -9.876 7.915 -1.248 0.216724
## TreatmentConNH:Timespring2019 -3.587 8.158 -0.440 0.661679
## TreatmentExclH:Timespring2019 -27.122 7.915 -3.427 0.001081 **
## TreatmentExclNH:Timespring2019 -8.396 7.915 -1.061 0.292829
## TreatmentConNH:Timesummer2018 3.178 7.915 0.402 0.689390
## TreatmentExclH:Timesummer2018 -23.144 7.915 -2.924 0.004794 **
## TreatmentExclNH:Timesummer2018 -5.330 7.915 -0.673 0.503140
## ---
## 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)
## Treatment 3 701.7 233.90 2.9870 0.03768 *
## Time 3 8549.0 2849.67 36.3923 9.07e-14 ***
## Treatment:Time 9 1705.2 189.46 2.4196 0.01984 *
## Residuals 63 4933.2 78.30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model1aov<-aov(NetNitr_OM~Treatment*Time, data=Deer_projectdata)
posthoc1 <- TukeyHSD(x=model1aov,conf.level=0.95)
posthoc1
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = NetNitr_OM ~ Treatment * Time, data = Deer_projectdata)
##
## $Treatment
## diff lwr upr p adj
## ConNH-ConH 0.2586316 -7.222455 7.73971806 0.9997249
## ExclH-ConH 0.8615000 -6.523052 8.24605228 0.9897809
## ExclNH-ConH -6.4405000 -13.825052 0.94405228 0.1085466
## ExclH-ConNH 0.6028684 -6.878218 8.08395490 0.9965661
## ExclNH-ConNH -6.6991316 -14.180218 0.78195490 0.0949339
## ExclNH-ExclH -7.3020000 -14.686552 0.08255228 0.0537128
##
## $Time
## diff lwr upr p adj
## Spring2018-fall2018 -1.55100 -8.935552 5.833552 0.9450517
## spring2019-fall2018 -20.16712 -27.648206 -12.686033 0.0000000
## summer2018-fall2018 8.65400 1.269448 16.038552 0.0152687
## spring2019-Spring2018 -18.61612 -26.097206 -11.135033 0.0000001
## summer2018-Spring2018 10.20500 2.820448 17.589552 0.0029623
## summer2018-spring2019 28.82112 21.340033 36.302206 0.0000000
##
## $`Treatment:Time`
## diff lwr upr
## ConNH:fall2018-ConH:fall2018 1.292 -18.67526739 21.2592674
## ExclH:fall2018-ConH:fall2018 20.472 0.50473261 40.4392674
## ExclNH:fall2018-ConH:fall2018 -0.540 -20.50726739 19.4272674
## ConH:Spring2018-ConH:fall2018 9.718 -10.24926739 29.6852674
## ConNH:Spring2018-ConH:fall2018 3.986 -15.98126739 23.9532674
## ExclH:Spring2018-ConH:fall2018 2.014 -17.95326739 21.9812674
## ExclNH:Spring2018-ConH:fall2018 -0.698 -20.66526739 19.2692674
## ConH:spring2019-ConH:fall2018 -10.360 -30.32726739 9.6072674
## ConNH:spring2019-ConH:fall2018 -12.655 -33.83348526 8.5234853
## ExclH:spring2019-ConH:fall2018 -17.010 -36.97726739 2.9572674
## ExclNH:spring2019-ConH:fall2018 -19.296 -39.26326739 0.6712674
## ConH:summer2018-ConH:fall2018 14.978 -4.98926739 34.9452674
## ConNH:summer2018-ConH:fall2018 19.448 -0.51926739 39.4152674
## ExclH:summer2018-ConH:fall2018 12.306 -7.66126739 32.2732674
## ExclNH:summer2018-ConH:fall2018 9.108 -10.85926739 29.0752674
## ExclH:fall2018-ConNH:fall2018 19.180 -0.78726739 39.1472674
## ExclNH:fall2018-ConNH:fall2018 -1.832 -21.79926739 18.1352674
## ConH:Spring2018-ConNH:fall2018 8.426 -11.54126739 28.3932674
## ConNH:Spring2018-ConNH:fall2018 2.694 -17.27326739 22.6612674
## ExclH:Spring2018-ConNH:fall2018 0.722 -19.24526739 20.6892674
## ExclNH:Spring2018-ConNH:fall2018 -1.990 -21.95726739 17.9772674
## ConH:spring2019-ConNH:fall2018 -11.652 -31.61926739 8.3152674
## ConNH:spring2019-ConNH:fall2018 -13.947 -35.12548526 7.2314853
## ExclH:spring2019-ConNH:fall2018 -18.302 -38.26926739 1.6652674
## ExclNH:spring2019-ConNH:fall2018 -20.588 -40.55526739 -0.6207326
## ConH:summer2018-ConNH:fall2018 13.686 -6.28126739 33.6532674
## ConNH:summer2018-ConNH:fall2018 18.156 -1.81126739 38.1232674
## ExclH:summer2018-ConNH:fall2018 11.014 -8.95326739 30.9812674
## ExclNH:summer2018-ConNH:fall2018 7.816 -12.15126739 27.7832674
## ExclNH:fall2018-ExclH:fall2018 -21.012 -40.97926739 -1.0447326
## ConH:Spring2018-ExclH:fall2018 -10.754 -30.72126739 9.2132674
## ConNH:Spring2018-ExclH:fall2018 -16.486 -36.45326739 3.4812674
## ExclH:Spring2018-ExclH:fall2018 -18.458 -38.42526739 1.5092674
## ExclNH:Spring2018-ExclH:fall2018 -21.170 -41.13726739 -1.2027326
## ConH:spring2019-ExclH:fall2018 -30.832 -50.79926739 -10.8647326
## ConNH:spring2019-ExclH:fall2018 -33.127 -54.30548526 -11.9485147
## ExclH:spring2019-ExclH:fall2018 -37.482 -57.44926739 -17.5147326
## ExclNH:spring2019-ExclH:fall2018 -39.768 -59.73526739 -19.8007326
## ConH:summer2018-ExclH:fall2018 -5.494 -25.46126739 14.4732674
## ConNH:summer2018-ExclH:fall2018 -1.024 -20.99126739 18.9432674
## ExclH:summer2018-ExclH:fall2018 -8.166 -28.13326739 11.8012674
## ExclNH:summer2018-ExclH:fall2018 -11.364 -31.33126739 8.6032674
## ConH:Spring2018-ExclNH:fall2018 10.258 -9.70926739 30.2252674
## ConNH:Spring2018-ExclNH:fall2018 4.526 -15.44126739 24.4932674
## ExclH:Spring2018-ExclNH:fall2018 2.554 -17.41326739 22.5212674
## ExclNH:Spring2018-ExclNH:fall2018 -0.158 -20.12526739 19.8092674
## ConH:spring2019-ExclNH:fall2018 -9.820 -29.78726739 10.1472674
## ConNH:spring2019-ExclNH:fall2018 -12.115 -33.29348526 9.0634853
## ExclH:spring2019-ExclNH:fall2018 -16.470 -36.43726739 3.4972674
## ExclNH:spring2019-ExclNH:fall2018 -18.756 -38.72326739 1.2112674
## ConH:summer2018-ExclNH:fall2018 15.518 -4.44926739 35.4852674
## ConNH:summer2018-ExclNH:fall2018 19.988 0.02073261 39.9552674
## ExclH:summer2018-ExclNH:fall2018 12.846 -7.12126739 32.8132674
## ExclNH:summer2018-ExclNH:fall2018 9.648 -10.31926739 29.6152674
## ConNH:Spring2018-ConH:Spring2018 -5.732 -25.69926739 14.2352674
## ExclH:Spring2018-ConH:Spring2018 -7.704 -27.67126739 12.2632674
## ExclNH:Spring2018-ConH:Spring2018 -10.416 -30.38326739 9.5512674
## ConH:spring2019-ConH:Spring2018 -20.078 -40.04526739 -0.1107326
## ConNH:spring2019-ConH:Spring2018 -22.373 -43.55148526 -1.1945147
## ExclH:spring2019-ConH:Spring2018 -26.728 -46.69526739 -6.7607326
## ExclNH:spring2019-ConH:Spring2018 -29.014 -48.98126739 -9.0467326
## ConH:summer2018-ConH:Spring2018 5.260 -14.70726739 25.2272674
## ConNH:summer2018-ConH:Spring2018 9.730 -10.23726739 29.6972674
## ExclH:summer2018-ConH:Spring2018 2.588 -17.37926739 22.5552674
## ExclNH:summer2018-ConH:Spring2018 -0.610 -20.57726739 19.3572674
## ExclH:Spring2018-ConNH:Spring2018 -1.972 -21.93926739 17.9952674
## ExclNH:Spring2018-ConNH:Spring2018 -4.684 -24.65126739 15.2832674
## ConH:spring2019-ConNH:Spring2018 -14.346 -34.31326739 5.6212674
## ConNH:spring2019-ConNH:Spring2018 -16.641 -37.81948526 4.5374853
## ExclH:spring2019-ConNH:Spring2018 -20.996 -40.96326739 -1.0287326
## ExclNH:spring2019-ConNH:Spring2018 -23.282 -43.24926739 -3.3147326
## ConH:summer2018-ConNH:Spring2018 10.992 -8.97526739 30.9592674
## ConNH:summer2018-ConNH:Spring2018 15.462 -4.50526739 35.4292674
## ExclH:summer2018-ConNH:Spring2018 8.320 -11.64726739 28.2872674
## ExclNH:summer2018-ConNH:Spring2018 5.122 -14.84526739 25.0892674
## ExclNH:Spring2018-ExclH:Spring2018 -2.712 -22.67926739 17.2552674
## ConH:spring2019-ExclH:Spring2018 -12.374 -32.34126739 7.5932674
## ConNH:spring2019-ExclH:Spring2018 -14.669 -35.84748526 6.5094853
## ExclH:spring2019-ExclH:Spring2018 -19.024 -38.99126739 0.9432674
## ExclNH:spring2019-ExclH:Spring2018 -21.310 -41.27726739 -1.3427326
## ConH:summer2018-ExclH:Spring2018 12.964 -7.00326739 32.9312674
## ConNH:summer2018-ExclH:Spring2018 17.434 -2.53326739 37.4012674
## ExclH:summer2018-ExclH:Spring2018 10.292 -9.67526739 30.2592674
## ExclNH:summer2018-ExclH:Spring2018 7.094 -12.87326739 27.0612674
## ConH:spring2019-ExclNH:Spring2018 -9.662 -29.62926739 10.3052674
## ConNH:spring2019-ExclNH:Spring2018 -11.957 -33.13548526 9.2214853
## ExclH:spring2019-ExclNH:Spring2018 -16.312 -36.27926739 3.6552674
## ExclNH:spring2019-ExclNH:Spring2018 -18.598 -38.56526739 1.3692674
## ConH:summer2018-ExclNH:Spring2018 15.676 -4.29126739 35.6432674
## ConNH:summer2018-ExclNH:Spring2018 20.146 0.17873261 40.1132674
## ExclH:summer2018-ExclNH:Spring2018 13.004 -6.96326739 32.9712674
## ExclNH:summer2018-ExclNH:Spring2018 9.806 -10.16126739 29.7732674
## ConNH:spring2019-ConH:spring2019 -2.295 -23.47348526 18.8834853
## ExclH:spring2019-ConH:spring2019 -6.650 -26.61726739 13.3172674
## ExclNH:spring2019-ConH:spring2019 -8.936 -28.90326739 11.0312674
## ConH:summer2018-ConH:spring2019 25.338 5.37073261 45.3052674
## ConNH:summer2018-ConH:spring2019 29.808 9.84073261 49.7752674
## ExclH:summer2018-ConH:spring2019 22.666 2.69873261 42.6332674
## ExclNH:summer2018-ConH:spring2019 19.468 -0.49926739 39.4352674
## ExclH:spring2019-ConNH:spring2019 -4.355 -25.53348526 16.8234853
## ExclNH:spring2019-ConNH:spring2019 -6.641 -27.81948526 14.5374853
## ConH:summer2018-ConNH:spring2019 27.633 6.45451474 48.8114853
## ConNH:summer2018-ConNH:spring2019 32.103 10.92451474 53.2814853
## ExclH:summer2018-ConNH:spring2019 24.961 3.78251474 46.1394853
## ExclNH:summer2018-ConNH:spring2019 21.763 0.58451474 42.9414853
## ExclNH:spring2019-ExclH:spring2019 -2.286 -22.25326739 17.6812674
## ConH:summer2018-ExclH:spring2019 31.988 12.02073261 51.9552674
## ConNH:summer2018-ExclH:spring2019 36.458 16.49073261 56.4252674
## ExclH:summer2018-ExclH:spring2019 29.316 9.34873261 49.2832674
## ExclNH:summer2018-ExclH:spring2019 26.118 6.15073261 46.0852674
## ConH:summer2018-ExclNH:spring2019 34.274 14.30673261 54.2412674
## ConNH:summer2018-ExclNH:spring2019 38.744 18.77673261 58.7112674
## ExclH:summer2018-ExclNH:spring2019 31.602 11.63473261 51.5692674
## ExclNH:summer2018-ExclNH:spring2019 28.404 8.43673261 48.3712674
## ConNH:summer2018-ConH:summer2018 4.470 -15.49726739 24.4372674
## ExclH:summer2018-ConH:summer2018 -2.672 -22.63926739 17.2952674
## ExclNH:summer2018-ConH:summer2018 -5.870 -25.83726739 14.0972674
## ExclH:summer2018-ConNH:summer2018 -7.142 -27.10926739 12.8252674
## ExclNH:summer2018-ConNH:summer2018 -10.340 -30.30726739 9.6272674
## ExclNH:summer2018-ExclH:summer2018 -3.198 -23.16526739 16.7692674
## p adj
## ConNH:fall2018-ConH:fall2018 1.0000000
## ExclH:fall2018-ConH:fall2018 0.0389927
## ExclNH:fall2018-ConH:fall2018 1.0000000
## ConH:Spring2018-ConH:fall2018 0.9285022
## ConNH:Spring2018-ConH:fall2018 0.9999960
## ExclH:Spring2018-ConH:fall2018 1.0000000
## ExclNH:Spring2018-ConH:fall2018 1.0000000
## ConH:spring2019-ConH:fall2018 0.8872718
## ConNH:spring2019-ConH:fall2018 0.7380404
## ExclH:spring2019-ConH:fall2018 0.1828202
## ExclNH:spring2019-ConH:fall2018 0.0688067
## ConH:summer2018-ConH:fall2018 0.3684197
## ConNH:summer2018-ConH:fall2018 0.0640835
## ExclH:summer2018-ConH:fall2018 0.6943878
## ExclNH:summer2018-ConH:fall2018 0.9568794
## ExclH:fall2018-ConNH:fall2018 0.0726094
## ExclNH:fall2018-ConNH:fall2018 1.0000000
## ConH:Spring2018-ConNH:fall2018 0.9777761
## ConNH:Spring2018-ConNH:fall2018 1.0000000
## ExclH:Spring2018-ConNH:fall2018 1.0000000
## ExclNH:Spring2018-ConNH:fall2018 1.0000000
## ConH:spring2019-ConNH:fall2018 0.7689380
## ConNH:spring2019-ConNH:fall2018 0.5902619
## ExclH:spring2019-ConNH:fall2018 0.1075871
## ExclNH:spring2019-ConNH:fall2018 0.0367902
## ConH:summer2018-ConNH:fall2018 0.5226596
## ConNH:summer2018-ConNH:fall2018 0.1145718
## ExclH:summer2018-ConNH:fall2018 0.8329376
## ExclNH:summer2018-ConNH:fall2018 0.9888961
## ExclNH:fall2018-ExclH:fall2018 0.0296560
## ConH:Spring2018-ExclH:fall2018 0.8559940
## ConNH:Spring2018-ExclH:fall2018 0.2226552
## ExclH:Spring2018-ExclH:fall2018 0.1005133
## ExclNH:Spring2018-ExclH:fall2018 0.0273342
## ConH:spring2019-ExclH:fall2018 0.0000773
## ConNH:spring2019-ExclH:fall2018 0.0000590
## ExclH:spring2019-ExclH:fall2018 0.0000008
## ExclNH:spring2019-ExclH:fall2018 0.0000002
## ConH:summer2018-ExclH:fall2018 0.9997607
## ConNH:summer2018-ExclH:fall2018 1.0000000
## ExclH:summer2018-ExclH:fall2018 0.9832572
## ExclNH:summer2018-ExclH:fall2018 0.7990551
## ConH:Spring2018-ExclNH:fall2018 0.8946306
## ConNH:Spring2018-ExclNH:fall2018 0.9999787
## ExclH:Spring2018-ExclNH:fall2018 1.0000000
## ExclNH:Spring2018-ExclNH:fall2018 1.0000000
## ConH:spring2019-ExclNH:fall2018 0.9227538
## ConNH:spring2019-ExclNH:fall2018 0.7931307
## ExclH:spring2019-ExclNH:fall2018 0.2239608
## ExclNH:spring2019-ExclNH:fall2018 0.0880685
## ConH:summer2018-ExclNH:fall2018 0.3111279
## ConNH:summer2018-ExclNH:fall2018 0.0494989
## ExclH:summer2018-ExclNH:fall2018 0.6283210
## ExclNH:summer2018-ExclNH:fall2018 0.9322756
## ConNH:Spring2018-ConH:Spring2018 0.9996038
## ExclH:Spring2018-ConH:Spring2018 0.9903368
## ExclNH:Spring2018-ConH:Spring2018 0.8831017
## ConH:spring2019-ConH:Spring2018 0.0473746
## ConNH:spring2019-ConH:Spring2018 0.0284375
## ExclH:spring2019-ConH:Spring2018 0.0011097
## ExclNH:spring2019-ConH:Spring2018 0.0002574
## ConH:summer2018-ConH:Spring2018 0.9998588
## ConNH:summer2018-ConH:Spring2018 0.9278414
## ExclH:summer2018-ConH:Spring2018 1.0000000
## ExclNH:summer2018-ConH:Spring2018 1.0000000
## ExclH:Spring2018-ConNH:Spring2018 1.0000000
## ExclNH:Spring2018-ConNH:Spring2018 0.9999669
## ConH:spring2019-ConNH:Spring2018 0.4414874
## ConNH:spring2019-ConNH:Spring2018 0.2940943
## ExclH:spring2019-ConNH:Spring2018 0.0299008
## ExclNH:spring2019-ConNH:Spring2018 0.0086845
## ConH:summer2018-ConNH:Spring2018 0.8349604
## ConNH:summer2018-ConNH:Spring2018 0.3168202
## ExclH:summer2018-ConNH:Spring2018 0.9801567
## ExclNH:summer2018-ConNH:Spring2018 0.9998982
## ExclNH:Spring2018-ExclH:Spring2018 1.0000000
## ConH:spring2019-ExclH:Spring2018 0.6862429
## ConNH:spring2019-ExclH:Spring2018 0.5046572
## ExclH:spring2019-ExclH:Spring2018 0.0780052
## ExclNH:spring2019-ExclH:Spring2018 0.0254158
## ConH:summer2018-ExclH:Spring2018 0.6135650
## ConNH:summer2018-ExclH:Spring2018 0.1546598
## ExclH:summer2018-ExclH:Spring2018 0.8922117
## ExclNH:summer2018-ExclH:Spring2018 0.9957751
## ConH:spring2019-ExclNH:Spring2018 0.9315320
## ConNH:spring2019-ExclNH:Spring2018 0.8081477
## ExclH:spring2019-ExclNH:Spring2018 0.2371400
## ExclNH:spring2019-ExclNH:Spring2018 0.0944962
## ConH:summer2018-ExclNH:Spring2018 0.2953949
## ConNH:summer2018-ExclNH:Spring2018 0.0458232
## ExclH:summer2018-ExclNH:Spring2018 0.6085470
## ExclNH:summer2018-ExclNH:Spring2018 0.9235605
## ConNH:spring2019-ConH:spring2019 1.0000000
## ExclH:spring2019-ConH:spring2019 0.9978703
## ExclNH:spring2019-ConH:spring2019 0.9631402
## ConH:summer2018-ConH:spring2019 0.0026062
## ConNH:summer2018-ConH:spring2019 0.0001528
## ExclH:summer2018-ConH:spring2019 0.0122591
## ExclNH:summer2018-ConH:spring2019 0.0634833
## ExclH:spring2019-ConNH:spring2019 0.9999940
## ExclNH:spring2019-ConNH:spring2019 0.9989048
## ConH:summer2018-ConNH:spring2019 0.0016863
## ConNH:summer2018-ConNH:spring2019 0.0001127
## ExclH:summer2018-ConNH:spring2019 0.0075276
## ExclNH:summer2018-ConNH:spring2019 0.0380989
## ExclNH:spring2019-ExclH:spring2019 1.0000000
## ConH:summer2018-ExclH:spring2019 0.0000354
## ConNH:summer2018-ExclH:spring2019 0.0000016
## ExclH:summer2018-ExclH:spring2019 0.0002112
## ExclNH:summer2018-ExclH:spring2019 0.0016199
## ConH:summer2018-ExclNH:spring2019 0.0000074
## ConNH:summer2018-ExclNH:spring2019 0.0000003
## ExclH:summer2018-ExclNH:spring2019 0.0000460
## ExclNH:summer2018-ExclNH:spring2019 0.0003825
## ConNH:summer2018-ConH:summer2018 0.9999819
## ExclH:summer2018-ConH:summer2018 1.0000000
## ExclNH:summer2018-ConH:summer2018 0.9994765
## ExclH:summer2018-ConNH:summer2018 0.9954703
## ExclNH:summer2018-ConNH:summer2018 0.8887389
## ExclNH:summer2018-ExclH:summer2018 0.9999998
plot(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.3
model1b<-lmer(NetNitr_OM~Treatment + (1|Site),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model1b)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: NetNitr_OM ~ Treatment + (1 | Site)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 651.6 665.9 -319.8 639.6 73
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.95420 -0.63641 -0.04861 0.50641 2.87370
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.0 0.00
## Residual 192.2 13.87
## Number of obs: 79, groups: Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 27.9840 3.1004 9.026
## TreatmentConNH 0.2586 4.4419 0.058
## TreatmentExclH 0.8615 4.3846 0.196
## TreatmentExclNH -6.4405 4.3846 -1.469
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH
## TretmntCnNH -0.698
## TrtmntExclH -0.707 0.494
## TrtmntExcNH -0.707 0.494 0.500
anova(model1b)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 701.69 233.9 1.2167
resid1b<-resid(model1b)
hist(resid1b)
model1ba<-lmer(NetNitr_OM~Treatment + (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 ~ Treatment + (1 | Site/Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 625.6 642.2 -305.8 611.6 72
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5967 -0.5954 -0.1883 0.4220 3.7259
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time:Site (Intercept) 109.58 10.468
## Site (Intercept) 0.00 0.000
## Residual 85.44 9.243
## Number of obs: 79, groups: Time:Site, 20; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 27.9840 3.1226 8.962
## TreatmentConNH -0.6386 2.9712 -0.215
## TreatmentExclH 0.8615 2.9230 0.295
## TreatmentExclNH -6.4405 2.9230 -2.203
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH
## TretmntCnNH -0.460
## TrtmntExclH -0.468 0.492
## TrtmntExcNH -0.468 0.492 0.500
anova(model1ba)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 658.25 219.41 2.5681
resid1ba<-resid(model1ba)
hist(resid1ba)
model1bb<-lmer(NetNitr_OM~Treatment + (1|Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model1bb)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: NetNitr_OM ~ Treatment + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 603.2 617.4 -295.6 591.2 73
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4882 -0.5273 -0.1337 0.4128 3.9019
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 105.95 10.293
## Residual 88.52 9.408
## Number of obs: 79, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 27.9840 5.5600 5.033
## TreatmentConNH -0.5952 3.0156 -0.197
## TreatmentExclH 0.8615 2.9752 0.290
## TreatmentExclNH -6.4405 2.9752 -2.165
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH
## TretmntCnNH -0.264
## TrtmntExclH -0.268 0.493
## TrtmntExcNH -0.268 0.493 0.500
anova(model1bb)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 659.98 219.99 2.4854
resid1bb<-resid(model1bb)
hist(resid1bb)
model1bc<-lmer(NetNitr_OM~Treatment + 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 ~ Treatment + Time + (1 | Site)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 590.8 612.1 -286.4 572.8 70
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7471 -0.5597 -0.0707 0.2644 4.1281
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 5.13 2.265
## Residual 78.87 8.881
## Number of obs: 79, groups: Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 31.2601 2.8179 11.093
## TreatmentConNH -0.6373 2.8476 -0.224
## TreatmentExclH 0.8615 2.8083 0.307
## TreatmentExclNH -6.4405 2.8083 -2.293
## TimeSpring2018 -1.5510 2.8083 -0.552
## Timespring2019 -20.2073 2.8476 -7.096
## Timesummer2018 8.6540 2.8083 3.082
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH TrtENH TS2018 Tm2019
## TretmntCnNH -0.498
## TrtmntExclH -0.498 0.493
## TrtmntExcNH -0.498 0.493 0.500
## TmSprng2018 -0.498 0.000 0.000 0.000
## Tmsprng2019 -0.498 0.027 0.000 0.000 0.493
## Timsmmr2018 -0.498 0.000 0.000 0.000 0.500 0.493
anova(model1bc)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 702.8 234.28 2.9705
## Time 3 8544.0 2847.99 36.1111
anova(model1b,model1ba,model1bb,model1bc)
## Data: Deer_projectdata
## Models:
## model1b: NetNitr_OM ~ Treatment + (1 | Site)
## model1bb: NetNitr_OM ~ Treatment + (1 | Time)
## model1ba: NetNitr_OM ~ Treatment + (1 | Site/Time)
## model1bc: NetNitr_OM ~ Treatment + Time + (1 | Site)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model1b 6 651.64 665.85 -319.82 639.64
## model1bb 6 603.18 617.40 -295.59 591.18 48.455 0 < 2.2e-16 ***
## model1ba 7 625.60 642.18 -305.80 611.60 0.000 1 1
## model1bc 9 590.78 612.11 -286.39 572.78 38.818 2 3.722e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
resid1bc<-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~Treatment* Time, data=Deer_projectdata,na.action=na.omit)
summary(model2)
##
## Call:
## lm(formula = Nmin_OM ~ Treatment * 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 ***
## TreatmentConNH 2.686 6.251 0.430 0.668876
## TreatmentExclH 13.302 6.251 2.128 0.037203 *
## TreatmentExclNH 3.140 6.251 0.502 0.617181
## 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 **
## TreatmentConNH:TimeSpring2018 -11.552 8.841 -1.307 0.195992
## TreatmentExclH:TimeSpring2018 -22.200 8.841 -2.511 0.014570 *
## TreatmentExclNH:TimeSpring2018 -14.774 8.841 -1.671 0.099575 .
## TreatmentConNH:Timespring2019 -7.914 8.841 -0.895 0.374041
## TreatmentExclH:Timespring2019 -18.982 8.841 -2.147 0.035577 *
## TreatmentExclNH:Timespring2019 -11.224 8.841 -1.270 0.208826
## TreatmentConNH:Timesummer2018 -5.278 8.841 -0.597 0.552602
## TreatmentExclH:Timesummer2018 -18.372 8.841 -2.078 0.041712 *
## TreatmentExclNH:Timesummer2018 -9.988 8.841 -1.130 0.262783
## ---
## 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
anova(model2)
## Analysis of Variance Table
##
## Response: Nmin_OM
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 382.6 127.52 1.3053 0.2804
## Time 3 8584.7 2861.56 29.2907 4.858e-12 ***
## Treatment:Time 9 822.8 91.42 0.9358 0.5008
## Residuals 64 6252.5 97.70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
resid2<-resid(model2)
hist(resid2)
model2aov<-aov(Nmin_OM~Treatment*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 ~ Treatment * 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.0000000
plot(posthoc2)
#Run linear mixed effects model setting "Site and Time" as random effects
#Also compare the AICs
library(lme4)
model2b<-lmer(Nmin_OM~Treatment + (1|Site),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model2b)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Nmin_OM ~ Treatment + (1 | Site)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 661.2 675.5 -324.6 649.2 74
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9558 -0.6375 -0.1031 0.5501 2.8230
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.9988 0.9994
## Residual 194.7511 13.9553
## Number of obs: 80, groups: Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.745 3.152 8.484
## TreatmentConNH -3.500 4.413 -0.793
## TreatmentExclH -1.587 4.413 -0.360
## TreatmentExclNH -5.857 4.413 -1.327
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH
## TretmntCnNH -0.700
## TrtmntExclH -0.700 0.500
## TrtmntExcNH -0.700 0.500 0.500
anova(model2b)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 382.57 127.52 0.6548
resid2b<-resid(model2b)
plot(resid2b)
model2ba<-lmer(Nmin_OM~Treatment + (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 ~ Treatment + (1 | Site/Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 628.4 645.1 -307.2 614.4 73
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1112 -0.5855 -0.1180 0.3914 2.9520
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time:Site (Intercept) 118.00 10.863
## Site (Intercept) 0.00 0.000
## Residual 77.75 8.818
## Number of obs: 80, groups: Time:Site, 20; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.745 3.128 8.549
## TreatmentConNH -3.500 2.788 -1.255
## TreatmentExclH -1.587 2.788 -0.569
## TreatmentExclNH -5.857 2.788 -2.100
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH
## TretmntCnNH -0.446
## TrtmntExclH -0.446 0.500
## TrtmntExcNH -0.446 0.500 0.500
anova(model2ba)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 382.57 127.52 1.6402
resid2ba<-resid(model2ba)
hist(resid2ba)
model2bb<-lmer(Nmin_OM~Treatment + (1|Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model2bb)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Nmin_OM ~ Treatment + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 614.3 628.6 -301.1 602.3 74
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0602 -0.5932 -0.1909 0.5983 3.5562
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 102.7 10.132
## Residual 93.1 9.649
## Number of obs: 80, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.745 5.506 4.857
## TreatmentConNH -3.500 3.051 -1.147
## TreatmentExclH -1.587 3.051 -0.520
## TreatmentExclNH -5.857 3.051 -1.919
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH
## TretmntCnNH -0.277
## TrtmntExclH -0.277 0.500
## TrtmntExcNH -0.277 0.500 0.500
anova(model2bb)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 382.57 127.52 1.3698
resid2bb<-resid(model2bb)
hist(resid2bb)
model2bc<-lmer(Nmin_OM~Treatment + 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 ~ Treatment + Time + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 603.6 625.1 -292.8 585.6 71
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1682 -0.6116 -0.1303 0.5967 3.6243
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 0.00 0.000
## Residual 88.44 9.404
## Number of obs: 80, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 25.832 2.782 9.286
## TreatmentConNH -3.500 2.974 -1.177
## TreatmentExclH -1.587 2.974 -0.533
## TreatmentExclNH -5.857 2.974 -1.969
## TimeSpring2018 6.184 2.974 2.080
## Timespring2019 -15.250 2.974 -5.128
## Timesummer2018 12.718 2.974 4.277
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH TrtENH TS2018 Tm2019
## TretmntCnNH -0.535
## TrtmntExclH -0.535 0.500
## TrtmntExcNH -0.535 0.500 0.500
## TmSprng2018 -0.535 0.000 0.000 0.000
## Tmsprng2019 -0.535 0.000 0.000 0.000 0.500
## Timsmmr2018 -0.535 0.000 0.000 0.000 0.500 0.500
anova(model2bc)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 382.6 127.52 1.4419
## Time 3 8584.7 2861.56 32.3555
anova(model2b,model2ba,model2bb,model2bc)
## Data: Deer_projectdata
## Models:
## model2b: Nmin_OM ~ Treatment + (1 | Site)
## model2bb: Nmin_OM ~ Treatment + (1 | Time)
## model2ba: Nmin_OM ~ Treatment + (1 | Site/Time)
## model2bc: Nmin_OM ~ Treatment + Time + (1 | Time)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model2b 6 661.16 675.45 -324.58 649.16
## model2bb 6 614.27 628.56 -301.14 602.27 46.890 0 < 2.2e-16 ***
## model2ba 7 628.43 645.10 -307.21 614.43 0.000 1 1
## model2bc 9 603.62 625.06 -292.81 585.62 28.812 2 5.542e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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~Treatment*Time, data=Deer_projectdata,na.action=na.omit)
summary(model3)
##
## Call:
## lm(formula = lgSMR_OM ~ Treatment * 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 ***
## TreatmentConNH -0.03400 0.09414 -0.361 0.71915
## TreatmentExclH 0.09600 0.09414 1.020 0.31166
## TreatmentExclNH 0.03800 0.09414 0.404 0.68780
## 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
## TreatmentConNH:TimeSpring2018 -0.01200 0.13313 -0.090 0.92846
## TreatmentExclH:TimeSpring2018 -0.15400 0.13313 -1.157 0.25166
## TreatmentExclNH:TimeSpring2018 -0.04400 0.13313 -0.331 0.74209
## TreatmentConNH:Timespring2019 0.04600 0.13313 0.346 0.73083
## TreatmentExclH:Timespring2019 -0.20800 0.13313 -1.562 0.12312
## TreatmentExclNH:Timespring2019 -0.24600 0.13313 -1.848 0.06925 .
## TreatmentConNH:Timesummer2018 0.10800 0.13313 0.811 0.42023
## TreatmentExclH:Timesummer2018 -0.14600 0.13313 -1.097 0.27689
## TreatmentExclNH:Timesummer2018 -0.08200 0.13313 -0.616 0.54011
## ---
## 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)
## Treatment 3 0.04406 0.01469 0.6630 0.5778
## Time 3 1.63387 0.54462 24.5839 1.063e-10 ***
## Treatment:Time 9 0.22521 0.02502 1.1295 0.3556
## Residuals 64 1.41784 0.02215
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
resid3<-resid(model3)
hist(resid3)
model3aov <- aov(lgSMR_OM~Treatment*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 ~ Treatment * 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.0000000
plot(posthoc3)
#Run linear mixed effects model setting "Site and Time" as random effects
#Also compare the AICs
library(lme4)
model3b<-lmer(SMR_OM~Treatment + (1|Site),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model3b)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: SMR_OM ~ Treatment + (1 | Site)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 749.0 763.3 -368.5 737.0 74
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3251 -0.6143 -0.1590 0.2994 6.3516
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 24.26 4.925
## Residual 567.83 23.829
## Number of obs: 80, groups: Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 44.569 5.766 7.730
## TreatmentConNH 4.516 7.535 0.599
## TreatmentExclH -0.865 7.535 -0.115
## TreatmentExclNH -2.740 7.535 -0.364
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH
## TretmntCnNH -0.653
## TrtmntExclH -0.653 0.500
## TrtmntExcNH -0.653 0.500 0.500
resid3b<-resid(model3b)
hist(resid3b)
model3ba<-lmer(SMR_OM~Treatment + (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 ~ Treatment + (1 | Site/Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 731.0 747.7 -358.5 717.0 73
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6168 -0.3443 -0.0654 0.1186 6.5122
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time:Site (Intercept) 2.799e+02 1.673e+01
## Site (Intercept) 5.208e-14 2.282e-07
## Residual 3.122e+02 1.767e+01
## Number of obs: 80, groups: Time:Site, 20; Site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 44.569 5.441 8.191
## TreatmentConNH 4.516 5.588 0.808
## TreatmentExclH -0.865 5.588 -0.155
## TreatmentExclNH -2.740 5.588 -0.490
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH
## TretmntCnNH -0.513
## TrtmntExclH -0.513 0.500
## TrtmntExcNH -0.513 0.500 0.500
resid3ba<-resid(model3ba)
hist(resid3ba)
model3bb<-lmer(SMR_OM~Treatment + (1|Time),na.action=na.omit,REML="FALSE",data=Deer_projectdata)
summary(model3bb)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: SMR_OM ~ Treatment + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 732.9 747.2 -360.4 720.9 74
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8009 -0.4992 -0.1123 0.2611 6.9564
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 161.0 12.69
## Residual 431.1 20.76
## Number of obs: 80, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 44.569 7.862 5.669
## TreatmentConNH 4.516 6.566 0.688
## TreatmentExclH -0.865 6.566 -0.132
## TreatmentExclNH -2.740 6.566 -0.417
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH
## TretmntCnNH -0.418
## TrtmntExclH -0.418 0.500
## TrtmntExcNH -0.418 0.500 0.500
resid3bb<-resid(model3bb)
hist(resid3bb)
anova(model3bb)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 568.77 189.59 0.4398
model3bc<-lmer(SMR_OM~Treatment + 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 ~ Treatment + Time + (1 | Time)
## Data: Deer_projectdata
##
## AIC BIC logLik deviance df.resid
## 726.2 747.7 -354.1 708.2 71
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9172 -0.5260 -0.0193 0.3024 7.0591
##
## Random effects:
## Groups Name Variance Std.Dev.
## Time (Intercept) 0.0 0.00
## Residual 409.5 20.24
## Number of obs: 80, groups: Time, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 56.482 5.986 9.436
## TreatmentConNH 4.516 6.399 0.706
## TreatmentExclH -0.865 6.399 -0.135
## TreatmentExclNH -2.740 6.399 -0.428
## TimeSpring2018 -17.826 6.399 -2.785
## Timespring2019 -31.273 6.399 -4.887
## Timesummer2018 1.450 6.399 0.227
##
## Correlation of Fixed Effects:
## (Intr) TrtCNH TrtmEH TrtENH TS2018 Tm2019
## TretmntCnNH -0.535
## TrtmntExclH -0.535 0.500
## TrtmntExcNH -0.535 0.500 0.500
## TmSprng2018 -0.535 0.000 0.000 0.000
## Tmsprng2019 -0.535 0.000 0.000 0.000 0.500
## Timsmmr2018 -0.535 0.000 0.000 0.000 0.500 0.500
anova(model3bc)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Treatment 3 568.8 189.6 0.463
## Time 3 14605.3 4868.4 11.888
anova(model3b,model3ba,model3bb,model3bc)
## Data: Deer_projectdata
## Models:
## model3b: SMR_OM ~ Treatment + (1 | Site)
## model3bb: SMR_OM ~ Treatment + (1 | Time)
## model3ba: SMR_OM ~ Treatment + (1 | Site/Time)
## model3bc: SMR_OM ~ Treatment + Time + (1 | Time)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model3b 6 748.98 763.27 -368.49 736.98
## model3bb 6 732.88 747.17 -360.44 720.88 16.1012 0 < 2e-16 ***
## model3ba 7 730.99 747.66 -358.49 716.99 3.8936 1 0.04847 *
## model3bc 9 726.23 747.67 -354.11 708.23 8.7561 2 0.01255 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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