Winter 2020 Flight Trials: Mass Modeling

Flight Trials Winter 2020 Dataset was conducted from 2/17/2020 - 3/10/2020. Soapberry bugs were flight tested twice for multiple hours in the flight mill and observed from 8 AM to (5-8 PM) each day. Used multivariate (lm) and mixed effect modeling (lmer) to analyze the mass results.

All Data

Cleaning Data

lmer() host, sex, sym_dist | ID and trial_type

Morphology

Reading the data

rm(list=ls())
output_col = FALSE # Recommend changing this to TRUE if working in Base R or RStudio, and FALSE if generating an html
source("src/clean_flight_data.R") # Script that loads and cleans up the data
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
source("src/regression_output.R") # A script that cleans up regression outputs and prints in color or black and white
source("src/center_flight_data.R")
source("src/get_warnings.R")

data <- read_flight_data("data/all_flight_data-Winter2020.csv")
data_all <- data[[1]]
data_tested <- data[[2]]

Testing Models and Covariates

####Check for convergence without splitting by trial type
test_model<-lmer(mass~host_c*sex_c*sym_dist+ (1|trial_type) + (1|ID), data=data_tested, REML=FALSE) 
tidy_regression(test_model, is_color=output_col)
## lmer mass ~ host_c * sex_c * sym_dist + (1 | trial_type) + (1 | ID) data_tested FALSE 
## AIC:  -3504.559 -3504.559 
## (Intercept)              coeff:  0.0595158   tval:  25.87272
## host_c                   coeff:  -0.0046768  tval:  -2.35777
## sex_c                    coeff:  0.0240503   tval:  12.1226
## sym_dist                 coeff:  -0.0063641  tval:  -1.459875
## host_c:sex_c             coeff:  -0.0032198  tval:  -1.62324
## host_c:sym_dist          coeff:  0.0056225   tval:  1.289781
## sex_c:sym_dist           coeff:  -0.0107798  tval:  -2.472509
## host_c:sex_c:sym_dist    coeff:  0.0090185   tval:  2.068683
getME(test_model, "lower")
## [1] 0 0

Experimental Set-Up Effects

## lmer mass_c ~ chamber + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -3158.876 -3158.876 
## (Intercept)  coeff:  -6.5e-05    tval:  -0.02715385
## chamberA-2   coeff:  -2.78e-05   tval:  -0.01242178
## chamberA-3   coeff:  0.0014132   tval:  0.6492935
## chamberA-4   coeff:  -0.0011209  tval:  -0.5387897
## chamberB-1   coeff:  0.0007711   tval:  0.3244294
## chamberB-2   coeff:  -0.0005543  tval:  -0.2597971
## chamberB-3   coeff:  0.0008228   tval:  0.3685056
## chamberB-4   coeff:  0.0013812   tval:  0.6343744
## lmer mass_c ~ days_from_start_c + (1 | ID) data_tested FALSE 
## AIC:  -3172.693 -3172.693 
## (Intercept)          coeff:  8.98e-05    tval:  0.06844824
## days_from_start_c    coeff:  0.0002329   tval:  3.941263
## lmer mass_c ~ min_from_IncStart_c + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -3167.725 -3167.725 
## (Intercept)          coeff:  0.0001751   tval:  0.09777721
## min_from_IncStart_c  coeff:  2e-07   tval:  0.04538841

Biological Effects

## lmer mass_c ~ total_eggs + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -933.7505 -933.7505 
## (Intercept)  coeff:  0.0168041   tval:  5.759655
## total_eggs   coeff:  0.0001144   tval:  5.728754
## lmer mass_c ~ eggs_b + (1 | ID) data_tested FALSE 
## AIC:  -3377.796 -3377.796 
## (Intercept)  coeff:  -0.0070548  tval:  -7.199589
## eggs_b       coeff:  0.030596    tval:  18.52713

Morphology Effects

## lmer mass_c ~ beak_c + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -3476.712 -3476.712 
## (Intercept)  coeff:  0.0001451   tval:  0.1036325
## beak_c       coeff:  0.0181976   tval:  22.65519
## lmer mass_c ~ thorax_c + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -3453.797 -3453.797 
## (Intercept)  coeff:  8.78e-05    tval:  0.06098014
## thorax_c     coeff:  0.0596814   tval:  21.37375
## lmer mass_c ~ body_c + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -3431.141 -3431.141 
## (Intercept)  coeff:  0.0001935   tval:  0.1343313
## body_c       coeff:  0.0165333   tval:  20.04531
## lmer mass_c ~ wing_c + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -3341.774 -3341.774 
## (Intercept)  coeff:  0.0001803   tval:  0.1167289
## wing_c       coeff:  0.0185423   tval:  15.12271

All Data | lmer() A=host, B=sex, C=sym_dist, X=ID

data<-data.frame(R=data_tested$mass_c, 
                 A=data_tested$host_c, 
                 B=data_tested$sex_c, 
                 C=data_tested$sym_dist, 
                 X=data_tested$ID)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R"))
##        [,1]      [,2]      [,3]      [,4]      [,5]      [,6]       [,7]      
## AICs   -3494.901 -3493.398 -3493.123 -3492.966 -3492.913 -3491.608  -3491.415 
## models 13        16        8         4         15        10         17        
## probs  0.2882015 0.135925  0.1184743 0.1094847 0.1066492 0.05552537 0.05043566
## 
## m13  R ~ B * C + A + (1 | X)
## m16  R ~ A * C + B * C + (1 | X)
## m8   R ~ A * B + (1 | X)
## m4   R ~ A + B + (1 | X)
## m15  R ~ A * B + B * C + (1 | X)
## m10  R ~ B * C + (1 | X)
## m17  R ~ A * B + A * C + B * C + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings))
## Number of models that failed to converge:  0
anova(m13, m15, test="Chisq") # Adding A*B does not improve fit
## Data: data
## Models:
## m13: R ~ B * C + A + (1 | X)
## m15: R ~ A * B + B * C + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m13  7 -3494.9 -3464.0 1754.5  -3508.9                         
## m15  8 -3492.9 -3457.6 1754.5  -3508.9 0.0118      1     0.9136
model_alldata <- lmer(mass_c ~ sex_c * sym_dist + host_c + (1|ID), data=data_tested, REML=FALSE)
tidy_regression(model_alldata, is_color=output_col) 
## lmer mass_c ~ sex_c * sym_dist + host_c + (1 | ID) data_tested FALSE 
## AIC:  -3494.901 -3494.901 
## (Intercept)      coeff:  0.0049957   tval:  3.403682
## sex_c            coeff:  0.0212837   tval:  21.06267
## sym_dist         coeff:  -0.0010486  tval:  -0.9496991
## host_c           coeff:  -0.0027384  tval:  -2.308398
## sex_c:sym_dist   coeff:  -0.0021001  tval:  -2.43957
model_alldata <- lmer(mass_c ~ sex_c * sym_dist + host_c + (1|ID) + (1|trial_type), data=data_tested, REML=FALSE)
tidy_regression(model_alldata, is_color=output_col) # adding trial type improves fit
## lmer mass_c ~ sex_c * sym_dist + host_c + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -3505.776 -3505.776 
## (Intercept)      coeff:  0.005125    tval:  2.732317
## sex_c            coeff:  0.0212714   tval:  21.076
## sym_dist         coeff:  -0.0009623  tval:  -0.873093
## host_c           coeff:  -0.0027059  tval:  -2.283953
## sex_c:sym_dist   coeff:  -0.0020552  tval:  -2.392181

  • positive effect of sex, where if Female then weigh more.

  • negative weak effect of host where if from GRT then weigh less.

  • negative weak effect of sym_dist

  • negative weak effect of sex*sym_dist where if F and farther from Homestead, then weigh less

Morphology Modeling

data<-data.frame(R=data_tested$mass_c, 
                 A=data_tested$beak_c, 
                 B=data_tested$thorax_c, 
                 C=data_tested$body_c, 
                 D=data_tested$wing_c,
                 X=data_tested$ID)
