Binomial Modeling: flight_yes

Testing glmer() and Covariates

Without Mass

With Mass

Morphology

Morphology + Mass + Sex

Cleaning the Data

rm(list=ls())
output_col = FALSE # Recommend changing this to TRUE if working in Base R or RStudio, and FALSE if generating an html
source("src/clean_flight_data.R") # Script that loads and cleans up the data

## 
## Attaching package: 'lubridate'

## The following object is masked from 'package:base':
## 
##     date

source("src/regression_output.R") # A script that cleans up regression outputs and prints in color or black and white
source("src/center_flight_data.R")
source("get_warnings.R")

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

d <- data_tested %>%
   group_by(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, wing2body_s) %>%
   summarise_all(funs(list(na.omit(.))))

## Warning: funs() is soft deprecated as of 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))
## This warning is displayed once per session.

d$num_flew <- 0
d$num_notflew <- 0
d$average_mass <- 0

for(row in 1:length(d$flew_b)){
  
  n_flew <- sum(d$flew_b[[row]] == 1) # total number of times did fly among trails
  d$num_flew[[row]] <- n_flew 
  
  n_notflew <- sum(d$flew_b[[row]] == 0) # total number of times did not fly among trails
  d$num_notflew[[row]] <- n_notflew
  
  avg_mass <- mean(d$mass[[row]])
  d$average_mass[[row]] <- avg_mass

}

d <- select(d, -filename, -channel_letter, -set_number)
d <- center_data(d, is_not_binded = FALSE)
d

## # A tibble: 333 x 63
## # 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 [333]
##    ID    sex   population site  host_plant latitude longitude total_eggs  beak
##    <fct> <fct> <fct>      <fct> <fct>         <dbl>     <dbl>      <dbl> <dbl>
##  1 1     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  5.5 
##  2 2     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  5.64
##  3 3     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  5.75
##  4 4     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  6.21
##  5 5     F     Plantatio… Areg… C. corind…     25.0     -80.6        209  7.47
##  6 6     M     Plantatio… Foun… C. corind…     25.0     -80.6         NA  5.76
##  7 7     F     Plantatio… Foun… C. corind…     25.0     -80.6          8  8.19
##  8 8     F     Plantatio… Foun… C. corind…     25.0     -80.6         92  8.39
##  9 9     M     Plantatio… Foun… C. corind…     25.0     -80.6         NA  6.1 
## 10 10    M     Plantatio… Foun… C. corind…     25.0     -80.6         NA  6.09
## # … with 323 more rows, and 54 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>,
## #   total_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>, 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>, average_mass <dbl>, average_mass_c <dbl>,
## #   average_mass_s <dbl>, trial_type_og <list>, num_flew <dbl>,
## #   num_notflew <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 ?isSingular

tidy_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:  688.8968 688.8968 
## (Intercept)      coeff:  -0.0042496  Pr(>|t|):  0.9662782
## sex_c            coeff:  -0.6640022  Pr(>|t|):  3.955971e-11 *
## host_c           coeff:  -0.0679509  Pr(>|t|):  0.4990415
## sex_c:host_c     coeff:  0.1698073   Pr(>|t|):  0.09116195 .

model <- 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:  684.4669 684.4669 
## (Intercept)      coeff:  -0.0080055  Pr(>|t|):  0.9538386
## sex_c            coeff:  -0.7069343  Pr(>|t|):  3.332778e-11 *
## host_c           coeff:  -0.0475807  Pr(>|t|):  0.7311514
## sex_c:host_c     coeff:  0.1894487   Pr(>|t|):  0.07068904 .

getME(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) ~ average_mass_c binomial d 
## AIC:  690.0321 
## (Intercept)      coeff:  0.2314802   Pr(>|t|):  0.007290253 *
## average_mass_c   coeff:  -31.1658278 Pr(>|t|):  3.565029e-14 *

## glm cbind(num_flew, num_notflew) ~ total_eggs binomial d 
## AIC:  184.4609 
## (Intercept)  coeff:  0.0580002   Pr(>|t|):  0.8111338
## total_eggs   coeff:  -0.0117451  Pr(>|t|):  1.730675e-05 *

## glm cbind(num_flew, num_notflew) ~ beak_c binomial d 
## AIC:  734.2146 
## (Intercept)  coeff:  0.23902 Pr(>|t|):  0.003973372 *
## beak_c       coeff:  -0.4085962  Pr(>|t|):  9.565269e-07 *

## glm cbind(num_flew, num_notflew) ~ thorax_c binomial d 
## AIC:  739.1983 
## (Intercept)  coeff:  0.2385862   Pr(>|t|):  0.003890438 *
## thorax_c     coeff:  -1.2298779  Pr(>|t|):  1.111691e-05 *

## glm cbind(num_flew, num_notflew) ~ body_c binomial d 
## AIC:  750.4139 
## (Intercept)  coeff:  0.2389879   Pr(>|t|):  0.003520911 *
## body_c       coeff:  -0.2303836  Pr(>|t|):  0.003002383 *

## glm cbind(num_flew, num_notflew) ~ wing_c binomial d 
## AIC:  757.3505 
## (Intercept)  coeff:  0.2364567   Pr(>|t|):  0.003681042 *
## wing_c       coeff:  -0.1417783  Pr(>|t|):  0.1546876

Without Mass

R1 = d$num_flew
R2 = d$num_notflew
A = d$host_c
B = d$sex_c
C = d$sym_dist
D = d$average_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 average_mass (down below). 

data<-data.frame(R1, R2, A, B, C, D, X)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 3-FF.R"))

##        [,1]      [,2]      [,3]      [,4]       [,5]       [,6]      
## AICs   688.8968  688.95    689.707   690.5014   690.7554   691.6281  
## models 8         2         4         6          11         15        
## probs  0.2170513 0.2113509 0.1447477 0.09730022 0.08569733 0.05539398
##        [,7]      
## AICs   691.7067  
## models 7         
## probs  0.05325955
## 
## m8   cbind(R1, R2) ~ A * B + (1 | X)
## m2   cbind(R1, R2) ~ B + (1 | X)
## m4   cbind(R1, R2) ~ A + B + (1 | X)
## m6   cbind(R1, R2) ~ B + C + (1 | X)
## m11  cbind(R1, R2) ~ A * B + C + (1 | X)
## m15  cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m7   cbind(R1, R2) ~ A + B + C + (1 | X)

length(errors$warnings)

## [1] 0

anova(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)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m4  4 689.71 704.94 -340.85   681.71                           
## m8  5 688.90 707.94 -339.45   678.90 2.8103      1    0.09366 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(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)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2  3 688.95 700.37 -341.47   682.95                         
## m4  4 689.71 704.94 -340.85   681.71 1.2429      1     0.2649

anova(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)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2  3 688.95 700.37 -341.47   682.95                         
## m6  4 690.50 705.73 -341.25   682.50 0.4486      1      0.503

nomass_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:  688.95 688.95 
## (Intercept)  coeff:  -0.0019532  Pr(>|t|):  0.9872246
## sex_c        coeff:  -0.7440024  Pr(>|t|):  2.652377e-16 *

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:  684.2778 684.2778 
## (Intercept)  coeff:  0.0077662   Pr(>|t|):  0.9513932
## sex_c        coeff:  -0.7933747  Pr(>|t|):  3.26677e-16 *

strong negative effect of sex where if F then less number of times will fly

With Mass

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1-RF + 4-FF.R"))

##        [,1]       [,2]       [,3]      
## AICs   681.2078   681.549    682.5085  
## models 63         26         85        
## probs  0.09984206 0.08418323 0.05210509
## 
## m63  cbind(R1, R2) ~ A * D + C * D + B + (1 | X)
## m26  cbind(R1, R2) ~ A * D + B + (1 | X)
## m85  cbind(R1, R2) ~ A * D + B * D + C * D + (1 | X)

length(errors$warnings)

## [1] 1

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m63  8 681.21 711.67 -332.60   665.21                         
## m85  9 682.51 716.78 -332.25   664.51 0.6993      1      0.403

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m26  6 681.55 704.40 -334.77   669.55                         
## m36  7 683.36 710.02 -334.68   669.36 0.1877      1     0.6648

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m12  5 684.77 703.81 -337.39   674.77                           
## m26  6 681.55 704.40 -334.77   669.55 5.2227      1    0.02229 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wanted 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:  760.9329 760.9329 
## (Intercept)  coeff:  0.1891344   Pr(>|t|):  0.1982647
## sym_dist     coeff:  0.0156086   Pr(>|t|):  0.8614592

# 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 * average_mass_c +  sex_c + (1|population), data=d, family=binomial)

## boundary (singular) fit: see ?isSingular

tidy_regression(mass_model_all, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ host_c * average_mass_c + sex_c + (1 | population) d binomial 
## AIC:  681.549 681.549 
## (Intercept)              coeff:  0.0708971   Pr(>|t|):  0.493529
## host_c                   coeff:  -0.1284394  Pr(>|t|):  0.1898795
## average_mass_c           coeff:  -10.6860919 Pr(>|t|):  0.136257
## sex_c                    coeff:  -0.440198   Pr(>|t|):  0.004480143 *
## host_c:average_mass_c    coeff:  10.488122   Pr(>|t|):  0.01997165 *

mass_model_all <- glmer(cbind(num_flew, num_notflew) ~ host_c * average_mass_c + sym_dist * average_mass_c + sex_c + (1|population), data=d, family=binomial)
tidy_regression(mass_model_all, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ host_c * average_mass_c + sym_dist * average_mass_c + sex_c + (1 | population) d binomial 
## AIC:  681.2078 681.2078 
## (Intercept)              coeff:  0.0662598   Pr(>|t|):  0.6987177
## host_c                   coeff:  -0.1455773  Pr(>|t|):  0.3316777
## average_mass_c           coeff:  -2.5133936  Pr(>|t|):  0.7606884
## sym_dist                 coeff:  -0.0428036  Pr(>|t|):  0.7456766
## sex_c                    coeff:  -0.4042136  Pr(>|t|):  0.009820863 *
## host_c:average_mass_c    coeff:  16.0913743  Pr(>|t|):  0.002395689 *
## average_mass_c:sym_dist  coeff:  -11.9487833 Pr(>|t|):  0.06336045 .

Now there’s a conflicting interaction between average mass and host*average_mass…and there’s

no effect of host
no effect of average_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, average_mass), data=d , col=rangi2)

How does num_flew or num_notflew vary with host plant?

coplot(num_flew ~ average_mass_c | host_c, data=d)

coplot(num_notflew~ average_mass_c | host_c, data=d)

Plot the interaction between host*average_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 ~ average_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 ~ average_mass, col=~host_plant, data=d[d$average_mass_c < 0,], title="(a)   avgmass_c <0") + 
  geom_smooth(method='lm')
