########################################### ########## MIXED EFFECT LINEAR MODELS > summary (lizards1) Maternal.ID Egg.ID Oviposition.day Month Habitat Egg.Mass Snout.Vent.Length Tail.Length Hatchling.Mass Sex Survival.to.3.weeks weeks.at.racing Mass.at.racing Speed.over.25.cm Stops.while.racing Min. : 1.00 Min. : 4.00 Min. :45.00 April:32 Open :45 Min. :0.1276 Min. :17.00 Min. :26.00 Min. :0.1348 F:28 N:10 Min. :1.000 Min. :0.1062 Min. :0.0523 Min. : 0.400 1st Qu.:11.00 1st Qu.:21.00 1st Qu.:52.00 July :47 Shade:34 1st Qu.:0.1633 1st Qu.:18.00 1st Qu.:30.00 1st Qu.:0.1579 M:51 Y:69 1st Qu.:1.000 1st Qu.:0.1550 1st Qu.:0.2357 1st Qu.: 3.200 Median :32.00 Median :44.00 Median :66.00 Median :0.1760 Median :18.00 Median :31.00 Median :0.1720 Median :1.000 Median :0.1698 Median :0.3798 Median : 4.600 Mean :42.29 Mean :46.06 Mean :66.58 Mean :0.1816 Mean :18.42 Mean :30.92 Mean :0.1755 Mean :1.911 Mean :0.1772 Mean :0.5001 Mean : 5.256 3rd Qu.:84.00 3rd Qu.:69.50 3rd Qu.:78.50 3rd Qu.:0.2006 3rd Qu.:19.00 3rd Qu.:32.00 3rd Qu.:0.1924 3rd Qu.:3.000 3rd Qu.:0.2029 3rd Qu.:0.5301 3rd Qu.: 6.700 Max. :95.00 Max. :96.00 Max. :89.00 Max. :0.2480 Max. :21.00 Max. :36.00 Max. :0.2315 Max. :3.000 Max. :0.2945 Max. :3.2821 Max. :15.600 > lizards1 <- within(lizards1, { + Egg.ID <- as.factor(Egg.ID) + }) library(lme4, pos=28) > lizards1 <- within(lizards1, { + weeks.at.racing <- as.factor(weeks.at.racing) + }) > summary (lizards1) Maternal.ID Egg.ID Oviposition.day Month Habitat Egg.Mass Snout.Vent.Length Tail.Length Hatchling.Mass Sex Survival.to.3.weeks weeks.at.racing Mass.at.racing Speed.over.25.cm Stops.while.racing Min. : 1.00 4 : 2 Min. :45.00 April:32 Open :45 Min. :0.1276 Min. :17.00 Min. :26.00 Min. :0.1348 F:28 N:10 1:43 Min. :0.1062 Min. :0.0523 Min. : 0.400 1st Qu.:11.00 5 : 2 1st Qu.:52.00 July :47 Shade:34 1st Qu.:0.1633 1st Qu.:18.00 1st Qu.:30.00 1st Qu.:0.1579 M:51 Y:69 3:36 1st Qu.:0.1550 1st Qu.:0.2357 1st Qu.: 3.200 Median :32.00 7 : 2 Median :66.00 Median :0.1760 Median :18.00 Median :31.00 Median :0.1720 Median :0.1698 Median :0.3798 Median : 4.600 Mean :42.29 8 : 2 Mean :66.58 Mean :0.1816 Mean :18.42 Mean :30.92 Mean :0.1755 Mean :0.1772 Mean :0.5001 Mean : 5.256 3rd Qu.:84.00 9 : 2 3rd Qu.:78.50 3rd Qu.:0.2006 3rd Qu.:19.00 3rd Qu.:32.00 3rd Qu.:0.1924 3rd Qu.:0.2029 3rd Qu.:0.5301 3rd Qu.: 6.700 Max. :95.00 11 : 2 Max. :89.00 Max. :0.2480 Max. :21.00 Max. :36.00 Max. :0.2315 Max. :0.2945 Max. :3.2821 Max. :15.600 (Other):67 > speed1 <- lmer (Speed.over.25.cm ~ Mass.at.racing + Sex + Month + Sex:Month + (1|Egg.ID), lizards1) > summary(speed1) Linear mixed model fit by REML ['lmerMod'] Formula: Speed.over.25.cm ~ Mass.at.racing + Sex + Month + Sex:Month + (1 | Egg.ID) Data: lizards1 REML criterion at convergence: 107.8 Scaled residuals: Min 1Q Median 3Q Max -1.7356 -0.4326 -0.1597 0.1556 3.9141 Random effects: Groups Name Variance Std.Dev. Egg.ID (Intercept) 0.09044 0.3007 Residual 0.15017 0.3875 Number of obs: 79, groups: Egg.ID, 46 Fixed effects: Estimate Std. Error t value (Intercept) -0.33469 0.34974 -0.957 Mass.at.racing 3.72965 1.64414 2.268 Sex[T.M] -0.02216 0.20889 -0.106 Month[T.July] 0.23496 0.22285 1.054 Sex[T.M]:Month[T.July] 0.17631 0.27269 0.647 Correlation of Fixed Effects: (Intr) Mss.t. Sx[T.M] M[T.J] Mass.t.rcng -0.873 Sex[T.M] -0.433 0.039 Mnth[T.Jly] -0.526 0.174 0.632 S[T.M]:M[T. 0.428 -0.141 -0.770 -0.817 > ###### MODEL SELECTION IN AICcmodavg package (tedious!) > speed1a <- lmer (Speed.over.25.cm ~ Mass.at.racing + Sex + Month + Sex:Month + (1|Egg.ID), lizards1, REML=FALSE) > speed2 <- lmer (Speed.over.25.cm ~ Mass.at.racing + Sex + Month + (1|Egg.ID), lizards1, REML=FALSE) > ## EXCLUDE INTERACTION TERM > anova (speed1a, speed2, test = "LRT") Data: lizards1 Models: speed2: Speed.over.25.cm ~ Mass.at.racing + Sex + Month + (1 | Egg.ID) speed1a: Speed.over.25.cm ~ Mass.at.racing + Sex + Month + Sex:Month + speed1a: (1 | Egg.ID) Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq) speed2 6 113.24 127.46 -50.621 101.24 speed1a 7 114.86 131.44 -50.429 100.86 0.3839 1 0.5355 > ### INTERACTION term is insignificant and could be omitted > speed3 <- lmer (Speed.over.25.cm ~ Mass.at.racing + Month +(1|Egg.ID), lizards1, REML=FALSE) > speed4 <- lmer (Speed.over.25.cm ~ Mass.at.racing + Sex + (1|Egg.ID), lizards1, REML=FALSE) > speed5 <- lmer (Speed.over.25.cm ~ Sex + Month + (1|Egg.ID), lizards1, REML=FALSE) > speed6 <- lmer (Speed.over.25.cm ~ Mass.at.racing + (1|Egg.ID), lizards1, REML=FALSE) > speed7 <- lmer (Speed.over.25.cm ~ Sex + (1|Egg.ID), lizards1, REML=FALSE) > speed8 <- lmer (Speed.over.25.cm ~ Month + (1|Egg.ID), lizards1, REML=FALSE) > speed9 <- lmer (Speed.over.25.cm ~ 1 + (1|Egg.ID), lizards1, REML=FALSE) > Cand.models1 <- list(speed2, speed3, speed4, speed5, speed6, speed7, speed8, speed9) > Modnames <- c("speed2", "speed3", "speed4", "speed5", "speed6", "speed7", "speed8", "speed9") > print(aictab(cand.set = Cand.models1, modnames = Modnames, second.ord = TRUE), digits = 4) Model selection based on AICc: K AICc Delta_AICc AICcWt Cum.Wt LL speed3 5 112.2897 0.0000 0.6353 0.6353 -50.7339 speed2 6 114.4081 2.1184 0.2203 0.8555 -50.6207 speed8 4 117.0127 4.7230 0.0599 0.9154 -54.2361 speed6 4 118.3764 6.0867 0.0303 0.9457 -54.9179 speed5 5 118.5304 6.2407 0.0280 0.9737 -53.8542 speed4 5 120.3290 8.0393 0.0114 0.9851 -54.7536 speed9 3 120.5483 8.2586 0.0102 0.9954 -57.1142 speed7 4 122.1307 9.8410 0.0046 1.0000 -56.7951 > ###### MODEL SELECTION IN MuMIn package (easy!) > options(na.action = na.fail) > res5 <- dredge(speed1, rank="AICc", REML=F) > est5 <- model.avg(res5, subset=delta<7, revised.var=T) > summary(est5) Call: model.avg(object = res5, subset = delta < 7, revised.var = T) Component model call: lmer(formula = Speed.over.25.cm ~ <6 unique rhs>, data = lizards1) Component models: df logLik AICc delta weight 12 5 -52.72 112.65 0.00 0.60 123 6 -53.69 114.89 2.23 0.20 1234 7 -53.88 117.03 4.38 0.07 2 4 -57.14 117.07 4.42 0.07 1 4 -55.33 118.42 5.77 0.03 23 5 -57.81 118.65 6.00 0.03 Term codes: Mass.at.racing Month Sex Month:Sex 1 2 3 4 Model-averaged coefficients: (full average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) -0.33615 0.36574 0.36968 0.909 0.3632 Mass.at.racing 3.63039 1.92112 1.94084 1.871 0.0614 . Month[T.July] 0.32239 0.14573 0.14767 2.183 0.0290 * Sex[T.M] 0.01768 0.09095 0.09223 0.192 0.8480 Month[T.July]:Sex[T.M] 0.01192 0.08357 0.08460 0.141 0.8880 (conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) -0.33615 0.36574 0.36968 0.909 0.3632 Mass.at.racing 4.01848 1.58924 1.61555 2.487 0.0129 * Month[T.July] 0.33368 0.13497 0.13713 2.433 0.0150 * Sex[T.M] 0.05982 0.15959 0.16206 0.369 0.7120 Month[T.July]:Sex[T.M] 0.17631 0.27269 0.27735 0.636 0.5250 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Relative variable importance: Month Mass.at.racing Sex Month:Sex Importance: 0.97 0.90 0.30 0.07 N containing models: 5 4 3 1 > ############### LIKELIHOOD RATIO TESTS library(foreign, pos=23) library(survival, pos=23) library(MASS, pos=23) library(epiDisplay, pos=23) > ### SEX EFFECT > lrtest (speed2, speed3) Likelihood ratio test Model 1: Speed.over.25.cm ~ Mass.at.racing + Sex + Month + (1 | Egg.ID) Model 2: Speed.over.25.cm ~ Mass.at.racing + Month + (1 | Egg.ID) #Df LogLik Df Chisq Pr(>Chisq) 1 6 -50.621 2 5 -50.734 -1 0.2263 0.6343 > #### TUKEY POST HOC COMPARISON library(emmeans, pos=27) > lsmeans(speed2, pairwise ~ Sex*Month, adjust="tukey") $lsmeans Sex Month lsmean SE df lower.CL upper.CL F April 0.276 0.1288 79.0 0.0193 0.532 M April 0.327 0.1061 81.5 0.1165 0.538 F July 0.597 0.1164 79.1 0.3651 0.828 M July 0.649 0.0917 80.4 0.4662 0.831 Degrees-of-freedom method: kenward-roger Confidence level used: 0.95 $contrasts contrast estimate SE df t.ratio p.value F,April - M,April -0.0519 0.130 79.6 -0.399 0.9783 F,April - F,July -0.3211 0.125 80.8 -2.559 0.0585 F,April - M,July -0.3730 0.181 79.6 -2.060 0.1754 M,April - F,July -0.2693 0.180 80.8 -1.496 0.4445 M,April - M,July -0.3211 0.125 80.8 -2.559 0.0585 F,July - M,July -0.0519 0.130 79.6 -0.399 0.9783 Degrees-of-freedom method: kenward-roger P value adjustment: tukey method for comparing a family of 4 estimates ########################################################### ######################SCRIPT############################## ########################################### ########## MIXED EFFECT LINEAR MODELS summary (lizards1) lizards1 <- within(lizards1, { Egg.ID <- as.factor(Egg.ID) }) library(lme4, pos=28) lizards1 <- within(lizards1, { weeks.at.racing <- as.factor(weeks.at.racing) }) speed1 <- lmer (Speed.over.25.cm ~ Mass.at.racing + Sex + Month + Sex:Month + (1|Egg.ID), lizards1) summary(speed1) library(AICcmodavg, pos=30) ###### MODEL SELECTION IN AICcmodavg package (tedious!) speed1a <- lmer (Speed.over.25.cm ~ Mass.at.racing + Sex + Month + Sex:Month + (1|Egg.ID), lizards1, REML=FALSE) speed2 <- lmer (Speed.over.25.cm ~ Mass.at.racing + Sex + Month + (1|Egg.ID), lizards1, REML=FALSE) ## EXCLUDE INTERACTION TERM anova (speed1a, speed2, test = "LRT") ### INTERACTION term is insignificant and could be omitted speed3 <- lmer (Speed.over.25.cm ~ Mass.at.racing + Month +(1|Egg.ID), lizards1, REML=FALSE) speed4 <- lmer (Speed.over.25.cm ~ Mass.at.racing + Sex + (1|Egg.ID), lizards1, REML=FALSE) speed5 <- lmer (Speed.over.25.cm ~ Sex + Month + (1|Egg.ID), lizards1, REML=FALSE) speed6 <- lmer (Speed.over.25.cm ~ Mass.at.racing + (1|Egg.ID), lizards1, REML=FALSE) speed7 <- lmer (Speed.over.25.cm ~ Sex + (1|Egg.ID), lizards1, REML=FALSE) speed8 <- lmer (Speed.over.25.cm ~ Month + (1|Egg.ID), lizards1, REML=FALSE) speed9 <- lmer (Speed.over.25.cm ~ 1 + (1|Egg.ID), lizards1, REML=FALSE) Cand.models1 <- list(speed2, speed3, speed4, speed5, speed6, speed7, speed8, speed9) Modnames <- c("speed2", "speed3", "speed4", "speed5", "speed6", "speed7", "speed8", "speed9") print(aictab(cand.set = Cand.models1, modnames = Modnames, second.ord = TRUE), digits = 4) ###### MODEL SELECTION IN MuMIn package (easy!) options(na.action = na.fail) res5 <- dredge(speed1, rank="AICc", REML=F) est5 <- model.avg(res5, subset=delta<7, revised.var=T) summary(est5) ############### LIKELIHOOD RATIO TESTS library(foreign, pos=23) library(survival, pos=23) library(MASS, pos=23) library(epiDisplay, pos=23) ### SEX EFFECT lrtest (speed2, speed3) ### Month effect lrtest (speed2, speed4) ### Mass effect lrtest (speed2, speed5) #### TUKEY POST HOC COMPARISON library(emmeans, pos=27) lsmeans(speed2, pairwise ~ Sex*Month, adjust="tukey")