data <- data %>%
  filter(!is.na(C))

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF REMLF.R"))
##        [,1]      [,2]      [,3]      [,4]     
## AICs   -3612.036 -3610.19  -3610.124 -3608.516
## models 91        107       106       112      
## probs  0.2922233 0.1161077 0.1123722 0.0502832
## 
## m91  R ~ A * B + A * C + A * D + B * C + (1 | X)
## m107     R ~ A * B + A * C + A * D + B * C + C * D + (1 | X)
## m106     R ~ A * B + A * C + A * D + B * C + B * D + (1 | X)
## m112     R ~ A * B + A * C + A * D + B * C + B * D + C * D + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings))
## Number of models that failed to converge:  0
anova(m91, m107, test="Chisq") # Adding C*D does not improve fit
## Data: data
## Models:
## m91: R ~ A * B + A * C + A * D + B * C + (1 | X)
## m107: R ~ A * B + A * C + A * D + B * C + C * D + (1 | X)
##      Df     AIC     BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m91  11 -3612.0 -3563.5 1817.0  -3634.0                        
## m107 12 -3610.2 -3557.2 1817.1  -3634.2 0.154      1     0.6947
allmorph_model <- lmer(mass_c ~ beak_c * thorax_c + beak_c * body_c + beak_c * wing_c + thorax_c * body_c + (1|ID), data=data_tested, REML=FALSE)
tidy_regression(allmorph_model, is_color=output_col) 
## lmer mass_c ~ beak_c * thorax_c + beak_c * body_c + beak_c * wing_c + thorax_c * body_c + (1 | ID) data_tested FALSE 
## AIC:  -3612.036 -3612.036 
## (Intercept)      coeff:  -0.0025323  tval:  -2.948129
## beak_c           coeff:  0.0096482   tval:  7.273435
## thorax_c         coeff:  0.0132373   tval:  2.150457
## body_c           coeff:  0.021193    tval:  5.300612
## wing_c           coeff:  -0.0202512  tval:  -5.093224
## beak_c:thorax_c  coeff:  -0.0143255  tval:  -2.446847
## beak_c:body_c    coeff:  0.0095108   tval:  2.96648
## beak_c:wing_c    coeff:  -0.0163177  tval:  -4.673686
## thorax_c:body_c  coeff:  0.02382 tval:  7.378229
allmorph_model <- lmer(mass_c ~ beak_c * thorax_c + beak_c * body_c + beak_c * wing_c + thorax_c * body_c + (1|ID) + (1|trial_type), data=data_tested, REML=FALSE)
tidy_regression(allmorph_model, is_color=output_col) 
## lmer mass_c ~ beak_c * thorax_c + beak_c * body_c + beak_c * wing_c + thorax_c * body_c + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -3621.147 -3621.147 
## (Intercept)      coeff:  -0.0023988  tval:  -1.744651
## beak_c           coeff:  0.0096144   tval:  7.243959
## thorax_c         coeff:  0.0135746   tval:  2.204269
## body_c           coeff:  0.0210139   tval:  5.252242
## wing_c           coeff:  -0.0201397  tval:  -5.061703
## beak_c:thorax_c  coeff:  -0.014353   tval:  -2.451471
## beak_c:body_c    coeff:  0.0094697   tval:  2.952669
## beak_c:wing_c    coeff:  -0.0161723  tval:  -4.629564
## thorax_c:body_c  coeff:  0.0237478   tval:  7.351297
data<-data.frame(R=data_tested$mass_c, 
                 A=data_tested$beak_c, 
                 B=data_tested$thorax_c, 
                 C=data_tested$wing2body_c,
                 X=data_tested$ID)
data <- data %>%
  filter(!is.na(A))

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R"))
##        [,1]      [,2]      [,3]      [,4]     
## AICs   -3556.68  -3556.605 -3554.686 -3554.613
## models 14        12        17        16       
## probs  0.3660182 0.3524344 0.1350214 0.1301697
## 
## m14  R ~ A * B + A * C + (1 | X)
## m12  R ~ A * C + B + (1 | X)
## m17  R ~ A * B + A * C + B * C + (1 | X)
## m16  R ~ A * C + B * C + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings))
## Number of models that failed to converge:  0
anova(m12, m14, test="Chisq") # Adding A*B does not improve fit
## Data: data
## Models:
## m12: R ~ A * C + B + (1 | X)
## m14: R ~ A * B + A * C + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m12  7 -3556.6 -3525.7 1785.3  -3570.6                         
## m14  8 -3556.7 -3521.4 1786.3  -3572.7 2.0756      1     0.1497
allmorph_model <- lmer(mass_c ~  beak_c * wing2body_c + thorax_c + (1 | ID), data=data_tested, REML = FALSE)
tidy_regression(allmorph_model, is_color=output_col) 
## lmer mass_c ~ beak_c * wing2body_c + thorax_c + (1 | ID) data_tested FALSE 
## AIC:  -3556.605 -3556.605 
## (Intercept)          coeff:  -0.0006169  tval:  -0.8517232
## beak_c               coeff:  0.0086669   tval:  7.367686
## wing2body_c          coeff:  -0.2330116  tval:  -5.460915
## thorax_c             coeff:  0.0377287   tval:  9.715388
## beak_c:wing2body_c   coeff:  -0.1777679  tval:  -4.020089
allmorph_model <- lmer(mass_c ~  beak_c * wing2body_c + thorax_c + (1 | ID) + (1|trial_type), data=data_tested, REML = FALSE)
tidy_regression(allmorph_model, is_color=output_col) 
## lmer mass_c ~ beak_c * wing2body_c + thorax_c + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -3566.212 -3566.212 
## (Intercept)          coeff:  -0.0004484  tval:  -0.3390396
## beak_c               coeff:  0.0086493   tval:  7.351589
## wing2body_c          coeff:  -0.2322143  tval:  -5.440538
## thorax_c             coeff:  0.0377582   tval:  9.72031
## beak_c:wing2body_c   coeff:  -0.1752147  tval:  -3.961331

  • positive weak effect of beak, where the longer beak the heavier the bug

  • negative effect of wing2body, where the longer the wing to body the lighter the bug

  • positive effect of thorax, where the longer the thorax the heavier the bug

  • negative beak*wing2body interaction, where the longer the beak and longer the wing2body, the lighter the bug

All Data Plotting

Days from start by mass

source("src/plotting-lm.R") 
## Loading required package: rstan
## Loading required package: StanHeaders
## rstan (Version 2.19.3, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Loading required package: parallel
## Loading required package: dagitty
## rethinking (Version 2.00)
## 
## Attaching package: 'rethinking'
## The following objects are masked from 'package:mosaic':
## 
##     logit, resample
## The following object is masked from 'package:stats':
## 
##     rstudent
lm_plot(data_tested, "days_from_start_c","mass_c", "Days From Start", "Mass (g)", "Mass by Days From Start") 

tidy_regression(lmer(mass_c~days_from_start_c + (1|ID), data=data_tested, REML=FALSE), is_color=output_col)
## lmer mass_c ~ days_from_start_c + (1 | ID) data_tested FALSE 
## AIC:  -3172.693 -3172.693 
## (Intercept)          coeff:  8.98e-05    tval:  0.06844824
## days_from_start_c    coeff:  0.0002329   tval:  3.941263

Flight Duration by Mass

#pdf("massby-flightduration.pdf") # alternative way to save graph
flight_p1 <- as.grob(expression(
  plot(data_tested$mass,data_tested$minute_duration, 
          main = "All Trials",
          xlab = "Mass (g)", 
          ylab = "Flight duration (min)") %>% 
  # Seems like bugs with a mass between 0.025-0.075 g can hit longer flight durations
  abline(v=0.025, col="red") %>%
  abline(v=0.075, col="red")              
#dev.off() # alt way to save graph
))
yes_fly <- filter(data_tested, flew_b == 1)
flight_p2 <- as.grob(expression(
  plot(yes_fly$mass,yes_fly$minute_duration, 
                  main = "All Trials, only bugs that flew (yes_flew)",
                  xlab="Mass (g)",
                  ylab = "Flight duration(min)") %>%
  abline(v=0.025, col="red")%>%
  abline(v=0.075, col="red")
))

C_fly <- filter(yes_fly, flight_type =="C")
flight_p3 <- as.grob(expression(
  plot(C_fly$mass,C_fly$minute_duration,
                  main = "All Trials, only bugs that flew continuously",
                  xlab="Mass (g)",
                  ylab = "Flight duration(min)") %>%
  abline(v=0.025, col="red") %>%
  abline(v=0.075, col="red")
))

grob1 <- as.grob(expression(flight_p1))

grid.arrange(flight_p1, flight_p2, flight_p3, ncol=3)

Plotting mass vs. sym_dist * host * sex

# filter out missing masses
data_tested <- data_tested %>%
  filter(!is.na(mass_c))

p <- gf_point(mass ~ sym_dist, data = data_tested,
                  color = ~host_plant, 
                  alpha = ~sex, 
                  ylab="Mass (g)", 
                  xlab="Distance From Sympatric Zone")
p

#, fig.width=6, fig.height=2.5}
summary<-aggregate(mass~sex*host_plant*sym_dist, data=data_all, FUN=mean)

# Uncomment below to save graph

#pdf("lm-mass-model.pdf", width=9, height=6)

plot(summary$mass~summary$sym_dist, 
     col=c(1,2)[as.factor(summary$sex)], # Female = Black, Male = Red
     pch=c(19,22)[as.factor(summary$host_plant)],
     main="Observed Data",
     xlab = "Distance from Sympatric Zone (°)",
     ylab= "Mass (g)", # K. elegans = Squares C.corindum = circles
     ) 
legend("topright",
       legend = c("C.corindum and F","K.elegans and M"),
       col=c(1,2),
       pch = c(19,22))

#dev.off()

Plotting best-fit model (m38) | lmer(mass ~ sex)

#, fig.width=2.8, fig.height=2.3}
# Uncomment below to save graphs 
# For both graphs: need to add xaxt="n" inside plot function to make the customized x-axis

#pdf("lm-massbysex-model.pdf", width=6, height=6)
source("src/plotting-lm.R") 
source("src/plotting-lm2.R") 
lm_plot2(data_tested, "sex_c","mass_c", "Sex", "Mass (g)", "Mass by Sex") 
axis(1, at = seq(-1,1,.5),labels=c("M", " ", " ", " ", "F"))

#dev.off()

#pdf("lm-massbyhost-model.pdf", width=6, height=6)
lm_plot2(data_tested, "host_c","mass_c", "Host Plant", "Mass (g)", "Mass by Host Plant") 
axis(1, at = seq(-1,1,.5),labels=c("BV", " ", " ", " ", "GRT"))

#dev.off()

lm_plot(data_tested, "sym_dist","mass_c", "Distance From Sympatric Zone", "Mass (g)", "Mass by Sym_Dist")

Morphology Plots

The following plots are nonlinear - it’s worth exploring this nonlinear relationship with a polynomial regression.

p <- gf_point(mass ~ wing, data = data_tested,
                  color = ~thorax, 
                  alpha = ~body, 
                  size = ~beak,
                  ylab="Mass (g)", 
                  xlab="Wing (mm)")
p

#ggsave("MassbyMorphology-alldata_gf", plot = p, device="jpeg", width = 9, height = 7)

