Modeling Mass

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

source_path = "~/Desktop/git_repositories/SBB-dispersal/avbernat/Rsrc/"

script_names = c("center_flight_data.R", # Re-centers data 
                 "clean_flight_data.R", # Loads and cleans data
                 "get_warnings.R", 
                 "compare_models.R",
                 "regression_output.R", # Cleans regression outputs; prints in color or B&W
                 "AICprobabilities.R",
                 "get_Akaike_weights.R",
                 "plotting-lm.R",
                 "plotting-lm2.R",
                 "plotting-polyreg.R")

for (script in script_names) { 
  path = paste0(source_path, script)
  source(path) 
}

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

## Loading required package: rstan

## Loading required package: StanHeaders

## rstan (Version 2.21.2, 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

## rethinking (Version 2.13)

## 
## Attaching package: 'rethinking'

## The following objects are masked from 'package:mosaic':
## 
##     compare, logit, resample

## The following object is masked from 'package:stats':
## 
##     rstudent

output_col = FALSE # Change to TRUE if working in Base R or RStudio; FALSE if generating an HTML

data <- read_flight_data("data/all_flight_data-Winter2020.csv")

## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))

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.059512    tval:  25.89193
## host_c                   coeff:  -0.0046734  tval:  -2.358484
## sex_c                    coeff:  0.024044    tval:  12.13203
## sym_dist                 coeff:  -0.0063639  tval:  -1.459846
## host_c:sex_c             coeff:  -0.0032142  tval:  -1.622066
## host_c:sym_dist          coeff:  0.0056227   tval:  1.289829
## sex_c:sym_dist           coeff:  -0.0107795  tval:  -2.472457
## host_c:sex_c:sym_dist    coeff:  0.0090188   tval:  2.068768

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.02715618
## chamberA-2   coeff:  -2.78e-05   tval:  -0.01242455
## chamberA-3   coeff:  0.0014132   tval:  0.6493085
## chamberA-4   coeff:  -0.0011209  tval:  -0.5387918
## chamberB-1   coeff:  0.0007711   tval:  0.3244349
## chamberB-2   coeff:  -0.0005543  tval:  -0.2597874
## chamberB-3   coeff:  0.0008228   tval:  0.3685082
## chamberB-4   coeff:  0.0013812   tval:  0.6343803

## 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.09777646
## min_from_IncStart_c  coeff:  2e-07   tval:  0.04539345

Biological Effects

## lmer mass_c ~ total_eggs + (1 | ID) + (1 | trial_type) data_tested FALSE 
## AIC:  -933.7505 -933.7505 
## (Intercept)  coeff:  0.016804    tval:  5.759778
## total_eggs   coeff:  0.0001144   tval:  5.728782

## 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.1036336
## beak_c       coeff:  0.0181976   tval:  22.65525

## 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.06098019
## thorax_c     coeff:  0.0596814   tval:  21.37376

## 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.1167303
## wing_c       coeff:  0.0185423   tval:  15.12272

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)

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

##        [,1]      [,2]      [,3]      [,4]      [,5]      [,6]       [,7]      
## AICs   -3494.901 -3493.398 -3493.123 -3492.966 -3492.913 -3491.609  -3491.415 
## models 13        16        8         4         15        10         17        
## probs  0.2881694 0.1359125 0.1184829 0.1094926 0.10664   0.05555309 0.05043256
## 
## 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)
##     npar     AIC     BIC logLik deviance  Chisq 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.9134

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.0049949   tval:  3.404279
## sex_c            coeff:  0.0212822   tval:  21.06871
## sym_dist         coeff:  -0.001048   tval:  -0.9495004
## host_c           coeff:  -0.0027384  tval:  -2.30813
## sex_c:sym_dist   coeff:  -0.002099   tval:  -2.4395

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.0051243   tval:  2.732483
## sex_c            coeff:  0.0212699   tval:  21.08207
## sym_dist         coeff:  -0.0009617  tval:  -0.8729017
## host_c           coeff:  -0.0027059  tval:  -2.283697
## sex_c:sym_dist   coeff:  -0.0020542  tval:  -2.392127

