Data Visualizations: flight_yes_no
Anastasia Bernat
5/4/2020-1/6/2021
Winter 2020 Flight Trials
Flight Trials Winter 2020 Dataset was conducted from 2/17/2020 - 3/10/2020. Soapberry bugs were flight tested twice in the flight mill and observed from 8 AM to (5-8 PM) each day.
Loading, cleaning, and recoding the data
source_path = "~/Desktop/git_repositories/SBB-dispersal/avbernat/Rsrc/"
script_names = c("center_flight_data.R", # Re-centers data
"clean_flight_data.R", # Loads and cleans data
"unique_flight_data.R",
"get_warnings.R",
"compare_models.R",
"regression_output.R", # Cleans regression outputs; prints in color or B&W
"AICprobabilities.R",
"get_Akaike_weights.R",
"categorize_binomial_data.R",
"categorize_binomial_data-MF.R",
"categorize_binomial_data-MF2.R")
for (script in script_names) {
path = paste0(source_path, script)
source(path)
}
output_col = FALSE # Change to TRUE if working in Base R or RStudio; FALSE if generating an HTML
data <- read_flight_data("data/all_flight_data-Winter2020.csv")## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object length## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))## Warning in d$latitude - sym_zone: longer object length is not a multiple of
## shorter object lengthdata_all <- data[[1]]
data_tested <- data[[2]]Plotting and Categorizing Functions
categorize_data: It groups explanatory values into categories and calculates the sample of proportion of flight for each category
Experimental Set-Up
Incubator Start and Trial Duration
How did/does minutes since the incubator started and trial length (min) impact whether a bug flew or not? The answer is they don’t, which is good. That means that our experimental set-up did not impact if a bug flew or not.
cat_data1 <- categorize_data(data_tested, all_of("min_from_IncStart"), all_of("flew_b"), 30, 0, 30)## min_from_IncStart 30 Cases: 72 Sample Proportion: 0.5694444
## min_from_IncStart 60 Cases: 84 Sample Proportion: 0.6190476
## min_from_IncStart 90 Cases: 66 Sample Proportion: 0.4848485
## min_from_IncStart 120 Cases: 53 Sample Proportion: 0.5471698
## min_from_IncStart 150 Cases: 43 Sample Proportion: 0.3488372
## min_from_IncStart 180 Cases: 42 Sample Proportion: 0.547619
## min_from_IncStart 210 Cases: 43 Sample Proportion: 0.5348837
## min_from_IncStart 240 Cases: 31 Sample Proportion: 0.7096774
## min_from_IncStart 270 Cases: 17 Sample Proportion: 0.5882353
## min_from_IncStart 300 Cases: 19 Sample Proportion: 0.6315789
## min_from_IncStart 330 Cases: 22 Sample Proportion: 0.6363636
## min_from_IncStart 360 Cases: 26 Sample Proportion: 0.5769231
## min_from_IncStart 390 Cases: 21 Sample Proportion: 0.7142857
## min_from_IncStart 420 Cases: 28 Sample Proportion: 0.5
## min_from_IncStart 450 Cases: 14 Sample Proportion: 0.7142857
## min_from_IncStart 480 Cases: 22 Sample Proportion: 0.6363636
## min_from_IncStart 510 Cases: 7 Sample Proportion: 0.2857143cat_data2 <- categorize_data(data_tested, all_of("minute_duration"), all_of("flew_b"), 30, 0, 30)## minute_duration 30 Cases: 9 Sample Proportion: 0.1111111
## minute_duration 60 Cases: 431 Sample Proportion: 0.3897912
## minute_duration 90 Cases: 31 Sample Proportion: 1
## minute_duration 120 Cases: 25 Sample Proportion: 1
## minute_duration 150 Cases: 9 Sample Proportion: 1
## minute_duration 180 Cases: 14 Sample Proportion: 1
## minute_duration 210 Cases: 7 Sample Proportion: 1
## minute_duration 240 Cases: 6 Sample Proportion: 1
## minute_duration 270 Cases: 3 Sample Proportion: 1
## minute_duration 300 Cases: 6 Sample Proportion: 1
## minute_duration 330 Cases: 4 Sample Proportion: 1
## minute_duration 360 Cases: 3 Sample Proportion: 1
## minute_duration 390 Cases: 5 Sample Proportion: 1
## minute_duration 420 Cases: 8 Sample Proportion: 1
## minute_duration 450 Cases: 6 Sample Proportion: 1
## minute_duration 480 Cases: 6 Sample Proportion: 1
## minute_duration 510 Cases: 5 Sample Proportion: 1
## minute_duration 540 Cases: 2 Sample Proportion: 1
## minute_duration 570 Cases: 3 Sample Proportion: 1
## minute_duration 600 Cases: 2 Sample Proportion: 1
## minute_duration 630 Cases: 3 Sample Proportion: 1
## minute_duration 660 Cases: 4 Sample Proportion: 1
## minute_duration 690 Cases: 5 Sample Proportion: 1
## minute_duration 720 Cases: 2 Sample Proportion: 1
## minute_duration 750 Cases: 6 Sample Proportion: 1
## minute_duration 780 Cases: 1 Sample Proportion: 1
## minute_duration 810 Cases: 8 Sample Proportion: 1# Plotting observed proportions of yes flew by minutes from incubation start
p1 <- as.grob(expression(plot(sample_prop ~ min_from_IncStart, data=cat_data1,
ylab="Sample Proportion of Yes Flew", xlab="Minutes From Incubtation Start (min) ",
main="Observed proportions of yes flew by minutes from incubation start")
))
# Plotting observed proportions of yes flew by minutes from trial start
p2 <- as.grob(expression(
plot(sample_prop ~ minute_duration, data=cat_data2,
ylab="Sample Proportion of Yes Flew", xlab="Minutes From Trial Start (min)",
main="Observed proportions of yes flew by minutes from trial start") %>%
abline(v=120, col="red")
))
grid.arrange(p1,p2, ncol=2)
Testing Start Time
How did/does the time of day in which a bug was tested impact how long the bug flew for/whether the bug flew or not? When all bugs were considered, the effect of start time was marginal. When only bugs that flew were considered, the effect of start time was significant. From this we see that time was a limitation on our experimental set-up.
par(mfrow=c(1,2))
### Plot trial duration (min) vs. time_start | all bugs
x = chron(times = data_tested$time_start)
y = data_tested$recording_duration / (60*60)
fit <- lm(y~x, data=data_tested)
summary(fit) ##
## Call:
## lm(formula = y ~ x, data = data_tested)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8061 -1.4060 -1.1175 -0.5736 11.2788
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.9330 0.5573 5.263 1.97e-07 ***
## x -2.2654 1.1912 -1.902 0.0577 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.89 on 612 degrees of freedom
## Multiple R-squared: 0.005875, Adjusted R-squared: 0.004251
## F-statistic: 3.617 on 1 and 612 DF, p-value: 0.05767cf <- round(coef(fit), 2)
eq <- paste0("portion_flew = ", cf[1],
ifelse(sign(cf[2])==1, " + ", " - "), abs(cf[2]), "*time_start ")
plot(x,y, xlab= "Starting Time", ylab="Duration (hr)")
abline(coef(fit)[1:2], col="blue")
### Plot trial duration (min) vs. time_start | bugs that flew only
all_flew <- filter(data_all, flew_b == 1)
x = chron(times = all_flew$time_start)
y = all_flew$recording_duration / (60*60)
fit <- lm(y~x, data=all_flew)
summary(fit) ##
## Call:
## lm(formula = y ~ x, data = all_flew)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8203 -2.3285 -1.6447 0.4281 10.0744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8432 0.9162 5.286 2.23e-07 ***
## x -4.2197 1.9521 -2.162 0.0313 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.547 on 341 degrees of freedom
## Multiple R-squared: 0.01352, Adjusted R-squared: 0.01062
## F-statistic: 4.672 on 1 and 341 DF, p-value: 0.03135cf <- round(coef(fit), 2)
eq <- paste0("portion_flew = ", cf[1],
ifelse(sign(cf[2])==1, " + ", " - "), abs(cf[2]), "*time_start ")
plot(x,y, xlab= "Starting Time", ylab="Duration (hr)")
abline(coef(fit)[1:2], col="blue") 
c1=1.4
c2=1.2
c3=1.2
data_flew = data_tested %>%
filter(flew_b==1)
d = data_flew
data_flew$start_time = chron(times = data_flew$time_start)
data_flew$recording_dur_hr = data_flew$recording_duration / (60*60)
x=data_flew$start_time
y=data_flew$recording_dur_hr
m <- lm(recording_dur_hr ~ start_time, data=data_flew)
x.seq = seq(min(x) - sd(x), max(x) + sd(x), length.out=100)
m <- lm(y ~ x, data=data_flew)
prd <- data.frame(x=x.seq) # newdata
err <- predict(m, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit
mu_ci <- t(matrix(c(prd$lci,prd$uci), ncol=2))
colfunc <- colorRampPalette(c("red", "orange"))
plot(data_flew$recording_dur_hr ~ data_flew$start_time,
pch=16,
xlab = "Trial Start Time (HR:MIN)",
ylab= "Recording Duration (HR)",
cex=c2,
cex.lab=c3,
cex.axis=c2,
ylim=c(0,14),
#col = alpha("black", 0.7),
#bty="n"
)
abline(m, lty=2)
shade(mu_ci, lim = prd$x)
pval <- paste0("p = ", round(summary(m)$coefficients[8],3), "*")
text(chron(times = "14:15:15"), 4, pval, cex=c3)
#title("Time of Day Tested", adj = 0.05, line = 0, cex.main=c1)
#title("(a)", adj = 0.05, line = 0, cex.main=1.4)
summary(lm(recording_dur_hr ~ start_time, data=data_flew))##
## Call:
## lm(formula = recording_dur_hr ~ start_time, data = data_flew)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8203 -2.3285 -1.6447 0.4281 10.0744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8432 0.9162 5.286 2.23e-07 ***
## start_time -4.2197 1.9521 -2.162 0.0313 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.547 on 341 degrees of freedom
## Multiple R-squared: 0.01352, Adjusted R-squared: 0.01062
## F-statistic: 4.672 on 1 and 341 DF, p-value: 0.03135Days from Experiment Start
cat_data_days <- categorize_data(data_tested, all_of("days_from_start"), all_of("flew_b"), 1, 0, 1)## days_from_start 1 Cases: 37 Sample Proportion: 0.7027027
## days_from_start 2 Cases: 21 Sample Proportion: 0.7142857
## days_from_start 3 Cases: 35 Sample Proportion: 0.6857143
## days_from_start 4 Cases: 31 Sample Proportion: 0.5483871
## days_from_start 5 Cases: 22 Sample Proportion: 0.5
## days_from_start 8 Cases: 43 Sample Proportion: 0.4883721
## days_from_start 9 Cases: 41 Sample Proportion: 0.5365854
## days_from_start 10 Cases: 41 Sample Proportion: 0.7073171
## days_from_start 11 Cases: 30 Sample Proportion: 0.5666667
## days_from_start 12 Cases: 33 Sample Proportion: 0.6060606
## days_from_start 16 Cases: 47 Sample Proportion: 0.5106383
## days_from_start 17 Cases: 52 Sample Proportion: 0.5961538
## days_from_start 18 Cases: 52 Sample Proportion: 0.4807692
## days_from_start 19 Cases: 28 Sample Proportion: 0.6071429
## days_from_start 22 Cases: 63 Sample Proportion: 0.3968254x=cat_data_days$days_from_start
y=cat_data_days$sample_prop
m <- lm(sample_prop ~ days_from_start, data=cat_data_days)
x.seq = seq(min(x) - sd(x), max(x) + sd(x), length.out=100)
m <- lm(y ~ x, data=cat_data_days)
prd <- data.frame(x=x.seq) # newdata
err <- predict(m, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit
mu_ci <- t(matrix(c(prd$lci,prd$uci), ncol=2))
colfunc <- colorRampPalette(c("red", "orange"))
plot(cat_data_days$sample_prop ~ cat_data_days$days_from_start,
pch=16,
xlab = "Day From Start",
ylab= "Flight Probability",
cex=c2,
cex.lab=c3,
cex.axis=c2,
xlim=c(0,25),
ylim=c(0.4, 0.8),
#col = alpha("black", 0.7),
#bty="n"
)
abline(m, lty=2)
shade(mu_ci, lim = prd$x)
pval <- paste0("p = ", round(summary(m)$coefficients[8],3), "*")
text(15, 0.65, pval, cex=c3)
#title("Age", adj = 0.05, line = 0, cex.main=c1)
#title("(a)", adj = 0.05, line = 0, cex.main=1.4)
summary(glm(flew_b ~ days_from_start_c, data=data_tested, family="binomial"))##
## Call:
## glm(formula = flew_b ~ days_from_start_c, family = "binomial",
## data = data_tested)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.460 -1.226 0.932 1.071 1.235
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.23928 0.08192 2.921 0.00349 **
## days_from_start_c -0.03545 0.01184 -2.995 0.00274 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 842.72 on 613 degrees of freedom
## Residual deviance: 833.62 on 612 degrees of freedom
## AIC: 837.62
##
## Number of Fisher Scoring iterations: 4summary(lm(sample_prop ~ days_from_start, data=cat_data_days))##
## Call:
## lm(formula = sample_prop ~ days_from_start, data = cat_data_days)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.120275 -0.065780 -0.005571 0.060198 0.127073
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.660306 0.040058 16.484 4.29e-10 ***
## days_from_start -0.008006 0.003250 -2.464 0.0285 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08194 on 13 degrees of freedom
## Multiple R-squared: 0.3183, Adjusted R-squared: 0.2658
## F-statistic: 6.069 on 1 and 13 DF, p-value: 0.02848Biological Features
Host, Sex, Distance From Sympatric Zone
summary <-aggregate(flew_b~sex*host_plant*sym_dist, data=data_all, FUN=mean)
summary## sex host_plant sym_dist flew_b
## 1 F K. elegans 0.0000000 0.33333333
## 2 M K. elegans 0.0000000 0.30555556
## 3 F C. corindum 0.2046150 0.44444444
## 4 M C. corindum 0.2046150 0.50000000
## 5 F C. corindum 0.2047546 0.00000000
## 6 M C. corindum 0.2047546 0.00000000
## 7 F C. corindum 0.2178735 0.29629630
## 8 M C. corindum 0.2178735 0.76470588
## 9 F C. corindum 0.2178735 0.00000000
## 10 M C. corindum 0.2178735 0.00000000
## 11 F C. corindum 0.2348005 0.10000000
## 12 M C. corindum 0.2348005 0.95000000
## 13 F C. corindum 0.2638382 0.18181818
## 14 M C. corindum 0.2638382 0.52173913
## 15 F C. corindum 0.2962120 0.25000000
## 16 M C. corindum 0.2962120 0.59259259
## 17 F C. corindum 0.3092689 0.00000000
## 18 M C. corindum 0.3092689 0.61666667
## 19 F C. corindum 0.3157918 0.25000000
## 20 M C. corindum 0.3157918 0.51515152
## 21 F C. corindum 0.3629042 0.52631579
## 22 M C. corindum 0.3629042 0.63636364
## 23 F C. corindum 0.3913469 0.20000000
## 24 M C. corindum 0.3913469 0.50000000
## 25 F C. corindum 0.5177872 0.00000000
## 26 M C. corindum 0.5177872 0.58620690
## 27 F C. corindum 0.5271220 0.30303030
## 28 M C. corindum 0.5271220 0.70000000
## 29 F K. elegans 1.8072690 0.37500000
## 30 M K. elegans 1.8072690 0.76000000
## 31 F K. elegans 1.8072980 0.00000000
## 32 M K. elegans 1.8072980 0.00000000
## 33 F K. elegans 2.4053783 0.41666667
## 34 M K. elegans 2.4053783 0.62500000
## 35 F K. elegans 3.3049815 0.09090909
## 36 M K. elegans 3.3049815 0.46428571
## 37 F K. elegans 4.1706402 0.00000000
## 38 M K. elegans 4.1706402 0.50000000
## 39 F K. elegans 4.1706604 0.00000000
## 40 M K. elegans 4.1706604 0.00000000par(mar=c(6.5, 5.5, 5.5, 9.5), xpd=TRUE)
plot(summary$flew_b~summary$sym_dist,
col=c(1,2)[as.factor(summary$sex)], # Female = Black, Male = Red
pch=c(19,22)[as.factor(summary$host_plant)], # K. elegans = Squares C.corindum = circles
main="Observed Data",
xlab = "Distance from Sympatric Zone (°)",
ylab= "Proportion Flew")
legend("topright",
legend = c("C. corindum and F", "C. corindum and M", "K. elegans and F" , "K. elegans and M"),
inset=c(-0.23,0.2),
col=c(1,2),
pch = c(19,19, 22,22),
title="Groups")
Host, Sex, Distance From Sympatric Zone, Mass
data_mass <- data_all %>%
filter(!is.na(mass))
cat_data5 <- categorize_data_MF(data_mass, all_of("mass"), all_of("sex_c"),
all_of("host_c"), all_of("flew_b"), 0.010, 0.015, 0.025)## Mass 0.035 Sex: M Host Plant: 1 Cases: 2 Sample Proportion: 0
## Mass 0.045 Sex: M Host Plant: -1 Cases: 2 Sample Proportion: 0.5
## Mass 0.055 Sex: F Host Plant: 1 Cases: 8 Sample Proportion: 0.25
## Mass 0.055 Sex: M Host Plant: 1 Cases: 11 Sample Proportion: 0.3636364
## Mass 0.055 Sex: F Host Plant: -1 Cases: 22 Sample Proportion: 0.4545455
## Mass 0.055 Sex: M Host Plant: -1 Cases: 51 Sample Proportion: 0.6666667
## Mass 0.065 Sex: F Host Plant: 1 Cases: 10 Sample Proportion: 0.5
## Mass 0.065 Sex: M Host Plant: 1 Cases: 1 Sample Proportion: 0
## Mass 0.065 Sex: F Host Plant: -1 Cases: 20 Sample Proportion: 0.4
## Mass 0.065 Sex: M Host Plant: -1 Cases: 10 Sample Proportion: 0.5
## Mass 0.075 Sex: F Host Plant: 1 Cases: 12 Sample Proportion: 0.5833333
## Mass 0.085 Sex: F Host Plant: -1 Cases: 23 Sample Proportion: 0.4782609
## Mass 0.105 Sex: F Host Plant: 1 Cases: 1 Sample Proportion: 1
## Mass 0.115 Sex: F Host Plant: -1 Cases: 19 Sample Proportion: 0.1578947
## Mass 0.135 Sex: F Host Plant: 1 Cases: 1 Sample Proportion: 1
## Mass 0.145 Sex: F Host Plant: -1 Cases: 5 Sample Proportion: 0
## Mass 0.205 Sex: F Host Plant: -1 Cases: 1 Sample Proportion: 0par(mar=c(6.5, 5.5, 5.5, 9.5), xpd=TRUE)
plot(cat_data5$sample_prop~cat_data5$mass,
col=c(1,2)[as.factor(cat_data5$sex_c)], # Female = Black, Male = Red
pch=c(19,22)[as.factor(cat_data5$host_c)], # K. elegans = Squares, C.corindum = circles
main="Observed Data",
xlab = "Mass (g)",
ylab= "Proportion Flew")
legend("topright",
legend = c("C. corindum and F", "C. corindum and M", "K. elegans and F" , "K. elegans and M"),
inset=c(-0.23,0.2),
col=c(1,2),
pch = c(19,19, 22,22),
title="Groups")
Mass
# Plot All
data_mass <- data_all %>%
filter(!is.na(mass))
cat_data3 <- categorize_data(data_mass, all_of("mass"), all_of("flew_b"), 0.010, 0.015, 0.025) ## mass 0.025 Cases: 4 Sample Proportion: 0.25
## mass 0.035 Cases: 123 Sample Proportion: 0.6747967
## mass 0.045 Cases: 206 Sample Proportion: 0.7281553
## mass 0.055 Cases: 92 Sample Proportion: 0.5434783
## mass 0.065 Cases: 41 Sample Proportion: 0.4390244
## mass 0.075 Cases: 35 Sample Proportion: 0.5142857
## mass 0.085 Cases: 31 Sample Proportion: 0.1935484
## mass 0.095 Cases: 28 Sample Proportion: 0.2142857
## mass 0.105 Cases: 20 Sample Proportion: 0.2
## mass 0.115 Cases: 13 Sample Proportion: 0.3076923
## mass 0.125 Cases: 7 Sample Proportion: 0
## mass 0.135 Cases: 6 Sample Proportion: 0.1666667
## mass 0.145 Cases: 3 Sample Proportion: 0.3333333
## mass 0.155 Cases: 1 Sample Proportion: 0
## mass 0.185 Cases: 1 Sample Proportion: 0# Plot Trial 1
data_T1 <-data_all[data_all$trial_type=="T1",]
data_mass_T1 <- data_T1 %>%
filter(!is.na(mass))
cat_data4 <- categorize_data(data_mass_T1, all_of("mass"), all_of("flew_b"), 0.010, 0.015, 0.025) ## mass 0.025 Cases: 4 Sample Proportion: 0.25
## mass 0.035 Cases: 79 Sample Proportion: 0.7468354
## mass 0.045 Cases: 104 Sample Proportion: 0.7596154
## mass 0.055 Cases: 52 Sample Proportion: 0.5384615
## mass 0.065 Cases: 21 Sample Proportion: 0.5238095
## mass 0.075 Cases: 19 Sample Proportion: 0.5789474
## mass 0.085 Cases: 12 Sample Proportion: 0.25
## mass 0.095 Cases: 11 Sample Proportion: 0.1818182
## mass 0.105 Cases: 9 Sample Proportion: 0.2222222
## mass 0.115 Cases: 8 Sample Proportion: 0.375
## mass 0.125 Cases: 4 Sample Proportion: 0
## mass 0.135 Cases: 4 Sample Proportion: 0.25
## mass 0.155 Cases: 1 Sample Proportion: 0
## mass 0.185 Cases: 1 Sample Proportion: 0# Plot Trial 2
data_T2 <-data_all[data_all$trial_type=="T2",]
data_mass_T2 <- data_T2 %>%
filter(!is.na(mass))
cat_data5 <- categorize_data(data_mass_T2, all_of("mass"), all_of("flew_b"), 0.010, 0.015, 0.025) ## mass 0.035 Cases: 44 Sample Proportion: 0.5454545
## mass 0.045 Cases: 102 Sample Proportion: 0.6960784
## mass 0.055 Cases: 40 Sample Proportion: 0.55
## mass 0.065 Cases: 20 Sample Proportion: 0.35
## mass 0.075 Cases: 16 Sample Proportion: 0.4375
## mass 0.085 Cases: 19 Sample Proportion: 0.1578947
## mass 0.095 Cases: 17 Sample Proportion: 0.2352941
## mass 0.105 Cases: 11 Sample Proportion: 0.1818182
## mass 0.115 Cases: 5 Sample Proportion: 0.2
## mass 0.125 Cases: 3 Sample Proportion: 0
## mass 0.135 Cases: 2 Sample Proportion: 0
## mass 0.145 Cases: 3 Sample Proportion: 0.3333333par(mfrow=c(1,2))
# Plot Trial 1
fit<-lm(flew_b~mass, data=data_T1)
coeff <- coefficients(summary(fit))
mass_seq <- seq(from=0.0, to=0.18, length.out=30)
mu <- link(fit, data=data.frame(mass=mass_seq)) ## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.mu.PI <- apply(mu,2,PI, prob=0.89)
eq <- paste0("prop flew = ", round(coeff[1],2),
ifelse(sign(coeff[2])==1, " + ", " - "),
abs(round(coeff[2],2)), "*mass")
plot(as.matrix(cat_data4[1]),
as.matrix(cat_data4[3]),
ylab="Sample Proportion of Yes Flew",
xlab="Mass (g)",
main="Trial 1", pch=16)
shade(mu.PI, mass_seq)
abline(coeff[1], coeff[2], lty=2)
text(0.13, 0.68, eq)
# Plot Trial 2
fit<-lm(flew_b~mass, data=data_T2)
coeff <- coefficients(summary(fit))
mass_seq <- seq(from=0.0, to=0.18, length.out=30)
mu <- link(fit, data=data.frame(mass=mass_seq))## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.mu.PI <- apply(mu,2,PI, prob=0.89)
eq <- paste0("prop flew = ", round(coeff[1],2),
ifelse(sign(coeff[2])==1, " + ", " - "),
abs(round(coeff[2],2)), "*mass")
plot(as.matrix(cat_data5[1]),
as.matrix(cat_data5[3]),
ylab="Sample Proportion of Yes Flew",
xlab="Mass (g)",
main="Trial 2", pch=16)
shade(mu.PI, mass_seq)
abline(coeff[1], coeff[2], lty=2)
text(0.105, 0.6, eq)
Females
data_fem <- data_all[data_all$sex=="F",]Mass, Eggs
pf1 <- as.grob(expression(
plot(data_fem$flew_b, data_fem$total_eggs, xlab="flight response", ylab="total eggs")))
# seems like if you generally laid less eggs you flew, but a lot of females laid few and also didn't fly
pf2 <- as.grob(expression(
plot(data_fem$mass, data_fem$total_eggs, xlab="mass (g)", ylab="total eggs")))
pf3 <- as.grob(expression(
plot(data_fem$sym_dist,data_fem$total_eggs, xlab="total eggs", ylab="distance from sympatric zone")))
grid.arrange(pf1,pf2,pf3, ncol=3)
# remove missing mass
data_fem_temp = data_fem %>%
filter(!is.na(mass))
x=data_fem_temp$mass
y=data_fem_temp$total_eggs
m <- lm(total_eggs ~ mass, data=data_fem_temp)
x.seq = seq(min(x) - sd(x), max(x) + sd(x), length.out=100)
m <- lm(y ~ x, data=data_fem_temp)
prd <- data.frame(x=x.seq) # newdata
err <- predict(m, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit
mu_ci <- t(matrix(c(prd$lci,prd$uci), ncol=2))
colfunc <- colorRampPalette(c("red", "orange"))
plot(data_fem_temp$total_eggs ~ data_fem_temp$mass,
pch=16,
xlab = "mass (g)",
ylab= "total eggs",
cex=c2,
cex.lab=c3,
cex.axis=c2,
#ylim=c(0,14),
#col = alpha("black", 0.7),
#bty="n"
)
abline(m, lty=2)
shade(mu_ci, lim = prd$x)
pval <- paste0("p = ", round(summary(m)$coefficients[8],10), "*")
text(0.16, 150, pval, cex=c3)
#title("Time of Day Tested", adj = 0.05, line = 0, cex.main=c1)
#title("(a)", adj = 0.05, line = 0, cex.main=1.4)
summary(lm(total_eggs ~ mass, data=data_fem_temp))##
## Call:
## lm(formula = total_eggs ~ mass, data = data_fem_temp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -198.66 -61.16 -14.70 61.67 236.08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -34.45 22.75 -1.515 0.132
## mass 1721.31 266.48 6.459 9.02e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 85.47 on 185 degrees of freedom
## (26 observations deleted due to missingness)
## Multiple R-squared: 0.184, Adjusted R-squared: 0.1796
## F-statistic: 41.72 on 1 and 185 DF, p-value: 9.022e-10path = paste0(source_path, "delta-females.R")
source(path)## Warning in if (egg_data_T1$ID %in% egg_data_T2$ID) {: the condition has length >
## 1 and only the first element will be used## ID flew_delta mass_diff flew_count egg_diff
## 21 53 0 -0.039 2 -1
## 34 214 0 -0.038 2 1
## 17 319 0 -0.034 2 1
## 32 101 0 -0.026 0 1
## 14 155 0 -0.026 0 -1
## 26 407 0 -0.020 0 -1
## 7 147 0 -0.014 2 1
## 11 151 0 -0.014 2 1
## 5 130 -1 -0.008 1 1
## 13 132 0 -0.006 0 1
## 2 37 0 -0.004 0 1
## 8 109 0 -0.002 0 1
## 20 233 0 -0.002 0 -1
## 28 438 0 0.003 2 1
## 23 96 0 0.006 0 -1
## 12 149 0 0.008 0 -1
## 25 140 0 0.010 2 1
## 15 375 -1 0.015 1 1
## 6 164 0 0.016 2 1
## 24 211 0 0.017 0 1
## 4 63 -1 0.018 1 1
## 33 116 0 0.019 2 1
## 10 269 -1 0.021 1 1
## 22 161 -1 0.022 1 1
## 1 318 -1 0.023 1 1
## 19 47 -1 0.027 1 1
## 16 271 -1 0.030 1 1
## 29 424 -1 0.030 1 1
## 31 97 0 0.031 0 1
## 30 20 -1 0.039 1 1
## 18 7 0 0.041 2 1
## 27 26 0 0.042 2 1
## 3 92 0 0.054 2 1
## 9 12 -1 0.064 1 1Three things noticed:
12 of the 34 rows flew twice (in both trials).
12 of the 34 did not flew at all (in either trial).
there were no cases in which the flight diff = 1, so no female bugs flew in T2 but not T1. But 10 cases of the 34 that flew in T1 but not T2.
cat_data10 <- categorize_data_MF2(egg_data_final, "mass_diff", "flew_count", "flew_delta", 0.005, 0.010, -0.040) ## Mass -0.03 Flight Count: 2 Cases: 3 Sample Proportion: 0
## Mass -0.02 Flight Count: 0 Cases: 2 Sample Proportion: 0
## Mass -0.01 Flight Count: 0 Cases: 1 Sample Proportion: 0
## Mass -0.01 Flight Count: 2 Cases: 2 Sample Proportion: 0
## Mass 0.00 Flight Count: 1 Cases: 5 Sample Proportion: -0.2
## Mass 0.01 Flight Count: 0 Cases: 2 Sample Proportion: 0
## Mass 0.01 Flight Count: 2 Cases: 1 Sample Proportion: 0
## Mass 0.02 Flight Count: 1 Cases: 3 Sample Proportion: -0.6666667
## Mass 0.02 Flight Count: 2 Cases: 3 Sample Proportion: 0
## Mass 0.03 Flight Count: 1 Cases: 6 Sample Proportion: -1
## Mass 0.04 Flight Count: 0 Cases: 2 Sample Proportion: -0.5
## Mass 0.05 Flight Count: 2 Cases: 2 Sample Proportion: 0
## Mass 0.06 Flight Count: 2 Cases: 1 Sample Proportion: 0
## Mass 0.07 Flight Count: 1 Cases: 1 Sample Proportion: -1gf_point(sample_prop ~ mass_diff, data=cat_data10, col=~flew_count,
xlab = "Mass Differences between T1 and T2 (g)", ylab= "Proportion Flew Over Trials",
title= "Observed proportions of yes flew by mass differences between trials")
fit <- lm(sample_prop~mass_diff, data=cat_data10)
coeff <- coefficients(summary(fit))
eq <- paste0("prop_flew = ", round(coeff[1],3),
ifelse(sign(coeff[2])==1, " + ", " - "),
abs(round(coeff[2],3)), "*mass_diff")
mass_diff.seq <- seq(from=-0.04, to=0.08, length.out=30)
mu <- link(fit, data=data.frame(mass_diff=mass_diff.seq)) # gives warning## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.mu.PI <- apply(mu,2,PI, prob=0.89) # maybe change to 95 later
plot(as.matrix(cat_data10[1]),
as.matrix(cat_data10[3]),
#col=c(3,2,1)[as.factor(cat_data10$flew_count)], # green = 2, black = -1
#col=c(1,2)[as.factor(summary$sex)], # Female = Black, Male = Red
pch=c(21,22,19)[as.factor(cat_data10$flew_count)],
#col=c(1,2)[as.factor(cat_data11$egg_diff)], # there would be
ylim = c(-1,1),
ylab="Proportion Flew Over Trials",
xlab="Mass Differences between T1 to T2 (g)",
main="Observed proportions of yes flew by mass differences between trials",
cex=0.9)
abline(coeff[1], coeff[2], col="darkgrey")
mtext(eq, side=4, cex=0.9)
shade(mu.PI, mass_diff.seq)
abline(v=0, col="mediumpurple3", lty=2)
text(0.011, 0.95, c("gained mass"),col="mediumpurple4")
text(-0.009, 0.963, c("lost mass"), col="mediumpurple4")
abline(h=0, col="darkolivegreen4", lty=2)
text(0.060, 0.127, c("flew in T2 but not T1"),col="darkolivegreen4")
text(0.060, -0.127, c("flew in T1 but not T2"),col="darkolivegreen4")
legend("topright",
#legend = c("Both Trials","One Trial", "No Trial"),
legend = c("2", "1", "0"),
#inset=c(-0.03,0.02),
pch = c(19,22,21),
title="Flight Count")
egg_data_summary## ID flew_delta mass_diff flew_count egg_diff
## 21 53 0 -0.039 2 -1
## 34 214 0 -0.038 2 1
## 17 319 0 -0.034 2 1
## 32 101 0 -0.026 0 1
## 14 155 0 -0.026 0 -1
## 26 407 0 -0.020 0 -1
## 7 147 0 -0.014 2 1
## 11 151 0 -0.014 2 1
## 5 130 -1 -0.008 1 1
## 13 132 0 -0.006 0 1
## 2 37 0 -0.004 0 1
## 8 109 0 -0.002 0 1
## 20 233 0 -0.002 0 -1
## 28 438 0 0.003 2 1
## 23 96 0 0.006 0 -1
## 12 149 0 0.008 0 -1
## 25 140 0 0.010 2 1
## 15 375 -1 0.015 1 1
## 6 164 0 0.016 2 1
## 24 211 0 0.017 0 1
## 4 63 -1 0.018 1 1
## 33 116 0 0.019 2 1
## 10 269 -1 0.021 1 1
## 22 161 -1 0.022 1 1
## 1 318 -1 0.023 1 1
## 19 47 -1 0.027 1 1
## 16 271 -1 0.030 1 1
## 29 424 -1 0.030 1 1
## 31 97 0 0.031 0 1
## 30 20 -1 0.039 1 1
## 18 7 0 0.041 2 1
## 27 26 0 0.042 2 1
## 3 92 0 0.054 2 1
## 9 12 -1 0.064 1 1plot(egg_diff ~ mass_diff, data=egg_data_summary)
cat_data11 <- categorize_data_MF2(egg_data_final, "mass_diff", "egg_diff", "flew_delta", 0.005, 0.010, -0.040) ## Mass -0.03 Flight Count: -1 Cases: 3 Sample Proportion: 0
## Mass -0.03 Flight Count: NA Cases: 3 Sample Proportion: 0
## Mass -0.02 Flight Count: 1 Cases: 2 Sample Proportion: 0
## Mass -0.02 Flight Count: NA Cases: 2 Sample Proportion: 0
## Mass -0.01 Flight Count: -1 Cases: 3 Sample Proportion: 0
## Mass -0.01 Flight Count: NA Cases: 3 Sample Proportion: 0
## Mass 0.00 Flight Count: 1 Cases: 5 Sample Proportion: -0.2
## Mass 0.00 Flight Count: NA Cases: 5 Sample Proportion: 0
## Mass 0.01 Flight Count: 1 Cases: 3 Sample Proportion: 0
## Mass 0.01 Flight Count: NA Cases: 3 Sample Proportion: 0
## Mass 0.02 Flight Count: 1 Cases: 6 Sample Proportion: -0.3333333
## Mass 0.02 Flight Count: NA Cases: 6 Sample Proportion: 0
## Mass 0.03 Flight Count: 1 Cases: 6 Sample Proportion: -1
## Mass 0.03 Flight Count: NA Cases: 6 Sample Proportion: 0
## Mass 0.04 Flight Count: 1 Cases: 2 Sample Proportion: -0.5
## Mass 0.04 Flight Count: NA Cases: 2 Sample Proportion: 0
## Mass 0.05 Flight Count: 1 Cases: 2 Sample Proportion: 0
## Mass 0.05 Flight Count: NA Cases: 2 Sample Proportion: 0
## Mass 0.06 Flight Count: 1 Cases: 1 Sample Proportion: 0
## Mass 0.06 Flight Count: NA Cases: 1 Sample Proportion: 0
## Mass 0.07 Flight Count: 1 Cases: 1 Sample Proportion: -1
## Mass 0.07 Flight Count: NA Cases: 1 Sample Proportion: 0gf_point(sample_prop ~ mass_diff, data=cat_data11, col=~egg_diff,
xlab = "Mass Differences between T1 and T2 (g)", ylab= "Proportion Flew Over Trials",
title= "Observed proportions of yes flew by mass differences between trials")
fit <- lm(sample_prop~mass_diff, data=cat_data11)
coeff <- coefficients(summary(fit))
eq <- paste0("portion_flew = ", round(coeff[1],3),
ifelse(sign(coeff[2])==1, " + ", " - "),
abs(round(coeff[2],3)), "*mass_diff")
mass_diff.seq <- seq(from=-0.04, to=0.08, length.out=30)
mu <- link(fit, data=data.frame(mass_diff=mass_diff.seq)) # gives warning## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.mu.PI <- apply(mu,2,PI, prob=0.89) # maybe change to 95 later
plot(as.matrix(cat_data11[1]),
as.matrix(cat_data11[3]),
#col=c(3,2,1)[as.factor(cat_data11$flew_count)], # green = 2, black = -1
#col=c(1,2)[as.factor(summary$sex)], # Female = Black, Male = Red
pch=c(21,22)[as.factor(cat_data11$egg_diff)],
ylim = c(-1,1),
ylab="Proportion Flew Over Trials",
xlab="Mass Differences between T1 to T2 (g)",
main="Observed proportions of yes flew by mass differences between trials",
cex=0.9)
abline(coeff[1], coeff[2], col="blue")
mtext(eq, side=4, cex=0.9)
shade(mu.PI, mass_diff.seq)
abline(v=0, col="mediumpurple3", lty=2)
text(0.01, 0.95, c("gained mass"),col="mediumpurple4")
text(-0.009, 0.963, c("lost mass"), col="mediumpurple4")
abline(h=0, col="darkolivegreen4", lty=2)
text(0.060, 0.087, c("flew in T2 but not T1"),col="darkolivegreen4")
text(0.060, -0.087, c("flew in T1 but not T2"),col="darkolivegreen4")
legend("topright",
#legend = c("Both Trials","One Trial", "No Trial"),
legend = c("T2 not T1", "T1 not T2"),
#inset=c(-0.03,0.02),
pch = c(22,21),
title="Egg Event")
fit <- lm(flew_delta~mass_diff, data=egg_data_final)
p <- paste0("p = ", round(summary(fit)[[4]][8],3), "*")
coeff <- coefficients(summary(fit))
eq <- paste0("portion_flew = ", round(coeff[1],3),
ifelse(sign(coeff[2])==1, " + ", " - "),
abs(round(coeff[2],3)), "*mass_diff")
mass_diff.seq <- seq(from=-0.045, to=0.07, length.out=30)
mu <- link(fit, data=data.frame(mass_diff=mass_diff.seq)) # gives warning## Warning in .local(fit, data, n, ...): The link method for lm fits is not
## complete and not tested. Proceed with abundant caution.mu.PI <- apply(mu,2,PI, prob=0.89)summary <- aggregate(flew_delta~mass_diff*flew_count, data=egg_data_final, FUN=mean)
#par(mar=c(5.5, 9.5, 4.5, 9.5), xpd=TRUE) # Add extra space to right of plot area; change clipping to figure
plot(egg_data_final$flew_delta~egg_data_final$mass_diff,
col=c(3,4,1)[as.factor(egg_data_final$flew_count)], # Green = Flew 0, Blue = Flew 1, Black = Flew 2,
main="Observed proportions of yes flew by mass differences between trials",
xlab = "Mass Differences between T1 and T2 (g)",
ylab= "Proportion Flew",
ylim=c(-1,0.2),
pch=c(16,1)[as.factor(egg_data_final$egg_diff)] # filled in means laid eggs in T1 but not T2 (not laying eggs in T2),
)
abline(coeff[1], coeff[2], col="black")
mtext(eq, side=4, cex=0.9)
shade(mu.PI, mass_diff.seq)
abline(v=0, col="mediumpurple3", lty=2)
text(0.01, 0.15, c("gained mass"),col="mediumpurple4")
text(-0.009, 0.155, c("lost mass"), col="mediumpurple4")
abline(h=0, col="darkolivegreen4", lty=2)
text(0.045, 0.087, c("flew in T2 but not T1"),col="darkolivegreen4")
text(0.045, -0.087, c("flew in T1 but not T2"),col="darkolivegreen4")
text(0.02, -0.5, p, col= "black", cex=0.9) #pval
legend(-0.04,-0.4,
legend = c("2", "1", "0"),
col= c(1,4,3),
pch = c(1),
title="Flew Count")
# for(r in 1:nrow(egg_data_final)){
# delta_egg <- egg_data_final$egg_diff[r]
# if(delta_egg < 1) {
# lwd = 3
# points(egg_data_final$mass_diff[r], egg_data_final$flew_delta[r], lwd=lwd, col=c(3,4,1)[as.factor(egg_data_final$flew_count[r])])
# #pch=16,
# #col=c(3,4,1)[as.factor(egg_data_final$flew_count[r])])
# }
# }Filled in circles mean that the bug laid eggs in T1, but not T2 so can show that really what matters more is the mass difference than the egg difference.
Now let’s do this for everyone.