Per Household

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>

Models

All prey

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

Mammals

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

Birds

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

Rate Ratios and Forest Plot

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()