Binomial Modeling: flight_yes_no
Anastasia Bernat
3/30/2020 - 8/5/2020
Winter 2020 Flight Trials: Binomial 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.
All Trials
Testing glmer() and Covariates
Cleaning 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
"unique_flight_data.R",
"get_warnings.R",
"compare_models.R",
"regression_output.R", # Cleans regression outputs; prints in color or B&W
"AICprobabilities.R",
"get_Akaike_weights.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##
## Attaching package: 'chron'## The following objects are masked from 'package:lubridate':
##
## days, hours, minutes, seconds, yearsoutput_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 in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object length## 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))## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_all <- data[[1]]
data_tested <- data[[2]]
d <- create_delta_data(data_tested)d <- center_data(d, is_not_unique_data = FALSE)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthd## # A tibble: 281 x 86
## # Groups: ID, sex, population, site, host_plant, latitude, longitude,
## # total_eggs, beak, thorax, wing, body, w_morph, morph_notes, tested, host_c,
## # sex_c, w_morph_c, lat_c, sym_dist, sym_dist_s, total_eggs_c, beak_c,
## # thorax_c, thorax_s, body_c, wing_c, wing2body, wing2body_c [281]
## ID sex population site host_plant latitude longitude total_eggs beak
## <fct> <fct> <fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1 M Plantation… Arego… C. corind… 25.0 -80.6 NA 5.5
## 2 2 M Plantation… Arego… C. corind… 25.0 -80.6 NA 5.64
## 3 3 M Plantation… Arego… C. corind… 25.0 -80.6 NA 5.75
## 4 4 M Plantation… Arego… C. corind… 25.0 -80.6 NA 6.21
## 5 5 F Plantation… Arego… C. corind… 25.0 -80.6 209 7.47
## 6 6 M Plantation… Found… C. corind… 25.0 -80.6 NA 5.76
## 7 7 F Plantation… Found… C. corind… 25.0 -80.6 8 8.19
## 8 9 M Plantation… Found… C. corind… 25.0 -80.6 NA 6.1
## 9 10 M Plantation… Found… C. corind… 25.0 -80.6 NA 6.09
## 10 11 M Key Largo KLMRL C. corind… 25.1 -80.4 NA 5.88
## # … with 271 more rows, and 77 more variables: thorax <dbl>, wing <dbl>,
## # body <dbl>, w_morph <fct>, morph_notes <fct>, tested <fct>, host_c <dbl>,
## # sex_c <dbl>, w_morph_c <dbl>, lat_c <dbl>, sym_dist <dbl>,
## # sym_dist_s <dbl>, total_eggs_c <dbl>, beak_c <dbl>, thorax_c <dbl>,
## # thorax_s <dbl>, body_c <dbl>, wing_c <dbl>, wing2body <dbl>,
## # wing2body_c <dbl>, wing2body_s <dbl>, trial_type <list>, chamber <list>,
## # average_speed <list>, total_flight_time <list>, distance <list>,
## # shortest_flying_bout <list>, longest_flying_bout <list>,
## # recording_duration <list>, max_speed <list>, test_date <list>,
## # time_start <list>, time_end <list>, flew <list>, flight_type <list>,
## # mass <list>, EWM <list>, flew_b <list>, host_temp <list>, eggs_b <list>,
## # minute_duration <list>, minute_duration_c <list>, min_from_IncStart <dbl>,
## # min_from_IncStart_c <dbl>, min_from_IncStart_s <dbl>,
## # days_from_start <list>, days_from_start_c <list>, mass_c <list>,
## # mass_s <list>, avg_mass <dbl>, avg_mass_c <dbl>, avg_mass_s <dbl>,
## # trial_type_og <list>, dt_start_c <list>, dt_end_c <list>,
## # datetime_start <list>, datetime_end <list>, hr_start <list>, hr_end <list>,
## # num_flew <dbl>, num_notflew <dbl>, avg_days <dbl>, num_egg <dbl>,
## # egg_diff <dbl>, mass_diff <dbl>, flew_diff <dbl>, dist_diff <dbl>,
## # speed_diff <dbl>, mass_per <dbl>, avg_rec_dur <dbl>, avg_time_start <dbl>,
## # flight_case <fct>, egg_case <dbl>, flew_prob <dbl>, avg_days_c <dbl>,
## # avg_mass_logsqrt <dbl>, wing2body_logsqrt_i <dbl>Note
Modeling the dataset as it was before led to lots of converging issues that couldn’t be easily resolved. Now we are using a modified dataset, d, to model the data.
Testing Models and Covariates
# test mixed effect model
model <- glmer(cbind(num_flew, num_notflew)~sex_c*host_c + (1|population), data=d, family=binomial)## boundary (singular) fit: see ?isSingulartidy_regression(model, is_color=output_col) ####MLC: updated model, but won't converge with population as a RF ## glmer cbind(num_flew, num_notflew) ~ sex_c * host_c + (1 | population) d binomial
## AIC: 627.2934 627.2934
## (Intercept) coeff: 0.1066936 Pr(>|t|): 0.3161209
## sex_c coeff: -0.6143089 Pr(>|t|): 7.839931e-09 *
## host_c coeff: -0.0478021 Pr(>|t|): 0.6533351
## sex_c:host_c coeff: 0.1499523 Pr(>|t|): 0.1588626model <- glmer(cbind(num_flew, num_notflew)~sex_c*host_c + (1|site), data=d, family=binomial)
tidy_regression(model, is_color=output_col) ## glmer cbind(num_flew, num_notflew) ~ sex_c * host_c + (1 | site) d binomial
## AIC: 622.8293 622.8293
## (Intercept) coeff: 0.1133526 Pr(>|t|): 0.4310746
## sex_c coeff: -0.6459361 Pr(>|t|): 1.026995e-08 *
## host_c coeff: -0.0265116 Pr(>|t|): 0.8540307
## sex_c:host_c coeff: 0.1805152 Pr(>|t|): 0.1077054getME(model, "lower")## [1] 0## AB: What is a random factor here exactly? B/c it seems to matter which we use.
## AB: I tried this and it didn't work but is there a way to cbind() random factors so I can put in chamber, test date, or test time?Experimental, Biological, and Morphological Effects
## glm cbind(num_flew, num_notflew) ~ avg_mass_c binomial d
## AIC: 617.6888
## (Intercept) coeff: 0.3362491 Pr(>|t|): 0.0002083159 *
## avg_mass_c coeff: -30.9034972 Pr(>|t|): 2.466421e-13 *## glm cbind(num_flew, num_notflew) ~ total_eggs binomial d
## AIC: 167.0929
## (Intercept) coeff: 0.3465927 Pr(>|t|): 0.1799292
## total_eggs coeff: -0.0131785 Pr(>|t|): 3.659423e-06 *## glm cbind(num_flew, num_notflew) ~ beak_c binomial d
## AIC: 663.4421
## (Intercept) coeff: 0.337106 Pr(>|t|): 0.0001046852 *
## beak_c coeff: -0.3564658 Pr(>|t|): 4.393565e-05 *## glm cbind(num_flew, num_notflew) ~ thorax_c binomial d
## AIC: 665.7521
## (Intercept) coeff: 0.3376217 Pr(>|t|): 9.915633e-05 *
## thorax_c coeff: -1.1155472 Pr(>|t|): 0.000148351 *## glm cbind(num_flew, num_notflew) ~ body_c binomial d
## AIC: 674.9982
## (Intercept) coeff: 0.3334611 Pr(>|t|): 0.0001052822 *
## body_c coeff: -0.1930026 Pr(>|t|): 0.01824023 *## glm cbind(num_flew, num_notflew) ~ wing_c binomial d
## AIC: 679.9601
## (Intercept) coeff: 0.3307729 Pr(>|t|): 0.0001109575 *
## wing_c coeff: -0.086374 Pr(>|t|): 0.4125933Without Mass
R1 = d$num_flew
R2 = d$num_notflew
A = d$host_c
B = d$sex_c
C = d$sym_dist
D = d$avg_mass_c
X = d$population # or population get the same top model. Variance is 0 if use pop, not 0 if use site. However, variance is 0 if use site and avg_mass (down below).
data<-data.frame(R1, R2, A, B, C, D, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs 625.9441 627.2552 627.2934 627.9081 628.5906 629.0372
## models 2 4 8 6 11 7
## probs 0.282125 0.1464644 0.143692 0.1056744 0.07512103 0.06008558
##
## m2 cbind(R1, R2) ~ B + (1 | X)
## m4 cbind(R1, R2) ~ A + B + (1 | X)
## m8 cbind(R1, R2) ~ A * B + (1 | X)
## m6 cbind(R1, R2) ~ B + C + (1 | X)
## m11 cbind(R1, R2) ~ A * B + C + (1 | X)
## m7 cbind(R1, R2) ~ A + B + C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m4, m8, test="Chisq") # Adding A*B marginally improves fit if use pop as RF | marginally improves fit if use site as RF## Data: data
## Models:
## m4: cbind(R1, R2) ~ A + B + (1 | X)
## m8: cbind(R1, R2) ~ A * B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m4 4 627.26 641.81 -309.63 619.26
## m8 5 627.29 645.49 -308.65 617.29 1.9618 1 0.1613anova(m2, m4, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m2: cbind(R1, R2) ~ B + (1 | X)
## m4: cbind(R1, R2) ~ A + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 3 625.94 636.86 -309.97 619.94
## m4 4 627.26 641.81 -309.63 619.26 0.6889 1 0.4066anova(m2, m6, test="Chisq") # Adding C does not improve fit## Data: data
## Models:
## m2: cbind(R1, R2) ~ B + (1 | X)
## m6: cbind(R1, R2) ~ B + C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 3 625.94 636.86 -309.97 619.94
## m6 4 627.91 642.46 -309.95 619.91 0.036 1 0.8495nomass_model_all <- glmer(cbind(num_flew, num_notflew) ~ sex_c + (1|population), data=d, family=binomial)
tidy_regression(nomass_model_all, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ sex_c + (1 | population) d binomial
## AIC: 625.9441 625.9441
## (Intercept) coeff: 0.1272439 Pr(>|t|): 0.3043372
## sex_c coeff: -0.6885895 Pr(>|t|): 2.314013e-13 *nomass_model_all <- glmer(cbind(num_flew, num_notflew) ~ sex_c + (1|site), data=d, family=binomial)
tidy_regression(nomass_model_all, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ sex_c + (1 | site) d binomial
## AIC: 621.6971 621.6971
## (Intercept) coeff: 0.123132 Pr(>|t|): 0.338327
## sex_c coeff: -0.7326913 Pr(>|t|): 3.493679e-13 *- strong negative effect of sex where if F then less number of times will fly
With Mass
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1-RF + 4-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2]
## AICs 617.1318 617.1369
## models 26 51
## probs 0.05023765 0.05011018
##
## m26 cbind(R1, R2) ~ A * D + B + (1 | X)
## m51 cbind(R1, R2) ~ A * D + C * D + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m63, m85, test="Chisq") # Adding B*D does not improve fit## Data: data
## Models:
## m63: cbind(R1, R2) ~ A * D + C * D + B + (1 | X)
## m85: cbind(R1, R2) ~ A * D + B * D + C * D + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m63 8 617.54 646.65 -300.77 601.54
## m85 9 617.85 650.59 -299.92 599.85 1.6959 1 0.1928anova(m26, m36, test="Chisq") # Adding C does not improve fit | is sym dist important?## Data: data
## Models:
## m26: cbind(R1, R2) ~ A * D + B + (1 | X)
## m36: cbind(R1, R2) ~ A * D + B + C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m26 6 617.13 638.96 -302.57 605.13
## m36 7 618.25 643.72 -302.13 604.25 0.8803 1 0.3481anova(m12, m26, test="Chisq") # Adding A*D does improve fit## Data: data
## Models:
## m12: cbind(R1, R2) ~ A + B + D + (1 | X)
## m26: cbind(R1, R2) ~ A * D + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m12 5 619.44 637.64 -304.72 609.44
## m26 6 617.13 638.96 -302.57 605.13 4.3129 1 0.03782 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wanted to check whether sym_dist is making much of an impact overall and you can see it’s far from significant.
tidy_regression(glmer(cbind(num_flew, num_notflew) ~ sym_dist + (1|population), data=d, family=binomial), is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ sym_dist + (1 | population) d binomial
## AIC: 681.6667 681.6667
## (Intercept) coeff: 0.2488054 Pr(>|t|): 0.1197266
## sym_dist coeff: 0.0599719 Pr(>|t|): 0.5496975# model m26 might be the better model than this because even adding C did not improve fit
# and the difference in AIC between the two isn't huge.
mass_model_all <- glmer(cbind(num_flew, num_notflew) ~ host_c * avg_mass_c + sex_c + (1|population), data=d, family=binomial)## boundary (singular) fit: see ?isSingulartidy_regression(mass_model_all, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ host_c * avg_mass_c + sex_c + (1 | population) d binomial
## AIC: 617.1318 617.1318
## (Intercept) coeff: 0.2123599 Pr(>|t|): 0.05627584 .
## host_c coeff: -0.1174916 Pr(>|t|): 0.2583448
## avg_mass_c coeff: -16.4082233 Pr(>|t|): 0.03718048 *
## sex_c coeff: -0.2657803 Pr(>|t|): 0.1209117
## host_c:avg_mass_c coeff: 9.9767364 Pr(>|t|): 0.03448337 *mass_model_all <- glmer(cbind(num_flew, num_notflew) ~ host_c * avg_mass_c + sym_dist * avg_mass_c + sex_c + (1|population), data=d, family=binomial)## boundary (singular) fit: see ?isSingulartidy_regression(mass_model_all, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ host_c * avg_mass_c + sym_dist * avg_mass_c + sex_c + (1 | population) d binomial
## AIC: 617.5419 617.5419
## (Intercept) coeff: 0.1290136 Pr(>|t|): 0.4444862
## host_c coeff: -0.1982126 Pr(>|t|): 0.1446491
## avg_mass_c coeff: -9.6293052 Pr(>|t|): 0.2700795
## sym_dist coeff: 0.043203 Pr(>|t|): 0.7390353
## sex_c coeff: -0.219919 Pr(>|t|): 0.2054395
## host_c:avg_mass_c coeff: 15.1162124 Pr(>|t|): 0.005333247 *
## avg_mass_c:sym_dist coeff: -10.0438672 Pr(>|t|): 0.1098937Now there’s a conflicting interaction between average mass and host*avg_mass…and there’s
- no effect of host
- no effect of avg_mass_c
- no effect of sym_dist
but effect of sex and interactions but the interactions are conflicting…
How do you graph these models? How would you graph cbind(num_flew, num_notflew)?
plot( cbind(num_flew, num_notflew) ~ cbind(host_c, avg_mass), data=d , col=rangi2)
How does num_flew or num_notflew vary with host plant?
coplot(num_flew ~ avg_mass_c | host_c, data=d)
coplot(num_notflew~ avg_mass_c | host_c, data=d)
Plot the interaction between host*avg_mass:
Consistent with the model and graphs: If GRT and had an average mass below the mean of the population, then you flew less times. But if were GRT and had an average mass above the mean of the population, then you flew more times.
Overall mass vs. num_flew trend:
gf_point(num_flew ~ avg_mass, data=d) +
geom_smooth(method='lm')## `geom_smooth()` using formula 'y ~ x'
Then break it between host and average mass:
p1 <- gf_point(num_flew ~ avg_mass, col=~host_plant, data=d[d$avg_mass_c < 0,], title="(a) avgmass_c <0") +
geom_smooth(method='lm')
p2 <- gf_point(num_flew ~ avg_mass, col=~host_plant, data=d[d$avg_mass_c > 0,], title="(b) avgmass_c >0") +
geom_smooth(method='lm')
grid.arrange(p1,p2, nrow=2)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
p3 <- gf_point(num_notflew ~ avg_mass, col=~host_plant, data=d[d$avg_mass_c < 0,], title="(a) avgmass_c <0") +
geom_smooth(method='lm')
p4 <- gf_point(num_notflew ~ avg_mass, col=~host_plant, data=d[d$avg_mass_c > 0,], title="(b) avgmass_c >0") +
geom_smooth(method='lm')
grid.arrange(p3,p4, nrow=2)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
GRT flyers: Even when they weighed more than the average mass of the population, they flew more times, which breaks the overall trend between number times bug flew vs. average mass which is typically negative.
Morphology
d <- d %>%
filter(!is.na(body))
d$thorax_c <- d$thorax - mean(d$thorax)
d$wing_c <- d$wing - mean(d$wing)
d$body_c <- d$body - mean(d$body)
R1 = d$num_flew
R2 = d$num_notflew
A = d$thorax_c
B = d$body_c
C = d$wing_c
X = d$population
data<-data.frame(R1, R2, A, B, C, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3]
## AICs 625.3875 625.6398 627.1639
## models 16 15 17
## probs 0.3972829 0.3501824 0.1634401
##
## m16 cbind(R1, R2) ~ A * C + B * C + (1 | X)
## m15 cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m17 cbind(R1, R2) ~ 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: 0anova(m15, m17, test="Chisq") # Adding A*C does not improve fit## Data: data
## Models:
## m15: cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m17: cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m15 7 625.64 651.11 -305.82 611.64
## m17 8 627.16 656.27 -305.58 611.16 0.476 1 0.4902anova(m13, m15, test="Chisq") # Adding A*B does improve fit## Data: data
## Models:
## m13: cbind(R1, R2) ~ B * C + A + (1 | X)
## m15: cbind(R1, R2) ~ A * B + B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m13 6 629.74 651.57 -308.87 617.74
## m15 7 625.64 651.11 -305.82 611.64 6.0983 1 0.01353 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1morph_model_all <- glmer(cbind(num_flew, num_notflew) ~ thorax_c * body_c + body_c * wing_c + (1|population), data=d, family=binomial)
tidy_regression(morph_model_all, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ thorax_c * body_c + body_c * wing_c + (1 | population) d binomial
## AIC: 625.6398 625.6398
## (Intercept) coeff: 0.3930088 Pr(>|t|): 0.01062241 *
## thorax_c coeff: -2.6397658 Pr(>|t|): 0.004835852 *
## body_c coeff: -1.2315294 Pr(>|t|): 0.0120045 *
## wing_c coeff: 2.3022068 Pr(>|t|): 1.136889e-05 *
## thorax_c:body_c coeff: 1.5450568 Pr(>|t|): 0.01672995 *
## body_c:wing_c coeff: -0.6778824 Pr(>|t|): 0.003735883 *strong negative effect of thorax where the wider the thorax, the less times the bug flies
strong negative effect of body where the longer the body, the less times a bug flies
strong positive effect of wing where the longer the wing, the more times the bug flies
strong positive effect of thorax*body where the longer the body and wider the thorax, the more times a bug flies (hm seems conflicting to the single effects)
negative effect of of body*wing where the longer the body and wing, the less times the bug flies
R1 = d$num_flew
R2 = d$num_notflew
A = d$thorax_c
B = d$wing2body_c
X = d$population
data<-data.frame(R1, R2, A, B, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1-RF + 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2]
## AICs 629.6603 630.291
## models 3 4
## probs 0.5780482 0.4217172
##
## m3 cbind(R1, R2) ~ A + B + (1 | X)
## m4 cbind(R1, R2) ~ A * B + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m3, m4, test="Chisq") # Adding A*B improves fit## Data: data
## Models:
## m3: cbind(R1, R2) ~ A + B + (1 | X)
## m4: cbind(R1, R2) ~ A * B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m3 4 629.66 644.21 -310.83 621.66
## m4 5 630.29 648.48 -310.14 620.29 1.3694 1 0.2419morph_model_all2 <- glmer(cbind(num_flew, num_notflew) ~ wing2body_c * thorax_c + (1|population), data=d, family=binomial)
tidy_regression(morph_model_all2, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing2body_c * thorax_c + (1 | population) d binomial
## AIC: 630.291 630.291
## (Intercept) coeff: 0.2909256 Pr(>|t|): 0.05135432 .
## wing2body_c coeff: 31.3397207 Pr(>|t|): 4.251342e-08 *
## thorax_c coeff: -1.2838472 Pr(>|t|): 2.764461e-05 *
## wing2body_c:thorax_c coeff: -23.6737268 Pr(>|t|): 0.2432869- strong positive effect of wing2body where the greater the ratio the more times a bug flies
- negative effect of thorax where the the wider the thorax the less times a bug flies
- marginal strong negative effect of wing2body*thorax where the wider the thorax and large the ratio, the less times a bug flies
# check for collinearity
covariates <- data.frame(d$thorax, d$wing2body, d$avg_mass)
pairs(covariates)
Exponential relationships between mass and thorax. Thorax can only get so big - mass can increase but thorax caps out (top right).
covariates <- data.frame(d$thorax, d$wing, d$body, d$avg_mass)
pairs(covariates)
Seems like once morphology caps out so does mass not stretch too far.
Morphology + Mass + Sex
R1 = d$num_flew
R2 = d$num_notflew
A = d$thorax_c
B = d$wing2body_c
C = d$avg_mass_c
X = d$population
data<-data.frame(R1, R2, A, B, C, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5]
## AICs 592.4714 594.4679 594.5592 595.9454 596.2326
## models 15 17 11 7 14
## probs 0.4465144 0.1645513 0.1572123 0.07860943 0.06809484
##
## m15 cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m17 cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)
## m11 cbind(R1, R2) ~ A * B + C + (1 | X)
## m7 cbind(R1, R2) ~ A + B + C + (1 | X)
## m14 cbind(R1, R2) ~ A * B + A * C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 11anova(m11, m15, test="Chisq") # Adding B*C improves fit## Data: data
## Models:
## m11: cbind(R1, R2) ~ A * B + C + (1 | X)
## m15: cbind(R1, R2) ~ A * B + B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m11 6 594.56 616.39 -291.28 582.56
## m15 7 592.47 617.94 -289.24 578.47 4.0877 1 0.0432 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m15, m17, test="Chisq") # Adding A*C does not improve fit## Data: data
## Models:
## m15: cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m17: cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m15 7 592.47 617.94 -289.24 578.47
## m17 8 594.47 623.57 -289.23 578.47 0.0035 1 0.9528multi_model_all <- glmer(cbind(num_flew, num_notflew) ~ thorax_c * wing2body_c + wing2body_c * avg_mass + (1|population), data=d, family=binomial) ## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluestidy_regression(multi_model_all, is_color=output_col)## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## glmer cbind(num_flew, num_notflew) ~ thorax_c * wing2body_c + wing2body_c * avg_mass + (1 | population) d binomial
## AIC: 592.4714 592.4714
## (Intercept) coeff: 2.4526127 Pr(>|t|): 1.145317e-08 *
## thorax_c coeff: 1.2176489 Pr(>|t|): 0.03003484 *
## wing2body_c coeff: -30.493285 Pr(>|t|): 0.2231488
## avg_mass coeff: -40.650045 Pr(>|t|): 2.257211e-07 *
## thorax_c:wing2body_c coeff: -86.1556666 Pr(>|t|): 0.008399759 *
## wing2body_c:avg_mass coeff: 876.2448761 Pr(>|t|): 0.05366202 .coefficients are huge. Probably a collinearity issue between morph measurements * mass.
R1 = d$num_flew
R2 = d$num_notflew
A = d$thorax_c
B = d$wing2body_c
C = d$sex_c
X = d$population
data<-data.frame(R1, R2, A, B, C, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs 612.8582 614.6533 614.8412 614.8563 615.865 616.4857
## models 6 11 10 7 15 14
## probs 0.3388222 0.1380988 0.1257107 0.1247697 0.0753452 0.05524289
##
## m6 cbind(R1, R2) ~ B + C + (1 | X)
## m11 cbind(R1, R2) ~ A * B + C + (1 | X)
## m10 cbind(R1, R2) ~ B * C + (1 | X)
## m7 cbind(R1, R2) ~ A + B + C + (1 | X)
## m15 cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m14 cbind(R1, R2) ~ A * B + A * C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m11, m15, test="Chisq") # Adding B*C does not improve fit## Data: data
## Models:
## m11: cbind(R1, R2) ~ A * B + C + (1 | X)
## m15: cbind(R1, R2) ~ A * B + B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m11 6 614.65 636.48 -301.33 602.65
## m15 7 615.87 641.33 -300.93 601.87 0.7882 1 0.3746anova(m11, m14, test="Chisq") # Adding A*C does not improve fit## Data: data
## Models:
## m11: cbind(R1, R2) ~ A * B + C + (1 | X)
## m14: cbind(R1, R2) ~ A * B + A * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m11 6 614.65 636.48 -301.33 602.65
## m14 7 616.49 641.95 -301.24 602.49 0.1675 1 0.6823anova(m7, m11, test="Chisq") # Adding A*B does improve fit## Data: data
## Models:
## m7: cbind(R1, R2) ~ A + B + C + (1 | X)
## m11: cbind(R1, R2) ~ A * B + C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m7 5 614.86 633.05 -302.43 604.86
## m11 6 614.65 636.48 -301.33 602.65 2.203 1 0.1377multi_model_all2 <- glmer(cbind(num_flew, num_notflew) ~ thorax_c * wing2body_c + sex_c + (1|population), data=d, family=binomial)
tidy_regression(multi_model_all2, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ thorax_c * wing2body_c + sex_c + (1 | population) d binomial
## AIC: 614.6533 614.6533
## (Intercept) coeff: 0.1267477 Pr(>|t|): 0.3645612
## thorax_c coeff: 0.0340484 Pr(>|t|): 0.937718
## wing2body_c coeff: 19.1759863 Pr(>|t|): 0.002452555 *
## sex_c coeff: -0.6001514 Pr(>|t|): 3.22563e-05 *
## thorax_c:wing2body_c coeff: -29.0188446 Pr(>|t|): 0.136408no effect of thorax
strong positve effect of ratio where the larger the ratio the more times a bug flew
negative effect of sex where if F then less times bug will fly
strong negative effect of thorax*ratio where the larger the ratio and wider the thorax, the less times a bug flew.
# weak positive relationship between thorax vs. wing2body
p1 <- gf_point(thorax ~ wing2body, col=~num_flew, data=d) +
geom_smooth(method='lm')
p2 <- gf_point(thorax ~ wing2body, data=d[d$num_flew == 2,], title="num_flew=2") +
geom_smooth(method='lm')
p3 <- gf_point(thorax ~ wing2body, data=d[d$num_flew == 1,], title="num_flew=1") +
geom_smooth(method='lm')
p4 <- gf_point(thorax ~ wing2body, data=d[d$num_flew == 0,], title="num_flew=0") +
geom_smooth(method='lm')
grid.arrange(p1,p2,p3,p4, ncol=2)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
As wing2body and thorax increases, num_flew decreases too, but its really the thorax shrinking that seems to be driving this.
gf_point(num_flew ~ thorax, data=d) +
geom_smooth(method='lm')## `geom_smooth()` using formula 'y ~ x'
gf_point(num_flew ~ wing2body, data=d) +
geom_smooth(method='lm')## `geom_smooth()` using formula 'y ~ x'
R1 = d$num_flew
R2 = d$num_notflew
A = d$thorax_s
B = d$wing2body_s
C = d$sex_c
D = d$avg_mass_s
X = d$population
data<-data.frame(R1, R2, A, B, C, D, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1-RF + 4-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5]
## AICs 592.4714 593.9312 594.1685 594.4679 594.5592
## models 43 74 58 72 23
## probs 0.1607423 0.07747159 0.06880296 0.05923741 0.05659541
##
## m43 cbind(R1, R2) ~ A * B + B * D + (1 | X)
## m74 cbind(R1, R2) ~ A * B + B * C + B * D + (1 | X)
## m58 cbind(R1, R2) ~ A * B + B * D + C + (1 | X)
## m72 cbind(R1, R2) ~ A * B + A * D + B * D + (1 | X)
## m23 cbind(R1, R2) ~ A * B + D + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 2# Less fail if use standardized variables. Also, all the top models have multiple interactions b/c I'm guessing, they're correlated in some ways.anova(m58, m74, test="Chisq") # Adding B*C marginally improves fit## Data: data
## Models:
## m58: cbind(R1, R2) ~ A * B + B * D + C + (1 | X)
## m74: cbind(R1, R2) ~ A * B + B * C + B * D + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m58 8 594.17 623.28 -289.08 578.17
## m74 9 593.93 626.68 -287.97 575.93 2.2373 1 0.1347anova(m74, m94, test="Chisq") # Adding A*C does not improve fit## Data: data
## Models:
## m74: cbind(R1, R2) ~ A * B + B * C + B * D + (1 | X)
## m94: cbind(R1, R2) ~ A * B + A * C + B * C + B * D + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m74 9 593.93 626.68 -287.97 575.93
## m94 10 595.42 631.80 -287.71 575.42 0.5107 1 0.4748multi_model <- glmer(cbind(num_flew, num_notflew) ~ thorax_s * wing2body_s + wing2body_s * avg_mass_s + sex_c + (1|population), data=d, family=binomial)
tidy_regression(multi_model, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ thorax_s * wing2body_s + wing2body_s * avg_mass_s + sex_c + (1 | population) d binomial
## AIC: 594.1685 594.1685
## (Intercept) coeff: 0.2729026 Pr(>|t|): 0.08862536 .
## thorax_s coeff: 0.373077 Pr(>|t|): 0.02697113 *
## wing2body_s coeff: 0.272076 Pr(>|t|): 0.02646234 *
## avg_mass_s coeff: -0.872786 Pr(>|t|): 0.0001663526 *
## sex_c coeff: -0.1048923 Pr(>|t|): 0.5810486
## thorax_s:wing2body_s coeff: -0.4727165 Pr(>|t|): 0.00840624 *
## wing2body_s:avg_mass_s coeff: 0.3742382 Pr(>|t|): 0.05279777 .Summary
tidy_regression(nomass_model_all, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ sex_c + (1 | site) d binomial
## AIC: 621.6971 621.6971
## (Intercept) coeff: 0.123132 Pr(>|t|): 0.338327
## sex_c coeff: -0.7326913 Pr(>|t|): 3.493679e-13 *tidy_regression(mass_model_all, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ host_c * avg_mass_c + sym_dist * avg_mass_c + sex_c + (1 | population) d binomial
## AIC: 617.5419 617.5419
## (Intercept) coeff: 0.1290136 Pr(>|t|): 0.4444862
## host_c coeff: -0.1982126 Pr(>|t|): 0.1446491
## avg_mass_c coeff: -9.6293052 Pr(>|t|): 0.2700795
## sym_dist coeff: 0.043203 Pr(>|t|): 0.7390353
## sex_c coeff: -0.219919 Pr(>|t|): 0.2054395
## host_c:avg_mass_c coeff: 15.1162124 Pr(>|t|): 0.005333247 *
## avg_mass_c:sym_dist coeff: -10.0438672 Pr(>|t|): 0.1098937Why this host*avg_mass interaction? GRT flyers: Even when they weighed more than the average mass of the population, they flew more times, which breaks the overall trend between number times bug flew vs. average mass which is typically negative. (Refer to plots in the “All Trials” tab)
tidy_regression(morph_model_all, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ thorax_c * body_c + body_c * wing_c + (1 | population) d binomial
## AIC: 625.6398 625.6398
## (Intercept) coeff: 0.3930088 Pr(>|t|): 0.01062241 *
## thorax_c coeff: -2.6397658 Pr(>|t|): 0.004835852 *
## body_c coeff: -1.2315294 Pr(>|t|): 0.0120045 *
## wing_c coeff: 2.3022068 Pr(>|t|): 1.136889e-05 *
## thorax_c:body_c coeff: 1.5450568 Pr(>|t|): 0.01672995 *
## body_c:wing_c coeff: -0.6778824 Pr(>|t|): 0.003735883 *tidy_regression(morph_model_all2, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing2body_c * thorax_c + (1 | population) d binomial
## AIC: 630.291 630.291
## (Intercept) coeff: 0.2909256 Pr(>|t|): 0.05135432 .
## wing2body_c coeff: 31.3397207 Pr(>|t|): 4.251342e-08 *
## thorax_c coeff: -1.2838472 Pr(>|t|): 2.764461e-05 *
## wing2body_c:thorax_c coeff: -23.6737268 Pr(>|t|): 0.2432869Why this wing2body*thorax interaction? As wing2body and thorax increases, num_flew decreases too, but its really the thorax shrinking that seems to be driving this.
tidy_regression(multi_model_all, is_color=output_col)## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## glmer cbind(num_flew, num_notflew) ~ thorax_c * wing2body_c + wing2body_c * avg_mass + (1 | population) d binomial
## AIC: 592.4714 592.4714
## (Intercept) coeff: 2.4526127 Pr(>|t|): 1.145317e-08 *
## thorax_c coeff: 1.2176489 Pr(>|t|): 0.03003484 *
## wing2body_c coeff: -30.493285 Pr(>|t|): 0.2231488
## avg_mass coeff: -40.650045 Pr(>|t|): 2.257211e-07 *
## thorax_c:wing2body_c coeff: -86.1556666 Pr(>|t|): 0.008399759 *
## wing2body_c:avg_mass coeff: 876.2448761 Pr(>|t|): 0.05366202 .tidy_regression(multi_model_all2, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ thorax_c * wing2body_c + sex_c + (1 | population) d binomial
## AIC: 614.6533 614.6533
## (Intercept) coeff: 0.1267477 Pr(>|t|): 0.3645612
## thorax_c coeff: 0.0340484 Pr(>|t|): 0.937718
## wing2body_c coeff: 19.1759863 Pr(>|t|): 0.002452555 *
## sex_c coeff: -0.6001514 Pr(>|t|): 3.22563e-05 *
## thorax_c:wing2body_c coeff: -29.0188446 Pr(>|t|): 0.136408Huge coefficients! I’m sure this is all due to collinearity. Thoughts?
All Trials - Old Method
Testing glmer() and Covariates
Cleaning the Data
data <- read_flight_data("data/all_flight_data-Winter2020.csv")## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object length
## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_all <- data[[1]]
data_tested <- data[[2]]Testing Models and Covariates
# test mixed effect model
model <- glmer(flew_b~sex_c*host_c + (1|ID) + (1|trial_type), data=data_tested, family=binomial)
summary(model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: flew_b ~ sex_c * host_c + (1 | ID) + (1 | trial_type)
## Data: data_tested
##
## AIC BIC logLik deviance df.resid
## 712.1 738.6 -350.0 700.1 608
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6403 -0.2638 0.1487 0.2755 1.2783
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 10.8325 3.2913
## trial_type (Intercept) 0.4311 0.6566
## Number of obs: 614, groups: ID, 333; trial_type, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3013 0.5557 -0.542 0.5877
## sex_c -1.9269 0.5903 -3.264 0.0011 **
## host_c -0.2048 0.2907 -0.704 0.4812
## sex_c:host_c 0.5212 0.3180 1.639 0.1012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sex_c host_c
## sex_c 0.232
## host_c 0.264 0.206
## sex_c:hst_c -0.032 -0.179 0.151getME(model, "lower")## [1] 0 0Experimental Set-Up Effects
## glm flew_b ~ chamber binomial data_tested
## AIC: 851.7226
## chamberA-2 coeff: 0.2113091 Pr(>|t|): 0.5455023
## chamberA-3 coeff: 0.5039052 Pr(>|t|): 0.1528972
## chamberA-4 coeff: 0.1625189 Pr(>|t|): 0.6163908
## chamberB-1 coeff: 0.5007753 Pr(>|t|): 0.1949202
## chamberB-2 coeff: 0.504556 Pr(>|t|): 0.1437077
## chamberB-3 coeff: -0.1000835 Pr(>|t|): 0.7718034
## chamberB-4 coeff: 0.1625189 Pr(>|t|): 0.6435895## glm flew_b ~ days_from_start_c binomial data_tested
## AIC: 837.6227
## (Intercept) coeff: 0.2392821 Pr(>|t|): 0.003488978 *
## days_from_start_c coeff: -0.035453 Pr(>|t|): 0.002742396 *## glm flew_b ~ min_from_IncStart_c binomial data_tested
## AIC: 846.5487
## (Intercept) coeff: 0.2356811 Pr(>|t|): 0.003738474 *
## min_from_IncStart_c coeff: 0.0002404 Pr(>|t|): 0.6770982Biological Effects
## glm flew_b ~ mass_c binomial data_tested
## AIC: 760.7559
## (Intercept) coeff: 0.2222918 Pr(>|t|): 0.01123772 *
## mass_c coeff: -33.9300431 Pr(>|t|): 1.462452e-15 *## glm flew_b ~ total_eggs binomial data_tested
## AIC: 209.4142
## (Intercept) coeff: 0.0580002 Pr(>|t|): 0.8111337
## total_eggs coeff: -0.0117451 Pr(>|t|): 1.730662e-05 *## glm flew_b ~ eggs_b binomial data_tested
## AIC: 723.9316
## (Intercept) coeff: 0.7344885 Pr(>|t|): 6.73571e-14 *
## eggs_b coeff: -2.40562 Pr(>|t|): 1.45431e-21 *Morphology Effects
## glm flew_b ~ beak_c binomial data_tested
## AIC: 821.5512
## (Intercept) coeff: 0.2400935 Pr(>|t|): 0.003814044 *
## beak_c coeff: -0.4085962 Pr(>|t|): 9.565263e-07 *## glm flew_b ~ thorax_c binomial data_tested
## AIC: 826.5349
## (Intercept) coeff: 0.2413545 Pr(>|t|): 0.003501085 *
## thorax_c coeff: -1.2298779 Pr(>|t|): 1.11169e-05 *## glm flew_b ~ body_c binomial data_tested
## AIC: 837.7504
## (Intercept) coeff: 0.2386438 Pr(>|t|): 0.003567537 *
## body_c coeff: -0.2303836 Pr(>|t|): 0.003002383 *## glm flew_b ~ wing_c binomial data_tested
## AIC: 844.687
## (Intercept) coeff: 0.2363842 Pr(>|t|): 0.003691238 *
## wing_c coeff: -0.1417783 Pr(>|t|): 0.1546876Without Mass (Test)
# Remove any missing masses
data_tested <- data_tested %>%
filter(!is.na(mass))
data_tested <- center_data(data_tested)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthR = data_tested$flew_b
A = data_tested$host_c
B = data_tested$sex_c
C = data_tested$sym_dist
D = data_tested$mass_c
X = data_tested$ID
Y = data_tested$trial_type
data<-data.frame(R, A, B, C, D, X, Y)
model_script = paste0(source_path,"generic models-binomial glmer 2 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs 707.3685 708.3419 709.0323 709.3669 710.1521 710.8463
## models 46 40 42 49 44 45
## probs 0.2925964 0.1798529 0.127345 0.107727 0.07274823 0.05141327
##
## m46 R ~ A * B + (1 | X) + (1 | Y)
## m40 R ~ B + (1 | X) + (1 | Y)
## m42 R ~ A + B + (1 | X) + (1 | Y)
## m49 R ~ A * B + C + (1 | X) + (1 | Y)
## m44 R ~ B + C + (1 | X) + (1 | Y)
## m45 R ~ A + B + C + (1 | X) + (1 | Y)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 16anova(m42, m46, test="Chisq") # Adding A*B marginally improves fit## Data: data
## Models:
## m42: R ~ A + B + (1 | X) + (1 | Y)
## m46: R ~ A * B + (1 | X) + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m42 5 709.03 731.11 -349.52 699.03
## m46 6 707.37 733.86 -347.68 695.37 3.6638 1 0.05561 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m40, m42, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m40: R ~ B + (1 | X) + (1 | Y)
## m42: R ~ A + B + (1 | X) + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m40 4 708.34 726.00 -350.17 700.34
## m42 5 709.03 731.11 -349.52 699.03 1.3095 1 0.2525nomass_model_all_old<-glmer(flew_b~sex_c*host_c + (1|ID) + (1|trial_type), family=binomial, data=data_tested)
tidy_regression(nomass_model_all_old, is_color=output_col)## glmer flew_b ~ sex_c * host_c + (1 | ID) + (1 | trial_type) data_tested binomial
## AIC: 707.3685 707.3685
## (Intercept) coeff: -0.3215517 Pr(>|t|): 0.5754789
## sex_c coeff: -2.0261213 Pr(>|t|): 0.003378218 *
## host_c coeff: -0.190949 Pr(>|t|): 0.5291606
## sex_c:host_c coeff: 0.570523 Pr(>|t|): 0.09505378 .- strong negative effect of sex where if F then less likely to fly
- no effect of host
- marginal effect of host*sex
Probbaly females from GRT are flyers becuse there are the colonizers of the population.
With Mass (Test)
# model_script = paste0(source_path,"generic models-binomial glmer 2-RF + 4-FF.R")
# errors = catch_warnings(model_comparisonsAIC(model_script))
# cat("Number of models that failed to converge: ", length(errors$warnings))
# too many models failed to even keep runningMorphology
data_tested <- data_tested %>%
filter(!is.na(body))
data_tested <- center_data(data_tested)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthR = data_tested$flew_b
A = data_tested$thorax_c
B = data_tested$body_c
C = data_tested$wing_c
X = data_tested$ID
Y = data_tested$trial_type
data<-data.frame(R, A, B, C, X, Y)
model_script = paste0(source_path,"generic models-binomial glmer 2 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4]
## AICs 723.4441 724.4474 725.3581 725.4187
## models 53 54 51 55
## probs 0.3625276 0.2195278 0.1392247 0.1350739
##
## m53 R ~ A * B + B * C + (1 | X) + (1 | Y)
## m54 R ~ A * C + B * C + (1 | X) + (1 | Y)
## m51 R ~ B * C + A + (1 | X) + (1 | Y)
## m55 R ~ A * B + A * C + B * C + (1 | X) + (1 | Y)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 6anova(m51, m53, test="Chisq") # Adding A*B improves fit## Data: data
## Models:
## m51: R ~ B * C + A + (1 | X) + (1 | Y)
## m53: R ~ A * B + B * C + (1 | X) + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m51 7 725.36 756.26 -355.68 711.36
## m53 8 723.44 758.76 -353.72 707.44 3.914 1 0.04788 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1morph_model_all_old<-glmer(flew_b~ thorax_c * body_c + body_c * wing_c + (1|ID), family=binomial, data=data_tested)
tidy_regression(morph_model_all_old, is_color=output_col) ## glmer flew_b ~ thorax_c * body_c + body_c * wing_c + (1 | ID) data_tested binomial
## AIC: 734.9448 734.9448
## (Intercept) coeff: 0.6537184 Pr(>|t|): 0.004835281 *
## thorax_c coeff: -4.0530068 Pr(>|t|): 0.02222891 *
## body_c coeff: -2.5044461 Pr(>|t|): 0.01314011 *
## wing_c coeff: 4.1255961 Pr(>|t|): 0.0002256748 *
## thorax_c:body_c coeff: 2.1987885 Pr(>|t|): 0.06018357 .
## body_c:wing_c coeff: -1.0819163 Pr(>|t|): 0.01165561 *R = data_tested$flew_b
A = data_tested$thorax_c
B = data_tested$wing2body_c
X = data_tested$ID
Y = data_tested$trial_type
data<-data.frame(R, A, B, C, X, Y)
model_script = paste0(source_path,"generic models-binomial glmer 2-RF + 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2]
## AICs 724.7296 725.3385
## models 14 13
## probs 0.5735094 0.4229865
##
## m14 R ~ A * B + (1 | X) + (1 | Y)
## m13 R ~ A + B + (1 | X) + (1 | Y)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 2anova(m13, m14, test="Chisq") # Adding A*B does not improve fit## Data: data
## Models:
## m13: R ~ A + B + (1 | X) + (1 | Y)
## m14: R ~ A * B + (1 | X) + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m13 5 725.34 747.41 -357.67 715.34
## m14 6 724.73 751.22 -356.36 712.73 2.6089 1 0.1063morph_model_all_old2<-glmer(flew_b~ thorax_c + wing2body_c + (1|trial_type) + (1|ID), family=binomial, data=data_tested)
tidy_regression(morph_model_all_old2, is_color=output_col) # coefficients are huge...## glmer flew_b ~ thorax_c + wing2body_c + (1 | trial_type) + (1 | ID) data_tested binomial
## AIC: 725.3385 725.3385
## (Intercept) coeff: 0.40038 Pr(>|t|): 0.3800688
## thorax_c coeff: -3.3927389 Pr(>|t|): 0.0002105466 *
## wing2body_c coeff: 72.7414372 Pr(>|t|): 1.569864e-05 *Summary
tidy_regression(nomass_model_all_old, is_color=output_col)## glmer flew_b ~ sex_c * host_c + (1 | ID) + (1 | trial_type) data_tested binomial
## AIC: 707.3685 707.3685
## (Intercept) coeff: -0.3215517 Pr(>|t|): 0.5754789
## sex_c coeff: -2.0261213 Pr(>|t|): 0.003378218 *
## host_c coeff: -0.190949 Pr(>|t|): 0.5291606
## sex_c:host_c coeff: 0.570523 Pr(>|t|): 0.09505378 .It’s possible that females from GRT are more likley to fly because they are colonizers.
# There isn't a mass_model_all_old because too many models failed to converge tidy_regression(morph_model_all_old, is_color=output_col) # same morph model as the morph_model_all## glmer flew_b ~ thorax_c * body_c + body_c * wing_c + (1 | ID) data_tested binomial
## AIC: 734.9448 734.9448
## (Intercept) coeff: 0.6537184 Pr(>|t|): 0.004835281 *
## thorax_c coeff: -4.0530068 Pr(>|t|): 0.02222891 *
## body_c coeff: -2.5044461 Pr(>|t|): 0.01314011 *
## wing_c coeff: 4.1255961 Pr(>|t|): 0.0002256748 *
## thorax_c:body_c coeff: 2.1987885 Pr(>|t|): 0.06018357 .
## body_c:wing_c coeff: -1.0819163 Pr(>|t|): 0.01165561 *tidy_regression(morph_model_all_old2, is_color=output_col) ## glmer flew_b ~ thorax_c + wing2body_c + (1 | trial_type) + (1 | ID) data_tested binomial
## AIC: 725.3385 725.3385
## (Intercept) coeff: 0.40038 Pr(>|t|): 0.3800688
## thorax_c coeff: -3.3927389 Pr(>|t|): 0.0002105466 *
## wing2body_c coeff: 72.7414372 Pr(>|t|): 1.569864e-05 *This morph model is a better fit model than the morph_model_all_old model, but HUGE coefficient for wing2body_c. Why?
Trial 1
Testing glmer() with pop and chamber
Cleaning the Data
data <- read_flight_data("data/all_flight_data-Winter2020.csv")## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object length
## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_all <- data[[1]]
data_tested <- data[[2]]
data_T1 <-data_tested[data_tested$trial_type=="T1",]
data_T1 <- center_data(data_T1)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthExperimental Set-Up Effects
## glm flew_b ~ chamber binomial data_T1
## AIC: 458.9055
## (Intercept) coeff: 0.2876821 Pr(>|t|): 0.5942526
## chamberA-2 coeff: -0.0966268 Pr(>|t|): 0.8766875
## chamberA-3 coeff: 0.117783 Pr(>|t|): 0.8474768
## chamberA-4 coeff: 0.4054651 Pr(>|t|): 0.5053931
## chamberB-1 coeff: 0.0918075 Pr(>|t|): 0.8875092
## chamberB-2 coeff: 0.2130932 Pr(>|t|): 0.7267933
## chamberB-3 coeff: -0.074108 Pr(>|t|): 0.9040261
## chamberB-4 coeff: 0.3254224 Pr(>|t|): 0.6114074## glm flew_b ~ days_from_start_c binomial data_T1
## AIC: 448.1122
## (Intercept) coeff: 0.4298329 Pr(>|t|): 0.0001338893 *
## days_from_start_c coeff: -0.0339142 Pr(>|t|): 0.2615186## glm flew_b ~ min_from_IncStart binomial data_T1
## AIC: 449.1128
## (Intercept) coeff: 0.360952 Pr(>|t|): 0.03534992 *
## min_from_IncStart coeff: 0.0004162 Pr(>|t|): 0.6064583Biological Effects
## glm flew_b ~ mass_c binomial data_T1
## AIC: 401.7318
## (Intercept) coeff: 0.4447371 Pr(>|t|): 0.0002377227 *
## mass_c coeff: -33.2650205 Pr(>|t|): 6.742316e-09 *Morphology Effects
## glm flew_b ~ beak_c binomial data_T1
## AIC: 438.1289
## (Intercept) coeff: 0.4391226 Pr(>|t|): 0.0001235556 *
## beak_c coeff: -0.368997 Pr(>|t|): 0.0009363177 *## glm flew_b ~ thorax_c binomial data_T1
## AIC: 441.4424
## (Intercept) coeff: 0.4374127 Pr(>|t|): 0.0001215489 *
## thorax_c coeff: -1.0468489 Pr(>|t|): 0.005435159 *## glm flew_b ~ body_c binomial data_T1
## AIC: 447.0143
## (Intercept) coeff: 0.4311034 Pr(>|t|): 0.000131644 *
## body_c coeff: -0.1613906 Pr(>|t|): 0.1253027## glm flew_b ~ wing_c binomial data_T1
## AIC: 449.2264
## (Intercept) coeff: 0.4283101 Pr(>|t|): 0.0001371065 *
## wing_c coeff: -0.053208 Pr(>|t|): 0.6957995Glmer() Testing Notes
(1|chamber) as a random effect does not need to be included in the Trial 1 glmer() analyses because in all cases the variance is essentially 0. What pops up if you do try glmer on the best fit model? A singular fit error is raised: “When you obtain a singular fit, this is often indicating that the model is overfitted – that is, the random effects structure is too complex to be supported by the data, which naturally leads to the advice to remove the most complex part of the random effects structure (usually random slopes). The benefit of this approach is that it leads to a more parsimonious (conservative) model that is not over-fitted.” Source: https://stats.stackexchange.com/questions/378939/dealing-with-singular-fit-in-mixed-models
However, (1|population) and (1|days_from_start_c) does matter when doing the morphology analyses in Trial 1, and will not throw an error because the variance is large enough to change the fit of a model.
Without Mass
# Remove any missing masses
data_T1 <- data_T1 %>%
filter(!is.na(mass))
data_T1 <- center_data(data_T1)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthR = data_T1$flew_b
A = data_T1$host_c
B = data_T1$sex_c
C = data_T1$sym_dist
D = data_T1$mass_c
data<-data.frame(R, A, B, C, D)
model_script = paste0(source_path,"generic models-binomial glm 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs 403.0853 405.0821 405.4319 405.5377 406.4848 406.5945
## models 8 11 2 15 4 6
## probs 0.3595574 0.1324877 0.1112301 0.1054958 0.06570071 0.06219381
##
## m8 glm(formula = R ~ A * B, family = binomial, data = data)
## m11 glm(formula = R ~ A * B + C, family = binomial, data = data)
## m2 glm(formula = R ~ B, family = binomial, data = data)
## m15 glm(formula = R ~ A * B + B * C, family = binomial, data = data)
## m4 glm(formula = R ~ A + B, family = binomial, data = data)
## m6 glm(formula = R ~ B + C, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m8, m11, test="Chisq") ## Adding C does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A * B
## Model 2: R ~ A * B + C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 325 395.09
## 2 324 395.08 1 0.0032319 0.9547anova(m11, m15, test="Chisq") ## Adding B*C interaction does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A * B + C
## Model 2: R ~ A * B + B * C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 324 395.08
## 2 323 393.54 1 1.5444 0.214anova(m11, m14, test="Chisq") ## Adding A*C interaction does not improve fit ## Analysis of Deviance Table
##
## Model 1: R ~ A * B + C
## Model 2: R ~ A * B + A * C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 324 395.08
## 2 323 395.08 1 0.0041801 0.9484nomass_T1<-glm(flew_b~sex_c*host_c, family=binomial, data=data_T1)
tidy_regression(nomass_T1, is_color=output_col)## glm flew_b ~ sex_c * host_c binomial data_T1
## AIC: 403.0853
## (Intercept) coeff: 0.2398549 Pr(>|t|): 0.07395882 .
## sex_c coeff: -0.6224552 Pr(>|t|): 3.532318e-06 *
## host_c coeff: -0.0561982 Pr(>|t|): 0.675461
## sex_c:host_c coeff: 0.3136354 Pr(>|t|): 0.01946434 *- strong negative effect of sex, where if you are female you are less likely to fly.
- No direct effect of host plant
- strong positive effect of sex and host where if female and from GRT, then more likely to fly.
With Mass
model_script = paste0(source_path,"generic models-binomial glm 4-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3]
## AICs 395.3859 395.7653 395.8829
## models 63 26 51
## probs 0.08435482 0.06977781 0.06579118
##
## m63 glm(formula = R ~ A * D + C * D + B, family = binomial, data = data)
## m26 glm(formula = R ~ A * D + B, family = binomial, data = data)
## m51 glm(formula = R ~ A * D + C * D, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0R = data_T1$flew_b
A = data_T1$host_c
B = data_T1$sex_c
C = data_T1$mass_c
data<-data.frame(R, A, B, C)
model_script = paste0(source_path,"generic models-binomial glm 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5]
## AICs 395.7653 397.0146 397.2769 397.5315 397.6811
## models 12 9 11 14 16
## probs 0.2751383 0.1473246 0.1292128 0.1137668 0.10557
##
## m12 glm(formula = R ~ A * C + B, family = binomial, data = data)
## m9 glm(formula = R ~ A * C, family = binomial, data = data)
## m11 glm(formula = R ~ A * B + C, family = binomial, data = data)
## m14 glm(formula = R ~ A * B + A * C, family = binomial, data = data)
## m16 glm(formula = R ~ A * C + B * C, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m9, m12, test="Chisq") # Adding B does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A * C
## Model 2: R ~ A * C + B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 325 389.01
## 2 324 385.77 1 3.2493 0.07146 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m51, m63, test="Chisq") # Adding B does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A * D + C * D
## Model 2: R ~ A * D + C * D + B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 323 383.88
## 2 322 381.39 1 2.4971 0.1141anova(m27, m51, test="Chisq") # Adding C*D does improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A * D + C
## Model 2: R ~ A * D + C * D
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 324 389.01
## 2 323 383.88 1 5.1241 0.0236 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# top model if don't include sym_dist
mass_T1_nosym <- glm(flew_b~host_c*mass_c, family=binomial, data=data_T1)
tidy_regression(mass_T1_nosym, is_color=output_col)## glm flew_b ~ host_c * mass_c binomial data_T1
## AIC: 397.0146
## (Intercept) coeff: 0.3972549 Pr(>|t|): 0.002534805 *
## host_c coeff: -0.1665608 Pr(>|t|): 0.2055597
## mass_c coeff: -28.0261972 Pr(>|t|): 4.143597e-06 *
## host_c:mass_c coeff: 15.6704317 Pr(>|t|): 0.01004477 *# top model if you do include sym_dist
mass_T1_sym<- glm(flew_b~host_c*mass_c + sym_dist*mass_c, family=binomial, data=data_T1)
tidy_regression(mass_T1_sym, is_color=output_col)## glm flew_b ~ host_c * mass_c + sym_dist * mass_c binomial data_T1
## AIC: 395.8829
## (Intercept) coeff: 0.4350653 Pr(>|t|): 0.04305314 *
## host_c coeff: -0.1640846 Pr(>|t|): 0.350775
## mass_c coeff: -15.1574096 Pr(>|t|): 0.07927326 .
## sym_dist coeff: -0.1083302 Pr(>|t|): 0.5155587
## host_c:mass_c coeff: 22.7816838 Pr(>|t|): 0.001430077 *
## mass_c:sym_dist coeff: -17.3520072 Pr(>|t|): 0.04407444 *Morphology
data_T1 <-data_tested[data_tested$trial_type=="T1",]
data_T1 <- data_T1 %>%
filter(!is.na(body))
data_T1 <- center_data(data_T1)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthR = data_T1$flew_b
A = data_T1$thorax_c
B = data_T1$body_c
C = data_T1$wing_c
data<-data.frame(R, A, B, C)
model_script = paste0(source_path,"generic models-binomial glm 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 420.4885 421.1253 421.4696 421.5025 421.8375 421.9204 422.082
## models 13 15 7 12 14 16 11
## probs 0.1864269 0.1355936 0.1141476 0.1122827 0.09496689 0.09111367 0.08404131
## [,8]
## AICs 422.6621
## models 17
## probs 0.06288207
##
## m13 glm(formula = R ~ B * C + A, family = binomial, data = data)
## m15 glm(formula = R ~ A * B + B * C, family = binomial, data = data)
## m7 glm(formula = R ~ A + B + C, family = binomial, data = data)
## m12 glm(formula = R ~ A * C + B, family = binomial, data = data)
## m14 glm(formula = R ~ A * B + A * C, family = binomial, data = data)
## m16 glm(formula = R ~ A * C + B * C, family = binomial, data = data)
## m11 glm(formula = R ~ A * B + C, family = binomial, data = data)
## m17 glm(formula = R ~ A * B + A * C + B * C, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m7, m13, test="Chisq") # Adding B*C marginally improves fit## Analysis of Deviance Table
##
## Model 1: R ~ A + B + C
## Model 2: R ~ B * C + A
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 328 413.47
## 2 327 410.49 1 2.9811 0.08424 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m13, m16, test="Chisq") # Adding A*C does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ B * C + A
## Model 2: R ~ A * C + B * C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 327 410.49
## 2 326 409.92 1 0.56814 0.451morph_T1 <- glm(flew_b ~ body_c * wing_c + thorax_c, family = binomial, data = data_T1)
tidy_regression(morph_T1, is_color=output_col)## glm flew_b ~ body_c * wing_c + thorax_c binomial data_T1
## AIC: 420.4885
## (Intercept) coeff: 0.6086248 Pr(>|t|): 2.775207e-05 *
## body_c coeff: -1.1762506 Pr(>|t|): 0.06266837 .
## wing_c coeff: 2.350309 Pr(>|t|): 0.0004808155 *
## thorax_c coeff: -2.6772762 Pr(>|t|): 0.01747959 *
## body_c:wing_c coeff: -0.1698641 Pr(>|t|): 0.08756481 .maringal negative effect of body wher ethe larger the body the less likely to fly
strong positive effect of wing length where the longer the wing, the more likely to fly
storng negative effect of thorax where the larger the thorasxs, the less likely to fly
no effect of body*wing (we know body and wing are closely related so they could be canceling each other out (b/c body is a proxi for mass))
morph_T1glmer <- glmer(formula = flew_b ~ body_c * wing_c + thorax_c + (1| population), family = binomial, data = data_T1)
tidy_regression(morph_T1glmer, is_color=output_col) # did change the coefficients slightly, but not a better fit## glmer flew_b ~ body_c * wing_c + thorax_c + (1 | population) data_T1 binomial
## AIC: 421.8532 421.8532
## (Intercept) coeff: 0.5335423 Pr(>|t|): 0.005217949 *
## body_c coeff: -1.2372535 Pr(>|t|): 0.05374932 .
## wing_c coeff: 2.4704559 Pr(>|t|): 0.0003676219 *
## thorax_c coeff: -2.792138 Pr(>|t|): 0.01517357 *
## body_c:wing_c coeff: -0.1651854 Pr(>|t|): 0.1002469R = data_T1$flew_b
A = data_T1$thorax_c
B = data_T1$wing2body_c
data<-data.frame(R, A, B)
model_script = paste0(source_path,"generic models-binomial glm 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2]
## AICs 419.9404 420.2258
## models 4 3
## probs 0.531513 0.4608213
##
## m4 glm(formula = R ~ A * B, family = binomial, data = data)
## m3 glm(formula = R ~ A + B, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m3, m4, test="Chisq") # Adding A*B does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A + B
## Model 2: R ~ A * B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 329 414.23
## 2 328 411.94 1 2.2854 0.1306morph_T1_ratio <- glm(flew_b ~ wing2body_c + thorax_c, family = binomial, data = data_T1)
tidy_regression(morph_T1_ratio, is_color=output_col)## glm flew_b ~ wing2body_c + thorax_c binomial data_T1
## AIC: 420.2258
## (Intercept) coeff: 0.4670867 Pr(>|t|): 8.08552e-05 *
## wing2body_c coeff: 32.6671543 Pr(>|t|): 6.925235e-06 *
## thorax_c coeff: -1.2219033 Pr(>|t|): 0.001695651 *Summary
tidy_regression(nomass_T1, is_color=output_col)## glm flew_b ~ sex_c * host_c binomial data_T1
## AIC: 403.0853
## (Intercept) coeff: 0.2398549 Pr(>|t|): 0.07395882 .
## sex_c coeff: -0.6224552 Pr(>|t|): 3.532318e-06 *
## host_c coeff: -0.0561982 Pr(>|t|): 0.675461
## sex_c:host_c coeff: 0.3136354 Pr(>|t|): 0.01946434 *tidy_regression(mass_T1_sym, is_color=output_col)## glm flew_b ~ host_c * mass_c + sym_dist * mass_c binomial data_T1
## AIC: 395.8829
## (Intercept) coeff: 0.4350653 Pr(>|t|): 0.04305314 *
## host_c coeff: -0.1640846 Pr(>|t|): 0.350775
## mass_c coeff: -15.1574096 Pr(>|t|): 0.07927326 .
## sym_dist coeff: -0.1083302 Pr(>|t|): 0.5155587
## host_c:mass_c coeff: 22.7816838 Pr(>|t|): 0.001430077 *
## mass_c:sym_dist coeff: -17.3520072 Pr(>|t|): 0.04407444 *tidy_regression(mass_T1_nosym, is_color=output_col)## glm flew_b ~ host_c * mass_c binomial data_T1
## AIC: 397.0146
## (Intercept) coeff: 0.3972549 Pr(>|t|): 0.002534805 *
## host_c coeff: -0.1665608 Pr(>|t|): 0.2055597
## mass_c coeff: -28.0261972 Pr(>|t|): 4.143597e-06 *
## host_c:mass_c coeff: 15.6704317 Pr(>|t|): 0.01004477 *tidy_regression(morph_T1, is_color=output_col) # glmer did not improve fit## glm flew_b ~ body_c * wing_c + thorax_c binomial data_T1
## AIC: 420.4885
## (Intercept) coeff: 0.6086248 Pr(>|t|): 2.775207e-05 *
## body_c coeff: -1.1762506 Pr(>|t|): 0.06266837 .
## wing_c coeff: 2.350309 Pr(>|t|): 0.0004808155 *
## thorax_c coeff: -2.6772762 Pr(>|t|): 0.01747959 *
## body_c:wing_c coeff: -0.1698641 Pr(>|t|): 0.08756481 .tidy_regression(morph_T1_ratio, is_color=output_col)## glm flew_b ~ wing2body_c + thorax_c binomial data_T1
## AIC: 420.2258
## (Intercept) coeff: 0.4670867 Pr(>|t|): 8.08552e-05 *
## wing2body_c coeff: 32.6671543 Pr(>|t|): 6.925235e-06 *
## thorax_c coeff: -1.2219033 Pr(>|t|): 0.001695651 *Trial 2
Testing glmer() with pop and chamber
Cleaning the Data
data <- read_flight_data("data/all_flight_data-Winter2020.csv")## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object length
## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_all <- data[[1]]
data_tested <- data[[2]]
data_T2 <-data_tested[data_tested$trial_type=="T2",]
data_T2 <- center_data(data_T2)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthExperimental Set-Up Effects
## glm flew_b ~ chamber binomial data_T2
## AIC: 393.888
## (Intercept) coeff: -0.0909718 Pr(>|t|): 0.7631039
## chamberA-2 coeff: 0.3273606 Pr(>|t|): 0.4754205
## chamberA-3 coeff: 0.784119 Pr(>|t|): 0.1224819
## chamberA-4 coeff: -0.2837217 Pr(>|t|): 0.488549
## chamberB-1 coeff: 0.784119 Pr(>|t|): 0.1559208
## chamberB-2 coeff: 0.6017974 Pr(>|t|): 0.2039747
## chamberB-3 coeff: -0.468644 Pr(>|t|): 0.3199648
## chamberB-4 coeff: -0.1809619 Pr(>|t|): 0.6866405## glm flew_b ~ days_from_start_c binomial data_T2
## AIC: 393.3056
## (Intercept) coeff: 0.014415 Pr(>|t|): 0.9039395
## days_from_start_c coeff: -0.0540142 Pr(>|t|): 0.2070645## glm flew_b ~ min_from_IncStart binomial data_T2
## AIC: 394.6235
## (Intercept) coeff: -0.0779168 Pr(>|t|): 0.7062456
## min_from_IncStart coeff: 0.0004647 Pr(>|t|): 0.5857749Biological Effects
## glm flew_b ~ mass_c binomial data_T2
## AIC: 358.9575
## (Intercept) coeff: -0.0288154 Pr(>|t|): 0.8227357
## mass_c coeff: -33.6659527 Pr(>|t|): 1.030843e-07 *Morphology Effects
## glm flew_b ~ beak_c binomial data_T2
## AIC: 380.4646
## (Intercept) coeff: 0.0069373 Pr(>|t|): 0.9547796
## beak_c coeff: -0.468637 Pr(>|t|): 0.0002555501 *## glm flew_b ~ thorax_c binomial data_T2
## AIC: 381.7558
## (Intercept) coeff: 0.0109082 Pr(>|t|): 0.9287235
## thorax_c coeff: -1.49497 Pr(>|t|): 0.0004698825 *## glm flew_b ~ body_c binomial data_T2
## AIC: 387.3314
## (Intercept) coeff: 0.0133229 Pr(>|t|): 0.9121223
## body_c coeff: -0.3168907 Pr(>|t|): 0.006959575 *## glm flew_b ~ wing_c binomial data_T2
## AIC: 392.1204
## (Intercept) coeff: 0.014232 Pr(>|t|): 0.905352
## wing_c coeff: -0.2472874 Pr(>|t|): 0.09699422 .Glmer() Testing Notes
(1|population) + (1|chamber) as random effects do not need to be included in the Trial 2 glmer() mass and no mass analyses because in all cases the variance is essentially 0. What pops up if you do try glmer on the best fit model? A singular fit error is raised. Refer to the reasoning in Trial 1.
However, no error is thrown when included in the morphology analyses.
Without Mass
R = data_T2$flew_b
A = data_T2$host_c
B = data_T2$sex_c
C = data_T2$sym_dist
D = data_T2$mass_c # as a covariate
data<-data.frame(R, A, B, C, D)
model_script = paste0(source_path,"generic models-binomial glm 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs 361.1644 362.5404 363.037 364.4845 364.5336 364.8157
## models 2 4 6 7 8 10
## probs 0.3641586 0.1830186 0.1427814 0.06923818 0.06755975 0.05867098
##
## m2 glm(formula = R ~ B, family = binomial, data = data)
## m4 glm(formula = R ~ A + B, family = binomial, data = data)
## m6 glm(formula = R ~ B + C, family = binomial, data = data)
## m7 glm(formula = R ~ A + B + C, family = binomial, data = data)
## m8 glm(formula = R ~ A * B, family = binomial, data = data)
## m10 glm(formula = R ~ B * C, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m2, m4, test="Chisq") # Adding A does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ B
## Model 2: R ~ A + B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 280 357.16
## 2 279 356.54 1 0.624 0.4296anova(m2, m6, test="Chisq") # Adding C does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ B
## Model 2: R ~ B + C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 280 357.16
## 2 279 357.04 1 0.12745 0.7211nomass_T2 <-glm(flew_b~sex_c, family=binomial, data=data_T2)
tidy_regression(nomass_T2, is_color=output_col)## glm flew_b ~ sex_c binomial data_T2
## AIC: 361.1644
## (Intercept) coeff: -0.2425498 Pr(>|t|): 0.07779975 .
## sex_c coeff: -0.7620335 Pr(>|t|): 3.010899e-08 *- sex is the only significant and strong effect, so that if you are female you are much less likely to disperse/fly
With Mass
model_script = paste0(source_path,"generic models-binomial glm 4-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4]
## AICs 358.8454 358.9575 359.2153 359.563
## models 9 4 7 12
## probs 0.0889176 0.08406912 0.07390225 0.06210914
##
## m9 glm(formula = R ~ B + D, family = binomial, data = data)
## m4 glm(formula = R ~ D, family = binomial, data = data)
## m7 glm(formula = R ~ A + D, family = binomial, data = data)
## m12 glm(formula = R ~ A + B + D, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m4, m7, test="Chisq") # Adding A does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ D
## Model 2: R ~ A + D
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 280 354.96
## 2 279 353.22 1 1.7422 0.1869anova(m7, m12, test="Chisq") # Adding B does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A + D
## Model 2: R ~ A + B + D
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 279 353.22
## 2 278 351.56 1 1.6523 0.1986anova(m4, m9, test="Chisq") # Adding B does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ D
## Model 2: R ~ B + D
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 280 354.96
## 2 279 352.85 1 2.1121 0.1461mass_T2 <-glm(flew_b~mass_c, family=binomial, data=data_T2)
tidy_regression(mass_T2, is_color=output_col)## glm flew_b ~ mass_c binomial data_T2
## AIC: 358.9575
## (Intercept) coeff: -0.0288154 Pr(>|t|): 0.8227357
## mass_c coeff: -33.6659527 Pr(>|t|): 1.030843e-07 *- mass is the only significant and extremely strong effect, so that if you are heavy you are much less likely to disperse/fly
Mass is really dominating everything, so let’s split by sex after looking at Trial 2.
Morphology
# Remove any missing masses and morphology measurements
data_T2 <- data_T2 %>%
filter(!is.na(mass)) %>%
filter(!is.na(body))
data_T2 <- center_data(data_T2)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthR = data_T2$flew_b
A = data_T2$thorax_c
B = data_T2$body_c
C = data_T2$wing_c
data<-data.frame(R, A, B, C)
model_script = paste0(source_path,"generic models-binomial glm 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4]
## AICs 364.316 365.2947 366.2914 368.3788
## models 16 15 17 13
## probs 0.4230217 0.2593179 0.1575417 0.05547893
##
## m16 glm(formula = R ~ A * C + B * C, family = binomial, data = data)
## m15 glm(formula = R ~ A * B + B * C, family = binomial, data = data)
## m17 glm(formula = R ~ A * B + A * C + B * C, family = binomial, data = data)
## m13 glm(formula = R ~ B * C + A, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m13, m16, test="Chisq") # Adding A*C improves fit## Analysis of Deviance Table
##
## Model 1: R ~ B * C + A
## Model 2: R ~ A * C + B * C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 277 358.38
## 2 276 352.32 1 6.0628 0.01381 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1morph_T2 <- glm(flew_b ~ thorax_c * wing_c + body_c * wing_c, family = binomial, data = data_T2)
tidy_regression(morph_T2, is_color=output_col)## glm flew_b ~ thorax_c * wing_c + body_c * wing_c binomial data_T2
## AIC: 364.316
## (Intercept) coeff: 0.1914188 Pr(>|t|): 0.2140046
## thorax_c coeff: -2.3082152 Pr(>|t|): 0.0828982 .
## wing_c coeff: 1.9063552 Pr(>|t|): 0.008518188 *
## body_c coeff: -1.1472691 Pr(>|t|): 0.09990423 .
## thorax_c:wing_c coeff: 3.692847 Pr(>|t|): 0.0199121 *
## wing_c:body_c coeff: -1.216138 Pr(>|t|): 0.006636944 *negative marginal effect of thorax, where the wider the thorax, the then bug less likely to fly
positive effect of wing length, where if wing longer then bug more likely to fly
negative effect of body, where the longer the body the less likely the bug is to fly
positive effect of thorax:wing where the wider the thorax and wing the more likely the bug is to fly (this encapsulates overall flight power in a way b/c thorax muscles and wing length)
negative effect of wing and body where the longer the wing and body the less likely the bug is to fly (this encapsulates mass limitations)
R = data_T2$flew_b
A = data_T2$thorax_c
B = data_T2$wing2body_c
data<-data.frame(R, A, B)
model_script = paste0(source_path,"generic models-binomial glm 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2]
## AICs 367.6219 368.6156
## models 3 4
## probs 0.6206869 0.3776612
##
## m3 glm(formula = R ~ A + B, family = binomial, data = data)
## m4 glm(formula = R ~ A * B, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m3, m4, test="Chisq") # Adding A*B does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A + B
## Model 2: R ~ A * B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 279 361.62
## 2 278 360.62 1 1.0063 0.3158morph_T2_ratio <- glm(flew_b ~ thorax_c + wing2body_c, family = binomial, data = data_T2)
tidy_regression(morph_T2_ratio, is_color=output_col)## glm flew_b ~ thorax_c + wing2body_c binomial data_T2
## AIC: 367.6219
## (Intercept) coeff: 0.0109675 Pr(>|t|): 0.9303663
## thorax_c coeff: -1.6132431 Pr(>|t|): 0.0002366193 *
## wing2body_c coeff: 28.5706239 Pr(>|t|): 0.0001653885 *Summary
tidy_regression(nomass_T2, is_color=output_col)## glm flew_b ~ sex_c binomial data_T2
## AIC: 361.1644
## (Intercept) coeff: -0.2425498 Pr(>|t|): 0.07779975 .
## sex_c coeff: -0.7620335 Pr(>|t|): 3.010899e-08 *tidy_regression(mass_T2, is_color=output_col)## glm flew_b ~ mass_c binomial data_T2
## AIC: 358.9575
## (Intercept) coeff: -0.0288154 Pr(>|t|): 0.8227357
## mass_c coeff: -33.6659527 Pr(>|t|): 1.030843e-07 *tidy_regression(morph_T2, is_color=output_col) # glmer did not improve fit## glm flew_b ~ thorax_c * wing_c + body_c * wing_c binomial data_T2
## AIC: 364.316
## (Intercept) coeff: 0.1914188 Pr(>|t|): 0.2140046
## thorax_c coeff: -2.3082152 Pr(>|t|): 0.0828982 .
## wing_c coeff: 1.9063552 Pr(>|t|): 0.008518188 *
## body_c coeff: -1.1472691 Pr(>|t|): 0.09990423 .
## thorax_c:wing_c coeff: 3.692847 Pr(>|t|): 0.0199121 *
## wing_c:body_c coeff: -1.216138 Pr(>|t|): 0.006636944 *tidy_regression(morph_T2_ratio, is_color=output_col)## glm flew_b ~ thorax_c + wing2body_c binomial data_T2
## AIC: 367.6219
## (Intercept) coeff: 0.0109675 Pr(>|t|): 0.9303663
## thorax_c coeff: -1.6132431 Pr(>|t|): 0.0002366193 *
## wing2body_c coeff: 28.5706239 Pr(>|t|): 0.0001653885 *Females
Cleaning the Data
data <- read_flight_data("data/all_flight_data-Winter2020.csv")## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object length
## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_all <- data[[1]]
data_tested <- data[[2]]
d = create_delta_data(data_tested)data_fem <- d[d$sex=="F",]
data_fem <- center_data(data_fem, is_not_unique_data = FALSE)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_fem## # A tibble: 96 x 86
## # Groups: ID, sex, population, site, host_plant, latitude, longitude,
## # total_eggs, beak, thorax, wing, body, w_morph, morph_notes, tested, host_c,
## # sex_c, w_morph_c, lat_c, sym_dist, sym_dist_s, total_eggs_c, beak_c,
## # thorax_c, thorax_s, body_c, wing_c, wing2body, wing2body_c [96]
## ID sex population site host_plant latitude longitude total_eggs beak
## <fct> <fct> <fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 5 F Plantation… Arego… C. corind… 25.0 -80.6 209 7.47
## 2 7 F Plantation… Found… C. corind… 25.0 -80.6 8 8.19
## 3 12 F Key Largo KLMRL C. corind… 25.1 -80.4 329 8.79
## 4 13 F Key Largo KLMRL C. corind… 25.1 -80.4 296 6.99
## 5 19 F Key Largo JP C. corind… 25.1 -80.4 NA 7.54
## 6 20 F Key Largo JP C. corind… 25.1 -80.4 55 8.73
## 7 25 F North Key … DD fr… C. corind… 25.3 -80.3 NA 8.03
## 8 26 F North Key … Charl… C. corind… 25.2 -80.3 33 7.15
## 9 27 F North Key … Charl… C. corind… 25.2 -80.3 158 7.58
## 10 29 F North Key … Charl… C. corind… 25.2 -80.3 67 6.24
## # … with 86 more rows, and 77 more variables: thorax <dbl>, wing <dbl>,
## # body <dbl>, w_morph <fct>, morph_notes <fct>, tested <fct>, host_c <dbl>,
## # sex_c <dbl>, w_morph_c <dbl>, lat_c <dbl>, sym_dist <dbl>,
## # sym_dist_s <dbl>, total_eggs_c <dbl>, beak_c <dbl>, thorax_c <dbl>,
## # thorax_s <dbl>, body_c <dbl>, wing_c <dbl>, wing2body <dbl>,
## # wing2body_c <dbl>, wing2body_s <dbl>, trial_type <list>, chamber <list>,
## # average_speed <list>, total_flight_time <list>, distance <list>,
## # shortest_flying_bout <list>, longest_flying_bout <list>,
## # recording_duration <list>, max_speed <list>, test_date <list>,
## # time_start <list>, time_end <list>, flew <list>, flight_type <list>,
## # mass <list>, EWM <list>, flew_b <list>, host_temp <list>, eggs_b <list>,
## # minute_duration <list>, minute_duration_c <list>, min_from_IncStart <dbl>,
## # min_from_IncStart_c <dbl>, min_from_IncStart_s <dbl>,
## # days_from_start <list>, days_from_start_c <list>, mass_c <list>,
## # mass_s <list>, avg_mass <dbl>, avg_mass_c <dbl>, avg_mass_s <dbl>,
## # trial_type_og <list>, dt_start_c <list>, dt_end_c <list>,
## # datetime_start <list>, datetime_end <list>, hr_start <list>, hr_end <list>,
## # num_flew <dbl>, num_notflew <dbl>, avg_days <dbl>, num_egg <dbl>,
## # egg_diff <dbl>, mass_diff <dbl>, flew_diff <dbl>, dist_diff <dbl>,
## # speed_diff <dbl>, mass_per <dbl>, avg_rec_dur <dbl>, avg_time_start <dbl>,
## # flight_case <fct>, egg_case <dbl>, flew_prob <dbl>, avg_days_c <dbl>,
## # avg_mass_logsqrt <dbl>, wing2body_logsqrt_i <dbl>Experimental, Biological, and Morphological Effects
## glm cbind(num_flew, num_notflew) ~ avg_mass binomial data_fem
## AIC: 219.2622
## (Intercept) coeff: 1.5788786 Pr(>|t|): 0.01741761 *
## avg_mass coeff: -27.2115254 Pr(>|t|): 0.001204519 *## glm cbind(num_flew, num_notflew) ~ total_eggs binomial data_fem
## AIC: 163.3111
## (Intercept) coeff: 0.2982785 Pr(>|t|): 0.2614142
## total_eggs coeff: -0.0128174 Pr(>|t|): 7.35417e-06 *## glm cbind(num_flew, num_notflew) ~ beak_c binomial data_fem
## AIC: 220.9202
## (Intercept) coeff: -0.5860233 Pr(>|t|): 0.0001583204 *
## beak_c coeff: 0.5964896 Pr(>|t|): 0.002221452 *## glm cbind(num_flew, num_notflew) ~ thorax_c binomial data_fem
## AIC: 229.0582
## (Intercept) coeff: -0.5611857 Pr(>|t|): 0.0002001936 *
## thorax_c coeff: 0.7483512 Pr(>|t|): 0.1706597## glm cbind(num_flew, num_notflew) ~ body_c binomial data_fem
## AIC: 225.1188
## (Intercept) coeff: -0.5753868 Pr(>|t|): 0.0001721078 *
## body_c coeff: 0.3628557 Pr(>|t|): 0.01896725 *## glm cbind(num_flew, num_notflew) ~ wing_c binomial data_fem
## AIC: 224.6506
## (Intercept) coeff: -0.5778662 Pr(>|t|): 0.0001669733 *
## wing_c coeff: 0.4469276 Pr(>|t|): 0.01544936 *Without Mass
R1 = data_fem$num_flew
R2 = data_fem$num_notflew
A = data_fem$host_c
B = data_fem$wing_c
C = data_fem$avg_mass
X = data_fem$population # can't do trial type anymore
data<-data.frame(R1, R2, A, B, C, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4]
## AICs 226.6506 227.6224 228.5451 230.9545
## models 2 4 3 5
## probs 0.4607286 0.2834078 0.1786729 0.05356317
##
## m2 cbind(R1, R2) ~ B + (1 | X)
## m4 cbind(R1, R2) ~ A * B + (1 | X)
## m3 cbind(R1, R2) ~ A + B + (1 | X)
## m0 cbind(R1, R2) ~ 1 + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m0, m2, test="Chisq") # Adding B improves fit## Data: data
## Models:
## m0: cbind(R1, R2) ~ 1 + (1 | X)
## m2: cbind(R1, R2) ~ B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 2 230.95 236.08 -113.48 226.95
## m2 3 226.65 234.34 -110.33 220.65 6.3039 1 0.01205 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2, m3, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m2: cbind(R1, R2) ~ B + (1 | X)
## m3: cbind(R1, R2) ~ A + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 3 226.65 234.34 -110.33 220.65
## m3 4 228.54 238.80 -110.27 220.54 0.1055 1 0.7453nomass_fem <- glmer(cbind(num_flew, num_notflew) ~ wing_c + (1|population), data=data_fem, family=binomial)## boundary (singular) fit: see ?isSingulartidy_regression(nomass_fem, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing_c + (1 | population) data_fem binomial
## AIC: 226.6506 226.6506
## (Intercept) coeff: -0.5778662 Pr(>|t|): 0.0001669843 *
## wing_c coeff: 0.4469276 Pr(>|t|): 0.01545033 *- positive effect of wing length where the longer the wing the more times the bug will fly
With Mass
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 203.03 203.7663 203.9601 204.6943 204.8636 205.4558 205.7597
## models 12 11 14 6 16 15 17
## probs 0.2396203 0.1658205 0.1505049 0.1042629 0.09580079 0.07124868 0.06120348
##
## m12 cbind(R1, R2) ~ A * C + B + (1 | X)
## m11 cbind(R1, R2) ~ A * B + C + (1 | X)
## m14 cbind(R1, R2) ~ A * B + A * C + (1 | X)
## m6 cbind(R1, R2) ~ B + C + (1 | X)
## m16 cbind(R1, R2) ~ A * C + B * C + (1 | X)
## m15 cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m17 cbind(R1, R2) ~ 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: 0anova(m7, m12, test="Chisq") # Adding A*C improves fit## Data: data
## Models:
## m7: cbind(R1, R2) ~ A + B + C + (1 | X)
## m12: cbind(R1, R2) ~ A * C + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m7 5 206.24 219.06 -98.121 196.24
## m12 6 203.03 218.42 -95.515 191.03 5.2121 1 0.02243 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m12, m14, test="Chisq") # Adding A*B does not improve fit## Data: data
## Models:
## m12: cbind(R1, R2) ~ A * C + B + (1 | X)
## m14: cbind(R1, R2) ~ A * B + A * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m12 6 203.03 218.42 -95.515 191.03
## m14 7 203.96 221.91 -94.980 189.96 1.0699 1 0.301mass_fem <- glmer(cbind(num_flew, num_notflew) ~ host_c * avg_mass + wing_c + (1|population), data=data_fem, family=binomial)## boundary (singular) fit: see ?isSingulartidy_regression(mass_fem, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ host_c * avg_mass + wing_c + (1 | population) data_fem binomial
## AIC: 203.03 203.03
## (Intercept) coeff: 1.9138191 Pr(>|t|): 0.03019205 *
## host_c coeff: -1.9373234 Pr(>|t|): 0.01803064 *
## avg_mass coeff: -31.826984 Pr(>|t|): 0.005588096 *
## wing_c coeff: 0.9189461 Pr(>|t|): 7.071493e-05 *
## host_c:avg_mass coeff: 24.0687046 Pr(>|t|): 0.02131547 *negative effect of host where if from GRT the less times a bug will fly
strong negative effect of average mass where the heavier the less times a bug will fly
positive effect of wing length where the longer the wing the more times the bug will fly
strong positive effect of host*average mass where the heavier the bug and from GRT the more times the bug will fly (I’ve seen these strong conflicting interactions a lot with host and something else)
# check for collinearity between mass and wing length
data<-data.frame(data_fem$host_c, data_fem$wing_c, data_fem$avg_mass )
colnames(data) <- c("Host Plant", "Wing Length", "Average Mass")
pairs(data)
GRT bugs weigh less.
Eggs
R1 = data_fem$num_flew
R2 = data_fem$num_notflew
A = data_fem$host_c
B = data_fem$wing_c
C = data_fem$avg_mass
D = data_fem$total_eggs_c # can't use eggs_b anymore
X = data_fem$population
data<-data.frame(R1, R2, A, B, C, D, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1-RF + 4-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2]
## AICs 159.6743 159.792
## models 9 20
## probs 0.09002536 0.08487926
##
## m9 cbind(R1, R2) ~ B + D + (1 | X)
## m20 cbind(R1, R2) ~ B * D + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 83anova(m9, m20, test="Chisq") # Adding B*D does not improve fit## Data: data
## Models:
## m9: cbind(R1, R2) ~ B + D + (1 | X)
## m20: cbind(R1, R2) ~ B * D + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m9 4 159.67 169.54 -75.837 151.67
## m20 5 159.79 172.12 -74.896 149.79 1.8823 1 0.1701egg_model <- glmer(cbind(num_flew, num_notflew) ~ wing_c + total_eggs_c + (1|population), data=data_fem, family=binomial)## boundary (singular) fit: see ?isSingulartidy_regression(egg_model, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing_c + total_eggs_c + (1 | population) data_fem binomial
## AIC: 159.6743 159.6743
## (Intercept) coeff: -1.1574516 Pr(>|t|): 1.796822e-07 *
## wing_c coeff: 0.5954488 Pr(>|t|): 0.00931373 *
## total_eggs_c coeff: -0.0125051 Pr(>|t|): 9.730972e-06 *- positive effect of wing length
- negative effect of total eggs
Morphology
d <- data_fem %>%
filter(!is.na(body))
d$thorax_c <- d$thorax - mean(d$thorax)
d$wing_c <- d$wing - mean(d$wing)
d$body_c <- d$body - mean(d$body)
R1 = d$num_flew
R2 = d$num_notflew
A = d$thorax_c
B = d$body_c
C = d$wing_c
X = d$population
data<-data.frame(R1, R2, A, B, C, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs 224.2721 224.8072 225.6219 226.0386 226.2219 226.6506
## models 16 15 4 5 17 3
## probs 0.1927991 0.1475452 0.09817542 0.07971243 0.07273028 0.05870019
## [,7] [,8]
## AICs 226.9197 226.9665
## models 10 13
## probs 0.05130975 0.05012388
##
## m16 cbind(R1, R2) ~ A * C + B * C + (1 | X)
## m15 cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m4 cbind(R1, R2) ~ A + B + (1 | X)
## m5 cbind(R1, R2) ~ A + C + (1 | X)
## m17 cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)
## m3 cbind(R1, R2) ~ C + (1 | X)
## m10 cbind(R1, R2) ~ B * C + (1 | X)
## m13 cbind(R1, R2) ~ B * C + A + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m15, m16, test="Chisq") # replacing A*B with A*C improves fit## Data: data
## Models:
## m15: cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m16: cbind(R1, R2) ~ A * C + B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m15 7 224.81 242.76 -105.40 210.81
## m16 7 224.27 242.22 -105.14 210.27 0.535 0anova(m13, m16, test="Chisq") # Adding B*C improves fit## Data: data
## Models:
## m13: cbind(R1, R2) ~ B * C + A + (1 | X)
## m16: cbind(R1, R2) ~ A * C + B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m13 6 226.97 242.35 -107.48 214.97
## m16 7 224.27 242.22 -105.14 210.27 4.6943 1 0.03026 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1morph_fem <- glmer(cbind(num_flew, num_notflew) ~ thorax_c * wing_c + body_c * wing_c + (1 | population), data=d, family=binomial)## boundary (singular) fit: see ?isSingulartidy_regression(morph_fem, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ thorax_c * wing_c + body_c * wing_c + (1 | population) d binomial
## AIC: 224.2721 224.2721
## (Intercept) coeff: -0.3931431 Pr(>|t|): 0.05371586 .
## thorax_c coeff: -3.7184292 Pr(>|t|): 0.03426514 *
## wing_c coeff: 0.2200509 Pr(>|t|): 0.8215141
## body_c coeff: 1.2013836 Pr(>|t|): 0.2274205
## thorax_c:wing_c coeff: 4.418756 Pr(>|t|): 0.04558103 *
## wing_c:body_c coeff: -1.5087329 Pr(>|t|): 0.0272959 *strong negative effect of thorax where the wider the thorax, the less times the bug flies
no effect of wing length
no effect of body length
strong positive (marginal) effect of thorax*body where the longer the body and wider the thorax, the more times a bug flies (seems to conflict the single effects)
negative effect of of body*wing where the longer the body and wing, the less times the bug flies
R1 = d$num_flew
R2 = d$num_notflew
A = d$thorax_c
B = d$wing2body_c
X = d$population
data<-data.frame(R1, R2, A, B, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5]
## AICs 227.9456 228.9807 229.7947 230.9545 231.0582
## models 2 4 3 5 1
## probs 0.4122398 0.2456901 0.163544 0.09157709 0.08694903
##
## m2 cbind(R1, R2) ~ B + (1 | X)
## m4 cbind(R1, R2) ~ A * B + (1 | X)
## m3 cbind(R1, R2) ~ A + B + (1 | X)
## m0 cbind(R1, R2) ~ 1 + (1 | X)
## m1 cbind(R1, R2) ~ A + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m0, m2, test="Chisq") # Adding B marginally improves fit## Data: data
## Models:
## m0: cbind(R1, R2) ~ 1 + (1 | X)
## m2: cbind(R1, R2) ~ B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 2 230.95 236.08 -113.48 226.95
## m2 3 227.95 235.64 -110.97 221.95 5.0088 1 0.02522 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2, m3, test="Chisq") # Adding B does not improve fit## Data: data
## Models:
## m2: cbind(R1, R2) ~ B + (1 | X)
## m3: cbind(R1, R2) ~ A + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 3 227.95 235.64 -110.97 221.95
## m3 4 229.79 240.05 -110.90 221.79 0.151 1 0.6976morph_fem2 <- glmer(cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population), data=d, family=binomial)## boundary (singular) fit: see ?isSingulartidy_regression(morph_fem2, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population) d binomial
## AIC: 227.9456 227.9456
## (Intercept) coeff: -0.57421 Pr(>|t|): 0.0001722135 *
## wing2body_c coeff: 19.6346273 Pr(>|t|): 0.03293398 *- marginal strong positive effect of wing:body where the longer the body and heavier the body, the more times a female bug flew. (hm)
Summary
tidy_regression(nomass_fem, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing_c + (1 | population) data_fem binomial
## AIC: 226.6506 226.6506
## (Intercept) coeff: -0.5778662 Pr(>|t|): 0.0001669843 *
## wing_c coeff: 0.4469276 Pr(>|t|): 0.01545033 *tidy_regression(mass_fem, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ host_c * avg_mass + wing_c + (1 | population) data_fem binomial
## AIC: 203.03 203.03
## (Intercept) coeff: 1.9138191 Pr(>|t|): 0.03019205 *
## host_c coeff: -1.9373234 Pr(>|t|): 0.01803064 *
## avg_mass coeff: -31.826984 Pr(>|t|): 0.005588096 *
## wing_c coeff: 0.9189461 Pr(>|t|): 7.071493e-05 *
## host_c:avg_mass coeff: 24.0687046 Pr(>|t|): 0.02131547 *tidy_regression(egg_model, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing_c + total_eggs_c + (1 | population) data_fem binomial
## AIC: 159.6743 159.6743
## (Intercept) coeff: -1.1574516 Pr(>|t|): 1.796822e-07 *
## wing_c coeff: 0.5954488 Pr(>|t|): 0.00931373 *
## total_eggs_c coeff: -0.0125051 Pr(>|t|): 9.730972e-06 *tidy_regression(morph_fem, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ thorax_c * wing_c + body_c * wing_c + (1 | population) d binomial
## AIC: 224.2721 224.2721
## (Intercept) coeff: -0.3931431 Pr(>|t|): 0.05371586 .
## thorax_c coeff: -3.7184292 Pr(>|t|): 0.03426514 *
## wing_c coeff: 0.2200509 Pr(>|t|): 0.8215141
## body_c coeff: 1.2013836 Pr(>|t|): 0.2274205
## thorax_c:wing_c coeff: 4.418756 Pr(>|t|): 0.04558103 *
## wing_c:body_c coeff: -1.5087329 Pr(>|t|): 0.0272959 *tidy_regression(morph_fem2, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population) d binomial
## AIC: 227.9456 227.9456
## (Intercept) coeff: -0.57421 Pr(>|t|): 0.0001722135 *
## wing2body_c coeff: 19.6346273 Pr(>|t|): 0.03293398 *Females - Old Method
Cleaning the Data
data <- read_flight_data("data/all_flight_data-Winter2020.csv")## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object length
## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_all <- data[[1]]
data_tested <- data[[2]]
data_fem <- data_tested[data_tested$sex=="F",]
data_fem <- center_data(data_fem)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthExperimental Set-Up Effects
## glm flew_b ~ chamber binomial data_fem
## AIC: 283.3421
## (Intercept) coeff: -1.6094379 Pr(>|t|): 0.01093576 *
## chamberA-2 coeff: 0.8109302 Pr(>|t|): 0.2789959
## chamberA-3 coeff: 1.0498221 Pr(>|t|): 0.1740305
## chamberA-4 coeff: 1.0388931 Pr(>|t|): 0.1498307
## chamberB-1 coeff: 1.2417131 Pr(>|t|): 0.1053886
## chamberB-2 coeff: 0.8556661 Pr(>|t|): 0.2627736
## chamberB-3 coeff: 0.4307829 Pr(>|t|): 0.5660442
## chamberB-4 coeff: 1.3411739 Pr(>|t|): 0.06690121 .## glm flew_b ~ days_from_start binomial data_fem
## AIC: 274.8378
## (Intercept) coeff: -0.3276578 Pr(>|t|): 0.228011
## days_from_start coeff: -0.0345995 Pr(>|t|): 0.1007042## glm flew_b ~ min_from_IncStart binomial data_fem
## AIC: 277.2395
## (Intercept) coeff: -0.6067624 Pr(>|t|): 0.01010641 *
## min_from_IncStart coeff: -0.0005694 Pr(>|t|): 0.5687423Biological Effects
## glm flew_b ~ mass_c binomial data_fem
## AIC: 259.8795
## (Intercept) coeff: -0.7797155 Pr(>|t|): 4.942528e-07 *
## mass_c coeff: -25.14884 Pr(>|t|): 0.0005365364 *## glm flew_b ~ total_eggs_c binomial data_fem
## AIC: 203.7404
## (Intercept) coeff: -1.1924462 Pr(>|t|): 4.467475e-09 *
## total_eggs_c coeff: -0.0113041 Pr(>|t|): 3.560418e-05 *## glm flew_b ~ eggs_b binomial data_fem
## AIC: 226.3734
## (Intercept) coeff: 0.5596158 Pr(>|t|): 0.0181655 *
## eggs_b coeff: -2.2307473 Pr(>|t|): 1.790399e-11 *Morphology Effects
## glm flew_b ~ beak_c binomial data_fem
## AIC: 269.434
## (Intercept) coeff: -0.743354 Pr(>|t|): 6.111363e-07 *
## beak_c coeff: 0.5097352 Pr(>|t|): 0.005387442 *## glm flew_b ~ thorax_c binomial data_fem
## AIC: 275.7036
## (Intercept) coeff: -0.7206805 Pr(>|t|): 7.691419e-07 *
## thorax_c coeff: 0.7177388 Pr(>|t|): 0.1741144## glm flew_b ~ body_c binomial data_fem
## AIC: 271.6644
## (Intercept) coeff: -0.7373394 Pr(>|t|): 6.431673e-07 *
## body_c coeff: 0.3567498 Pr(>|t|): 0.01823304 *## glm flew_b ~ wing_c binomial data_fem
## AIC: 271.5651
## (Intercept) coeff: -0.7384629 Pr(>|t|): 6.342982e-07 *
## wing_c coeff: 0.4261769 Pr(>|t|): 0.01788508 *Without Mass
R = data_fem$flew_b
A = data_fem$host_c
B = data_fem$wing_c # before had sym_dist but now leave out because did not have significant effects or interactions within the total model dataset.
C = data_fem$mass_c
data<-data.frame(R, A, B, C)
model_script = paste0(source_path,"generic models-binomial glm 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4]
## AICs 271.5651 272.7081 273.4338 275.5669
## models 2 4 3 5
## probs 0.4643959 0.2622274 0.1824284 0.06279057
##
## m2 glm(formula = R ~ B, family = binomial, data = data)
## m4 glm(formula = R ~ A * B, family = binomial, data = data)
## m3 glm(formula = R ~ A + B, family = binomial, data = data)
## m0 glm(formula = R ~ 1, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m0, m2, test="Chisq") # Adding B does improve fit## Analysis of Deviance Table
##
## Model 1: R ~ 1
## Model 2: R ~ B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 215 273.57
## 2 214 267.56 1 6.0019 0.01429 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2, m3, test="Chisq") # Adding A does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ B
## Model 2: R ~ A + B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 214 267.56
## 2 213 267.43 1 0.13124 0.7171## glm flew_b ~ wing_c binomial data_fem
## AIC: 271.5651
## (Intercept) coeff: -0.7384629 Pr(>|t|): 6.342982e-07 *
## wing_c coeff: 0.4261769 Pr(>|t|): 0.01788508 *- positive effect of wing length where the longer the wing then female more likely to fly.
With Mass
data_fem <- data_fem %>%
filter(!is.na(mass))
data_fem <- center_data(data_fem)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthR = data_fem$flew_b
A = data_fem$host_c
B = data_fem$wing_c # used to be sym_dis
C = data_fem$mass_c
data<-data.frame(R, A, B, C)
model_script = paste0(source_path,"generic models-binomial glm 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 243.8633 244.0385 245.1296 245.4473 245.6227 245.8632 246.0379
## models 6 11 12 14 7 10 15
## probs 0.22037 0.2018842 0.1169977 0.09981397 0.09143357 0.08107111 0.07429147
##
## m6 glm(formula = R ~ B + C, family = binomial, data = data)
## m11 glm(formula = R ~ A * B + C, family = binomial, data = data)
## m12 glm(formula = R ~ A * C + B, family = binomial, data = data)
## m14 glm(formula = R ~ A * B + A * C, family = binomial, data = data)
## m7 glm(formula = R ~ A + B + C, family = binomial, data = data)
## m10 glm(formula = R ~ B * C, family = binomial, data = data)
## m15 glm(formula = R ~ A * B + B * C, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m3, m6, test="Chisq") # Adding B does improve fit## Analysis of Deviance Table
##
## Model 1: R ~ C
## Model 2: R ~ B + C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 211 255.88
## 2 210 237.86 1 18.016 2.19e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m3, m5, test="Chisq") # Adding A does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ C
## Model 2: R ~ A + C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 211 255.88
## 2 210 255.81 1 0.066494 0.7965anova(m6, m7, test="Chisq") # Adding A does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ B + C
## Model 2: R ~ A + B + C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 210 237.86
## 2 209 237.62 1 0.24061 0.6238anova(m6, m10, test="Chisq") # Adding B*C does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ B + C
## Model 2: R ~ B * C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 210 237.86
## 2 209 237.86 1 3.7622e-05 0.9951## glm flew_b ~ wing_c + mass_c binomial data_fem
## AIC: 243.8633
## (Intercept) coeff: -0.8643395 Pr(>|t|): 2.342868e-07 *
## wing_c coeff: 0.8582809 Pr(>|t|): 9.339491e-05 *
## mass_c coeff: -37.2787841 Pr(>|t|): 5.071396e-06 *- large negative effect of mass, where the heavier the bug is the less likely the bus will fly/disperse.
- positve effect of wing morph where the longer the wing then more likely to disperse.
Consider for females: A=host, B=mass, C=wing, X=ID
R = data_fem$flew_b
A = data_fem$host_c
B = data_fem$mass_c
C = data_fem$wing_c # replaced sym_dist
X = data_fem$ID # variance is zero if use trial_type | 3 models fail to converge if use ID
data<-data.frame(R, A, B, C, X)
model_script = paste0(source_path,"generic models-binomial glmer 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4]
## AICs 203.132 204.4579 206.663 206.7247
## models 14 17 11 6
## probs 0.4219636 0.2174506 0.07219793 0.07000463
##
## m14 R ~ A * B + A * C + (1 | X)
## m17 R ~ A * B + A * C + B * C + (1 | X)
## m11 R ~ A * B + C + (1 | X)
## m6 R ~ B + C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 4# top model if use ID:
multi_glmer_fem <-glmer(flew_b~host_c * mass_c + host_c * wing_c + (1 | ID), family=binomial, data=data_fem)
tidy_regression(multi_glmer_fem, is_color=output_col) ## glmer flew_b ~ host_c * mass_c + host_c * wing_c + (1 | ID) data_fem binomial
## AIC: 203.132 203.132
## (Intercept) coeff: -14.5424939 Pr(>|t|): 0.0001255127 *
## host_c coeff: -3.4870166 Pr(>|t|): 0.1065177
## mass_c coeff: -358.5446649 Pr(>|t|): 0.0006948939 *
## wing_c coeff: 8.6005138 Pr(>|t|): 0.002515753 *
## host_c:mass_c coeff: -194.6760791 Pr(>|t|): 0.01740641 *
## host_c:wing_c coeff: 6.8070032 Pr(>|t|): 0.01524722 *# top model if use trial_type:
multi_glmer_fem2 <-glmer(flew_b~mass_c + wing_c + (1|trial_type), family=binomial, data=data_fem)## boundary (singular) fit: see ?isSingular# if you replace trial_type with ID then you'll see the coefficients are way too huge and at different scales. Why?
tidy_regression(multi_glmer_fem2, is_color=output_col) ## glmer flew_b ~ mass_c + wing_c + (1 | trial_type) data_fem binomial
## AIC: 245.8633 245.8633
## (Intercept) coeff: -0.8643395 Pr(>|t|): 2.348399e-07 *
## mass_c coeff: -37.2787841 Pr(>|t|): 5.170195e-06 *
## wing_c coeff: 0.8582809 Pr(>|t|): 9.362589e-05 *Eggs
R = data_fem$flew_b
A = data_fem$host_c
B = data_fem$wing_c
C = data_fem$mass_c
D = data_fem$eggs_b
data<-data.frame(R, A, B, C, D)
model_script = paste0(source_path,"generic models-binomial glm 4-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3]
## AICs 217.2239 217.6946 217.9508
## models 14 34 35
## probs 0.07362638 0.05818669 0.05118926
##
## m14 glm(formula = R ~ B + C + D, family = binomial, data = data)
## m34 glm(formula = R ~ A * B + C + D, family = binomial, data = data)
## m35 glm(formula = R ~ A * C + B + D, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0egg_model2 <-glm(flew_b~wing_c + mass_c + eggs_b, family=binomial, data=data_fem)
tidy_regression(egg_model2, is_color=output_col)## glm flew_b ~ wing_c + mass_c + eggs_b binomial data_fem
## AIC: 217.2239
## (Intercept) coeff: 0.3324593 Pr(>|t|): 0.2344807
## wing_c coeff: 0.6277736 Pr(>|t|): 0.008255178 *
## mass_c coeff: -16.9840966 Pr(>|t|): 0.06369226 .
## eggs_b coeff: -1.9388382 Pr(>|t|): 2.680918e-07 *- Strong negative effect if laid eggs that day
- strong marginal negative effect of mass
- positive effect of wing length
egg_glmer <-glmer(flew_b~wing_c + mass_c + eggs_b + (1|ID), family=binomial, data=data_fem)
tidy_regression(egg_glmer, is_color=output_col) # model fits better## glmer flew_b ~ wing_c + mass_c + eggs_b + (1 | ID) data_fem binomial
## AIC: 175.24 175.24
## (Intercept) coeff: 5.779837 Pr(>|t|): 0.01174729 *
## wing_c coeff: 1.8027952 Pr(>|t|): 0.2580401
## mass_c coeff: -81.8894826 Pr(>|t|): 0.1449163
## eggs_b coeff: -15.887739 Pr(>|t|): 1.17845e-07 *still a strong negative effect if laid eggs that day but mass and wing length are no longer significant…
glmer() A=host, B=wing_c, C=mass, D=eggs_b, X=ID
R = data_fem$flew_b
A = data_fem$host_c
B = data_fem$wing_c
C = data_fem$mass_c
D = data_fem$eggs_b
X = data_fem$ID
data<-data.frame(R, A, B, C, D, X)
model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 4-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4]
## AICs 195.2918 197.7208 198.2143 198.4093
## models 53 21 33 86
## probs 0.4395504 0.1304878 0.1019551 0.09248245
##
## m53 R ~ B * C + C * D + (1 | X)
## m21 R ~ C * D + (1 | X)
## m33 R ~ C * D + B + (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: 0anova(m21, m33, test="Chisq") # Adding B improves fits## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m21: R ~ C * D + (1 | X)
## m33: R ~ C * D + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m21 6 194.90 215.07 -91.452 182.90
## m33 7 190.78 214.31 -88.392 176.78 6.1191 1 0.01337 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m33, m53, test="Chisq") # Adding B*C does not improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m33: R ~ C * D + B + (1 | X)
## m53: R ~ B * C + C * D + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m33 7 190.78 214.31 -88.392 176.78
## m53 8 190.52 217.41 -87.262 174.52 2.2609 1 0.1327egg_glmer2<-glmer(flew_b~mass_c *eggs_b + wing_c + (1|ID), family=binomial, data=data_fem)
tidy_regression(egg_glmer2, is_color=output_col)## glmer flew_b ~ mass_c * eggs_b + wing_c + (1 | ID) data_fem binomial
## AIC: 171.3691 171.3691
## (Intercept) coeff: -1.3889397 Pr(>|t|): 0.6520714
## mass_c coeff: -423.980237 Pr(>|t|): 0.01574343 *
## eggs_b coeff: -9.2371298 Pr(>|t|): 0.0002705239 *
## wing_c coeff: 2.5592644 Pr(>|t|): 0.3315875
## mass_c:eggs_b coeff: 361.4810625 Pr(>|t|): 0.02710765 *coefficients are extremely out of scale. * negative effect of mass * negative effect of eggs positive effect of wing positivie effect mass*eggs? doesn’t really make sense. –> collinearity
Morphology
R = data_fem$flew_b
A = data_fem$thorax_c
B = data_fem$body_c
C = data_fem$wing_c
data<-data.frame(R, A, B, C)
model_script = paste0(source_path,"generic models-binomial glm 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 265.5347 265.6369 266.4738 266.9023 267.3331 267.3598 267.8101
## models 16 15 4 5 17 3 2
## probs 0.163925 0.1557537 0.1024978 0.08272912 0.06669857 0.0658151 0.05254679
##
## m16 glm(formula = R ~ A * C + B * C, family = binomial, data = data)
## m15 glm(formula = R ~ A * B + B * C, family = binomial, data = data)
## m4 glm(formula = R ~ A + B, family = binomial, data = data)
## m5 glm(formula = R ~ A + C, family = binomial, data = data)
## m17 glm(formula = R ~ A * B + A * C + B * C, family = binomial, data = data)
## m3 glm(formula = R ~ C, family = binomial, data = data)
## m2 glm(formula = R ~ B, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m13, m16, test="Chisq") # Adding A*C improves fit## Analysis of Deviance Table
##
## Model 1: R ~ B * C + A
## Model 2: R ~ A * C + B * C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 208 257.96
## 2 207 253.53 1 4.428 0.03535 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1morph_fem2 <- glm(flew_b ~ thorax_c * wing_c + body_c * wing_c, data=data_fem)
tidy_regression(morph_fem2, is_color=output_col)## glm flew_b ~ thorax_c * wing_c + body_c * wing_c data_fem
## AIC: 286.9107
## (Intercept) coeff: 0.3431652 Pr(>|t|): 8.437584e-15 *
## thorax_c coeff: -0.3812667 Pr(>|t|): 0.1847385
## wing_c coeff: 0.0074331 Pr(>|t|): 0.9708583
## body_c coeff: 0.1484849 Pr(>|t|): 0.4500281
## thorax_c:wing_c coeff: 0.1167622 Pr(>|t|): 0.6453306
## wing_c:body_c coeff: -0.0440551 Pr(>|t|): 0.4384809- no significant effects
morph_fem_glmer <- glmer(flew_b ~ thorax_c * wing_c + body_c * wing_c + (1|ID), data=data_fem, family=binomial)
tidy_regression(morph_fem_glmer, is_color=output_col) #even better fitted model; effects no longer significant and coefficients changed. ## glmer flew_b ~ thorax_c * wing_c + body_c * wing_c + (1 | ID) data_fem binomial
## AIC: 237.7825 237.7825
## (Intercept) coeff: -6.416729 Pr(>|t|): 3.219458e-05 *
## thorax_c coeff: -7.1858445 Pr(>|t|): 0.3667043
## wing_c coeff: 0.1655431 Pr(>|t|): 0.9636254
## body_c coeff: 2.5171036 Pr(>|t|): 0.5259
## thorax_c:wing_c coeff: 8.0921535 Pr(>|t|): 0.4434568
## wing_c:body_c coeff: -2.698719 Pr(>|t|): 0.4019887Trying wing2body
R = data_fem$flew_b
A = data_fem$thorax_c
B = data_fem$wing2body_c
data<-data.frame(R, A, B)
model_script = paste0(source_path,"generic models-binomial glm 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5]
## AICs 268.7935 269.1631 270.5935 271.7234 271.7476
## models 2 4 3 1 5
## probs 0.37075 0.3081929 0.1507324 0.08567727 0.08464737
##
## m2 glm(formula = R ~ B, family = binomial, data = data)
## m4 glm(formula = R ~ A * B, family = binomial, data = data)
## m3 glm(formula = R ~ A + B, family = binomial, data = data)
## m1 glm(formula = R ~ A, family = binomial, data = data)
## m0 glm(formula = R ~ 1, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m0, m2, test="Chisq") # Adding B improves fit## Analysis of Deviance Table
##
## Model 1: R ~ 1
## Model 2: R ~ B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 212 269.75
## 2 211 264.79 1 4.9541 0.02603 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2, m3, test="Chisq") # Adding A does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ B
## Model 2: R ~ A + B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 211 264.79
## 2 210 264.59 1 0.19996 0.6548morph_fem3 <- glm(flew_b ~ wing2body_c, data=data_fem, family=binomial) # variance of population and days_from_start_c = 0
tidy_regression(morph_fem3, is_color=output_col)## glm flew_b ~ wing2body_c binomial data_fem
## AIC: 268.7935
## (Intercept) coeff: -0.7354108 Pr(>|t|): 7.728476e-07 *
## wing2body_c coeff: 19.29705 Pr(>|t|): 0.03322314 *- strong positive effect of wing2body where the larger the ration (the more wing than body), the more likely the female bug disperses
Worth running a glmer() script.
Consider for females morph: A=thorax, B=body, C=wing X=trial_type
R = data_fem$flew_b
A = data_fem$thorax_c
B = data_fem$body_c
C = data_fem$wing_c
X = data_fem$ID # you get a different set of top models if you use trial_type or ID, you get the null model if ID and the usual model if trial_type
data<-data.frame(R, A, B, C, X)
model_script = paste0(source_path,"generic models-binomial glmer 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 230.1878 231.2014 231.218 231.8429 232.7172 232.8936 233.1978
## models 19 3 2 1 4 5 6
## probs 0.2501252 0.1506785 0.149432 0.1093346 0.0706141 0.06465224 0.05553038
##
## m0 R ~ 1 + (1 | X)
## m3 R ~ C + (1 | X)
## m2 R ~ B + (1 | X)
## m1 R ~ A + (1 | X)
## m4 R ~ A + B + (1 | X)
## m5 R ~ A + C + (1 | X)
## m6 R ~ B + C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 1anova(m0, m3, test="Chisq") # Adding C does not improve fit## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m3: R ~ C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 2 230.19 236.91 -113.09 226.19
## m3 3 231.20 241.28 -112.60 225.20 0.9864 1 0.3206anova(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 2 230.19 236.91 -113.09 226.19
## m2 3 231.22 241.30 -112.61 225.22 0.9698 1 0.3247anova(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 2 230.19 236.91 -113.09 226.19
## m1 3 231.84 241.93 -112.92 225.84 0.3449 1 0.557# when use ID and trial type as RFs get the null model as the top model
morph_fem_glmer3 <- glmer(flew_b ~ 1 + (1|ID), data=data_fem, family=binomial)
tidy_regression(morph_fem_glmer3 , is_color=output_col) ## glmer flew_b ~ 1 + (1 | ID) data_fem binomial
## AIC: 230.1878 230.1878
## 1 coeff: -7.3543181 Pr(>|t|): 1.49237e-11 *# when use only trial type as RF:
morph_fem_glmer4 <- glmer(flew_b ~ thorax_c * wing_c + body_c * wing_c + (1|trial_type), data=data_fem, family=binomial)
tidy_regression(morph_fem_glmer4 , is_color=output_col) ## glmer flew_b ~ thorax_c * wing_c + body_c * wing_c + (1 | trial_type) data_fem binomial
## AIC: 267.2253 267.2253
## (Intercept) coeff: -0.5902431 Pr(>|t|): 0.01223559 *
## thorax_c coeff: -3.8030181 Pr(>|t|): 0.02886067 *
## wing_c coeff: 0.0956551 Pr(>|t|): 0.9205474
## body_c coeff: 1.3297507 Pr(>|t|): 0.178274
## thorax_c:wing_c coeff: 4.371407 Pr(>|t|): 0.04916682 *
## wing_c:body_c coeff: -1.4756477 Pr(>|t|): 0.02992012 *Trying wing2body
R = data_fem$flew_b
A = data_fem$thorax_c
B = data_fem$wing2body_c
X = data_fem$trial_type # you get a different set of top models if you use trial_type or ID
data<-data.frame(R, A, B, X)
model_script = paste0(source_path,"generic models-binomial glmer 1-RF + 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5]
## AICs 270.687 271.0019 272.4749 273.5265 273.5722
## models 2 4 3 1 5
## probs 0.3647736 0.3116245 0.1492076 0.08819381 0.08620041
##
## m2 R ~ B + (1 | X)
## m4 R ~ A * B + (1 | X)
## m3 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: 0anova(m0, m2, test="Chisq") # Adding A marginally imprvoes fit## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 2 273.57 280.30 -134.79 269.57
## m2 3 270.69 280.77 -132.34 264.69 4.8852 1 0.02709 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m0, m1, test="Chisq")## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 2 273.57 280.30 -134.79 269.57
## m1 3 273.53 283.61 -133.76 267.53 2.0457 1 0.1526# best fit model if only use trial type as RF:
morph_fem_glmer5 <- glmer(flew_b ~ wing2body_c + (1|trial_type), data=data_fem, family=binomial)
tidy_regression(morph_fem_glmer5 , is_color=output_col) ## glmer flew_b ~ wing2body_c + (1 | trial_type) data_fem binomial
## AIC: 270.687 270.687
## (Intercept) coeff: -0.7440755 Pr(>|t|): 2.431814e-05 *
## wing2body_c coeff: 19.20888 Pr(>|t|): 0.03433405 *# best fit model is the null if use both trial type and ID as RFs:
morph_fem_glmer6 <- glmer(flew_b ~ 1 + (1|ID), data=data_fem, family=binomial)
tidy_regression(morph_fem_glmer6, is_color=output_col) ## glmer flew_b ~ 1 + (1 | ID) data_fem binomial
## AIC: 230.1878 230.1878
## 1 coeff: -7.3543181 Pr(>|t|): 1.49237e-11 *Summary
tidy_regression(nomass_fem2, is_color=output_col) ## glm flew_b ~ wing_c binomial data_fem
## AIC: 271.5651
## (Intercept) coeff: -0.7384629 Pr(>|t|): 6.342982e-07 *
## wing_c coeff: 0.4261769 Pr(>|t|): 0.01788508 *tidy_regression(mass_fem2, is_color=output_col) ## glm flew_b ~ wing_c + mass_c binomial data_fem
## AIC: 243.8633
## (Intercept) coeff: -0.8643395 Pr(>|t|): 2.342868e-07 *
## wing_c coeff: 0.8582809 Pr(>|t|): 9.339491e-05 *
## mass_c coeff: -37.2787841 Pr(>|t|): 5.071396e-06 *tidy_regression(multi_glmer_fem, is_color=output_col) ## glmer flew_b ~ host_c * mass_c + host_c * wing_c + (1 | ID) data_fem binomial
## AIC: 203.132 203.132
## (Intercept) coeff: -14.5424939 Pr(>|t|): 0.0001255127 *
## host_c coeff: -3.4870166 Pr(>|t|): 0.1065177
## mass_c coeff: -358.5446649 Pr(>|t|): 0.0006948939 *
## wing_c coeff: 8.6005138 Pr(>|t|): 0.002515753 *
## host_c:mass_c coeff: -194.6760791 Pr(>|t|): 0.01740641 *
## host_c:wing_c coeff: 6.8070032 Pr(>|t|): 0.01524722 *tidy_regression(multi_glmer_fem2, is_color=output_col) ## glmer flew_b ~ mass_c + wing_c + (1 | trial_type) data_fem binomial
## AIC: 245.8633 245.8633
## (Intercept) coeff: -0.8643395 Pr(>|t|): 2.348399e-07 *
## mass_c coeff: -37.2787841 Pr(>|t|): 5.170195e-06 *
## wing_c coeff: 0.8582809 Pr(>|t|): 9.362589e-05 *tidy_regression(egg_model2, is_color=output_col)## glm flew_b ~ wing_c + mass_c + eggs_b binomial data_fem
## AIC: 217.2239
## (Intercept) coeff: 0.3324593 Pr(>|t|): 0.2344807
## wing_c coeff: 0.6277736 Pr(>|t|): 0.008255178 *
## mass_c coeff: -16.9840966 Pr(>|t|): 0.06369226 .
## eggs_b coeff: -1.9388382 Pr(>|t|): 2.680918e-07 *tidy_regression(egg_glmer, is_color=output_col)## glmer flew_b ~ wing_c + mass_c + eggs_b + (1 | ID) data_fem binomial
## AIC: 175.24 175.24
## (Intercept) coeff: 5.779837 Pr(>|t|): 0.01174729 *
## wing_c coeff: 1.8027952 Pr(>|t|): 0.2580401
## mass_c coeff: -81.8894826 Pr(>|t|): 0.1449163
## eggs_b coeff: -15.887739 Pr(>|t|): 1.17845e-07 *tidy_regression(egg_glmer2, is_color=output_col)## glmer flew_b ~ mass_c * eggs_b + wing_c + (1 | ID) data_fem binomial
## AIC: 171.3691 171.3691
## (Intercept) coeff: -1.3889397 Pr(>|t|): 0.6520714
## mass_c coeff: -423.980237 Pr(>|t|): 0.01574343 *
## eggs_b coeff: -9.2371298 Pr(>|t|): 0.0002705239 *
## wing_c coeff: 2.5592644 Pr(>|t|): 0.3315875
## mass_c:eggs_b coeff: 361.4810625 Pr(>|t|): 0.02710765 *All the morph models had an issue between using ID vs. trial type as the RF. When I used trial_type I got reasonable results. When I used ID I always got the null model. Why is this?
tidy_regression(morph_fem2, is_color=output_col)## glm flew_b ~ thorax_c * wing_c + body_c * wing_c data_fem
## AIC: 286.9107
## (Intercept) coeff: 0.3431652 Pr(>|t|): 8.437584e-15 *
## thorax_c coeff: -0.3812667 Pr(>|t|): 0.1847385
## wing_c coeff: 0.0074331 Pr(>|t|): 0.9708583
## body_c coeff: 0.1484849 Pr(>|t|): 0.4500281
## thorax_c:wing_c coeff: 0.1167622 Pr(>|t|): 0.6453306
## wing_c:body_c coeff: -0.0440551 Pr(>|t|): 0.4384809tidy_regression(morph_fem_glmer, is_color=output_col) ## glmer flew_b ~ thorax_c * wing_c + body_c * wing_c + (1 | ID) data_fem binomial
## AIC: 237.7825 237.7825
## (Intercept) coeff: -6.416729 Pr(>|t|): 3.219458e-05 *
## thorax_c coeff: -7.1858445 Pr(>|t|): 0.3667043
## wing_c coeff: 0.1655431 Pr(>|t|): 0.9636254
## body_c coeff: 2.5171036 Pr(>|t|): 0.5259
## thorax_c:wing_c coeff: 8.0921535 Pr(>|t|): 0.4434568
## wing_c:body_c coeff: -2.698719 Pr(>|t|): 0.4019887tidy_regression(morph_fem3, is_color=output_col)## glm flew_b ~ wing2body_c binomial data_fem
## AIC: 268.7935
## (Intercept) coeff: -0.7354108 Pr(>|t|): 7.728476e-07 *
## wing2body_c coeff: 19.29705 Pr(>|t|): 0.03322314 *tidy_regression(morph_fem_glmer3 , is_color=output_col) ## glmer flew_b ~ 1 + (1 | ID) data_fem binomial
## AIC: 230.1878 230.1878
## 1 coeff: -7.3543181 Pr(>|t|): 1.49237e-11 *tidy_regression(morph_fem_glmer4 , is_color=output_col) ## glmer flew_b ~ thorax_c * wing_c + body_c * wing_c + (1 | trial_type) data_fem binomial
## AIC: 267.2253 267.2253
## (Intercept) coeff: -0.5902431 Pr(>|t|): 0.01223559 *
## thorax_c coeff: -3.8030181 Pr(>|t|): 0.02886067 *
## wing_c coeff: 0.0956551 Pr(>|t|): 0.9205474
## body_c coeff: 1.3297507 Pr(>|t|): 0.178274
## thorax_c:wing_c coeff: 4.371407 Pr(>|t|): 0.04916682 *
## wing_c:body_c coeff: -1.4756477 Pr(>|t|): 0.02992012 *tidy_regression(morph_fem_glmer5 , is_color=output_col) ## glmer flew_b ~ wing2body_c + (1 | trial_type) data_fem binomial
## AIC: 270.687 270.687
## (Intercept) coeff: -0.7440755 Pr(>|t|): 2.431814e-05 *
## wing2body_c coeff: 19.20888 Pr(>|t|): 0.03433405 *tidy_regression(morph_fem_glmer6, is_color=output_col)## glmer flew_b ~ 1 + (1 | ID) data_fem binomial
## AIC: 230.1878 230.1878
## 1 coeff: -7.3543181 Pr(>|t|): 1.49237e-11 *Males
Cleaning the Data
data <- read_flight_data("data/all_flight_data-Winter2020.csv")## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object length
## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_all <- data[[1]]
data_tested <- data[[2]]
d = create_delta_data(data_tested)
data_male <- d[d$sex=="M",]
data_male <- center_data(data_male, is_not_unique_data = FALSE) ## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_male## # A tibble: 185 x 86
## # Groups: ID, sex, population, site, host_plant, latitude, longitude,
## # total_eggs, beak, thorax, wing, body, w_morph, morph_notes, tested, host_c,
## # sex_c, w_morph_c, lat_c, sym_dist, sym_dist_s, total_eggs_c, beak_c,
## # thorax_c, thorax_s, body_c, wing_c, wing2body, wing2body_c [185]
## ID sex population site host_plant latitude longitude total_eggs beak
## <fct> <fct> <fct> <fct> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1 M Plantation… Arego… C. corind… 25.0 -80.6 NA 5.5
## 2 2 M Plantation… Arego… C. corind… 25.0 -80.6 NA 5.64
## 3 3 M Plantation… Arego… C. corind… 25.0 -80.6 NA 5.75
## 4 4 M Plantation… Arego… C. corind… 25.0 -80.6 NA 6.21
## 5 6 M Plantation… Found… C. corind… 25.0 -80.6 NA 5.76
## 6 9 M Plantation… Found… C. corind… 25.0 -80.6 NA 6.1
## 7 10 M Plantation… Found… C. corind… 25.0 -80.6 NA 6.09
## 8 11 M Key Largo KLMRL C. corind… 25.1 -80.4 NA 5.88
## 9 14 M Key Largo KLMRL C. corind… 25.1 -80.4 NA 5.25
## 10 15 M Key Largo JP C. corind… 25.1 -80.4 NA 6.53
## # … with 175 more rows, and 77 more variables: thorax <dbl>, wing <dbl>,
## # body <dbl>, w_morph <fct>, morph_notes <fct>, tested <fct>, host_c <dbl>,
## # sex_c <dbl>, w_morph_c <dbl>, lat_c <dbl>, sym_dist <dbl>,
## # sym_dist_s <dbl>, total_eggs_c <dbl>, beak_c <dbl>, thorax_c <dbl>,
## # thorax_s <dbl>, body_c <dbl>, wing_c <dbl>, wing2body <dbl>,
## # wing2body_c <dbl>, wing2body_s <dbl>, trial_type <list>, chamber <list>,
## # average_speed <list>, total_flight_time <list>, distance <list>,
## # shortest_flying_bout <list>, longest_flying_bout <list>,
## # recording_duration <list>, max_speed <list>, test_date <list>,
## # time_start <list>, time_end <list>, flew <list>, flight_type <list>,
## # mass <list>, EWM <list>, flew_b <list>, host_temp <list>, eggs_b <list>,
## # minute_duration <list>, minute_duration_c <list>, min_from_IncStart <dbl>,
## # min_from_IncStart_c <dbl>, min_from_IncStart_s <dbl>,
## # days_from_start <list>, days_from_start_c <list>, mass_c <list>,
## # mass_s <list>, avg_mass <dbl>, avg_mass_c <dbl>, avg_mass_s <dbl>,
## # trial_type_og <list>, dt_start_c <list>, dt_end_c <list>,
## # datetime_start <list>, datetime_end <list>, hr_start <list>, hr_end <list>,
## # num_flew <dbl>, num_notflew <dbl>, avg_days <dbl>, num_egg <dbl>,
## # egg_diff <dbl>, mass_diff <dbl>, flew_diff <dbl>, dist_diff <dbl>,
## # speed_diff <dbl>, mass_per <dbl>, avg_rec_dur <dbl>, avg_time_start <dbl>,
## # flight_case <fct>, egg_case <dbl>, flew_prob <dbl>, avg_days_c <dbl>,
## # avg_mass_logsqrt <dbl>, wing2body_logsqrt_i <dbl>Experimental, Biological, and Morphological Effects
## glm cbind(num_flew, num_notflew) ~ avg_mass binomial data_male
## AIC: 396.8874
## (Intercept) coeff: 0.6024 Pr(>|t|): 0.3759086
## avg_mass coeff: 5.6296207 Pr(>|t|): 0.7440755## glm cbind(num_flew, num_notflew) ~ beak_c binomial data_male
## AIC: 393.8903
## (Intercept) coeff: 0.8296418 Pr(>|t|): 2.964899e-13 *
## beak_c coeff: 0.5133225 Pr(>|t|): 0.08098797 .## glm cbind(num_flew, num_notflew) ~ thorax_c binomial data_male
## AIC: 396.6049
## (Intercept) coeff: 0.8226484 Pr(>|t|): 3.297608e-13 *
## thorax_c coeff: 0.3574224 Pr(>|t|): 0.5326183## glm cbind(num_flew, num_notflew) ~ body_c binomial data_male
## AIC: 392.488
## (Intercept) coeff: 0.8322146 Pr(>|t|): 2.85788e-13 *
## body_c coeff: 0.3502639 Pr(>|t|): 0.03396248 *## glm cbind(num_flew, num_notflew) ~ wing_c binomial data_male
## AIC: 389.9123
## (Intercept) coeff: 0.837782 Pr(>|t|): 2.635901e-13 *
## wing_c coeff: 0.5108897 Pr(>|t|): 0.008114818 *Without Mass
R1 = data_male$num_flew
R2 = data_male$num_notflew
A = data_male$host_c
B = data_male$wing_c
C = data_male$avg_mass
X = data_male$population # can't do trial type anymore
data<-data.frame(R1, R2, A, B, C, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3]
## AICs 385.9057 387.7504 388.6232
## models 2 3 4
## probs 0.5757326 0.2288993 0.1479542
##
## m2 cbind(R1, R2) ~ B + (1 | X)
## m3 cbind(R1, R2) ~ A + B + (1 | X)
## m4 cbind(R1, R2) ~ A * B + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m2, m3, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m2: cbind(R1, R2) ~ B + (1 | X)
## m3: cbind(R1, R2) ~ A + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 3 385.91 395.57 -189.95 379.91
## m3 4 387.75 400.63 -189.88 379.75 0.1553 1 0.6935nomass_male <- glmer(cbind(num_flew, num_notflew) ~ wing_c + (1|population), data=data_male, family=binomial)
tidy_regression(nomass_male, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing_c + (1 | population) data_male binomial
## AIC: 385.9057 385.9057
## (Intercept) coeff: 0.7365366 Pr(>|t|): 0.0004586381 *
## wing_c coeff: 0.5742953 Pr(>|t|): 0.005717746 *- positive effect of wing
With Mass
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 382.6366 383.1203 383.5774 383.6552 384.0872 385.1117 385.2434
## models 6 14 10 17 7 13 12
## probs 0.1991648 0.1563822 0.1244322 0.1196797 0.09643019 0.05777622 0.05409489
## [,8]
## AICs 385.2604
## models 11
## probs 0.05363823
##
## m6 cbind(R1, R2) ~ B + C + (1 | X)
## m14 cbind(R1, R2) ~ A * B + A * C + (1 | X)
## m10 cbind(R1, R2) ~ B * C + (1 | X)
## m17 cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)
## m7 cbind(R1, R2) ~ A + B + C + (1 | X)
## m13 cbind(R1, R2) ~ B * C + A + (1 | X)
## m12 cbind(R1, R2) ~ A * C + B + (1 | X)
## m11 cbind(R1, R2) ~ A * B + C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m6, m10, test="Chisq") # Adding B*C does not improve fit## Data: data
## Models:
## m6: cbind(R1, R2) ~ B + C + (1 | X)
## m10: cbind(R1, R2) ~ B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m6 4 382.64 395.52 -187.32 374.64
## m10 5 383.58 399.68 -186.79 373.58 1.0593 1 0.3034anova(m6, m7, test="Chisq") # Adding A does not improve fot## Data: data
## Models:
## m6: cbind(R1, R2) ~ B + C + (1 | X)
## m7: cbind(R1, R2) ~ A + B + C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m6 4 382.64 395.52 -187.32 374.64
## m7 5 384.09 400.19 -187.04 374.09 0.5494 1 0.4586mass_male <- glmer(cbind(num_flew, num_notflew) ~ wing_c + avg_mass + (1|population), data=data_male, family=binomial)
tidy_regression(mass_male, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing_c + avg_mass + (1 | population) data_male binomial
## AIC: 382.6366 382.6366
## (Intercept) coeff: 2.876638 Pr(>|t|): 0.002539525 *
## wing_c coeff: 1.0077334 Pr(>|t|): 0.0004203606 *
## avg_mass coeff: -55.6507412 Pr(>|t|): 0.02060649 *- strong positive effect of wing
- strong negative effect of average mass
# check for collinearity between mass and wing length
data<-data.frame(data_male$host_c, data_male$wing_c, data_male$avg_mass)
colnames(data) <- c("Host Plant", "Wing Length", "Average Mass")
pairs(data)
Morphology
d <- data_male %>%
filter(!is.na(body))
d$thorax_c <- d$thorax - mean(d$thorax)
d$wing_c <- d$wing - mean(d$wing)
d$body_c <- d$body - mean(d$body)
R1 = d$num_flew
R2 = d$num_notflew
A = d$thorax_c
B = d$body_c
C = d$wing_c
X = d$population
data<-data.frame(R1, R2, A, B, C, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 379.0646 379.3479 379.9851 380.1507 380.2589 380.6564 380.7822
## models 16 13 10 5 15 17 9
## probs 0.1835943 0.1593475 0.1158705 0.106662 0.1010466 0.08283175 0.07778032
##
## m16 cbind(R1, R2) ~ A * C + B * C + (1 | X)
## m13 cbind(R1, R2) ~ B * C + A + (1 | X)
## m10 cbind(R1, R2) ~ B * C + (1 | X)
## m5 cbind(R1, R2) ~ A + C + (1 | X)
## m15 cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m17 cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)
## m9 cbind(R1, R2) ~ A * C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m13, m16, test="Chisq") # Adding B*C does not improve fit## Data: data
## Models:
## m13: cbind(R1, R2) ~ B * C + A + (1 | X)
## m16: cbind(R1, R2) ~ A * C + B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m13 6 379.35 398.67 -183.67 367.35
## m16 7 379.06 401.61 -182.53 365.06 2.2833 1 0.1308anova(m10, m13, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m10: cbind(R1, R2) ~ B * C + (1 | X)
## m13: cbind(R1, R2) ~ B * C + A + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m10 5 379.99 396.09 -184.99 369.99
## m13 6 379.35 398.67 -183.67 367.35 2.6372 1 0.1044anova(m6, m10, test="Chisq") # Adding B*C improves fit## Data: data
## Models:
## m6: cbind(R1, R2) ~ B + C + (1 | X)
## m10: cbind(R1, R2) ~ B * C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m6 4 383.77 396.65 -187.88 375.77
## m10 5 379.99 396.09 -184.99 369.99 5.7832 1 0.01618 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1morph_male <- glmer(cbind(num_flew, num_notflew) ~ body_c * wing_c + (1 | population), data=d, family=binomial)
tidy_regression(morph_male, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ body_c * wing_c + (1 | population) d binomial
## AIC: 379.9851 379.9851
## (Intercept) coeff: 0.9549736 Pr(>|t|): 7.787216e-05 *
## body_c coeff: -1.2383232 Pr(>|t|): 0.05485939 .
## wing_c coeff: 1.7059399 Pr(>|t|): 0.02543455 *
## body_c:wing_c coeff: -0.5575577 Pr(>|t|): 0.02037575 *- marginal negative effect of body
- positive effect of wing
- negative wing*body where the longer the wing and body the less times a male bug flies.
R1 = d$num_flew
R2 = d$num_notflew
A = d$thorax_c
B = d$wing2body_c
X = d$population
data<-data.frame(R1, R2, A, B, X)
model_script = paste0(source_path,"generic models-binomial glmer 2R ~ 1 RF + 2-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3]
## AICs 381.251 383.1217 384.2777
## models 2 3 4
## probs 0.6170651 0.2421659 0.1358595
##
## m2 cbind(R1, R2) ~ B + (1 | X)
## m3 cbind(R1, R2) ~ A + B + (1 | X)
## m4 cbind(R1, R2) ~ A * B + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m2, m3, test="Chisq") # Adding B does not improve fit## Data: data
## Models:
## m2: cbind(R1, R2) ~ B + (1 | X)
## m3: cbind(R1, R2) ~ A + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 3 381.25 390.91 -187.63 375.25
## m3 4 383.12 396.00 -187.56 375.12 0.1293 1 0.7192morph_male2 <- glmer(cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population), data=d, family=binomial)
tidy_regression(morph_male2, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population) d binomial
## AIC: 381.251 381.251
## (Intercept) coeff: 0.7315339 Pr(>|t|): 0.002148333 *
## wing2body_c coeff: 25.3265477 Pr(>|t|): 0.000655473 *- strong positive effect of wing2body where the more wing2body the more times a male bug flies.
# check for collinearity between mass and morhpology
data<-data.frame(data_male$thorax_c, data_male$wing_c, data_male$body_c, data_male$avg_mass)
colnames(data) <- c("Thorax Length", "Wing Length", "Body Length", "Average Mass")
pairs(data)
data_male <- data_tested[data_tested$sex=="M",]
data_male <- center_data(data_male) ## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata<-data.frame(data_male$thorax_c, data_male$wing_c, data_male$body_c, data_male$mass_c, data_male$days_from_start_c)
colnames(data) <- c("Thorax Length", "Wing Length", "Body Length", "Average Mass", "Days from Start")
pairs(data)
m <- lm(mass ~ days_from_start, data=data_male)
summary(m)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0364684619 6.846321e-04 53.267242 1.153823e-182
## days_from_start 0.0001997694 5.118653e-05 3.902772 1.117406e-04plot(mass ~ days_from_start, data=data_male)
abline(m, col="red")
#plot(body ~ days_from_start, data=data_male)
text(mass ~ days_from_start, labels=data_male$ID, data=data_male, cex=0.5, font=2, pos=4) # maybe the smaller males died faster? that's why it looks like.
Summary
tidy_regression(nomass_male, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing_c + (1 | population) data_male binomial
## AIC: 385.9057 385.9057
## (Intercept) coeff: 0.7365366 Pr(>|t|): 0.0004586381 *
## wing_c coeff: 0.5742953 Pr(>|t|): 0.005717746 *tidy_regression(mass_male, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing_c + avg_mass + (1 | population) data_male binomial
## AIC: 382.6366 382.6366
## (Intercept) coeff: 2.876638 Pr(>|t|): 0.002539525 *
## wing_c coeff: 1.0077334 Pr(>|t|): 0.0004203606 *
## avg_mass coeff: -55.6507412 Pr(>|t|): 0.02060649 *tidy_regression(morph_male, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ body_c * wing_c + (1 | population) d binomial
## AIC: 379.9851 379.9851
## (Intercept) coeff: 0.9549736 Pr(>|t|): 7.787216e-05 *
## body_c coeff: -1.2383232 Pr(>|t|): 0.05485939 .
## wing_c coeff: 1.7059399 Pr(>|t|): 0.02543455 *
## body_c:wing_c coeff: -0.5575577 Pr(>|t|): 0.02037575 *tidy_regression(morph_male2, is_color=output_col)## glmer cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population) d binomial
## AIC: 381.251 381.251
## (Intercept) coeff: 0.7315339 Pr(>|t|): 0.002148333 *
## wing2body_c coeff: 25.3265477 Pr(>|t|): 0.000655473 *Males - Old Method
Without Mass Modeling + Days From Start
With Mass Modeling + Days From Start
glmer() A=thorax, B=body, C=wing, D=days_from_start, X=ID, Y=trial_type
glmer() A=wing2body, B=thorax, C=days_from_start, X=ID, Y=trial_type
Cleaning Data
data <- read_flight_data("data/all_flight_data-Winter2020.csv")## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object length
## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_all <- data[[1]]
data_tested <- data[[2]]
data_male <- data_tested[data_tested$sex=="M",]
data_male <- center_data(data_male)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthExperimental Set-Up Effects
## glm flew_b ~ chamber binomial data_male
## AIC: 507.3498
## (Intercept) coeff: 0.6190392 Pr(>|t|): 0.06184483 .
## chamberA-2 coeff: 0.238411 Pr(>|t|): 0.604278
## chamberA-3 coeff: 0.36179 Pr(>|t|): 0.4203252
## chamberA-4 coeff: -0.1010961 Pr(>|t|): 0.8045408
## chamberB-1 coeff: 0.6131045 Pr(>|t|): 0.2585014
## chamberB-2 coeff: 0.4795731 Pr(>|t|): 0.2820852
## chamberB-3 coeff: 0.0095695 Pr(>|t|): 0.9831673
## chamberB-4 coeff: -0.1564157 Pr(>|t|): 0.7302194## glm flew_b ~ days_from_start_c binomial data_male
## AIC: 493.3651
## (Intercept) coeff: 0.7861825 Pr(>|t|): 6.92517e-13 *
## days_from_start_c coeff: -0.0431183 Pr(>|t|): 0.006662926 *## glm flew_b ~ min_from_IncStart binomial data_male
## AIC: 497.8801
## (Intercept) coeff: 0.53576 Pr(>|t|): 0.001679447 *
## min_from_IncStart coeff: 0.0013963 Pr(>|t|): 0.08561738 .Biological Effects
## glm flew_b ~ mass_c binomial data_male
## AIC: 499.9643
## (Intercept) coeff: 0.7715063 Pr(>|t|): 8.956991e-13 *
## mass_c coeff: -14.4801559 Pr(>|t|): 0.3287944Morphology Effects
## glm flew_b ~ beak_c binomial data_male
## AIC: 496.7797
## (Intercept) coeff: 0.77863 Pr(>|t|): 7.884113e-13 *
## beak_c coeff: 0.559369 Pr(>|t|): 0.04417121 *## glm flew_b ~ thorax_c binomial data_male
## AIC: 500.4994
## (Intercept) coeff: 0.7704035 Pr(>|t|): 9.132743e-13 *
## thorax_c coeff: 0.3470788 Pr(>|t|): 0.5183982## glm flew_b ~ body_c binomial data_male
## AIC: 496.0185
## (Intercept) coeff: 0.7792955 Pr(>|t|): 7.83718e-13 *
## body_c coeff: 0.3439638 Pr(>|t|): 0.02711686 *## glm flew_b ~ wing_c binomial data_male
## AIC: 493.4751
## (Intercept) coeff: 0.7840123 Pr(>|t|): 7.251728e-13 *
## wing_c coeff: 0.495168 Pr(>|t|): 0.006665089 *Without Mass + Days From Start
# Remove any missing masses
data_male <- data_male %>%
filter(!is.na(mass), !is.na(wing))
data_male <- center_data(data_male)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthR = data_male$flew_b
A = data_male$host_c
B = data_male$days_from_start_c
C = data_male$wing_c
D = data_male$mass_c
data<-data.frame(R, A, B, C, D)
model_script = paste0(source_path,"generic models-binomial glm 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 484.7966 484.9322 485.7839 485.8981 485.9905 486.145 486.7854
## models 12 7 13 16 11 14 6
## probs 0.1786059 0.166895 0.1090184 0.1029697 0.09832031 0.09100878 0.06607356
## [,8]
## AICs 486.9662
## models 15
## probs 0.06036154
##
## m12 glm(formula = R ~ A * C + B, family = binomial, data = data)
## m7 glm(formula = R ~ A + B + C, family = binomial, data = data)
## m13 glm(formula = R ~ B * C + A, family = binomial, data = data)
## m16 glm(formula = R ~ A * C + B * C, family = binomial, data = data)
## m11 glm(formula = R ~ A * B + C, family = binomial, data = data)
## m14 glm(formula = R ~ A * B + A * C, family = binomial, data = data)
## m6 glm(formula = R ~ B + C, family = binomial, data = data)
## m15 glm(formula = R ~ A * B + B * C, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m7, m12, test="Chisq") # Adding A*C does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A + B + C
## Model 2: R ~ A * C + B
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 394 476.93
## 2 393 474.80 1 2.1356 0.1439anova(m6, m7, test="Chisq") # Adding A improves fit## Analysis of Deviance Table
##
## Model 1: R ~ B + C
## Model 2: R ~ A + B + C
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 395 480.79
## 2 394 476.93 1 3.8532 0.04965 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## glm flew_b ~ host_c + days_from_start_c + wing_c binomial data_male
## AIC: 484.9322
## (Intercept) coeff: 0.690237 Pr(>|t|): 2.552291e-08 *
## host_c coeff: -0.2450768 Pr(>|t|): 0.04811069 *
## days_from_start_c coeff: -0.0497262 Pr(>|t|): 0.002335063 *
## wing_c coeff: 0.5077837 Pr(>|t|): 0.006220775 *- Negative effect if from GRT
- Negative effect of days from start where the farther out the day from the start of the experiment (e.g. maybe a proxi for age), the less likely the male bug is to fly
- Positive effect wing length where the longer the wing the more likely the male bug is to fly
## Consider covariates
model_male_final <-glmer(flew_b~host_c + days_from_start_c + wing_c + (1|ID), family=binomial, data=data_male)
tidy_regression(model_male_final, is_color=output_col)## glmer flew_b ~ host_c + days_from_start_c + wing_c + (1 | ID) data_male binomial
## AIC: 461.6402 461.6402
## (Intercept) coeff: 1.1896769 Pr(>|t|): 8.63996e-05 *
## host_c coeff: -0.4623157 Pr(>|t|): 0.07102692 .
## days_from_start_c coeff: -0.0827028 Pr(>|t|): 0.0007555347 *
## wing_c coeff: 0.8533121 Pr(>|t|): 0.02535222 *## Changed the effect slightly and improved the model- Marginal negative effect if from GRT
- Negative effect of days from start where the farther out the day from the start of the experiment (e.g. maybe a proxi for age), the less likely the male bug is to fly
- Positive effect wing length where the longer the wing the more likely the male bug is to fly
With Mass + Days From Start
model_script = paste0(source_path,"generic models-binomial glm 4-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1]
## AICs 475.8103
## models 15
## probs 0.06605644
##
## m15 glm(formula = R ~ A + B + C + D, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m15, m35, test="Chisq") # Adding A*C does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ A + B + C + D
## Model 2: R ~ A * C + B + D
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 393 465.81
## 2 392 464.38 1 1.4288 0.232## glm flew_b ~ host_c + mass_c + days_from_start_c + wing_c binomial data_male
## AIC: 475.8103
## (Intercept) coeff: 0.6685714 Pr(>|t|): 1.159895e-07 *
## host_c coeff: -0.3331418 Pr(>|t|): 0.01025867 *
## mass_c coeff: -67.4777477 Pr(>|t|): 0.001017975 *
## days_from_start_c coeff: -0.040922 Pr(>|t|): 0.01456519 *
## wing_c coeff: 1.0036593 Pr(>|t|): 4.388538e-05 *- Strong negative effect if from GRT
- Strong negative effect of mass, that if weigh more then less likely to disperse
- negative effect of days from start where the farther out the day from the start of the experiment (e.g. maybe a proxi for age), the less likely the male bug is to fly
- positive effect of wing length where the longer th wing, the more likely the male bug is to fly
## Consider covariates
mass_male_glmer <-glmer(flew_b~host_c + mass_c + days_from_start_c + wing_c + (1|population) + (1|trial_type), family=binomial, data=data_male) # no error, but singular fit error ## boundary (singular) fit: see ?isSingulartidy_regression(mass_male_glmer, is_color=output_col)## glmer flew_b ~ host_c + mass_c + days_from_start_c + wing_c + (1 | population) + (1 | trial_type) data_male binomial
## AIC: 475.2368 475.2368
## (Intercept) coeff: 0.6347795 Pr(>|t|): 0.001638181 *
## host_c coeff: -0.2540031 Pr(>|t|): 0.2179676
## mass_c coeff: -69.4180012 Pr(>|t|): 0.0009737074 *
## days_from_start_c coeff: -0.0423158 Pr(>|t|): 0.01287823 *
## wing_c coeff: 1.1066491 Pr(>|t|): 2.8196e-05 *## Changed the effects slightly, a better fitting model, and made host not significantMorphology
data_male <- data_tested[data_tested$sex=="M",]
data_male <- data_male %>%
filter(!is.na(mass)) %>%
filter(!is.na(body))
data_male <- center_data(data_male)## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthR = data_male$flew_b
A = data_male$thorax_c
B = data_male$body_c
C = data_male$wing_c
data<-data.frame(R, A, B, C)
model_script = paste0(source_path,"generic models-binomial glm 3-FF.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 483.6585 484.2478 485.1861 485.4292 485.6202 486.4286 486.7925
## models 13 10 16 9 15 17 11
## probs 0.2532581 0.1886292 0.1179935 0.1044878 0.09497182 0.06339317 0.05284828
##
## m13 glm(formula = R ~ B * C + A, family = binomial, data = data)
## m10 glm(formula = R ~ B * C, family = binomial, data = data)
## m16 glm(formula = R ~ A * C + B * C, family = binomial, data = data)
## m9 glm(formula = R ~ A * C, family = binomial, data = data)
## m15 glm(formula = R ~ A * B + B * C, family = binomial, data = data)
## m17 glm(formula = R ~ A * B + A * C + B * C, family = binomial, data = data)
## m11 glm(formula = R ~ A * B + C, family = binomial, data = data)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m10, m13, test="Chisq") # Adding A does not improve fit## Analysis of Deviance Table
##
## Model 1: R ~ B * C
## Model 2: R ~ B * C + A
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 394 476.25
## 2 393 473.66 1 2.5892 0.1076morph_male <-glm(flew_b~body_c* wing_c, family=binomial, data=data_male)
tidy_regression(morph_male, is_color=output_col) ## glm flew_b ~ body_c * wing_c binomial data_male
## AIC: 484.2478
## (Intercept) coeff: 1.0621697 Pr(>|t|): 6.603942e-14 *
## body_c coeff: -0.9138665 Pr(>|t|): 0.1217577
## wing_c coeff: 1.2195143 Pr(>|t|): 0.08033274 .
## body_c:wing_c coeff: -0.7043743 Pr(>|t|): 0.0016459 *no effect of body length
marginal positive effect of wing length, where if longer wing then more likely to fly
negative effect of body*wing where if longer body and wing then less likely to fly
morph_male_glmer <- glmer(flew_b ~ body_c* wing_c + (1|ID), family=binomial, data=data_male)
tidy_regression(morph_male_glmer, is_color=output_col) # model improves fit and it converges## glmer flew_b ~ body_c * wing_c + (1 | ID) data_male binomial
## AIC: 465.3772 465.3772
## (Intercept) coeff: 1.7666191 Pr(>|t|): 7.532794e-07 *
## body_c coeff: -1.4734131 Pr(>|t|): 0.1642051
## wing_c coeff: 1.9167017 Pr(>|t|): 0.1248792
## body_c:wing_c coeff: -1.2218112 Pr(>|t|): 0.004369603 *no effect of body length
no effect of wing length
negative effect of body*wing where if longer body and wing then less likely to fly
glmer() A=thorax, B=body, C=wing, D=days_from_start, X=ID, Y=trial_type
R = data_male$flew_b
A = data_male$thorax_c
B = data_male$body_c
C = data_male$wing_c
D = data_male$days_from_start_c
X = data_male$ID # 3 models failed to converge when ID as RF
#Y = data_male$trial_type # more models fail if both ID and trial as RF
data<-data.frame(R, A, B, C, D, X)
model_script = paste0(source_path,"generic models-binomial glmer 1-RF + 4-FF-noD-interactions.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs 455.1421 455.2078 455.9558 457.041 457.1342 457.4628 457.8314
## models 27 24 22 33 32 25 26
## probs 0.210914 0.2040983 0.1404132 0.08161466 0.07789699 0.06609419 0.05497175
## [,8]
## AICs 457.9353
## models 20
## probs 0.05218774
##
## m27 R ~ B * C + A + D + (1 | X)
## m24 R ~ B * C + D + (1 | X)
## m22 R ~ A * C + D + (1 | X)
## m33 R ~ A * C + B * C + D + (1 | X)
## m32 R ~ A * B + B * C + D + (1 | X)
## m25 R ~ A * B + C + D + (1 | X)
## m26 R ~ A * C + B + D + (1 | X)
## m20 R ~ A * B + D + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m24, m27, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m24: R ~ B * C + D + (1 | X)
## m27: R ~ B * C + A + D + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m24 6 455.21 479.13 -221.60 443.21
## m27 7 455.14 483.05 -220.57 441.14 2.0657 1 0.1506anova(m22, m26, test="Chisq") # Adding B does not improve fit## Data: data
## Models:
## m22: R ~ A * C + D + (1 | X)
## m26: R ~ A * C + B + D + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m22 6 455.96 479.87 -221.98 443.96
## m26 7 457.83 485.74 -221.92 443.83 0.1245 1 0.7242morph_male_glmer2 <- glmer(flew_b ~ thorax_c * wing_c + days_from_start_c + (1 | ID) , family=binomial, data=data_male)
tidy_regression(morph_male_glmer2, is_color=output_col) ## glmer flew_b ~ thorax_c * wing_c + days_from_start_c + (1 | ID) data_male binomial
## AIC: 455.9558 455.9558
## (Intercept) coeff: 1.833075 Pr(>|t|): 7.571767e-06 *
## thorax_c coeff: -3.8151257 Pr(>|t|): 0.04417442 *
## wing_c coeff: 1.653975 Pr(>|t|): 0.01484067 *
## days_from_start_c coeff: -0.0814431 Pr(>|t|): 0.001113034 *
## thorax_c:wing_c coeff: -4.0709097 Pr(>|t|): 0.02784661 *glmer() A=wing2body, B=thorax, C=days_from_start, X=ID, Y=trial_type
# left off here
R = data_male$flew_b
A = data_male$wing2body_c
B = data_male$thorax_c
C = data_male$days_from_start_c
X = data_male$ID
Y = data_male$trial_type
data<-data.frame(R, A, B, C, X, Y)
model_script = paste0(source_path,"generic models-binomial glmer 2 RF + 3-FF-noCinteractions.R")
errors = catch_warnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5]
## AICs 461.9652 463.9597 463.9652 464.2025 465.9597
## models 5 7 24 9 26
## probs 0.4062747 0.1498759 0.1494601 0.1327392 0.05513628
##
## m5 R ~ A + C + (1 | X)
## m7 R ~ A + B + C + (1 | X)
## m24 R ~ A + C + (1 | X) + (1 | Y)
## m9 R ~ A * B + C + (1 | X)
## m26 R ~ A + B + C + (1 | X) + (1 | Y)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 6anova(m5, m7, test="Chisq") # Adding B does not improve fit## Data: data
## Models:
## m5: R ~ A + C + (1 | X)
## m7: R ~ A + B + C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m5 4 461.97 477.91 -226.98 453.97
## m7 5 463.96 483.89 -226.98 453.96 0.0056 1 0.9406anova(m5, m24, test="Chisq") # Adding Y does not improve fit## Data: data
## Models:
## m5: R ~ A + C + (1 | X)
## m24: R ~ A + C + (1 | X) + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m5 4 461.97 477.91 -226.98 453.97
## m24 5 463.97 483.90 -226.98 453.97 0 1 1morph_male_glmer3 <- glmer(flew_b ~ wing2body_c + days_from_start_c + (1 | ID), family=binomial, data=data_male)
tidy_regression(morph_male_glmer3, is_color=output_col) ## glmer flew_b ~ wing2body_c + days_from_start_c + (1 | ID) data_male binomial
## AIC: 461.9652 461.9652
## (Intercept) coeff: 1.42291 Pr(>|t|): 4.609553e-06 *
## wing2body_c coeff: 36.4422713 Pr(>|t|): 0.01117388 *
## days_from_start_c coeff: -0.0777844 Pr(>|t|): 0.00146109 *- strong positive effect of wing2body, where if you have more wing to body, then more likely to fly
- negative effect of days from start, where farther out in days from the start day, the less likely the bug will fly.
Summary
tidy_regression(nomass_male, is_color=output_col)## glm flew_b ~ host_c + days_from_start_c + wing_c binomial data_male
## AIC: 484.9322
## (Intercept) coeff: 0.690237 Pr(>|t|): 2.552291e-08 *
## host_c coeff: -0.2450768 Pr(>|t|): 0.04811069 *
## days_from_start_c coeff: -0.0497262 Pr(>|t|): 0.002335063 *
## wing_c coeff: 0.5077837 Pr(>|t|): 0.006220775 *tidy_regression(mass_male, is_color=output_col)## glm flew_b ~ host_c + mass_c + days_from_start_c + wing_c binomial data_male
## AIC: 475.8103
## (Intercept) coeff: 0.6685714 Pr(>|t|): 1.159895e-07 *
## host_c coeff: -0.3331418 Pr(>|t|): 0.01025867 *
## mass_c coeff: -67.4777477 Pr(>|t|): 0.001017975 *
## days_from_start_c coeff: -0.040922 Pr(>|t|): 0.01456519 *
## wing_c coeff: 1.0036593 Pr(>|t|): 4.388538e-05 *tidy_regression(mass_male_glmer, is_color=output_col)## glmer flew_b ~ host_c + mass_c + days_from_start_c + wing_c + (1 | population) + (1 | trial_type) data_male binomial
## AIC: 475.2368 475.2368
## (Intercept) coeff: 0.6347795 Pr(>|t|): 0.001638181 *
## host_c coeff: -0.2540031 Pr(>|t|): 0.2179676
## mass_c coeff: -69.4180012 Pr(>|t|): 0.0009737074 *
## days_from_start_c coeff: -0.0423158 Pr(>|t|): 0.01287823 *
## wing_c coeff: 1.1066491 Pr(>|t|): 2.8196e-05 *tidy_regression(morph_male, is_color=output_col) ## glm flew_b ~ body_c * wing_c binomial data_male
## AIC: 484.2478
## (Intercept) coeff: 1.0621697 Pr(>|t|): 6.603942e-14 *
## body_c coeff: -0.9138665 Pr(>|t|): 0.1217577
## wing_c coeff: 1.2195143 Pr(>|t|): 0.08033274 .
## body_c:wing_c coeff: -0.7043743 Pr(>|t|): 0.0016459 *tidy_regression(morph_male_glmer, is_color=output_col) ## glmer flew_b ~ body_c * wing_c + (1 | ID) data_male binomial
## AIC: 465.3772 465.3772
## (Intercept) coeff: 1.7666191 Pr(>|t|): 7.532794e-07 *
## body_c coeff: -1.4734131 Pr(>|t|): 0.1642051
## wing_c coeff: 1.9167017 Pr(>|t|): 0.1248792
## body_c:wing_c coeff: -1.2218112 Pr(>|t|): 0.004369603 *tidy_regression(morph_male_glmer2, is_color=output_col)## glmer flew_b ~ thorax_c * wing_c + days_from_start_c + (1 | ID) data_male binomial
## AIC: 455.9558 455.9558
## (Intercept) coeff: 1.833075 Pr(>|t|): 7.571767e-06 *
## thorax_c coeff: -3.8151257 Pr(>|t|): 0.04417442 *
## wing_c coeff: 1.653975 Pr(>|t|): 0.01484067 *
## days_from_start_c coeff: -0.0814431 Pr(>|t|): 0.001113034 *
## thorax_c:wing_c coeff: -4.0709097 Pr(>|t|): 0.02784661 *tidy_regression(morph_male_glmer3, is_color=output_col) ## glmer flew_b ~ wing2body_c + days_from_start_c + (1 | ID) data_male binomial
## AIC: 461.9652 461.9652
## (Intercept) coeff: 1.42291 Pr(>|t|): 4.609553e-06 *
## wing2body_c coeff: 36.4422713 Pr(>|t|): 0.01117388 *
## days_from_start_c coeff: -0.0777844 Pr(>|t|): 0.00146109 *