Modeling Flight Speed
Anastasia Bernat
5/4/2020
Winter 2020 Flight Trials: Speed Flight Modeling
Flight Trials Winter 2020 Dataset was conducted from 2/17/2020 - 3/10/2020. Soapberry bugs were flight tested twice for multiple hours in the flight mill and observed from 8 AM to (5-8 PM) each day. Used multivariate (glm) and mixed effect modeling (glmer) to analyze the flight results.
Problems Noticing
qqplots of all models are not normal. Might need to retransform the speed_diff data either by standardizing it or performing another transformation.
errors “refitting model(s) with ML (instead of REML)” and “boundary singular…” showing up a lot
All Flyers
Cleaning the Data
source_path = "~/Desktop/git_repositories/SBB-dispersal/avbernat_working_on/Rsrc/"
script_names = c("center_flight_data.R", # Re-centers data
"clean_flight_data.R", # Loads and cleans data
"unique_flight_data.R",
"get_warnings.R",
"compare_models.R",
"regression_output.R", # Cleans regression outputs; prints in color or B&W
"AICprobabilities.R",
"plotting-lm.R",
"plotting-lm2.R")
for (script in script_names) {
path = paste0(source_path, script)
source(path)
}##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionoutput_col = FALSE # Change to TRUE if working in Base R or RStudio; FALSE if generating an HTML
data <- read_flight_data("data/all_flight_data-Winter2020.csv")
data_all <- data[[1]]
data_tested <- data[[2]]
#d <- create_delta_data(data_tested)### Remove everyone who didn't fly
data_flew_all <- data_tested[data_tested$flew_b == 1, ]
### Check for low speeds
low_speeds <- data_flew_all %>%
filter(average_speed <0.05)
### Check for high speeds
high_speeds <- data_flew_all %>%
filter(average_speed >0.65)
low_speeds$flight_type# have 15 bugs with average_speed = 0 but were marked as bursters (this could be just something very short (second burst) - not enough to grant a calculation) - I decided to remove them. But one bug was continuous and had 0 distance and 0 speeds - that was bug 196 T2 set011-3-03-2020-A3_196.txt## [1] B B B B B B B B B B C B B B B
## Levels: B BC C CBhigh_speeds$flight_type # 3 bugs - also bursters. Could also be short explosive bursts but not true to the biology of these bugs (more like us blowing on them).## factor(0)
## Levels: B BC C CB### Remove those have 0 m/s speeds
data_flew <- data_flew_all %>%
filter(average_speed > 0.05)
data_flew <- center_data(data_flew)Testing Models and Covariates
test_model<-lmer(average_speed~host_c*sex_c*sym_dist + (1|ID), data=data_flew) #it converges - Run one set of models with ID as a random factor.
getME(test_model, "lower")## [1] 0## lmer average_speed ~ chamber + (1 | ID) data_flew
## AIC: -621.057 -641.057
## (Intercept) coeff: 0.346498 tval: 25.65336
## chamberA-1 coeff: 0.0141921 tval: 0.6776612
## chamberA-3 coeff: -0.0387153 tval: -2.119609
## chamberA-4 coeff: -0.0520144 tval: -2.978639
## chamberB-1 coeff: -0.009754 tval: -0.4836775
## chamberB-2 coeff: -0.0640839 tval: -3.57309
## chamberB-3 coeff: -0.0830768 tval: -4.211742
## chamberB-4 coeff: -0.0947142 tval: -4.947998## lmer average_speed ~ days_from_start + (1 | chamber) data_flew
## AIC: -642.687 -650.687
## (Intercept) coeff: 0.3063359 tval: 19.7795
## days_from_start coeff: -0.0001106 tval: -0.1613477## boundary (singular) fit: see ?isSingular## lmer average_speed ~ min_from_IncStart + (1 | trial_type) data_flew
## AIC: -611.892 -619.892
## (Intercept) coeff: 0.3121393 tval: 38.26806
## min_from_IncStart coeff: -5.28e-05 tval: -1.509154Without Mass
data<-data.frame(R=data_flew$average_speed,
A=data_flew$host_c,
B=data_flew$sex_c,
C=data_flew$sym_dist,
X=data_flew$ID)
model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
model_comparisonsAIC(model_script)## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -647.2402 -645.3635 -645.2482 -643.6194 -643.469 -643.2526
## models 2 4 6 7 8 10
## probs 0.3990449 0.1561342 0.1473938 0.06527899 0.06055064 0.05434196
##
## m2 R ~ B + (1 | X)
## m4 R ~ A + B + (1 | X)
## m6 R ~ B + C + (1 | X)
## m7 R ~ A + B + C + (1 | X)
## m8 R ~ A * B + (1 | X)
## m10 R ~ B * C + (1 | X)anova(m0, m2, test="Chisq") # Adding B improves fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ (1 | X)
## m2: R ~ B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 -639.12 -627.75 322.56 -645.12
## m2 4 -647.24 -632.07 327.62 -655.24 10.116 1 0.00147 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2, m4, test="Chisq") # Adding Adoes not improve fit## Data: data
## Models:
## m2: R ~ B + (1 | X)
## m4: R ~ A + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 4 -647.24 -632.07 327.62 -655.24
## m4 5 -645.36 -626.40 327.68 -655.36 0.1233 1 0.7255anova(m2, m6, test="Chisq") # Adding C does not improve fit## Data: data
## Models:
## m2: R ~ B + (1 | X)
## m6: R ~ B + C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 4 -647.24 -632.07 327.62 -655.24
## m6 5 -645.25 -626.28 327.62 -655.25 0.0081 1 0.9284speed_model <- lmer(average_speed ~ sex_c + (1 | ID), data=data_flew, REML=FALSE)
tidy_regression(speed_model, is_color=output_col)## lmer average_speed ~ sex_c + (1 | ID) data_flew FALSE
## AIC: -647.2402 -647.2402
## (Intercept) coeff: 0.3138262 tval: 51.2686
## sex_c coeff: 0.01962 tval: 3.205254- marginal effect of sex
s.test <- paste("pval: ", shapiro.test(residuals(speed_model))$p.value)
qqnorm(resid(speed_model))
qqline(resid(speed_model))
text(-2, 0.1, s.test)
With Mass
data_flew <- data_flew %>%
filter(!is.na(mass))
data_flew <- center_data(data_flew)
data<-data.frame(R=data_flew$average_speed,
A=data_flew$host_c,
B=data_flew$sex_c,
C=data_flew$sym_dist,
D=data_flew$mass_c,
X=data_flew$chamber)
model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 4-FF REMLF.R")
model_comparisonsAIC(model_script)## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -680.856 -679.928 -679.0715 -678.2061 -678.1845 -678.1173
## models 20 30 31 50 43 38
## probs 0.2065858 0.1298893 0.08464489 0.05491386 0.05432315 0.05252862
##
## m20 R ~ B * D + (1 | X)
## m30 R ~ B * D + A + (1 | X)
## m31 R ~ B * D + C + (1 | X)
## m50 R ~ A * D + B * D + (1 | X)
## m43 R ~ A * B + B * D + (1 | X)
## m38 R ~ B * D + A + C + (1 | X)anova(m0, m4, test="Chisq") # Adding D does improve fit## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m4: R ~ D + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 -661.40 -650.03 333.70 -667.40
## m4 4 -672.32 -657.16 340.16 -680.32 12.915 1 0.000326 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m0, m2, test="Chisq") # Adding B does improve fit## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m2: R ~ B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 -661.4 -650.03 333.70 -667.4
## m2 4 -674.5 -659.34 341.25 -682.5 15.096 1 0.0001022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1speed_model <- lmer(average_speed ~ sex_c + (1 | chamber), data=data_flew, REML=FALSE)
tidy_regression(speed_model, is_color=output_col)## lmer average_speed ~ sex_c + (1 | chamber) data_flew FALSE
## AIC: -674.4998 -674.4998
## (Intercept) coeff: 0.318611 tval: 23.69559
## sex_c coeff: 0.0226141 tval: 3.937881- positive effect of sex, so that if F then faster
s.test <- paste("pval: ", shapiro.test(residuals(speed_model))$p.value)
qqnorm(resid(speed_model))
qqline(resid(speed_model))
text(-2, 0.1, s.test)
Morphology
data_flew <- data_flew %>%
filter(!is.na(mass)) %>%
filter(!is.na(body))
data_flew <- center_data(data_flew)
data<-data.frame(R=data_flew$average_speed,
A=data_flew$thorax_c,
B=data_flew$body_c,
C=data_flew$wing_c,
X=data_flew$ID)
model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
model_comparisonsAIC(model_script)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs -658.6189 -657.5618 -657.4346 -657.0179 -656.7503 -656.7312 -656.3556
## models 3 1 5 14 6 7 2
## probs 0.1911358 0.1126632 0.1057212 0.08584012 0.07509004 0.07437402 0.0616421
##
## m3 R ~ C + (1 | X)
## m1 R ~ A + (1 | X)
## m5 R ~ A + C + (1 | X)
## m14 R ~ A * B + A * C + (1 | X)
## m6 R ~ B + C + (1 | X)
## m7 R ~ A + B + C + (1 | X)
## m2 R ~ B + (1 | X)anova(m0, m1, test="Chisq") # Adding A does improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ (1 | X)
## m1: R ~ A + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 -636.99 -625.62 321.49 -642.99
## m1 4 -657.56 -642.40 332.78 -665.56 22.573 1 2.023e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m0, m2, test="Chisq") # Adding B does improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ (1 | X)
## m2: R ~ B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 -636.99 -625.62 321.49 -642.99
## m2 4 -656.36 -641.20 332.18 -664.36 21.367 1 3.792e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m0, m3, test="Chisq") # Adding C does improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ (1 | X)
## m3: R ~ C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 -636.99 -625.62 321.49 -642.99
## m3 4 -658.62 -643.46 333.31 -666.62 23.631 1 1.167e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m1, m5, test="Chisq") # Addihng C does not improve fit## Data: data
## Models:
## m1: R ~ A + (1 | X)
## m5: R ~ A + C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 4 -657.56 -642.40 332.78 -665.56
## m5 5 -657.43 -638.48 333.72 -667.43 1.8728 1 0.1712anova(m1, m4, test="Chisq") # Adding B does not improve fit## Data: data
## Models:
## m1: R ~ A + (1 | X)
## m4: R ~ A + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 4 -657.56 -642.40 332.78 -665.56
## m4 5 -655.84 -636.89 332.92 -665.84 0.2757 1 0.5995speed_morph_model <- lmer(average_speed ~ thorax_c + (1|ID), data=data_flew, REML=FALSE)## boundary (singular) fit: see ?isSingulartidy_regression(speed_morph_model, is_color=output_col)## lmer average_speed ~ thorax_c + (1 | ID) data_flew FALSE
## AIC: -657.5618 -657.5618
## (Intercept) coeff: 0.3026606 tval: 62.5805
## thorax_c coeff: 0.0864877 tval: 4.841675- positive effect of thorax, where the longer the thorax, the faster the speed
s.test <- paste("pval: ", shapiro.test(residuals(speed_morph_model))$p.value)
qqnorm(resid(speed_morph_model))
qqline(resid(speed_morph_model))
text(-2, 0.1, s.test)
All Flyers Plots
gf_histogram(~average_speed,data=data_flew_all, col=~flight_type) # before filtering outliers
gf_histogram(~average_speed,data=data_flew, col=~flight_type) # after filtering outliers
Continuous Flyers
Cleaning the data
### Remove everyone who didn't fly (then remove distances = 0, if that didn't work fully)
data_flew <- data_tested[data_tested$flew_b == 1, ]
data_flew <- data_flew[data_flew$average_speed > 0, ]
data_flew <- center_data(data_flew)
### Break up by flight type
dC<-data_flew[data_flew$flight_type=="C" | data_flew$flight_type=="BC" | data_flew$flight_type=="CB",] # this includes BC or CB flyers
dC <- center_data(dC)Testing Models and Covariates
test_model<-lmer(average_speed~host_c*sex_c*sym_dist + (1|ID), data=dC) #it converges - Run one set of models with ID as a random factor.
getME(test_model, "lower")## [1] 0## lmer average_speed ~ chamber + (1 | ID) dC
## AIC: -341.6616 -361.6616
## (Intercept) coeff: 0.3425476 tval: 22.17328
## chamberA-1 coeff: 0.0294355 tval: 1.145375
## chamberA-3 coeff: -0.0333228 tval: -1.557068
## chamberA-4 coeff: -0.0540571 tval: -2.451812
## chamberB-1 coeff: -0.0051492 tval: -0.2315663
## chamberB-2 coeff: -0.0810492 tval: -3.671022
## chamberB-3 coeff: -0.0925548 tval: -3.712048
## chamberB-4 coeff: -0.1073299 tval: -4.7833## lmer average_speed ~ days_from_start + (1 | ID) dC
## AIC: -336.1318 -344.1318
## (Intercept) coeff: 0.2958053 tval: 25.24525
## days_from_start coeff: 0.0002718 tval: 0.3020569## lmer average_speed ~ min_from_IncStart + (1 | ID) dC
## AIC: -330.7818 -338.7818
## (Intercept) coeff: 0.3045469 tval: 28.48669
## min_from_IncStart coeff: -3.63e-05 tval: -0.7317197Without Mass
data<-data.frame(R=dC$average_speed,
A=dC$host_c,
B=dC$sex_c,
C=dC$sym_dist,
X=dC$ID)
model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
model_comparisonsAIC(model_script)## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -359.6489 -358.499 -357.9787 -356.9865 -356.7949 -356.6058
## models 2 4 6 1 3 7
## probs 0.2780126 0.1564441 0.1206103 0.07343791 0.06673105 0.06070952
## [,7]
## AICs -356.5292
## models 8
## probs 0.05842998
##
## m2 R ~ B + (1 | X)
## m4 R ~ A + B + (1 | X)
## m6 R ~ B + C + (1 | X)
## m1 R ~ A + (1 | X)
## m3 R ~ C + (1 | X)
## m7 R ~ A + B + C + (1 | X)
## m8 R ~ A * B + (1 | X)anova(m2, m4, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m2: R ~ B + (1 | X)
## m4: R ~ A + B + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 4 -359.65 -346.79 183.82 -367.65
## m4 5 -358.50 -342.42 184.25 -368.50 0.8501 1 0.3565anova(m2, m6, test="Chisq") # Adding B does not improve fit## Data: data
## Models:
## m2: R ~ B + (1 | X)
## m6: R ~ B + C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 4 -359.65 -346.79 183.82 -367.65
## m6 5 -357.98 -341.90 183.99 -367.98 0.3298 1 0.5658continuous_model <- lmer(average_speed ~ sex_c + (1|ID), data=dC, REML=FALSE)
tidy_regression(continuous_model, is_color=output_col)## lmer average_speed ~ sex_c + (1 | ID) dC FALSE
## AIC: -359.6489 -359.6489
## (Intercept) coeff: 0.3075654 tval: 36.37051
## sex_c coeff: 0.015495 tval: 1.832329continuous_model <- lmer(average_speed ~ sex_c + (1|chamber) + (1|ID), data=dC, REML=FALSE)
tidy_regression(continuous_model, is_color=output_col) # Adding chamber improves fit## lmer average_speed ~ sex_c + (1 | chamber) + (1 | ID) dC FALSE
## AIC: -384.2756 -384.2756
## (Intercept) coeff: 0.3104449 tval: 18.11553
## sex_c coeff: 0.0190961 tval: 2.438286- positive effect of sex_c, where if a female continuous-flyer bug then more likely to fly faster on average
s.test <- paste("pval: ", shapiro.test(residuals(continuous_model))$p.value)
qqnorm(resid(continuous_model))
qqline(resid(continuous_model))
text(-1, 0.1, s.test)
With Mass
data<-data.frame(R=dC$average_speed,
A=dC$host_c,
B=dC$sex_c,
C=dC$sym_dist,
D=dC$mass_c,
X=dC$ID)
model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 4-FF REMLF.R")
errors = withWarnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3]
## AICs -363.1026 -361.3765 -361.177
## models 20 30 31
## probs 0.2148525 0.09063804 0.08203591
##
## m20 R ~ B * D + (1 | X)
## m30 R ~ B * D + A + (1 | X)
## m31 R ~ B * D + C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 1anova(m20, m30, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m20: R ~ B * D + (1 | X)
## m30: R ~ B * D + A + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m20 6 -363.10 -343.85 187.55 -375.10
## m30 7 -361.38 -338.91 187.69 -375.38 0.2738 1 0.6008anova(m20, m31, test="Chisq") # Adding C does not improve fit## Data: data
## Models:
## m20: R ~ B * D + (1 | X)
## m31: R ~ B * D + C + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m20 6 -363.10 -343.85 187.55 -375.10
## m31 7 -361.18 -338.71 187.59 -375.18 0.0744 1 0.785continuous_model <- lmer(average_speed ~ sex_c*mass_c + (1|ID), data=dC)
tidy_regression(continuous_model, is_color=output_col)## lmer average_speed ~ sex_c * mass_c + (1 | ID) dC
## AIC: -349.7477 -361.7477
## (Intercept) coeff: 0.326688 tval: 24.59279
## sex_c coeff: 0.0165075 tval: 1.242668
## mass_c coeff: 1.3044507 tval: 1.705021
## sex_c:mass_c coeff: -2.1438024 tval: -2.80212continuous_model <- lmer(average_speed ~ sex_c*mass_c + (1|chamber) + (1|ID), data=dC)
tidy_regression(continuous_model, is_color=output_col) # Adding chamber improves model## lmer average_speed ~ sex_c * mass_c + (1 | chamber) + (1 | ID) dC
## AIC: -376.4037 -390.4037
## (Intercept) coeff: 0.3308877 tval: 16.20645
## sex_c coeff: 0.0208851 tval: 1.736473
## mass_c coeff: 1.2220835 tval: 1.766986
## sex_c:mass_c coeff: -2.220121 tval: -3.145823- no effect of sex_c
- no effect of mass_c
- negative effect of sex_c*mass_c, where if female and heavier, then slower avg speed (make sense because probably due to egg-laying)
s.test <- paste("pval: ", shapiro.test(residuals(continuous_model))$p.value)
qqnorm(resid(continuous_model))
qqline(resid(continuous_model))
text(-1, 0.1, s.test)
Morphology
dC<-data_flew[data_flew$flight_type=="C",]
dC <- dC %>%
filter(!is.na(body))
dC <- center_data(dC)
data<-data.frame(R=dC$average_speed,
A=dC$thorax_c,
B=dC$body_c,
C=dC$wing_c,
X=dC$ID)
model_script = paste0(source_path,"generic models-gaussian glmer 1-RF + 3-FF REMLF.R")
errors = withWarnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -286.5804 -286.5032 -286.1151 -286.0869 -286.0455 -285.7361
## models 10 9 2 14 1 13
## probs 0.1066223 0.1025859 0.08449057 0.08330788 0.08160233 0.06990504
## [,7] [,8] [,9] [,10]
## AICs -285.7337 -285.6573 -285.6435 -285.1613
## models 17 3 8 15
## probs 0.06982295 0.06720644 0.06674425 0.05244422
##
## m10 R ~ B * C + (1 | X)
## m9 R ~ A * C + (1 | X)
## m2 R ~ B + (1 | X)
## m14 R ~ A * B + A * C + (1 | X)
## m1 R ~ A + (1 | X)
## m13 R ~ B * C + A + (1 | X)
## m17 R ~ A * B + A * C + B * C + (1 | X)
## m3 R ~ C + (1 | X)
## m8 R ~ A * B + (1 | X)
## m15 R ~ A * B + B * C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m10, m13, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m10: R ~ B * C + (1 | X)
## m13: R ~ B * C + A + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m10 6 -286.58 -268.8 149.29 -298.58
## m13 7 -285.74 -265.0 149.87 -299.74 1.1557 1 0.2824continuous_morph_model <- lmer(average_speed ~ body_c*wing_c + (1|ID), data=dC)
tidy_regression(continuous_morph_model, is_color=output_col)## lmer average_speed ~ body_c * wing_c + (1 | ID) dC
## AIC: -258.8789 -270.8789
## (Intercept) coeff: 0.3136288 tval: 34.68586
## body_c coeff: 0.0490269 tval: 1.349415
## wing_c coeff: -0.0196597 tval: -0.4202444
## body_c:wing_c coeff: -0.0191163 tval: -2.10057continuous_morph_model <- lmer(average_speed ~ body_c*wing_c + (1|ID) + (1|chamber), data=dC)
tidy_regression(continuous_morph_model, is_color=output_col) # improves model## lmer average_speed ~ body_c * wing_c + (1 | ID) + (1 | chamber) dC
## AIC: -268.9601 -282.9601
## (Intercept) coeff: 0.3117589 tval: 19.42156
## body_c coeff: 0.0500163 tval: 1.44706
## wing_c coeff: -0.0223555 tval: -0.5024798
## body_c:wing_c coeff: -0.0173923 tval: -2.031037- positive effect of body length, where the longer the body the faster the avg speed
- negative effect of wing length, where the longer the wing the slower the avg speed
- neg eggect of body* wing, where the longer the wing *body the slower the avg speed
s.test <- paste("pval: ", shapiro.test(residuals(continuous_morph_model))$p.value)
qqnorm(resid(continuous_morph_model))
qqline(resid(continuous_morph_model))
text(-1, 0.1, s.test)
Continuous Flyers Plots
gf_histogram(~average_speed,data=dC)
Bursters
Cleaning the data
### Remove everyone who didn't fly (then remove distances = 0, if that didn't work fully)
data_flew <- data_tested[data_tested$flew_b == 1, ]
data_flew <- data_flew %>%
filter(average_speed > 0) %>%
filter(average_speed < 0.65)
data_flew <- center_data(data_flew)
### Break up by flight type
dB<-data_flew[data_flew$flight_type=="B",] # this exclusdes BC or CB flyers
dB <- center_data(dB)Testing Models and Covariates
test_model<-lmer(average_speed~host_c*sex_c*sym_dist + (1|ID), data=dB) #it converges - Run one set of models with ID as a random factor.
getME(test_model, "lower")## [1] 0## lmer average_speed ~ chamber + (1 | ID) dB
## AIC: -225.0329 -245.0329
## (Intercept) coeff: 0.3549856 tval: 13.75944
## chamberA-1 coeff: -0.002658 tval: -0.0746513
## chamberA-3 coeff: -0.0411172 tval: -1.239861
## chamberA-4 coeff: -0.058236 tval: -1.942705
## chamberB-1 coeff: -0.0353766 tval: -0.8706123
## chamberB-2 coeff: -0.0497021 tval: -1.61028
## chamberB-3 coeff: -0.0840925 tval: -2.498602
## chamberB-4 coeff: -0.0757869 tval: -2.182173## lmer average_speed ~ days_from_start + (1 | ID) dB
## AIC: -252.2945 -260.2945
## (Intercept) coeff: 0.3005704 tval: 21.3898
## days_from_start coeff: 0.000591 tval: 0.5517944## boundary (singular) fit: see ?isSingular## lmer average_speed ~ min_from_IncStart + (1 | trial_type) dB
## AIC: -248.2631 -256.2631
## (Intercept) coeff: 0.3272916 tval: 24.90064
## min_from_IncStart coeff: -9.34e-05 tval: -1.832417Without Mass
data<-data.frame(R=dB$average_speed,
A=dB$host_c,
B=dB$sex_c,
C=dB$sym_dist,
X=dB$ID, Y=dB$chamber)
model_script = paste0(source_path,"generic models-gaussian glmer 2-RF + 3-FF.R")
errors = withWarnings(model_comparisonsAIC(model_script))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -265.832 -265.4378 -264.3919 -263.7995 -263.2714 -262.3944
## models 54 18 36 20 2 38
## probs 0.3131483 0.2571244 0.152417 0.1133411 0.08703753 0.05613952
##
## m0 R ~ (1 | X)
## m18 R ~ (1 | Y)
## m36 R ~ (1 | X) + (1 | Y)
## m20 R ~ B + (1 | Y)
## m2 R ~ B + (1 | X)
## m38 R ~ B + (1 | X) + (1 | Y)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m18, m0, test="Chisq") # Replacing Y with X improves the fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m18: R ~ (1 | Y)
## m0: R ~ (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m18 3 -272.95 -264.08 139.48 -278.95
## m0 3 -273.65 -264.78 139.82 -279.65 0.6978 0anova(m0, m36, test="Chisq") # Adding Y Marginally improves the fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ (1 | X)
## m36: R ~ (1 | X) + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 -273.65 -264.78 139.82 -279.65
## m36 4 -271.84 -260.02 139.92 -279.84 0.1937 1 0.6599burster_model <- lmer(average_speed ~ (1|ID), data=dB) # null model the best fit model
tidy_regression(burster_model, is_color=output_col)## lmer average_speed ~ (1 | ID) dB
## AIC: -265.832 -271.832
## 1 coeff: 0.3069246 tval: 38.21524s.test <- paste("pval: ", shapiro.test(residuals(burster_model))$p.value)
qqnorm(resid(burster_model))
qqline(resid(burster_model))
text(-1, 0.1, s.test)
With Mass
data<-data.frame(R=dB$average_speed,
A=dB$host_c,
B=dB$sex_c,
C=dB$sym_dist,
D=dB$mass_c,
X=dB$trial_type, Y=dB$chamber)
model_script = paste0(source_path,"generic models-gaussian glmer 2-RF + 4-FF.R")
errors = withWarnings(model_comparisonsAIC(model_script))## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## [,1] [,2] [,3] [,4] [,5]
## AICs -268.4246 -267.842 -266.4246 -265.4378 -264.8565
## models 117 4 230 113 339
## probs 0.3179288 0.2375868 0.1169595 0.07140806 0.05339739
##
## m117 R ~ D + (1 | Y)
## m4 R ~ D + (1 | X)
## m230 R ~ D + (1 | X) + (1 | Y)
## m113 R ~ 1 + (1 | Y)
## m0 R ~ 1 + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m113, m117, test="Chisq") # Adding D does not improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m113: R ~ 1 + (1 | Y)
## m117: R ~ D + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m113 3 -272.95 -264.08 139.48 -278.95
## m117 4 -275.98 -264.16 141.99 -283.98 5.0296 1 0.02492 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m113, m226, test="Chisq") # Adidng X does not improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m113: R ~ 1 + (1 | Y)
## m226: R ~ 1 + (1 | X) + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m113 3 -272.95 -264.08 139.48 -278.95
## m226 4 -270.95 -259.13 139.48 -278.95 0 1 1anova(m0, m4, test="Chisq") # Adding D does not improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ 1 + (1 | X)
## m4: R ~ D + (1 | X)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 -272.77 -263.90 139.39 -278.77
## m4 4 -275.79 -263.96 141.89 -283.79 5.0165 1 0.02511 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m226, m230, test="Chisq") # Adding D does not improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m226: R ~ 1 + (1 | X) + (1 | Y)
## m230: R ~ D + (1 | X) + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m226 4 -270.95 -259.13 139.48 -278.95
## m230 5 -273.98 -259.20 141.99 -283.98 5.0296 1 0.02492 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1burster_model <- lmer(average_speed ~ (1|chamber), data=dB) # null model the best fit model
tidy_regression(burster_model, is_color=output_col)## lmer average_speed ~ (1 | chamber) dB
## AIC: -265.4378 -271.4378
## 1 coeff: 0.3093318 tval: 31.46028s.test <- paste("pval: ", shapiro.test(residuals(burster_model))$p.value)
qqnorm(resid(burster_model))
qqline(resid(burster_model))
text(-1, 0.1, s.test)
Morphology
dB<-data_flew[data_flew$flight_type=="B",]
dB <- dB %>%
filter(!is.na(body))
dB <- center_data(dB)
data<-data.frame(R=dB$average_speed,
A=dB$thorax_c,
B=dB$body_c,
C=dB$wing_c,
X=dB$trial_type, Y=dB$chamber)
model_script = paste0(source_path,"generic models-gaussian glmer 2-RF + 3-FF.R")
errors = withWarnings(model_comparisonsAIC(model_script))## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs -273.1931 -272.7015 -272.0438 -271.368 -271.1931 -271.0022 -270.296
## models 21 3 19 1 39 20 2
## probs 0.2318323 0.1813148 0.1304977 0.09307913 0.08528632 0.0775225 0.05446066
##
## m21 R ~ C + (1 | Y)
## m3 R ~ C + (1 | X)
## m19 R ~ A + (1 | Y)
## m1 R ~ A + (1 | X)
## m39 R ~ C + (1 | X) + (1 | Y)
## m20 R ~ B + (1 | Y)
## m2 R ~ B + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 0anova(m18, m19, test="Chisq") # Adding A does marginally improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m18: R ~ (1 | Y)
## m19: R ~ A + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m18 3 -272.95 -264.08 139.48 -278.95
## m19 4 -285.09 -273.27 146.54 -293.09 14.138 1 0.0001699 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m0, m36, test="Chisq") # Adding Y does not impove fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ (1 | X)
## m36: R ~ (1 | X) + (1 | Y)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m0 3 -272.77 -263.90 139.39 -278.77
## m36 4 -270.95 -259.13 139.48 -278.95 0.1791 1 0.6722burster_morph_model <- lmer(average_speed ~ thorax_c + (1|chamber), data=dB)
tidy_regression(burster_morph_model, is_color=output_col) # thorax not significant - null model best fit ## lmer average_speed ~ thorax_c + (1 | chamber) dB
## AIC: -272.0438 -280.0438
## (Intercept) coeff: 0.3096333 tval: 32.57902
## thorax_c coeff: 0.0995045 tval: 3.833736burster_morph_model <- lmer(average_speed ~ (1|chamber), data=dB)
tidy_regression(burster_morph_model, is_color=output_col) ## lmer average_speed ~ (1 | chamber) dB
## AIC: -265.4378 -271.4378
## 1 coeff: 0.3093318 tval: 31.46028s.test <- paste("pval: ", shapiro.test(residuals(burster_morph_model))$p.value)
qqnorm(resid(burster_morph_model))
qqline(resid(burster_morph_model))
text(-1, 0.1, s.test)
Burster Plots
gf_histogram(~average_speed,data=dB)