First-Order Regression

p1 <- gf_point(mass ~ wing2body, data = data_tested,
                  ylab="Mass (g)", 
                  xlab="Wing-to-Body Ratio") %>%
  gf_lm(alpha=0.2)
p2 <- gf_point(mass ~ wing2body, data = data_tested,
                  color = ~thorax, 
                  alpha = ~beak, 
                  ylab="Mass (g)", 
                  xlab="Wing-to-Body Ratio") %>%
  gf_lm(alpha=0.2)
p3 <- gf_point(mass ~ thorax, data = data_tested,
                  color = ~wing2body, 
                  alpha = ~beak, 
                  ylab="Mass (g)", 
                  xlab="Thorax (mm)") %>%
  gf_lm(alpha=0.2)
p4 <- gf_point(mass ~ beak, data = data_tested,
                  color = ~wing2body, 
                  alpha = ~thorax, 
                  ylab="Mass (g)", 
                  xlab="Beak (mm)") %>%
  gf_lm(alpha=0.2)
g <- grid.arrange(p1,p2,p3,p4, ncol=2, top="All Data")

#ggsave("MassbyMorphology-alldata2", plot = g, device="jpeg", width = 9, height = 7)
#ggsave("MassbyMorphology-alldata2-lm", plot = g, device="jpeg", width = 12, height = 7)
p5 <- as.grob(expression(
  plot(data_tested$body, data_tested$mass, ylab="Mass (g)", xlab="Body (mm)")))
p6 <- as.grob(expression(
  plot(data_tested$beak, data_tested$mass, ylab="Mass (g)", xlab="Beak (mm)")))
p7 <- as.grob(expression(
  plot(data_tested$thorax, data_tested$mass, ylab="Mass (g)", xlab="Throax (mm)")))
p8 <- as.grob(expression(
  plot(data_tested$wing, data_tested$mass, ylab="Mass (g)", xlab="Wing (mm)")))

g2 <- grid.arrange(p5,p6,p7,p8, ncol=2, top="All Data")

# Uncomment below to save graph:
#ggsave("MassbyMorphology-alldata", plot = g2, device="jpeg", width = 9, height = 7)

p9 <- plot(data_tested$wing2body, data_tested$mass, ylab="Mass (g)", xlab="Wing-to-Body Ratio", main="All Data")

d <- data_tested %>%
  filter(!is.na(mass)) %>%
  filter(!is.na(body))
#, fig.width=5, fig.height=4}
source("src/plotting-lm.R") 
lm_plot(d, "wing2body_c","mass_c", "Wing-to-Body Ratio", "Mass (g)", "Mass by Wing:Body") 

par(mfrow=c(2,2))
lm_plot(d, "body_c","mass_c", "Body Length (mm)", "Mass (g)", "Mass by Body Length") 
lm_plot(d, "wing_c","mass_c", "Wing Length (mm)", "Mass (g)", "Mass by Wing Length") 
lm_plot(d, "thorax_c","mass_c", "Thorax Length (mm)", "Mass (g)", "Mass by Thorax Length") 
lm_plot(d, "beak_c","mass_c", "Beak Length (mm)", "Mass (g)", "Mass by Beak Length")

Second-Order Polynomial Regression

#, fig.width=5, fig.height=4}
source("src/plotting-polyreg.R")
lm_polyplot(d, "wing2body_c","mass_c", "Wing-to-Body Ratio", "Mass (g)", "Mass by Wing:Body") 

#par(mfrow=c(2,2))
lm_polyplot(d, "body", "mass", "Body Length (mm)", "Mass (g)", "Mass by Body Length")

lm_polyplot(d, "thorax", "mass", "Thorax Length (mm)", "Mass (g)", "Mass by Thorax Length") # check some bugs 

lm_polyplot(d, "wing", "mass", "Wing Length (mm)", "Mass (g)", "Mass by Wing Length")

lm_polyplot(d, "beak", "mass", "Beak Length (mm)", "Mass (g)", "Mass by Beak Length")

More gf_point explanations: https://cran.r-project.org/web/packages/ggformula/vignettes/ggformula-blog.html

Try an exponential line next: https://www.theanalysisfactor.com/r-tutorial-5/

Female vs. Male Histograms

# Had twice as many males as females.
data_fem <- filter(data_tested, sex == "F") 
data_male <- filter(data_tested, sex == "M") 

data_fem1 <- filter(data_tested, sex == "F", trial_type=="T1") 
data_fem2 <- filter(data_tested, sex == "F", trial_type=="T2")  
data_male1 <- filter(data_tested, sex == "M", trial_type=="T1")  
data_male2 <- filter(data_tested, sex == "M", trial_type=="T2") 

h3 <- as.grob(expression(hist(data_fem1$mass, main="Females Trial 1", xlab= "Mass (g)")))
h4 <- as.grob(expression(hist(data_fem2$mass, main="Females Trial 2", xlab= "Mass (g)")))
h5 <- as.grob(expression(hist(data_male1$mass, main="Males Trial 1", xlab = "Mass (g)")))
h6 <- as.grob(expression(hist(data_male2$mass, main="Males Trial 2", xlab = "Mass (g)")))

h7 <- as.grob(expression(hist(data_fem$mass, main="Females", xlab= "Mass (g)")))
h8 <- as.grob(expression(hist(data_male$mass, main="Males", xlab= "Mass (g)")))

grid.arrange(h7,h8,ncol=2)

grid.arrange(h3,h4,h5,h6,ncol=4)

Females

Cleaning Data

Testing Covariates

lmer() A=host, B=sex, C=sym_dist, X=ID, Y=trial_type, Z=chamber

Morphology

Egg and Morphology

Reading the data

rm(list=ls())
output_col = FALSE 
source("src/clean_flight_data.R") 
source("src/regression_output.R") 
source("src/center_flight_data.R")
source("src/get_warnings.R")

data <- read_flight_data("data/all_flight_data-Winter2020.csv")
data_all <- data[[1]]
data_tested <- data[[2]]

# filter out missing masses
data_tested <- data_tested %>%
  filter(!is.na(mass))

data_fem <- data_tested[data_tested$sex=="F",]
data_fem <- center_data(data_fem)

Testing Covariates

Experimental Set-Up Effects

## lmer mass_c ~ chamber + (1 | ID) + (1 | trial_type) data_fem FALSE 
## AIC:  -1005.932 -1005.932 
## (Intercept)  coeff:  0.0020167   tval:  0.4078058
## chamberA-2   coeff:  -0.0024752  tval:  -0.4380517
## chamberA-3   coeff:  0.001937    tval:  0.3217407
## chamberA-4   coeff:  -0.0048502  tval:  -0.8921571
## chamberB-1   coeff:  -0.0015052  tval:  -0.2414935
## chamberB-2   coeff:  -0.0070746  tval:  -1.192351
## chamberB-3   coeff:  -0.0033839  tval:  -0.6096036
## chamberB-4   coeff:  -0.0002183  tval:  -0.03900598
## lmer mass_c ~ days_from_start_c + (1 | ID) data_fem FALSE 
## AIC:  -1019.432 -1019.432 
## (Intercept)          coeff:  -0.0005159  tval:  -0.2532834
## days_from_start_c    coeff:  0.0003761   tval:  2.328514
## lmer mass_c ~ min_from_IncStart_c + (1 | ID) + (1 | trial_type) data_fem FALSE 
## AIC:  -1014.35 -1014.35 
## (Intercept)          coeff:  -0.000414   tval:  -0.1610887
## min_from_IncStart_c  coeff:  5.5e-06 tval:  0.5800358

Biological Effects

## lmer mass_c ~ total_eggs + (1 | ID) + (1 | trial_type) data_fem FALSE 
## AIC:  -915.1807 -915.1807 
## (Intercept)  coeff:  -0.0081024  tval:  -2.778514
## total_eggs   coeff:  0.0001081   tval:  5.462349
## lmer mass_c ~ eggs_b + (1 | ID) data_fem FALSE 
## AIC:  -1048.447 -1048.447 
## (Intercept)  coeff:  -0.0115389  tval:  -4.449116
## eggs_b       coeff:  0.0169169   tval:  6.110421

Morphology Effects

## lmer mass_c ~ beak_c + (1 | ID) + (1 | trial_type) data_fem FALSE 
## AIC:  -1028.029 -1028.029 
## (Intercept)  coeff:  -0.0003137  tval:  -0.1243193
## beak_c       coeff:  0.0089663   tval:  3.845114
## lmer mass_c ~ thorax_c + (1 | ID) + (1 | trial_type) data_fem FALSE 
## AIC:  -1050.582 -1050.582 
## (Intercept)  coeff:  -0.0002928  tval:  -0.122799
## thorax_c     coeff:  0.0413919   tval:  6.546765
## lmer mass_c ~ body_c + (1 | ID) data_fem FALSE 
## AIC:  -1043.769 -1043.769 
## (Intercept)  coeff:  -0.0004848  tval:  -0.2670899
## body_c       coeff:  0.0102721   tval:  5.784586
## lmer mass_c ~ wing_c + (1 | ID) data_fem FALSE 
## AIC:  -1034.742 -1034.742 
## (Intercept)  coeff:  -0.0005514  tval:  -0.2930797
## wing_c       coeff:  0.0102625   tval:  4.74619

Female Data | lmer() A=host, B=sex, C=sym_dist, X=ID

data<-data.frame(R=data_fem$mass_c, 
                 A=data_fem$host_c, 
                 B=data_fem$sym_dist,
                 X=data_fem$ID)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 2-FF REMLF.R"))
