#### Human shields mediate sexual conflict in a top predator ### load the required packages library(MuMIn) library(lme4) ### load the data procdata <- read.table('HS_data_24052016.txt', header = TRUE, sep = '\t', dec = '.') str(procdata) procdata$young <- as.factor(procdata$young) procdata$treerichbo <- as.factor(procdata$treerichbo) procdata$old <- as.factor(procdata$old) procdata$midage <- as.factor(procdata$midage) procdata$bogs <- as.factor(procdata$bogs) procdata$clearcut <- as.factor(procdata$clearcut) procdata$SURV <- as.factor(procdata$SURV) ### Run the candidate models mfullr <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + young * SURV + treerichbo * SURV + old * SURV + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) mfullr1 <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + young + treerichbo + old * SURV + midage + bogs + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) mfullr2 <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc + d_habitation_sc * SURV + young * SURV + treerichbo + old * SURV + midage * SURV + bogs + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) mfullr3 <- glmer(RESP~ d_road_sc + d_forest_r_sc + d_habitation_sc * SURV + young + treerichbo + old * SURV + midage + bogs + clearcut + NDVI2007_r_sc + (1|bear_coded) +(1|year), data = procdata, family = binomial) mfullr4 <- glmer(RESP~ d_road_sc + d_forest_r_sc + d_habitation_sc + young + treerichbo + old + midage + bogs + clearcut + NDVI2007_r_sc + (1|bear_coded) +(1|year), data = procdata, family = binomial) mhab <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) mhab2 <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc + d_habitation_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) mhab3 <- glmer(RESP~ d_road_sc + d_forest_r_sc + d_habitation_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) mhab4 <- glmer(RESP~ d_road_sc + d_forest_r_sc + d_habitation_sc + (1|bear_coded) +(1|year), data = procdata, family = binomial) mLC <- glmer(RESP~ young * SURV + treerichbo * SURV + old * SURV + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) mLC2 <- glmer(RESP~ young * SURV + treerichbo * SURV + old * SURV + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc + (1|bear_coded) +(1|year), data = procdata, family = binomial) mLC3 <- glmer(RESP~ young * SURV + treerichbo + old * SURV + midage * SURV + bogs + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) mLC4 <- glmer(RESP~ young + treerichbo + old * SURV + midage + bogs + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) mLC5 <- glmer(RESP~ young + treerichbo + old + midage + bogs + clearcut + NDVI2007_r_sc + (1|bear_coded) +(1|year), data = procdata, family = binomial) null <- mLC <- glmer(RESP~ 1 + (1|bear_coded) +(1|year), data = procdata, family = binomial) ### use the model.sel function to obtain the model selection table model.sel(mfullr, mfullr1, mfullr2, mfullr3, mfullr4, mhab, mhab2, mhab3, mhab4, mLC, mLC2, mLC3, mLC4, mLC5, null) ### the full model has the lowest AICc score and all other candidates are inconclusive summary(mfullr) ### Model output (note that we reversed the sign of the 'distance to...' variables in the paper to facilitate interpretation) # AIC BIC logLik deviance df.resid # 55806.5 56014.8 -27879.2 55758.5 43458 #Scaled residuals: # Min 1Q Median 3Q Max #-3.7345 -0.9250 0.1638 0.8816 4.6459 #Random effects: # Groups Name Variance Std.Dev. # bear_coded (Intercept) 0.151230 0.38888 # year (Intercept) 0.002164 0.04652 #Number of obs: 43482, groups: bear_coded, 22; year, 7 #Fixed effects: # Estimate Std. Error z value Pr(>|z|) #(Intercept) -0.32320 0.10863 -2.975 0.002927 ** #d_road_sc 0.55723 0.03098 17.985 < 2e-16 *** #SURV1 0.02292 0.08693 0.264 0.792069 #d_forest_r_sc 0.42022 0.02547 16.496 < 2e-16 *** #d_habitation_sc 0.20985 0.02178 9.636 < 2e-16 *** #young1 0.50383 0.07435 6.776 1.23e-11 *** #treerichbo1 0.45246 0.13467 3.360 0.000780 *** #old1 0.67485 0.07332 9.204 < 2e-16 *** #midage1 0.34764 0.06609 5.260 1.44e-07 *** #bogs1 -0.76125 0.09049 -8.413 < 2e-16 *** #clearcut1 -0.50123 0.10089 -4.968 6.76e-07 *** #NDVI2007_r_sc 0.18833 0.01957 9.625 < 2e-16 *** #d_road_sc:SURV1 -0.21938 0.03598 -6.098 1.07e-09 *** #SURV1:d_forest_r_sc -0.07133 0.02998 -2.379 0.017348 * #SURV1:d_habitation_sc -0.73405 0.02685 -27.344 < 2e-16 *** #SURV1:young1 -0.06706 0.09040 -0.742 0.458191 #SURV1:treerichbo1 -0.34129 0.16440 -2.076 0.037891 * #SURV1:old1 -0.48461 0.09001 -5.384 7.28e-08 *** #SURV1:midage1 0.20978 0.08072 2.599 0.009352 ** #SURV1:bogs1 -0.39779 0.11360 -3.502 0.000462 *** #SURV1:clearcut1 1.06127 0.11658 9.103 < 2e-16 *** #SURV1:NDVI2007_r_sc -0.00650 0.02385 -0.273 0.785223 ### AICcdiff scores, compare the full model with the models that have the survival interaction ### term on a specific model term excluded # d_roads mfullroads <- glmer(RESP~ d_road_sc + d_forest_r_sc * SURV + d_habitation_sc * SURV + young * SURV + treerichbo * SURV + old * SURV + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) # d_forest_rds mfullforestroads <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc + d_habitation_sc * SURV + young * SURV + treerichbo * SURV + old * SURV + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) # d_habitation mfullhabitation <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc + young * SURV + treerichbo * SURV + old * SURV + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) # young forest mfullyoung <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + young + treerichbo * SURV + old * SURV + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) # Tree rich bog mfulltrb <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + young * SURV + treerichbo + old * SURV + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) # old mfullold <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + young * SURV + treerichbo * SURV + old + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) # midage mfullmidage <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + young * SURV + treerichbo * SURV + old * SURV + midage + bogs * SURV + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) # bog mfullbog <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + young * SURV + treerichbo * SURV + old * SURV + midage * SURV + bogs + clearcut * SURV + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) # clearcut mfullclearcut <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + young * SURV + treerichbo * SURV + old * SURV + midage * SURV + bogs * SURV + clearcut + NDVI2007_r_sc * SURV + (1|bear_coded) +(1|year), data = procdata, family = binomial) # ndvi mfullndvi <- glmer(RESP~ d_road_sc * SURV + d_forest_r_sc * SURV + d_habitation_sc * SURV + young * SURV + treerichbo * SURV + old * SURV + midage * SURV + bogs * SURV + clearcut * SURV + NDVI2007_r_sc + (1|bear_coded) +(1|year), data = procdata, family = binomial) ### the model.sel table reveals the AICcdiff scores for each model term in relation to the full model # the terms for which the interaction with litter fate improved the model (note that some terms here have delta AICc values below 4 and were thus considered inconclusive in the paper) model.sel(mfullr, mfullroads, mfullhabitation, mfullold, mfullmidage, mfullbog, mfullclearcut) # the terms for which the interaction with litter fate was not or less important in relation to the full model model.sel(mfullr, mfullyoung, mfullndvi)