• 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))

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 4-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,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)
##      npar     AIC     BIC logLik deviance Chisq 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.745286
## beak_c           coeff:  0.0096144   tval:  7.243925
## thorax_c         coeff:  0.0135745   tval:  2.204248
## body_c           coeff:  0.0210139   tval:  5.252228
## wing_c           coeff:  -0.0201397  tval:  -5.061687
## beak_c:thorax_c  coeff:  -0.014353   tval:  -2.451462
## beak_c:body_c    coeff:  0.0094698   tval:  2.952663
## beak_c:wing_c    coeff:  -0.0161723  tval:  -4.629552
## thorax_c:body_c  coeff:  0.0237478   tval:  7.351266

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))

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,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)
##     npar     AIC     BIC logLik deviance  Chisq 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

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)
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}
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}
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

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.0020164   tval:  0.4077486
## chamberA-2   coeff:  -0.0024748  tval:  -0.4379859
## chamberA-3   coeff:  0.0019373   tval:  0.3217941
## chamberA-4   coeff:  -0.0048499  tval:  -0.8920994
## chamberB-1   coeff:  -0.0015048  tval:  -0.2414317
## chamberB-2   coeff:  -0.0070743  tval:  -1.192299
## chamberB-3   coeff:  -0.0033837  tval:  -0.6095594
## chamberB-4   coeff:  -0.0002181  tval:  -0.03897331

## 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.1610815
## min_from_IncStart_c  coeff:  5.5e-06 tval:  0.5800283

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.778524
## total_eggs   coeff:  0.0001081   tval:  5.462321

## 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.1243353
## beak_c       coeff:  0.0089663   tval:  3.845077

## 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.2930789
## wing_c       coeff:  0.0102625   tval:  4.746197

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)

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 2-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,1]      [,2]      [,3]      [,4]      [,5]     
## AICs   -1018.297 -1017.899 -1017.107 -1016.95  -1016.082
## models 2         1         3         4         5        
## probs  0.3113937 0.2552317 0.1717411 0.1587412 0.1028923
## 
## 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)
##    npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## m2    4 -1018.3 -1004.9 513.15  -1026.3                     
## m3    5 -1017.1 -1000.3 513.55  -1027.1 0.8099  1     0.3682

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.002563    tval:  1.00972
## sym_dist     coeff:  -0.0048643  tval:  -2.069924

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.0026618   tval:  0.894951
## sym_dist     coeff:  -0.0047002  tval:  -2.001503

• 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)

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,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.2128116 0.1683993 0.1186873 0.1012688 0.08791961 0.07907835 0.06342447
## 
## 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)
##     npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## m6     5 -1049.6 -1032.8 529.81  -1059.6                     
## m10    6 -1047.6 -1027.5 529.82  -1059.6 0.0201  1     0.8874

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.0087942  tval:  -2.928927
## sym_dist     coeff:  -0.0038776  tval:  -1.790074
## eggs_b       coeff:  0.0166019   tval:  6.007508

• 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'

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,1]      [,2]      [,3]      [,4]    
## AICs   -1018.297 -1017.899 -1017.107 -1016.95
## models 2         1         4         8       
## probs  0.346903  0.2843367 0.1913254 0.176843
## 
## 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)
##    npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## m2    4 -1018.3 -1004.9 513.15  -1026.3                     
## m4    5 -1017.1 -1000.3 513.55  -1027.1 0.8099  1     0.3682

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.002563    tval:  1.00972
## sym_dist     coeff:  -0.0048643  tval:  -2.069924

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.0026618   tval:  0.894951
## sym_dist     coeff:  -0.0047002  tval:  -2.001503

• 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))

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 4-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,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_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 4-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,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

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

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))

lm_plot(data_eggs, "total_eggs_c", "mass_c", "Total Eggs Laid", "Mass (g)", "Mass by Total Eggs Laid")

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}

lm_plot(data_fem, "wing2body_c", "mass_c", "Wing:Body", "Mass (g)", "Mass by Wing:Body")

lm_plot(data_male, "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")

lm_polyplot(data_fem, "wing2body_c", "mass_c", "Wing:Body", "Mass (g)", "Mass by Wing:Body")

lm_polyplot(data_male, "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

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.01018054
## chamberA-3   coeff:  -0.0009262  tval:  -0.8482126
## chamberA-4   coeff:  0.000511    tval:  0.4739869
## 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.8764587

## 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.05730676
## days_from_start_c    coeff:  8.06e-05    tval:  1.304994

## 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.09158187
## beak_c       coeff:  0.0095844   tval:  10.17196

## 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.07795958
## body_c       coeff:  0.0069936   tval:  15.40859

## 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)

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 2-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,1]      [,2]      [,3]     
## AICs   -2912.775 -2912.531 -2910.353
## models 4         3         1        
## probs  0.4574888 0.404919  0.136296 
## 
## 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)
##    npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## m3    5 -2912.5 -2892.6 1461.3  -2922.5                     
## m4    6 -2912.8 -2888.9 1462.4  -2924.8 2.2441  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.0024281  tval:  -2.817716
## host_c       coeff:  -0.0028508  tval:  -3.982546
## sym_dist     coeff:  0.0011125   tval:  2.052889

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.0022514  tval:  -1.792416
## host_c       coeff:  -0.0027559  tval:  -3.864916
## sym_dist     coeff:  0.0011121   tval:  2.061372

