Winter 2020 Flight Trials: Distance Flight 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 (glm) and mixed effect modeling (glmer) to analyze the flight results.

Problems Noticing

  • qqplots of models are not normal and most best fit models are the null model. Might need to retransform the the distance data.

All Data

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

script_names = c("center_flight_data.R", # Re-centers data 
                 "clean_flight_data.R", # Loads and cleans data
                 "unique_flight_data.R",
                 "get_warnings.R", 
                 "compare_models.R",
                 "regression_output.R", # Cleans regression outputs; prints in color or B&W
                 "AICprobabilities.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
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")
data_all <- data[[1]]
data_tested <- data[[2]]

Cleaning the Data

### Remove everyone who didn't fly (then remove distances = 0, if that didn't work fully)
data_flew <- data_tested[data_tested$flew_b == 1, ]
data_flew <- data_flew[data_flew$distance > 0, ]
data_flew <- center_data(data_flew)

### Break up by trial type
d_T1 <-data_flew[data_flew$trial_type=="T1",] 
d_T1 <- center_data(d_T1)
d_T2 <-data_flew[data_flew$trial_type=="T2",]
d_T2 <- center_data(d_T2)

Plots

h1 <-  as.grob(expression(
  hist(data_flew$distance)))
p1 <- as.grob(expression(
  plot(distance ~ mass_c, data=data_flew)))
p2 <- as.grob(expression(
  plot(distance ~ days_from_start_c, data=data_flew)))
p3 <- as.grob(expression(
  plot(distance ~ min_from_IncStart, data=data_flew)))
p4 <- as.grob(expression(
  plot(distance ~ wing_c, data=data_flew)))
p5 <- as.grob(expression(
  plot(distance ~ beak_c, data=data_flew)))
p6 <- as.grob(expression(
  plot(distance ~ thorax_c, data=data_flew)))
p7 <- as.grob(expression(
  plot(distance ~ body_c, data=data_flew)))
p8 <- as.grob(expression(
  plot(distance ~ wing2body_c, data=data_flew)))

grid.arrange(h1,p1,p2,p3,p4,p5,p6,p7,p8, ncol=3)

gf_point(distance~sym_dist, col=~host_plant, alpha=~sex_c, data=data_flew) 

summary <- aggregate(distance~sym_dist*sex*host_plant, data=data_flew, FUN=mean)
plot(summary$distance~summary$sym_dist, 
     col=c(1,2)[as.factor(summary$sex)], # Female = Black, Male = Red
     pch=c(19,22)[as.factor(summary$host_plant)], # Filled circle is GRT, Open square is BV
     main="Observed Data: distance ~ sex*host_plant*sym_dist",
     xlab = "Distance from Sympatric Zone (°)",
     ylab= "Distance Flew", 
     #sub=eq_glmer
     ) 
legend("topright",
       legend = c("F and K.elegans", "M and C.corindum", "F and C.corindum","M and K.elegans"),
       #inset=c(-0.27,0.2),
       col= 1:2,
       pch = c(0,16,19),
       title="Groups")

p <- ggplot(data_flew, aes(x=sym_dist, y=distance, color=host_plant)) + 
  geom_violin()
p + stat_summary(fun=mean, geom="point", shape=23, size=2)
## Warning: position_dodge requires non-overlapping x intervals

Trial 1

Testing link and covariates

Without Mass

With Mass

Min From Inc Start

Morphology

Testing link functions

model_test<-glm(distance~chamber, data=d_T1, family=Gamma(link="log")) #equivalent to using log link function in Gamma is to log transform distance, except it won't fuss about non-0 values
summary(model_test)
## 
## Call:
## glm(formula = distance ~ chamber, family = Gamma(link = "log"), 
##     data = d_T1)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -3.861  -2.851  -1.153   0.302   3.708  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.1002     0.5888  15.456  < 2e-16 ***
## chamberA-2   -0.9188     0.6835  -1.344  0.18049    
## chamberA-3   -1.5270     0.6626  -2.304  0.02228 *  
## chamberA-4   -2.1010     0.6493  -3.236  0.00143 ** 
## chamberB-1   -0.4164     0.7019  -0.593  0.55366    
## chamberB-2   -1.6692     0.6604  -2.528  0.01230 *  
## chamberB-3   -1.5418     0.6733  -2.290  0.02312 *  
## chamberB-4   -1.4030     0.6798  -2.064  0.04041 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 2.773147)
## 
##     Null deviance: 960.77  on 197  degrees of freedom
## Residual deviance: 890.43  on 190  degrees of freedom
## AIC: 3152.6
## 
## Number of Fisher Scoring iterations: 19
plot(model_test)

Testing experimental covariates

####### Effect of chamber A-4
d_T1$chamber <- relevel(d_T1$chamber, ref="A-2")
tidy_regression(glm(distance~chamber, data=d_T1, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ chamber Gamma(link = "log") d_T1 
## AIC:  3152.602 
## (Intercept)  coeff:  8.1813766   Pr(>|t|):  2.778239e-58 *
## chamberA-1   coeff:  0.9187927   Pr(>|t|):  0.1804913
## chamberA-3   coeff:  -0.6082159  Pr(>|t|):  0.1891488
## chamberA-4   coeff:  -1.1822543  Pr(>|t|):  0.008155016 *
## chamberB-1   coeff:  0.5023509   Pr(>|t|):  0.3317638
## chamberB-2   coeff:  -0.7503904  Pr(>|t|):  0.1032069
## chamberB-3   coeff:  -0.6230025  Pr(>|t|):  0.1928124
## chamberB-4   coeff:  -0.48418    Pr(>|t|):  0.3203151
####### No effect of test date
tidy_regression(glm(distance~days_from_start_c, data=d_T1, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ days_from_start_c Gamma(link = "log") d_T1 
## AIC:  3162.359 
## (Intercept)          coeff:  7.871885    Pr(>|t|):  4.737044e-138 *
## days_from_start_c    coeff:  0.0133802   Pr(>|t|):  0.655233
####### Effect of Minutes from IncStart time
tidy_regression(glm(distance~min_from_IncStart_c, data=d_T1, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ min_from_IncStart_c Gamma(link = "log") d_T1 
## AIC:  3158.304 
## (Intercept)          coeff:  7.8381821   Pr(>|t|):  6.743288e-136 *
## min_from_IncStart_c  coeff:  -0.00203    Pr(>|t|):  0.01733534 *

Without Mass

data<-data.frame(R=d_T1$distance, 
                 A=d_T1$host_c, 
                 B=d_T1$sex_c, 
                 C=d_T1$wing2body,
                 D=d_T1$mass_c,
                 X=d_T1$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 3-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      [,4]      [,5]       [,6]       [,7]     
## AICs   3166.985  3167.458  3167.509  3167.581  3168.401   3168.697   3168.993 
## models 2         18        4         1         6          12         7        
## probs  0.1480341 0.1169016 0.1139401 0.1099037 0.07292471 0.06291987 0.0542636
## 
## m2   R ~ B + (1 | X)
## m0   R ~ 1 + (1 | X)
## m4   R ~ A + B + (1 | X)
## m1   R ~ A + (1 | X)
## m6   R ~ B + C + (1 | X)
## m12  R ~ A * C + B + (1 | X)
## m7   R ~ A + B + C + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 9
anova(m0, m2, test="Chisq") # Adding B does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 3167.5 3177.3 -1580.7   3161.5                     
## m2    4 3167.0 3180.1 -1579.5   3159.0 2.4722  1     0.1159
anova(m0, m1, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 3167.5 3177.3 -1580.7   3161.5                     
## m1    4 3167.6 3180.7 -1579.8   3159.6 1.8765  1     0.1707
best.fit <- glm(distance ~ 1, family = Gamma(link = "log"), data = d_T1)
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ 1 Gamma(link = "log") d_T1 
## AIC:  3160.521 
## 1    coeff:  7.8732445   Pr(>|t|):  3.084544e-138 *
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

Null is the best and model not normal.

With Mass

d_T1 <- d_T1 %>%
  filter(!is.na(mass))
d_T1 <- center_data(d_T1)

data<-data.frame(R=d_T1$distance, 
                 A=d_T1$host_c, 
                 B=d_T1$sex_c, 
                 C=d_T1$mass_c, 
                 X=d_T1$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 3-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      [,4]       [,5]     
## AICs   3149.738  3149.947  3151.309  3151.59    3152.335 
## models 13        10        16        15         2        
## probs  0.2248449 0.2025211 0.1024962 0.08908694 0.0613667
## 
## m13  R ~ B * C + A + (1 | X)
## m10  R ~ B * C + (1 | X)
## m16  R ~ A * C + B * C + (1 | X)
## m15  R ~ A * B + B * C + (1 | X)
## m2   R ~ B + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 5
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 3149.9 3169.6 -1569.0   3137.9                     
## m13    7 3149.7 3172.7 -1567.9   3135.7 2.2091  1     0.1372
anova(m6, m10, test="Chisq") # Adding B*C does 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 3153.6 3170.0 -1571.8   3143.6                       
## m10    6 3149.9 3169.6 -1569.0   3137.9 5.6775  1    0.01718 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m2, m6, test="Chisq") # Adding C does not improve fit
## Data: data
## Models:
## m2: R ~ B + (1 | X)
## m6: R ~ B + C + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m2    4 3152.3 3165.5 -1572.2   3144.3                     
## m6    5 3153.6 3170.0 -1571.8   3143.6 0.7104  1     0.3993
best.fit <- glm(distance ~ sex_c * mass_c, family = Gamma(link = "log"), data = d_T1) # model failed to converge using glmer and (1|chamber) OR when using Gamma and standardize mass
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ sex_c * mass_c Gamma(link = "log") d_T1 
## AIC:  3145.635 
## (Intercept)      coeff:  8.4558379   Pr(>|t|):  1.118312e-92 *
## sex_c            coeff:  0.60635 Pr(>|t|):  0.00620706 *
## mass_c           coeff:  -14.0584879 Pr(>|t|):  0.2450473
## sex_c:mass_c     coeff:  -28.0921617 Pr(>|t|):  0.02084149 *
  • positive effect of sex where if F then more likely to disperse farther than M
  • no effect of mass
  • negative effect of sex*mass where if female and heavy then less likely to disperse far

Summary of model: Light females disperse farther than males. But the model is not normal:

s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

Min From Inc Start

# sym dist and host isn't really showing up before so left them out
# less models fail if use standardized variables and it fixes the scaling issues.
data<-data.frame(R=d_T1$distance, 
                 A=d_T1$min_from_IncStart_s, 
                 B=d_T1$sex_c, 
                 C=d_T1$mass_s, 
                 X=d_T1$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 3-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      [,4]      [,5]     
## AICs   3148.933  3149.271  3149.947  3150.241  3150.9   
## models 15        13        10        16        17       
## probs  0.2218554 0.1873472 0.1336146 0.1153273 0.0829801
## 
## m15  R ~ A * B + B * C + (1 | X)
## m13  R ~ B * C + A + (1 | X)
## m10  R ~ B * 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: 5
best.fit <- glmer(distance ~ min_from_IncStart_s * sex_c + sex_c * mass_s  + (1|chamber), family = Gamma(link = "log"), data = d_T1) # converges! 
tidy_regression(best.fit, is_color=output_col)
## glmer distance ~ min_from_IncStart_s * sex_c + sex_c * mass_s + (1 | chamber) d_T1 Gamma(link = "log") 
## AIC:  3148.933 3148.933 
## (Intercept)                  coeff:  8.658505    Pr(>|t|):  3.753127e-131 *
## min_from_IncStart_s          coeff:  -0.4161222  Pr(>|t|):  0.02019857 *
## sex_c                        coeff:  0.7552641   Pr(>|t|):  0.01404303 *
## mass_s                       coeff:  -0.3058487  Pr(>|t|):  0.2585672
## min_from_IncStart_s:sex_c    coeff:  -0.2694879  Pr(>|t|):  0.1144714
## sex_c:mass_s                 coeff:  -0.7933401  Pr(>|t|):  0.004200295 *

  • negative effect of min from start, where the later in the day the bug is tested, the less distance the bug will fly

  • positive effect from sex where if F then more likely to disperse far

  • no effect of mass

  • negative effect of min from start*sex, where if F and tested late in the day then less likely to disperse far

  • negative effect of sex*mass where the more mass and F then less likely to disperse far

s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

Morphology

d_T1 <- d_T1 %>%
  filter(!is.na(body))
d_T1 <- center_data(d_T1)

data<-data.frame(R=d_T1$distance,
                 A=d_T1$sex_c,
                 B=d_T1$thorax_c,
                 X=d_T1$chamber) # singular fit if run with min_from_IncStart meaning the variance is near 0

model_script = paste0(source_path,"generic models-Gamma glmer 2-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]     [,2]      [,3]      [,4]       [,5]      
## AICs   3149.004 3150.829  3152.094  3152.335   3153.046  
## models 2        3         4         1          5         
## probs  0.516395 0.2073736 0.1101474 0.09764641 0.06843761
## 
## m2   R ~ B + (1 | X)
## m3   R ~ A + B + (1 | X)
## m4   R ~ A * B + (1 | X)
## m1   R ~ A + (1 | X)
## m0   R ~ 1 + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 1
anova(m0, m2, test="Chisq") # adding B improves fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m0    3 3153.1 3162.9 -1573.5   3147.1                       
## m2    4 3149.0 3162.1 -1570.5   3141.0 6.0419  1    0.01397 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m0, m1, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m0    3 3153.1 3162.9 -1573.5   3147.1                       
## m1    4 3152.3 3165.5 -1572.2   3144.3 2.7109  1    0.09967 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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 3149.0 3162.1 -1570.5   3141.0                     
## m3    5 3150.8 3167.2 -1570.4   3140.8 0.1753  1     0.6754
best.fit <- glmer(distance ~ thorax_c + (1|chamber), family = Gamma(link = "log"), data = d_T1) 
tidy_regression(best.fit, is_color=output_col)
## glmer distance ~ thorax_c + (1 | chamber) d_T1 Gamma(link = "log") 
## AIC:  3149.004 3149.004 
## (Intercept)  coeff:  7.8194398   Pr(>|t|):  1.04196e-271 *
## thorax_c     coeff:  1.1856293   Pr(>|t|):  0.01900159 *
  • positive effect of thorax length where the larger the thorax the more likely the bug disperse farther
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

thorax_means<-aggregate(distance~thorax*sex, data=d_T1, FUN="mean")
d_T1$min_from_IncStart_b <- 0
d_T1$min_from_IncStart_b[d_T1$min_from_IncStart > 240] = 1

thorax_means<-aggregate(distance~thorax*min_from_IncStart_b*sex, data=d_T1, FUN="mean")
mass_meansT1<-aggregate(distance~mass*min_from_IncStart_b*sex, data=d_T1, FUN="mean")

Trial_1_Plots = function() {par(mfrow=c(1,2))
  plot(thorax_means$thorax, thorax_means$distance, 
       col=c(19,21)[as.factor(thorax_means$min_from_IncStart_b)], # green is 0 and 1 is blue
       pch=c(19,21)[as.factor(thorax_means$sex)], # open circles is male and closed is female
       xlab="thorax length (mm)",
       ylab="distance (m)") 
  mtext(expression(bold("A) Trial 1 | n = 197")), side=3, adj=0.01, line=0.5, cex=1.3)
  
  plot(mass_meansT1$mass, mass_meansT1$distance, 
       col=c(19,21)[as.factor(mass_meansT1$min_from_IncStart_b)], # green is 0 and 1 is blue
       pch=c(19,21)[as.factor(mass_meansT1$sex)], # open circles is male and closed is female
       xlab="mass (g)",
       ylab="distance (m)") 
  legend(0.097, 18000, 
         legend=c(expression(italic("Morning & F")), 
                  expression(italic("Afternoon & F")), 
                  expression(italic("Morning & M")), 
                  expression(italic("Afternoon & M"))), 
         pch=c(19,19,21,21), col=c(19,21,19,21), cex=1)
  mtext(expression(bold("B) Trial 1 | n = 197")), side=3, adj=0.01, line=0.5, cex=1.3)
}
Trial_1_Plots()

Trial 1 Summary: For females getting that morning start is important, for males not so much. Thorax matters because of the sex differences but wing2body does not matter for either sex this time around.

For mass, it seems like there isn’t a relationship for males but for females, if have extreme masses or gain more mass then don’t go far in distance.

Trial 2

Testing link and covariates

Without Mass

With Mass

Morphology

Testing link functions

model_test<-glm(distance~chamber, data=d_T2, family=Gamma(link="log")) 
## Warning: glm.fit: algorithm did not converge
summary(model_test)
## 
## Call:
## glm(formula = distance ~ chamber, family = Gamma(link = "log"), 
##     data = d_T2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.8630  -2.9530  -1.1910   0.1881   3.0812  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.659161   0.365132  20.976   <2e-16 ***
## chamberA-2   0.009787   0.530524   0.018    0.985    
## chamberA-3   0.607590   0.538674   1.128    0.261    
## chamberA-4  -0.731418   0.504502  -1.450    0.150    
## chamberB-1   0.754300   0.569017   1.326    0.187    
## chamberB-2   0.057136   0.516374   0.111    0.912    
## chamberB-3   0.339896   0.596258   0.570    0.570    
## chamberB-4  -0.567144   0.547698  -1.036    0.302    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 2.666424)
## 
##     Null deviance: 680.38  on 138  degrees of freedom
## Residual deviance: 649.25  on 131  degrees of freedom
## AIC: 2193.3
## 
## Number of Fisher Scoring iterations: 25
plot(model_test)

Testing experimental Covariates

####### No effect of chamber
tidy_regression(glm(distance~chamber, data=d_T2, family=Gamma(link="log")), is_color=output_col)
## Warning: glm.fit: algorithm did not converge
## glm distance ~ chamber Gamma(link = "log") d_T2 
## AIC:  2193.29 
## (Intercept)  coeff:  7.6591608   Pr(>|t|):  1.045132e-43 *
## chamberA-2   coeff:  0.0097871   Pr(>|t|):  0.9853096
## chamberA-3   coeff:  0.6075901   Pr(>|t|):  0.2614083
## chamberA-4   coeff:  -0.7314184  Pr(>|t|):  0.1495089
## chamberB-1   coeff:  0.7543003   Pr(>|t|):  0.1872725
## chamberB-2   coeff:  0.0571363   Pr(>|t|):  0.912064
## chamberB-3   coeff:  0.3398956   Pr(>|t|):  0.5696213
## chamberB-4   coeff:  -0.5671436  Pr(>|t|):  0.3023408
####### No effect of test date
tidy_regression(glm(distance~days_from_start, data=d_T2, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ days_from_start Gamma(link = "log") d_T2 
## AIC:  2190.528 
## (Intercept)      coeff:  7.2989157   Pr(>|t|):  3.40326e-18 *
## days_from_start  coeff:  0.0340738   Pr(>|t|):  0.5061641
####### No effect of Minutes from IncStart time
tidy_regression(glm(distance~min_from_IncStart, data=d_T2, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ min_from_IncStart Gamma(link = "log") d_T2 
## AIC:  2189.8 
## (Intercept)          coeff:  8.0188723   Pr(>|t|):  8.902866e-67 *
## min_from_IncStart    coeff:  -0.0012527  Pr(>|t|):  0.2056701

Without Mass

data<-data.frame(R=d_T2$distance, 
                 A=d_T2$host_c, 
                 B=d_T2$sex_c, 
                 C=d_T2$sym_dist, 
                 X=d_T2$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 3-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      [,4]      [,5]       [,6]      
## AICs   2192.804  2194.723  2194.772  2194.797  2196.668   2196.691  
## models 18        1         2         3         5          4         
## probs  0.3505176 0.1342952 0.1310719 0.1294131 0.05077559 0.05020805
## 
## m0   R ~ 1 + (1 | X)
## m1   R ~ A + (1 | X)
## m2   R ~ B + (1 | X)
## m3   R ~ C + (1 | X)
## m5   R ~ A + C + (1 | X)
## m4   R ~ A + B + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 0
anova(m0, m1, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 2192.8 2201.6 -1093.4   2186.8                     
## m1    4 2194.7 2206.5 -1093.4   2186.7 0.0813  1     0.7756
anova(m1, m4, test="Chisq") # Adding B does not improve fit
## Data: data
## Models:
## m1: R ~ A + (1 | X)
## m4: R ~ A + B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m1    4 2194.7 2206.5 -1093.4   2186.7                     
## m4    5 2196.7 2211.4 -1093.3   2186.7 0.0323  1     0.8574
anova(m1, m5, test="Chisq") # Adding C does not improve fit
## Data: data
## Models:
## m1: R ~ A + (1 | X)
## m5: R ~ A + C + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m1    4 2194.7 2206.5 -1093.4   2186.7                     
## m5    5 2196.7 2211.3 -1093.3   2186.7 0.0548  1      0.815
best.fit <- glm(distance ~ 1, family = Gamma(link = "log"), data = d_T2) 
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ 1 Gamma(link = "log") d_T2 
## AIC:  2188.804 
## 1    coeff:  7.776082    Pr(>|t|):  2.869045e-96 *
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

With Mass

d_T2 <- d_T2 %>%
  filter(!is.na(mass))
d_T2 <- center_data(d_T2)

data<-data.frame(R=d_T2$distance, 
                 A=d_T2$host_c, 
                 B=d_T2$sex_c, 
                 C=d_T2$wing2body_c, # if use sym_dist the null model is the best model
                 D=d_T2$mass_c, X=d_T2$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 4-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]       [,2]      [,3]       [,4]     
## AICs   2192.213   2192.315  2192.804   2193.355 
## models 104        112       113        86       
## probs  0.09703104 0.0922427 0.07220904 0.0548264
## 
## m104     R ~ A * C + B * C + B * D + C * D + (1 | X)
## m112     R ~ A * B + A * C + A * D + B * C + B * D + C * D + (1 | X)
## m0   R ~ 1 + (1 | X)
## m86  R ~ 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: 41
best.fit <- glm(distance ~ host_c*wing2body_c + sex_c*wing2body_c + sex_c*mass_c + wing2body*mass_c, 
                family = Gamma(link = "log"), data = d_T2) 
## Warning: glm.fit: algorithm did not converge
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ host_c * wing2body_c + sex_c * wing2body_c + sex_c * mass_c + wing2body * mass_c Gamma(link = "log") d_T2 
## AIC:  2188.213 
## (Intercept)          coeff:  8.7901879   Pr(>|t|):  5.841391e-47 *
## host_c               coeff:  0.1109322   Pr(>|t|):  0.5268982
## wing2body_c          coeff:  49.6252618  Pr(>|t|):  0.01985163 *
## sex_c                coeff:  0.5696769   Pr(>|t|):  0.1246919
## mass_c               coeff:  3844.4895006    Pr(>|t|):  2.289974e-07 *
## host_c:wing2body_c   coeff:  -30.7301024 Pr(>|t|):  0.004474874 *
## wing2body_c:sex_c    coeff:  91.269555   Pr(>|t|):  0.0004919267 *
## sex_c:mass_c         coeff:  -56.7178313 Pr(>|t|):  0.0001810646 *
## mass_c:wing2body     coeff:  -5308.9820089   Pr(>|t|):  2.482135e-07 *
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

Min From Start

d_T2 <- d_T2 %>%
  filter(!is.na(mass))
d_T2 <- center_data(d_T2)

data<-data.frame(R=d_T2$distance, 
                 A=d_T2$min_from_IncStart, 
                 B=d_T2$sex_c, 
                 C=d_T2$wing2body_c, # if use sym_dist the null model is the best model
                 D=d_T2$mass_c, 
                 X=d_T2$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 4-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      
## AICs   2192.804  2193.355  2193.8    
## models 113       86        1         
## probs  0.0931372 0.0707166 0.05662476
## 
## m0   R ~ 1 + (1 | X)
## m86  R ~ B * C + B * D + C * D + (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: 97
anova(m0, m1, test="Chisq") # Adding A does not improve fit - null model the best fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 2192.8 2201.6 -1093.4   2186.8                     
## m1    4 2193.8 2205.5 -1092.9   2185.8 1.0047  1     0.3162

Morphology

d_T2 <- d_T2 %>%
  filter(!is.na(thorax))
d_T2 <- center_data(d_T2)

data<-data.frame(R=d_T2$distance,
                 A=d_T2$thorax_c, 
                 B=d_T2$wing2body_c, 
                 X=d_T2$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 2-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      [,4]      [,5]      
## AICs   2192.804  2193.92   2194.462  2195.203  2196.308  
## models 5         1         2         3         4         
## probs  0.4025872 0.2304731 0.1757752 0.1213282 0.06983625
## 
## m0   R ~ 1 + (1 | X)
## m1   R ~ A + (1 | X)
## m2   R ~ B + (1 | X)
## m3   R ~ A + B + (1 | X)
## m4   R ~ A * B + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 0
anova(m0, m2, test="Chisq") # Adding B does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 2192.8 2201.6 -1093.4   2186.8                     
## m2    4 2194.5 2206.2 -1093.2   2186.5 0.3426  1     0.5583
anova(m0, m3, test="Chisq") # Adding C does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m3: R ~ A + B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 2192.8 2201.6 -1093.4   2186.8                     
## m3    5 2195.2 2209.9 -1092.6   2185.2 1.6012  2     0.4491
best.fit <- glm(distance ~ 1, family = Gamma(link = "log"), data = d_T2) 
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ 1 Gamma(link = "log") d_T2 
## AIC:  2188.804 
## 1    coeff:  7.776082    Pr(>|t|):  2.869045e-96 *
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

d_T2$min_from_IncStart_b <- 0
d_T2$min_from_IncStart_b[d_T2$min_from_IncStart > 240] = 1

wing_means<-aggregate(distance~wing2body*sex*min_from_IncStart_b*mass, data=d_T2, FUN="mean")
mass_means<-aggregate(distance~min_from_IncStart_b*sex*mass, data=d_T2, FUN="mean")

Trial_2_Plots = function() {
  par(mfrow=c(1,2))
  plot(wing_means$wing2body, wing_means$distance, 
       col=c(19,21)[as.factor(wing_means$min_from_IncStart_b)], # green is 0 and 1 is blue
       pch=c(19,21)[as.factor(wing_means$sex)], # open circles is male and closed is female
       xlab="wing-to-body ratio",
       ylab="distance (m)") 
  mtext(expression(bold("A) Trial 2 | n = 139")), side=3, adj=0.01, line=0.5, cex=1.3)
  
  plot(mass_means$mass, mass_means$distance,
       col=c(19,21)[as.factor(mass_means$min_from_IncStart_b)], # green is 0 and 1 is blue
       pch=c(19,21)[as.factor(mass_means$sex)], # open circles is male and closed is female
       xlab="mass (g)",
       ylab="distance (m)")
  legend(0.097, 9000, 
         legend=c(expression(italic("Morning & F")), 
                  expression(italic("Afternoon & F")), 
                  expression(italic("Morning & M")), 
                  expression(italic("Afternoon & M"))), 
         pch=c(19,19,21,21), col=c(19,21,19,21), cex=1)
  mtext(expression(bold("B) Trial 1 | n = 197")), side=3, adj=0.01, line=0.5, cex=1.3)
}
Trial_2_Plots()

A) Males have larger wing to body ratios but it doesn’t necessarily determine distance (large spread). More important is whether the males flew in the morning (green). For females no factor seems to stand out as influences its distance flown.

  1. Except for some outliers it seems like females that weigh more are less likely to fly whereas with males it’s not affecting them.

Females

Testing Link Functions

Without Mass

With Mass

Minutes from IncStart

Morphology

Eggs

Testing link functions

d_fem = data_flew[data_flew$sex =="F",]
model_test<-glm(distance~chamber, data=d_fem, family=Gamma(link="log")) 
## Warning: glm.fit: algorithm did not converge
summary(model_test)
## 
## Call:
## glm(formula = distance ~ chamber, family = Gamma(link = "log"), 
##     data = d_fem)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.0731  -2.9133  -1.3246   0.3074   2.7411  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.0202     0.8218   4.892 7.24e-06 ***
## chamberA-2    4.3581     0.9489   4.593 2.15e-05 ***
## chamberA-3    3.7835     0.9637   3.926 0.000216 ***
## chamberA-4    3.4060     0.9117   3.736 0.000405 ***
## chamberB-1    4.9465     0.9489   5.213 2.19e-06 ***
## chamberB-2   10.4795     0.9637  10.875 4.29e-16 ***
## chamberB-3    4.8171     0.9637   4.999 4.88e-06 ***
## chamberB-4    3.2252     0.9117   3.538 0.000764 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 2.026096)
## 
##     Null deviance: 342.61  on 70  degrees of freedom
## Residual deviance: 393.35  on 63  degrees of freedom
## AIC: 1189.4
## 
## Number of Fisher Scoring iterations: 25
plot(model_test)

Testing covariates

####### Effect of chamber A-1, B-2, and marginal effect of B-4 (also algorithm did not converge)
d_fem$chamber <- relevel(d_fem$chamber, ref="A-2")
tidy_regression(glm(distance~chamber, data=d_fem, family=Gamma(link="log")), is_color=output_col)
## Warning: glm.fit: algorithm did not converge
## glm distance ~ chamber Gamma(link = "log") d_fem 
## AIC:  1189.373 
## (Intercept)  coeff:  8.3783004   Pr(>|t|):  4.397746e-26 *
## chamberA-1   coeff:  -4.3581407  Pr(>|t|):  2.148564e-05 *
## chamberA-3   coeff:  -0.5745949  Pr(>|t|):  0.4092489
## chamberA-4   coeff:  -0.9521374  Pr(>|t|):  0.1279382
## chamberB-1   coeff:  0.5883994   Pr(>|t|):  0.3838743
## chamberB-2   coeff:  6.1213245   Pr(>|t|):  1.17643e-12 *
## chamberB-3   coeff:  0.4589939   Pr(>|t|):  0.5093556
## chamberB-4   coeff:  -1.1329356  Pr(>|t|):  0.07115085 .
####### No effect of test date
tidy_regression(glm(distance~days_from_start, data=d_fem, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ days_from_start Gamma(link = "log") d_fem 
## AIC:  1161.236 
## (Intercept)      coeff:  8.3841161   Pr(>|t|):  3.165135e-37 *
## days_from_start  coeff:  -0.0369398  Pr(>|t|):  0.1576866
####### Marginal effect of Minutes from IncStart time
tidy_regression(glm(distance~min_from_IncStart, data=d_fem, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ min_from_IncStart Gamma(link = "log") d_fem 
## AIC:  1161.236 
## (Intercept)          coeff:  8.4039362   Pr(>|t|):  1.15342e-39 *
## min_from_IncStart    coeff:  -0.002202   Pr(>|t|):  0.09089658 .
#### No effect of eggs laid
tidy_regression(glm(distance~eggs_b, data=d_fem, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ eggs_b Gamma(link = "log") d_fem 
## AIC:  1162.838 
## (Intercept)  coeff:  8.0161517   Pr(>|t|):  1.301713e-46 *
## eggs_b       coeff:  0.0743005   Pr(>|t|):  0.8523971
### No effect of total eggs laid
tidy_regression(glm(distance~total_eggs, data=d_fem, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ total_eggs Gamma(link = "log") d_fem 
## AIC:  882.7586 
## (Intercept)  coeff:  8.2965485   Pr(>|t|):  2.359389e-36 *
## total_eggs   coeff:  -0.0020423  Pr(>|t|):  0.4713187

Without Mass

data<-data.frame(R=d_fem$distance,
                 A=d_fem$host_c, 
                 B=d_fem$sym_dist, 
                 X=d_fem$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 2-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      [,4]      
## AICs   1164.863  1166.824  1166.831  1168.655  
## models 5         1         2         3         
## probs  0.5115039 0.1918707 0.1912365 0.07682723
## 
## m0   R ~ 1 + (1 | X)
## m1   R ~ A + (1 | X)
## m2   R ~ B + (1 | X)
## m3   R ~ A + B + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 0
anova(m0, m2, test="Chisq") # Adding B does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 1164.9 1171.7 -579.43   1158.9                     
## m2    4 1166.8 1175.9 -579.42   1158.8 0.0323  1     0.8573
anova(m0, m1, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 1164.9 1171.7 -579.43   1158.9                     
## m1    4 1166.8 1175.9 -579.41   1158.8 0.0389  1     0.8436
anova(m0, m3, test="Chisq") # Adding C does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m3: R ~ A + B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 1164.9 1171.7 -579.43   1158.9                     
## m3    5 1168.7 1180.0 -579.33   1158.7 0.2084  2      0.901
best.fit <- glm(distance ~ 1, family = Gamma(link = "log"), data = d_fem) 
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ 1 Gamma(link = "log") d_fem 
## AIC:  1160.863 
## 1    coeff:  8.0397701   Pr(>|t|):  7.154653e-53 *
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

With Mass

d_fem <- d_fem %>%
  filter(!is.na(mass))
d_fem <- center_data(d_fem)

data<-data.frame(R=d_fem$distance, 
                 A=d_fem$host_c, 
                 B=d_fem$wing2body_c, 
                 C=d_fem$mass_c, 
                 X=d_fem$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 2-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      [,4]     
## AICs   1150.318  1151.516  1152.262  1153.261 
## models 5         2         1         3        
## probs  0.4458972 0.2450045 0.1687157 0.1023804
## 
## m0   R ~ 1 + (1 | X)
## m2   R ~ B + (1 | X)
## m1   R ~ A + (1 | X)
## m3   R ~ A + B + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 0
anova(m0, m2, test="Chisq") # Adding B does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 1150.3 1157.1 -572.16   1144.3                     
## m2    4 1151.5 1160.5 -571.76   1143.5 0.8024  1     0.3704
anova(m0, m1, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 1150.3 1157.1 -572.16   1144.3                     
## m1    4 1152.3 1161.3 -572.13   1144.3 0.0563  1     0.8125
best.fit <- glm(distance ~ 1, family = Gamma(link = "log"), data = d_fem) 
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ 1 Gamma(link = "log") d_fem 
## AIC:  1146.318 
## 1    coeff:  8.0534097   Pr(>|t|):  1.84817e-52 *
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

Min From Inc Start

data<-data.frame(R=d_fem$distance, 
                 A=d_fem$min_from_IncStart_s, 
                 B=d_fem$wing2body_c, # eggs_b doesn't converge
                 C=d_fem$mass_s, 
                 X=d_fem$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 3-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]     [,2]      [,3]      [,4]       [,5]       [,6]       [,7]      
## AICs   1150.318 1151.042  1151.077  1151.206   1151.516   1151.584   1151.931  
## models 18       1         3         13         2          10         5         
## probs  0.151147 0.1052397 0.1034211 0.09695245 0.08304985 0.08025896 0.06748582
##        [,8]       [,9]     
## AICs   1152.346   1152.396 
## models 6          4        
## probs  0.05484872 0.0534821
## 
## m0   R ~ 1 + (1 | X)
## m1   R ~ A + (1 | X)
## m3   R ~ C + (1 | X)
## m13  R ~ B * C + A + (1 | X)
## m2   R ~ B + (1 | X)
## m10  R ~ B * C + (1 | X)
## m5   R ~ A + C + (1 | X)
## m6   R ~ B + C + (1 | X)
## m4   R ~ A + B + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 0
anova(m0, m1, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## m0    3 1150.3 1157.1 -572.16   1144.3                    
## m1    4 1151.0 1160.0 -571.52   1143.0 1.276  1     0.2586
anova(m0, m3, test="Chisq") # Adding C does not impove fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m3: R ~ C + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 1150.3 1157.1 -572.16   1144.3                     
## m3    4 1151.1 1160.1 -571.54   1143.1 1.2411  1     0.2653
anova(m0, m2, test="Chisq") # Adding B does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 1150.3 1157.1 -572.16   1144.3                     
## m2    4 1151.5 1160.5 -571.76   1143.5 0.8024  1     0.3704
best.fit <- glm(distance ~ 1, family = Gamma(link = "log"), data = d_fem) 
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ 1 Gamma(link = "log") d_fem 
## AIC:  1146.318 
## 1    coeff:  8.0534097   Pr(>|t|):  1.84817e-52 *
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

Morphology

d_fem <- d_fem %>%
  filter(!is.na(mass))
d_fem <- center_data(d_fem)

data<-data.frame(R=d_fem$distance, 
                 A=d_fem$wing2body, 
                 B=d_fem$thorax_c,
                 X=d_fem$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 2-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]     
## AICs   1141.397 
## models 4        
## probs  0.9421605
## 
## m4   R ~ A * B + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 0
best.fit <- glmer(distance ~ wing2body*thorax_c + (1|chamber), family = Gamma(link = "log"), data = d_fem) 
## boundary (singular) fit: see ?isSingular
tidy_regression(best.fit, is_color=output_col) # the coeffients are crazy large
## glmer distance ~ wing2body * thorax_c + (1 | chamber) d_fem Gamma(link = "log") 
## AIC:  1141.397 1141.397 
## (Intercept)          coeff:  -4.8717952  Pr(>|t|):  0.7516128
## wing2body            coeff:  18.1205325  Pr(>|t|):  0.3997316
## thorax_c             coeff:  245.3248807 Pr(>|t|):  0.0004739619 *
## wing2body:thorax_c   coeff:  -337.9088657    Pr(>|t|):  0.0005213684 *
  • no effect of wing2body
  • positive effect of thorax where the wider the thorax the more likely the bug flies far
  • negative effect of wing2body*thorax where the wider the thorax and longer the wings in comparison to the body, the less likely the bugs fly far (contradicting)
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

Eggs

d_fem <- d_fem %>%
  filter(!is.na(mass))
d_fem <- center_data(d_fem)

data<-data.frame(R=d_fem$distance, 
                 A=d_fem$eggs_b, 
                 B=d_fem$mass_c, 
                 C=d_fem$thorax_c, 
                 X=d_fem$chamber) # adding D=d_fem$wing2body_c leads to a massive converging error

model_script = paste0(source_path,"generic models-Gamma glmer 3-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      [,4]     
## AICs   1138.718  1139.228  1139.731  1140.537 
## models 16        13        15        17       
## probs  0.3433065 0.2660775 0.2068931 0.1382829
## 
## m16  R ~ A * C + B * C + (1 | X)
## m13  R ~ B * C + A + (1 | X)
## m15  R ~ A * B + 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, m16, test="Chisq") # Adding A*C does not improve fit
## Data: data
## Models:
## m13: R ~ B * C + A + (1 | X)
## m16: R ~ A * C + B * C + (1 | X)
##     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m13    7 1139.2 1155.0 -562.61   1125.2                     
## m16    8 1138.7 1156.7 -561.36   1122.7 2.5097  1     0.1131
anova(m10, m13, test="Chisq") # Adding A improves 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 1143.9 1157.4 -565.95   1131.9                        
## m13    7 1139.2 1155.0 -562.61   1125.2 6.6709  1     0.0098 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
best.fit <- glmer(distance ~  mass_c * thorax_c + eggs_b +  (1|chamber), family = Gamma(link = "log"), data = d_fem) 
## boundary (singular) fit: see ?isSingular
tidy_regression(best.fit, is_color=output_col) # the coeffients are crazy large
## glmer distance ~ mass_c * thorax_c + eggs_b + (1 | chamber) d_fem Gamma(link = "log") 
## AIC:  1139.228 1139.228 
## (Intercept)      coeff:  7.6876454   Pr(>|t|):  1.667819e-170 *
## mass_c           coeff:  -37.7315159 Pr(>|t|):  0.006364513 *
## thorax_c         coeff:  4.8026494   Pr(>|t|):  7.085296e-05 *
## eggs_b           coeff:  1.4719465   Pr(>|t|):  0.01622924 *
## mass_c:thorax_c  coeff:  -229.6775316    Pr(>|t|):  0.000277704 *
  • negative effect of mass where the larger the mass the less likely the bug will disperse far
  • positive effect of thorax where the wider the thorax the more likely the bug will disperse far
  • positive effect of eggs laid on test day where if female laid eggs then more likely to disperse far
  • negative effect of mass*thorax where if heavy and have a wide thorax then less likley to disperse far (conflicting)
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

d_fem$min_from_IncStart_b <- 0
d_fem$min_from_IncStart_b[d_fem$min_from_IncStart > 240] = 1

means<-aggregate(distance~mass*thorax*eggs_b*min_from_IncStart_b, data=d_fem, FUN="mean")

Fem_Plots = function() {par(mfrow=c(1,2))
  plot(means$thorax, means$distance, 
       col=c(19,21)[as.factor(means$min_from_IncStart_b)], # green is 0 and 1 is blue
       pch=c(19,21)[as.factor(means$eggs_b)], # open circles is male and closed is female
       xlab="thorax length (mm)",
       ylab="distance (m)") 
  mtext(expression(bold("A) Females | n = 70")), side=3, adj=0.01, line=0.5, cex=1.3)
  
  plot(means$mass, means$distance, 
       col=c(19,21)[as.factor(means$min_from_IncStart_b)], # green is 0 and 1 is blue
       pch=c(19,21)[as.factor(means$eggs_b)], # solid is 0 and open is 1
       xlab="mass (g)",
       ylab="distance (m)") 
  legend(0.103, 18000, 
         legend=c(expression(italic("Morning & No Eggs")), 
                  expression(italic("Afternoon & No Eggs")), 
                  expression(italic("Morning & Eggs")), 
                  expression(italic("Afternoon & Eggs"))), 
         pch=c(19,19,21,21), col=c(19,21,19,21), cex=0.85)
  mtext(expression(bold("B) Females | n = 70")), side=3, adj=0.01, line=0.5, cex=1.3)
}
Fem_Plots()

Females who started flying in the morning were especially able to hit 15km and most who flew had low mass and no eggs and a larger thorax. Those with extreme masses, small thoraxes, and flew in the afternoon had less of a chance of lying far.

Some factors: eggs, mass, thorax, wing2body

Males

Testing Link Functions

Without Mass

With Mass

Minutes from IncStart

Morphology

Testing link functions

d_male = data_flew[data_flew$sex == "M",]
model_test<-glm(distance~chamber, data=d_male, family=Gamma(link="log"))
summary(model_test)
## 
## Call:
## glm(formula = distance ~ chamber, family = Gamma(link = "log"), 
##     data = d_male)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.7579  -2.8903  -1.2139   0.1589   3.6639  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.42409    0.33751  24.960  < 2e-16 ***
## chamberA-2  -0.58064    0.45045  -1.289 0.198544    
## chamberA-3  -0.52660    0.43235  -1.218 0.224344    
## chamberA-4  -1.62769    0.41931  -3.882 0.000132 ***
## chamberB-1  -0.04254    0.48225  -0.088 0.929775    
## chamberB-2  -0.95771    0.42443  -2.256 0.024877 *  
## chamberB-3  -0.89942    0.45699  -1.968 0.050121 .  
## chamberB-4  -0.82562    0.46838  -1.763 0.079136 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 2.847771)
## 
##     Null deviance: 1295.1  on 265  degrees of freedom
## Residual deviance: 1231.3  on 258  degrees of freedom
## AIC: 4181.8
## 
## Number of Fisher Scoring iterations: 19
plot(model_test)

Testing covariates

####### Effect of chamber A-4
d_male$chamber <- relevel(d_male$chamber, ref="A-2")
tidy_regression(glm(distance~chamber, data=d_male, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ chamber Gamma(link = "log") d_male 
## AIC:  4181.812 
## (Intercept)  coeff:  7.8434462   Pr(>|t|):  6.006056e-75 *
## chamberA-1   coeff:  0.580639    Pr(>|t|):  0.1985441
## chamberA-3   coeff:  0.0540395   Pr(>|t|):  0.893304
## chamberA-4   coeff:  -1.0470501  Pr(>|t|):  0.007492747 *
## chamberB-1   coeff:  0.5380972   Pr(>|t|):  0.2387497
## chamberB-2   coeff:  -0.3770749  Pr(>|t|):  0.3394169
## chamberB-3   coeff:  -0.3187812  Pr(>|t|):  0.4579601
## chamberB-4   coeff:  -0.2449808  Pr(>|t|):  0.5790111
####### No effect of test date
tidy_regression(glm(distance~days_from_start, data=d_male, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ days_from_start Gamma(link = "log") d_male 
## AIC:  4188.748 
## (Intercept)      coeff:  7.6194576   Pr(>|t|):  9.885506e-115 *
## days_from_start  coeff:  0.0137231   Pr(>|t|):  0.3534092
####### Effect of Minutes from IncStart time
tidy_regression(glm(distance~min_from_IncStart, data=d_male, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ min_from_IncStart Gamma(link = "log") d_male 
## AIC:  4185.843 
## (Intercept)          coeff:  8.0396056   Pr(>|t|):  1.062798e-132 *
## min_from_IncStart    coeff:  -0.0016051  Pr(>|t|):  0.02797034 *
#### No effect of eggs laid
tidy_regression(glm(distance~eggs_b, data=d_male, family=Gamma(link="log")), is_color=output_col)
## glm distance ~ eggs_b Gamma(link = "log") d_male 
## AIC:  4187.417 
## 1    coeff:  7.7714908   Pr(>|t|):  1.699765e-181 *

Without Mass

data<-data.frame(R=d_male$distance, 
                 A=d_male$host_c, 
                 B=d_male$sym_dist, 
                 X=d_male$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 2-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]     [,4]     
## AICs   4195.531  4197.118  4197.474 4198.595 
## models 5         1         2        3        
## probs  0.4686086 0.2118835 0.177363 0.1012803
## 
## m0   R ~ 1 + (1 | X)
## m1   R ~ A + (1 | X)
## m2   R ~ B + (1 | X)
## m3   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(m0, m2, test="Chisq") # Adding B does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 4195.5 4206.3 -2094.8   4189.5                     
## m2    4 4197.5 4211.8 -2094.7   4189.5 0.0569  1     0.8115
anova(m0, m1, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 4195.5 4206.3 -2094.8   4189.5                     
## m1    4 4197.1 4211.5 -2094.6   4189.1 0.4125  1     0.5207
anova(m0, m3, test="Chisq") # Adding C does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m3: R ~ A + B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 4195.5 4206.3 -2094.8   4189.5                     
## m3    5 4198.6 4216.5 -2094.3   4188.6 0.9362  2     0.6262
best.fit <- glm(distance ~ 1, family = Gamma(link = "log"), data = d_male) 
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ 1 Gamma(link = "log") d_male 
## AIC:  4187.417 
## 1    coeff:  7.7714908   Pr(>|t|):  1.699765e-181 *
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

With Mass

d_male <- d_male %>%
  filter(!is.na(mass))
d_male <- center_data(d_male)

data<-data.frame(R=d_male$distance, 
                 A=d_male$host_c, 
                 B=d_male$wing2body_c, 
                 C=d_male$mass_c, 
                 X=d_male$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 3-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]      [,3]      [,4]       [,5]      [,6]       [,7]      
## AICs   4193.593  4193.958  4194.143  4194.17    4194.864  4195.053   4195.374  
## models 3         11        15        5          17        6          14        
## probs  0.1326704 0.1105416 0.1007816 0.09943557 0.0702871 0.06395294 0.05446059
##        [,8]      
## AICs   4195.531  
## models 18        
## probs  0.05035235
## 
## m3   R ~ C + (1 | X)
## m11  R ~ A * B + C + (1 | X)
## m15  R ~ A * B + B * C + (1 | X)
## m5   R ~ A + C + (1 | X)
## m17  R ~ A * B + A * C + B * C + (1 | X)
## m6   R ~ B + C + (1 | X)
## m14  R ~ A * B + A * C + (1 | X)
## m0   R ~ 1 + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 6
anova(m0, m3, test="Chisq") # Adding C improves fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m3: R ~ C + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m0    3 4195.5 4206.3 -2094.8   4189.5                       
## m3    4 4193.6 4207.9 -2092.8   4185.6 3.9376  1    0.04722 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m3, m5, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m3: R ~ C + (1 | X)
## m5: R ~ A + C + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m3    4 4193.6 4207.9 -2092.8   4185.6                     
## m5    5 4194.2 4212.1 -2092.1   4184.2 1.4233  1     0.2329
best.fit <- glm(distance ~ mass_c, family = Gamma(link = "log"), data = d_male) 
tidy_regression(best.fit, is_color=output_col)
## glm distance ~ mass_c Gamma(link = "log") d_male 
## AIC:  4187.136 
## (Intercept)  coeff:  7.7572867   Pr(>|t|):  1.252219e-179 *
## mass_c       coeff:  27.0665642  Pr(>|t|):  0.08062868 .
  • marginal positive effect of mass, where the heavier the male bug, the more likely the male bug flies far (this could be because of thorax wider so has more muscle power?)
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

Min From Inc Start

data<-data.frame(R=d_male$distance, 
                 A=d_male$min_from_IncStart_s, 
                 B=d_male$wing2body_c, 
                 C=d_male$mass_s, 
                 X=d_male$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 3-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]    [,2]      [,3]      [,4]       [,5]       [,6]       [,7]      
## AICs   4193.14 4193.199  4193.593  4194.866   4195.037   4195.053   4195.101  
## models 5       9         3         7          12         6          1         
## probs  0.15505 0.1505744 0.1235998 0.06540127 0.06004464 0.05958049 0.05817052
## 
## m5   R ~ A + C + (1 | X)
## m9   R ~ A * C + (1 | X)
## m3   R ~ C + (1 | X)
## m7   R ~ A + B + C + (1 | X)
## m12  R ~ A * C + B + (1 | X)
## m6   R ~ B + C + (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: 1
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 4193.1 4211.1 -2091.6   4183.1                     
## m9    6 4193.2 4214.7 -2090.6   4181.2 1.9414  1     0.1635
anova(m3, m5, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m3: R ~ C + (1 | X)
## m5: R ~ A + C + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m3    4 4193.6 4207.9 -2092.8   4185.6                     
## m5    5 4193.1 4211.1 -2091.6   4183.1 2.4534  1     0.1173
anova(m0, m3, test="Chisq") # Adding C improves fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m3: R ~ C + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m0    3 4195.5 4206.3 -2094.8   4189.5                       
## m3    4 4193.6 4207.9 -2092.8   4185.6 3.9376  1    0.04722 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
best.fit <- glmer(distance ~ mass_s + (1|chamber), family = Gamma(link = "log"), data = d_fem) 
## boundary (singular) fit: see ?isSingular
tidy_regression(best.fit, is_color=output_col)
## glmer distance ~ mass_s + (1 | chamber) d_fem Gamma(link = "log") 
## AIC:  1151.077 1151.077 
## (Intercept)  coeff:  8.0243169   Pr(>|t|):  3.773659e-302 *
## mass_s       coeff:  -0.2839124  Pr(>|t|):  0.2492655
  • no effect of mass
s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

Morphology

data<-data.frame(R=d_male$distance, A=d_male$wing2body, B=d_male$thorax_c, X=d_male$chamber)

model_script = paste0(source_path,"generic models-Gamma glmer 2-FF log link.R")
errors <- withWarnings(model_comparisonsAIC(model_script))
##        [,1]      [,2]     [,3]      [,4]      [,5]      
## AICs   4194.25   4195.531 4195.91   4196.567  4197.78   
## models 2         5        3         1         4         
## probs  0.4084932 0.215251 0.1781139 0.1282408 0.06990104
## 
## m2   R ~ B + (1 | X)
## m0   R ~ 1 + (1 | X)
## m3   R ~ A + B + (1 | X)
## m1   R ~ A + (1 | X)
## m4   R ~ A * B + (1 | X)
cat("Number of models that failed to converge:", length(errors$warnings))
## Number of models that failed to converge: 2
anova(m0, m2, test="Chisq") # Adding B marginally improves fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m0    3 4195.5 4206.3 -2094.8   4189.5                       
## m2    4 4194.2 4208.6 -2093.1   4186.2 3.2813  1    0.07007 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m0, m1, test="Chisq") # Adding A does not improve fit
## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m0    3 4195.5 4206.3 -2094.8   4189.5                     
## m1    4 4196.6 4210.9 -2094.3   4188.6 0.9642  1     0.3261
best.fit <- glmer(distance ~ thorax_c + (1|chamber), family = Gamma(link = "log"), data = d_male) 
tidy_regression(best.fit, is_color=output_col)
## glmer distance ~ thorax_c + (1 | chamber) d_male Gamma(link = "log") 
## AIC:  4194.25 4194.25 
## (Intercept)  coeff:  7.705892    Pr(>|t|):  0 *
## thorax_c     coeff:  1.3249771   Pr(>|t|):  0.06795104 .

*marginal positive effect of thorax where the wider the thorax the more likely the bug will disperse far.

s.test <- paste("pval: ", shapiro.test(residuals(best.fit))$p.value)
qqnorm(resid(best.fit))
qqline(resid(best.fit))
text(-1, 0.1, s.test)

d_male$min_from_IncStart_b <- 0
d_male$min_from_IncStart_b[d_male$min_from_IncStart > 240] = 1

Mmeans<-aggregate(distance~mass*wing2body*min_from_IncStart_b, data=d_male, FUN="mean")

Male_Plots = function() {par(mfrow=c(1,2))
  plot(Mmeans$wing2body, Mmeans$distance, 
       col=c(19,21)[as.factor(Mmeans$min_from_IncStart_b)], # green is 0 and 1 is blue
       pch=19,
       xlab="wing-to-body ratio",
       ylab="distance (m)") 
  mtext(expression(bold("A) Males | n = 266")), side=3, adj=0.01, line=0.5, cex=1.3)
  
  plot(Mmeans$mass, Mmeans$distance, 
       col=c(19,21)[as.factor(Mmeans$min_from_IncStart_b)], # green is 0 and 1 is blue
       pch=19,
       xlab="mass (g)",
       ylab="distance (m)") 
  legend(0.053, 21000, 
         legend=c(expression(italic("Morning")), 
                  expression(italic("Afternoon"))), 
         pch=c(19,19,21,21), col=c(19,21,19,21), cex=1)
  mtext(expression(bold("B) Males | n = 266")), side=3, adj=0.01, line=0.5, cex=1.3)
}
Male_Plots()

Mass doesn’t seem to matter, neither does thorax length, so plotted with wing2body ratio. Hard to see a pattern here.

Summary

  • (For bugs that flew some distance) Distance or delta distance has nothing to do with the wing2body ratio. It rarely has anything to do with sex or mass. But it seems to be heavily dependent on time of day.
Trial_1_Plots()

Trial_2_Plots()

Fem_Plots()

Male_Plots()