Phenology ABM: Analyses for Revisions
Shannon Carter
July 10, 2020
2-Species Models
Load data
This has already been processed and cleaned in another script. This data has one row per individual and gives info on treatment parameters (11 levels of synchrony, 3 mean treatments), results of the simulation (number of survivors, biomass), and stats for the individual (hatching date, survival status, development time, final size).
## X run_num sync_fish sync_dfly mean_fish mean_dfly surv_fish surv_dfly
## 1 1 1 0 5 early 60 40 29
## 2 2 1 0 5 early 60 40 29
## 3 3 1 0 5 early 60 40 29
## 4 4 1 0 5 early 60 40 29
## 5 5 1 0 5 early 60 40 29
## 6 6 1 0 5 early 60 40 29
## biom_fish biom_dfly biom_tot sp id meta hatch_date end_date growth_time
## 1 287.3 192.7 480 d 1 1 63 109 46
## 2 287.3 192.7 480 d 10 1 64 126 62
## 3 287.3 192.7 480 d 11 1 70 158 88
## 4 287.3 192.7 480 d 12 1 66 158 92
## 5 287.3 192.7 480 d 13 1 66 130 64
## 6 287.3 192.7 480 d 14 1 70 148 78
## growth_rate max_size
## 1 0.13369565 6.15
## 2 0.11048387 6.85
## 3 0.07500000 6.60
## 4 0.07065217 6.50
## 5 0.09687500 6.20
## 6 0.08333333 6.50Survival model
Prepare data
Essentially, pool replicate treatments, filter for focal species, and calcuate standardized metrics for each treatment.
## Prepare data
survival_model_data <- ind %>%
# remove run_num to pool all replicates of a treatment. also remove others we don't need
select(-c(X, run_num, biom_fish, biom_dfly, biom_tot, end_date, growth_rate, id, sync_dfly, mean_dfly)) %>%
# "fish" is focal species that we manipulated
filter(sp == "f") %>%
# only a few of these... I think the simulation ended before they "hatched"
filter(!is.na(max_size)) %>%
# group by treatment to do the relativazation calculations
group_by(sync_fish, mean_fish) %>%
# standardize arrival date for each treatment combination
mutate(trt_max_hatch_date = max(hatch_date),
trt_min_hatch_date = min(hatch_date),
trt_mean_hatch_date = mean(hatch_date),
hatch_percentile = round((hatch_date - trt_mean_hatch_date)/(trt_max_hatch_date - trt_min_hatch_date), 3))
head(survival_model_data)## # A tibble: 6 x 13
## # Groups: sync_fish, mean_fish [1]
## sync_fish mean_fish surv_fish surv_dfly sp meta hatch_date
## <int> <fct> <int> <int> <fct> <int> <int>
## 1 0 early 40 29 f 1 51
## 2 0 early 40 29 f 1 52
## 3 0 early 40 29 f 1 52
## 4 0 early 40 29 f 1 52
## 5 0 early 40 29 f 1 52
## 6 0 early 40 29 f 1 52
## # … with 6 more variables: growth_time <int>, max_size <dbl>,
## # trt_max_hatch_date <int>, trt_min_hatch_date <int>,
## # trt_mean_hatch_date <dbl>, hatch_percentile <dbl>Model survival
Logistic regression model with survival as a binary outcome
## Fitness = survival model
m_surv <- glm(data = survival_model_data, meta ~ hatch_percentile + sync_fish + mean_fish +
hatch_percentile*sync_fish +
hatch_percentile*mean_fish +
sync_fish*mean_fish +
hatch_percentile*sync_fish*mean_fish,
family = binomial(link = "logit"))
summary(m_surv)##
## Call:
## glm(formula = meta ~ hatch_percentile + sync_fish + mean_fish +
## hatch_percentile * sync_fish + hatch_percentile * mean_fish +
## sync_fish * mean_fish + hatch_percentile * sync_fish * mean_fish,
## family = binomial(link = "logit"), data = survival_model_data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0793 -0.1754 0.0019 0.1594 3.2142
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 5.442283 0.527421 10.319
## hatch_percentile -2.818863 2.635637 -1.070
## sync_fish -0.012956 0.029140 -0.445
## mean_fishlate -7.282795 0.543728 -13.394
## mean_fishsame -4.311525 0.541471 -7.963
## hatch_percentile:sync_fish -2.472857 0.262704 -9.413
## hatch_percentile:mean_fishlate 1.365869 2.743686 0.498
## hatch_percentile:mean_fishsame 0.339938 2.935901 0.116
## sync_fish:mean_fishlate 0.031611 0.030563 1.034
## sync_fish:mean_fishsame 0.007751 0.030684 0.253
## hatch_percentile:sync_fish:mean_fishlate 0.857482 0.285706 3.001
## hatch_percentile:sync_fish:mean_fishsame 0.474967 0.316789 1.499
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## hatch_percentile 0.28484
## sync_fish 0.65659
## mean_fishlate < 2e-16 ***
## mean_fishsame 1.68e-15 ***
## hatch_percentile:sync_fish < 2e-16 ***
## hatch_percentile:mean_fishlate 0.61861
## hatch_percentile:mean_fishsame 0.90782
## sync_fish:mean_fishlate 0.30100
## sync_fish:mean_fishsame 0.80057
## hatch_percentile:sync_fish:mean_fishlate 0.00269 **
## hatch_percentile:sync_fish:mean_fishsame 0.13379
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 9012.9 on 6599 degrees of freedom
## Residual deviance: 2678.6 on 6588 degrees of freedom
## AIC: 2702.6
##
## Number of Fisher Scoring iterations: 9## (Intercept)
## 5.442282864
## hatch_percentile
## -2.818863209
## sync_fish
## -0.012956338
## mean_fishlate
## -7.282795294
## mean_fishsame
## -4.311524613
## hatch_percentile:sync_fish
## -2.472857174
## hatch_percentile:mean_fishlate
## 1.365868662
## hatch_percentile:mean_fishsame
## 0.339938417
## sync_fish:mean_fishlate
## 0.031610736
## sync_fish:mean_fishsame
## 0.007751155
## hatch_percentile:sync_fish:mean_fishlate
## 0.857481560
## hatch_percentile:sync_fish:mean_fishsame
## 0.474967443# make a data frame of individuals to add predicted survival to
predict_df <- data.frame(hatch_percentile = survival_model_data$hatch_percentile,
sync_fish = survival_model_data$sync_fish,
mean_fish = survival_model_data$mean_fish)
# calculate predicted survival + new value for each individual
m2 <- predict_df %>%
mutate(predicted_survival = predict(m_surv, type = "response")) %>%
mutate(coef = predicted_survival * (1 - predicted_survival))
# histogram of predicted survival
ggplot(m2, aes(x = predicted_survival)) + geom_histogram() +
facet_grid(sync_fish ~ mean_fish)
# calculate selection coefficients -- absolute and standardized
selection_coefficients <- m2 %>%
group_by(sync_fish, mean_fish) %>%
summarize(selection_coef = mean(coef), # absolute: 0.0628
standardized_selection_coef = selection_coef/mean(survival_model_data$meta)) # standardized: 0.110
selection_coefficients## # A tibble: 33 x 4
## # Groups: sync_fish [11]
## sync_fish mean_fish selection_coef standardized_selection_coef
## <int> <fct> <dbl> <dbl>
## 1 0 early 0.00449 0.00786
## 2 0 late 0.118 0.207
## 3 0 same 0.186 0.326
## 4 3 early 0.0167 0.0291
## 5 3 late 0.129 0.225
## 6 3 same 0.146 0.255
## 7 6 early 0.0420 0.0734
## 8 6 late 0.113 0.197
## 9 6 same 0.100 0.175
## 10 9 early 0.0460 0.0805
## # … with 23 more rowsMass model
Prepare data
Almost the same as before, except we also filter out individuals that did not survive, and calculate relative size in each treatment
## Prepare data
size_model_data <- ind %>%
# remove run_num to pool all replicates of a treatment. also remove others we don't need
select(-c(X, run_num, biom_fish, biom_dfly, biom_tot, end_date, growth_rate, id, sync_dfly, mean_dfly)) %>%
# "fish" is focal species that we manipulated
filter(sp == "f") %>%
# filter only those that survived
filter(meta == 1) %>%
# only a few of these... I think the simulation ended before they "hatched"
filter(!is.na(max_size)) %>%
# group by treatment to do the relativazation calculations
group_by(sync_fish, mean_fish) %>%
# standardize fitness by dividing mass by treatment mean
mutate(size_relative = max_size / mean(max_size)) %>%
# standardize arrival date for each treatment combination
mutate(trt_max_hatch_date = max(hatch_date),
trt_min_hatch_date = min(hatch_date),
trt_mean_hatch_date = mean(hatch_date),
hatch_percentile = round((hatch_date - trt_mean_hatch_date)/(trt_max_hatch_date - trt_min_hatch_date), 3))
head(size_model_data)## # A tibble: 6 x 14
## # Groups: sync_fish, mean_fish [1]
## sync_fish mean_fish surv_fish surv_dfly sp meta hatch_date
## <int> <fct> <int> <int> <fct> <int> <int>
## 1 0 early 40 29 f 1 51
## 2 0 early 40 29 f 1 52
## 3 0 early 40 29 f 1 52
## 4 0 early 40 29 f 1 52
## 5 0 early 40 29 f 1 52
## 6 0 early 40 29 f 1 52
## # … with 7 more variables: growth_time <int>, max_size <dbl>,
## # size_relative <dbl>, trt_max_hatch_date <int>,
## # trt_min_hatch_date <int>, trt_mean_hatch_date <dbl>,
## # hatch_percentile <dbl># View distribution of final sizes - pretty normal
ggplot(size_model_data, aes(x = size_relative)) + geom_histogram() + mytheme
Model mass
## Fitness = mass model
m_size <- glm(data = size_model_data,
size_relative ~
hatch_percentile + sync_fish + mean_fish +
hatch_percentile*sync_fish +
hatch_percentile*mean_fish +
sync_fish*mean_fish +
hatch_percentile*sync_fish*mean_fish)
summary(m_size)##
## Call:
## glm(formula = size_relative ~ hatch_percentile + sync_fish +
## mean_fish + hatch_percentile * sync_fish + hatch_percentile *
## mean_fish + sync_fish * mean_fish + hatch_percentile * sync_fish *
## mean_fish, data = size_model_data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.48132 -0.04248 -0.00595 0.03783 0.48299
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.000e+00 3.058e-03 327.049
## hatch_percentile -4.101e-01 1.687e-02 -24.316
## sync_fish 9.646e-07 1.864e-04 0.005
## mean_fishlate -1.619e-06 7.366e-03 0.000
## mean_fishsame 1.499e-05 4.731e-03 0.003
## hatch_percentile:sync_fish -1.625e-02 8.659e-04 -18.771
## hatch_percentile:mean_fishlate 5.176e-01 3.613e-02 14.324
## hatch_percentile:mean_fishsame 2.086e-01 2.540e-02 8.210
## sync_fish:mean_fishlate -1.120e-06 3.843e-04 -0.003
## sync_fish:mean_fishsame -1.283e-06 2.836e-04 -0.005
## hatch_percentile:sync_fish:mean_fishlate -2.946e-02 1.813e-03 -16.243
## hatch_percentile:sync_fish:mean_fishsame -1.433e-02 1.323e-03 -10.839
## Pr(>|t|)
## (Intercept) <2e-16 ***
## hatch_percentile <2e-16 ***
## sync_fish 0.996
## mean_fishlate 1.000
## mean_fishsame 0.997
## hatch_percentile:sync_fish <2e-16 ***
## hatch_percentile:mean_fishlate <2e-16 ***
## hatch_percentile:mean_fishsame 3e-16 ***
## sync_fish:mean_fishlate 0.998
## sync_fish:mean_fishsame 0.996
## hatch_percentile:sync_fish:mean_fishlate <2e-16 ***
## hatch_percentile:sync_fish:mean_fishsame <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.00539373)
##
## Null deviance: 109.445 on 3748 degrees of freedom
## Residual deviance: 20.156 on 3737 degrees of freedom
## (25 observations deleted due to missingness)
## AIC: -8926
##
## Number of Fisher Scoring iterations: 2## (Intercept)
## 9.999893e-01
## hatch_percentile
## -4.101455e-01
## sync_fish
## 9.645720e-07
## mean_fishlate
## -1.618530e-06
## mean_fishsame
## 1.499360e-05
## hatch_percentile:sync_fish
## -1.625421e-02
## hatch_percentile:mean_fishlate
## 5.175879e-01
## hatch_percentile:mean_fishsame
## 2.085683e-01
## sync_fish:mean_fishlate
## -1.119798e-06
## sync_fish:mean_fishsame
## -1.282986e-06
## hatch_percentile:sync_fish:mean_fishlate
## -2.945621e-02
## hatch_percentile:sync_fish:mean_fishsame
## -1.433445e-02Plotting
# adjust data for plotting-- remove some sync levels and put factors in order
m2_fewsync <- m2 %>%
filter(sync_fish == 3 | sync_fish == 15 | sync_fish == 27) %>%
mutate(synchrony_fish = case_when(sync_fish == 3 ~ "high synchrony",
sync_fish == 15 ~ "medium synchrony",
sync_fish == 27 ~ "low synchrony")) %>%
mutate(synchrony_fish = factor(synchrony_fish, levels = c("low", "medium", "high")),
mean_fish = factor(mean_fish, levels = c("early", "same", "late")))
# or use y predicted_survival
ggplot(data = m2_fewsync, aes(x = hatch_percentile, y = coef, color = mean_fish, fill = mean_fish, shape = mean_fish)) +
geom_point() +
geom_smooth(se = F) + #method = "lm", formula = y ~ poly(x, 3)) +
facet_wrap(~synchrony_fish) +
#ylim(0, 1) +
mytheme +
labs(color = "relative mean arrival\nof focal species",
fill = "relative mean arrival\nof focal species",
shape = "relative mean arrival\nof focal species",
x = "relative hatching date",
y = "selection coefficient") +
scale_color_manual(values = c("#01A08A", "#F4AC00", "#FF2401")) +
scale_fill_manual(values = c("#01A08A", "#F4AC00", "#FF2401")) +
scale_shape_manual(values = c(21, 24, 23))
1-Species Models
Load data
Very similar to the 2-species dataset. This data has one row per individual and gives info on treatment parameters (11 levels of synchrony, 3 levels of competitive asymmetry), results of the simulation (number of survivors, biomass), and stats for the individual (hatching date, survival status, development time, final size).
## X run_num id sync asym synchrony asymmetry n_surv n_dead biomass
## 1 1 1 1 0 0 0 none 65 15 409.7
## 2 2 1 2 0 0 0 none 65 15 409.7
## 3 3 1 3 0 0 0 none 65 15 409.7
## 4 4 1 4 0 0 0 none 65 15 409.7
## 5 5 1 5 0 0 0 none 65 15 409.7
## 6 6 1 6 0 0 0 none 65 15 409.7
## mean_size meta final_size hatch_date end_date growth_time growth_rate
## 1 0 1 6.10 46 124 78 0.07820513
## 2 0 1 6.35 47 110 63 0.10079365
## 3 0 0 5.80 47 112 65 0.08923077
## 4 0 1 6.30 47 123 76 0.08289474
## 5 0 1 6.20 47 125 78 0.07948718
## 6 0 0 6.00 47 122 75 0.08000000Survival model
Prepare data
Essentially, pool replicate treatments and calcuate standardized metrics for each treatment.
## Prepare data
survival_model_data_1 <- ind_1 %>%
# remove run_num to pool all replicates of a treatment. also remove others we don't need
select(-c(X, run_num, synchrony, asymmetry, n_surv, n_dead, biomass, mean_size, end_date, growth_time, growth_rate)) %>%
# group by treatment to do the relativazation calculations
group_by(sync, asym) %>%
# standardize arrival date for each treatment combination
mutate(trt_max_hatch_date = max(hatch_date),
trt_min_hatch_date = min(hatch_date),
trt_mean_hatch_date = mean(hatch_date),
hatch_percentile = round((hatch_date - trt_mean_hatch_date)/(trt_max_hatch_date - trt_min_hatch_date), 3))
head(survival_model_data_1)## # A tibble: 6 x 10
## # Groups: sync, asym [1]
## id sync asym meta final_size hatch_date trt_max_hatch_d…
## <int> <int> <dbl> <int> <dbl> <dbl> <dbl>
## 1 1 0 0 1 6.1 46 47
## 2 2 0 0 1 6.35 47 47
## 3 3 0 0 0 5.8 47 47
## 4 4 0 0 1 6.3 47 47
## 5 5 0 0 1 6.2 47 47
## 6 6 0 0 0 6 47 47
## # … with 3 more variables: trt_min_hatch_date <dbl>,
## # trt_mean_hatch_date <dbl>, hatch_percentile <dbl>Model survival
Logistic regression model with survival as a binary outcome
m_surv1 <- glm(data = survival_model_data_1,
meta ~
hatch_percentile + sync + asym +
hatch_percentile*sync +
hatch_percentile*asym +
sync*asym +
hatch_percentile*sync*asym,
family=binomial(link="logit"))
print(summary(m_surv1))##
## Call:
## glm(formula = meta ~ hatch_percentile + sync + asym + hatch_percentile *
## sync + hatch_percentile * asym + sync * asym + hatch_percentile *
## sync * asym, family = binomial(link = "logit"), data = survival_model_data_1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.6025 0.0000 0.0969 0.3713 3.7849
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.721959 0.071314 24.146 < 2e-16 ***
## hatch_percentile -4.799948 0.563571 -8.517 < 2e-16 ***
## sync 0.065950 0.006319 10.436 < 2e-16 ***
## asym 2.042836 0.151396 13.493 < 2e-16 ***
## hatch_percentile:sync -0.726147 0.047994 -15.130 < 2e-16 ***
## hatch_percentile:asym -6.303161 1.176633 -5.357 8.46e-08 ***
## sync:asym -0.109889 0.010679 -10.290 < 2e-16 ***
## hatch_percentile:sync:asym -1.023680 0.108833 -9.406 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27753 on 24799 degrees of freedom
## Residual deviance: 10858 on 24792 degrees of freedom
## (1600 observations deleted due to missingness)
## AIC: 10874
##
## Number of Fisher Scoring iterations: 8## (Intercept) hatch_percentile
## 1.72195870 -4.79994791
## sync asym
## 0.06595003 2.04283629
## hatch_percentile:sync hatch_percentile:asym
## -0.72614682 -6.30316078
## sync:asym hatch_percentile:sync:asym
## -0.10988852 -1.02367964# make a data frame of individuals to add predicted survival to
predict_df1 <- data.frame(hatch_percentile = survival_model_data_1$hatch_percentile,
sync = survival_model_data_1$sync,
asym = survival_model_data_1$asym)
# calculate predicted survival + new value for each individual
m1 <- as_tibble(predict(m_surv1, predict_df1, type = "response")) %>%
mutate(coef = value * (1 - value))
m1 <- cbind(predict_df1, m1)
# calculate selection coefficients -- absolute and standardized
selection_coefficients_1sp <- m1 %>%
group_by(sync, asym) %>%
summarize(selection_coef = mean(coef, na.rm = T), # absolute: 0.06667
standardized_selection_coef = selection_coef/mean(survival_model_data_1$meta, na.rm = T)) # standardized: 0.0889
selection_coefficients_1sp## # A tibble: 33 x 4
## # Groups: sync [11]
## sync asym selection_coef standardized_selection_coef
## <int> <dbl> <dbl> <dbl>
## 1 0 0 0.133 0.177
## 2 0 0.5 0.0614 0.0819
## 3 0 1 0.0251 0.0334
## 4 3 0 0.123 0.164
## 5 3 0.5 0.0895 0.119
## 6 3 1 0.0656 0.0875
## 7 6 0 0.111 0.148
## 8 6 0.5 0.0843 0.112
## 9 6 1 0.0698 0.0931
## 10 9 0 0.103 0.137
## # … with 23 more rowsMass model
Prepare data
Almost the same as before, except we also filter out individuals that did not survive, and calculate relative size in each treatment
# Prepare data
size_model_data_1 <- ind_1 %>%
# remove run_num to pool all replicates of a treatment. also remove others we don't need
select(-c(X, run_num, synchrony, asymmetry, n_surv, n_dead, biomass, mean_size, end_date, growth_time, growth_rate)) %>%
# group by treatment to do the relativazation calculations
group_by(sync, asym) %>%
# standardize fitness by dividing mass by treatment mean
mutate(size_relative = final_size / mean(final_size)) %>%
# standardize arrival date for each treatment combination
mutate(trt_max_hatch_date = max(hatch_date),
trt_min_hatch_date = min(hatch_date),
trt_mean_hatch_date = mean(hatch_date),
hatch_percentile = round((hatch_date - trt_mean_hatch_date)/(trt_max_hatch_date - trt_min_hatch_date), 3))
head(size_model_data_1)## # A tibble: 6 x 11
## # Groups: sync, asym [1]
## id sync asym meta final_size hatch_date size_relative
## <int> <int> <dbl> <int> <dbl> <dbl> <dbl>
## 1 1 0 0 1 6.1 46 0.978
## 2 2 0 0 1 6.35 47 1.02
## 3 3 0 0 0 5.8 47 0.930
## 4 4 0 0 1 6.3 47 1.01
## 5 5 0 0 1 6.2 47 0.994
## 6 6 0 0 0 6 47 0.962
## # … with 4 more variables: trt_max_hatch_date <dbl>,
## # trt_min_hatch_date <dbl>, trt_mean_hatch_date <dbl>,
## # hatch_percentile <dbl># View distribution of final sizes - pretty normal
ggplot(size_model_data_1, aes(x = size_relative)) + geom_histogram() + mytheme
Model mass
## Fitness = mass model
m_size_1 <- glm(data = size_model_data_1,
size_relative ~
hatch_percentile + sync + asym +
hatch_percentile*sync +
hatch_percentile*asym +
sync*asym +
hatch_percentile*sync*asym)
summary(m_size_1)##
## Call:
## glm(formula = size_relative ~ hatch_percentile + sync + asym +
## hatch_percentile * sync + hatch_percentile * asym + sync *
## asym + hatch_percentile * sync * asym, data = size_model_data_1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.88888 -0.05303 0.00365 0.05285 1.16072
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.000e+00 2.071e-03 482.888 <2e-16 ***
## hatch_percentile -1.557e-01 1.342e-02 -11.605 <2e-16 ***
## sync -1.688e-07 1.283e-04 -0.001 0.999
## asym 1.093e-04 3.157e-03 0.035 0.972
## hatch_percentile:sync -2.342e-02 7.422e-04 -31.551 <2e-16 ***
## hatch_percentile:asym -7.195e-01 2.046e-02 -35.163 <2e-16 ***
## sync:asym -6.571e-06 1.875e-04 -0.035 0.972
## hatch_percentile:sync:asym -3.710e-02 1.102e-03 -33.680 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.01204348)
##
## Null deviance: 1852.80 on 24799 degrees of freedom
## Residual deviance: 298.58 on 24792 degrees of freedom
## (1600 observations deleted due to missingness)
## AIC: -39208
##
## Number of Fisher Scoring iterations: 2## (Intercept) hatch_percentile
## 1.000031e+00 -1.557391e-01
## sync asym
## -1.688253e-07 1.093462e-04
## hatch_percentile:sync hatch_percentile:asym
## -2.341789e-02 -7.195091e-01
## sync:asym hatch_percentile:sync:asym
## -6.571219e-06 -3.710324e-02Plotting
# adjust data for plotting-- remove some sync levels and put factors in order
m1_fewsync <- m1 %>%
filter(sync == 3 | sync == 15 | sync == 27) %>%
mutate(synchrony = case_when(sync == 3 ~ "high",
sync == 15 ~ "medium",
sync == 27 ~ "low")) %>%
mutate(asymmetry = case_when(asym == 0 ~ "none",
asym == 0.5 ~ "weak",
asym == 1 ~ "strong")) %>%
mutate(synchrony = factor(synchrony, levels = c("low", "medium", "high")),
asymmetry = factor(asymmetry, levels = c("none", "weak", "strong")))
ggplot(data = m1_fewsync, aes(x = hatch_percentile, y = value, color = asymmetry, fill = asymmetry)) +
geom_point() +
geom_smooth() + #method = "lm", formula = y ~ poly(x, 3)) +
facet_wrap(~synchrony) +
ylim(0, 1) +
mytheme +
labs(color = "competitive asymmetry",
fill ="competitive asymmetry",
x = "relative hatching date",
y = "predicted survival") +
scale_color_manual(values = c("#acbab3", "#fbd364", "#446353")) +
scale_fill_manual(values = c("#acbab3", "#fbd364", "#446353")) 