Eclosion analyses for Gmm incubator experiments
DB
2025-08-29
Does temperature impact adult eclosion?
YES – SIG DIF **
simple way first. Drop with less than 5 pupae. get mean eclosion time
for each cage. compare to all other cages by genotype and by
temperature.
| temperature | mean_eclosion | sd | SE |
|---|---|---|---|
| 25 | 29.02680 | 0.8148064 | 0.1698989 |
| 27 | 24.32199 | 1.1138832 | 0.1180714 |
| 28.5 | 22.00108 | 1.0906557 | 0.3148452 |
| .y. | n | statistic | df | p | method |
|---|---|---|---|---|---|
| mean_eclosion_per_cage | 124 | 71.2685 | 2 | 0 | Kruskal-Wallis |
| .y. | group1 | group2 | n1 | n2 | statistic | p | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|
| mean_eclosion_per_cage | 25 | 27 | 23 | 89 | -6.761411 | 0.00e+00 | 0.00e+00 | **** |
| mean_eclosion_per_cage | 25 | 28.5 | 23 | 12 | -7.840872 | 0.00e+00 | 0.00e+00 | **** |
| mean_eclosion_per_cage | 27 | 28.5 | 89 | 12 | -3.936692 | 8.26e-05 | 8.26e-05 | **** |
more complex models. mixed linear models
model <- lmer(days_from_pupate_to_eclose ~ temperature + (1 | name), data = mod_dat)
days needed to elcose is the response with temperature as the predictor. dataframe has 1 row per pupae. grouping variable is the eclosion cup. maybe try this with a beta distribution and not gaussian?
should also try??: model is perhaps too complex…
glmmTMB(totalEclosedAdults ~ temperature * poly(days_from_pupate_to_eclose, 2) + (1 | name) + offset(log(numpupae)),
here we would have one row per eclosion cup. total elcosed adults is the response variable with temperature and days from pupate to eclose as the predictors. use cage as the grouping variable DO I drop smaller cages still? or do I keep them? keep them here..
## prepare vars
#mod_dat = long_eclosion_time
mod_dat = eclosion_dat %>% uncount(total_eclosed_adults)
str(mod_dat)## tibble [2,764 × 17] (S3: tbl_df/tbl/data.frame)
## $ treatment : chr [1:2764] "Inc28.5" "Inc28.5" "Inc28.5" "Inc28.5" ...
## $ Gono : chr [1:2764] "F1" "F1" "F1" "F1" ...
## $ name : chr [1:2764] "28.5_F1_8/4/23-8/11/23" "28.5_F1_8/4/23-8/11/23" "28.5_F1_8/4/23-8/11/23" "28.5_F1_8/4/23-8/11/23" ...
## $ genotype : chr [1:2764] "28.5_F1" "28.5_F1" "28.5_F1" "28.5_F1" ...
## $ generation : chr [1:2764] "F1" "F1" "F1" "F1" ...
## $ trayName : chr [1:2764] "blueG2F1" "blueG2F1" "blueG2F1" "blueG2F1" ...
## $ dateEclose : Date[1:2764], format: "2023-08-22" "2023-08-22" ...
## $ t0_start : Date[1:2764], format: "2023-08-04" "2023-08-04" ...
## $ t0_end : Date[1:2764], format: "2023-08-11" "2023-08-11" ...
## $ datePupate : chr [1:2764] "8/4/23-8/11/23" "8/4/23-8/11/23" "8/4/23-8/11/23" "8/4/23-8/11/23" ...
## $ numPupae : int [1:2764] 50 50 50 50 50 50 50 50 50 50 ...
## $ numFemaleEclose : int [1:2764] 4 4 4 4 0 0 7 7 7 7 ...
## $ numMaleEclose : int [1:2764] 0 0 0 0 1 1 4 4 4 4 ...
## $ mean.pupate.date : Date[1:2764], format: "2023-08-07" "2023-08-07" ...
## $ temperature : chr [1:2764] "28.5" "28.5" "28.5" "28.5" ...
## $ days_from_pupate_to_eclose: num [1:2764] 14.5 14.5 14.5 14.5 16.5 17.5 20.5 20.5 20.5 20.5 ...
## $ generation_numeric : chr [1:2764] "1" "1" "1" "1" ...mod_dat$genotype = as.factor(mod_dat$genotype)
mod_dat$temperature = as.factor(mod_dat$temperature)
mod_dat$generation_numeric = as.numeric(mod_dat$generation_numeric)
mod_dat$generation = as.factor(mod_dat$generation_numeric)
mod_dat$name = as.factor(mod_dat$name)
mod_dat$trayName = as.factor(mod_dat$trayName)
## lmer model
ggplot(mod_dat, aes(x = temperature, y = days_from_pupate_to_eclose, fill = temperature)) +
geom_boxplot() +
labs(title = "Eclosion Time by Rearing Temperature",
x = "Temperature (°C)",
y = "Days from Pupation to Eclosion") +
theme_minimal()
model <- lmer(days_from_pupate_to_eclose ~ temperature + (1 | name), data = mod_dat)
# Get the summary of the model, including p-values
summary(model)## Linear mixed model fit by REML ['lmerMod']
## Formula: days_from_pupate_to_eclose ~ temperature + (1 | name)
## Data: mod_dat
##
## REML criterion at convergence: 11422.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1191 -0.6028 -0.0410 0.5425 3.5100
##
## Random effects:
## Groups Name Variance Std.Dev.
## name (Intercept) 0.888 0.9423
## Residual 3.368 1.8351
## Number of obs: 2764, groups: name, 155
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 29.0028 0.2139 135.58
## temperature27 -4.6952 0.2392 -19.63
## temperature28.5 -7.0199 0.3743 -18.76
##
## Correlation of Fixed Effects:
## (Intr) tmpr27
## temperatr27 -0.894
## temprtr28.5 -0.572 0.511# Simulate residuals
simulationOutput <- simulateResiduals(fittedModel = model)
# Plot the residuals
plot(simulationOutput)
emmeans_temp <- emmeans(model, pairwise ~ temperature)
# Print the results
print(emmeans_temp)## $emmeans
## temperature emmean SE df lower.CL upper.CL
## 25 29.0 0.214 131 28.6 29.4
## 27 24.3 0.107 141 24.1 24.5
## 28.5 22.0 0.307 122 21.4 22.6
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## temperature25 - temperature27 4.70 0.239 133 19.618 <.0001
## temperature25 - temperature28.5 7.02 0.374 124 18.749 <.0001
## temperature27 - temperature28.5 2.32 0.325 123 7.145 <.0001
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesis there evidence for adaptation in the time it takes incubator flies to eclose?
contrast 27 incubator flies.
NO, NOT SIG
incEclose_time = long_eclosion_time
## select only 27 flies
incEclose_time27 = incEclose_time %>% filter(treatment == "Inc27")
# prep vars
incEclose_time27$genotype = as.factor(incEclose_time27$genotype)
incEclose_time27$temperature = as.factor(incEclose_time27$temperature)
incEclose_time27$generation_numeric = as.numeric(incEclose_time27$generation_numeric)
incEclose_time27$generation = as.factor(incEclose_time27$generation_numeric)
str(incEclose_time27)## tibble [2,000 × 17] (S3: tbl_df/tbl/data.frame)
## $ treatment : chr [1:2000] "Inc27" "Inc27" "Inc27" "Inc27" ...
## $ Gono : chr [1:2000] "" "" "" "" ...
## $ name : chr [1:2000] "27_F1_12/26/23" "27_F1_12/26/23" "27_F1_12/26/23" "27_F1_12/26/23" ...
## $ genotype : Factor w/ 4 levels "27_F1","27_F2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ generation : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ trayName : chr [1:2000] "24G2-1" "24G2-1" "24G2-1" "24G2-1" ...
## $ dateEclose : Date[1:2000], format: "2024-01-17" "2024-01-17" ...
## $ t0_start : Date[1:2000], format: "2023-12-26" "2023-12-26" ...
## $ t0_end : Date[1:2000], format: "2023-12-26" "2023-12-26" ...
## $ datePupate : chr [1:2000] "12/26/23" "12/26/23" "12/26/23" "12/26/23" ...
## $ numPupae : int [1:2000] 20 20 20 20 20 20 20 20 20 20 ...
## $ numFemaleEclose : int [1:2000] 13 13 13 13 13 13 13 13 13 13 ...
## $ numMaleEclose : int [1:2000] 0 0 0 0 0 0 0 0 0 0 ...
## $ mean.pupate.date : Date[1:2000], format: "2023-12-26" "2023-12-26" ...
## $ temperature : Factor w/ 1 level "27": 1 1 1 1 1 1 1 1 1 1 ...
## $ days_from_pupate_to_eclose: num [1:2000] 22 22 22 22 22 22 22 22 22 22 ...
## $ generation_numeric : num [1:2000] 1 1 1 1 1 1 1 1 1 1 ...sum_daysEclose27 = incEclose_time27 %>% group_by(name) %>%
summarise(mean_eclosion_per_cage = mean(days_from_pupate_to_eclose), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature))
eclose_days_table_gen27 = sum_daysEclose27 %>% group_by(genotype) %>%
summarise(mean_eclosion = mean(mean_eclosion_per_cage), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature),
sd = sd(mean_eclosion_per_cage), SE = sd/sqrt(length(mean_eclosion_per_cage)))
eclose_days_table_gen27## # A tibble: 4 × 6
## genotype mean_eclosion generation temperature sd SE
## <fct> <dbl> <fct> <fct> <dbl> <dbl>
## 1 27_F1 24.7 1 27 1.21 0.199
## 2 27_F2 24.1 2 27 1.11 0.207
## 3 27_F3 24.0 3 27 0.826 0.180
## 4 27_F4 24.3 4 27 0.990 0.700kw_genotype = sum_daysEclose27 %>% kruskal_test(mean_eclosion_per_cage ~ genotype)
dunn_genotype = sum_daysEclose27 %>% dunn_test(mean_eclosion_per_cage ~ genotype, p.adjust.method = "BH")
kable(kw_genotype, caption = "KW test eclosion by genotype")| .y. | n | statistic | df | p | method |
|---|---|---|---|---|---|
| mean_eclosion_per_cage | 89 | 6.569586 | 3 | 0.087 | Kruskal-Wallis |
sum_daysEclose27 %>%
ggplot(aes(y = mean_eclosion_per_cage, x = genotype, color = genotype)) +
stat_summary(fun = mean, geom = "point", position = position_dodge(0.5)) +
stat_summary(fun.max = function(x) mean(x) + (sd(x) / sqrt(length(x))),
fun.min = function(x) mean(x) - (sd(x) / sqrt(length(x))),
geom = "errorbar", width = 0.2, position = position_dodge(0.5), linewidth = 1) +
theme_set(theme_bw(base_size = 12)) +
theme(legend.position = "top") +
labs(x = "treatment", y = "mean time to eclosion per cage") +
ggtitle("mean time to eclosion by genotype") +
scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) +
scale_fill_viridis(discrete = TRUE) +
coord_cartesian(clip = 'off') +
geom_jitter(height = 0, width = .2, alpha = .4) + theme_bw()
#### more complex model
incEclose_time = eclosion_dat %>% uncount(total_eclosed_adults)
## select only 27 flies
incEclose_time27 = incEclose_time %>% filter(treatment == "Inc27")
# prep vars
incEclose_time27$genotype = as.factor(incEclose_time27$genotype)
incEclose_time27$temperature = as.factor(incEclose_time27$temperature)
incEclose_time27$generation_numeric = as.numeric(incEclose_time27$generation_numeric)
incEclose_time27$generation = as.factor(incEclose_time27$generation_numeric)
model <- lmer(days_from_pupate_to_eclose ~ generation + (1 | trayName), data = incEclose_time27)
# Get the summary of the model, including p-values
summary(model)## Linear mixed model fit by REML ['lmerMod']
## Formula: days_from_pupate_to_eclose ~ generation + (1 | trayName)
## Data: incEclose_time27
##
## REML criterion at convergence: 8429.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4439 -0.6121 -0.0106 0.5437 3.2798
##
## Random effects:
## Groups Name Variance Std.Dev.
## trayName (Intercept) 1.003 1.001
## Residual 3.255 1.804
## Number of obs: 2055, groups: trayName, 117
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 24.5661 0.1783 137.768
## generation2 -0.2907 0.2639 -1.102
## generation3 -0.5058 0.2920 -1.732
## generation4 -1.0028 0.5490 -1.826
##
## Correlation of Fixed Effects:
## (Intr) gnrtn2 gnrtn3
## generation2 -0.676
## generation3 -0.611 0.413
## generation4 -0.325 0.219 0.198# Simulate residuals
simulationOutput <- simulateResiduals(fittedModel = model)
# Plot the residuals
plot(simulationOutput)
emmeans_temp <- emmeans(model, pairwise ~ generation)
# Print the results
print(emmeans_temp)## $emmeans
## generation emmean SE df lower.CL upper.CL
## 1 24.6 0.178 77.3 24.2 24.9
## 2 24.3 0.195 102.1 23.9 24.7
## 3 24.1 0.231 147.6 23.6 24.5
## 4 23.6 0.520 261.5 22.5 24.6
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## generation1 - generation2 0.291 0.264 89.5 1.101 0.6900
## generation1 - generation3 0.506 0.292 113.3 1.731 0.3127
## generation1 - generation4 1.003 0.549 221.6 1.825 0.2643
## generation2 - generation3 0.215 0.303 125.7 0.711 0.8926
## generation2 - generation4 0.712 0.555 227.4 1.283 0.5747
## generation3 - generation4 0.497 0.569 235.4 0.874 0.8185
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimateswhat about for 28 treatment?
NO, NO SIG DIFF between generations
incEclose_time = long_eclosion_time
## select only 28.5 flies
incEclose_time28.5 = incEclose_time %>% filter(treatment == "Inc28.5")
# prep vars
incEclose_time28.5$genotype = as.factor(incEclose_time28.5$genotype)
incEclose_time28.5$temperature = as.factor(incEclose_time28.5$temperature)
incEclose_time28.5$generation_numeric = as.numeric(incEclose_time28.5$generation_numeric)
incEclose_time28.5$generation = as.factor(incEclose_time28.5$generation_numeric)
str(incEclose_time28.5)## tibble [294 × 17] (S3: tbl_df/tbl/data.frame)
## $ treatment : chr [1:294] "Inc28.5" "Inc28.5" "Inc28.5" "Inc28.5" ...
## $ Gono : chr [1:294] "F1" "F1" "F1" "F1" ...
## $ name : chr [1:294] "28.5_F1_8/4/23-8/11/23" "28.5_F1_8/4/23-8/11/23" "28.5_F1_8/4/23-8/11/23" "28.5_F1_8/4/23-8/11/23" ...
## $ genotype : Factor w/ 2 levels "28.5_F1","28.5_F2": 1 1 1 1 1 1 1 1 1 1 ...
## $ generation : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ trayName : chr [1:294] "blueG2F1" "blueG2F1" "blueG2F1" "blueG2F1" ...
## $ dateEclose : Date[1:294], format: "2023-08-22" "2023-08-22" ...
## $ t0_start : Date[1:294], format: "2023-08-04" "2023-08-04" ...
## $ t0_end : Date[1:294], format: "2023-08-11" "2023-08-11" ...
## $ datePupate : chr [1:294] "8/4/23-8/11/23" "8/4/23-8/11/23" "8/4/23-8/11/23" "8/4/23-8/11/23" ...
## $ numPupae : int [1:294] 50 50 50 50 50 50 50 50 50 50 ...
## $ numFemaleEclose : int [1:294] 4 4 4 4 0 0 7 7 7 7 ...
## $ numMaleEclose : int [1:294] 0 0 0 0 1 1 4 4 4 4 ...
## $ mean.pupate.date : Date[1:294], format: "2023-08-07" "2023-08-07" ...
## $ temperature : Factor w/ 1 level "28.5": 1 1 1 1 1 1 1 1 1 1 ...
## $ days_from_pupate_to_eclose: num [1:294] 14.5 14.5 14.5 14.5 16.5 17.5 20.5 20.5 20.5 20.5 ...
## $ generation_numeric : num [1:294] 1 1 1 1 1 1 1 1 1 1 ...sum_daysEclose28 = incEclose_time28.5 %>% group_by(name) %>%
summarise(mean_eclosion_per_cage = mean(days_from_pupate_to_eclose), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature))
eclose_days_table_gen28 = sum_daysEclose28 %>% group_by(genotype) %>%
summarise(mean_eclosion = mean(mean_eclosion_per_cage), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature),
sd = sd(mean_eclosion_per_cage), SE = sd/sqrt(length(mean_eclosion_per_cage)))
eclose_days_table_gen28## # A tibble: 2 × 6
## genotype mean_eclosion generation temperature sd SE
## <fct> <dbl> <fct> <fct> <dbl> <dbl>
## 1 28.5_F1 21.6 1 28.5 1.07 0.436
## 2 28.5_F2 22.4 2 28.5 1.07 0.439kw_genotype = sum_daysEclose28 %>% kruskal_test(mean_eclosion_per_cage ~ genotype)
dunn_genotype = sum_daysEclose28 %>% dunn_test(mean_eclosion_per_cage ~ genotype, p.adjust.method = "BH")
kable(kw_genotype, caption = "KW test eclosion by genotype")| .y. | n | statistic | df | p | method |
|---|---|---|---|---|---|
| mean_eclosion_per_cage | 12 | 0.6410256 | 1 | 0.423 | Kruskal-Wallis |
sum_daysEclose28 %>%
ggplot(aes(y = mean_eclosion_per_cage, x = genotype, color = genotype)) +
stat_summary(fun = mean, geom = "point", position = position_dodge(0.5)) +
stat_summary(fun.max = function(x) mean(x) + (sd(x) / sqrt(length(x))),
fun.min = function(x) mean(x) - (sd(x) / sqrt(length(x))),
geom = "errorbar", width = 0.2, position = position_dodge(0.5), linewidth = 1) +
theme_set(theme_bw(base_size = 12)) +
theme(legend.position = "top") +
labs(x = "treatment", y = "mean time to eclosion per cage") +
ggtitle("mean time to eclosion by genotype") +
scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) +
scale_fill_viridis(discrete = TRUE) +
coord_cartesian(clip = 'off') +
geom_jitter(height = 0, width = .2, alpha = .4) + theme_bw()
#### more complex model
## select only 28.5 flies
incEclose_time = eclosion_dat %>% uncount(total_eclosed_adults)
incEclose_time28.5 = incEclose_time %>% filter(treatment == "Inc28.5")
# prep vars
incEclose_time28.5$genotype = as.factor(incEclose_time28.5$genotype)
incEclose_time28.5$temperature = as.factor(incEclose_time28.5$temperature)
incEclose_time28.5$generation_numeric = as.numeric(incEclose_time28.5$generation_numeric)
incEclose_time28.5$generation = as.factor(incEclose_time28.5$generation_numeric)
model <- lmer(days_from_pupate_to_eclose ~ generation + (1 | trayName), data = incEclose_time28.5)
# Get the summary of the model, including p-values
summary(model)## Linear mixed model fit by REML ['lmerMod']
## Formula: days_from_pupate_to_eclose ~ generation + (1 | trayName)
## Data: incEclose_time28.5
##
## REML criterion at convergence: 1334.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2846 -0.6112 -0.1757 0.6415 2.8302
##
## Random effects:
## Groups Name Variance Std.Dev.
## trayName (Intercept) 0.703 0.8385
## Residual 5.242 2.2895
## Number of obs: 294, groups: trayName, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 21.7421 0.3884 55.982
## generation2 0.5677 0.6024 0.942
##
## Correlation of Fixed Effects:
## (Intr)
## generation2 -0.645# Simulate residuals
simulationOutput <- simulateResiduals(fittedModel = model)
# Plot the residuals
plot(simulationOutput)
emmeans_temp <- emmeans(model, pairwise ~ generation)
# Print the results
print(emmeans_temp)## $emmeans
## generation emmean SE df lower.CL upper.CL
## 1 21.7 0.391 7.26 20.8 22.7
## 2 22.3 0.463 13.93 21.3 23.3
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## generation1 - generation2 -0.568 0.606 10.4 -0.937 0.3700
##
## Degrees-of-freedom method: kenward-rogerdoes temperature impact adult eclosion success?
YES, SIG diff
eclose_rate = eclosion_dat %>% filter(numPupae >= 5)
eclose_rate_bycage = eclose_rate %>% group_by(name) %>%
summarise(treatment = unique(treatment), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature), trayName = unique(trayName), generation_numeric = unique(generation_numeric),
pupaeInCage = unique(numPupae), total_female_eclosed = sum(numFemaleEclose), total_male_eclosed = sum(numMaleEclose))
eclose_rate_bycage = eclose_rate_bycage %>% mutate(total_adult_eclosed = total_female_eclosed + total_male_eclosed) %>%
mutate(total_NOT_eclosed_adults = pupaeInCage - total_adult_eclosed) %>%
mutate(percent_adult_eclosed = total_adult_eclosed / pupaeInCage)
## summary tables
eclose_rate_table_gen = eclose_rate_bycage %>% group_by(genotype) %>%
summarise(mean_adult_eclosion_rate = mean(percent_adult_eclosed), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature),
sd = sd(percent_adult_eclosed), SE = sd/sqrt(length(percent_adult_eclosed)))
eclose_rate_table_temp = eclose_rate_bycage %>% group_by(temperature) %>%
summarise(mean_adult_eclosion_rate = mean(percent_adult_eclosed), temperature = unique(temperature),
sd = sd(percent_adult_eclosed), SE = sd/sqrt(length(percent_adult_eclosed)))
## simple stats
kw_temp= eclose_rate_bycage %>% kruskal_test(percent_adult_eclosed ~ temperature)
dunn_temp = eclose_rate_bycage %>% dunn_test(percent_adult_eclosed ~ temperature, p.adjust.method = "BH")
kable(kw_temp, caption = "KW test adult eclosion rate by temperature")| .y. | n | statistic | df | p | method |
|---|---|---|---|---|---|
| percent_adult_eclosed | 124 | 29.50406 | 2 | 4e-07 | Kruskal-Wallis |
kable(dunn_temp, caption = "Dunn test adult eclosion rate by temperature")| .y. | group1 | group2 | n1 | n2 | statistic | p | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|
| percent_adult_eclosed | 25 | 27 | 23 | 89 | -5.4178931 | 0.0000001 | 0.0000002 | **** |
| percent_adult_eclosed | 25 | 28.5 | 23 | 12 | -3.1588598 | 0.0015839 | 0.0023758 | ** |
| percent_adult_eclosed | 27 | 28.5 | 89 | 12 | 0.4631023 | 0.6432911 | 0.6432911 | ns |
### graph
eclose_rate_bycage %>%
ggplot(aes(y = percent_adult_eclosed, x = temperature, color = temperature)) +
stat_summary(fun = mean, geom = "point", position = position_dodge(0.5)) +
stat_summary(fun.max = function(x) mean(x) + (sd(x) / sqrt(length(x))),
fun.min = function(x) mean(x) - (sd(x) / sqrt(length(x))),
geom = "errorbar", width = 0.2, position = position_dodge(0.5), linewidth = 1) +
theme_set(theme_bw(base_size = 12)) +
theme(legend.position = "top") +
labs(x = "treatment", y = "mean eclosion ratet per cage") +
ggtitle("mean eclosion rate by temperature treatment") +
scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) +
scale_fill_viridis(discrete = TRUE) +
coord_cartesian(clip = 'off') +
geom_jitter(height = 0, width = .2, alpha = .4) + theme_bw()
more complex model
binary glm
mod_dat = eclosion_dat
mod_dat_bycage = mod_dat %>% group_by(name) %>%
summarise(treatment = unique(treatment), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature), trayName = unique(trayName), generation_numeric = unique(generation_numeric),
pupaeInCage = unique(numPupae), total_female_eclosed = sum(numFemaleEclose), total_male_eclosed = sum(numMaleEclose))
mod_dat_bycage = mod_dat_bycage %>% mutate(total_adult_eclosed = total_female_eclosed + total_male_eclosed) %>%
mutate(total_NOT_eclosed_adults = pupaeInCage - total_adult_eclosed) %>%
mutate(percent_adult_eclosed = total_adult_eclosed / pupaeInCage)
# Now, transform the data
mod_dat_bycage_long <- mod_dat_bycage %>%
# Step 1: Pivot the three count columns into a long format
pivot_longer(
cols = c(total_female_eclosed, total_male_eclosed, total_NOT_eclosed_adults),
names_to = "fly_type",
values_to = "count") %>%
# Step 2: Expand the rows based on the 'count' for each fly type
uncount(count) %>%
# Step 3: Create the 'eclosed' and 'sex' columns based on the 'fly_type'
mutate(
eclosed = if_else(fly_type == "total_NOT_eclosed_adults", 0, 1),
sex = case_when(
fly_type == "total_female_eclosed" ~ "female",
fly_type == "total_male_eclosed" ~ "male",
TRUE ~ NA_character_ # Assign NA for not eclosed
)
) %>%
# Step 4: Clean up by removing unnecessary columns
select(
-fly_type,
-total_adult_eclosed,
-percent_adult_eclosed
)
# prep vars
mod_dat_bycage_long$genotype = as.factor(mod_dat_bycage_long$genotype)
mod_dat_bycage_long$temperature = as.factor(mod_dat_bycage_long$temperature)
mod_dat_bycage_long$generation_numeric = as.numeric(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$generation = as.factor(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$name = as.factor(mod_dat_bycage_long$name)
mod_dat_bycage_long$sex = as.factor(mod_dat_bycage_long$sex)
str(mod_dat_bycage_long)## tibble [3,169 × 10] (S3: tbl_df/tbl/data.frame)
## $ name : Factor w/ 157 levels "25_wt_10/01/23",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ treatment : chr [1:3169] "Amb25" "Amb25" "Amb25" "Amb25" ...
## $ genotype : Factor w/ 7 levels "25_wt","27_F1",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ generation : Factor w/ 5 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ temperature : Factor w/ 3 levels "25","27","28.5": 1 1 1 1 1 1 1 1 1 1 ...
## $ trayName : chr [1:3169] "ColCage15" "ColCage15" "ColCage15" "ColCage15" ...
## $ generation_numeric: num [1:3169] 0 0 0 0 0 0 0 0 0 0 ...
## $ pupaeInCage : int [1:3169] 15 15 15 15 15 15 15 15 15 15 ...
## $ eclosed : num [1:3169] 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 2 2 ...model <- glmer(eclosed ~ temperature + (1|name),
data = mod_dat_bycage_long,
family = binomial(link = "logit"))
summary(model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: eclosed ~ temperature + (1 | name)
## Data: mod_dat_bycage_long
##
## AIC BIC logLik -2*log(L) df.resid
## 2271.1 2295.4 -1131.6 2263.1 3165
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5443 0.1674 0.3289 0.3833 1.1136
##
## Random effects:
## Groups Name Variance Std.Dev.
## name (Intercept) 0.6208 0.7879
## Number of obs: 3169, groups: name, 157
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.9531 0.3748 10.548 < 2e-16 ***
## temperature27 -2.0176 0.3864 -5.221 1.78e-07 ***
## temperature28.5 -2.0556 0.4787 -4.294 1.75e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tmpr27
## temperatr27 -0.957
## temprtr28.5 -0.766 0.738# Simulate residuals
simulationOutput <- simulateResiduals(fittedModel = model)
# Plot the residuals
plot(simulationOutput)
summary(model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: eclosed ~ temperature + (1 | name)
## Data: mod_dat_bycage_long
##
## AIC BIC logLik -2*log(L) df.resid
## 2271.1 2295.4 -1131.6 2263.1 3165
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5443 0.1674 0.3289 0.3833 1.1136
##
## Random effects:
## Groups Name Variance Std.Dev.
## name (Intercept) 0.6208 0.7879
## Number of obs: 3169, groups: name, 157
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.9531 0.3748 10.548 < 2e-16 ***
## temperature27 -2.0176 0.3864 -5.221 1.78e-07 ***
## temperature28.5 -2.0556 0.4787 -4.294 1.75e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tmpr27
## temperatr27 -0.957
## temprtr28.5 -0.766 0.738emmeans_temp <- emmeans(model, pairwise ~ temperature)
emmeans_temp## $emmeans
## temperature emmean SE df asymp.LCL asymp.UCL
## 25 3.95 0.375 Inf 3.22 4.69
## 27 1.94 0.112 Inf 1.72 2.15
## 28.5 1.90 0.308 Inf 1.29 2.50
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## temperature25 - temperature27 2.018 0.386 Inf 5.221 <.0001
## temperature25 - temperature28.5 2.056 0.479 Inf 4.294 0.0001
## temperature27 - temperature28.5 0.038 0.325 Inf 0.117 0.9925
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatesis there evidence for adaptation in the incubator treatments in eclosion rate?
NO, NO evidence for 27
eclose_rate = eclosion_dat %>% filter(numPupae >= 5) %>% filter(temperature == "27")
eclose_rate_bycage = eclose_rate %>% group_by(name) %>%
summarise(treatment = unique(treatment), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature), trayName = unique(trayName), generation_numeric = unique(generation_numeric),
pupaeInCage = unique(numPupae), total_female_eclosed = sum(numFemaleEclose), total_male_eclosed = sum(numMaleEclose))
eclose_rate_bycage = eclose_rate_bycage %>% mutate(total_adult_eclosed = total_female_eclosed + total_male_eclosed) %>%
mutate(total_NOT_eclosed_adults = pupaeInCage - total_adult_eclosed) %>%
mutate(percent_adult_eclosed = total_adult_eclosed / pupaeInCage)
## summary tables
eclose_rate_table_gen = eclose_rate_bycage %>% group_by(genotype) %>%
summarise(mean_adult_eclosion_rate = mean(percent_adult_eclosed), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature),
sd = sd(percent_adult_eclosed), SE = sd/sqrt(length(percent_adult_eclosed)))
eclose_rate_table_temp = eclose_rate_bycage %>% group_by(generation) %>%
summarise(mean_adult_eclosion_rate = mean(percent_adult_eclosed), generation = unique(generation),
sd = sd(percent_adult_eclosed), SE = sd/sqrt(length(percent_adult_eclosed)))
## simple stats
kw_temp= eclose_rate_bycage %>% kruskal_test(percent_adult_eclosed ~ generation)
dunn_temp = eclose_rate_bycage %>% dunn_test(percent_adult_eclosed ~ generation, p.adjust.method = "BH")
kable(kw_temp, caption = "KW test adult eclosion rate by temperature")| .y. | n | statistic | df | p | method |
|---|---|---|---|---|---|
| percent_adult_eclosed | 89 | 1.901981 | 3 | 0.593 | Kruskal-Wallis |
### graph
eclose_rate_bycage %>%
ggplot(aes(y = percent_adult_eclosed, x = generation, color = generation)) +
stat_summary(fun = mean, geom = "point", position = position_dodge(0.5)) +
stat_summary(fun.max = function(x) mean(x) + (sd(x) / sqrt(length(x))),
fun.min = function(x) mean(x) - (sd(x) / sqrt(length(x))),
geom = "errorbar", width = 0.2, position = position_dodge(0.5), linewidth = 1) +
theme_set(theme_bw(base_size = 12)) +
theme(legend.position = "top") +
labs(x = "treatment", y = "mean eclosion ratet per cage") +
ggtitle("mean eclosion rate by temperature treatment") +
scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) +
scale_fill_viridis(discrete = TRUE) +
coord_cartesian(clip = 'off') +
geom_jitter(height = 0, width = .2, alpha = .4) + theme_bw()
more complex model for 27
mod_dat = eclosion_dat %>% filter(temperature == "27")
mod_dat_bycage = mod_dat %>% group_by(name) %>%
summarise(treatment = unique(treatment), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature), trayName = unique(trayName), generation_numeric = unique(generation_numeric),
pupaeInCage = unique(numPupae), total_female_eclosed = sum(numFemaleEclose), total_male_eclosed = sum(numMaleEclose))
mod_dat_bycage = mod_dat_bycage %>% mutate(total_adult_eclosed = total_female_eclosed + total_male_eclosed) %>%
mutate(total_NOT_eclosed_adults = pupaeInCage - total_adult_eclosed) %>%
mutate(percent_adult_eclosed = total_adult_eclosed / pupaeInCage)
# Now, transform the data
mod_dat_bycage_long <- mod_dat_bycage %>%
# Step 1: Pivot the three count columns into a long format
pivot_longer(
cols = c(total_female_eclosed, total_male_eclosed, total_NOT_eclosed_adults),
names_to = "fly_type",
values_to = "count") %>%
# Step 2: Expand the rows based on the 'count' for each fly type
uncount(count) %>%
# Step 3: Create the 'eclosed' and 'sex' columns based on the 'fly_type'
mutate(
eclosed = if_else(fly_type == "total_NOT_eclosed_adults", 0, 1),
sex = case_when(
fly_type == "total_female_eclosed" ~ "female",
fly_type == "total_male_eclosed" ~ "male",
TRUE ~ NA_character_ # Assign NA for not eclosed
)
) %>%
# Step 4: Clean up by removing unnecessary columns
select(
-fly_type,
-total_adult_eclosed,
-percent_adult_eclosed
)
# prep vars
mod_dat_bycage_long$genotype = as.factor(mod_dat_bycage_long$genotype)
mod_dat_bycage_long$temperature = as.factor(mod_dat_bycage_long$temperature)
mod_dat_bycage_long$generation_numeric = as.numeric(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$generation = as.factor(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$name = as.factor(mod_dat_bycage_long$name)
mod_dat_bycage_long$sex = as.factor(mod_dat_bycage_long$sex)
str(mod_dat_bycage_long)## tibble [2,388 × 10] (S3: tbl_df/tbl/data.frame)
## $ name : Factor w/ 119 levels "27_F1_1/10/24",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ treatment : chr [1:2388] "Inc27" "Inc27" "Inc27" "Inc27" ...
## $ genotype : Factor w/ 4 levels "27_F1","27_F2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ generation : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ temperature : Factor w/ 1 level "27": 1 1 1 1 1 1 1 1 1 1 ...
## $ trayName : chr [1:2388] "24G2-7" "24G2-7" "24G2-7" "24G2-7" ...
## $ generation_numeric: num [1:2388] 1 1 1 1 1 1 1 1 1 1 ...
## $ pupaeInCage : int [1:2388] 38 38 38 38 38 38 38 38 38 38 ...
## $ eclosed : num [1:2388] 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 1 1 ...model <- glmer(eclosed ~ generation + (1|name),
data = mod_dat_bycage_long,
family = binomial(link = "logit"))
summary(model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: eclosed ~ generation + (1 | name)
## Data: mod_dat_bycage_long
##
## AIC BIC logLik -2*log(L) df.resid
## 1852.3 1881.1 -921.1 1842.3 2383
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7608 0.2938 0.3473 0.3875 1.1376
##
## Random effects:
## Groups Name Variance Std.Dev.
## name (Intercept) 0.6238 0.7898
## Number of obs: 2388, groups: name, 119
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.83762 0.15907 11.552 <2e-16 ***
## generation2 0.07602 0.24927 0.305 0.760
## generation3 0.35499 0.31324 1.133 0.257
## generation4 0.45482 0.71691 0.634 0.526
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) gnrtn2 gnrtn3
## generation2 -0.631
## generation3 -0.500 0.329
## generation4 -0.218 0.144 0.116# Simulate residuals
simulationOutput <- simulateResiduals(fittedModel = model)
# Plot the residuals
plot(simulationOutput)
emmeans_temp <- emmeans(model, pairwise ~ generation)
emmeans_temp## $emmeans
## generation emmean SE df asymp.LCL asymp.UCL
## 1 1.84 0.159 Inf 1.526 2.15
## 2 1.91 0.193 Inf 1.535 2.29
## 3 2.19 0.271 Inf 1.661 2.72
## 4 2.29 0.700 Inf 0.921 3.66
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## generation1 - generation2 -0.0760 0.249 Inf -0.305 0.9902
## generation1 - generation3 -0.3550 0.313 Inf -1.133 0.6689
## generation1 - generation4 -0.4548 0.717 Inf -0.634 0.9209
## generation2 - generation3 -0.2790 0.330 Inf -0.845 0.8327
## generation2 - generation4 -0.3788 0.724 Inf -0.523 0.9536
## generation3 - generation4 -0.0998 0.748 Inf -0.133 0.9992
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimateswhat about for eclosion rate for 28.5 C
NOPE not sig
eclose_rate = eclosion_dat %>% filter(numPupae >= 5) %>% filter(temperature == "28.5")
eclose_rate_bycage = eclose_rate %>% group_by(name) %>%
summarise(treatment = unique(treatment), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature), trayName = unique(trayName), generation_numeric = unique(generation_numeric),
pupaeInCage = unique(numPupae), total_female_eclosed = sum(numFemaleEclose), total_male_eclosed = sum(numMaleEclose))
eclose_rate_bycage = eclose_rate_bycage %>% mutate(total_adult_eclosed = total_female_eclosed + total_male_eclosed) %>%
mutate(total_NOT_eclosed_adults = pupaeInCage - total_adult_eclosed) %>%
mutate(percent_adult_eclosed = total_adult_eclosed / pupaeInCage)
## summary tables
eclose_rate_table_gen = eclose_rate_bycage %>% group_by(genotype) %>%
summarise(mean_adult_eclosion_rate = mean(percent_adult_eclosed), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature),
sd = sd(percent_adult_eclosed), SE = sd/sqrt(length(percent_adult_eclosed)))
eclose_rate_table_temp = eclose_rate_bycage %>% group_by(generation) %>%
summarise(mean_adult_eclosion_rate = mean(percent_adult_eclosed), generation = unique(generation),
sd = sd(percent_adult_eclosed), SE = sd/sqrt(length(percent_adult_eclosed)))
## simple stats
kw_temp= eclose_rate_bycage %>% kruskal_test(percent_adult_eclosed ~ generation)
dunn_temp = eclose_rate_bycage %>% dunn_test(percent_adult_eclosed ~ generation, p.adjust.method = "BH")
kable(kw_temp, caption = "KW test adult eclosion rate by temperature")| .y. | n | statistic | df | p | method |
|---|---|---|---|---|---|
| percent_adult_eclosed | 12 | 0.6524318 | 1 | 0.419 | Kruskal-Wallis |
### graph
eclose_rate_bycage %>%
ggplot(aes(y = percent_adult_eclosed, x = generation, color = generation)) +
stat_summary(fun = mean, geom = "point", position = position_dodge(0.5)) +
stat_summary(fun.max = function(x) mean(x) + (sd(x) / sqrt(length(x))),
fun.min = function(x) mean(x) - (sd(x) / sqrt(length(x))),
geom = "errorbar", width = 0.2, position = position_dodge(0.5), linewidth = 1) +
theme_set(theme_bw(base_size = 12)) +
theme(legend.position = "top") +
labs(x = "treatment", y = "mean eclosion ratet per cage") +
ggtitle("mean eclosion rate by temperature treatment") +
scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) +
scale_fill_viridis(discrete = TRUE) +
coord_cartesian(clip = 'off') +
geom_jitter(height = 0, width = .2, alpha = .4) + theme_bw()
more complex model for 28.5
nope not sig
mod_dat = eclosion_dat %>% filter(temperature == "28.5")
mod_dat_bycage = mod_dat %>% group_by(name) %>%
summarise(treatment = unique(treatment), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature), trayName = unique(trayName), generation_numeric = unique(generation_numeric),
pupaeInCage = unique(numPupae), total_female_eclosed = sum(numFemaleEclose), total_male_eclosed = sum(numMaleEclose))
mod_dat_bycage = mod_dat_bycage %>% mutate(total_adult_eclosed = total_female_eclosed + total_male_eclosed) %>%
mutate(total_NOT_eclosed_adults = pupaeInCage - total_adult_eclosed) %>%
mutate(percent_adult_eclosed = total_adult_eclosed / pupaeInCage)
# Now, transform the data
mod_dat_bycage_long <- mod_dat_bycage %>%
pivot_longer(
cols = c(total_female_eclosed, total_male_eclosed, total_NOT_eclosed_adults),
names_to = "fly_type", values_to = "count") %>%
uncount(count) %>%
mutate(eclosed = if_else(fly_type == "total_NOT_eclosed_adults", 0, 1),
sex = case_when(
fly_type == "total_female_eclosed" ~ "female",
fly_type == "total_male_eclosed" ~ "male",
TRUE ~ NA_character_ # Assign NA for not eclosed
)
) %>%select(
-fly_type,
-total_adult_eclosed,
-percent_adult_eclosed
)
# prep vars
mod_dat_bycage_long$genotype = as.factor(mod_dat_bycage_long$genotype)
mod_dat_bycage_long$temperature = as.factor(mod_dat_bycage_long$temperature)
mod_dat_bycage_long$generation_numeric = as.numeric(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$generation = as.factor(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$name = as.factor(mod_dat_bycage_long$name)
mod_dat_bycage_long$sex = as.factor(mod_dat_bycage_long$sex)
str(mod_dat_bycage_long)## tibble [355 × 10] (S3: tbl_df/tbl/data.frame)
## $ name : Factor w/ 12 levels "28.5_F1_8/11/23-8/18/23",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ treatment : chr [1:355] "Inc28.5" "Inc28.5" "Inc28.5" "Inc28.5" ...
## $ genotype : Factor w/ 2 levels "28.5_F1","28.5_F2": 1 1 1 1 1 1 1 1 1 1 ...
## $ generation : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ temperature : Factor w/ 1 level "28.5": 1 1 1 1 1 1 1 1 1 1 ...
## $ trayName : chr [1:355] "pinkG2F1" "pinkG2F1" "pinkG2F1" "pinkG2F1" ...
## $ generation_numeric: num [1:355] 1 1 1 1 1 1 1 1 1 1 ...
## $ pupaeInCage : int [1:355] 116 116 116 116 116 116 116 116 116 116 ...
## $ eclosed : num [1:355] 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 1 1 ...model <- glmer(eclosed ~ generation + (1|name),
data = mod_dat_bycage_long,
family = binomial(link = "logit"))
summary(model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: eclosed ~ generation + (1 | name)
## Data: mod_dat_bycage_long
##
## AIC BIC logLik -2*log(L) df.resid
## 324.8 336.4 -159.4 318.8 352
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3220 0.3145 0.3393 0.5367 0.5819
##
## Random effects:
## Groups Name Variance Std.Dev.
## name (Intercept) 0.4131 0.6427
## Number of obs: 355, groups: name, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6810 0.3372 4.986 6.18e-07 ***
## generation2 0.4468 0.5735 0.779 0.436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## generation2 -0.530# Simulate residuals
simulationOutput <- simulateResiduals(fittedModel = model)
# Plot the residuals
plot(simulationOutput)
emmeans_temp <- emmeans(model, pairwise ~ generation)
emmeans_temp## $emmeans
## generation emmean SE df asymp.LCL asymp.UCL
## 1 1.68 0.337 Inf 1.02 2.34
## 2 2.13 0.488 Inf 1.17 3.08
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## generation1 - generation2 -0.447 0.573 Inf -0.779 0.4360
##
## Results are given on the log odds ratio (not the response) scale.is there a sex bias in eclosion at different temperatures?
NO
mod_dat = eclosion_dat
mod_dat_bycage = mod_dat %>% group_by(name) %>%
summarise(treatment = unique(treatment), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature), trayName = unique(trayName), generation_numeric = unique(generation_numeric),
pupaeInCage = unique(numPupae), total_female_eclosed = sum(numFemaleEclose), total_male_eclosed = sum(numMaleEclose))
mod_dat_bycage = mod_dat_bycage %>% mutate(total_adult_eclosed = total_female_eclosed + total_male_eclosed) %>%
mutate(total_NOT_eclosed_adults = pupaeInCage - total_adult_eclosed) %>%
mutate(percent_adult_eclosed = total_adult_eclosed / pupaeInCage)
# Now, transform the data
mod_dat_bycage_long <- mod_dat_bycage %>%
# Step 1: Pivot the three count columns into a long format
pivot_longer(
cols = c(total_female_eclosed, total_male_eclosed, total_NOT_eclosed_adults),
names_to = "fly_type",
values_to = "count") %>%
# Step 2: Expand the rows based on the 'count' for each fly type
uncount(count) %>%
# Step 3: Create the 'eclosed' and 'sex' columns based on the 'fly_type'
mutate(
eclosed = if_else(fly_type == "total_NOT_eclosed_adults", 0, 1),
sex = case_when(
fly_type == "total_female_eclosed" ~ "female",
fly_type == "total_male_eclosed" ~ "male",
TRUE ~ NA_character_ # Assign NA for not eclosed
)
) %>%
# Step 4: Clean up by removing unnecessary columns
select(
-fly_type,
-total_adult_eclosed,
-percent_adult_eclosed
)
mod_dat_bycage_long = mod_dat_bycage_long %>% filter(eclosed ==1)
# prep vars
mod_dat_bycage_long$genotype = as.factor(mod_dat_bycage_long$genotype)
mod_dat_bycage_long$temperature = as.factor(mod_dat_bycage_long$temperature)
mod_dat_bycage_long$generation_numeric = as.numeric(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$generation = as.factor(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$name = as.factor(mod_dat_bycage_long$name)
mod_dat_bycage_long$sex = as.factor(mod_dat_bycage_long$sex)
str(mod_dat_bycage_long)## tibble [2,764 × 10] (S3: tbl_df/tbl/data.frame)
## $ name : Factor w/ 155 levels "25_wt_10/01/23",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ treatment : chr [1:2764] "Amb25" "Amb25" "Amb25" "Amb25" ...
## $ genotype : Factor w/ 7 levels "25_wt","27_F1",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ generation : Factor w/ 5 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ temperature : Factor w/ 3 levels "25","27","28.5": 1 1 1 1 1 1 1 1 1 1 ...
## $ trayName : chr [1:2764] "ColCage15" "ColCage15" "ColCage15" "ColCage15" ...
## $ generation_numeric: num [1:2764] 0 0 0 0 0 0 0 0 0 0 ...
## $ pupaeInCage : int [1:2764] 15 15 15 15 15 15 15 15 15 15 ...
## $ eclosed : num [1:2764] 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 2 2 ...# We use family = binomial because our outcome (sex) has two levels.
sex_ratio_model <- glm(sex ~ temperature,
data = mod_dat_bycage_long,
family = binomial)
# Get the detailed summary of the model
summary(sex_ratio_model)##
## Call:
## glm(formula = sex ~ temperature, family = binomial, data = mod_dat_bycage_long)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.004819 0.098176 -0.049 0.961
## temperature27 0.025258 0.107635 0.235 0.814
## temperature28.5 -0.117783 0.152628 -0.772 0.440
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 3831.7 on 2763 degrees of freedom
## Residual deviance: 3830.4 on 2761 degrees of freedom
## AIC: 3836.4
##
## Number of Fisher Scoring iterations: 3# Simulate residuals
simulationOutput <- simulateResiduals(fittedModel = sex_ratio_model)
# Plot the residuals
plot(simulationOutput)
prob_male_25C <- plogis(coef(sex_ratio_model)["(Intercept)"])
print(paste("Predicted probability of being male at 25°C:", round(prob_male_25C, 3)))## [1] "Predicted probability of being male at 25°C: 0.499"prob_male_27C <- plogis(coef(sex_ratio_model)["(Intercept)"] + coef(sex_ratio_model)["temperature27"])
print(paste("Predicted probability of being male at 27°C:", round(prob_male_27C, 3)))## [1] "Predicted probability of being male at 27°C: 0.505"prob_male_28.5C <- plogis(coef(sex_ratio_model)["(Intercept)"] + coef(sex_ratio_model)["temperature28.5"])
print(paste("Predicted probability of being male at 28.5°C:", round(prob_male_28.5C, 3)))## [1] "Predicted probability of being male at 28.5°C: 0.469"### but model shows differences thus these are NOT significant. Did adaptation cause sex bias in 27 flies?
NO, not sig
mod_dat = eclosion_dat %>% filter(temperature == "27")
mod_dat_bycage = mod_dat %>% group_by(name) %>%
summarise(treatment = unique(treatment), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature), trayName = unique(trayName), generation_numeric = unique(generation_numeric),
pupaeInCage = unique(numPupae), total_female_eclosed = sum(numFemaleEclose), total_male_eclosed = sum(numMaleEclose))
mod_dat_bycage = mod_dat_bycage %>% mutate(total_adult_eclosed = total_female_eclosed + total_male_eclosed) %>%
mutate(total_NOT_eclosed_adults = pupaeInCage - total_adult_eclosed) %>%
mutate(percent_adult_eclosed = total_adult_eclosed / pupaeInCage)
# Now, transform the data
mod_dat_bycage_long <- mod_dat_bycage %>%
# Step 1: Pivot the three count columns into a long format
pivot_longer(
cols = c(total_female_eclosed, total_male_eclosed, total_NOT_eclosed_adults),
names_to = "fly_type",
values_to = "count") %>%
# Step 2: Expand the rows based on the 'count' for each fly type
uncount(count) %>%
# Step 3: Create the 'eclosed' and 'sex' columns based on the 'fly_type'
mutate(
eclosed = if_else(fly_type == "total_NOT_eclosed_adults", 0, 1),
sex = case_when(
fly_type == "total_female_eclosed" ~ "female",
fly_type == "total_male_eclosed" ~ "male",
TRUE ~ NA_character_ # Assign NA for not eclosed
)
) %>%
# Step 4: Clean up by removing unnecessary columns
select(
-fly_type,
-total_adult_eclosed,
-percent_adult_eclosed
)
mod_dat_bycage_long = mod_dat_bycage_long %>% filter(eclosed ==1)
# prep vars
mod_dat_bycage_long$genotype = as.factor(mod_dat_bycage_long$genotype)
mod_dat_bycage_long$temperature = as.factor(mod_dat_bycage_long$temperature)
mod_dat_bycage_long$generation_numeric = as.numeric(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$generation = as.factor(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$name = as.factor(mod_dat_bycage_long$name)
mod_dat_bycage_long$sex = as.factor(mod_dat_bycage_long$sex)
str(mod_dat_bycage_long)## tibble [2,055 × 10] (S3: tbl_df/tbl/data.frame)
## $ name : Factor w/ 117 levels "27_F1_1/10/24",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ treatment : chr [1:2055] "Inc27" "Inc27" "Inc27" "Inc27" ...
## $ genotype : Factor w/ 4 levels "27_F1","27_F2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ generation : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ temperature : Factor w/ 1 level "27": 1 1 1 1 1 1 1 1 1 1 ...
## $ trayName : chr [1:2055] "24G2-7" "24G2-7" "24G2-7" "24G2-7" ...
## $ generation_numeric: num [1:2055] 1 1 1 1 1 1 1 1 1 1 ...
## $ pupaeInCage : int [1:2055] 38 38 38 38 38 38 38 38 38 38 ...
## $ eclosed : num [1:2055] 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 1 1 ...# We use family = binomial because our outcome (sex) has two levels.
sex_ratio_model <- glm(sex ~ generation,
data = mod_dat_bycage_long,
family = binomial)
# Get the detailed summary of the model
summary(sex_ratio_model)##
## Call:
## glm(formula = sex ~ generation, family = binomial, data = mod_dat_bycage_long)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.007918 0.056277 0.141 0.888
## generation2 0.052277 0.101258 0.516 0.606
## generation3 -0.007918 0.151556 -0.052 0.958
## generation4 -0.087960 0.404257 -0.218 0.828
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2848.6 on 2054 degrees of freedom
## Residual deviance: 2848.3 on 2051 degrees of freedom
## AIC: 2856.3
##
## Number of Fisher Scoring iterations: 3# Simulate residuals
simulationOutput <- simulateResiduals(fittedModel = sex_ratio_model)
# Plot the residuals
plot(simulationOutput)
prob_F1 <- plogis(coef(sex_ratio_model)["(Intercept)"])
print(paste("Predicted probability of being male at 27°C F1:", round(prob_F1, 3)))## [1] "Predicted probability of being male at 27°C F1: 0.502"prob_F2 <- plogis(coef(sex_ratio_model)["(Intercept)"] + coef(sex_ratio_model)["generation2"])
print(paste("Predicted probability of being male at 27°C F2:", round(prob_F2, 3)))## [1] "Predicted probability of being male at 27°C F2: 0.515"prob_F3 <- plogis(coef(sex_ratio_model)["(Intercept)"] + coef(sex_ratio_model)["generation3"])
print(paste("Predicted probability of being male at 27°C F2:", round(prob_F3, 3)))## [1] "Predicted probability of being male at 27°C F2: 0.5"prob_F4 <- plogis(coef(sex_ratio_model)["(Intercept)"] + coef(sex_ratio_model)["generation4"])
print(paste("Predicted probability of being male at 27°C F2:", round(prob_F4, 3)))## [1] "Predicted probability of being male at 27°C F2: 0.48"### but model shows differences are NOT signifigant. Did adaptation cause sex bias in 28.5 flies?
YES, sig diff
mod_dat = eclosion_dat %>% filter(temperature == "28.5")
mod_dat_bycage = mod_dat %>% group_by(name) %>%
summarise(treatment = unique(treatment), genotype = unique(genotype), generation = unique(generation), temperature = unique(temperature), trayName = unique(trayName), generation_numeric = unique(generation_numeric),
pupaeInCage = unique(numPupae), total_female_eclosed = sum(numFemaleEclose), total_male_eclosed = sum(numMaleEclose))
mod_dat_bycage = mod_dat_bycage %>% mutate(total_adult_eclosed = total_female_eclosed + total_male_eclosed) %>%
mutate(total_NOT_eclosed_adults = pupaeInCage - total_adult_eclosed) %>%
mutate(percent_adult_eclosed = total_adult_eclosed / pupaeInCage)
# Now, transform the data
mod_dat_bycage_long <- mod_dat_bycage %>%
# Step 1: Pivot the three count columns into a long format
pivot_longer(
cols = c(total_female_eclosed, total_male_eclosed, total_NOT_eclosed_adults),
names_to = "fly_type",
values_to = "count") %>%
# Step 2: Expand the rows based on the 'count' for each fly type
uncount(count) %>%
# Step 3: Create the 'eclosed' and 'sex' columns based on the 'fly_type'
mutate(
eclosed = if_else(fly_type == "total_NOT_eclosed_adults", 0, 1),
sex = case_when(
fly_type == "total_female_eclosed" ~ "female",
fly_type == "total_male_eclosed" ~ "male",
TRUE ~ NA_character_ # Assign NA for not eclosed
)
) %>%
# Step 4: Clean up by removing unnecessary columns
select(
-fly_type,
-total_adult_eclosed,
-percent_adult_eclosed
)
mod_dat_bycage_long = mod_dat_bycage_long %>% filter(eclosed ==1)
# prep vars
mod_dat_bycage_long$genotype = as.factor(mod_dat_bycage_long$genotype)
mod_dat_bycage_long$temperature = as.factor(mod_dat_bycage_long$temperature)
mod_dat_bycage_long$generation_numeric = as.numeric(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$generation = as.factor(mod_dat_bycage_long$generation_numeric)
mod_dat_bycage_long$name = as.factor(mod_dat_bycage_long$name)
mod_dat_bycage_long$sex = as.factor(mod_dat_bycage_long$sex)
str(mod_dat_bycage_long)## tibble [294 × 10] (S3: tbl_df/tbl/data.frame)
## $ name : Factor w/ 12 levels "28.5_F1_8/11/23-8/18/23",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ treatment : chr [1:294] "Inc28.5" "Inc28.5" "Inc28.5" "Inc28.5" ...
## $ genotype : Factor w/ 2 levels "28.5_F1","28.5_F2": 1 1 1 1 1 1 1 1 1 1 ...
## $ generation : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ temperature : Factor w/ 1 level "28.5": 1 1 1 1 1 1 1 1 1 1 ...
## $ trayName : chr [1:294] "pinkG2F1" "pinkG2F1" "pinkG2F1" "pinkG2F1" ...
## $ generation_numeric: num [1:294] 1 1 1 1 1 1 1 1 1 1 ...
## $ pupaeInCage : int [1:294] 116 116 116 116 116 116 116 116 116 116 ...
## $ eclosed : num [1:294] 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 1 1 ...# We use family = binomial because our outcome (sex) has two levels.
sex_ratio_model <- glm(sex ~ generation,
data = mod_dat_bycage_long,
family = binomial)
# Get the detailed summary of the model
summary(sex_ratio_model)##
## Call:
## glm(formula = sex ~ generation, family = binomial, data = mod_dat_bycage_long)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2776 0.1326 -2.094 0.0362 *
## generation2 0.7372 0.2925 2.520 0.0117 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 406.47 on 293 degrees of freedom
## Residual deviance: 399.95 on 292 degrees of freedom
## AIC: 403.95
##
## Number of Fisher Scoring iterations: 4# Simulate residuals
simulationOutput <- simulateResiduals(fittedModel = sex_ratio_model)
# Plot the residuals
plot(simulationOutput)
prob_F1 <- plogis(coef(sex_ratio_model)["(Intercept)"])
print(paste("Predicted probability of being male at 28.5°C F1:", round(prob_F1, 3)))## [1] "Predicted probability of being male at 28.5°C F1: 0.431"prob_F2 <- plogis(coef(sex_ratio_model)["(Intercept)"] + coef(sex_ratio_model)["generation2"])
print(paste("Predicted probability of being male at 28.5°C F2:", round(prob_F2, 3)))## [1] "Predicted probability of being male at 28.5°C F2: 0.613"prob_F3 <- plogis(coef(sex_ratio_model)["(Intercept)"] + coef(sex_ratio_model)["generation3"])
print(paste("Predicted probability of being male at 27°C F2:", round(prob_F3, 3)))## [1] "Predicted probability of being male at 27°C F2: NA"prob_F4 <- plogis(coef(sex_ratio_model)["(Intercept)"] + coef(sex_ratio_model)["generation4"])
print(paste("Predicted probability of being male at 27°C F2:", round(prob_F4, 3)))## [1] "Predicted probability of being male at 27°C F2: NA"### but model shows differences are NOT signifigant.