• 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))

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 4-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,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)
##     npar     AIC     BIC logLik deviance  Chisq 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)
##     npar     AIC     BIC logLik deviance  Chisq 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.7212858
## thorax_c         coeff:  0.0124757   tval:  4.682477
## wing_c           coeff:  -0.0026086  tval:  -1.525296
## body_c           coeff:  0.0054843   tval:  3.177489
## beak_c           coeff:  0.0039302   tval:  4.906169
## thorax_c:wing_c  coeff:  -0.0290863  tval:  -3.156054
## thorax_c:body_c  coeff:  0.0242504   tval:  3.652633
## wing_c:body_c    coeff:  0.0015995   tval:  1.474573

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}
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

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.2101054
## chamberA-2   coeff:  0.0063281   tval:  0.8408199
## chamberA-3   coeff:  -0.0043664  tval:  -0.59317
## chamberA-4   coeff:  -0.0018981  tval:  -0.2599025
## chamberB-1   coeff:  0.0011979   tval:  0.153893
## chamberB-2   coeff:  -0.0066758  tval:  -0.9129257
## chamberB-3   coeff:  -0.0001199  tval:  -0.01620094
## chamberB-4   coeff:  -0.0023396  tval:  -0.3043398

## lmer mass_c ~ days_from_start_c + (1 | chamber) data_T1 FALSE 
## AIC:  -1497.218 -1497.218 
## (Intercept)          coeff:  4.94e-05    tval:  0.03384751
## 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.07e-05    tval:  0.0285217
## min_from_IncStart_c  coeff:  6.3e-06 tval:  0.6418641

Morphology Effects

## lmer mass_c ~ beak_c + (1 | population) data_T1 FALSE 
## AIC:  -1736.24 -1736.24 
## (Intercept)  coeff:  0.0003952   tval:  0.3474807
## 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.896248
## body_c       coeff:  0.0164097   tval:  18.27008

## lmer mass_c ~ wing_c + (1 | population) data_T1 FALSE 
## AIC:  -1646.291 -1646.291 
## (Intercept)  coeff:  -0.001967   tval:  -0.9477481
## 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)

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,1]      [,2]      [,3]       [,4]      [,5]       [,6]      
## AICs   -1771.646 -1771.277 -1769.663  -1769.578 -1769.577  -1768.546 
## models 10        13        16         15        5          9         
## probs  0.2608932 0.2169678 0.09682692 0.0928092 0.09276491 0.05537529
## 
## 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)
##     npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## m10    6 -1771.7 -1748.9 891.82  -1783.7                     
## m13    7 -1771.3 -1744.7 892.64  -1785.3 1.6313  1     0.2015

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.646 -1771.646 
## (Intercept)      coeff:  0.0076014   tval:  6.536645
## sex_c            coeff:  0.020893    tval:  17.96643
## sym_dist         coeff:  -0.0027453  tval:  -2.854718
## sex_c:sym_dist   coeff:  -0.0023604  tval:  -2.454449

• 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))

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 4-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,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.0020693  tval:  -1.651817
## beak_c           coeff:  0.0088408   tval:  5.462607
## thorax_c         coeff:  0.0165342   tval:  2.248405
## body_c           coeff:  0.0220665   tval:  4.579658
## wing_c           coeff:  -0.0220187  tval:  -4.606201
## beak_c:thorax_c  coeff:  -0.0154746  tval:  -2.254761
## beak_c:body_c    coeff:  0.0094591   tval:  2.47639
## beak_c:wing_c    coeff:  -0.0168317  tval:  -4.023263
## thorax_c:body_c  coeff:  0.0247601   tval:  6.458338

• 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))

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

## 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)
##     npar     AIC     BIC logLik deviance  Chisq 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)
##     npar     AIC     BIC logLik deviance  Chisq 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

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.01455641
## days_from_start_c    coeff:  0.0002657   tval:  0.5213987

## 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)

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,1]      [,2]      [,3]      [,4]      [,5]       [,6]       [,7]     
## AICs   -1613.161 -1612.278 -1611.799 -1611.402 -1610.574  -1610.439  -1609.916
## models 5         9         7         13        12         11         15       
## probs  0.2615762 0.1682551 0.1323713 0.108583  0.07176833 0.06708562 0.0516501
## 
## 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)
##    npar     AIC     BIC logLik deviance  Chisq 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))

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 4-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,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)
##     npar     AIC     BIC logLik deviance  Chisq 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))

model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))

##        [,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)
##     npar     AIC     BIC logLik deviance  Chisq 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)
##     npar     AIC     BIC logLik deviance  Chisq 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)