p2 <- gf_point(num_flew ~ average_mass, col=~host_plant, data=d[d$average_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 ~ average_mass, col=~host_plant, data=d[d$average_mass_c < 0,], title="(a)   avgmass_c <0") + 
  geom_smooth(method='lm')
p4 <- gf_point(num_notflew ~ average_mass, col=~host_plant, data=d[d$average_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)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 3-FF.R"))

##        [,1]      [,2]      [,3]     [,4]      
## AICs   696.4932  697.7959  698.4702 699.7618  
## models 15        16        17       13        
## probs  0.4471192 0.2331006 0.166388 0.08722792
## 
## m15  cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m16  cbind(R1, R2) ~ A * C + B * C + (1 | X)
## m17  cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)
## m13  cbind(R1, R2) ~ B * C + A + (1 | X)

length(errors$warnings)

## [1] 1

anova(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)
##     Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## m15  7 696.49 723.15 -341.25   682.49                        
## m17  8 698.47 728.94 -341.24   682.47 0.023      1     0.8795

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m13  6 699.76 722.61 -343.88   687.76                           
## m15  7 696.49 723.15 -341.25   682.49 5.2686      1    0.02171 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

morph_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:  696.4932 696.4932 
## (Intercept)      coeff:  0.2804294   Pr(>|t|):  0.07044802 .
## thorax_c         coeff:  -2.549434   Pr(>|t|):  0.003596569 *
## body_c           coeff:  -1.2548354  Pr(>|t|):  0.007228665 *
## wing_c           coeff:  2.2612361   Pr(>|t|):  6.326959e-06 *
## thorax_c:body_c  coeff:  1.3630407   Pr(>|t|):  0.02490653 *
## body_c:wing_c    coeff:  -0.6384384  Pr(>|t|):  0.003529522 *

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)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1-RF + 2-FF.R"))

##        [,1]      [,2]     
## AICs   699.1221  700.8564 
## models 4         3        
## probs  0.7041422 0.2958542
## 
## m4   cbind(R1, R2) ~ A * B + (1 | X)
## m3   cbind(R1, R2) ~ A + B + (1 | X)

length(errors$warnings)

## [1] 0

anova(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)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m3  4 700.86 716.09 -346.43   692.86                           
## m4  5 699.12 718.16 -344.56   689.12 3.7342      1    0.05331 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

morph_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:  699.1221 699.1221 
## (Intercept)              coeff:  0.1644896   Pr(>|t|):  0.2796688
## wing2body_c              coeff:  30.4358638  Pr(>|t|):  3.405945e-08 *
## thorax_c                 coeff:  -1.4513019  Pr(>|t|):  8.423364e-07 *
## wing2body_c:thorax_c     coeff:  -37.3187734 Pr(>|t|):  0.05638751 .

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$average_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$average_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$average_mass_c
X = d$population

data<-data.frame(R1, R2, A, B, C, X)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 3-FF.R"))

##        [,1]      [,2]     [,3]      
## AICs   660.5244  662.4846 665.1787  
## models 15        17       11        
## probs  0.6458782 0.242385 0.06302092
## 
## 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)

length(errors$warnings) # 13 models failed to converge

## [1] 13

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## m11  6 665.18 688.03 -326.59   653.18                            
## m15  7 660.52 687.18 -323.26   646.52 6.6543      1   0.009892 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m15  7 660.52 687.18 -323.26   646.52                         
## m17  8 662.48 692.95 -323.24   646.48 0.0398      1     0.8418

multi_model_all <- glmer(cbind(num_flew, num_notflew) ~ thorax_c * wing2body_c + wing2body_c * average_mass + (1|population), data=d, family=binomial) 
tidy_regression(multi_model_all, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ thorax_c * wing2body_c + wing2body_c * average_mass + (1 | population) d binomial 
## AIC:  660.5244 660.5244 
## (Intercept)                  coeff:  2.0903463   Pr(>|t|):  2.37968e-07 *
## thorax_c                     coeff:  0.742338    Pr(>|t|):  0.1527069
## wing2body_c                  coeff:  -39.7258582 Pr(>|t|):  9.264061e-07 *
## average_mass                 coeff:  -36.1155189 Pr(>|t|):  7.888075e-07 *
## thorax_c:wing2body_c         coeff:  -108.5821295    Pr(>|t|):  2.376879e-07 *
## wing2body_c:average_mass     coeff:  1062.7542387    Pr(>|t|):  7.274993e-27 *

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)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 3-FF.R"))

##        [,1]      [,2]      [,3]      [,4]      [,5]       [,6]       [,7]      
## AICs   676.149   676.8094  677.0927  678.0144  678.6671   678.7985   679.0892  
## models 11        6         15        14        10         7          17        
## probs  0.2652567 0.1906682 0.1654839 0.1043776 0.07531222 0.07052579 0.06098305
## 
## m11  cbind(R1, R2) ~ A * B + C + (1 | X)
## m6   cbind(R1, R2) ~ B + C + (1 | X)
## m15  cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m14  cbind(R1, R2) ~ A * B + A * C + (1 | X)
## m10  cbind(R1, R2) ~ B * C + (1 | X)
## m7   cbind(R1, R2) ~ A + B + C + (1 | X)
## m17  cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)

length(errors$warnings) # no models failed

## [1] 0

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m11  6 676.15 699.00 -332.07   664.15                         
## m15  7 677.09 703.75 -331.55   663.09 1.0564      1      0.304

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m11  6 676.15 699.00 -332.07   664.15                         
## m14  7 678.01 704.67 -332.01   664.01 0.1346      1     0.7137

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m7   5 678.80 697.84 -334.40   668.80                           
## m11  6 676.15 699.00 -332.07   664.15 4.6494      1    0.03106 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

multi_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:  676.149 676.149 
## (Intercept)              coeff:  -0.0074703  Pr(>|t|):  0.9588266
## thorax_c                 coeff:  0.0512288   Pr(>|t|):  0.902632
## wing2body_c              coeff:  16.9221328  Pr(>|t|):  0.00567254 *
## sex_c                    coeff:  -0.6784392  Pr(>|t|):  9.704959e-07 *
## thorax_c:wing2body_c     coeff:  -40.4953172 Pr(>|t|):  0.03486006 *

no 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$average_mass_s
X = d$population

data<-data.frame(R1, R2, A, B, C, D, X)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1-RF + 4-FF.R"))

##        [,1]      [,2]      [,3]       [,4]       [,5]       [,6]      
## AICs   657.3179  658.2776  658.5987   658.8669   659.0584   659.2736  
## models 74        106       58         94         100        97        
## probs  0.1766309 0.1093142 0.09310048 0.08141795 0.07398394 0.06643595
##        [,7]      
## AICs   659.4118  
## models 112       
## probs  0.06200052
## 
## m74  cbind(R1, R2) ~ A * B + B * C + B * D + (1 | X)
## m106     cbind(R1, R2) ~ A * B + A * C + A * D + B * C + B * D + (1 | 
##     X)
## m58  cbind(R1, R2) ~ A * B + B * D + C + (1 | X)
## m94  cbind(R1, R2) ~ A * B + A * C + B * C + B * D + (1 | X)
## m100     cbind(R1, R2) ~ A * B + B * C + B * D + C * D + (1 | X)
## m97  cbind(R1, R2) ~ A * B + A * D + B * C + B * D + (1 | X)
## m112     cbind(R1, R2) ~ A * B + A * C + A * D + B * C + B * D + C * D + 
##     (1 | X)

length(errors$warnings) # 108 models failed if use centered variables but only 7 fail if use standardized variables. Also, all the top models have multiple interactions b/c I'm guessing, they're correlated in some ways.

## [1] 7

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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m58  8 658.60 689.06 -321.30   642.60                           
## m74  9 657.32 691.59 -319.66   639.32 3.2808      1     0.0701 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m74  9 657.32 691.59 -319.66   639.32                         
## m94 10 658.87 696.95 -319.43   638.87 0.4511      1     0.5018

multi_model <- glmer(cbind(num_flew, num_notflew) ~ thorax_s * wing2body_s + wing2body_s * average_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 * average_mass_s + sex_c + (1 | population) d binomial 
## AIC:  658.5987 658.5987 
## (Intercept)                  coeff:  0.1135378   Pr(>|t|):  0.4745605
## thorax_s                     coeff:  0.2865295   Pr(>|t|):  0.07625582 .
## wing2body_s                  coeff:  0.224649    Pr(>|t|):  0.05382945 .
## average_mass_s               coeff:  -0.6417996  Pr(>|t|):  0.001958101 *
## sex_c                        coeff:  -0.3372454  Pr(>|t|):  0.04717871 *
## thorax_s:wing2body_s         coeff:  -0.5894181  Pr(>|t|):  0.000514089 *
## wing2body_s:average_mass_s   coeff:  0.4501488   Pr(>|t|):  0.01478315 *

Summary

tidy_regression(nomass_model_all, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ sex_c + (1 | site) d binomial 
## AIC:  684.2778 684.2778 
## (Intercept)  coeff:  0.0077662   Pr(>|t|):  0.9513932
## sex_c        coeff:  -0.7933747  Pr(>|t|):  3.26677e-16 *

tidy_regression(mass_model_all, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ host_c * average_mass_c + sym_dist * average_mass_c + sex_c + (1 | population) d binomial 
## AIC:  681.2078 681.2078 
## (Intercept)              coeff:  0.0662598   Pr(>|t|):  0.6987177
## host_c                   coeff:  -0.1455773  Pr(>|t|):  0.3316777
## average_mass_c           coeff:  -2.5133936  Pr(>|t|):  0.7606884
## sym_dist                 coeff:  -0.0428036  Pr(>|t|):  0.7456766
## sex_c                    coeff:  -0.4042136  Pr(>|t|):  0.009820863 *
## host_c:average_mass_c    coeff:  16.0913743  Pr(>|t|):  0.002395689 *
## average_mass_c:sym_dist  coeff:  -11.9487833 Pr(>|t|):  0.06336045 .

Why this host*average_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:  696.4932 696.4932 
## (Intercept)      coeff:  0.2804294   Pr(>|t|):  0.07044802 .
## thorax_c         coeff:  -2.549434   Pr(>|t|):  0.003596569 *
## body_c           coeff:  -1.2548354  Pr(>|t|):  0.007228665 *
## wing_c           coeff:  2.2612361   Pr(>|t|):  6.326959e-06 *
## thorax_c:body_c  coeff:  1.3630407   Pr(>|t|):  0.02490653 *
## body_c:wing_c    coeff:  -0.6384384  Pr(>|t|):  0.003529522 *

tidy_regression(morph_model_all2, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ wing2body_c * thorax_c + (1 | population) d binomial 
## AIC:  699.1221 699.1221 
## (Intercept)              coeff:  0.1644896   Pr(>|t|):  0.2796688
## wing2body_c              coeff:  30.4358638  Pr(>|t|):  3.405945e-08 *
## thorax_c                 coeff:  -1.4513019  Pr(>|t|):  8.423364e-07 *
## wing2body_c:thorax_c     coeff:  -37.3187734 Pr(>|t|):  0.05638751 .

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

## glmer cbind(num_flew, num_notflew) ~ thorax_c * wing2body_c + wing2body_c * average_mass + (1 | population) d binomial 
## AIC:  660.5244 660.5244 
## (Intercept)                  coeff:  2.0903463   Pr(>|t|):  2.37968e-07 *
## thorax_c                     coeff:  0.742338    Pr(>|t|):  0.1527069
## wing2body_c                  coeff:  -39.7258582 Pr(>|t|):  9.264061e-07 *
## average_mass                 coeff:  -36.1155189 Pr(>|t|):  7.888075e-07 *
## thorax_c:wing2body_c         coeff:  -108.5821295    Pr(>|t|):  2.376879e-07 *
## wing2body_c:average_mass     coeff:  1062.7542387    Pr(>|t|):  7.274993e-27 *

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:  676.149 676.149 
## (Intercept)              coeff:  -0.0074703  Pr(>|t|):  0.9588266
## thorax_c                 coeff:  0.0512288   Pr(>|t|):  0.902632
## wing2body_c              coeff:  16.9221328  Pr(>|t|):  0.00567254 *
## sex_c                    coeff:  -0.6784392  Pr(>|t|):  9.704959e-07 *
## thorax_c:wing2body_c     coeff:  -40.4953172 Pr(>|t|):  0.03486006 *

Huge coefficients! I’m sure this is all due to collinearity. Thoughts?

All Trials - Old Method

Testing glmer() and Covariates

Without Mass

With Mass

Morphology

Cleaning the Data

rm(list=ls())
output_col = FALSE # Recommend changing this to TRUE if working in Base R or RStudio, and FALSE if generating an html
source("src/clean_flight_data.R") # Script that loads and cleans up the data
source("src/regression_output.R") # A script that cleans up regression outputs and prints in color or black and white
source("src/center_flight_data.R")
source("get_warnings.R")

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

Testing Models and Covariates

# 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.151

getME(model, "lower")

## [1] 0 0

Experimental 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.6770982

Biological 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.1546876

Without Mass (Test)

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

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

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2 RF + 3-FF.R"))

##        [,1]      [,2]      [,3]      [,4]      [,5]       [,6]      
## AICs   707.3685  708.3419  709.0323  709.1583  709.3669   710.2419  
## models 46        40        42        44        49         53        
## probs  0.2655926 0.1632542 0.1155923 0.1085372 0.09778638 0.06313657
## 
## m46  R ~ A * B + (1 | X) + (1 | Y)
## m40  R ~ B + (1 | X) + (1 | Y)
## m42  R ~ A + B + (1 | X) + (1 | Y)
## m44  R ~ B + C + (1 | X) + (1 | Y)
## m49  R ~ A * B + C + (1 | X) + (1 | Y)
## m53  R ~ A * B + B * C + (1 | X) + (1 | Y)

length(errors$warnings) # 11 models did not converge

## [1] 11

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi 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 ' ' 1

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi 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.2525

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

# source("src/compare_models.R")
# model_comparisonsAIC("src/generic models-binomial glmer 2-RF + 4-FF.R")
# too many models failed to even keep running

Morphology

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

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

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2 RF + 3-FF.R"))

##        [,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)) # 7 models did not converge

## Number of models that failed to converge: 7

anova(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)
##     Df    AIC    BIC  logLik deviance Chisq Chi 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 ' ' 1

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

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2-RF + 2-FF.R"))

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

anova(m13, m14, test="Chisq") # Adding A*B doe snot improve fit

## Data: data
## Models:
## m13: R ~ A + B + (1 | X) + (1 | Y)
## m14: R ~ A * B + (1 | X) + (1 | Y)
##     Df    AIC    BIC  logLik deviance  Chisq Chi 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.1063

morph_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?

Delta Flight Response (T1 vs. T2)

Introduction

Clean the Data

Flew Diff ~ Mass Diff

Flew Diff ~ Egg Diff

Introduction

Some graphs I generated illustrated how changes in mass or egg-laying lead to changes in flight response, speed, and distance. Below, I run analyses on changes in flight response due to changes in mass or egg-laying that will help interpret and describe those graphs. Other graphs such as speed and distance changes will be found in the speed and distance scripts.

To perfom these analyses, a new variable will be created called “flew_diff” which calculates the flight response differences between T1 and T2. If flew_diff = -1, then the bug flew in T1 but not T2. If flew_diff = 0, then the bug either did not fly in either or flew in both T1 and T2 (no flight response). If flew_diff = 1, then the bug flew in T2 but not T1. Since this response variable is no longer binomial, I will need to use multicategorical logit models (refer to Chapter 6 of Introduction to Categorical Data Analysis by Alan Agresti). Below I’ve written out notes on performing these analyses.

What we are interested in is a multicategorical, nomial response variable. When the response variable is nominal, there is no natural order among the response variable categories (unordered categories). When the response variable is multicategorical, its multicategory models assume that the counts in the categories of \(Y\) have a multinomial distribution.

Let J = number of categories for \(Y\).

\(\pi_1,...\pi_J\) = the response probabilities where \(\sum_j \pi_j = 1\).

With \(n\) independent observations, the probability distribution for the number of outcomes of the J types is the multinomial. It specifices the probability for each possible way the \(n\) observations can fall in the J categories. Multicategory logit models simultaneously use all pairs of categories by specifying the odds of outocme in one category instead of another.

Logit models for nomial response variables pair each category with a baseline category. The choice of the baseline category is arbitrary. When the last category (J) is the baseline, the baseline-category logits are

\(\log (\frac{\pi_j}{\pi_J}), j = 1, ..., J - 1\).

Given that the response falls in category j or category J, this is the log odds that the response is j. For J = 3, for instance, the model uses \(\log(\pi_1/\pi_3)\) and \(\log(\pi_2/\pi_3)\). The baseline-category logit model with a predictor x is

\(\log \frac{\pi_j}{\pi_J} = \alpha_J + \beta_j x, j = 1, ..., J - 1\)

The model has J - 1 equations, with separate parameters for each. The effects vary according to the category paried with the baseline. When J = 2, this model simplifies to a single equaiton for \(\log(\pi_1/\pi_2) = logit(\pi_1)\), resulting in ordinary logistic regression for binary responses.

So how do these pair of categories determine equations for all other pairs of categories? Here is an arbitrary pair of categories a and b that follows the general equation above,

\(\log (\frac{\pi_a}{\pi_b}) = \log ( \frac{\pi_a/pi_J}{\pi_b/pi_J}) = \log (\frac{\pi_a}{\pi_J}) - \log (\frac{\pi_b}{\pi_J})\)

\(= (\alpha_a + \beta_a x) - (\alpha_b + \beta_b x)\)

\(= (\alpha_a - \alpha_b) + (\beta_a - \beta_b) x\)

So, the equaiton for categories a and b has the form \(\alpha + \beta x\) with intercept parameter \(\alpha = (\alpha_a - \alpha_b)\) and with slope parameter \(\beta = (\beta_a - \beta_b)\).

Let’s apply it to our data now:

Cleaning the Data

rm(list=ls())
output_col = FALSE # Recommend changing this to TRUE if working in Base R or RStudio, and FALSE if generating an html
source("src/clean_flight_data.R") # Script that loads and cleans up the data
source("src/regression_output.R") # A script that cleans up regression outputs and prints in color or black and white
source("src/center_flight_data.R")
source("get_warnings.R")

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

d <- data_tested %>%
   group_by(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, wing2body_s) %>%
   summarise_all(funs(list(na.omit(.))))

d$num_flew <- 0
d$num_notflew <- 0
d$average_mass <- 0

d$mass_diff <- 0
d$flew_diff <- 0
d$dist_diff <- 0
d$speed_diff <- 0

for(row in 1:length(d$flew_b)){
  
  n_flew <- sum(d$flew_b[[row]] == 1) # total number of times did fly among trails
  d$num_flew[[row]] <- n_flew 
  
  n_notflew <- sum(d$flew_b[[row]] == 0) # total number of times did not fly among trails
  d$num_notflew[[row]] <- n_notflew
  
  avg_mass <- mean(d$mass[[row]])
  d$average_mass[[row]] <- avg_mass
  
  # mass, flight, distance, and speed changes between T1 and T2
  
  d$mass_diff[[row]] <- d$mass[[row]][2] - d$mass[[row]][1]  # T2 - T1
  d$flew_diff[[row]] <- d$flew_b[[row]][2] - d$flew_b[[row]][1]  # T2 - T1
  d$dist_diff[[row]] <- d$distance[[row]][2] - d$distance[[row]][1]  # T2 - T1
  d$speed_diff[[row]] <- d$average_speed[[row]][2] - d$average_speed[[row]][1]  # T2 - T1
  d$egg_diff[[row]] <- d$eggs_b[[row]][2] - d$eggs_b[[row]][1]  # T2 - T1

}

## Warning: Unknown or uninitialised column: `egg_diff`.

d <- select(d, -filename, -channel_letter, -set_number)

d # NA's generated are good because that means it's accounted only for bugs that flew in both trials

## # A tibble: 333 x 68
## # 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 [333]
##    ID    sex   population site  host_plant latitude longitude total_eggs  beak
##    <fct> <fct> <fct>      <fct> <fct>         <dbl>     <dbl>      <dbl> <dbl>
##  1 1     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  5.5 
##  2 2     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  5.64
##  3 3     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  5.75
##  4 4     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  6.21
##  5 5     F     Plantatio… Areg… C. corind…     25.0     -80.6        209  7.47
##  6 6     M     Plantatio… Foun… C. corind…     25.0     -80.6         NA  5.76
##  7 7     F     Plantatio… Foun… C. corind…     25.0     -80.6          8  8.19
##  8 8     F     Plantatio… Foun… C. corind…     25.0     -80.6         92  8.39
##  9 9     M     Plantatio… Foun… C. corind…     25.0     -80.6         NA  6.1 
## 10 10    M     Plantatio… Foun… C. corind…     25.0     -80.6         NA  6.09
## # … with 323 more rows, and 59 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>,
## #   total_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>, eggs_b <list>,
## #   minute_duration <list>, minute_duration_c <list>, min_from_IncStart <list>,
## #   min_from_IncStart_c <list>, min_from_IncStart_s <list>,
## #   days_from_start <list>, days_from_start_c <list>, mass_c <list>,
## #   mass_s <list>, average_mass <dbl>, average_mass_c <list>,
## #   average_mass_s <list>, trial_type_og <list>, num_flew <dbl>,
## #   num_notflew <dbl>, mass_diff <dbl>, flew_diff <dbl>, dist_diff <dbl>,
## #   speed_diff <dbl>, egg_diff <dbl>

### AB: NOTE this includes also bugs that only flew in T1 so if we want to constrict it to bugs that flew in T1 and T2 we would need to adjust for that. If want to do that just uncomment the following:

# Filter out bugs that ONLY flew in T1:
# rows_remove <- c()
# for (row in 1:nrow(d)){
#   if (length(d$trial_type[[row]]) < 2) {
#     rows_remove <- c(rows_remove, row)
#   }
# }
# d <- d[-rows_remove, ]
# d

Below I used the multinom function from the nnet package to estimate a multinomial logistic regression model. There are other functions in other R packages capable of multinomial regression. I chose the multinom function because it does not require the data to be reshaped (as the mlogit package does) and to mirror the example code found in Hilbe’s Logistic Regression Models.

First, I need to choose the level of our outcome that we wish to use as our baseline and specify this in the relevel function. Then, I ran our model using multinom. The multinom package does not include p-value calculation for the regression coefficients, so I calculate p-values using Wald tests (here z-tests).

Flew Diff ~ Mass Diff

df <- d %>%
  filter(!is.na(mass_diff), !is.na(flew_diff)) 
  
df <- df[with(df, order(mass_diff)),]

df$flew_diff <- relevel(as.factor(df$flew_diff), ref = "0")

Null model:

null <- multinom(flew_diff ~ 1, data = df)

## # weights:  6 (2 variable)
## initial  value 305.414216 
## iter  10 value 174.928010
## iter  10 value 174.928009
## iter  10 value 174.928009
## final  value 174.928009 
## converged

Mass_diff model:

model <- multinom(flew_diff ~ mass_diff, data = df)

## # weights:  9 (4 variable)
## initial  value 305.414216 
## iter  10 value 166.730466
## final  value 166.154287 
## converged

s <- summary(model) # multinom doesn't compute pvalues so need to do it by hand.
z <- s$coefficients/s$standard.errors  #  calculate p-values using Wald tests (here z-tests)
wald <- z^2
p <- (1 - pnorm(abs(z), 0, 1)) * 2

table <- cbind(s$coefficients, c(s$edf, s$edf), s$standard.errors[,2], z[,2], wald[,2], p[,2])
colnames(table) <- c("(Intercept)"," Estimate","DF", "Std. Err.", "z", "wald", "P > |z|")
cat("\n")

cat("AIC: ", s$AIC, "\n")

## AIC:  340.3086

table

##    (Intercept)  Estimate DF Std. Err.          z        wald      P > |z|
## -1   -1.699907  52.84829  4  13.64735  3.8724203 14.99563911 0.0001077599
## 1    -3.066996  -4.57556  4  27.85556 -0.1642602  0.02698141 0.8695263283

anova(null, model, test="Chisq") # adding mass_diff improves fit

## Likelihood ratio tests of Multinomial Models
## 
## Response: flew_diff
##       Model Resid. df Resid. Dev   Test    Df LR stat.      Pr(Chi)
## 1         1       554   349.8560                                   
## 2 mass_diff       552   332.3086 1 vs 2     2 17.54744 0.0001547465

The ML prediction equations are:

\(\log(\hat{\pi_{-1}}/\hat{\pi_0}) = -1.70 + 52.85 x\) and

\(\log(\hat{\pi_1}/\hat{\pi_0}) = -3.07 - 4.58 x\)

The estimated log odds that the response is -1 (Flew in T1, not T2) rather than 1 (flew in T2, not T1) equals:

\(\log(\hat{\pi_{-1}}/\hat{\pi_1}) = (-1.70 - (-3.07)) + (52.85 - (- 4.58)) x = 1.37 + 57.43 x\) where x = mass_diff

So, the more positive the mass difference (positive (+) mass means gained mass from T1 to T2, and negative (-) mass means lost mass), the more likely to select flew_diff = -1, meaning that the bug flew in T1 but not T2. While the smaller/more negative the mass difference, the more likely to select flew_diff = 1, meaning that the bug flew in T2 but not T1.

For bugs of mass_diff x + 1/50 (0.02) g, the estimated odds that flew_diff = -1 rather than flew_diff = 1 equals

exp(52.85/50)

## [1] 2.877725

times the estimated odds at mass_diff x (g).

exp(coef(model)/50)

##    (Intercept) mass_diff
## -1   0.9665733 2.8776264
## 1    0.9405035 0.9125511

The relative risk ratio for a 1/50-unit increase in the variable mass_diff is 2.90 for flew_diff = -1 (T1 flew only) rather than flew_diff = 1 (T2 flew only).

Predicted probabilities

You can also use predicted probabilities to help you understand the model. The multicategory logit model has an alternative expression in terms of the response probabilities. This is

\(\pi_j = \frac{e^{\alpha_j + \beta_j x}}{\sum_he^{\alpha_h + \beta_h x}}, j=1,...,J\)

The denominator is the same for each probability, and the numerators for various j sum to the denominator where \(\sum_j \pi_j = 1\). The parameters (\(\alpha and \beta\)) equal zero for whichever category is the baseline in the logit expressions. So these would be the equations for the model above:

\(\hat{\pi_-1} = \frac{e^{-1.70 + 50.85 x}}{1 + e^{-1.70 + 50.85x} + e^{-3.07 - 4.58 x}}, j=1,...,J\)

\(\hat{\pi_1} = \frac{e^{-3.07 - 4.58 x}}{1 + e^{-1.70 + 50.85x} + e^{-3.07 - 4.58 x}}, j=1,...,J\)

\(\hat{\pi_0} = \frac{e^{0 + 0 x} = 1}{1 + e^{-1.70 + 50.85x} + e^{-3.07 - 4.58x}}, j=1,...,J\)

x = mass_diff which can give you the estimated probabilities.

You can calculate these predicted probabilities for each of our outcome levels using the fitted function. You can start by generating the predicted probabilities for the observations in our dataset and viewing the first few rows. Then you can plot them.

head(pp <- fitted(model))

##           0         -1          1
## 1 0.9309721 0.01577063 0.05325723
## 2 0.9268548 0.02155929 0.05158587
## 3 0.9259895 0.02270809 0.05130243
## 4 0.9219192 0.02793021 0.05015061
## 5 0.9194762 0.03096174 0.04956208
## 6 0.9167257 0.03431055 0.04896370

plot(df$mass_diff, pp[,1], ylim=c(0,1), col="red", type="l", ylab="Predicted Probability", xlab="Changes in Mass From T1 to T2 (g)") # no flight_diff
points(df$mass_diff, pp[,2], col="blue", type="l") # flew in T1 but not T2
points(df$mass_diff, pp[,3], col="green", type="l") # flew in T2 but not T1
text(0.02,0.8, labels="No Flight Diff", col="red")
text(0.0,0.35, labels="Flew in T1 but not T2", col="blue")
text(0.02,0.1, labels="Flew in T2 but not T1", col="green")
abline(v=0.035, lty=2)
abline(v=-0.027, lty=2)

Summary: Those who gain mass (at least more than 0.035 g) are more likely to fly in T1 but not T2 and less likely to have no change in flight response. Those who loose mass are more likely to fly in T2 but not T1 (but need to loose significant mass - more than 0.03 (g)).

Now let’s compare some models:

data <- data.frame(R = df$flew_diff, 
         A = df$mass_diff, 
         B = df$host_c, 
         C = df$sex_c)
source("src/compare_models.R")
model_comparisonsAIC("src/generic multinomial models- multinom 1RF + 3 FF.R")

##        [,1]      [,2]      [,3]      [,4]       [,5]       [,6]      
## AICs   333.1306  333.6353  334.6408  335.332    335.5504   335.8177  
## models 5         11        13        16         7          14        
## probs  0.2715071 0.2109556 0.1275994 0.09031437 0.08097221 0.07083936
##        [,7]      
## AICs   335.9544  
## models 9         
## probs  0.06615921
## 
## m5   multinom(formula = R ~ A + C, data = data, trace = FALSE)
## m11  multinom(formula = R ~ A * B + C, data = data, trace = FALSE)
## m13  multinom(formula = R ~ B * C + A, data = data, trace = FALSE)
## m16  multinom(formula = R ~ B * C + A * B, data = data, trace = FALSE)
## m7   multinom(formula = R ~ A + B + C, data = data, trace = FALSE)
## m14  multinom(formula = R ~ A * B + A * C, data = data, trace = FALSE)
## m9   multinom(formula = R ~ A * C, data = data, trace = FALSE)

anova(m5, m7, test="Chisq") # Adding B does not improve fit

## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1     A + C       550   321.1306                                
## 2 A + B + C       548   319.5504 1 vs 2     2 1.580235 0.4537914

anova(m5, m9, test="Chisq") # Adding A*C does not improve fit

## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##   Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1 A + C       550   321.1306                                
## 2 A * C       548   319.9544 1 vs 2     2 1.176152 0.5553949

delta_mass_model <- multinom(flew_diff ~ mass_diff + sex_c, data = df)

## # weights:  12 (6 variable)
## initial  value 305.414216 
## iter  10 value 167.318267
## iter  20 value 161.019802
## iter  30 value 160.852946
## iter  40 value 160.741183
## iter  50 value 160.658113
## iter  60 value 160.603172
## iter  70 value 160.587839
## iter  80 value 160.572910
## iter  90 value 160.565452
## final  value 160.565295 
## converged

s <- summary(delta_mass_model) # multinom doesn't compute pvalues so need to do it by hand.
z <- s$coefficients/s$standard.errors  #  calculate p-values using Wald tests (here z-tests)
wald <- z^2
p <- (1 - pnorm(abs(z), 0, 1)) * 2

table <- cbind(s$coefficients, c(s$edf, s$edf), s$standard.errors[,2:3], z[,2:3], wald[,2:3], p[,2:3])
colnames(table) <- c("(Intercept)","coeff mass_diff", "coeff sex_c","DF", "MD SE", "sex SE", "MD z", "sex z", 
                     "MD wald","sex wald","MD P > |z|", "sex P > |z|")
cat("\n", "AIC: ", s$AIC, "\n", "MD = mass_diff", "\n")

## 
##  AIC:  333.1306 
##  MD = mass_diff

table

##    (Intercept) coeff mass_diff coeff sex_c DF    MD SE     sex SE       MD z
## -1   -1.863449        61.25330  -0.2881294  6 15.99454  0.1991477  3.8296384
## 1    -8.900205       -55.86612  -6.3472957  6 67.78557 89.6299371 -0.8241596
##          sex z   MD wald    sex wald   MD P > |z| sex P > |z|
## -1 -1.44681204 14.666130 2.093265087 0.0001283317   0.1479496
## 1  -0.07081669  0.679239 0.005015004 0.4098489056   0.9435436

The general ML prediction equaiton is:

\(\log(\hat{\pi_{j}}/\hat{\pi_0}) =\alpha_{j} + \beta_{-1}^{\delta mass} x_1 + \beta_{j}^{sex} x_2\) where j = -1 or 1

The individual ML prediction equations are:

\(\log(\hat{\pi_{-1}}/\hat{\pi_0}) = -1.86 + 61.25 (\delta mass) - 0.29 (sex)\) and

\(\log(\hat{\pi_1}/\hat{\pi_0}) = -8.90 - 55.87 (\delta mass) - 6.35 (sex)\)

The estimated log odds that the response is -1 (Flew in T1, not T2) rather than 1 (flew in T2, not T1) equals:

\(\log(\hat{\pi_{-1}}/\hat{\pi_1}) = (-1.86 - (-8.90)) + (61.25 - (- 55.87)) \delta mass + (-0.29 - (-6.35)) sex\)

\(\log(\hat{\pi_{-1}}/\hat{\pi_1}) = 7.04 + 117.12 \delta mass + 6.06 sex\)

For bugs of mass_diff x + 1/50 (0.02) g and female, the estimated odds that flew_diff = -1 rather than flew_diff = 1 equals

exp((117.12+6.06)/50)

## [1] 11.74702

times (almost 12 times) the estimated odds at mass_diff x (g).

For bugs of mass_diff x + 1/50 (0.02) g and male, the estimated odds that flew_diff = -1 rather than flew_diff = 1 equals

exp((117.12-6.06)/50)

## [1] 9.218386

times the estimated odds at mass_diff x (g).

exp(coef(delta_mass_model)/50)

##    (Intercept) mass_diff     sex_c
## -1    0.963417 3.4043905 0.9942540
## 1     0.836939 0.3271546 0.8807813

head(pp <- fitted(delta_mass_model))

##           0          -1            1
## 1 0.9926638 0.007333222 2.928886e-06
## 2 0.9894420 0.010555888 2.087933e-06
## 3 0.9887828 0.011215207 1.973171e-06
## 4 0.9857139 0.014284503 1.573132e-06
## 5 0.9838824 0.016116154 1.404213e-06
## 6 0.9818204 0.018178338 1.253133e-06

plot(df$mass_diff, pp[,1], ylim=c(0,1), col="red", type="b", ylab="Predicted Probability", xlab="Changes in Mass From T1 to T2 (g)") # no flight_diff
points(df$mass_diff, pp[,2], col="blue", type="b") # flew in T1 but not T2
points(df$mass_diff, pp[,3], col="green", type="b") # flew in T2 but not T1
text(0.031,0.8, labels="No Flight Diff", col="red")
text(0.008,0.5, labels="Flew in T1 but not T2", col="blue")
text(0.052,0.1, labels="Flew in T2 but not T1", col="green")
# female line is always higher
text(0.0,0.96, labels="F", col="red")
text(0.0,0.7, labels="M", col="red")
text(0.0,0.25, labels="F", col="blue")
text(0.02,0.2, labels="M", col="blue")
text(-0.03,0.25, labels="F", col="green")
text(0.025,0.05, labels="M", col="green")
abline(v=0.03, lty=2, col="slategrey")
abline(v=0.035, lty=2, col="slategrey")
abline(v=0.04, lty=2, col="slategrey")
abline(v=-0.01, lty=2, col="slategrey")
abline(v=-0.045, lty=2, col="slategrey")

\(\log(\hat{\pi_{-1}}/\hat{\pi_1}) = 6.67 + 95.21 \delta mass + 5.45 sex\)

As females gain mass (need to gain more than males - 0.04 vs. 0.03 g) they are more likely to fly in T1 than T2. As females loose mass they are much more likely than males to fly in T2 than T1. But both are mor likely to have no mass flight response change.

More on multinomial plotting for when have more than 1 predictor variables: https://stats.idre.ucla.edu/r/dae/multinomial-logistic-regression/

Flew Diff ~ Egg Diff (females only)

df <- d %>%
  filter(!is.na(egg_diff), !is.na(flew_diff), sex_c == 1)
  
df$flew_diff <- relevel(as.factor(df$flew_diff), ref = "0")

Null model:

null <- multinom(flew_diff ~ 1, data = df)

## # weights:  2 (1 variable)
## initial  value 66.542129 
## final  value 46.327446 
## converged

Egg_diff model:

model <- multinom(flew_diff ~ egg_diff, data = df)

## # weights:  3 (2 variable)
## initial  value 66.542129 
## iter  10 value 40.203864
## iter  10 value 40.203864
## iter  10 value 40.203864
## final  value 40.203864 
## converged

s <- summary(model) # multinom doesn't compute pvalues so need to do it by hand.
z <- s$coefficients/s$standard.errors  #  calculate p-values using Wald tests (here z-tests)
wald <- z^2
p <- (1 - pnorm(abs(z), 0, 1)) * 2

table <- cbind(s$coefficients, c(s$edf, s$edf), s$standard.errors[2], z[2], wald[2], p[2])
colnames(table) <- c(" Estimate","DF", "Std. Err.", "z", "wald", "P > |z|")
cat("\n")

cat("AIC: ", s$AIC, "\n")

## AIC:  84.40773

table

##              Estimate DF Std. Err.        z     wald    P > |z|
## (Intercept) -2.219471  2 0.5486809 3.234434 10.46156 0.00121884
## egg_diff     1.774672  2 0.5486809 3.234434 10.46156 0.00121884

anova(null, model, test="Chisq") # adding egg_diff improves fit

## Likelihood ratio tests of Multinomial Models
## 
## Response: flew_diff
##      Model Resid. df Resid. Dev   Test    Df LR stat.      Pr(Chi)
## 1        1        -1   92.65489                                   
## 2 egg_diff        -2   80.40773 1 vs 2     1 12.24716 0.0004659658

The ML prediction equations are:

\(\log(\hat{\pi_{-1}}/\hat{\pi_0}) = -2.22 + 1.78 x\) and that’s it. Why? Because these comparisons have been reduced to a binomial situation. There are no flight responses = 1. None of the female bugs Fly in T2 but not T1. So we can do our regular binomial comparisons.

# stopped here
# Remove any missing masses
df <- df[!is.na(df$mass_diff),]

# Change the -1's to 1's, but still need to keep in mind 1 will then mean that all the bugs flew in T1 but not T2
df$flew_diff <- abs(as.integer(df$flew_diff)) - 1
df$flew_diff <- as.factor(df$flew_diff)

R = df$flew_diff
A = df$egg_diff
B = df$mass_diff
#C = df$sym_dist
C = df$host_c

data<-data.frame(R, A, B, C)

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 3-FF.R")

##        [,1]     [,2]      [,3]      [,4]       [,5]       [,6]       [,7]      
## AICs   73.28097 73.75454  74.48026  75.28029   75.42404   75.67354   75.73601  
## models 7        12        13        11         16         6          14        
## probs  0.210355 0.1660043 0.1154862 0.07741166 0.07204284 0.06359357 0.06163784
## 
## m7   glm(formula = R ~ A + B + C, family = binomial, data = data)
## m12  glm(formula = R ~ A * C + B, family = binomial, data = data)
## m13  glm(formula = R ~ B * C + A, family = binomial, data = data)
## m11  glm(formula = R ~ A * B + C, family = binomial, data = data)
## m16  glm(formula = R ~ A * C + B * C, family = binomial, data = data)
## m6   glm(formula = R ~ B + C, family = binomial, data = data)
## m14  glm(formula = R ~ A * B + A * C, family = binomial, data = data)

anova(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        89     65.281                     
## 2        88     63.755  1   1.5264   0.2166

anova(m6, m7, test="Chisq") # Addin A does 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        90     69.674                       
## 2        89     65.281  1   4.3926   0.0361 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

delta_egg_model <- glm(flew_diff ~ egg_diff + mass_diff + host_c, family = binomial, data = df)
tidy_regression(delta_egg_model, is_color = output_col)

## glm flew_diff ~ egg_diff + mass_diff + host_c binomial df 
## AIC:  73.28097 
## (Intercept)  coeff:  -2.1794225  Pr(>|t|):  2.381079e-06 *
## egg_diff     coeff:  1.3447603   Pr(>|t|):  0.03846687 *
## mass_diff    coeff:  40.2176192  Pr(>|t|):  0.02047179 *
## host_c       coeff:  0.7246695   Pr(>|t|):  0.02731218 *

strong effect of egg_diff where the more eggs a female bug laid the more likely a female bug flew in T1 but not T2
strong effect of mass where higher the mass gain the more likely the female bug flew in T1 but not T2. The smaller the mass_diff or if lost mass then the more likely the female bug did not change its flight response.
positive effect of host where if from GRT then more likely to have flown in T1 but not T2

gf_point(flew_diff ~ mass_diff, data=df, col=~egg_diff)

Summary

delta_mass_model

## Call:
## multinom(formula = flew_diff ~ mass_diff + sex_c, data = df)
## 
## Coefficients:
##    (Intercept) mass_diff      sex_c
## -1   -1.863449  61.25330 -0.2881294
## 1    -8.900205 -55.86612 -6.3472957
## 
## Residual Deviance: 321.1306 
## AIC: 333.1306

Look within the ‘Delta flight response’ tab to see the making of the graphs generated from this model. Here is the graph:

tidy_regression(delta_egg_model, is_color = output_col)

## glm flew_diff ~ egg_diff + mass_diff + host_c binomial df 
## AIC:  73.28097 
## (Intercept)  coeff:  -2.1794225  Pr(>|t|):  2.381079e-06 *
## egg_diff     coeff:  1.3447603   Pr(>|t|):  0.03846687 *
## mass_diff    coeff:  40.2176192  Pr(>|t|):  0.02047179 *
## host_c       coeff:  0.7246695   Pr(>|t|):  0.02731218 *

Trial 1

Testing glmer() with pop and chamber

Cleaning the Data

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

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

data_T1 <-data_tested[data_tested$trial_type=="T1",]
data_T1 <- center_data(data_T1)

Experimental 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.6064583

Biological 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.6957995

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

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 3-FF.R")

##        [,1]      [,2]      [,3]     [,4]     [,5]       [,6]      
## AICs   403.0853  405.0821  405.4319 405.5382 406.4848   406.594   
## models 8         11        2        15       4          6         
## probs  0.3595604 0.1324857 0.111231 0.105473 0.06570126 0.06221179
## 
## 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)

anova(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.0031848    0.955

anova(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.544    0.214

anova(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.0041746   0.9485

nomass_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 less likely to fly.

With Mass

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 4-FF.R")

##        [,1]       [,2]       [,3]    
## AICs   395.3856   395.7653   395.8828
## models 63         26         51      
## probs  0.08435731 0.06977224 0.06579 
## 
## 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)

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 3-FF.R")

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

anova(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 ' ' 1

anova(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.4972    0.114

anova(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.1242  0.02359 *
## ---
## 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.8828 
## (Intercept)      coeff:  0.4352287   Pr(>|t|):  0.04304544 *
## host_c           coeff:  -0.1640368  Pr(>|t|):  0.3508711
## mass_c           coeff:  -15.1464735 Pr(>|t|):  0.07961529 .
## sym_dist         coeff:  -0.1084774  Pr(>|t|):  0.5151625
## host_c:mass_c    coeff:  22.7792912  Pr(>|t|):  0.001430864 *
## mass_c:sym_dist  coeff:  -17.3605309 Pr(>|t|):  0.04408037 *

Morphology

data_T1 <-data_tested[data_tested$trial_type=="T1",]
data_T1 <- data_T1 %>%
  filter(!is.na(body))
data_T1 <- center_data(data_T1)

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 3-FF.R")

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

anova(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 ' ' 1

anova(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.451

morph_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.1002469

R = data_T1$flew_b
A = data_T1$thorax_c
B = data_T1$wing2body_c

data<-data.frame(R, A, B)

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 2-FF.R")

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

anova(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.1306

morph_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.8828 
## (Intercept)      coeff:  0.4352287   Pr(>|t|):  0.04304544 *
## host_c           coeff:  -0.1640368  Pr(>|t|):  0.3508711
## mass_c           coeff:  -15.1464735 Pr(>|t|):  0.07961529 .
## sym_dist         coeff:  -0.1084774  Pr(>|t|):  0.5151625
## host_c:mass_c    coeff:  22.7792912  Pr(>|t|):  0.001430864 *
## mass_c:sym_dist  coeff:  -17.3605309 Pr(>|t|):  0.04408037 *

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

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

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

data_T2 <-data_tested[data_tested$trial_type=="T2",]
data_T2 <- center_data(data_T2)

Experimental 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.5857749

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 3-FF.R")

##        [,1]      [,2]      [,3]     [,4]       [,5]       [,6]      
## AICs   361.1644  362.5404  363.0369 364.4845   364.5336   364.8156  
## models 2         4         6        7          8          10        
## probs  0.3641601 0.1830193 0.142783 0.06923625 0.06756002 0.05867243
## 
## 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)

anova(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.4296

anova(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.12746   0.7211

nomass_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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 4-FF.R")

##        [,1]       [,2]       [,3]      [,4]      
## AICs   358.8454   358.9575   359.2153  359.563   
## models 9          4          7         12        
## probs  0.08891803 0.08406952 0.0739026 0.06210944
## 
## 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)

anova(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.1869

anova(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.1986

anova(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.1461

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

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 3-FF.R")

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

anova(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 ' ' 1

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 2-FF.R")

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

anova(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.3158

morph_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

Eggs Modeling

Cleaning the Data

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

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

d <- data_tested %>%
   group_by(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, wing2body_s) %>%
   summarise_all(funs(list(na.omit(.))))

d$num_flew <- 0
d$num_notflew <- 0
d$average_mass <- 0

for(row in 1:length(d$flew_b)){
  
  n_flew <- sum(d$flew_b[[row]] == 1) # total number of times did fly among trails
  d$num_flew[[row]] <- n_flew 
  
  n_notflew <- sum(d$flew_b[[row]] == 0) # total number of times did not fly among trails
  d$num_notflew[[row]] <- n_notflew
  
  avg_mass <- mean(d$mass[[row]])
  d$average_mass[[row]] <- avg_mass

}

d <- select(d, -filename, -channel_letter, -set_number)

data_fem <- d[d$sex=="F",]
data_fem <- center_data(data_fem, is_not_binded = FALSE)
data_fem

## # A tibble: 120 x 63
## # 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 [120]
##    ID    sex   population site  host_plant latitude longitude total_eggs  beak
##    <fct> <fct> <fct>      <fct> <fct>         <dbl>     <dbl>      <dbl> <dbl>
##  1 5     F     Plantatio… Areg… C. corind…     25.0     -80.6        209  7.47
##  2 7     F     Plantatio… Foun… C. corind…     25.0     -80.6          8  8.19
##  3 8     F     Plantatio… Foun… C. corind…     25.0     -80.6         92  8.39
##  4 12    F     Key Largo  KLMRL C. corind…     25.1     -80.4        329  8.79
##  5 13    F     Key Largo  KLMRL C. corind…     25.1     -80.4        296  6.99
##  6 16    F     Key Largo  JP    C. corind…     25.1     -80.4          1  6.33
##  7 19    F     Key Largo  JP    C. corind…     25.1     -80.4         NA  7.54
##  8 20    F     Key Largo  JP    C. corind…     25.1     -80.4         55  8.73
##  9 25    F     North Key… DD f… C. corind…     25.3     -80.3         NA  8.03
## 10 26    F     North Key… Char… C. corind…     25.2     -80.3         33  7.15
## # … with 110 more rows, and 54 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>,
## #   total_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>, 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>, average_mass <dbl>, average_mass_c <dbl>,
## #   average_mass_s <dbl>, trial_type_og <list>, num_flew <dbl>,
## #   num_notflew <dbl>

Experimental, Biological, and Morphological Effects

## glm cbind(num_flew, num_notflew) ~ average_mass binomial data_fem 
## AIC:  245.2278 
## (Intercept)      coeff:  0.7621848   Pr(>|t|):  0.1847196
## average_mass     coeff:  -19.1280942 Pr(>|t|):  0.00927343 *

## glm cbind(num_flew, num_notflew) ~ total_eggs binomial data_fem 
## AIC:  180.1734 
## (Intercept)  coeff:  -0.000878   Pr(>|t|):  0.9971877
## total_eggs   coeff:  -0.0113041  Pr(>|t|):  3.560437e-05 *

## glm cbind(num_flew, num_notflew) ~ beak_c binomial data_fem 
## AIC:  244.4807 
## (Intercept)  coeff:  -0.7462472  Pr(>|t|):  5.668239e-07 *
## beak_c       coeff:  0.5097352   Pr(>|t|):  0.005387361 *

## glm cbind(num_flew, num_notflew) ~ thorax_c binomial data_fem 
## AIC:  250.7503 
## (Intercept)  coeff:  -0.7208467  Pr(>|t|):  7.657079e-07 *
## thorax_c     coeff:  0.7177388   Pr(>|t|):  0.1741137

## glm cbind(num_flew, num_notflew) ~ body_c binomial data_fem 
## AIC:  246.7111 
## (Intercept)  coeff:  -0.7397474  Pr(>|t|):  6.056371e-07 *
## body_c       coeff:  0.3567498   Pr(>|t|):  0.0182318 *

## glm cbind(num_flew, num_notflew) ~ wing_c binomial data_fem 
## AIC:  246.6118 
## (Intercept)  coeff:  -0.7379341  Pr(>|t|):  6.425547e-07 *
## wing_c       coeff:  0.4261769   Pr(>|t|):  0.01788375 *

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$average_mass 
X = data_fem$population # can't do trial type anymore

data<-data.frame(R1, R2, A, B, C, X)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 2-FF.R"))

##        [,1]      [,2]      [,3]      [,4]      
## AICs   248.6118  249.7548  250.4805  252.6136  
## models 2         4         3         5         
## probs  0.4643959 0.2622274 0.1824284 0.06279057
## 
## 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)

length(errors$warnings)

## [1] 0

anova(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)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m0  2 252.61 258.19 -124.31   248.61                           
## m2  3 248.61 256.97 -121.31   242.61 6.0019      1    0.01429 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m2, m3, test="Chisq") # Adding A does not improve fit

## Data: data
## Models:
## m2: cbind(R1, R2) ~ B + (1 | X)
## m3: cbind(R1, R2) ~ A + B + (1 | X)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2  3 248.61 256.97 -121.31   242.61                         
## m3  4 250.48 261.63 -121.24   242.48 0.1312      1     0.7171

nomass_fem <- glmer(cbind(num_flew, num_notflew) ~ wing_c + (1|population), data=data_fem, family=binomial)

## boundary (singular) fit: see ?isSingular

tidy_regression(nomass_fem, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ wing_c + (1 | population) data_fem binomial 
## AIC:  248.6118 248.6118 
## (Intercept)  coeff:  -0.7379341  Pr(>|t|):  6.426616e-07 *
## wing_c       coeff:  0.4261769   Pr(>|t|):  0.01788511 *

positive effect of wing length where the longer the wing the more times the bug will fly

With Mass

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 3-FF.R"))

##        [,1]      [,2]      [,3]      [,4]      [,5]       [,6]       [,7]      
## AICs   232.7316  233.6511  233.7866  233.975   234.7259   235.6114   235.6516  
## models 12        11        14        6         16         15         7         
## probs  0.2437271 0.1538936 0.1438118 0.1308843 0.08991714 0.05774999 0.05660236
##        [,8]      
## AICs   235.7758  
## models 17        
## probs  0.05319211
## 
## 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)
## m7   cbind(R1, R2) ~ A + B + C + (1 | X)
## m17  cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)

length(errors$warnings)

## [1] 0

anova(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)
##     Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)  
## m7   5 235.65 249.59 -112.83   225.65                          
## m12  6 232.73 249.46 -110.37   220.73  4.92      1    0.02655 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m12  6 232.73 249.46 -110.37   220.73                         
## m14  7 233.79 253.30 -109.89   219.79 0.9449      1      0.331

mass_fem <- glmer(cbind(num_flew, num_notflew) ~ host_c * average_mass + wing_c + (1|population), data=data_fem, family=binomial)

## boundary (singular) fit: see ?isSingular

tidy_regression(mass_fem, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ host_c * average_mass + wing_c + (1 | population) data_fem binomial 
## AIC:  232.7316 232.7316 
## (Intercept)          coeff:  0.961684    Pr(>|t|):  0.1989489
## host_c               coeff:  -1.589422   Pr(>|t|):  0.02368153 *
## average_mass         coeff:  -21.9064605 Pr(>|t|):  0.02794326 *
## wing_c               coeff:  0.7992536   Pr(>|t|):  0.0002548813 *
## host_c:average_mass  coeff:  20.3124168  Pr(>|t|):  0.02653981 *

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

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1-RF + 4-FF.R"))

##        [,1]       [,2]      
## AICs   177.1075   177.3402  
## models 9          20        
## probs  0.09166895 0.08160042
## 
## m9   cbind(R1, R2) ~ B + D + (1 | X)
## m20  cbind(R1, R2) ~ B * D + (1 | X)

length(errors$warnings) # 83 models failed

## [1] 83

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m9   4 177.11 187.65 -84.554   169.11                         
## m20  5 177.34 190.51 -83.670   167.34 1.7673      1     0.1837

egg_model <- glmer(cbind(num_flew, num_notflew) ~ wing_c + total_eggs_c + (1|population), data=data_fem, family=binomial)

## boundary (singular) fit: see ?isSingular

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:  177.1075 177.1075 
## (Intercept)      coeff:  -1.188912   Pr(>|t|):  5.122316e-09 *
## wing_c           coeff:  0.5576031   Pr(>|t|):  0.01181403 *
## total_eggs_c     coeff:  -0.0110275  Pr(>|t|):  4.517168e-05 *

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)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 3-FF.R"))

##        [,1]     [,2]      [,3]      [,4]       [,5]      
## AICs   245.4993 246.243   246.9547  247.4702   248.2304  
## models 16       15        4         17         5         
## probs  0.224862 0.1550342 0.1086168 0.08393821 0.05739459
## 
## 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)
## m17  cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)
## m5   cbind(R1, R2) ~ A + C + (1 | X)

length(errors$warnings)

## [1] 0

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m15  7 246.24 265.75 -116.12   232.24                             
## m16  7 245.50 265.01 -115.75   231.50 0.7437      0  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m13  6 248.52 265.25 -118.26   236.52                           
## m16  7 245.50 265.01 -115.75   231.50 5.0242      1      0.025 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

morph_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 ?isSingular

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:  245.4993 245.4993 
## (Intercept)      coeff:  -0.5878878  Pr(>|t|):  0.002748832 *
## thorax_c         coeff:  -4.0650614  Pr(>|t|):  0.0185575 *
## wing_c           coeff:  -0.1386652  Pr(>|t|):  0.8813351
## body_c           coeff:  1.5980987   Pr(>|t|):  0.09451309 .
## thorax_c:wing_c  coeff:  4.5899997   Pr(>|t|):  0.03769363 *
## wing_c:body_c    coeff:  -1.5360787  Pr(>|t|):  0.02326468 *

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)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 2-FF.R"))

##        [,1]      [,2]      [,3]      [,4]      [,5]     
## AICs   250.9253  251.4847  252.6136  252.6343  252.7503 
## models 2         4         5         3         1        
## probs  0.3319029 0.2509184 0.1426887 0.1412225 0.1332674
## 
## m2   cbind(R1, R2) ~ B + (1 | X)
## m4   cbind(R1, R2) ~ A * B + (1 | X)
## m0   cbind(R1, R2) ~ 1 + (1 | X)
## m3   cbind(R1, R2) ~ A + B + (1 | X)
## m1   cbind(R1, R2) ~ A + (1 | X)

length(errors$warnings)

## [1] 0

anova(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)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m0  2 252.61 258.19 -124.31   248.61                           
## m2  3 250.93 259.29 -122.46   244.93 3.6884      1    0.05479 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(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)
##    Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## m2  3 250.93 259.29 -122.46   244.93                        
## m3  4 252.63 263.78 -122.32   244.63 0.291      1     0.5896

morph_fem2 <- glmer(cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population), data=d, family=binomial)

## boundary (singular) fit: see ?isSingular

tidy_regression(morph_fem2, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population) d binomial 
## AIC:  250.9253 250.9253 
## (Intercept)  coeff:  -0.7209243  Pr(>|t|):  8.661145e-07 *
## wing2body_c  coeff:  16.2277047  Pr(>|t|):  0.06379128 .

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:  248.6118 248.6118 
## (Intercept)  coeff:  -0.7379341  Pr(>|t|):  6.426616e-07 *
## wing_c       coeff:  0.4261769   Pr(>|t|):  0.01788511 *

tidy_regression(mass_fem, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ host_c * average_mass + wing_c + (1 | population) data_fem binomial 
## AIC:  232.7316 232.7316 
## (Intercept)          coeff:  0.961684    Pr(>|t|):  0.1989489
## host_c               coeff:  -1.589422   Pr(>|t|):  0.02368153 *
## average_mass         coeff:  -21.9064605 Pr(>|t|):  0.02794326 *
## wing_c               coeff:  0.7992536   Pr(>|t|):  0.0002548813 *
## host_c:average_mass  coeff:  20.3124168  Pr(>|t|):  0.02653981 *

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:  177.1075 177.1075 
## (Intercept)      coeff:  -1.188912   Pr(>|t|):  5.122316e-09 *
## wing_c           coeff:  0.5576031   Pr(>|t|):  0.01181403 *
## total_eggs_c     coeff:  -0.0110275  Pr(>|t|):  4.517168e-05 *

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:  245.4993 245.4993 
## (Intercept)      coeff:  -0.5878878  Pr(>|t|):  0.002748832 *
## thorax_c         coeff:  -4.0650614  Pr(>|t|):  0.0185575 *
## wing_c           coeff:  -0.1386652  Pr(>|t|):  0.8813351
## body_c           coeff:  1.5980987   Pr(>|t|):  0.09451309 .
## thorax_c:wing_c  coeff:  4.5899997   Pr(>|t|):  0.03769363 *
## wing_c:body_c    coeff:  -1.5360787  Pr(>|t|):  0.02326468 *

tidy_regression(morph_fem2, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population) d binomial 
## AIC:  250.9253 250.9253 
## (Intercept)  coeff:  -0.7209243  Pr(>|t|):  8.661145e-07 *
## wing2body_c  coeff:  16.2277047  Pr(>|t|):  0.06379128 .

Females - Old Method

glmer() A=host, B=mass, C=wing, X=ID

Eggs Modeling

glmer() A=host, B=wing_c, C=mass, D=eggs_b, X=ID

glmer() A=thorax, B=body, C=wing X=trial_type

Cleaning the Data

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

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

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

Experimental 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.5687423

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 2-FF.R")

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

anova(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 ' ' 1

anova(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_tested[data_tested$sex=="F",]
data_fem <- data_fem %>%
  filter(!is.na(mass))
data_fem <- center_data(data_fem)

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 3-FF.R")

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

anova(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 ' ' 1

anova(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.7965

anova(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.6238

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glmer 1 RF + 3-FF.R")

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0513385 (tol = 0.001, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0731019 (tol = 0.001, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.085792 (tol = 0.001, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

##        [,1]      [,2]      [,3]       [,4]      
## AICs   203.132   204.311   206.663    206.7247  
## models 14        17        11         6         
## probs  0.4150825 0.2302117 0.07102064 0.06886306
## 
## 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)

# 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.5424898 Pr(>|t|):  0.0001373015 *
## host_c           coeff:  -3.4870152  Pr(>|t|):  0.1081201
## mass_c           coeff:  -358.5443417    Pr(>|t|):  0.0007737767 *
## wing_c           coeff:  8.6005095   Pr(>|t|):  0.002577142 *
## host_c:mass_c    coeff:  -194.6761103    Pr(>|t|):  0.01841898 *
## host_c:wing_c    coeff:  6.807007    Pr(>|t|):  0.01546067 *

# 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.350303e-07 *
## mass_c       coeff:  -37.2787841 Pr(>|t|):  5.206041e-06 *
## wing_c       coeff:  0.8582809   Pr(>|t|):  9.370685e-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)

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 4-FF.R")

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

egg_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.779839    Pr(>|t|):  0.01153589 *
## wing_c       coeff:  1.8027943   Pr(>|t|):  0.2569706
## mass_c       coeff:  -81.8893648 Pr(>|t|):  0.1419278
## eggs_b       coeff:  -15.8877372 Pr(>|t|):  1.173023e-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)

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF.R")

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

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi 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 ' ' 1

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi 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.1327

egg_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.3889394  Pr(>|t|):  0.7755781
## mass_c           coeff:  -423.9801397    Pr(>|t|):  0.2711374
## eggs_b           coeff:  -9.237128   Pr(>|t|):  0.0003012895 *
## wing_c           coeff:  2.559265    Pr(>|t|):  0.4402197
## mass_c:eggs_b    coeff:  361.4809036 Pr(>|t|):  0.3137307

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)

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 3-FF.R")

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

anova(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 ' ' 1

morph_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.4019887

Trying wing2body

R = data_fem$flew_b
A = data_fem$thorax_c
B = data_fem$wing2body_c

data<-data.frame(R, A, B)

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glm 2-FF.R")

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

anova(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 ' ' 1

anova(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.6548

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

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glmer 1 RF + 3-FF.R")

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0240326 (tol = 0.001, component 1)

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

anova(m0, m3, test="Chisq") # Adding C does not improve fit

## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m3: R ~ C + (1 | X)
##    Df    AIC    BIC  logLik deviance  Chisq Chi 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.3206

anova(m0, m2, test="Chisq") # Adding B does not improve fit

## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    Df    AIC    BIC  logLik deviance  Chisq Chi 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.3247

anova(m0, m1, test="Chisq") # Adding A does not improve fit

## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    Df    AIC    BIC  logLik deviance  Chisq Chi 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)

source("src/compare_models.R")
model_comparisonsAIC("src/generic models-binomial glmer 1-RF + 2-FF.R")

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

anova(m0, m2, test="Chisq") # Adding A marginally imprvoes fit

## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
##    Df    AIC    BIC  logLik deviance  Chisq Chi 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 ' ' 1

anova(m0, m1, test="Chisq")

## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m1: R ~ A + (1 | X)
##    Df    AIC    BIC  logLik deviance  Chisq Chi 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.5424898 Pr(>|t|):  0.0001373015 *
## host_c           coeff:  -3.4870152  Pr(>|t|):  0.1081201
## mass_c           coeff:  -358.5443417    Pr(>|t|):  0.0007737767 *
## wing_c           coeff:  8.6005095   Pr(>|t|):  0.002577142 *
## host_c:mass_c    coeff:  -194.6761103    Pr(>|t|):  0.01841898 *
## host_c:wing_c    coeff:  6.807007    Pr(>|t|):  0.01546067 *

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.350303e-07 *
## mass_c       coeff:  -37.2787841 Pr(>|t|):  5.206041e-06 *
## wing_c       coeff:  0.8582809   Pr(>|t|):  9.370685e-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.779839    Pr(>|t|):  0.01153589 *
## wing_c       coeff:  1.8027943   Pr(>|t|):  0.2569706
## mass_c       coeff:  -81.8893648 Pr(>|t|):  0.1419278
## eggs_b       coeff:  -15.8877372 Pr(>|t|):  1.173023e-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.3889394  Pr(>|t|):  0.7755781
## mass_c           coeff:  -423.9801397    Pr(>|t|):  0.2711374
## eggs_b           coeff:  -9.237128   Pr(>|t|):  0.0003012895 *
## wing_c           coeff:  2.559265    Pr(>|t|):  0.4402197
## mass_c:eggs_b    coeff:  361.4809036 Pr(>|t|):  0.3137307

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.4384809

tidy_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.4019887

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 *

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

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

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

d <- data_tested %>%
   group_by(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, wing2body_s) %>%
   summarise_all(funs(list(na.omit(.))))

d$num_flew <- 0
d$num_notflew <- 0
d$average_mass <- 0

for(row in 1:length(d$flew_b)){
  
  n_flew <- sum(d$flew_b[[row]] == 1) # total number of times did fly among trails
  d$num_flew[[row]] <- n_flew 
  
  n_notflew <- sum(d$flew_b[[row]] == 0) # total number of times did not fly among trails
  d$num_notflew[[row]] <- n_notflew
  
  avg_mass <- mean(d$mass[[row]])
  d$average_mass[[row]] <- avg_mass

}

d <- select(d, -filename, -channel_letter, -set_number)

data_male <- d[d$sex=="M",]
data_male <- center_data(data_male, is_not_binded = FALSE) 
data_male

## # A tibble: 213 x 63
## # 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 [213]
##    ID    sex   population site  host_plant latitude longitude total_eggs  beak
##    <fct> <fct> <fct>      <fct> <fct>         <dbl>     <dbl>      <dbl> <dbl>
##  1 1     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  5.5 
##  2 2     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  5.64
##  3 3     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  5.75
##  4 4     M     Plantatio… Areg… C. corind…     25.0     -80.6         NA  6.21
##  5 6     M     Plantatio… Foun… C. corind…     25.0     -80.6         NA  5.76
##  6 9     M     Plantatio… Foun… C. corind…     25.0     -80.6         NA  6.1 
##  7 10    M     Plantatio… Foun… 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 203 more rows, and 54 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>,
## #   total_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>, 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>, average_mass <dbl>, average_mass_c <dbl>,
## #   average_mass_s <dbl>, trial_type_og <list>, num_flew <dbl>,
## #   num_notflew <dbl>

Experimental, Biological, and Morphological Effects

## glm cbind(num_flew, num_notflew) ~ average_mass binomial data_male 
## AIC:  438.4795 
## (Intercept)      coeff:  0.6245134   Pr(>|t|):  0.3254026
## average_mass     coeff:  3.7451961   Pr(>|t|):  0.8169126

## glm cbind(num_flew, num_notflew) ~ beak_c binomial data_male 
## AIC:  434.3964 
## (Intercept)  coeff:  0.7697254   Pr(>|t|):  1.253049e-12 *
## beak_c       coeff:  0.559369    Pr(>|t|):  0.04417133 *

## glm cbind(num_flew, num_notflew) ~ thorax_c binomial data_male 
## AIC:  438.1161 
## (Intercept)  coeff:  0.769756    Pr(>|t|):  9.445937e-13 *
## thorax_c     coeff:  0.3470788   Pr(>|t|):  0.5183983

## glm cbind(num_flew, num_notflew) ~ body_c binomial data_male 
## AIC:  433.6353 
## (Intercept)  coeff:  0.773278    Pr(>|t|):  1.075558e-12 *
## body_c       coeff:  0.3439638   Pr(>|t|):  0.02711692 *

## glm cbind(num_flew, num_notflew) ~ wing_c binomial data_male 
## AIC:  431.0918 
## (Intercept)  coeff:  0.7774098   Pr(>|t|):  1.024409e-12 *
## wing_c       coeff:  0.495168    Pr(>|t|):  0.006665111 *

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$average_mass 
X = data_male$population # can't do trial type anymore

data<-data.frame(R1, R2, A, B, C, X)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 2-FF.R"))

##        [,1]      [,2]      [,3]     
## AICs   426.3587  427.9881  429.0797 
## models 2         3         4        
## probs  0.5634561 0.2494909 0.1445465
## 
## m2   cbind(R1, R2) ~ B + (1 | X)
## m3   cbind(R1, R2) ~ A + B + (1 | X)
## m4   cbind(R1, R2) ~ A * B + (1 | X)

length(errors$warnings)

## [1] 0

anova(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)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2  3 426.36 436.44 -210.18   420.36                         
## m3  4 427.99 441.43 -209.99   419.99 0.3707      1     0.5426

nomass_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:  426.3587 426.3587 
## (Intercept)  coeff:  0.6575878   Pr(>|t|):  0.001200224 *
## wing_c       coeff:  0.5535136   Pr(>|t|):  0.00490563 *

positive effect of wing

With Mass

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 3-FF.R"))

##        [,1]      [,2]      [,3]      [,4]       [,5]       [,6]      
## AICs   422.1946  422.6649  423.2442  423.8514   423.8748   424.0915  
## models 6         10        7         13         14         17        
## probs  0.2034319 0.1607984 0.1203663 0.08884698 0.08781608 0.07879646
##        [,7]       [,8]       [,9]      
## AICs   424.62     424.6987   424.7506  
## models 11         12         16        
## probs  0.06049864 0.05816318 0.05667366
## 
## m6   cbind(R1, R2) ~ B + C + (1 | X)
## m10  cbind(R1, R2) ~ B * C + (1 | X)
## m7   cbind(R1, R2) ~ A + B + C + (1 | X)
## m13  cbind(R1, R2) ~ B * C + A + (1 | X)
## m14  cbind(R1, R2) ~ A * B + A * C + (1 | X)
## m17  cbind(R1, R2) ~ A * B + A * C + B * C + (1 | X)
## m11  cbind(R1, R2) ~ A * B + C + (1 | X)
## m12  cbind(R1, R2) ~ A * C + B + (1 | X)
## m16  cbind(R1, R2) ~ A * C + B * C + (1 | X)

length(errors$warnings)

## [1] 0

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m6   4 422.19 435.64 -207.10   414.19                         
## m10  5 422.66 439.47 -206.33   412.66 1.5296      1     0.2162

anova(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)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m6  4 422.19 435.64 -207.10   414.19                         
## m7  5 423.24 440.05 -206.62   413.24 0.9504      1     0.3296

mass_male <- glmer(cbind(num_flew, num_notflew) ~ wing_c + average_mass + (1|population), data=data_male, family=binomial)
tidy_regression(mass_male, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ wing_c + average_mass + (1 | population) data_male binomial 
## AIC:  422.1946 422.1946 
## (Intercept)      coeff:  2.8003774   Pr(>|t|):  0.001610106 *
## wing_c           coeff:  0.9875339   Pr(>|t|):  0.0002309966 *
## average_mass     coeff:  -56.2532673 Pr(>|t|):  0.0126619 *

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

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 3-FF.R"))

##        [,1]      [,2]      [,3]     [,4]     [,5]      [,6]      
## AICs   416.449   416.901   417.9681 418.0949 418.2036  419.539   
## models 13        10        16       9        15        12        
## probs  0.2463597 0.1965326 0.115268 0.10819  0.1024653 0.05255385
## 
## m13  cbind(R1, R2) ~ B * C + A + (1 | X)
## m10  cbind(R1, R2) ~ B * C + (1 | X)
## m16  cbind(R1, R2) ~ A * C + B * C + (1 | X)
## m9   cbind(R1, R2) ~ A * C + (1 | X)
## m15  cbind(R1, R2) ~ A * B + B * C + (1 | X)
## m12  cbind(R1, R2) ~ A * C + B + (1 | X)

length(errors$warnings)

## [1] 0

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m13  6 416.45 436.62 -202.22   404.45                         
## m16  7 417.97 441.50 -201.98   403.97 0.4809      1      0.488

anova(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)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m10  5 416.90 433.71 -203.45   406.90                         
## m13  6 416.45 436.62 -202.22   404.45 2.4519      1     0.1174

anova(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)
##     Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)   
## m6   4 424.38 437.83 -208.19   416.38                           
## m10  5 416.90 433.71 -203.45   406.90 9.481      1   0.002076 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

morph_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:  416.901 416.901 
## (Intercept)      coeff:  0.9329442   Pr(>|t|):  6.186581e-05 *
## body_c           coeff:  -1.1785281  Pr(>|t|):  0.05570531 .
## wing_c           coeff:  1.6028889   Pr(>|t|):  0.02741834 *
## body_c:wing_c    coeff:  -0.676228   Pr(>|t|):  0.00342541 *

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)

source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-binomial glmer 2R ~ 1 RF + 2-FF.R"))

##        [,1]      [,2]      [,3]     
## AICs   421.839   423.0562  423.7015 
## models 2         4         3        
## probs  0.5140051 0.2796796 0.2025521
## 
## m2   cbind(R1, R2) ~ B + (1 | X)
## m4   cbind(R1, R2) ~ A * B + (1 | X)
## m3   cbind(R1, R2) ~ A + B + (1 | X)

length(errors$warnings)

## [1] 0

anova(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)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2  3 421.84 431.92 -207.92   415.84                         
## m3  4 423.70 437.15 -207.85   415.70 0.1375      1     0.7108

morph_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:  421.839 421.839 
## (Intercept)  coeff:  0.6512969   Pr(>|t|):  0.004022727 *
## wing2body_c  coeff:  24.3255246  Pr(>|t|):  0.0005769884 *

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

data<-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-04

plot(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:  426.3587 426.3587 
## (Intercept)  coeff:  0.6575878   Pr(>|t|):  0.001200224 *
## wing_c       coeff:  0.5535136   Pr(>|t|):  0.00490563 *

tidy_regression(mass_male, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ wing_c + average_mass + (1 | population) data_male binomial 
## AIC:  422.1946 422.1946 
## (Intercept)      coeff:  2.8003774   Pr(>|t|):  0.001610106 *
## wing_c           coeff:  0.9875339   Pr(>|t|):  0.0002309966 *
## average_mass     coeff:  -56.2532673 Pr(>|t|):  0.0126619 *

tidy_regression(morph_male, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ body_c * wing_c + (1 | population) d binomial 
## AIC:  416.901 416.901 
## (Intercept)      coeff:  0.9329442   Pr(>|t|):  6.186581e-05 *
## body_c           coeff:  -1.1785281  Pr(>|t|):  0.05570531 .
## wing_c           coeff:  1.6028889   Pr(>|t|):  0.02741834 *
## body_c:wing_c    coeff:  -0.676228   Pr(>|t|):  0.00342541 *

tidy_regression(morph_male2, is_color=output_col)

## glmer cbind(num_flew, num_notflew) ~ wing2body_c + (1 | population) d binomial 
## AIC:  421.839 421.839 
## (Intercept)  coeff:  0.6512969   Pr(>|t|):  0.004022727 *
## wing2body_c  coeff:  24.3255246  Pr(>|t|):  0.0005769884 *

Males - Old Method

Without Mass Modeling + Days From Start

With Mass Modeling + Days From Start