FCI and Selection
Megan Wanchuk
2021-08-27
All Grazing Systems
Patch-Burn Grazing
Season Long Grazing
Modified Twice-over Rest Rotational Grazing
glmmTMB Models
Family Selection
##glmmTMB
###Fit family
nb1glmm <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
ziformula=~1,
family=nbinom1(),
data=FFQA)
nb1glmmM <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
ziformula=~MDS,
family=nbinom1(),
data=FFQA)
nb2glmm <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~1,
family=nbinom2())
nb2glmmM <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~MDS,
family=nbinom2())
tnb1glmm <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom1)
tnb1glmmz1 <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~1,
family=truncated_nbinom1)
tnb2glmm <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmmz1 <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~1,
family=truncated_nbinom2)
# Compare models with information criteria
bbmle::AICtab(nb1glmm, nb2glmm,
nb1glmmM, nb2glmmM,
tnb1glmm, tnb2glmm,
tnb1glmmz1, tnb2glmmz1)## dAIC df
## tnb2glmm 0.0 7
## tnb1glmm 57.7 7
## tnb2glmmz1 328.6 5
## nb2glmm 384.8 5
## tnb1glmmz1 386.2 5
## nb2glmmM 386.5 6
## nb1glmm 415.1 5
## nb1glmmM 416.7 6# Are top two models different?
anova(tnb2glmm, tnb1glmm)## Data: FFQA
## Models:
## tnb2glmm: PerFecalR ~ MDS + (1 | plot:Year:Month), zi=~., disp=~1
## tnb1glmm: PerFecalR ~ MDS + (1 | plot:Year:Month), zi=~., disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## tnb2glmm 7 14272 14318 -7129.2 14258
## tnb1glmm 7 14330 14376 -7158.1 14316 0 0 1#Model Summary for top two models
summary(tnb2glmm)## Family: truncated_nbinom2 ( log )
## Formula: PerFecalR ~ MDS + (1 | plot:Year:Month)
## Zero inflation: ~.
## Data: FFQA
##
## AIC BIC logLik deviance df.resid
## 14272.5 14318.3 -7129.2 14258.5 5158
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot:Year:Month (Intercept) 0.3604 0.6004
## Number of obs: 5165, groups: plot:Year:Month, 180
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## plot:Year:Month (Intercept) 0.8893 0.943
## Number of obs: 5165, groups: plot:Year:Month, 180
##
## Dispersion parameter for truncated_nbinom2 family (): 3.72
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.69046 0.05016 53.64 < 2e-16 ***
## MDS 0.65537 0.12000 5.46 4.73e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.25015 0.08059 15.513 <2e-16 ***
## MDS -1.51329 0.24627 -6.145 8e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(tnb2glmm)## 2.5 % 97.5 % Estimate
## cond.(Intercept) 2.5921579 2.7887631 2.6904605
## cond.MDS 0.4201668 0.8905640 0.6553654
## zi.(Intercept) 1.0922007 1.4081022 1.2501514
## zi.MDS -1.9959677 -1.0306189 -1.5132933
## plot:Year:Month.cond.Std.Dev.(Intercept) 0.5269861 0.6839453 0.6003579
## plot:Year:Month.zi.Std.Dev.(Intercept) 0.8083845 1.1001283 0.9430412summary(tnb1glmm)## Family: truncated_nbinom1 ( log )
## Formula: PerFecalR ~ MDS + (1 | plot:Year:Month)
## Zero inflation: ~.
## Data: FFQA
##
## AIC BIC logLik deviance df.resid
## 14330.1 14376.0 -7158.1 14316.1 5158
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot:Year:Month (Intercept) 0.4104 0.6407
## Number of obs: 5165, groups: plot:Year:Month, 180
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## plot:Year:Month (Intercept) 0.8893 0.9431
## Number of obs: 5165, groups: plot:Year:Month, 180
##
## Dispersion parameter for truncated_nbinom1 family (): 3.2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.69925 0.05215 51.76 < 2e-16 ***
## MDS 0.53842 0.11272 4.78 1.78e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.25014 0.08059 15.513 < 2e-16 ***
## MDS -1.51326 0.24627 -6.145 8.01e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(tnb1glmm)## 2.5 % 97.5 % Estimate
## cond.(Intercept) 2.5970385 2.8014658 2.6992522
## cond.MDS 0.3174847 0.7593495 0.5384171
## zi.(Intercept) 1.0921904 1.4080933 1.2501419
## zi.MDS -1.9959388 -1.0305907 -1.5132647
## plot:Year:Month.cond.Std.Dev.(Intercept) 0.5677556 0.7229183 0.6406566
## plot:Year:Month.zi.Std.Dev.(Intercept) 0.8083920 1.1001390 0.9430502Model Selection with All Treatments
Family truncated_nbinom2 was the best fit. - Test relationship between selection and forage quality in different grazing systems. Is MTRR a significant term? - Expect that the MTRR will be different than PBG and SL, and not a significant term in the model.
###Fit models with TNB2
tnb2glmm0 <- glmmTMB(PerFecalR ~ 1 + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmm1 <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmm2 <- glmmTMB(PerFecalR ~ MDS + Treatment + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom2)
bbmle::AICtab(tnb2glmm0, tnb2glmm1, tnb2glmm2)## dAIC df
## tnb2glmm2 0.0 11
## tnb2glmm1 65.6 7
## tnb2glmm0 130.5 5anova(tnb2glmm2, tnb2glmm1)## Data: FFQA
## Models:
## tnb2glmm1: PerFecalR ~ MDS + (1 | plot:Year:Month), zi=~., disp=~1
## tnb2glmm2: PerFecalR ~ MDS + Treatment + (1 | plot:Year:Month), zi=~., disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## tnb2glmm1 7 14272 14318 -7129.2 14258
## tnb2glmm2 11 14207 14279 -7092.4 14185 73.64 4 3.863e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(tnb2glmm2)## Family: truncated_nbinom2 ( log )
## Formula: PerFecalR ~ MDS + Treatment + (1 | plot:Year:Month)
## Zero inflation: ~.
## Data: FFQA
##
## AIC BIC logLik deviance df.resid
## 14206.8 14278.9 -7092.4 14184.8 5154
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot:Year:Month (Intercept) 0.2081 0.4562
## Number of obs: 5165, groups: plot:Year:Month, 180
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## plot:Year:Month (Intercept) 0.8349 0.9137
## Number of obs: 5165, groups: plot:Year:Month, 180
##
## Dispersion parameter for truncated_nbinom2 family (): 3.67
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.35396 0.06105 38.56 < 2e-16 ***
## MDS 0.81511 0.12119 6.73 1.75e-11 ***
## TreatmentSL 0.22014 0.09723 2.26 0.0236 *
## TreatmentMTRR 0.88843 0.09735 9.13 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0382 0.1183 8.774 < 2e-16 ***
## MDS -1.4239 0.2498 -5.700 1.19e-08 ***
## TreatmentSL 0.3580 0.1855 1.930 0.0536 .
## TreatmentMTRR 0.3766 0.1910 1.972 0.0486 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1glht_tnb2glmm2<-multcomp::glht(tnb2glmm2, linfct = mcp(Treatment =
"Tukey"))
summary(glht_tnb2glmm2)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmmTMB(formula = PerFecalR ~ MDS + Treatment + (1 | plot:Year:Month),
## data = FFQA, family = truncated_nbinom2, ziformula = ~.,
## dispformula = ~1)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## SL - PBG == 0 0.22014 0.09723 2.264 0.0606 .
## MTRR - PBG == 0 0.88843 0.09735 9.126 <1e-04 ***
## MTRR - SL == 0 0.66829 0.10606 6.301 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(glht_tnb2glmm2)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmmTMB(formula = PerFecalR ~ MDS + Treatment + (1 | plot:Year:Month),
## data = FFQA, family = truncated_nbinom2, ziformula = ~.,
## dispformula = ~1)
##
## Quantile = 2.3421
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## SL - PBG == 0 0.220137 -0.007595 0.447869
## MTRR - PBG == 0 0.888427 0.660423 1.116432
## MTRR - SL == 0 0.668290 0.419884 0.916696Model Selection- PBG and SL
- Issues in sampling with including MTRR
- Test relationship between poop counts and FCI for PBG and SL
##No MTRR?
FFQPS<- FFQA %>%
dplyr::filter(Treatment!="MTRR")
tnb2glmm0PS <- glmmTMB(PerFecalR ~ 1 + (1|plot:Year:Month),
data=FFQPS,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmm1PS <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQPS,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmm2PS <- glmmTMB(PerFecalR ~ MDS + Treatment + (1|plot:Year:Month),
data=FFQPS,
ziformula=~.,
family=truncated_nbinom2)
bbmle::AICtab(tnb2glmm0PS, tnb2glmm1PS, tnb2glmm2PS)## dAIC df
## tnb2glmm2PS 0.0 9
## tnb2glmm1PS 1.9 7
## tnb2glmm0PS 114.4 5anova(tnb2glmm2PS, tnb2glmm1PS)## Data: FFQPS
## Models:
## tnb2glmm1PS: PerFecalR ~ MDS + (1 | plot:Year:Month), zi=~., disp=~1
## tnb2glmm2PS: PerFecalR ~ MDS + Treatment + (1 | plot:Year:Month), zi=~., disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## tnb2glmm1PS 7 11312 11356 -5649.0 11298
## tnb2glmm2PS 9 11310 11367 -5646.1 11292 5.8871 2 0.05268 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(tnb2glmm2PS)## Family: truncated_nbinom2 ( log )
## Formula: PerFecalR ~ MDS + Treatment + (1 | plot:Year:Month)
## Zero inflation: ~.
## Data: FFQPS
##
## AIC BIC logLik deviance df.resid
## 11310.1 11367.0 -5646.1 11292.1 4117
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot:Year:Month (Intercept) 0.245 0.495
## Number of obs: 4126, groups: plot:Year:Month, 128
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## plot:Year:Month (Intercept) 1.343 1.159
## Number of obs: 4126, groups: plot:Year:Month, 128
##
## Dispersion parameter for truncated_nbinom2 family (): 3.76
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.36048 0.06532 36.14 < 2e-16 ***
## MDS 0.80428 0.13002 6.19 6.18e-10 ***
## TreatmentSL 0.21940 0.10400 2.11 0.0349 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1234 0.1454 7.728 1.09e-14 ***
## MDS -2.4145 0.2845 -8.487 < 2e-16 ***
## TreatmentSL 0.2844 0.2275 1.250 0.211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1glht_tnb2glmm2PS<-multcomp::glht(tnb2glmm2PS, linfct = mcp(Treatment =
"Tukey"))
summary(glht_tnb2glmm2PS)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmmTMB(formula = PerFecalR ~ MDS + Treatment + (1 | plot:Year:Month),
## data = FFQPS, family = truncated_nbinom2, ziformula = ~.,
## dispformula = ~1)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## SL - PBG == 0 0.2194 0.1040 2.11 0.0349 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(glht_tnb2glmm2PS)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glmmTMB(formula = PerFecalR ~ MDS + Treatment + (1 | plot:Year:Month),
## data = FFQPS, family = truncated_nbinom2, ziformula = ~.,
## dispformula = ~1)
##
## Quantile = 1.96
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## SL - PBG == 0 0.21940 0.01556 0.42324Model Prediction
All grazing systems
SL and PBG
Script
if (!require("pacman")) install.packages("pacman")
pacman::p_load(knitr)
opts_chunk$set(message = FALSE, warning=FALSE,
echo=FALSE, eval=TRUE, fig.path= 'Figs/')
pacman::p_load(tidyverse, glmmTMB, multcomp, ggstance)
##17 and 18 Samp Sheets
URL3 <- url("https://docs.google.com/spreadsheets/d/e/2PACX-1vSCzVK8I-x8Hr6zKD2oXg9dBcD4DqHb8gEMwi0wLPZ2Vi3Df8AzBZl4HKWgG3E3BRqRiA0Js25j_fHH/pub?output=csv")
Samps17 <- read.csv(URL3)
Samps17 <- mutate_at(Samps17,vars(ID, point, Year, Biomass, BurnSize, TSF, Treatment), as.character)
URL4 <- url("https://docs.google.com/spreadsheets/d/e/2PACX-1vThwrjzlDD6GA5d237StJ0zbxB5XeevZWpHwjc_waiObIDdFxndyVpIMHU3fLHn2q_lfZkHlGEih_9q/pub?output=csv")
Samps18 <- read.csv(URL4)
Samps18 <- mutate_at(Samps18,vars(ID, point, Year, Biomass, BurnSize, TSF, Treatment), as.character)
#Merge two sheets 2019and2020
URL <- url("https://docs.google.com/spreadsheets/d/e/2PACX-1vQrev81cM0-kEdPgZq3xT5-OjE-_Qv943JzemCpNxfhwx6MwwFb4S917-CPiBJY2M5h7rjz3Sg6oJox/pub?output=csv")
Samps19 <- read.csv(URL)
Samps19 <- mutate_at(Samps19,vars(ID, point, Year, Biomass, BurnSize, TSF, Treatment), as.character)
URL2 <- url("https://docs.google.com/spreadsheets/d/e/2PACX-1vSaDJebDNlThVT0-TIPfHXxUFDbR4AixfjLxRu8WjThBLB7Y36YBkhY8dnZzrhWd7Rx1_VooUQkYydS/pub?output=csv")
Samps20 <- read.csv(URL2)
Samps20 <- mutate_at(Samps20,vars(ID, point, Year, Biomass, BurnSize, TSF, Treatment), as.character)
Samps.17181920<-rbind(Samps17, Samps18, Samps19, Samps20)
Samps.17181920$Added_Column <- "S"
Samps.17181920 <- Samps.17181920 %>%
unite(col="ID",
c("Added_Column", "ID"), sep="", remove = T)%>%
mutate(Month=recode(Month, "August"="Aug"))
#add in Quality sheet
xl_file = "C:/Users/megan.wanchuk/Documents/CGREC/Forage/ForageSheetwithFixedSampleNo.xlsx"
xl_file2= "C:/Users/megan.wanchuk/Documents/CGREC/Forage/AllCGREC17-18.xlsx"
FQ1 <- readxl::read_excel(xl_file)
FQ1<-FQ1%>%
dplyr::select("Sample Number", ADF:N)
FQ2 <- readxl::read_excel(xl_file2)
FQ<- rbind(FQ1,FQ2)
FQ<- FQ%>%
rename( "ID"="Sample Number")
FFQ<-right_join(Samps.17181920, FQ, by ="ID")
FFQ<-FFQ%>%
mutate(TDNg= 98.625-(1.048*ADF)) %>%
mutate(TDNm= 92.62-(0.9093*ADF))
#Select the spring only PBG
FFQPB<- FFQ %>%
dplyr::filter(Treatment=="Spring Only")%>%
mutate( TSF= recode(TSF,
"2YSB"="2-3YSB",
"3YSB"="2-3YSB"
))%>%
mutate(Month= fct_relevel(Month,
"May", "June", "July", "Aug", "Sept"))
FFQA<- FFQ %>%
filter(Treatment!="Spring+Summer")%>%
mutate( TSF= recode(TSF,
"2YSB"="2-3YSB",
"3YSB"="2-3YSB"))%>%
mutate(Month= fct_relevel(Month,
"May", "June", "July", "Aug", "Sept"))%>%
mutate(Treatment= recode(Treatment,
"Spring Only"= "PBG", "MTR"="MTRR"))%>%
mutate(Treatment = fct_relevel(Treatment,
"PBG", "SL", "MTRR"))
##Calculate MDS1
FFQA1<- unite(data= FFQ, col="YP",
c( "Year", "plot"), sep=" ", remove = FALSE) %>%
filter(Treatment!="Spring+Summer")%>%
mutate(Month= fct_relevel(Month,
"May", "June", "July", "Aug", "Sept"))%>%
mutate( Treatment= recode(Treatment,
"Spring Only"="PBG",
"MTR"="MTRR"))
FFQA.2<- FFQA1%>%
dplyr::select('CP','ADF','NDF','ADL', 'TDNg')%>%
scale()
FQA_pca <- vegan::capscale(FFQA.2 ~ 1, "euc")
##Creating a dataframe of MDS1
MDS_1 <-
vegan::scores(FQA_pca, display = "site") %>%
as.data.frame %>%
dplyr::select('MDS1')
MDS_1<-MDS_1%>%
mutate(MDS=MDS1*1)%>%
dplyr::select(MDS)
FFQA<-cbind(FFQA, MDS_1)
##Summarize patch points and fecal counts
ForageAll <- FFQA%>%
group_by(Year, Month, plot)%>%
summarize(
TCount=sum(Fecal))%>%
ungroup()
ForageAll<- ForageAll%>%
drop_na()
FFQA<-left_join(FFQA, ForageAll, by=c("Year", "Month", "plot"))
##Calculuate SI
FFQA <- FFQA%>%
mutate(PFecal= Fecal/ TCount)
FFQA$PFecal[is.na(FFQA$PFecal)]<-0
#calc precentages
FFQA <- FFQA%>%
mutate(PerFecal= PFecal*100)
FFQA <- FFQA%>%
mutate(PerFecalR = round(PerFecal, 0))
##All Grazing systems
FCIxF<-lm(MDS ~ PFecal, data = FFQA)
summary(FCIxF)
ggplot(FFQA, aes(x = MDS, y = PFecal, colour=TSF)) +
geom_point() +
theme_bw()+
stat_smooth(method = "lm", col = "red")
ggplot(FFQA, aes(x = MDS, y = PFecal, colour=Month)) +
geom_point() +
theme_bw()+
stat_smooth(method = "lm", col = "red")
##PBG
FFQPBG<- FFQA %>%
dplyr::filter(Treatment=="PBG")
FCIxFPBG<-lm(MDS ~ PFecal, data = FFQPBG)
summary(FCIxFPBG)
ggplot(FFQPBG, aes(x = MDS, y = PFecal, colour=TSF)) +
geom_point() +
theme_bw()+
stat_smooth(method = "lm", col = "red")
ggplot(FFQPBG, aes(x = MDS, y = PFecal, colour=Month)) +
geom_point() +
theme_bw()+
stat_smooth(method = "lm", col = "red")
##SL
FFQSL<- FFQA %>%
dplyr::filter(Treatment=="SL")
FCIxFSL<-lm(MDS ~ PFecal, data = FFQSL)
summary(FCIxFSL)
ggplot(FFQSL, aes(x = MDS, y = PFecal, colour=Month)) +
geom_point() +
theme_bw()+
stat_smooth(method = "lm", col = "red")
#MTRR
FFQMTRR<- FFQA %>%
dplyr::filter(Treatment=="MTRR")
FCIxFMTRR<-lm(MDS ~ PFecal, data = FFQMTRR)
summary(FCIxFMTRR)
ggplot(FFQMTRR, aes(x = MDS, y = PFecal, colour=Month)) +
geom_point() +
theme_bw()+
stat_smooth(method = "lm", col = "red")
##glmmTMB
###Fit family
nb1glmm <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
ziformula=~1,
family=nbinom1(),
data=FFQA)
nb1glmmM <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
ziformula=~MDS,
family=nbinom1(),
data=FFQA)
nb2glmm <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~1,
family=nbinom2())
nb2glmmM <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~MDS,
family=nbinom2())
tnb1glmm <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom1)
tnb1glmmz1 <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~1,
family=truncated_nbinom1)
tnb2glmm <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmmz1 <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~1,
family=truncated_nbinom2)
# Compare models with information criteria
bbmle::AICtab(nb1glmm, nb2glmm,
nb1glmmM, nb2glmmM,
tnb1glmm, tnb2glmm,
tnb1glmmz1, tnb2glmmz1)
# Are top two models different?
anova(tnb2glmm, tnb1glmm)
#Model Summary for top two models
summary(tnb2glmm)
confint(tnb2glmm)
summary(tnb1glmm)
confint(tnb1glmm)
###Fit models with TNB2
tnb2glmm0 <- glmmTMB(PerFecalR ~ 1 + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmm1 <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmm2 <- glmmTMB(PerFecalR ~ MDS + Treatment + (1|plot:Year:Month),
data=FFQA,
ziformula=~.,
family=truncated_nbinom2)
bbmle::AICtab(tnb2glmm0, tnb2glmm1, tnb2glmm2)
anova(tnb2glmm2, tnb2glmm1)
summary(tnb2glmm2)
glht_tnb2glmm2<-multcomp::glht(tnb2glmm2, linfct = mcp(Treatment =
"Tukey"))
summary(glht_tnb2glmm2)
confint(glht_tnb2glmm2)
##No MTRR?
FFQPS<- FFQA %>%
dplyr::filter(Treatment!="MTRR")
tnb2glmm0PS <- glmmTMB(PerFecalR ~ 1 + (1|plot:Year:Month),
data=FFQPS,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmm1PS <- glmmTMB(PerFecalR ~ MDS + (1|plot:Year:Month),
data=FFQPS,
ziformula=~.,
family=truncated_nbinom2)
tnb2glmm2PS <- glmmTMB(PerFecalR ~ MDS + Treatment + (1|plot:Year:Month),
data=FFQPS,
ziformula=~.,
family=truncated_nbinom2)
bbmle::AICtab(tnb2glmm0PS, tnb2glmm1PS, tnb2glmm2PS)
anova(tnb2glmm2PS, tnb2glmm1PS)
summary(tnb2glmm2PS)
glht_tnb2glmm2PS<-multcomp::glht(tnb2glmm2PS, linfct = mcp(Treatment =
"Tukey"))
summary(glht_tnb2glmm2PS)
confint(glht_tnb2glmm2PS)
prediction1 <- ggeffects::ggpredict(tnb2glmm2,c( "MDS[all]" ,
"Treatment"))
ggplot() +
geom_boxploth(data = FFQA,
aes(y = 85, x = MDS, color=Treatment),
size=1, width= 5)+
geom_point( data = FFQA,
aes(x= MDS, y= PerFecalR, colour=Treatment),
alpha=0.4, shape= 15)+
geom_line(data= prediction1,
aes(x = x, y = predicted, colour= group), size=1) +
geom_ribbon( data= prediction1,
aes(x = x, y = predicted, ymin = conf.low,
ymax = conf.high,fill = group),
alpha = .2, colour = NA) +
theme_bw()+
guides(fill = "none")+
labs(x="FCI",
y="Predicted Poop Count %")
predictions <- ggeffects::ggpredict(tnb2glmm2PS,c( "MDS[all]" ,
"Treatment"))
ggplot() +
geom_point( data = FFQPS,
aes(x= MDS, y= PerFecalR, colour=Treatment),alpha=0.6)+
geom_line(data= predictions,
aes(x = x, y = predicted, colour= group) ,size=1) +
geom_ribbon( data= predictions,
aes(x = x, y = predicted, ymin = conf.low,
ymax = conf.high, fill = group),alpha = .2, colour =
NA) +
geom_boxploth(data = FFQPS, aes(y = 85,
x = MDS, colour=Treatment ),size=1,
width= 5)+
theme_bw()+
guides(fill="none")+
labs(x="FCI",
y="Predicted Poop Count %")