##        [,1]      [,2]      [,3]      [,4]     [,5]     
## AICs   -1018.297 -1017.899 -1017.107 -1016.95 -1016.082
## models 2         1         3         4        5        
## probs  0.3113621 0.2552403 0.1717577 0.158744 0.1028958
## 
## m2   R ~ B + (1 | X)
## m1   R ~ A + (1 | X)
## m3   R ~ A + B + (1 | X)
## m4   R ~ A * B + (1 | X)
## m0   R ~ (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings))
## Number of models that failed to converge:  0
anova(m2, m3, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m2: R ~ B + (1 | X)
## m3: R ~ A + B + (1 | X)
##    Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2  4 -1018.3 -1004.9 513.15  -1026.3                         
## m3  5 -1017.1 -1000.3 513.55  -1027.1 0.8103      1      0.368
femm_model <- lmer(mass_c ~ sym_dist + (1|ID), data=data_fem, REML=FALSE)
tidy_regression(femm_model, is_color=output_col)
## lmer mass_c ~ sym_dist + (1 | ID) data_fem FALSE 
## AIC:  -1018.297 -1018.297 
## (Intercept)  coeff:  0.0025664   tval:  1.010654
## sym_dist     coeff:  -0.0048665  tval:  -2.069857
femm_model <- lmer(mass_c ~ sym_dist + (1|ID) + (1|trial_type), data=data_fem, REML=FALSE)
tidy_regression(femm_model, is_color=output_col)
## lmer mass_c ~ sym_dist + (1 | ID) + (1 | trial_type) data_fem FALSE 
## AIC:  -1017.961 -1017.961 
## (Intercept)  coeff:  0.0026651   tval:  0.8957712
## sym_dist     coeff:  -0.0047023  tval:  -2.001425
  • marginal negative weak effect of sym_dist where the farther from homestead the less GRT female bugs weigh

Let’s add eggs:

data<-data.frame(R=data_fem$mass_c, 
                 A=data_fem$host_c, 
                 B=data_fem$sym_dist,
                 C=data_fem$eggs_b,
                 X=data_fem$ID)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R"))
##        [,1]      [,2]      [,3]      [,4]      [,5]       [,6]       [,7]      
## AICs   -1049.615 -1049.147 -1048.447 -1048.13  -1047.847  -1047.635  -1047.194 
## models 6         5         3         7         11         10         9         
## probs  0.2128016 0.1683949 0.1186843 0.1012758 0.08791912 0.07907821 0.06342284
## 
## m6   R ~ B + C + (1 | X)
## m5   R ~ A + C + (1 | X)
## m3   R ~ C + (1 | X)
## m7   R ~ A + B + C + (1 | X)
## m11  R ~ A * B + C + (1 | X)
## m10  R ~ B * C + (1 | X)
## m9   R ~ A * C + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings))
## Number of models that failed to converge:  0
anova(m6, m10, test="Chisq") # Adding B*C does not improve fit
## Data: data
## Models:
## m6: R ~ B + C + (1 | X)
## m10: R ~ B * C + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m6   5 -1049.6 -1032.8 529.81  -1059.6                         
## m10  6 -1047.6 -1027.5 529.82  -1059.6 0.0202      1     0.8871
femm_model <- lmer(mass_c ~ sym_dist + eggs_b + (1|ID), data=data_fem, REML=FALSE)
tidy_regression(femm_model, is_color=output_col) # When addd trial type as RF, it has zero variance so not needed
## lmer mass_c ~ sym_dist + eggs_b + (1 | ID) data_fem FALSE 
## AIC:  -1049.615 -1049.615 
## (Intercept)  coeff:  -0.0087915  tval:  -2.927234
## sym_dist     coeff:  -0.0038795  tval:  -1.790062
## eggs_b       coeff:  0.016602    tval:  6.00753
  • no efffect of sym dist
  • Positive effect of eggs laid on test day where if a female bug laid eggs on her test day, she weighed more.
data<-data.frame(R=data_fem$mass_c, 
                 A=data_fem$host_c, 
                 B=data_fem$sym_dist,
                 C=data_fem$total_eggs_c,
                 X=data_fem$ID)
data <- data %>%
  filter(!is.na(D))
## Warning in is.na(D): is.na() applied to non-(list or vector) of type 'closure'
source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R"))
##        [,1]      [,2]      [,3]      [,4]     
## AICs   -1018.297 -1017.899 -1017.107 -1016.95 
## models 2         1         4         8        
## probs  0.3468692 0.2843474 0.1913446 0.1768469
## 
## m2   R ~ B + (1 | X)
## m1   R ~ A + (1 | X)
## m4   R ~ A + B + (1 | X)
## m8   R ~ A * B + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings))
## Number of models that failed to converge:  1
anova(m2, m4, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m2: R ~ B + (1 | X)
## m4: R ~ A + B + (1 | X)
##    Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2  4 -1018.3 -1004.9 513.15  -1026.3                         
## m4  5 -1017.1 -1000.3 513.55  -1027.1 0.8103      1      0.368
femm_model <- lmer(mass_c ~ sym_dist + (1|ID), data=data_fem, REML=FALSE)
tidy_regression(femm_model, is_color=output_col)
## lmer mass_c ~ sym_dist + (1 | ID) data_fem FALSE 
## AIC:  -1018.297 -1018.297 
## (Intercept)  coeff:  0.0025664   tval:  1.010654
## sym_dist     coeff:  -0.0048665  tval:  -2.069857
femm_model <- lmer(mass_c ~ sym_dist + (1|ID) + (1 | trial_type), data=data_fem, REML=FALSE)
tidy_regression(femm_model, is_color=output_col)
## lmer mass_c ~ sym_dist + (1 | ID) + (1 | trial_type) data_fem FALSE 
## AIC:  -1017.961 -1017.961 
## (Intercept)  coeff:  0.0026651   tval:  0.8957712
## sym_dist     coeff:  -0.0047023  tval:  -2.001425
  • negative marginal effect of sym_dist where the farther the bug is from Homestead, the less it weighed

Morphology Modeling

data<-data.frame(R=data_fem$mass_c, 
                 A=data_fem$beak_c, 
                 B=data_fem$thorax_c,
                 C=data_fem$wing_c,
                 D=data_fem$body_c,
                 X=data_fem$ID) 
data <- data %>%
  filter(!is.na(A))

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF REMLF.R"))
##        [,1]      
## AICs   -1063.81  
## models 58        
## probs  0.08469531
## 
## m58  R ~ A * B + B * D + C + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings))
## Number of models that failed to converge:  0
femm_morph_model <- lmer(mass_c ~ beak_c * thorax_c + thorax_c * body_c + wing_c  + (1|ID), data=data_fem, REML=FALSE)
tidy_regression(femm_morph_model, is_color=output_col)
## lmer mass_c ~ beak_c * thorax_c + thorax_c * body_c + wing_c + (1 | ID) data_fem FALSE 
## AIC:  -1063.81 -1063.81 
## (Intercept)      coeff:  -0.0031752  tval:  -1.530044
## beak_c           coeff:  0.0007551   tval:  0.2877745
## thorax_c         coeff:  0.016121    tval:  1.100156
## body_c           coeff:  0.0324877   tval:  3.195056
## wing_c           coeff:  -0.0303894  tval:  -2.873455
## beak_c:thorax_c  coeff:  -0.0238294  tval:  -2.562146
## thorax_c:body_c  coeff:  0.0240542   tval:  3.123057
femm_morph_model <- lmer(mass_c ~ beak_c * thorax_c + thorax_c * body_c + wing_c  + (1|ID) + (1|trial_type), data=data_fem, REML=FALSE)
tidy_regression(femm_morph_model, is_color=output_col)
## lmer mass_c ~ beak_c * thorax_c + thorax_c * body_c + wing_c + (1 | ID) + (1 | trial_type) data_fem FALSE 
## AIC:  -1063.408 -1063.408 
## (Intercept)      coeff:  -0.0030159  tval:  -1.188039
## beak_c           coeff:  0.000773    tval:  0.2949609
## thorax_c         coeff:  0.0162314   tval:  1.108989
## body_c           coeff:  0.0318749   tval:  3.136705
## wing_c           coeff:  -0.0297012  tval:  -2.809744
## beak_c:thorax_c  coeff:  -0.0237685  tval:  -2.557551
## thorax_c:body_c  coeff:  0.0241126   tval:  3.134129

  • no effect of beak

  • no effect of thorax

  • positive effect of body

  • negative effect of beak

  • negative effect of wing

  • beak*thorax negative effect where the longer the beak and throax the less the bug weighs

  • positive effect of thorax*body where the larger the thorax and body, more less the bug weighs

Eggs and Morphology

data<-data.frame(R=data_fem$mass_c, 
                 A=data_fem$eggs_b, 
                 B=data_fem$thorax_c,
                 C=data_fem$wing_c,
                 D=data_fem$body_c,
                 X=data_fem$ID) 
data <- data %>%
  filter(!is.na(A), !is.na(B)) 

## model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF.R") ## AB: wide difference between model selection, uncomment to compare.

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF REMLF.R"))
##        [,1]      
## AICs   -1106.339 
## models 62        
## probs  0.04304825
## 
## m62  R ~ A * D + B * D + C + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings)) 
## Number of models that failed to converge:  0
femm_eggmorph <- lmer(mass_c ~ eggs_b * body_c + thorax_c * body_c + wing_c  + (1|ID), data=data_fem, REML=TRUE)
tidy_regression(femm_eggmorph, is_color=output_col) # Adding (1|trial_type) does not improve fit
## lmer mass_c ~ eggs_b * body_c + thorax_c * body_c + wing_c + (1 | ID) data_fem TRUE 
## AIC:  -1039.952 -1057.952 
## (Intercept)      coeff:  -0.0134964  tval:  -5.494663
## eggs_b           coeff:  0.0176765   tval:  6.912355
## body_c           coeff:  0.0267177   tval:  2.959187
## thorax_c         coeff:  0.0241161   tval:  1.876735
## wing_c           coeff:  -0.0294202  tval:  -3.248083
## eggs_b:body_c    coeff:  0.0043137   tval:  1.626413
## body_c:thorax_c  coeff:  0.0073997   tval:  1.627234

  • if female laid eggs that day, then mass is larger

  • the larger the body, the heavier the female

  • marginal? effect of thorax

  • negative effect of wing where the larger the wing the lighter the female

  • no effect of eggs*body

  • no effect of body*thorax

Plotting Female Data

Sym Distance

p <- gf_point(mass ~ sym_dist, data = data_fem, 
                  color = ~host_plant, 
                  ylab="Mass (g)", 
                  xlab="Distance From Sympatric Zone")
p

Eggs

#, fig.width=2.8, fig.height=2.3}
data_eggs <- data_fem %>%
  filter(!is.na(total_eggs))
source("src/plotting-lm.R") 
lm_plot(data_eggs, "total_eggs_c", "mass_c", "Total Eggs Laid", "Mass (g)", "Mass by Total Eggs Laid")

source("src/plotting-lm2.R") 
lm_plot2(data_fem, "eggs_b", "mass_c", "Eggs Laid", "Mass (g)", "Mass by Egg Laid on Test Day")
axis(1, at = seq(-1,1,.5),labels=c(" ", " ", "M ", " ", "F"))

lm_plot(data_fem, "sym_dist", "mass_c", "Distance From Sympatric Zone", "Mass (g)", "Mass by Distance from Sympatric Zone")

Morphology

#, fig.width=2.8, fig.height=2.3}

source("src/plotting-lm.R") 
lm_plot(data_fem, "wing2body_c", "mass_c", "Wing:Body", "Mass (g)", "Mass by Wing:Body")

lm_plot(data_fem, "thorax_c", "mass_c", "Thorax Length (mm)", "Mass (g)", "Mass by Thorax Length")

lm_plot(data_fem, "wing_c", "mass_c", "Wing Length (mm)", "Mass (g)", "Mass by Wing Length")

lm_plot(data_fem, "body_c", "mass_c", "Body Length (mm)", "Mass (g)", "Mass by Body Length")

lm_plot(data_fem, "beak_c", "mass_c", "Beak Length (mm)", "Mass (g)", "Mass by Beak Length")

source("src/plotting-polyreg.R")
lm_polyplot(data_fem, "wing2body_c", "mass_c", "Wing:Body", "Mass (g)", "Mass by Wing:Body")

lm_polyplot(data_fem, "thorax_c", "mass_c", "Thorax Length (mm)", "Mass (g)", "Mass by Thorax Length")

lm_polyplot(data_fem, "wing_c", "mass_c", "Wing Length (mm)", "Mass (g)", "Mass by Wing Length")

lm_polyplot(data_fem, "body_c", "mass_c", "Body Length (mm)", "Mass (g)", "Mass by Body Length")

lm_polyplot(data_fem, "beak_c", "mass_c", "Beak Length (mm)", "Mass (g)", "Mass by Beak Length")

Males

Cleaning Data

Testing Covariates

lmer() A=host, B=sym_dist, X=ID, Y=trial_type

Morphology1

Reading the data

rm(list=ls())
output_col = FALSE 
source("src/clean_flight_data.R") 
source("src/regression_output.R") 
source("src/center_flight_data.R")
source("src/get_warnings.R")

data <- read_flight_data("data/all_flight_data-Winter2020.csv")
data_all <- data[[1]]
data_tested <- data[[2]]

data_male <- data_tested[data_tested$sex=="M",]
data_male <- center_data(data_male)

Testing Covariates

Experimental Set-Up Effects

## lmer mass_c ~ chamber + (1 | ID) data_male FALSE 
## AIC:  -2892.297 -2892.297 
## (Intercept)  coeff:  0.0003444   tval:  0.3731964
## chamberA-2   coeff:  -1.21e-05   tval:  -0.0101808
## chamberA-3   coeff:  -0.0009262  tval:  -0.8482127
## chamberA-4   coeff:  0.000511    tval:  0.4739871
## chamberB-1   coeff:  -0.0014291  tval:  -1.153582
## chamberB-2   coeff:  -0.0004232  tval:  -0.3988591
## chamberB-3   coeff:  -0.0015228  tval:  -1.303421
## chamberB-4   coeff:  -0.0010191  tval:  -0.8764588
## lmer mass_c ~ days_from_start_c + (1 | ID) + (1 | trial_type) data_male FALSE 
## AIC:  -2931.511 -2931.511 
## (Intercept)          coeff:  -4.29e-05   tval:  -0.05729814
## days_from_start_c    coeff:  8.06e-05    tval:  1.304836
## lmer mass_c ~ min_from_IncStart_c + (1 | ID) data_male FALSE 
## AIC:  -2898.075 -2898.075 
## (Intercept)          coeff:  -0.0001518  tval:  -0.3297338
## min_from_IncStart_c  coeff:  -4e-07  tval:  -0.2178442

Biological Effects - none b/c males don’t lay eggs

Morphology Effects

## lmer mass_c ~ beak_c + (1 | ID) + (1 | trial_type) data_male FALSE 
## AIC:  -3015.466 -3015.466 
## (Intercept)  coeff:  8.82e-05    tval:  0.09159024
## beak_c       coeff:  0.0095844   tval:  10.17221
## lmer mass_c ~ thorax_c + (1 | ID) + (1 | trial_type) data_male FALSE 
## AIC:  -3105.049 -3105.049 
## (Intercept)  coeff:  4.25e-05    tval:  0.04404148
## thorax_c     coeff:  0.024729    tval:  16.51937
## lmer mass_c ~ body_c + (1 | ID) + (1 | trial_type) data_male FALSE 
## AIC:  -3089.733 -3089.733 
## (Intercept)  coeff:  7.33e-05    tval:  0.0779836
## body_c       coeff:  0.0069936   tval:  15.40888
## lmer mass_c ~ wing_c + (1 | ID) + (1 | trial_type) data_male FALSE 
## AIC:  -3051.742 -3051.742 
## (Intercept)  coeff:  5.95e-05    tval:  0.0618658
## wing_c       coeff:  0.0074195   tval:  12.79243

Male Data | lmer() A=host, B=sym_dist, X=ID, Y=trial_type

data<-data.frame(R=data_male$mass_c, 
                 A=data_male$host_c, 
                 B=data_male$sym_dist,
                 X=data_male$ID)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 2-FF REMLF.R"))
##        [,1]      [,2]      [,3]     
## AICs   -2912.775 -2912.531 -2910.353
## models 4         3         1        
## probs  0.4574682 0.4049383 0.1362975
## 
## m4   R ~ A * B + (1 | X)
## m3   R ~ A + B + (1 | X)
## m1   R ~ A + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings)) 
## Number of models that failed to converge:  0
anova(m3, m4, test="Chisq") # Adding A*B does not improve fit
## Data: data
## Models:
## m3: R ~ A + B + (1 | X)
## m4: R ~ A * B + (1 | X)
##    Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m3  5 -2912.5 -2892.6 1461.3  -2922.5                         
## m4  6 -2912.8 -2888.9 1462.4  -2924.8 2.2439      1     0.1341
male_model <- lmer(mass_c ~ host_c + sym_dist + (1 | ID), data=data_male, REML=FALSE)
tidy_regression(male_model, is_color=output_col)
## lmer mass_c ~ host_c + sym_dist + (1 | ID) data_male FALSE 
## AIC:  -2912.531 -2912.531 
## (Intercept)  coeff:  -0.0024287  tval:  -2.817653
## host_c       coeff:  -0.0028506  tval:  -3.982671
## sym_dist     coeff:  0.0011128   tval:  2.052907
male_model <- lmer(mass_c ~ host_c + sym_dist + (1 | ID) + (1 | trial_type), data=data_male, REML=FALSE)
tidy_regression(male_model, is_color=output_col)
## lmer mass_c ~ host_c + sym_dist + (1 | ID) + (1 | trial_type) data_male FALSE 
## AIC:  -2944.146 -2944.146 
## (Intercept)  coeff:  -0.002252   tval:  -1.792665
## host_c       coeff:  -0.0027557  tval:  -3.86503
## sym_dist     coeff:  0.0011124   tval:  2.061392
  • significant but weak negative effect of host
  • no effect of sym_dist

Morphology Modeling

data<-data.frame(R=data_male$mass_c, 
                 A=data_male$beak_c, 
                 B=data_male$thorax_c,
                 C=data_male$wing_c,
                 D=data_male$body_c,
                 X=data_male$ID) # multiple models failed if add Y=data_male$trialtype
data <- data %>%
  filter(!is.na(A))

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF REMLF.R"))
##        [,1]      [,2]      [,3]       [,4]       [,5]       [,6]      
## AICs   -3116.99  -3116.818 -3115.772  -3115.759  -3115.604  -3115.585 
## models 90        64        80         104        83         105       
## probs  0.1305409 0.1198163 0.07101086 0.07054803 0.06529412 0.06468077
##        [,7]       [,8]      
## AICs   -3115.23   -3115.116 
## models 97         94        
## probs  0.05416005 0.05115706
## 
## m90  R ~ B * C + B * D + C * D + A + (1 | X)
## m64  R ~ B * C + B * D + A + (1 | X)
## m80  R ~ A * C + B * C + B * D + (1 | X)
## m104     R ~ A * C + B * C + B * D + C * D + (1 | X)
## m83  R ~ A * D + B * C + B * D + (1 | X)
## m105     R ~ A * D + B * C + B * D + C * D + (1 | X)
## m97  R ~ A * B + A * D + B * C + B * D + (1 | X)
## m94  R ~ A * B + A * C + B * C + B * D + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings)) 
## Number of models that failed to converge:  0
# this changed
anova(m64, m90, test="Chisq") # Adding C*D does not improve fit
## Data: data
## Models:
## m64: R ~ B * C + B * D + A + (1 | X)
## m90: R ~ B * C + B * D + C * D + A + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m64  9 -3116.8 -3080.9 1567.4  -3134.8                         
## m90 10 -3117.0 -3077.1 1568.5  -3137.0 2.1715      1     0.1406
anova(m64, m80, test="Chisq") # Adding A*C does not improve fit
## Data: data
## Models:
## m64: R ~ B * C + B * D + A + (1 | X)
## m80: R ~ A * C + B * C + B * D + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m64  9 -3116.8 -3080.9 1567.4  -3134.8                         
## m80 10 -3115.8 -3075.9 1567.9  -3135.8 0.9537      1     0.3288
male_morph_model <- lmer(mass_c ~ thorax_c * wing_c + thorax_c * body_c + wing_c * body_c + beak_c  + (1|ID), data=data_male, REML=FALSE)
tidy_regression(male_morph_model, is_color=output_col)
## lmer mass_c ~ thorax_c * wing_c + thorax_c * body_c + wing_c * body_c + beak_c + (1 | ID) data_male FALSE 
## AIC:  -3116.99 -3116.99 
## (Intercept)      coeff:  -0.0007847  tval:  -2.305689
## thorax_c         coeff:  0.012102    tval:  4.529882
## wing_c           coeff:  -0.0026275  tval:  -1.536531
## body_c           coeff:  0.0056071   tval:  3.247561
## beak_c           coeff:  0.0040832   tval:  5.097926
## thorax_c:wing_c  coeff:  -0.0288963  tval:  -3.126314
## thorax_c:body_c  coeff:  0.0240933   tval:  3.622992
## wing_c:body_c    coeff:  0.0016064   tval:  1.478191
male_morph_model <- lmer(mass_c ~ thorax_c * wing_c + thorax_c * body_c + wing_c * body_c + beak_c  + (1|ID) + (1|trial_type), data=data_male, REML=FALSE)
tidy_regression(male_morph_model, is_color=output_col)
## lmer mass_c ~ thorax_c * wing_c + thorax_c * body_c + wing_c * body_c + beak_c + (1 | ID) + (1 | trial_type) data_male FALSE 
## AIC:  -3147.511 -3147.511 
## (Intercept)      coeff:  -0.0006814  tval:  -0.7211702
## thorax_c         coeff:  0.0124757   tval:  4.682477
## wing_c           coeff:  -0.0026086  tval:  -1.525296
## body_c           coeff:  0.0054843   tval:  3.177486
## beak_c           coeff:  0.0039302   tval:  4.906163
## thorax_c:wing_c  coeff:  -0.0290863  tval:  -3.156053
## thorax_c:body_c  coeff:  0.0242504   tval:  3.652631
## wing_c:body_c    coeff:  0.0015995   tval:  1.474572

Plotting Male Data

data_mmale <- data_male %>%
  filter(!is.na(thorax))
p1 <- gf_point(mass ~ sym_dist, data = data_mmale, 
                  color = ~host_plant, 
                  ylab="Mass (g)", 
                  xlab="Distance From Sympatric Zone")
p1

#, fig.width=2.8, fig.height=2.3}
source("src/plotting-lm.R") 
lm_plot(data_mmale, "wing2body_c", "mass_c", "Wing:Body", "Mass (g)", "Mass by Wing:Body")
## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.

lm_plot(data_mmale, "thorax_c", "mass_c", "Thorax Length (mm)", "Mass (g)", "Mass by Thorax Length")
## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.

lm_plot(data_mmale, "wing_c", "mass_c", "Wing Length (mm)", "Mass (g)", "Mass by Wing Length")
## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.

lm_plot(data_mmale, "body_c", "mass_c", "Body Length (mm)", "Mass (g)", "Mass by Body Length")
## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.

lm_plot(data_mmale, "beak_c", "mass_c", "Beak Length (mm)", "Mass (g)", "Mass by Beak Length")
## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.

# , fig.width=2.8, fig.height=2.3}
p.morph <- gf_point(mass ~ thorax, data = data_mmale, 
                  alpha = ~beak, 
                  ylab="Mass (g)", 
                  xlab="Thoarx Length (mm)")
p.morph

Trial 1

Cleaning Data

Testing Covariates

Sex * Host * Sym_Dist Modeling

Morphology

Reading the data

rm(list=ls())
output_col = FALSE
source("src/clean_flight_data.R") 
source("src/regression_output.R") 
source("src/center_flight_data.R")
source("src/get_warnings.R")

data <- read_flight_data("data/all_flight_data-Winter2020.csv")
data_all <- data[[1]]
data_tested <- data[[2]]

# filter out missing masses
data_tested <- data_tested %>%
  filter(!is.na(mass_c))

data_T1 <- filter(data_tested, trial_type == "T1")
data_T2 <- filter(data_tested, trial_type == "T2")
data_T1 <- center_data(data_T1)

Testing Covariates

Experimental Set-Up Effects

## lmer mass_c ~ chamber + (1 | days_from_start_c) data_T1 FALSE 
## AIC:  -1491.671 -1491.671 
## (Intercept)  coeff:  0.0013799   tval:  0.2101058
## chamberA-2   coeff:  0.0063281   tval:  0.8408197
## chamberA-3   coeff:  -0.0043664  tval:  -0.5931701
## chamberA-4   coeff:  -0.0018981  tval:  -0.2599032
## chamberB-1   coeff:  0.0011979   tval:  0.1538927
## chamberB-2   coeff:  -0.0066758  tval:  -0.912926
## chamberB-3   coeff:  -0.0001199  tval:  -0.01620135
## chamberB-4   coeff:  -0.0023396  tval:  -0.3043403
## lmer mass_c ~ days_from_start_c + (1 | chamber) data_T1 FALSE 
## AIC:  -1497.218 -1497.218 
## (Intercept)          coeff:  4.94e-05    tval:  0.03384704
## days_from_start_c    coeff:  0.0005133   tval:  1.420627
## lmer mass_c ~ min_from_IncStart_c + (1 | chamber) data_T1 FALSE 
## AIC:  -1495.637 -1495.637 
## (Intercept)          coeff:  4.08e-05    tval:  0.0285751
## min_from_IncStart_c  coeff:  6.3e-06 tval:  0.641897

Morphology Effects

## lmer mass_c ~ beak_c + (1 | population) data_T1 FALSE 
## AIC:  -1736.24 -1736.24 
## (Intercept)  coeff:  0.0003952   tval:  0.3474711
## beak_c       coeff:  0.0175105   tval:  18.90378
## lmer mass_c ~ thorax_c + (1 | population) data_T1 FALSE 
## AIC:  -1740.169 -1740.169 
## (Intercept)  coeff:  -0.0011473  tval:  -0.5590813
## thorax_c     coeff:  0.0593184   tval:  19.43148
## lmer mass_c ~ body_c + (1 | population) data_T1 FALSE 
## AIC:  -1719.377 -1719.377 
## (Intercept)  coeff:  -0.0018591  tval:  -0.8962435
## body_c       coeff:  0.0164097   tval:  18.27009
## lmer mass_c ~ wing_c + (1 | population) data_T1 FALSE 
## AIC:  -1646.291 -1646.291 
## (Intercept)  coeff:  -0.001967   tval:  -0.9477477
## wing_c       coeff:  0.01847 tval:  14.13072

Sex * Host * Sym_Dist Modeling

data<-data.frame(R=data_T1$mass_c, 
                 A=data_T1$host_c, 
                 B=data_T1$sym_dist,
                 C=data_T1$sex_c,
                 X=data_T1$population)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R"))
##        [,1]      [,2]      [,3]       [,4]       [,5]       [,6]      
## AICs   -1771.645 -1771.277 -1769.663  -1769.579  -1769.577  -1768.546 
## models 10        13        16         15         5          9         
## probs  0.2608367 0.2170057 0.09682661 0.09282359 0.09277441 0.05538096
## 
## m10  R ~ B * C + (1 | X)
## m13  R ~ B * C + A + (1 | X)
## m16  R ~ A * C + B * C + (1 | X)
## m15  R ~ A * B + B * C + (1 | X)
## m5   R ~ A + C + (1 | X)
## m9   R ~ A * C + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings)) 
## Number of models that failed to converge:  0
anova(m10, m13, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m10: R ~ B * C + (1 | X)
## m13: R ~ B * C + A + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m10  6 -1771.6 -1748.9 891.82  -1783.6                         
## m13  7 -1771.3 -1744.7 892.64  -1785.3 1.6321      1     0.2014
T1_model <- lmer(mass_c ~ sex_c * sym_dist + (1|population), data=data_T1, REML=FALSE)
## boundary (singular) fit: see ?isSingular
tidy_regression(T1_model, is_color=output_col)
## lmer mass_c ~ sex_c * sym_dist + (1 | population) data_T1 FALSE 
## AIC:  -1771.645 -1771.645 
## (Intercept)      coeff:  0.0076033   tval:  6.535934
## sex_c            coeff:  0.0208947   tval:  17.96157
## sym_dist         coeff:  -0.0027465  tval:  -2.854617
## sex_c:sym_dist   coeff:  -0.0023616  tval:  -2.454528

  • positive effect of sex, where if female then weigh more

  • negative effect of sym_dist, where if from GRT then weigh less

  • negative effect of sym_dist*sex where if from GRT and F then weigh less

Morphology Modeling

data<-data.frame(R=data_T1$mass_c, 
                 A=data_T1$beak_c, 
                 B=data_T1$thorax_c,
                 C=data_T1$body_c,
                 D=data_T1$wing_c,
                 X=data_T1$population)
data <- data %>%
  filter(!is.na(A))

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF REMLF.R"))
##        [,1]      [,2]      [,3]       [,4]       [,5]       [,6]      
## AICs   -1852.027 -1851.993 -1851.434  -1851.061  -1850.618  -1850.463 
## models 91        105       110        97         106        107       
## probs  0.1297854 0.1275883 0.09649178 0.08008707 0.06417232 0.05937695
## 
## m91  R ~ A * B + A * C + A * D + B * C + (1 | X)
## m105     R ~ A * D + B * C + B * D + C * D + (1 | X)
## m110     R ~ A * B + A * D + B * C + B * D + C * D + (1 | X)
## m97  R ~ A * B + A * D + B * C + B * D + (1 | X)
## m106     R ~ A * B + A * C + A * D + B * C + B * D + (1 | X)
## m107     R ~ A * B + A * C + A * D + B * C + C * D + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings)) 
## Number of models that failed to converge:  0
T1_morph_model <- lmer(mass_c ~ beak_c * thorax_c + beak_c * body_c + beak_c * wing_c + thorax_c * body_c  + (1 | population), data=data_T1, REML=FALSE)
tidy_regression(T1_morph_model, is_color=output_col)
## lmer mass_c ~ beak_c * thorax_c + beak_c * body_c + beak_c * wing_c + thorax_c * body_c + (1 | population) data_T1 FALSE 
## AIC:  -1852.027 -1852.027 
## (Intercept)      coeff:  -0.0020692  tval:  -1.651501
## beak_c           coeff:  0.0088409   tval:  5.462576
## thorax_c         coeff:  0.0165347   tval:  2.248462
## body_c           coeff:  0.0220664   tval:  4.579629
## wing_c           coeff:  -0.0220187  tval:  -4.6062
## beak_c:thorax_c  coeff:  -0.0154745  tval:  -2.254743
## beak_c:body_c    coeff:  0.0094591   tval:  2.476381
## beak_c:wing_c    coeff:  -0.0168317  tval:  -4.023275
## thorax_c:body_c  coeff:  0.0247599   tval:  6.458283

  • positive effect of beak

  • positive effect of thorax

  • positive effect of wing

  • negative effect of wing

  • marginal negative effect of beak*thorax

  • positve effect of beak*body

  • negative effect of beak*wing

  • positive effect of thorax*body

data<-data.frame(R=data_T1$mass_c, 
                 A=data_T1$beak_c, 
                 B=data_T1$thorax_c,
                 C=data_T1$wing2body_c,
                 X=data_T1$population)
data <- data %>%
  filter(!is.na(A))
source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R"))
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
##        [,1]      [,2]      [,3]      [,4]     
## AICs   -1815.381 -1813.477 -1813.412 -1811.498
## models 12        14        16        17       
## probs  0.5086629 0.1962722 0.1900758 0.072979 
## 
## m12  R ~ A * C + B + (1 | X)
## m14  R ~ A * B + A * C + (1 | X)
## m16  R ~ A * C + B * C + (1 | X)
## m17  R ~ A * B + A * C + B * C + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings)) 
## Number of models that failed to converge:  0
anova(m12, m16, test="Chisq") # Adding B*C interaction does not improve fit
## Data: data
## Models:
## m12: R ~ A * C + B + (1 | X)
## m16: R ~ A * C + B * C + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m12  7 -1815.4 -1788.8 914.69  -1829.4                         
## m16  8 -1813.4 -1783.0 914.71  -1829.4 0.0313      1     0.8596
anova(m12, m14, test="Chisq") # Adding A*B interaction does not improve fit
## Data: data
## Models:
## m12: R ~ A * C + B + (1 | X)
## m14: R ~ A * B + A * C + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m12  7 -1815.4 -1788.8 914.69  -1829.4                         
## m14  8 -1813.5 -1783.1 914.74  -1829.5 0.0954      1     0.7574
T1_morph_model <- lmer(mass_c ~ beak_c * wing2body_c + thorax_c + (1 | population), data=data_T1, REML=FALSE)
tidy_regression(T1_morph_model, is_color=output_col)
## lmer mass_c ~ beak_c * wing2body_c + thorax_c + (1 | population) data_T1 FALSE 
## AIC:  -1815.381 -1815.381 
## (Intercept)          coeff:  -0.0003546  tval:  -0.2118321
## beak_c               coeff:  0.0070722   tval:  5.027234
## wing2body_c          coeff:  -0.2672837  tval:  -5.373914
## thorax_c             coeff:  0.041646    tval:  9.067206
## beak_c:wing2body_c   coeff:  -0.2023521  tval:  -3.911947

  • positive effect of beak where longer then weigh more

  • negative effect of wing2body where the more wing to body then weigh less

  • positive effect of thorax where the wider the thorax the more the bug weighs

  • negative effect of beak*wing2body wher ethe more wing to body and longer beak then weigh less

Trial 2

Cleaning Data

Testing Covariates

Sex * Host * Sym_Dist Modeling

Morphology

Reading the data

rm(list=ls())
output_col = FALSE 
source("src/clean_flight_data.R") 
source("src/regression_output.R") 
source("src/center_flight_data.R")
source("src/get_warnings.R")

data <- read_flight_data("data/all_flight_data-Winter2020.csv")
data_all <- data[[1]]
data_tested <- data[[2]]

# filter out missing masses
data_tested <- data_tested %>%
  filter(!is.na(mass_c))

data_T2 <- filter(data_tested, trial_type == "T2")
data_T2 <- center_data(data_T2)

Testing Experimental Effects

## lmer mass_c ~ chamber + (1 | site) data_T2 FALSE 
## AIC:  -1283.45 -1283.45 
## (Intercept)  coeff:  -0.0008941  tval:  -0.2427346
## chamberA-2   coeff:  -0.003167   tval:  -0.5807747
## chamberA-3   coeff:  -0.0018049  tval:  -0.30813
## chamberA-4   coeff:  0.0012498   tval:  0.256323
## chamberB-1   coeff:  0.0022032   tval:  0.3454594
## chamberB-2   coeff:  -0.0019955  tval:  -0.3579442
## chamberB-3   coeff:  0.0047892   tval:  0.8659809
## chamberB-4   coeff:  0.0056355   tval:  1.05502
## lmer mass_c ~ days_from_start_c + (1 | site) data_T2 FALSE 
## AIC:  -1291.515 -1291.515 
## (Intercept)          coeff:  2.39e-05    tval:  0.01455632
## days_from_start_c    coeff:  0.0002657   tval:  0.5213984
## lmer mass_c ~ min_from_IncStart_c + (1 | site) data_T2 FALSE 
## AIC:  -1292.561 -1292.561 
## (Intercept)          coeff:  2.58e-05    tval:  0.01586435
## min_from_IncStart_c  coeff:  1.18e-05    tval:  1.148711

Morphology Effects

## lmer mass_c ~ beak_c + (1 | population) data_T2 FALSE 
## AIC:  -1586.996 -1586.996 
## (Intercept)  coeff:  0.0009905   tval:  0.7722748
## beak_c       coeff:  0.0196729   tval:  22.98186
## lmer mass_c ~ thorax_c + (1 | population) data_T2 FALSE 
## AIC:  -1549.027 -1549.027 
## (Intercept)  coeff:  -0.001168   tval:  -0.7107158
## thorax_c     coeff:  0.0625115   tval:  20.71668
## lmer mass_c ~ body_c + (1 | population) data_T2 FALSE 
## AIC:  -1528.539 -1528.539 
## (Intercept)  coeff:  -0.002081   tval:  -1.06422
## body_c       coeff:  0.0174682   tval:  19.54523
## lmer mass_c ~ wing_c + (1 | population) data_T2 FALSE 
## AIC:  -1449.581 -1449.581 
## (Intercept)  coeff:  -0.0022619  tval:  -1.067528
## wing_c       coeff:  0.0198835   tval:  14.79646

Sex * Host * Sym_Dist Modeling

data<-data.frame(R=data_T2$mass_c, 
                 A=data_T2$host_c, 
                 B=data_T2$sym_dist,
                 C=data_T2$sex_c,
                 X=data_T2$population)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R"))
##        [,1]      [,2]      [,3]      [,4]      [,5]       [,6]       [,7]      
## AICs   -1613.161 -1612.278 -1611.798 -1611.402 -1610.574  -1610.439  -1609.916 
## models 5         9         7         13        12         11         15        
## probs  0.261577  0.1682556 0.1323672 0.1085882 0.07176605 0.06708515 0.05165231
## 
## m5   R ~ A + C + (1 | X)
## m9   R ~ A * C + (1 | X)
## m7   R ~ A + B + C + (1 | X)
## m13  R ~ B * C + A + (1 | X)
## m12  R ~ A * C + B + (1 | X)
## m11  R ~ A * B + C + (1 | X)
## m15  R ~ A * B + B * C + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings)) 
## Number of models that failed to converge:  0
anova(m5, m9, test="Chisq") # Adding A*C does not improve fit
## Data: data
## Models:
## m5: R ~ A + C + (1 | X)
## m9: R ~ A * C + (1 | X)
##    Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m5  5 -1613.2 -1595.0 811.58  -1623.2                         
## m9  6 -1612.3 -1590.4 812.14  -1624.3 1.1175      1     0.2905
T2_model <- lmer(mass_c ~ host_c + sex_c + (1|population), data=data_T2, REML=FALSE)
## boundary (singular) fit: see ?isSingular
tidy_regression(T2_model, is_color=output_col)
## lmer mass_c ~ host_c + sex_c + (1 | population) data_T2 FALSE 
## AIC:  -1613.161 -1613.161 
## (Intercept)  coeff:  0.0051739   tval:  5.276994
## host_c       coeff:  -0.0026546  tval:  -2.800963
## sex_c        coeff:  0.0209841   tval:  24.58715
  • positive effect of sex, where if female then weigh more.
  • negative effect of sex, where if from GRT then weigh less.

Morphology Modeling

data<-data.frame(R=data_T2$mass_c, 
                 A=data_T2$beak_c, 
                 B=data_T2$thorax_c,
                 C=data_T2$body_c,
                 D=data_T2$wing_c,
                 X=data_T2$population)
data <- data %>%
  filter(!is.na(A))

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF REMLF.R"))
##        [,1]      [,2]       [,3]       [,4]       [,5]       [,6]      
## AICs   -1678.367 -1677.537  -1677.364  -1676.806  -1676.622  -1676.6   
## models 49        78         77         91         92         71        
## probs  0.149204  0.09851999 0.09036947 0.06835842 0.06235193 0.06167388
##        [,7]       [,8]      
## AICs   -1676.464  -1676.373 
## models 83         84        
## probs  0.05762915 0.05506535
## 
## m49  R ~ A * D + B * C + (1 | X)
## m78  R ~ A * C + A * D + B * D + (1 | X)
## m77  R ~ A * C + A * D + B * C + (1 | X)
## m91  R ~ A * B + A * C + A * D + B * C + (1 | X)
## m92  R ~ A * B + A * C + A * D + B * D + (1 | X)
## m71  R ~ A * B + A * D + B * C + (1 | X)
## m83  R ~ A * D + B * C + B * D + (1 | X)
## m84  R ~ A * D + B * C + C * D + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings))
## Number of models that failed to converge:  0
anova(m49, m77, test="Chisq") # Adding A*C does not improve fit
## Data: data
## Models:
## m49: R ~ A * D + B * C + (1 | X)
## m77: R ~ A * C + A * D + B * C + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m49  9 -1678.4 -1645.6 848.18  -1696.4                         
## m77 10 -1677.4 -1640.9 848.68  -1697.4 0.9972      1      0.318
T2_morph_model <- lmer(mass_c ~ beak_c * wing_c + thorax_c * body_c  + (1|population), data=data_T2, REML=FALSE)
## boundary (singular) fit: see ?isSingular
tidy_regression(T2_morph_model, is_color=output_col)
## lmer mass_c ~ beak_c * wing_c + thorax_c * body_c + (1 | population) data_T2 FALSE 
## AIC:  -1678.367 -1678.367 
## (Intercept)      coeff:  -0.0017783  tval:  -2.020758
## beak_c           coeff:  0.0118304   tval:  8.493899
## wing_c           coeff:  -0.0153314  tval:  -3.587623
## thorax_c         coeff:  0.0140166   tval:  2.031179
## body_c           coeff:  0.0163878   tval:  3.760543
## beak_c:wing_c    coeff:  -0.0079698  tval:  -5.114734
## thorax_c:body_c  coeff:  0.0203836   tval:  6.13718
data<-data.frame(R=data_T2$mass_c, 
                 A=data_T2$beak_c, 
                 B=data_T2$thorax_c,
                 C=data_T2$wing2body_c,
                 X=data_T2$population)
data <- data %>%
  filter(!is.na(A))

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R"))
##        [,1]      [,2]      [,3]      [,4]     
## AICs   -1648.129 -1647.388 -1646.457 -1645.854
## models 12        14        16        17       
## probs  0.363354  0.2509301 0.1575115 0.1165296
## 
## m12  R ~ A * C + B + (1 | X)
## m14  R ~ A * B + A * C + (1 | X)
## m16  R ~ A * C + B * C + (1 | X)
## m17  R ~ A * B + A * C + B * C + (1 | X)
cat("Number of models that failed to converge: ", length(errors$warnings))
## Number of models that failed to converge:  0
anova(m7, m12, test="Chisq") # Adding A*C does improve the fit
## Data: data
## Models:
## m7: R ~ A + B + C + (1 | X)
## m12: R ~ A * C + B + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m7   6 -1643.5 -1621.7 827.75  -1655.5                           
## m12  7 -1648.1 -1622.6 831.06  -1662.1 6.6295      1    0.01003 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m12, m16, test="Chisq") # Adding B*C does not improve fit
## Data: data
## Models:
## m12: R ~ A * C + B + (1 | X)
## m16: R ~ A * C + B * C + (1 | X)
##     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m12  7 -1648.1 -1622.6 831.06  -1662.1                         
## m16  8 -1646.5 -1617.3 831.23  -1662.5 0.3282      1     0.5667
T2_morph_model <- lmer(mass_c ~ beak_c*wing2body_c + thorax_c + (1|population), data=data_T2, REML=TRUE)
tidy_regression(T2_morph_model, is_color=output_col)
## lmer mass_c ~ beak_c * wing2body_c + thorax_c + (1 | population) data_T2 TRUE 
## AIC:  -1606.206 -1620.206 
## (Intercept)          coeff:  -6.5e-06    tval:  -0.00493547
## beak_c               coeff:  0.0106175   tval:  8.070563
## wing2body_c          coeff:  -0.1715218  tval:  -3.789612
## thorax_c             coeff:  0.0353188   tval:  8.349534
## beak_c:wing2body_c   coeff:  -0.1234253  tval:  -2.59209
  • positive effect of beak length, where the longer beak the heavier the bug
  • negative effect of wing2body ratio, where the larger the ratio, the lighter the bug.
  • positive effect of thorax, where the longer the thorax the heavier the bug
  • negative effect of beak*wing2body ratio where the longer the beak and ratio, the lighter the bug.

T1 vs. T2 Plotting

Histograms

#, fig.width=6, fig.height=2.5}
data_T1 <- filter(data_tested, trial_type == "T1")
data_T2 <- filter(data_tested, trial_type == "T2")
h1 <- as.grob(expression(hist(data_T1$mass, main="Trial 1", xlab= "Mass (g)")))
h2 <- as.grob(expression(hist(data_T2$mass, main="Trial 2", xlab = "Mass (g)")))
hist(data_all$mass, xlab="Mass (g)", main="Soapberry Bug Mass Histogram")

grid.arrange(h1, h2, ncol=2)

h <- gf_histogram(~mass, data=data_tested, xlab="Mass (g)", ylab="Count",
             title="Soapberry Bug Mass Histogram", fill=~trial_type) %>%
    gf_theme(theme = theme_minimal() )
h

# Uncomment below to save graph:
#ggsave("Mass-Hist-splitbyTrial", plot = h, device="jpeg",width = 5, height = 3)

Trial 1: Mass by Sex, Host Plant and Distance from Sympatric Zone

p1 <- as.grob(expression(
  plot(data_T1$sex_c, data_T1$mass, ylab="Mass (g)", xlab="Sex")))
p2 <- as.grob(expression(
  plot(data_T1$host_c, data_T1$mass, ylab="Mass (g)", xlab="Host")))
p3 <- as.grob(expression(
  plot(data_T1$sym_dist, data_T1$mass, ylab="Mass (g)", 
       xlab="Distance From Sympatric Zone"))) 
p.all <- gf_point(mass ~ sym_dist, data = data_T1,
                  color = ~host_c, ylab="Mass (g)", 
                  xlab="Distance From Sympatric Zone", 
                  alpha = ~sex_c) # GRT = blue, BV = black
p.all

g <- grid.arrange(p1,p2,p3, ncol=3, top= "Trial 1")

# Uncomment below to save graph:
# ggsave("Massbysex,host,symp_zone", plot = g, device="jpeg", width = 9, height = 4)
# If Female and From GRT and farther from HS then have more mass.

Trial 1: Mass by Morphology

p4 <- as.grob(expression(
  plot(data_T1$body, data_T1$mass, ylab="Mass (g)", xlab="Body (mm)")))
p5 <- as.grob(expression(
  plot(data_T1$beak, data_T1$mass, ylab="Mass (g)", xlab="Beak (mm)")))
p6 <- as.grob(expression(
  plot(data_T1$thorax, data_T1$mass, ylab="Mass (g)", xlab="Throax (mm)")))
p7 <- as.grob(expression(
  plot(data_T1$wing, data_T1$mass, ylab="Mass (g)", xlab="Wing (mm)")))

g2 <- grid.arrange(p4,p5,p6,p7, ncol=2, top="Trial 1")

# Uncomment below to save graph:
#ggsave("MassbyMorphologyT2", plot = g2, device="jpeg", width = 9, height = 7)

Trial 2: Mass by Morphology

p8 <- as.grob(expression(
  plot(data_T2$body, data_T2$mass, ylab="Mass (g)", xlab="Body (mm)")))
p9 <- as.grob(expression(
  plot(data_T2$beak, data_T2$mass, ylab="Mass (g)", xlab="Beak (mm)")))
p10 <- as.grob(expression(
  plot(data_T2$thorax, data_T2$mass, ylab="Mass (g)", xlab="Throax (mm)")))
p11 <- as.grob(expression(
  plot(data_T2$wing, data_T2$mass, ylab="Mass (g)", xlab="Wing (mm)")))

g3 <- grid.arrange(p8,p9,p10,p11, ncol=2, top="Trial 2")

# Uncomment below to save graph:
#ggsave("MassbyMorphologyT2", plot = g3, device="jpeg", width = 9, height = 7)