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.
All Flyers
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':
##
## datesource("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]]
### Remove everyone who didn't fly (then remove distances = 0, if that didn't work fully)
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 7 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 C 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).## [1] B B B
## Levels: B BC C CB### Remove outliers
data_flew <- data_flew_all %>%
filter(average_speed > 0.05) %>%
filter(average_speed < 0.65)
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: -562.5465 -582.5465
## (Intercept) coeff: 0.3506719 tval: 23.79403
## chamberA-1 coeff: 0.0356284 tval: 1.557293
## chamberA-3 coeff: -0.0437178 tval: -2.187905
## chamberA-4 coeff: -0.0230032 tval: -1.197706
## chamberB-1 coeff: -0.0137146 tval: -0.6225656
## chamberB-2 coeff: -0.0632558 tval: -3.210209
## chamberB-3 coeff: -0.0402426 tval: -1.896561
## chamberB-4 coeff: -0.0978115 tval: -4.678375## lmer average_speed ~ days_from_start + (1 | chamber) data_flew
## AIC: -584.9133 -592.9133
## (Intercept) coeff: 0.3144849 tval: 19.82496
## days_from_start coeff: 0.0004536 tval: 0.604685## lmer average_speed ~ min_from_IncStart + (1 | trial_type) data_flew
## AIC: -557.1897 -565.1897
## (Intercept) coeff: 0.3252516 tval: 32.2628
## min_from_IncStart coeff: -4.34e-05 tval: -1.138938Without 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)
source("src/compare_models.R")
model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R")## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -587.268 -585.3789 -585.2915 -583.4716 -583.3959 -583.362
## models 2 4 6 8 7 10
## probs 0.3598186 0.1399207 0.1339365 0.05391452 0.05191336 0.05104095
##
## m2 R ~ B + (1 | X)
## m4 R ~ A + B + (1 | X)
## m6 R ~ B + C + (1 | X)
## m8 R ~ A * B + (1 | X)
## m7 R ~ A + B + C + (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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -584.96 -573.59 295.48 -590.96
## m2 4 -587.27 -572.11 297.63 -595.27 4.3055 1 0.03799 *
## ---
## 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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 4 -587.27 -572.11 297.63 -595.27
## m4 5 -585.38 -566.43 297.69 -595.38 0.111 1 0.7391anova(m2, m6, test="Chisq") # Adding C does not improve fit## Data: data
## Models:
## m2: R ~ B + (1 | X)
## m6: R ~ B + C + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 4 -587.27 -572.11 297.63 -595.27
## m6 5 -585.29 -566.34 297.65 -595.29 0.0235 1 0.8781speed_model <- lmer(average_speed ~ sex_c + (1 | ID), data=data_flew, REML=FALSE)## boundary (singular) fit: see ?isSingulartidy_regression(speed_model, is_color=output_col)## lmer average_speed ~ sex_c + (1 | ID) data_flew FALSE
## AIC: -587.268 -587.268
## (Intercept) coeff: 0.3248411 tval: 48.96227
## sex_c coeff: 0.0138119 tval: 2.081823- 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)
source("src/compare_models.R")
model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF REMLF.R")## [,1] [,2]
## AICs -608.5036 -607.7868
## models 2 4
## probs 0.08430075 0.05890849
##
## m2 R ~ B + (1 | X)
## m4 R ~ D + (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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -603.21 -591.85 304.61 -609.21
## m4 4 -607.79 -592.64 307.89 -615.79 6.5728 1 0.01035 *
## ---
## 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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -603.21 -591.85 304.61 -609.21
## m2 4 -608.50 -593.36 308.25 -616.50 7.2896 1 0.006935 **
## ---
## 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: -608.5036 -608.5036
## (Intercept) coeff: 0.3296762 tval: 23.99555
## sex_c coeff: 0.0172659 tval: 2.720721- 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)
source("src/compare_models.R")
model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R")## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -592.3049 -592.2102 -592.2089 -590.8658 -590.7851 -590.5642
## models 1 3 14 5 2 17
## probs 0.1451363 0.1384228 0.138331 0.07067493 0.06788193 0.06078307
## [,7] [,8] [,9]
## AICs -590.4638 -590.3663 -590.3064
## models 6 7 4
## probs 0.05780608 0.05505543 0.05343095
##
## m1 R ~ A + (1 | X)
## m3 R ~ C + (1 | X)
## m14 R ~ A * B + A * C + (1 | X)
## m5 R ~ A + C + (1 | X)
## m2 R ~ B + (1 | X)
## m17 R ~ A * B + A * C + B * C + (1 | X)
## m6 R ~ B + C + (1 | X)
## m7 R ~ A + B + C + (1 | X)
## m4 R ~ A + 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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -583.06 -571.70 294.53 -589.06
## m1 4 -592.30 -577.16 300.15 -600.30 11.243 1 0.0007991 ***
## ---
## 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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -583.06 -571.70 294.53 -589.06
## m2 4 -590.79 -575.64 299.39 -598.79 9.7236 1 0.001819 **
## ---
## 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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -583.06 -571.70 294.53 -589.06
## m3 4 -592.21 -577.06 300.11 -600.21 11.149 1 0.0008409 ***
## ---
## 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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 4 -592.30 -577.16 300.15 -600.30
## m5 5 -590.87 -571.93 300.43 -600.87 0.5608 1 0.4539anova(m1, m4, test="Chisq") # Adding B does not improve fit## Data: data
## Models:
## m1: R ~ A + (1 | X)
## m4: R ~ A + B + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 4 -592.30 -577.16 300.15 -600.30
## m4 5 -590.31 -571.37 300.15 -600.31 0.0014 1 0.9698speed_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: -592.3049 -592.3049
## (Intercept) coeff: 0.3170594 tval: 59.40831
## thorax_c coeff: 0.0667505 tval: 3.382239- 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
source("src/plotting-lm.R")
source("src/plotting-lm2.R")
Delta Speed
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]]
### Remove everyone who didn't fly (then remove distances = 0, if that didn't work fully)
data_flew_all <- data_tested[data_tested$flew_b == 1, ]
### Remove outliers
data_flew <- data_flew_all %>%
filter(average_speed > 0.05) %>%
filter(average_speed < 0.65)
data_flew <- center_data(data_flew)d <- data_flew %>%
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
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: 206 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 [206]
## 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 7 F Plantatio… Foun… C. corind… 25.0 -80.6 8 8.19
## 7 9 M Plantatio… Foun… C. corind… 25.0 -80.6 NA 6.1
## 8 10 M Plantatio… Foun… C. corind… 25.0 -80.6 NA 6.09
## 9 11 M Key Largo KLMRL C. corind… 25.1 -80.4 NA 5.88
## 10 12 F Key Largo KLMRL C. corind… 25.1 -80.4 329 8.79
## # … with 196 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># 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## # A tibble: 121 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 [121]
## 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 3 M Plantatio… Areg… C. corind… 25.0 -80.6 NA 5.75
## 3 6 M Plantatio… Foun… C. corind… 25.0 -80.6 NA 5.76
## 4 7 F Plantatio… Foun… C. corind… 25.0 -80.6 8 8.19
## 5 10 M Plantatio… Foun… C. corind… 25.0 -80.6 NA 6.09
## 6 11 M Key Largo KLMRL C. corind… 25.1 -80.4 NA 5.88
## 7 15 M Key Largo JP C. corind… 25.1 -80.4 NA 6.53
## 8 17 M Key Largo JP C. corind… 25.1 -80.4 NA 6.18
## 9 19 F Key Largo JP C. corind… 25.1 -80.4 NA 7.54
## 10 21 M North Key… DD f… C. corind… 25.3 -80.3 NA 6.03
## # … with 111 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>Speed Diff ~ Mass Diff
## Test
tidy_regression(lm(speed_diff ~ mass_diff, data=d), is_color = output_col) # no effect of mass_diff AFTER you filter out the outliers. If don't filter out the outliers then this relationship is significant and that data is residuals are not normal## lm speed_diff ~ mass_diff d
## AIC: -129.9329
## (Intercept) coeff: 0.0093927 Pr(>|t|): 0.4662632
## mass_diff coeff: -0.9686529 Pr(>|t|): 0.4027201hist(d$speed_diff) # much more normally distributed data and later will see the residuals pass the Shapiro-Wilk test
No mass_diff
R = d$speed_diff
A = d$host_c
B = d$sex_c
C = d$sym_dist
D = d$mass_diff
X = d$population
data<-data.frame(R, A, B, C, D, X)
source("src/compare_models.R")
model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R")## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -135.171 -134.456 -133.3357 -133.2134 -132.6768 -132.5989
## models 7 13 12 11 15 16
## probs 0.2366627 0.1655267 0.09453491 0.08892941 0.0680017 0.06540492
## [,7]
## AICs -132.3844
## models 5
## probs 0.05875271
##
## m7 R ~ A + B + C + (1 | X)
## m13 R ~ B * C + A + (1 | X)
## m12 R ~ A * C + B + (1 | X)
## m11 R ~ A * B + C + (1 | X)
## m15 R ~ A * B + B * C + (1 | X)
## m16 R ~ A * C + B * C + (1 | X)
## m5 R ~ A + C + (1 | X)anova(m7, m13, test="Chisq") # Adding B*C does not improve fit## Data: data
## Models:
## m7: R ~ A + B + C + (1 | X)
## m13: R ~ B * C + A + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m7 6 -135.17 -118.40 73.585 -147.17
## m13 7 -134.46 -114.89 74.228 -148.46 1.285 1 0.257speed_model <- lmer(speed_diff ~ host_c + sex_c + sym_dist + (1 | population), data=d, REML=FALSE)## boundary (singular) fit: see ?isSingulartidy_regression(speed_model, is_color=output_col)## lmer speed_diff ~ host_c + sex_c + sym_dist + (1 | population) d FALSE
## AIC: -135.171 -135.171
## (Intercept) coeff: 0.0585705 tval: 2.143171
## host_c coeff: 0.0509201 tval: 2.312272
## sex_c coeff: -0.0344581 tval: -2.209646
## sym_dist coeff: -0.060632 tval: -3.147974- (marginal) negative effect of sex where if F then more likely to fly slower in T2 than T1
- (marginal) negative effect of host where if from GRT then more likely to fly slower in T2 than T1
- negative effect of sym_dist where if farther from Homestead then more likely to fly slower in T2 than T1
s.test <- paste("pval: ", shapiro.test(residuals(speed_model))$p.value)
qqnorm(resid(speed_model))
qqline(resid(speed_model))
text(-1, 0.1, s.test)
Mass_diff
d <- d %>%
filter(!is.na(mass_diff))
R = d$speed_diff
A = d$host_c
B = d$sex_c
C = d$sym_dist
D = d$mass_diff
X = d$population
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 REMLF.R")## [,1] [,2]
## AICs -135.171 -134.456
## models 11 28
## probs 0.1107856 0.07748569
##
## m11 R ~ A + B + C + (1 | X)
## m28 R ~ B * C + A + (1 | X)anova(m11, m28, test="Chisq") # Adding B*C does not improve fit## Data: data
## Models:
## m11: R ~ A + B + C + (1 | X)
## m28: R ~ B * C + A + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m11 6 -135.17 -118.40 73.585 -147.17
## m28 7 -134.46 -114.89 74.228 -148.46 1.285 1 0.257speed_model <- lmer(speed_diff ~ host_c + sex_c + sym_dist + (1 | population), data=d, REML=FALSE)## boundary (singular) fit: see ?isSingulartidy_regression(speed_model, is_color=output_col)## lmer speed_diff ~ host_c + sex_c + sym_dist + (1 | population) d FALSE
## AIC: -135.171 -135.171
## (Intercept) coeff: 0.0585705 tval: 2.143171
## host_c coeff: 0.0509201 tval: 2.312272
## sex_c coeff: -0.0344581 tval: -2.209646
## sym_dist coeff: -0.060632 tval: -3.147974Same top model as before. All of these variables have weak effects.
So what if it was only females, since they had the greatest mass changes?
Delta Speed - Females Only
d_fem <- d[d$sex == "F",] # very low sample size No mass_diff
R = d_fem$speed_diff
A = d_fem$host_c
B = d_fem$sym_dist
C = d_fem$mass_diff
X = d_fem$population
data<-data.frame(R, A, B, C, X)
source("src/compare_models.R")
model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 2-FF REMLF.R")## [,1] [,2] [,3] [,4] [,5]
## AICs -18.93651 -18.19993 -17.16987 -16.80295 -16.20123
## models 5 1 2 4 3
## probs 0.3698021 0.2558721 0.1528792 0.1272548 0.09419171
##
## m0 R ~ (1 | X)
## m1 R ~ A + (1 | X)
## m2 R ~ B + (1 | X)
## m4 R ~ A * B + (1 | X)
## m3 R ~ A + B + (1 | X)anova(m0, m1, test="Chisq") # Adding A does not improve fit## Data: data
## Models:
## m0: R ~ (1 | X)
## m1: R ~ A + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -18.936 -15.530 12.468 -24.936
## m1 4 -18.200 -13.658 13.100 -26.200 1.2634 1 0.261anova(m0, m2, test="Chisq") # Adding B does not improve fit## Data: data
## Models:
## m0: R ~ (1 | X)
## m2: R ~ B + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -18.936 -15.530 12.468 -24.936
## m2 4 -17.170 -12.628 12.585 -25.170 0.2334 1 0.629speed_model <- lmer(speed_diff ~ 1 + (1 | population), data=d_fem, REML=FALSE)## boundary (singular) fit: see ?isSingulartidy_regression(speed_model, is_color=output_col)## lmer speed_diff ~ 1 + (1 | population) d_fem FALSE
## AIC: -18.93651 -18.93651
## 1 coeff: -0.0327732 tval: -1.116993mass_diff
source("src/compare_models.R")
model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R")## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -18.19993 -17.62841 -17.16987 -16.80295 -16.54578 -16.20123
## models 1 3 2 8 5 4
## probs 0.1780721 0.1338113 0.106395 0.08856192 0.07787564 0.06555192
## [,7] [,8]
## AICs -16.03881 -16.0132
## models 10 6
## probs 0.06043883 0.05966986
##
## m1 R ~ A + (1 | X)
## m3 R ~ C + (1 | X)
## m2 R ~ B + (1 | X)
## m8 R ~ A * B + (1 | X)
## m5 R ~ A + C + (1 | X)
## m4 R ~ A + B + (1 | X)
## m10 R ~ B * C + (1 | X)
## m6 R ~ B + C + (1 | X)anova(m0, m1, test="Chisq") # strange. get error "refitting model(s) with ML (instead of REML)" even though these have been refit to REML=FALSE...## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ (1 | X)
## m1: R ~ A + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -18.936 -15.530 12.468 -24.936
## m1 4 -18.200 -13.658 13.100 -26.200 1.2634 1 0.261anova(m0, m3, test="Chisq") # Adding C does not improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ (1 | X)
## m3: R ~ C + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -18.936 -15.530 12.468 -24.936
## m3 4 -17.628 -13.086 12.814 -25.628 0.6919 1 0.4055anova(m0, m2, test="Chisq") # Adding B does not improve fit## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## m0: R ~ (1 | X)
## m2: R ~ B + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -18.936 -15.530 12.468 -24.936
## m2 4 -17.170 -12.628 12.585 -25.170 0.2334 1 0.629speed_model <- lmer(speed_diff ~ 1 + (1 | population), data=d_fem, REML=FALSE)## boundary (singular) fit: see ?isSingulartidy_regression(speed_model, is_color=output_col)## lmer speed_diff ~ 1 + (1 | population) d_fem FALSE
## AIC: -18.93651 -18.93651
## 1 coeff: -0.0327732 tval: -1.116993s.test <- paste("pval: ", shapiro.test(residuals(speed_model))$p.value)
qqnorm(resid(speed_model))
qqline(resid(speed_model))
text(-1, 0.1, s.test)
Continuous Flyers
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]]
### 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.6909 -361.6909
## (Intercept) coeff: 0.3418703 tval: 22.2647
## chamberA-1 coeff: 0.0242972 tval: 0.9524348
## chamberA-3 coeff: -0.0327884 tval: -1.544873
## chamberA-4 coeff: -0.0460623 tval: -2.084857
## chamberB-1 coeff: -0.0041099 tval: -0.1862605
## chamberB-2 coeff: -0.0791824 tval: -3.615932
## chamberB-3 coeff: -0.0933522 tval: -3.764606
## chamberB-4 coeff: -0.10493 tval: -4.70953## lmer average_speed ~ days_from_start + (1 | ID) dC
## AIC: -338.2927 -346.2927
## (Intercept) coeff: 0.2964567 tval: 25.57216
## days_from_start coeff: 0.0002679 tval: 0.3024303## lmer average_speed ~ min_from_IncStart + (1 | ID) dC
## AIC: -332.9126 -340.9126
## (Intercept) coeff: 0.3049015 tval: 28.64662
## min_from_IncStart coeff: -3.45e-05 tval: -0.7039702Without Mass
data<-data.frame(R=dC$average_speed,
A=dC$host_c,
B=dC$sex_c,
C=dC$sym_dist,
X=dC$ID)
source("src/compare_models.R")
model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R")## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## AICs -361.7016 -360.4932 -360.0133 -359.1192 -358.9601 -358.5876 -358.5115
## models 2 4 6 1 3 7 8
## probs 0.2797239 0.1528701 0.1202616 0.07690534 0.0710261 0.0589576 0.05675413
##
## 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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 4 -361.70 -348.86 184.85 -369.70
## m4 5 -360.49 -344.45 185.25 -370.49 0.7916 1 0.3736anova(m2, m6, test="Chisq") # Adding B does not improve fit## Data: data
## Models:
## m2: R ~ B + (1 | X)
## m6: R ~ B + C + (1 | X)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 4 -361.70 -348.86 184.85 -369.70
## m6 5 -360.01 -343.97 185.01 -370.01 0.3117 1 0.5766continuous_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: -361.7016 -361.7016
## (Intercept) coeff: 0.3078765 tval: 36.61411
## sex_c coeff: 0.0150669 tval: 1.791825continuous_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.3798 -384.3798
## (Intercept) coeff: 0.3102547 tval: 18.70692
## sex_c coeff: 0.0182755 tval: 2.33547- 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)
source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 4-FF REMLF.R"))## [,1] [,2] [,3]
## AICs -365.1875 -363.4349 -363.2576
## models 20 30 31
## probs 0.2192431 0.09127204 0.08353253
##
## 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: 2anova(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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m20 6 -365.19 -345.96 188.59 -377.19
## m30 7 -363.43 -341.01 188.72 -377.43 0.2473 1 0.619anova(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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m20 6 -365.19 -345.96 188.59 -377.19
## m31 7 -363.26 -340.83 188.63 -377.26 0.0701 1 0.7912continuous_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: -351.7787 -363.7787
## (Intercept) coeff: 0.3272142 tval: 24.85364
## sex_c coeff: 0.0162801 tval: 1.236557
## mass_c coeff: 1.2832109 tval: 1.692175
## sex_c:mass_c coeff: -2.1490385 tval: -2.833945continuous_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.4728 -390.4728
## (Intercept) coeff: 0.3309742 tval: 16.65808
## sex_c coeff: 0.0202936 tval: 1.692638
## mass_c coeff: 1.2020658 tval: 1.744008
## sex_c:mass_c coeff: -2.2305168 tval: -3.172316- 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)
source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 1-RF + 3-FF REMLF.R"))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -287.9774 -287.7359 -287.6323 -287.5086 -287.2024 -287.0603
## models 10 9 17 14 13 1
## probs 0.1123919 0.09961121 0.09458204 0.08890941 0.07628609 0.0710534
## [,7] [,8] [,9] [,10]
## AICs -287.0587 -287.0125 -286.79 -286.6292
## models 2 15 8 3
## probs 0.07099847 0.06937735 0.0620739 0.05727863
##
## m10 R ~ B * C + (1 | X)
## m9 R ~ A * C + (1 | X)
## m17 R ~ A * B + A * C + B * C + (1 | X)
## m14 R ~ A * B + A * C + (1 | X)
## m13 R ~ B * C + A + (1 | X)
## m1 R ~ A + (1 | X)
## m2 R ~ B + (1 | X)
## m15 R ~ A * B + B * C + (1 | X)
## m8 R ~ A * B + (1 | X)
## m3 R ~ C + (1 | X)cat("Number of models that failed to converge: ", length(errors$warnings))## Number of models that failed to converge: 1anova(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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m10 6 -287.98 -270.24 149.99 -299.98
## m13 7 -287.20 -266.51 150.60 -301.20 1.225 1 0.2684continuous_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: -260.257 -272.257
## (Intercept) coeff: 0.3151502 tval: 34.78561
## body_c coeff: 0.0490906 tval: 1.351345
## wing_c coeff: -0.0205393 tval: -0.438782
## body_c:wing_c coeff: -0.0199209 tval: -2.205574continuous_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.9727 -282.9727
## (Intercept) coeff: 0.3131008 tval: 20.16806
## body_c coeff: 0.0485746 tval: 1.39824
## wing_c coeff: -0.0214165 tval: -0.4786581
## body_c:wing_c coeff: -0.0183365 tval: -2.146353- 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)
source("src/plotting-lm.R")
source("src/plotting-lm2.R")
Bursters
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]]
### 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: -196.4306 -216.4306
## (Intercept) coeff: 0.3735376 tval: 13.10017
## chamberA-1 coeff: 0.0411288 tval: 1.046403
## chamberA-3 coeff: -0.0610556 tval: -1.667595
## chamberA-4 coeff: -0.0327626 tval: -0.9814444
## chamberB-1 coeff: -0.0514394 tval: -1.146508
## chamberB-2 coeff: -0.0605026 tval: -1.764738
## chamberB-3 coeff: -0.0111004 tval: -0.3053705
## chamberB-4 coeff: -0.094332 tval: -2.458132## lmer average_speed ~ days_from_start + (1 | ID) dB
## AIC: -217.6996 -225.6996
## (Intercept) coeff: 0.3181063 tval: 19.84549
## days_from_start coeff: 0.0018077 tval: 1.437048## lmer average_speed ~ min_from_IncStart + (1 | trial_type) dB
## AIC: -214.1141 -222.1141
## (Intercept) coeff: 0.3624412 tval: 18.59292
## min_from_IncStart coeff: -0.0001127 tval: -1.980124Without 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)
source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 2-RF + 3-FF.R"))## [,1] [,2] [,3]
## AICs -232.1604 -230.8822 -229.1866
## models 18 36 54
## probs 0.5254394 0.2773097 0.1187823
##
## m18 R ~ (1 | Y)
## m36 R ~ (1 | X) + (1 | Y)
## m0 R ~ (1 | X)cat("Number of models that did not converge: ", length(errors$warnings))## Number of models that did not converge: 1anova(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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m18 3 -238.97 -230.12 122.48 -244.97
## m0 3 -236.84 -227.99 121.42 -242.84 0 0 1anova(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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -236.84 -227.99 121.42 -242.84
## m36 4 -237.56 -225.77 122.78 -245.56 2.727 1 0.09867 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1burster_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: -229.1866 -235.1866
## 1 coeff: 0.3375142 tval: 38.61569s.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)
source("src/compare_models.R")
errors <- withWarnings(model_comparisonsAIC("src/generic models-gaussian glmer 2-RF + 4-FF.R"))## [,1] [,2] [,3] [,4] [,5] [,6]
## AICs -232.1604 -230.6291 -230.3303 -229.9407 -228.8253 -228.6047
## models 113 117 226 339 230 4
## probs 0.3524708 0.1639065 0.1411633 0.1161731 0.06651277 0.05956686
##
## m113 R ~ 1 + (1 | Y)
## m117 R ~ D + (1 | Y)
## m226 R ~ 1 + (1 | X) + (1 | Y)
## m0 R ~ 1 + (1 | X)
## m230 R ~ D + (1 | X) + (1 | Y)
## m4 R ~ D + (1 | X)cat("Number of models that did not converge: ", length(errors$warnings))## Number of models that did not 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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m113 3 -238.97 -230.12 122.48 -244.97
## m117 4 -237.25 -225.45 122.62 -245.25 0.276 1 0.5994anova(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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m113 3 -238.97 -230.12 122.48 -244.97
## m226 4 -237.00 -225.21 122.50 -245.00 0.032 1 0.858anova(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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -236.94 -228.09 121.47 -242.94
## m4 4 -235.34 -223.55 121.67 -243.34 0.4 1 0.5271anova(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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m226 4 -237.00 -225.21 122.50 -245.00
## m230 5 -235.29 -220.55 122.64 -245.29 0.2873 1 0.592burster_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: -232.1604 -238.1604
## 1 coeff: 0.3392337 tval: 24.58075s.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)
source("src/compare_models.R")
model_comparisonsAIC("src/generic models-gaussian glmer 2-RF + 3-FF.R")## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00200157 (tol = 0.002, component 1)## [,1] [,2] [,3] [,4]
## AICs -232.1604 -230.3303 -229.9407 -228.2971
## models 18 36 54 19
## probs 0.4306301 0.1724659 0.1419341 0.06240025
##
## m18 R ~ (1 | Y)
## m36 R ~ (1 | X) + (1 | Y)
## m0 R ~ (1 | X)
## m19 R ~ A + (1 | Y)anova(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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m18 3 -238.97 -230.12 122.48 -244.97
## m19 4 -240.35 -228.56 124.18 -248.35 3.3842 1 0.06583 .
## ---
## 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)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m0 3 -236.94 -228.09 121.47 -242.94
## m36 4 -237.00 -225.21 122.50 -245.00 2.0608 1 0.1511burster_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: -228.2971 -236.2971
## (Intercept) coeff: 0.3393808 tval: 25.53237
## thorax_c coeff: 0.055011 tval: 1.825708burster_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: -232.1604 -238.1604
## 1 coeff: 0.3392337 tval: 24.58075s.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)
source("src/plotting-lm.R")
source("src/plotting-lm2.R")
