> #################### GENERALIZED MIXED EFFECT MODEL > summary (rabbits2) femaleID litterID weight parasites density disturbance color N_offspring habitat n1003 : 1 n79 : 4 Min. :115.0 Min. :2.000 Min. : 13.00 Min. : 15.03 black:29 Min. : 0.000 a:39 n1023 : 1 n105 : 3 1st Qu.:186.2 1st Qu.:4.000 1st Qu.: 30.73 1st Qu.: 68.71 grey :28 1st Qu.: 2.750 b:45 n1043 : 1 n5 : 3 Median :209.5 Median :5.000 Median : 42.02 Median :167.84 white:27 Median : 6.000 n1044 : 1 n54 : 3 Mean :209.4 Mean :4.845 Mean : 51.96 Mean :172.13 Mean : 6.262 n1045 : 1 n10 : 2 3rd Qu.:234.5 3rd Qu.:6.000 3rd Qu.: 65.32 3rd Qu.:287.47 3rd Qu.: 8.250 n1048 : 1 n11 : 2 Max. :319.0 Max. :9.000 Max. :138.40 Max. :381.08 Max. :26.000 (Other):78 (Other):67 > offspring1 <- glmer (N_offspring ~ scale(weight) + color + scale(disturbance) + (1|litterID), family = poisson, rabbits2) > summary (offspring1) Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod'] Family: poisson ( log ) Formula: N_offspring ~ scale(weight) + color + scale(disturbance) + (1 | litterID) Data: rabbits2 AIC BIC logLik deviance df.resid 542.3 556.9 -265.1 530.3 78 Scaled residuals: Min 1Q Median 3Q Max -2.6495 -1.2111 0.1404 0.6732 2.8998 Random effects: Groups Name Variance Std.Dev. litterID (Intercept) 0.2508 0.5008 Number of obs: 84, groups: litterID, 58 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.15565 0.14676 14.688 < 2e-16 *** scale(weight) 0.15675 0.08588 1.825 0.067965 . color[T.grey] -0.78080 0.23852 -3.274 0.001062 ** color[T.white] -0.81977 0.21222 -3.863 0.000112 *** scale(disturbance) -0.33236 0.09411 -3.532 0.000413 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) scl(w) clr[T.g] clr[T.w] scale(wght) 0.315 colr[T.gry] -0.712 -0.443 colr[T.wht] -0.671 -0.204 0.505 scl(dstrbn) -0.117 0.094 0.287 0.108 > ##### MODEL SELECTION IN MuMIn package > options(na.action = na.fail) > res7 <- dredge(offspring1, rank="AICc", REML=F) > est7 <- model.avg(res7, subset=delta<7, revised.var=T) > summary(est7) Call: model.avg(object = res7, subset = delta < 7, revised.var = T) Component model call: glmer(formula = N_offspring ~ <2 unique rhs>, data = rabbits2, family = poisson) Component models: df logLik AICc delta weight 123 6 -265.15 543.39 0.00 0.63 12 5 -266.83 544.43 1.04 0.37 Term codes: color scale(disturbance) scale(weight) 1 2 3 Model-averaged coefficients: (full average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 2.12173 0.15103 0.15320 13.850 < 2e-16 *** color[T.grey] -0.71147 0.24714 0.25051 2.840 0.004510 ** color[T.white] -0.79242 0.21402 0.21729 3.647 0.000266 *** scale(disturbance) -0.33876 0.09451 0.09598 3.530 0.000416 *** scale(weight) 0.09836 0.10184 0.10256 0.959 0.337517 (conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 2.12173 0.15103 0.15320 13.850 < 2e-16 *** color[T.grey] -0.71147 0.24714 0.25051 2.840 0.004510 ** color[T.white] -0.79242 0.21402 0.21729 3.647 0.000265 *** scale(disturbance) -0.33876 0.09451 0.09598 3.530 0.000416 *** scale(weight) 0.15675 0.08588 0.08723 1.797 0.072348 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Relative variable importance: color scale(disturbance) scale(weight) Importance: 1.00 1.00 0.63 N containing models: 2 2 1 > ############### LIKELIHOOD RATIO TESTS library(foreign, pos=23) library(survival, pos=23) library(MASS, pos=23) library(epiDisplay, pos=23) > ### WEIGHT EFFECT > offspring2 <- glmer (N_offspring ~ color + scale(disturbance) + (1|litterID), family = poisson, rabbits2) > lrtest (offspring1, offspring2) Likelihood ratio test Model 1: N_offspring ~ scale(weight) + color + scale(disturbance) + (1 | litterID) Model 2: N_offspring ~ color + scale(disturbance) + (1 | litterID) #Df LogLik Df Chisq Pr(>Chisq) 1 6 -265.15 2 5 -266.83 -1 3.3648 0.06661 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ### COLOR effect > offspring3 <- glmer (N_offspring ~ scale(weight) + scale(disturbance) + (1|litterID), family = poisson, rabbits2) > lrtest (offspring1, offspring3) Likelihood ratio test Model 1: N_offspring ~ scale(weight) + color + scale(disturbance) + (1 | litterID) Model 2: N_offspring ~ scale(weight) + scale(disturbance) + (1 | litterID) #Df LogLik Df Chisq Pr(>Chisq) 1 6 -265.15 2 4 -272.86 -2 15.427 0.0004467 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ### DISTURBANCE effect > offspring4 <- glmer (N_offspring ~ scale(weight) + color + (1|litterID), family = poisson, rabbits2) > lrtest (offspring1, offspring4) Likelihood ratio test Model 1: N_offspring ~ scale(weight) + color + scale(disturbance) + (1 | litterID) Model 2: N_offspring ~ scale(weight) + color + (1 | litterID) #Df LogLik Df Chisq Pr(>Chisq) 1 6 -265.15 2 5 -270.49 -1 10.674 0.001086 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > #### TUKEY POST HOC COMPARISON library(emmeans, pos=27) > lsmeans(offspring1, pairwise ~ color, adjust="tukey") $lsmeans color lsmean SE df asymp.LCL asymp.UCL black 2.16 0.147 Inf 1.87 2.44 grey 1.37 0.169 Inf 1.04 1.71 white 1.34 0.157 Inf 1.03 1.64 Results are given on the log (not the response) scale. Confidence level used: 0.95 $contrasts contrast estimate SE df z.ratio p.value black - grey 0.781 0.239 Inf 3.274 0.0030 black - white 0.820 0.212 Inf 3.863 0.0003 grey - white 0.039 0.225 Inf 0.173 0.9836 Results are given on the log (not the response) scale. P value adjustment: tukey method for comparing a family of 3 estimates ########################################################### #####################SCRIPT############################### #################### GENERALIZED MIXED EFFECT MODEL summary (rabbits2) offspring1 <- glmer (N_offspring ~ scale(weight) + color + scale(disturbance) + (1|litterID), family = poisson, rabbits2) summary (offspring1) ###### MODEL SELECTION IN MuMIn package options(na.action = na.fail) res7 <- dredge(offspring1, rank="AICc", REML=F) est7 <- model.avg(res7, subset=delta<7, revised.var=T) summary(est7) ############### LIKELIHOOD RATIO TESTS library(foreign, pos=23) library(survival, pos=23) library(MASS, pos=23) library(epiDisplay, pos=23) ### WEIGHT EFFECT offspring2 <- glmer (N_offspring ~ color + scale(disturbance) + (1|litterID), family = poisson, rabbits2) lrtest (offspring1, offspring2) ### COLOR effect offspring3 <- glmer (N_offspring ~ scale(weight) + scale(disturbance) + (1|litterID), family = poisson, rabbits2) lrtest (offspring1, offspring3) ### DISTURBANCE effect offspring4 <- glmer (N_offspring ~ scale(weight) + color + (1|litterID), family = poisson, rabbits2) lrtest (offspring1, offspring4) #### TUKEY POST HOC COMPARISON library(emmeans, pos=27) lsmeans(offspring1, pairwise ~ color, adjust="tukey")