# analysis of Desmognathus escape performance library(car) library(lme4) library(performance) library(emmeans) escape.dat <- read.csv("peerj-113480-QMescape.csv") names(escape.dat) # Species: species names except for the northern semi-aquatic form, which is listed as "bbnc" for 'black-bellied, northern clade' # Ecomorph: Aquatic or Semi-aquatic # Clade: Northern or Southern table(escape.dat$Species,escape.dat$Ecomorph,escape.dat$Clade) # Locality is "VA" for northern and "GA" for southern # ID is an arbitrary individual ID assigned at the time of capture # HL is head length # TrL is trunk length (groin to axilla) # TaL is tail length # VND and VNW are net velocities in the terrestrial (D) and aquatic (W) arenas, respectively # add body lengths escape.dat$BodyL <- escape.dat$HL+escape.dat$TrL escape.dat$TotL <- escape.dat$BodyL+escape.dat$TaL # compute relative velocity escape.dat$VNDr <- escape.dat$VND/escape.dat$BodyL escape.dat$VNWr <- escape.dat$VNW/escape.dat$BodyL # Basic test for trade-off (negative correlation) cor.test(escape.dat$VNDr,escape.dat$VNWr) # Pearson's product-moment correlation # data: escape.dat$VNDr and escape.dat$VNWr # t = -4.507, df = 42, p-value = 5.185e-05 # alternative hypothesis: true correlation is not equal to 0 # 95 percent confidence interval: # -0.7420477 -0.3300024 # sample estimates: # cor # -0.5709465 # breakdown by north vs. south cor.test(escape.dat$VNDr[escape.dat$Clade=="Northern"],escape.dat$VNWr[escape.dat$Clade=="Northern"]) # r = -0.5133669; p = 0.0103 # 95 percent confidence interval: # -0.7594821 -0.1386920 cor.test(escape.dat$VNDr[escape.dat$Clade=="Southern"],escape.dat$VNWr[escape.dat$Clade=="Southern"], exact=TRUE) # r = 0.6200049; p = 0.003544 # 95 percent confidence interval: # -0.8337687 -0.2445915 ### Scatterplot velocity in body lengths per sec quartz(width=8,height=5) par(mfrow=c(1,2), mar=c(5,5,3,2)) plot(escape.dat$VNDr,escape.dat$VNWr, xlab="Running velocity (BL/s)", ylab="Swimming velocity (BL/s)", pch="") points(escape.dat$VNDr[escape.dat$Species=="Dbbnc"],escape.dat$VNWr[escape.dat$Species=="Dbbnc"],pch=1) points(escape.dat$VNDr[escape.dat$Species=="Dmarmoratus"],escape.dat$VNWr[escape.dat$Species=="Dmarmoratus"],pch=16) mtext("(A) Northern ecomorphs", at=0.1, line=1) legend("bottomright", c("Semi-aquatic", "Aquatic"), pch=c(1,16), inset=c(0,1), xpd=TRUE, horiz=FALSE, bty="n" ) plot(escape.dat$VNDr,escape.dat$VNWr, xlab="Running velocity (BL/s)", ylab="Swimming velocity (BL/s)", pch="") points(escape.dat$VNDr[escape.dat$Species=="Damphileucus"],escape.dat$VNWr[escape.dat$Species=="Damphileucus"],pch=1) points(escape.dat$VNDr[escape.dat$Species=="Daureatus"],escape.dat$VNWr[escape.dat$Species=="Daureatus"],pch=16) mtext("(B) Southern ecomorphs", at=0.1, line=1) legend("bottomright", c("Semi-aquatic", "Aquatic"), pch=c(1,16), inset=c(0,1), xpd=TRUE, horiz=FALSE, bty="n" ) ###### New Analysis ## Single model: mixed effects, with species instead of ecomorph*clade (same thing, just simpler to interpret and report) # Test multicolinearity of main effects (interaction terms tend to increase VIFs) # Include species*bodyLength # Use model selection to simplify # Report coef of determination (some disagreement on the web regarding mixed models) # Extract marginal mean velocity for each species:habitat combination # Extract specific 'planned' comparisons: # aquatic vs. semi, northern clade, dry # aquatic vs. semi, northern clade, wet # aquatic vs. semi, southern clade, dry # aquatic vs. semi, southern clade, wet ###### Secondary Analysis ## Random effect above causes singularity warning because individual observations are negatively rather than positively correlated. This is probably fine, but just to be sure estimates are not distorted, perform separate analyses for wet vs dry tests. Consider adjusting significance cutoffs to maintain global alpha=0.05 # first combine aquatic and terrestrial Velocity <- c(escape.dat$VND,escape.dat$VNW) relVelocity <- Velocity/escape.dat$BodyL Individual <- as.factor(rep(escape.dat$ID,2)) Habitat <- rep(c("D","W"),each=nrow(escape.dat)) Covariates <- rbind(escape.dat,escape.dat) # Minimal models for diagnostics V.fit <- lmer(Velocity~Habitat+Species+BodyL+HL+TaL + (1|Individual), data=Covariates) # singularity warning because random effect (individual) variance is approximately zero ... as it should be if faster runners do NOT tend to be faster swimmers vif(V.fit) # morphometrics have multicolinearity Anova(V.fit) # drop HL and TaL V.fit <- lmer(Velocity~Habitat+Species+BodyL + (1|Individual), data=Covariates) vif(V.fit) # all good # residual diagnostics plot(V.fit) qqnorm(residuals(V.fit)) shapiro.test(residuals(V.fit)) # W = 0.98561, p-value = 0.4418 # All good: no need for transformations # Maximal model V.fit1 <- lmer(Velocity~Habitat*Species*BodyL + (1|Individual), data=Covariates) Anova(V.fit1) # Analysis of Deviance Table (Type II Wald chisquare tests) # Response: Velocity # Chisq Df Pr(>Chisq) # Habitat 18.9477 1 1.344e-05 *** # Species 3.7445 3 0.29041 # BodyL 28.2361 1 1.074e-07 *** # Habitat:Species 72.1850 3 1.453e-15 *** # Habitat:BodyL 0.4042 1 0.52494 # Species:BodyL 0.5772 3 0.90163 # Habitat:Species:BodyL 7.1427 3 0.06748 . # --- AIC(V.fit1) # 342.5413 AIC(V.fit) # 393.8661 AIC(lmer(Velocity~Habitat+Species+BodyL +Habitat:Species+Habitat:BodyL+Species:BodyL + (1|Individual), data=Covariates)) # 348.3492 # full model is much better r2_nakagawa(V.fit1) # Conditional R2: NA (fixed and random effects, NA because random effect ~ 0) # Marginal R2: 0.621 (only fixed effects) difftests <- emmeans(V.fit1, list(pairwise ~ Habitat:Species), adjust = "none") # marginal means difftests[1] # Habitat Species emmean SE df lower.CL upper.CL # D Damphileucus 5.76 0.546 72 4.67 6.84 # W Damphileucus 3.72 0.546 72 2.63 4.80 # D Daureatus 2.90 0.507 72 1.89 3.91 # W Daureatus 6.70 0.507 72 5.69 7.71 # D Dbbnc 5.10 0.614 72 3.88 6.33 # W Dbbnc 4.28 0.614 72 3.06 5.51 # D Dmarmoratus 3.25 0.475 72 2.30 4.20 # W Dmarmoratus 7.58 0.475 72 6.63 8.52 # Degrees-of-freedom method: kenward-roger # Confidence level used: 0.95 # pairs of interest difftests[2][[1]][c(2,9,24,27),] # 1 estimate SE df t.ratio p.value # D Damphileucus - D Daureatus 2.86 0.745 72 3.841 0.0003 # W Damphileucus - W Daureatus -2.98 0.745 72 -4.003 0.0002 # D Dbbnc - D Dmarmoratus 1.85 0.776 72 2.383 0.0198 # W Dbbnc - W Dmarmoratus -3.29 0.776 72 -4.240 0.0001 # note: these estimates are essentially unchanged if the bodyL interactions are dropped # # Secondary analysis ### separate by habitat dry.fit1 <- lm(VND~Species*BodyL, data=escape.dat) summary(dry.fit1)$adj.r.squared # 0.641615 drop1(dry.fit1) # Single term deletions # Model: # VND ~ Species * BodyL # Df Sum of Sq RSS AIC # 54.161 25.142 # Species:BodyL 3 5.6886 59.849 23.536 dry.fit <- lm(VND~Species+BodyL, data=escape.dat) Anova(dry.fit) # Response: VND # Sum Sq Df F value Pr(>F) # Species 61.637 3 13.388 3.738e-06 *** # BodyL 43.245 1 28.180 4.712e-06 *** # Residuals 59.849 39 summary(dry.fit)$adj.r.squared # 0.6344369 emmeans(dry.fit, list(pairwise~Species), adjust="none") # $`emmeans of Species` # Species emmean SE df lower.CL upper.CL # Damphileucus 5.92 0.408 39 5.09 6.74 # Daureatus 3.02 0.393 39 2.23 3.82 # Dbbnc 4.87 0.384 39 4.10 5.65 # Dmarmoratus 3.17 0.364 39 2.43 3.90 # Confidence level used: 0.95 # $`pairwise differences of Species` # 1 estimate SE df t.ratio p.value # Damphileucus - Daureatus 2.897 0.572 39 5.063 <.0001 # Damphileucus - Dbbnc 1.043 0.589 39 1.771 0.0844 # Damphileucus - Dmarmoratus 2.753 0.533 39 5.168 <.0001 # Daureatus - Dbbnc -1.854 0.542 39 -3.418 0.0015 # Daureatus - Dmarmoratus -0.144 0.539 39 -0.268 0.7905 # Dbbnc - Dmarmoratus 1.710 0.547 39 3.128 0.0033 emmeans(dry.fit1, list(pairwise~Species), adjust="none") # $`emmeans of Species` # Species emmean SE df lower.CL upper.CL # Damphileucus 5.76 0.428 36 4.89 6.63 # Daureatus 2.90 0.398 36 2.09 3.70 # Dbbnc 5.10 0.482 36 4.12 6.08 # Dmarmoratus 3.25 0.373 36 2.50 4.01 # Confidence level used: 0.95 # $`pairwise differences of Species` # 1 estimate SE df t.ratio p.value # Damphileucus - Daureatus 2.860 0.584 36 4.895 <.0001 # Damphileucus - Dbbnc 0.657 0.645 36 1.019 0.3152 # Damphileucus - Dmarmoratus 2.507 0.567 36 4.417 0.0001 # Daureatus - Dbbnc -2.204 0.625 36 -3.526 0.0012 # Daureatus - Dmarmoratus -0.354 0.545 36 -0.649 0.5205 # Dbbnc - Dmarmoratus 1.850 0.609 36 3.036 0.0044 wet.fit1 <- lm(VNW~Species*BodyL, data=escape.dat) summary(wet.fit1)$adj.r.squared # 0.5119936 drop1(wet.fit1) # Single term deletions # Model: # VNW ~ Species * BodyL # Df Sum of Sq RSS AIC # 121.77 60.789 # Species:BodyL 3 13.175 134.94 59.309 wet.fit <- lm(VNW~Species+BodyL, data=escape.dat) Anova(wet.fit) # Response: VNW # Sum Sq Df F value Pr(>F) # Species 123.893 3 11.9356 1.098e-05 *** # BodyL 26.736 1 7.7271 0.008329 ** # Residuals 134.942 39 # --- summary(wet.fit)$adj.r.squared # 0.5007941 emmeans(wet.fit, list(pairwise~Species), adjust="none") # $`emmeans of Species` # Species emmean SE df lower.CL upper.CL # Damphileucus 3.44 0.613 39 2.20 4.68 # Daureatus 6.58 0.590 39 5.39 7.78 # Dbbnc 4.67 0.577 39 3.50 5.83 # Dmarmoratus 7.79 0.546 39 6.68 8.89 # Confidence level used: 0.95 # $`pairwise differences of Species` # 1 estimate SE df t.ratio p.value # Damphileucus - Daureatus -3.14 0.859 39 -3.658 0.0008 # Damphileucus - Dbbnc -1.22 0.884 39 -1.385 0.1740 # Damphileucus - Dmarmoratus -4.35 0.800 39 -5.438 <.0001 # Daureatus - Dbbnc 1.92 0.814 39 2.354 0.0237 # Daureatus - Dmarmoratus -1.21 0.809 39 -1.492 0.1437 # Dbbnc - Dmarmoratus -3.12 0.821 39 -3.807 0.0005 emmeans(wet.fit1, list(pairwise~Species), adjust="none") # $`emmeans of Species` # Species emmean SE df lower.CL upper.CL # Damphileucus 3.72 0.642 36 2.41 5.02 # Daureatus 6.70 0.596 36 5.49 7.91 # Dbbnc 4.28 0.723 36 2.82 5.75 # Dmarmoratus 7.58 0.559 36 6.44 8.71 # Confidence level used: 0.95 # $`pairwise differences of Species` # 1 estimate SE df t.ratio p.value # Damphileucus - Daureatus -2.981 0.876 36 -3.403 0.0017 # Damphileucus - Dbbnc -0.568 0.967 36 -0.587 0.5608 # Damphileucus - Dmarmoratus -3.860 0.851 36 -4.536 0.0001 # Daureatus - Dbbnc 2.414 0.937 36 2.576 0.0143 # Daureatus - Dmarmoratus -0.879 0.817 36 -1.075 0.2894 # Dbbnc - Dmarmoratus -3.293 0.914 36 -3.604 0.0009 ###################################################################################### # To cite package ‘emmeans’ in publications use: # Lenth R (2025). _emmeans: Estimated Marginal Means, aka Least-Squares Means_. R package version 1.11.0, . ######################################################################### with relative velocity relVelocity <- Velocity/Covariates$BodyL shapiro.test(relVelocity) rel.fit1 <- lmer(relVelocity~Habitat*Species+BodyL + (1|Individual), data=Covariates) # singularity warning because random effect (individual) variance is approximately zero ... as it should be if faster runners do NOT tend to be faster swimmers # with head and tail length to see of body proportions have an effect rel.fit2 <- lmer(relVelocity~Habitat*Species+BodyL+HL+TaL + (1|Individual), data=Covariates) AIC(rel.fit1) # 59.53001 AIC(rel.fit2) # 68.46407 # simpler model is superior by > 2 AIC rel.fit0 <- lmer(relVelocity~Habitat*Species + (1|Individual), data=Covariates) AIC(rel.fit0) # 52.16361 # residual plots plot(rel.fit0) hist(residuals(rel.fit0)) qqnorm(residuals(rel.fit0)) boxplot(residuals(rel.fit0)~Covariates$Species:Habitat);abline(h=0,lty=2) shapiro.test(residuals(rel.fit0)) Anova(rel.fit0) # Analysis of Deviance Table (Type II Wald chisquare tests) # Response: relVelocity # Chisq Df Pr(>Chisq) # Habitat 21.4498 1 3.632e-06 *** # Species 3.5124 3 0.3192 # Habitat:Species 64.2895 3 7.117e-14 *** # --- rdifftests <- emmeans(rel.fit0, list(pairwise ~ Habitat:Species), adjust = "none") # marginal means rdifftests[1] # Habitat Species emmean SE df lower.CL upper.CL # D Damphileucus 0.994 0.0830 80 0.829 1.159 # W Damphileucus 0.610 0.0830 80 0.445 0.775 # D Daureatus 0.509 0.0830 80 0.344 0.674 # W Daureatus 1.130 0.0830 80 0.964 1.295 # D Dbbnc 0.803 0.0758 80 0.652 0.954 # W Dbbnc 0.821 0.0758 80 0.671 0.972 # D Dmarmoratus 0.561 0.0758 80 0.410 0.711 # W Dmarmoratus 1.296 0.0758 80 1.145 1.446 # Degrees-of-freedom method: kenward-roger # Confidence level used: 0.95 # pairs of interest rdifftests[2][[1]][c(2,9,24,27),] # 1 estimate SE df t.ratio p.value # D Damphileucus - D Daureatus 0.485 0.117 80 4.132 0.0001 # W Damphileucus - W Daureatus -0.520 0.117 80 -4.429 <.0001 # D Dbbnc - D Dmarmoratus 0.242 0.107 80 2.260 0.0265 # W Dbbnc - W Dmarmoratus -0.474 0.107 80 -4.426 <.0001 ######################################################################### ### some other alternative analyses ... all agree regarding conclusions # ######################################################################### ### breakdown Habitat x species interaction, accounting for body length dry.fitR <- lm(VNDr~BodyL+Ecomorph*Clade,data=escape.dat) summary(dry.fitR) Anova(dry.fitR) # Response: rVND # Sum Sq Df F value Pr(>F) Coef # BodyL 0.04871 1 0.9940 0.3249 0.03185 # Ecomorph 1.45023 1 29.5946 3.099e-06 *** 0.27912 (BB > SN by 0.28 body lengths per second) # Clade 0.02414 1 0.4926 0.4869 -0.03483 # Interaction 0.07290 1 1.4876 0.2299 0.18042 # Residuals 1.91113 39 # --- wet.fitR <- lm(VNWr~BodyL+Ecomorph*Clade,data=escape.dat) summary(wet.fitR) Anova(wet.fitR) # Response: rVNW # Sum Sq Df F value Pr(>F) Coef # BodyL 0.0245 1 0.2714 0.60535 -0.022608 # Ecomorph 2.7252 1 30.1343 2.648e-06 ***-0.500498 (BB < SN by 0.5 body lengths per second) # Clade 0.3350 1 3.7048 0.06157 . -0.177946 (N > S) # Interaction 0.0000 1 0.0000 0.99469 -0.001345 (BB*S) # Residuals 3.5270 39 # --- # same as above but Ecomorph*Clade captured as "Species" dry.fitRS <- lm(VNDr~BodyL+Species,data=escape.dat) wet.fitRS <- lm(VNWr~BodyL+Species,data=escape.dat) TukeyHSD(aov(dry.fitRS), which="Species") # Tukey multiple comparisons of means # 95% family-wise confidence level # Fit: aov(formula = dry.fitRS) # $Species # diff lwr upr p adj # Daureatus-Damphileucus -0.45740178 -0.72305043 -0.1917531 0.0002331 # Dbbnc-Damphileucus -0.14175685 -0.39609607 0.1125824 0.4498281 # Dmarmoratus-Damphileucus -0.42398368 -0.67832289 -0.1696445 0.0003655 # Dbbnc-Daureatus 0.31564493 0.06130572 0.5699841 0.0098770 # Dmarmoratus-Daureatus 0.03341811 -0.22092111 0.2877573 0.9847381 # Dmarmoratus-Dbbnc -0.28222682 -0.52472975 -0.0397239 0.0170219 TukeyHSD(aov(wet.fitRS), which="Species") # diff lwr upr p adj # Daureatus-Damphileucus 0.5096192 0.1487355 0.87050283 0.0027662 # Dbbnc-Damphileucus 0.1932340 -0.1522858 0.53875383 0.4468442 # Dmarmoratus-Damphileucus 0.6824323 0.3369125 1.02795214 0.0000280 # Dbbnc-Daureatus -0.3163851 -0.6619049 0.02913466 0.0830233 # Dmarmoratus-Daureatus 0.1728132 -0.1727066 0.51833297 0.5424172 # Dmarmoratus-Dbbnc 0.4891983 0.1597581 0.81863851 0.0015725 ### simple: take individual differences DWdiff <- escape.dat$VND-escape.dat$VNW DWrdiff <- escape.dat$VNDr-escape.dat$VNWr diff.fit <- aov(DWdiff~Species+BodyL, data=escape.dat) Anova(diff.fit) # Anova Table (Type II tests) # Response: DWdiff # Sum Sq Df F value Pr(>F) # Species 352.76 3 18.2456 1.478e-07 *** # BodyL 1.98 1 0.3065 0.583 # Residuals 251.34 39 # --- diff.tukey <- TukeyHSD(diff.fit,which="Species") # diff lwr upr p adj # Daureatus-Damphileucus -6.2005632 -9.2470203 -3.1541061 0.0000168 # Dbbnc-Damphileucus -2.5578504 -5.4746110 0.3589103 0.1035189 # Dmarmoratus-Damphileucus -7.1566735 -10.0734341 -4.2399129 0.0000005 # Dbbnc-Daureatus 3.6427128 0.7259522 6.5594734 0.0093347 # Dmarmoratus-Daureatus -0.9561104 -3.8728710 1.9606502 0.8153644 # Dmarmoratus-Dbbnc -4.5988232 -7.3798453 -1.8178011 0.0004076 rdiff.fit <- aov(DWrdiff~Species+BodyL, data=escape.dat) Anova(rdiff.fit) # Response: DWrdiff # Sum Sq Df F value Pr(>F) # Species 8.8624 3 17.1460 2.937e-07 *** # BodyL 0.1424 1 0.8265 0.3689 # Residuals 6.7194 39 # --- rdiff.tukey <- TukeyHSD(rdiff.fit,which="Species") # diff lwr upr p adj # Daureatus-Damphileucus -1.0048578 -1.5029718 -0.50674393 0.0000195 # Dbbnc-Damphileucus -0.4028347 -0.8797424 0.07407309 0.1236286 # Dmarmoratus-Damphileucus -1.1192791 -1.5961869 -0.64237135 0.0000012 # Dbbnc-Daureatus 0.6020232 0.1251154 1.07893093 0.0084694 # Dmarmoratus-Daureatus -0.1144213 -0.5913290 0.36248649 0.9170692 # Dmarmoratus-Dbbnc -0.7164444 -1.1711581 -0.26173073 0.0007661 quartz() par(mfrow=c(1,2)) plot(diff.tukey) plot(rdiff.tukey) # easier to interpret: par(mfrow=c(1,2),mar=c(8,7,2,3)) boxplot(DWdiff~Species,data=escape.dat, horizontal=TRUE, las=1, ylab="", ylim=c(-10,10), xlab="", main="(A) Difference in velocity (cm/s)") abline(v=0,lty=3) boxplot(DWrdiff~Species,data=escape.dat, horizontal=TRUE, las=1, ylab="", ylim=c(-1.6,1.6), xlab="", main="(B) Difference in relative velocity (BL/s)") abline(v=0,lty=3) quartz(); par(mar=c(6,8,3,2)) boxplot(DWdiff~Species,data=escape.dat, horizontal=TRUE, las=1, ylab="", ylim=c(-10,10), xlab="", main="Difference in velocity (cm/s)") abline(v=0,lty=3) mtext("Faster in water", side=1, line=2.5, at=-5) mtext("Faster on land", side=1, line=2.5, at=5) ############################################################################