The numbers of prey caught by multiple cats in the same household were combined, because prey could not in every case be confidently attributed to an individual cat.
Description of the data frame: p_id: participant identification number. period: 1 refers to pre-treatment period, 3 refers to treatment period. assigned: group of treatment, C=control, B=Bell, F=Food, P=Puzzle, S=Safe, T=Toy (Object play). trial_day: date of trial. daily_prey: daily total number of prey brought home. mammals: daily number of mammals brought home. birds: daily number of birds brought home.
reptiles:daily number of reptiles brought home. amphibians: daily number of amphibians brought home. insects: daily number of insects brought home. unsure_prey: daily number of unidentifiable brought home. int_cats: daily number of cats following the treatment.During treatment period,daily prey records were excluded from analysis if the cats in the house were not all following the treatment on that day. n_cats: total number of cats in the house.
library(dplyr)
library(tidyverse)
library(ggplot2)
library(DHARMa)
library(magrittr)
library(MuMIn)
library(lme4)
Read dataframe
per_house <- read.csv("~/Dropbox/Cat Project Team Folder/Intervention Study/Journal paper/R2/Cecchetti household data.csv", header = TRUE, stringsAsFactors = FALSE)
Work out the dates where not all cats in house were complying with treatment and exclude these. Major part of cat owners considered 20/06 the last day of recording, so exclude 21/06. Set Control group as reference one.
per_house$period <- as.factor(per_house$period)
exclude <- subset(per_house, per_house$period == '3' & per_house$int_cats< per_house$n_cats)
per_house <- per_house %>% anti_join(exclude, by= c('p_id', 'trial_date'))
per_house$trial_date<-as.Date(per_house$trial_date)
per_house<- filter(per_house, !trial_date > '2019-06-20')
per_house$assigned<- as.factor(per_house$assigned)
per_house$assigned <- relevel(per_house$assigned, ref= 'C')
Create a new column “active” for counting days of owner effort. Group counts of prey by treatment group, period and household
per_house$active<- as.integer(1)
per_house %>%
group_by(p_id, period,assigned)%>%
summarise(total_prey= sum(daily_prey, na.rm=TRUE),
mammals = sum(mammals, na.rm=TRUE),
birds = sum(birds, na.rm=TRUE),
reptiles = sum(reptiles, na.rm=TRUE),
amphibians= sum(amphibians, na.rm=TRUE),
insects = sum(insects, na.rm=TRUE),
unsure_prey = sum(unsure_prey, na.rm=TRUE),
eff= sum(active, na.rm= TRUE)) ->df
#add number of cats in the house
df$n_cats <- per_house$n_cats[match(df$p_id, per_house$p_id)]
df$p_id <- as.factor(df$p_id)
head(df)
## # A tibble: 6 x 12
## # Groups: p_id, period [6]
## p_id period assigned total_prey mammals birds reptiles amphibians
## <fct> <fct> <fct> <int> <int> <int> <int> <int>
## 1 3 1 F 6 3 3 0 0
## 2 3 3 F 6 1 5 0 0
## 3 6 1 P 7 4 3 0 0
## 4 6 3 P 14 13 1 0 0
## 5 8 1 B 23 4 2 0 0
## 6 8 3 B 8 0 7 0 0
## # … with 4 more variables: insects <int>, unsure_prey <int>, eff <int>,
## # n_cats <int>
Sample sizes:
df%>% group_by(assigned, period) %>%
summarise(p_id= length(unique(p_id)),
total_prey= sum(total_prey, na.rm=TRUE),
mammals= sum(mammals, na.rm=TRUE),
birds=sum(birds, na.rm=TRUE),
reptiles= sum(reptiles, na.rm=TRUE),
amphibians= sum(amphibians, na.rm=TRUE),
insects= sum(insects, na.rm=TRUE),
unsure_prey= sum(unsure_prey, na.rm=TRUE),
n_cats = sum(n_cats, na.rm=TRUE)) ->cats_per_group
head(cats_per_group)
## # A tibble: 6 x 11
## # Groups: assigned [3]
## assigned period p_id total_prey mammals birds reptiles amphibians
## <fct> <fct> <int> <int> <int> <int> <int> <int>
## 1 C 1 32 361 267 43 14 2
## 2 C 3 32 271 185 55 7 2
## 3 B 1 33 486 316 92 23 2
## 4 B 3 33 298 175 96 9 1
## 5 F 1 39 576 418 79 28 0
## 6 F 3 39 271 188 57 5 0
## # … with 3 more variables: insects <int>, unsure_prey <int>, n_cats <int>
GLMM with a Poisson error distribution and log link to analyse variation in total numbers of prey as a function of treatment.
m1<- glmer(total_prey ~ assigned*period+ (1|p_id) + offset(log(n_cats*eff)),
control=glmerControl(optimizer="bobyqa",optCtrl = list(maxfun=20000)),
family= poisson(link="log"), data= df)
summary(m1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## total_prey ~ assigned * period + (1 | p_id) + offset(log(n_cats *
## eff))
## Data: df
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
##
## AIC BIC logLik deviance df.resid
## 2678.7 2731.7 -1326.3 2652.7 425
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1951 -0.7058 -0.1333 0.5190 4.4242
##
## Random effects:
## Groups Name Variance Std.Dev.
## p_id (Intercept) 0.63 0.7938
## Number of obs: 438, groups: p_id, 219
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.32743 0.15344 -15.168 < 2e-16 ***
## assignedB 0.42180 0.21340 1.977 0.04809 *
## assignedF 0.44198 0.20502 2.156 0.03110 *
## assignedP 0.31854 0.20839 1.529 0.12636
## assignedS 0.34052 0.21736 1.567 0.11720
## assignedT 0.12420 0.20137 0.617 0.53737
## period3 0.07988 0.08014 0.997 0.31884
## assignedB:period3 -0.14365 0.10864 -1.322 0.18608
## assignedF:period3 -0.45126 0.10877 -4.149 3.34e-05 ***
## assignedP:period3 0.28459 0.10952 2.598 0.00936 **
## assignedS:period3 -0.14877 0.11609 -1.281 0.20003
## assignedT:period3 -0.28198 0.11747 -2.400 0.01638 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) assgnB assgnF assgnP assgnS assgnT perid3 assB:3 assF:3
## assignedB -0.718
## assignedF -0.747 0.537
## assignedP -0.735 0.528 0.550
## assignedS -0.705 0.507 0.527 0.519
## assignedT -0.761 0.547 0.569 0.560 0.537
## period3 -0.217 0.156 0.162 0.160 0.153 0.165
## assgndB:pr3 0.160 -0.199 -0.120 -0.118 -0.113 -0.122 -0.738
## assgndF:pr3 0.160 -0.115 -0.197 -0.118 -0.113 -0.122 -0.737 0.543
## assgndP:pr3 0.159 -0.114 -0.119 -0.223 -0.112 -0.121 -0.732 0.540 0.539
## assgndS:pr3 0.150 -0.108 -0.112 -0.110 -0.205 -0.114 -0.690 0.509 0.509
## assgndT:pr3 0.148 -0.106 -0.111 -0.109 -0.104 -0.216 -0.682 0.503 0.503
## assP:3 assS:3
## assignedB
## assignedF
## assignedP
## assignedS
## assignedT
## period3
## assgndB:pr3
## assgndF:pr3
## assgndP:pr3
## assgndS:pr3 0.505
## assgndT:pr3 0.499 0.471
#R squared
r.squaredGLMM(m1)
## R2m R2c
## delta 0.05943678 0.8690170
## lognormal 0.05981085 0.8744863
## trigamma 0.05902934 0.8630599
Model check
simulationOutput <- simulateResiduals(fittedModel = m1, n=10000)
plot(simulationOutput)
plotResiduals(df$assigned, simulationOutput$scaledResiduals)
## NULL
plotResiduals(df$period, simulationOutput$scaledResiduals)
## NULL
QQ plot residuals and tests normal. Some deviations between residual distribution at higher model predictions, but no patterning in residuals between treatments groups or time periods.
testResiduals(simulationOutput)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.032727, p-value = 0.7362
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.74822, p-value = 0.105
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test
##
## data: simulationOutput
## outLow = 0.000e+00, outHigh = 0.000e+00, nobs = 4.380e+02, freqH0
## = 9.999e-05, p-value = 1
## alternative hypothesis: two.sided
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.032727, p-value = 0.7362
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.74822, p-value = 0.105
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test
##
## data: simulationOutput
## outLow = 0.000e+00, outHigh = 0.000e+00, nobs = 4.380e+02, freqH0
## = 9.999e-05, p-value = 1
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.1552, p-value = 0.4674
## alternative hypothesis: two.sided
GLMM with a Poisson error distribution and log link to analyse variation in numbers of mammals as a function of treatment.
m_mammals<- glmer(mammals ~ assigned*period + (1|p_id) + offset(log(n_cats*eff)),
control=glmerControl(optimizer="bobyqa",optCtrl = list(maxfun=20000)),
family= poisson(link="log"), data= df)
summary(m_mammals)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: mammals ~ assigned * period + (1 | p_id) + offset(log(n_cats *
## eff))
## Data: df
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
##
## AIC BIC logLik deviance df.resid
## 2476.9 2530.0 -1225.5 2450.9 425
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0686 -0.7657 -0.1938 0.5409 5.3021
##
## Random effects:
## Groups Name Variance Std.Dev.
## p_id (Intercept) 0.8981 0.9477
## Number of obs: 438, groups: p_id, 219
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.72439 0.18398 -14.808 < 2e-16 ***
## assignedB 0.24093 0.25729 0.936 0.34906
## assignedF 0.42245 0.24565 1.720 0.08549 .
## assignedP 0.25387 0.25027 1.014 0.31041
## assignedS 0.30202 0.26073 1.158 0.24671
## assignedT 0.19999 0.24075 0.831 0.40613
## period3 -0.00597 0.09512 -0.063 0.94996
## assignedB:period3 -0.17453 0.13337 -1.309 0.19065
## assignedF:period3 -0.39726 0.12918 -3.075 0.00210 **
## assignedP:period3 0.39804 0.12890 3.088 0.00201 **
## assignedS:period3 -0.07249 0.13963 -0.519 0.60365
## assignedT:period3 -0.43315 0.14044 -3.084 0.00204 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) assgnB assgnF assgnP assgnS assgnT perid3 assB:3 assF:3
## assignedB -0.712
## assignedF -0.747 0.533
## assignedP -0.732 0.524 0.548
## assignedS -0.704 0.503 0.527 0.517
## assignedT -0.761 0.544 0.570 0.560 0.537
## period3 -0.204 0.146 0.153 0.150 0.144 0.156
## assgndB:pr3 0.146 -0.190 -0.109 -0.107 -0.103 -0.111 -0.713
## assgndF:pr3 0.150 -0.108 -0.188 -0.111 -0.106 -0.115 -0.736 0.525
## assgndP:pr3 0.151 -0.108 -0.113 -0.214 -0.106 -0.115 -0.738 0.526 0.543
## assgndS:pr3 0.139 -0.099 -0.104 -0.102 -0.199 -0.106 -0.681 0.486 0.502
## assgndT:pr3 0.139 -0.099 -0.104 -0.102 -0.098 -0.195 -0.677 0.483 0.499
## assP:3 assS:3
## assignedB
## assignedF
## assignedP
## assignedS
## assignedT
## period3
## assgndB:pr3
## assgndF:pr3
## assgndP:pr3
## assgndS:pr3 0.503
## assgndT:pr3 0.500 0.461
r.squaredGLMM(m_mammals)
## R2m R2c
## delta 0.04833047 0.8697199
## lognormal 0.04874392 0.8771601
## trigamma 0.04786534 0.8613497
Model check
simulationOutput2 <- simulateResiduals(fittedModel = m_mammals, n=100000)
plot(simulationOutput2)
plotResiduals(df$assigned, simulationOutput2$scaledResiduals)
## NULL
plotResiduals(df$period, simulationOutput2$scaledResiduals)
## NULL
QQ plot residuals normal. Some small deviations between higher residual distribution at higher model predictions, but no patterning in residuals between treatments groups or between time periods. Dispersion test suggests data slightly underdispersed, but on the border of significance/non-significance so not too concerning.
testResiduals(simulationOutput2)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.041454, p-value = 0.439
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.66185, p-value = 0.0463
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test
##
## data: simulationOutput
## outLow = 0.0000e+00, outHigh = 0.0000e+00, nobs = 4.3800e+02,
## freqH0 = 9.9999e-06, p-value = 1
## alternative hypothesis: two.sided
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.041454, p-value = 0.439
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.66185, p-value = 0.0463
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test
##
## data: simulationOutput
## outLow = 0.0000e+00, outHigh = 0.0000e+00, nobs = 4.3800e+02,
## freqH0 = 9.9999e-06, p-value = 1
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput2)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.1755, p-value = 0.2215
## alternative hypothesis: two.sided
GLMM with a Poisson error distribution and log link to analyse variation in numbers of birds as a function of treatment.
m_birds<- glmer(birds ~ assigned*period + (1|p_id) + offset(log(n_cats*eff)),
control=glmerControl(optimizer="bobyqa",optCtrl = list(maxfun=20000)),
family= poisson(link="log"), data= df)
summary(m_birds)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: birds ~ assigned * period + (1 | p_id) + offset(log(n_cats *
## eff))
## Data: df
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
##
## AIC BIC logLik deviance df.resid
## 1500.5 1553.6 -737.3 1474.5 425
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4431 -0.7242 -0.1898 0.5198 3.8098
##
## Random effects:
## Groups Name Variance Std.Dev.
## p_id (Intercept) 0.6663 0.8163
## Number of obs: 438, groups: p_id, 219
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.59792 0.22435 -20.495 < 2e-16 ***
## assignedB 0.89161 0.29113 3.063 0.00219 **
## assignedF 0.66784 0.28683 2.328 0.01989 *
## assignedP 0.80493 0.29170 2.759 0.00579 **
## assignedS 0.93798 0.29577 3.171 0.00152 **
## assignedT 0.19491 0.29631 0.658 0.51066
## period3 0.62946 0.20417 3.083 0.00205 **
## assignedB:period3 -0.12992 0.25177 -0.516 0.60585
## assignedF:period3 -0.57728 0.26892 -2.147 0.03182 *
## assignedP:period3 -0.27116 0.27586 -0.983 0.32563
## assignedS:period3 -0.54942 0.27591 -1.991 0.04645 *
## assignedT:period3 -0.04929 0.28259 -0.174 0.86155
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) assgnB assgnF assgnP assgnS assgnT perid3 assB:3 assF:3
## assignedB -0.757
## assignedF -0.771 0.593
## assignedP -0.760 0.583 0.593
## assignedS -0.750 0.575 0.584 0.575
## assignedT -0.742 0.575 0.583 0.573 0.565
## period3 -0.506 0.390 0.395 0.389 0.384 0.383
## assgndB:pr3 0.410 -0.463 -0.321 -0.315 -0.311 -0.310 -0.811
## assgndF:pr3 0.383 -0.296 -0.466 -0.295 -0.291 -0.291 -0.759 0.616
## assgndP:pr3 0.375 -0.288 -0.293 -0.478 -0.284 -0.283 -0.740 0.600 0.562
## assgndS:pr3 0.374 -0.288 -0.293 -0.288 -0.448 -0.283 -0.740 0.600 0.562
## assgndT:pr3 0.366 -0.281 -0.286 -0.281 -0.277 -0.511 -0.723 0.586 0.549
## assP:3 assS:3
## assignedB
## assignedF
## assignedP
## assignedS
## assignedT
## period3
## assgndB:pr3
## assgndF:pr3
## assgndP:pr3
## assgndS:pr3 0.548
## assgndT:pr3 0.535 0.535
r.squaredGLMM(m_birds)
## R2m R2c
## delta 0.09588806 0.5847098
## lognormal 0.10493360 0.6398680
## trigamma 0.08422070 0.5135641
Model check
simulationOutput3 <- simulateResiduals(fittedModel = m_birds, n=100000)
plot(simulationOutput3)
plotResiduals(df$assigned, simulationOutput3$scaledResiduals)
## NULL
plotResiduals(df$period, simulationOutput3$scaledResiduals)
## NULL
QQ plot residuals and tests normal. Expected distribution of Residuals at all model predictions.
testResiduals(simulationOutput3)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.027978, p-value = 0.8828
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.83872, p-value = 0.3414
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test
##
## data: simulationOutput
## outLow = 0.0000e+00, outHigh = 0.0000e+00, nobs = 4.3800e+02,
## freqH0 = 9.9999e-06, p-value = 1
## alternative hypothesis: two.sided
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.027978, p-value = 0.8828
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.83872, p-value = 0.3414
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test
##
## data: simulationOutput
## outLow = 0.0000e+00, outHigh = 0.0000e+00, nobs = 4.3800e+02,
## freqH0 = 9.9999e-06, p-value = 1
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput3)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99558, p-value = 0.9868
## alternative hypothesis: two.sided
Compute confidence intervals of the independent variable and interaction estimates for each model.
x <- confint(m1)[2:13 ,]
x1 <- confint(m_mammals)[2:13 ,]
x2 <- confint(m_birds)[2:13 ,]
Create three different dataframes and exponentiating values to obtain rate ratios with 95% confidence intervals (reported in table S2).
test_coefs <- cbind(fixef(m1), x) %>%
set_colnames(c("Estimate", "LCI", "UCI")) %>%
as.data.frame() %>%
rownames_to_column("Parameter") %>%
filter(Parameter != "(Intercept)")%>%
mutate_if(is.numeric, exp) %>%
mutate(model= "Total prey")
test_coefs2 <- cbind(fixef(m_mammals), x1) %>%
set_colnames(c("Estimate", "LCI", "UCI")) %>%
as.data.frame() %>%
rownames_to_column("Parameter") %>%
filter(Parameter != "(Intercept)")%>%
mutate_if(is.numeric, exp) %>%
mutate(model= "Mammals")
test_coefs3 <- cbind(fixef(m_birds), x2) %>%
set_colnames(c("Estimate", "LCI", "UCI")) %>%
as.data.frame() %>%
rownames_to_column("Parameter") %>%
filter(Parameter != "(Intercept)")%>%
mutate_if(is.numeric, exp) %>%
mutate(model= "Birds")
Bind the dataframes, rename variables and filter main effect of treatments.
rr <- rbind(test_coefs, test_coefs2, test_coefs3)
colnames(rr)[2] <- "RR"
colnames(rr)[5] <- "model"
rr <- filter(rr, !Parameter == 'assignedB')
rr <- filter(rr, !Parameter == 'assignedF')
rr <- filter(rr, !Parameter == 'assignedP')
rr <- filter(rr, !Parameter == 'assignedS')
rr <- filter(rr, !Parameter == 'assignedT')
rr[2:4] <- round(rr[2:4], 3)
knitr::kable(rr, format="markdown")
Parameter | RR | LCI | UCI | model |
---|---|---|---|---|
period3 | 1.083 | 0.924 | 1.269 | Total prey |
assignedB:period3 | 0.866 | 0.699 | 1.074 | Total prey |
assignedF:period3 | 0.637 | 0.513 | 0.790 | Total prey |
assignedP:period3 | 1.329 | 1.070 | 1.651 | Total prey |
assignedS:period3 | 0.862 | 0.684 | 1.084 | Total prey |
assignedT:period3 | 0.754 | 0.597 | 0.952 | Total prey |
period3 | 0.994 | 0.822 | 1.200 | Mammals |
assignedB:period3 | 0.840 | 0.644 | 1.094 | Mammals |
assignedF:period3 | 0.672 | 0.520 | 0.869 | Mammals |
assignedP:period3 | 1.489 | 1.153 | 1.925 | Mammals |
assignedS:period3 | 0.930 | 0.704 | 1.227 | Mammals |
assignedT:period3 | 0.648 | 0.490 | 0.857 | Mammals |
period3 | 1.877 | 1.260 | 2.813 | Birds |
assignedB:period3 | 0.878 | 0.535 | 1.437 | Birds |
assignedF:period3 | 0.561 | 0.330 | 0.949 | Birds |
assignedP:period3 | 0.762 | 0.443 | 1.307 | Birds |
assignedS:period3 | 0.577 | 0.335 | 0.988 | Birds |
assignedT:period3 | 0.952 | 0.546 | 1.656 | Birds |
Order model names and rename variables for the plot.
model <- rev(c('Birds', 'Mammals', 'Total_prey'))
varorder <- rev(c('period3', 'assignedB:period3', 'assignedF:period3', 'assignedP:period3',
'assignedS:period3', 'assignedT:period3'))
varnames <-rev(c('Treatment period', 'Bell:Treatment period', 'Food:Treatment period', 'Puzzle:Treatment period',
'Safe:Treatment period', 'Play:Treatment period'))
rr$model <- factor(rr$model, levels= c('Total prey', 'Mammals', 'Birds'))
Create plot showing rate ratios with 95% confidence intervals.
rr %>%
ggplot(aes(x = Parameter)) +
geom_point(aes(y = RR, color = model, shape = model, fill = model),
position=position_dodge(width=0.5), size = 1.5) +
geom_errorbar(aes(ymin = LCI, ymax = UCI, color = model), width = 0.4,
position=position_dodge(width=0.5), size = 1) +
geom_hline(yintercept = 1,linetype = "dashed") +
guides(colour = guide_legend(reverse = T), shape=FALSE, fill=FALSE)+
scale_y_continuous(breaks = c(0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8,2,2.2,2.4,2.6, 2.8,3))+
xlab('Explanatory variables\n')+
ylab('Rate Ratio\n')+
theme_classic()+
theme(axis.text = element_text(size=12, colour = 'black'),
axis.title = element_text(size=12),
axis.line = element_line(colour = "grey"),
legend.position = c(0.85,0.5),
legend.text = element_text(size= 12),
legend.title = element_blank())+
scale_x_discrete(limits = varorder, labels = varnames, expand=c(0,0.1))+
coord_flip()