Is there an interraction between social environment experienced and the susceptibility to challenges? We did two different factorial experiments that coupled manipulations targeted at changing the social environment (signal color manipulation) with manipulations targeted at imposing realistic environmental challenges (flight handicap or simulated predations).
In 2018, manipulations were ordered so that birds were first dulled and then later in the season experienced one of the challenges.
In 2019, the order was reversed so that birds experienced the challenge (only predator, no handicaping in this year) and then subsequently received a dulling manipulation.
For each year, we can ask whether dulling or challenges influenced a number of female measures (mass loss, cort response, provisioning rate, visitation & trips) and nestling outcomes (number fledged, nestling mass, nestling fledge age). Note that in these years we also cross fostered nestlings so that effects on nestlings can be disentangled from those caused by maternal effects vs. experienced nestling environment. Not sure what level of detail we want to go into with nestlings.
Data that we have for these birds but that at the moment I am thinking would not be included in this paper: adult/nestling microbiome, full paternity assignment, nestling/adult feather quality, glucose, some diet.
In this first code chunk I'm loading packages and data to use in later plotting.
# Load packages ----
pacman::p_load(tidyverse, here, lme4, scales)
## not loading plyr because it conflicts with here and tidyverse
# Load data ----
# the data files include similar data from each year but experiments will mainly be analyzed separately
# data by female: one row per female with measures throughout the breeding season
d_fem <- read.delim(here::here("1_raw_data", "data_by_female.txt"))
# Remove records that were noted as exclude. These are mostly cases where either a nest failed so early
# that no treatment was delivered or something weird happened (e.g., treatments mixed up).
d_fem <- subset(d_fem, d_fem$exclude != "Yes")
# data by nestling: one row per nestling with some more individual level nestling info
d_nestling <- read.delim(here::here("1_raw_data", "data_by_nestling.txt"))
# provisioning data: one row per nest per day with provisiong info from rfid sensors
d_provision <- read.delim(here::here("1_raw_data", "daily_provision_data.txt"))
# visitation data: one row per nest per day with social visits and trips for each nest
d_social <- read.delim(here::here("1_raw_data", "nest_visit_data.txt"))
# rfid reads: total rfid reads at each box on each day (used to filter out malfunctioning rfid boards)
d_tot_rfid <- read.delim(here::here("1_raw_data", "total_rfid_reads.txt"))
# Set colors for different treatment groups ----
col_dull <- "slateblue"
col_sham <- "orange"
ch_pred <- "chartreuse3"
ch_con <- "orange"
ch_tape <- "coral3"
# Full treatments
col_CC <- "gray50" # 2018 & 2019, full control
col_CS <- "coral3" # 2018, sham color then taping challenge
col_CP <- "orange" # 2018, sham color then predator
col_DC <- "slateblue" # 2018, dulling folowed by sham coloring
col_DS <- "violet" # 2018, dulling folowed by taping challenge
col_DP <- "chartreuse3" # 2018, dulling followed by predator
col_PC <- "coral3" # 2019, predator first then sham color
col_PD <- "chartreuse3" # 2019, predator first then dulling
col_CD <- "slateblue" # 2019, control first, then dulling
Females lose mass through the breeding season normally. The question here is whether the treatments we applied changed this mass loss trajectory.
## Female Mass Change ----
# Figures showing adult mass at each timepoint by treatment split by year (experiment)
long_mass <- d_fem %>%
pivot_longer(cols = c("fmass1", "fmass2", "fmass3"), names_to = "cap_num",
names_transform = list(
cap_num = ~ readr::parse_factor(.x, levels = c("fmass1", "fmass2", "fmass3"),
ordered = TRUE)),
values_to = "mass", values_drop_na = TRUE)
long_mass <- as.data.frame(long_mass)
long_mass$xpos <- as.numeric(long_mass$cap_num)
sum_mass <- d_fem %>%
pivot_longer(cols = c("fmass1", "fmass2", "fmass3"), names_to = "cap_num",
names_transform = list(
cap_num = ~ readr::parse_factor(.x, levels = c("fmass1", "fmass2", "fmass3"),
ordered = TRUE)),
values_to = "mass", values_drop_na = TRUE) %>%
group_by(as.factor(year), full_treatment, cap_num) %>%
summarise(n = n(), mu = mean(mass, na.rm = TRUE), se = sd(mass, na.rm = TRUE) / sqrt(n()))
sum_mass$x_pos <- as.numeric(sum_mass$cap_num)
for_plot <- data.frame(full_treatment = unique(sum_mass$full_treatment),
p_shape = c(21, 24, 22, 21, 24, 22, 21, 24, 24),
p_col = c(rep(col_sham, 3), rep(col_dull, 3), col_dull, col_sham, col_dull),
p_dodge = c(-.25, -.15, -.05, .05, .15, .25, .0833, -.0833, .25))
sum_mass2 <- as.data.frame(sum_mass)
colnames(sum_mass2)[1] <- "year"
sum_mass2 <- plyr::join(sum_mass2, for_plot, "full_treatment", "left", "first")
long_mass2 <- plyr::join(long_mass, for_plot, "full_treatment", "left", "first")
lm2_18 <- subset(long_mass2, long_mass2$year == "2018")
lm2_19 <- subset(long_mass2, long_mass2$year == "2019")
sm2_18 <- subset(sum_mass2, sum_mass2$year == "2018")
sm2_19 <- subset(sum_mass2, sum_mass2$year == "2019")
par(mfrow = c(1, 2))
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Dulling then Challenge", xlim = c(0.3, 3.7), ylim = c(15.5, 25),
xlab = "Capture Number", ylab = "Female Mass (g)")
axis(1, seq(0, 4, 1))
axis(2, seq(5, 30, 1), las = 2)
points(lm2_18$xpos + lm2_18$p_dodge, lm2_18$mass, pch = 16, col = alpha("black", 0.3))
for(i in 1:nrow(sm2_18)){
lines(rep(sm2_18$x_pos[i] + sm2_18$p_dodge[i], 2),
c(sm2_18$mu[i] - sm2_18$se[i], sm2_18$mu[i] + sm2_18$se[i]), lwd = 2)
}
points(sm2_18$x_pos + sm2_18$p_dodge, sm2_18$mu, pch = sm2_18$p_shape, cex = 1.4,
bg = as.character(sm2_18$p_col))
#legend("topright", c("No Color", "Dull Color", "Control", "Predator", "Handicap"),
# pch = c(21, 21, 21, 24, 22), pt.bg = c(col_sham, col_dull, rep("white", 3)), cex = 1.2, bty = "n")
rect(1.4, 15.7, 3.5, 16.2, col = alpha("chartreuse3", 0.3))
text(2.15, 15.95, "Dulling", pos = 4)
rect(2.2, 16.3, 3.5, 16.8, col = alpha("coral3", 0.2))
text(2.4, 16.55, "Challenge", pos = 4)
#2019
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Challenge then Dulling", xlim = c(0.3, 3.7), ylim = c(15.5, 25),
xlab = "Capture Number", ylab = "Female Mass (g)")
axis(1, seq(0, 4, 1))
axis(2, seq(5, 30, 1), las = 2)
points(lm2_19$xpos + lm2_19$p_dodge, lm2_19$mass, pch = 16, col = alpha("black", 0.3))
for(i in 1:nrow(sm2_19)){
lines(rep(sm2_19$x_pos[i] + sm2_19$p_dodge[i], 2),
c(sm2_19$mu[i] - sm2_19$se[i], sm2_19$mu[i] + sm2_19$se[i]), lwd = 2)
}
points(sm2_19$x_pos + sm2_19$p_dodge, sm2_19$mu, pch = sm2_19$p_shape, cex = 1.4,
bg = as.character(sm2_19$p_col))
legend("topright", c("Sham", "Dulled", "Control", "Predator", "Handicap"),
pch = c(21, 21, 21, 24, 22), pt.bg = c(col_sham, col_dull, rep("white", 3)), cex = 1, bty = "n")
rect(2.3, 15.7, 3.6, 16.2, col = alpha("chartreuse3", 0.3))
text(2.65, 15.95, "Dulling", pos = 4)
rect(1.4, 16.3, 2.5, 16.8, col = alpha("coral3", 0.2))
text(1.5, 16.55, "Challenge", pos = 4)
So the signature of mass loss across captures is very clear in all groups, but there is basically no suggestion that any of our treatments are resulting in a steaper mass loss than any others.
## Female Cort
long_cort <- d_fem %>%
pivot_longer(cols = c("fbase1", "fstr1", "fdex1", "fbase2", "fbase3", "fstr3", "fdex3"), names_to = "cort_type",
names_transform = list(
cap_num = ~ readr::parse_factor(.x, levels = c("fbase1", "fstr1", "fdex1", "fbase2", "fbase3", "fstr3", "fdex3"),
ordered = TRUE)),
values_to = "cort", values_drop_na = TRUE)
long_cort <- as.data.frame(long_cort)
cort_pos <- data.frame(cort_type = c("fbase1", "fstr1", "fdex1", "fbase2", "fbase3", "fstr3", "fdex3"),
xpos = c(1, 2, 3, 5, 7, 8, 9))
long_cort <- plyr::join(long_cort, cort_pos, "cort_type", "left", "first")
sum_cort <- d_fem %>%
pivot_longer(cols = c("fbase1", "fstr1", "fdex1", "fbase2", "fbase3", "fstr3", "fdex3"), names_to = "cort_type",
names_transform = list(
cap_num = ~ readr::parse_factor(.x, levels = c("fbase1", "fstr1", "fdex1", "fbase2", "fbase3", "fstr3", "fdex3"),
ordered = TRUE)),
values_to = "cort", values_drop_na = TRUE) %>%
group_by(as.factor(year), full_treatment, cort_type) %>%
summarise(n = n(), mu = mean(cort, na.rm = TRUE), se = sd(cort, na.rm = TRUE) / sqrt(n()))
sum_cort <- as.data.frame(sum_cort)
sum_cort <- plyr::join(sum_cort, cort_pos, "cort_type", "left", "first")
for_plot <- data.frame(full_treatment = unique(sum_cort$full_treatment),
p_shape = c(21, 24, 22, 21, 24, 22, 21, 24, 24),
p_col = c(rep(col_sham, 3), rep(col_dull, 3), col_dull, col_sham, col_dull),
p_dodge = c(-.25, -.15, -.05, .05, .15, .25, .0833, -.0833, .25))
sum_cort <- plyr::join(sum_cort, for_plot, "full_treatment")
long_cort <- plyr::join(long_cort, for_plot, "full_treatment")
colnames(sum_cort)[2] <- "year"
lc2_18 <- subset(long_cort, long_cort$year == "2018")
lc2_19 <- subset(long_cort, long_cort$year == "2019")
sc2_18 <- subset(sum_cort, sum_cort$year == "2018")
sc2_19 <- subset(sum_cort, sum_cort$year == "2019")
#2018 plot
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Dulling then Challenge", xlim = c(0.3, 9.7), ylim = c(-20, 80),
xlab = "", ylab = "Corticosterone (ng/ul)")
axis(1, c(-1, 1, 2, 3, 5, 7, 8, 9, 15), c("", "B 1", "S 1", "D 1", "B 2", "B 3", "S 3", "D 3", ""))
axis(2, c(-50, seq(0, 200, 10)), las = 2)
points(lc2_18$xpos + lc2_18$p_dodge, lc2_18$cort, pch = 16, col = alpha("black", 0.3))
abline(h = 0)
for(i in 1:nrow(sc2_18)){
lines(rep(sc2_18$xpos[i] + sc2_18$p_dodge[i], 2),
c(sc2_18$mu[i] - sc2_18$se[i], sc2_18$mu[i] + sc2_18$se[i]), lwd = 2)
}
points(sc2_18$xpos + sc2_18$p_dodge, sc2_18$mu, pch = sc2_18$p_shape, cex = 1.4,
bg = as.character(sc2_18$p_col))
rect(4.4, -15, 9.6, -10, col = alpha("chartreuse3", 0.3))
text(6.5, -12.5, "Dulling", pos = 4)
rect(6.4, -9, 9.6, -4, col = alpha("coral3", 0.2))
text(7.3, -6.5, "Challenge", pos = 4)
legend("bottomleft", c("No Color", "Dull Color", "Control", "Predator", "Handicap"),
pch = c(21, 21, 21, 24, 22), pt.bg = c(col_sham, col_dull, rep("white", 3)), cex =0.9, bty = "n")
#2019 plot
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Challenge then Dulling", xlim = c(0.3, 9.7), ylim = c(-20, 80),
xlab = "", ylab = "Corticosterone (ng/ul)")
axis(1, c(-1, 1, 2, 3, 5, 7, 8, 9, 15), c("", "B 1", "S 1", "D 1", "B 2", "B 3", "S 3", "A 3", ""))
axis(2, c(-50, seq(0, 200, 10)), las = 2)
points(lc2_19$xpos + lc2_19$p_dodge, lc2_19$cort, pch = 16, col = alpha("black", 0.3))
abline(h = 0)
for(i in 1:nrow(sc2_19)){
lines(rep(sc2_19$xpos[i] + sc2_19$p_dodge[i], 2),
c(sc2_19$mu[i] - sc2_19$se[i], sc2_19$mu[i] + sc2_19$se[i]), lwd = 2)
}
points(sc2_19$xpos + sc2_19$p_dodge, sc2_19$mu, pch = sc2_19$p_shape, cex = 1.4,
bg = as.character(sc2_19$p_col))
rect(6, -15, 9.6, -10, col = alpha("chartreuse3", 0.3))
text(7, -12.5, "Dulling", pos = 4)
rect(3.8, -9, 6.6, -4, col = alpha("coral3", 0.2))
text(4.5, -6.5, "Challenge", pos = 4)
legend("bottomleft", c("No Color", "Dull Color", "Control", "Predator"),
pch = c(21, 21, 21, 24), pt.bg = c(col_sham, col_dull, rep("white", 2)), cex =0.9, bty = "n")
There doesn't seem to be much going on with cort differences between groups. Maybe in the first year the predator groups have higher stress response and dex than the non-predator groups, but that certainly isn't very strong. In any case, there doesn't seem to be any interaction between coloring treatment and challenge treatment.
## Provisioning ----
# Figure showing nestling provisioning
d_fem$uby <- paste(d_fem$unitbox, d_fem$year, sep = "_")
d_prov2 <- plyr::join(d_provision, d_fem, "uby", "left", "first")
d_prov3 <- subset(d_prov2, d_prov2$offset < 16 & d_prov2$doy < d_prov2$end_soc_date & d_prov2$f_reads > 50)
d_prov3 <- d_prov3[, c("offset", "f_feed", "m_feed", "year", "full_treatment", "experiment")]
d_prov3_18 <- subset(d_prov3, d_prov3$year == 2018)
ggplot(d_prov3_18, aes(x = offset, y = f_feed, col = full_treatment)) +
geom_point() + geom_smooth() + xlim(0, 15) + xlab("Days after hatching") +
ylab("Daily female provisioning trips") + ggtitle("Dulling then Challenge") +
theme_bw()
ggplot(d_prov3_18, aes(x = offset, y = m_feed, col = full_treatment)) +
geom_point() + geom_smooth() + xlim(0, 15) + xlab("Days after hatching") +
ylab("Daily male provisioning trips") + ggtitle("Dulling then Challenge") +
theme_bw()
d_prov3_19 <- subset(d_prov3, d_prov3$year == 2019)
ggplot(d_prov3_19, aes(x = offset, y = f_feed, col = full_treatment)) +
geom_point() + geom_smooth() + xlim(0.6, 14.5) + xlab("Days after hatching") +
ylab("Daily female provisioning trips") + ggtitle("Challenge then Dulling") +
theme_bw()
ggplot(d_prov3_19, aes(x = offset, y = m_feed, col = full_treatment)) +
geom_point() + geom_smooth() + xlim(0.6, 14.5) + xlab("Days after hatching") +
ylab("Daily male provisioning trips") + ggtitle("Challenge then Dulling") +
theme_bw()
There isn't anything super obvious going on with provisioning rate in any of the treatment groups. Maybe some hints of lower overall rates in predator groups, but not much. I could model these in a better way looking at hourly feeding and controlling for some other things like temperature and time of day, etc.
## Visits to box
soc_vis <- data.frame(uby = unique(d_fem$uby))
for(i in 1:nrow(soc_vis)){
sub <- subset(d_fem, d_fem$uby == soc_vis$uby[i])
if(sub$end_soc_date[1] - sub$cap1[1] + 1 > 0){
dats <- seq(from = sub$cap1[1] + 1, to = sub$end_soc_date[1], by = 1)
temp <- data.frame(uby = rep(sub$uby, length(dats)),
doy = dats)
if(i == 1){
soc_visits <- temp
}
if(i > 1){
soc_visits <- rbind(soc_visits, temp)
}
}
}
soc_visits <- plyr::join(soc_visits, d_fem, "uby", "left", "first")
#for(i in 1:nrow(soc_visits)){
# sub1 <- subset(d_social, as.character(d_social$ubox) == as.character(soc_visits$unitbox[i])
# & d_social$doy == soc_visits$doy[i])
# subf <- subset(sub1, sub1$sex == "Female")
# subm <- subset(sub1, sub1$sex == "Male")
# soc_visits$tot_f_vis[i] <- nrow(subf)
# soc_visits$tot_m_vis[i] <- nrow(subm)
# soc_visits$uni_f_vis[i] <- length(unique(subf$rfid))
# soc_visits$uni_m_vis[i] <- length(unique(subm$rfid))
# print(paste(i, " of ", nrow(soc_visits), sep = ""))
#}
for(i in 1:nrow(d_fem)){
sub <- subset(d_social, as.character(d_social$ubox) == as.character(d_fem$unitbox[i]))
subf <- subset(sub, sub$sex == "Female")
subm <- subset(sub, sub$sex == "Male")
subf1 <- subset(subf, subf$doy >= d_fem$cap1[i] & subf$doy < d_fem$cap2[i])
subf2 <- subset(subf, subf$doy >= d_fem$cap2[i])
subm1 <- subset(subm, subm$doy >= d_fem$cap1[i] & subm$doy < d_fem$cap2[i])
subm2 <- subset(subm, subm$doy >= d_fem$cap2[i])
d_fem$tot_f_vis[i] <- nrow(subf)
d_fem$tot_m_vis[i] <- nrow(subm)
d_fem$uni_f_vis[i] <- length(unique(subf$rfid))
d_fem$uni_m_vis[i] <- length(unique(subm$rfid))
d_fem$tot_f_vis1[i] <- nrow(subf1)
d_fem$tot_m_vis1[i] <- nrow(subm1)
d_fem$uni_f_vis1[i] <- length(unique(subf1$rfid))
d_fem$uni_m_vis1[i] <- length(unique(subm1$rfid))
d_fem$tot_f_vis2[i] <- nrow(subf2)
d_fem$tot_m_vis2[i] <- nrow(subm2)
d_fem$uni_f_vis2[i] <- length(unique(subf2$rfid))
d_fem$uni_m_vis2[i] <- length(unique(subm2$rfid))
}
##visitors to the box
long_vis <- d_fem %>%
pivot_longer(cols = c("uni_f_vis1", "uni_m_vis1", "uni_f_vis2", "uni_m_vis2"), names_to = "type",
names_transform = list(
type = ~ readr::parse_factor(.x, levels = c("uni_f_vis1", "uni_m_vis1", "uni_f_vis2", "uni_m_vis2"),
ordered = TRUE)),
values_to = "count", values_drop_na = TRUE)
long_vis <- as.data.frame(long_vis)
vis_pos <- data.frame(type = c("uni_f_vis1", "uni_m_vis1", "uni_f_vis2", "uni_m_vis2"),
xpos = seq(1, 4, 1))
long_vis <- plyr::join(long_vis, vis_pos, "type", "left", "first")
sum_vis <- d_fem %>%
pivot_longer(cols = c("uni_f_vis1", "uni_m_vis1", "uni_f_vis2", "uni_m_vis2"), names_to = "type",
names_transform = list(
type = ~ readr::parse_factor(.x, levels = c("uni_f_vis1", "uni_m_vis1", "uni_f_vis2", "uni_m_vis2"),
ordered = TRUE)),
values_to = "count", values_drop_na = TRUE) %>%
group_by(as.factor(year), full_treatment, type) %>%
summarise(n = n(), mu = mean(count, na.rm = TRUE), se = sd(count, na.rm = TRUE) / sqrt(n()))
sum_vis <- as.data.frame(sum_vis)
sum_vis <- plyr::join(sum_vis, vis_pos, "type", "left", "first")
colnames(sum_vis)[2] <- "year"
for_plot <- data.frame(full_treatment = unique(sum_vis$full_treatment),
p_shape = c(21, 24, 22, 21, 24, 22, 21, 24, 24),
p_col = c(rep(col_sham, 3), rep(col_dull, 3), col_dull, col_sham, col_dull),
p_dodge = c(-.25, -.15, -.05, .05, .15, .25, .0833, -.0833, .25))
sum_vis <- plyr::join(sum_vis, for_plot, "full_treatment")
long_vis <- plyr::join(long_vis, for_plot, "full_treatment")
vis2_18 <- subset(long_vis, long_vis$year == "2018")
vis2_19 <- subset(long_vis, long_vis$year == "2019")
svis2_18 <- subset(sum_vis, sum_vis$year == "2018")
svis2_19 <- subset(sum_vis, sum_vis$year == "2019")
## Visitors
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Dulling then Challenge", xlim = c(0.5, 4.5), ylim = c(-2, 25),
xlab = "", ylab = "Number of Unique Visitors")
axis(1, c(-1, 1, 2, 3, 4, 10), c("", "Females 1", "Males 1", "Females 2", "Males 2", ""))
axis(2, seq(-30, 200, 5), las = 2)
abline(h = 0)
points(vis2_18$xpos + vis2_18$p_dodge, vis2_18$count, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(svis2_18)){
lines(rep(svis2_18$xpos[i] + svis2_18$p_dodge[i], 2),
c(svis2_18$mu[i] - svis2_18$se[i], svis2_18$mu[i] + svis2_18$se[i]), lwd = 2)
}
points(svis2_18$xpos + svis2_18$p_dodge, svis2_18$mu, pch = svis2_18$p_shape, cex = 1.4,
bg = as.character(svis2_18$p_col))
legend(.5, 23, c("No Color", "Dull Color", "Control", "Predator", "Handicap"),
pch = c(21, 21, 21, 24, 22), pt.bg = c(col_sham, col_dull, rep("white", 3)), cex =0.9, bty = "n")
rect(0.5, 23.5, 4.5, 25, col = alpha("chartreuse3", 0.3))
text(2.3, 24.25, "Dulling", pos = 4)
rect(2.5, 21.5, 4.5, 23, col = alpha("coral3", 0.2))
text(3.2, 22.25, "Challenge", pos = 4)
lines(c(2.5, 2.5), c(-5, 23.5), lty = 2)
text(1.5, -1, "First to second capture")
text(3.5, -1, "After second capture")
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Challenge then Dulling", xlim = c(0.5, 4.5), ylim = c(-2, 25),
xlab = "", ylab = "Number of Unique Visitors")
axis(1, c(-1, 1, 2, 3, 4, 10), c("", "Females 1", "Males 1", "Females 2", "Males 2", ""))
axis(2, seq(-30, 200, 5), las = 2)
abline(h = 0)
points(vis2_19$xpos + vis2_19$p_dodge, vis2_19$count, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(svis2_19)){
lines(rep(svis2_19$xpos[i] + svis2_19$p_dodge[i], 2),
c(svis2_19$mu[i] - svis2_19$se[i], svis2_19$mu[i] + svis2_19$se[i]), lwd = 2)
}
points(svis2_19$xpos + svis2_19$p_dodge, svis2_19$mu, pch = svis2_19$p_shape, cex = 1.4,
bg = as.character(svis2_19$p_col))
legend(.5, 22, c("No Color", "Dull Color", "Control", "Predator"),
pch = c(21, 21, 21, 24), pt.bg = c(col_sham, col_dull, rep("white", 2)), cex =0.9, bty = "n")
rect(2.5, 23.5, 4.5, 25, col = alpha("chartreuse3", 0.3))
text(3.3, 24.25, "Dulling", pos = 4)
rect(0.5, 21.5, 2.5, 23, col = alpha("coral3", 0.2))
text(1.2, 22.25, "Challenge", pos = 4)
lines(c(2.5, 2.5), c(-5, 23.5), lty = 2)
text(1.5, -1, "First to second capture")
text(3.5, -1, "After second capture")
Hey! An actual difference that makes sense. It seems that dulling in these two years pretty consistently resulted in more unique visitors (both male and female) coming to the box. This is done a little bit of a different (simpler) way than in the previous year experiment. I'm only looking at days 6-15 and I haven't (yet) thought about initial coloration or interactions with different breeding stages. We did that in the preprint for 2017, although it gets really messy to explain and this is already complicated... For the 'challenge then dulling' experiment I really should split up at least into two periods since those birds actually were not dulled until after the second capture and that might be resulting in what looks like a slightly smaller effect size than in 2018. In 2019, there looks like maybe a trend for an interaction with the predation treatment too (and remember that predation treatment happened during incubation so wouldn't be too surprising if effect on visitors was a bit different).
One thing that I want to check on is that the predator difference (especially in the second year) might be driven in part by earlier failure in those predator nests so that they have fewer days to accrue unique visitors. In the previous paper the model we used was uniqe visitors each day as the response with multiple days of observation from the same bird controlled by a random effect. That probably makes sense to do again here but I just need to work on the code to get it into shape for that.
### Trips to other boxes ----
## trips to other boxes
for(i in 1:nrow(d_fem)){
sub <- subset(d_social, as.character(d_social$rfid) == as.character(d_fem$fRFID[i]))
sub1 <- subset(sub, sub$doy >= d_fem$cap1[i] & sub$doy < d_fem$cap2[i])
sub2 <- subset(sub, sub$doy >= d_fem$cap2[i])
d_fem$tot_trip1[i] <- nrow(sub1)
d_fem$tot_trip2[i] <- nrow(sub2)
d_fem$uni_trip1[i] <- length(unique(sub1$ubox))
d_fem$uni_trip2[i] <- length(unique(sub2$ubox))
}
##trips to other boxes
long_vis <- d_fem %>%
pivot_longer(cols = c("uni_trip1", "uni_trip2"), names_to = "type",
names_transform = list(
type = ~ readr::parse_factor(.x, levels = c("uni_trip1", "uni_trip2"),
ordered = TRUE)),
values_to = "count", values_drop_na = TRUE)
long_vis <- as.data.frame(long_vis)
vis_pos <- data.frame(type = c("uni_trip1", "uni_trip2"),
xpos = seq(1, 4, 1))
long_vis <- plyr::join(long_vis, vis_pos, "type", "left", "first")
sum_vis <- d_fem %>%
pivot_longer(cols = c("uni_trip1", "uni_trip2"), names_to = "type",
names_transform = list(
type = ~ readr::parse_factor(.x, levels = c("uni_trip1", "uni_trip2"),
ordered = TRUE)),
values_to = "count", values_drop_na = TRUE) %>%
group_by(as.factor(year), full_treatment, type) %>%
summarise(n = n(), mu = mean(count, na.rm = TRUE), se = sd(count, na.rm = TRUE) / sqrt(n()))
sum_vis <- as.data.frame(sum_vis)
sum_vis <- plyr::join(sum_vis, vis_pos, "type", "left", "first")
colnames(sum_vis)[2] <- "year"
for_plot <- data.frame(full_treatment = unique(sum_vis$full_treatment),
p_shape = c(21, 24, 22, 21, 24, 22, 21, 24, 24),
p_col = c(rep(col_sham, 3), rep(col_dull, 3), col_dull, col_sham, col_dull),
p_dodge = c(-.25, -.15, -.05, .05, .15, .25, .0833, -.0833, .25))
sum_vis <- plyr::join(sum_vis, for_plot, "full_treatment")
long_vis <- plyr::join(long_vis, for_plot, "full_treatment")
vis2_18 <- subset(long_vis, long_vis$year == "2018")
vis2_19 <- subset(long_vis, long_vis$year == "2019")
svis2_18 <- subset(sum_vis, sum_vis$year == "2018")
svis2_19 <- subset(sum_vis, sum_vis$year == "2019")
## Trips
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Dulling then Challenge", xlim = c(0.5, 2.5), ylim = c(-.1, 25),
xlab = "", ylab = "Trips to Unique Boxes")
axis(1, c(-1, 1, 2, 10), c("", "Capture 1-2", "After Capture 2", ""))
axis(2, seq(-30, 200, 5), las = 2)
points(vis2_18$xpos + vis2_18$p_dodge, vis2_18$count, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(svis2_18)){
lines(rep(svis2_18$xpos[i] + svis2_18$p_dodge[i], 2),
c(svis2_18$mu[i] - svis2_18$se[i], svis2_18$mu[i] + svis2_18$se[i]), lwd = 2)
}
points(svis2_18$xpos + svis2_18$p_dodge, svis2_18$mu, pch = svis2_18$p_shape, cex = 1.4,
bg = as.character(svis2_18$p_col))
legend(.5, 23, c("No Color", "Dull Color", "Control", "Predator", "Handicap"),
pch = c(21, 21, 21, 24, 22), pt.bg = c(col_sham, col_dull, rep("white", 3)), cex =0.9, bty = "n")
rect(0.5, 23.5, 2.5, 25, col = alpha("chartreuse3", 0.3))
text(1.3, 24.25, "Dulling", pos = 4)
rect(1.5, 21.5, 2.5, 23, col = alpha("coral3", 0.2))
text(1.9, 22.25, "Challenge", pos = 4)
lines(c(1.5, 1.5), c(-5, 23.5), lty = 2)
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Challenge then Dulling", xlim = c(0.5, 2.5), ylim = c(-.1, 25),
xlab = "", ylab = "Trips to Unique Boxes")
axis(1, c(-1, 1, 2, 10), c("", "Capture 1-2", "After Capture 2", ""))
axis(2, seq(-30, 200, 5), las = 2)
points(vis2_19$xpos + vis2_19$p_dodge, vis2_19$count, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(svis2_19)){
lines(rep(svis2_19$xpos[i] + svis2_19$p_dodge[i], 2),
c(svis2_19$mu[i] - svis2_19$se[i], svis2_19$mu[i] + svis2_19$se[i]), lwd = 2)
}
points(svis2_19$xpos + svis2_19$p_dodge, svis2_19$mu, pch = svis2_19$p_shape, cex = 1.4,
bg = as.character(svis2_19$p_col))
legend(.5, 21.5, c("No Color", "Dull Color", "Control", "Predator"),
pch = c(21, 21, 21, 24), pt.bg = c(col_sham, col_dull, rep("white", 2)), cex =0.9, bty = "n")
rect(1.5, 23.5, 2.5, 25, col = alpha("chartreuse3", 0.3))
text(1.9, 24.25, "Dulling", pos = 4)
rect(0.5, 21.5, 1.5, 23, col = alpha("coral3", 0.2))
text(0.9, 22.25, "Challenge", pos = 4)
lines(c(1.5, 1.5), c(-5, 23.5), lty = 2)
So for the 2018 experiment it looks like dulled birds generally also made more trips to other boxes, but you can only really see this in the second time period when more trips are happening. In 2019 there isn't really much evidence for a difference, except that maybe the predator treatment birds made a bit fewer trips (although again those ones failed earlier so could be some effect of fewer days of observation, should model it like above).
I guess one interpretation for there being a trend towards social effects in 2019 but not as strong as in 2018 is just that the dulling wasn't applied until later and only applied twice instead of three times, so maybe it's expected to be a less obvious impact overall.
## Reproductive Success
long_rs <- d_fem %>%
pivot_longer(cols = c("clutch", "maxbrood", "numd6", "numband", "num_d15", "numfled"), names_to = "time_point",
names_transform = list(
time_point = ~ readr::parse_factor(.x, levels = c("clutch", "maxbrood", "numd6", "numband", "num_d15", "numfled"),
ordered = TRUE)),
values_to = "count", values_drop_na = TRUE)
long_rs <- as.data.frame(long_rs)
rs_pos <- data.frame(time_point = c("clutch", "maxbrood", "numd6", "numband", "num_d15", "numfled"),
xpos = seq(1, 6, 1))
long_rs <- plyr::join(long_rs, rs_pos, "time_point", "left", "first")
sum_rs <- d_fem %>%
pivot_longer(cols = c("clutch", "maxbrood", "numd6", "numband", "num_d15", "numfled"), names_to = "time_point",
names_transform = list(
time_point = ~ readr::parse_factor(.x, levels = c("clutch", "maxbrood", "numd6", "numband", "num_d15", "numfled"),
ordered = TRUE)),
values_to = "count", values_drop_na = TRUE) %>%
group_by(as.factor(year), full_treatment, time_point) %>%
summarise(n = n(), mu = mean(count, na.rm = TRUE), se = sd(count, na.rm = TRUE) / sqrt(n()))
sum_rs <- as.data.frame(sum_rs)
sum_rs <- plyr::join(sum_rs, rs_pos, "time_point", "left", "first")
colnames(sum_rs)[2] <- "year"
for_plot <- data.frame(full_treatment = unique(sum_rs$full_treatment),
p_shape = c(21, 24, 22, 21, 24, 22, 21, 24, 24),
p_col = c(rep(col_sham, 3), rep(col_dull, 3), col_dull, col_sham, col_dull),
p_dodge = c(-.25, -.15, -.05, .05, .15, .25, .0833, -.0833, .25))
sum_rs <- plyr::join(sum_rs, for_plot, "full_treatment")
long_rs <- plyr::join(long_rs, for_plot, "full_treatment")
rs2_18 <- subset(long_rs, long_rs$year == "2018")
rs2_19 <- subset(long_rs, long_rs$year == "2019")
srs2_18 <- subset(sum_rs, sum_rs$year == "2018")
srs2_19 <- subset(sum_rs, sum_rs$year == "2019")
#2018 plot
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Dulling then Challenge", xlim = c(0.5, 6.5), ylim = c(-0.1, 7.1),
xlab = "", ylab = "Count")
axis(1, c(-1, 1, 2, 3, 4, 5, 6, 10), c("", "Clutch", "Hatch", "Day 6", "Day 12", "Day 15", "Fledged", ""))
axis(2, seq(-1, 10, 1), las = 2)
points(rs2_18$xpos + rs2_18$p_dodge, rs2_18$count, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(srs2_18)){
lines(rep(srs2_18$xpos[i] + srs2_18$p_dodge[i], 2),
c(srs2_18$mu[i] - srs2_18$se[i], srs2_18$mu[i] + srs2_18$se[i]), lwd = 2)
}
points(srs2_18$xpos + srs2_18$p_dodge, srs2_18$mu, pch = srs2_18$p_shape, cex = 1.4,
bg = as.character(srs2_18$p_col))
abline(v = 1.5, lty = 2)
abline(v = 2.5, lty = 2)
abline(v = 3.5, lty = 2)
abline(v = 4.5, lty = 2)
abline(v = 5.5, lty = 2)
rect(0, 6.6, 3.5, 7, col = alpha("chartreuse3", 0.3))
text(1.7, 6.8, "Dulling", pos = 4)
rect(2, 6, 3, 6.4, col = alpha("coral3", 0.2))
text(2.1, 6.2, "Challenge", pos = 4)
legend("bottomleft", c("No Color", "Dull Color", "Control", "Predator", "Handicap"),
pch = c(21, 21, 21, 24, 22), pt.bg = c(col_sham, col_dull, rep("white", 3)), cex =0.9, bty = "n")
#2019 plot
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Challenge then Dulling", xlim = c(0.5, 6.5), ylim = c(-0.1, 7.1),
xlab = "", ylab = "Count")
axis(1, c(-1, 1, 2, 3, 4, 5, 6, 10), c("", "Clutch", "Hatch", "Day 6", "Day 12", "Day 15", "Fledged", ""))
axis(2, seq(-1, 10, 1), las = 2)
points(rs2_19$xpos + rs2_19$p_dodge, rs2_19$count, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(srs2_19)){
lines(rep(srs2_19$xpos[i] + srs2_19$p_dodge[i], 2),
c(srs2_19$mu[i] - srs2_19$se[i], srs2_19$mu[i] + srs2_19$se[i]), lwd = 2)
}
points(srs2_19$xpos + srs2_19$p_dodge, srs2_19$mu, pch = srs2_19$p_shape, cex = 1.4,
bg = as.character(srs2_19$p_col))
abline(v = 1.5, lty = 2)
abline(v = 2.5, lty = 2)
abline(v = 3.5, lty = 2)
abline(v = 4.5, lty = 2)
abline(v = 5.5, lty = 2)
rect(0, 6.6, 3.5, 7, col = alpha("chartreuse3", 0.3))
text(1.7, 6.8, "Dulling", pos = 4)
text(.9, 6.2, "Challenge", pos = 4)
arrows(0.5, 6.2, 0.9, 6.2, code = 1)
legend("bottomleft", c("No Color", "Dull Color", "Control", "Predator"),
pch = c(21, 21, 21, 24), pt.bg = c(col_sham, col_dull, rep("white", 2)), cex =0.9, bty = "n")
So there does seem to be a pretty clear effect of predator treatment on fledging success. In 2018 the pattern is very clear and appears right after the treatments were delivered (days 1-5 after hatching). It also looks like there could be a moderate decrease in fledging success from the handicap (wing taping) treatment, but not nearly as pronounced. In 2019, the same pattern is present though maybe a little bit weaker. Interestingly, in 2019 the challenges were applied during incubation but even though there was no difference in number hatched in the predator treatments, they had reduced numbers of nestlings later on (not sure how strong this evidence is).
## Nestling Mass ----
long_nm <- d_fem %>%
pivot_longer(cols = c("d6_av_mass", "d12_av_mass", "d15_av_mass"), names_to = "time_point",
names_transform = list(
time_point = ~ readr::parse_factor(.x, levels = c("d6_av_mass", "d12_av_mass", "d15_av_mass"),
ordered = TRUE)),
values_to = "mass", values_drop_na = TRUE)
long_nm <- as.data.frame(long_nm)
nm_pos <- data.frame(time_point = c("d6_av_mass", "d12_av_mass", "d15_av_mass"),
xpos = seq(1, 3, 1))
long_nm <- plyr::join(long_nm, nm_pos, "time_point", "left", "first")
sum_nm <- d_fem %>%
pivot_longer(cols = c("d6_av_mass", "d12_av_mass", "d15_av_mass"), names_to = "time_point",
names_transform = list(
time_point = ~ readr::parse_factor(.x, levels = c("d6_av_mass", "d12_av_mass", "d15_av_mass"),
ordered = TRUE)),
values_to = "mass", values_drop_na = TRUE) %>%
group_by(as.factor(year), full_treatment, time_point) %>%
summarise(n = n(), mu = mean(mass, na.rm = TRUE), se = sd(mass, na.rm = TRUE) / sqrt(n()))
sum_nm <- as.data.frame(sum_nm)
sum_nm <- plyr::join(sum_nm, nm_pos, "time_point", "left", "first")
colnames(sum_nm)[2] <- "year"
for_plot <- data.frame(full_treatment = unique(sum_nm$full_treatment),
p_shape = c(21, 24, 22, 21, 24, 22, 21, 24, 24),
p_col = c(rep(col_sham, 3), rep(col_dull, 3), col_dull, col_sham, col_dull),
p_dodge = c(-.25, -.15, -.05, .05, .15, .25, .0833, -.0833, .25))
sum_nm <- plyr::join(sum_nm, for_plot, "full_treatment")
long_nm <- plyr::join(long_nm, for_plot, "full_treatment")
nm2_18 <- subset(long_nm, long_nm$year == "2018")
nm2_19 <- subset(long_nm, long_nm$year == "2019")
snm2_18 <- subset(sum_nm, sum_nm$year == "2018")
snm2_19 <- subset(sum_nm, sum_nm$year == "2019")
##hbill and wing
long_hbw <- d_fem %>%
pivot_longer(cols = c("d12_av_hb", "d12_av_wing"), names_to = "measure",
names_transform = list(
measure = ~ readr::parse_factor(.x, levels = c("d12_av_hb", "d12_av_wing"),
ordered = TRUE)),
values_to = "length", values_drop_na = TRUE)
long_hbw <- as.data.frame(long_hbw)
hbw_pos <- data.frame(measure = c("d12_av_hb", "d12_av_wing"),
xpos = seq(1, 2, 1))
long_hbw <- plyr::join(long_hbw, hbw_pos, "measure", "left", "first")
sum_hbw <- d_fem %>%
pivot_longer(cols = c("d12_av_hb", "d12_av_wing"), names_to = "measure",
names_transform = list(
measure = ~ readr::parse_factor(.x, levels = c("d12_av_hb", "d12_av_wing"),
ordered = TRUE)),
values_to = "length", values_drop_na = TRUE) %>%
group_by(as.factor(year), full_treatment, measure) %>%
summarise(n = n(), mu = mean(length, na.rm = TRUE), se = sd(length, na.rm = TRUE) / sqrt(n()))
sum_hbw <- as.data.frame(sum_hbw)
sum_hbw <- plyr::join(sum_hbw, hbw_pos, "measure", "left", "first")
colnames(sum_hbw)[2] <- "year"
for_plot <- data.frame(full_treatment = unique(sum_hbw$full_treatment),
p_shape = c(21, 24, 22, 21, 24, 22, 21, 24, 24),
p_col = c(rep(col_sham, 3), rep(col_dull, 3), col_dull, col_sham, col_dull),
p_dodge = c(-.25, -.15, -.05, .05, .15, .25, .0833, -.0833, .25))
sum_hbw <- plyr::join(sum_hbw, for_plot, "full_treatment")
long_hbw <- plyr::join(long_hbw, for_plot, "full_treatment")
hbw2_18 <- subset(long_hbw, long_hbw$year == "2018")
hbw2_19 <- subset(long_hbw, long_hbw$year == "2019")
shbw2_18 <- subset(sum_hbw, sum_hbw$year == "2018")
shbw2_19 <- subset(sum_hbw, sum_hbw$year == "2019")
#2018 plot
par(mfrow = c(1, 1))
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Dulling then Challenge", xlim = c(0.5, 3.5), ylim = c(-0.1, 25),
xlab = "", ylab = "Average Nestling Mass (g)")
axis(1, c(-1, 1, 2, 3, 10), c("", "Day 6", "Day 12", "Day 15", ""))
axis(2, seq(-30, 30, 5), las = 2)
points(nm2_18$xpos + nm2_18$p_dodge, nm2_18$mass, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(snm2_18)){
lines(rep(snm2_18$xpos[i] + snm2_18$p_dodge[i], 2),
c(snm2_18$mu[i] - snm2_18$se[i], snm2_18$mu[i] + snm2_18$se[i]), lwd = 2)
}
points(snm2_18$xpos + snm2_18$p_dodge, snm2_18$mu, pch = snm2_18$p_shape, cex = 1.4,
bg = as.character(snm2_18$p_col))
abline(v = 1.5, lty = 2)
abline(v = 2.5, lty = 2)
legend("bottomright", c("No Color", "Dull Color", "Control", "Predator", "Handicap"),
pch = c(21, 21, 21, 24, 22), pt.bg = c(col_sham, col_dull, rep("white", 3)), cex =0.9, bty = "n")
## HB and wing
par(mfrow = c(1, 2))
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Dulling then Challenge", xlim = c(1.5, 2.5), ylim = c(15, 60),
xlab = "", ylab = "Length on Day 12 (mm)")
axis(1, c(-1, 1, 2, 10), c("", "Head + Bill", "Wing", ""))
axis(2, seq(-30, 200, 5), las = 2)
points(hbw2_18$xpos + hbw2_18$p_dodge, hbw2_18$length, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(shbw2_18)){
lines(rep(shbw2_18$xpos[i] + shbw2_18$p_dodge[i], 2),
c(shbw2_18$mu[i] - shbw2_18$se[i], shbw2_18$mu[i] + shbw2_18$se[i]), lwd = 2)
}
points(shbw2_18$xpos + shbw2_18$p_dodge, shbw2_18$mu, pch = shbw2_18$p_shape, cex = 1.4,
bg = as.character(shbw2_18$p_col))
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Dulling then Challenge", xlim = c(0.5, 1.5), ylim = c(19, 29),
xlab = "", ylab = "Length on Day 12 (mm)")
axis(1, c(-1, 1, 2, 10), c("", "Head + Bill", "Wing", ""))
axis(2, seq(-30, 200, 2), las = 2)
points(hbw2_18$xpos + hbw2_18$p_dodge, hbw2_18$length, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(shbw2_18)){
lines(rep(shbw2_18$xpos[i] + shbw2_18$p_dodge[i], 2),
c(shbw2_18$mu[i] - shbw2_18$se[i], shbw2_18$mu[i] + shbw2_18$se[i]), lwd = 2)
}
points(shbw2_18$xpos + shbw2_18$p_dodge, shbw2_18$mu, pch = shbw2_18$p_shape, cex = 1.4,
bg = as.character(shbw2_18$p_col))
legend("bottomright", c("No Color", "Dull Color", "Control", "Predator", "Handicap"),
pch = c(21, 21, 21, 24, 22), pt.bg = c(col_sham, col_dull, rep("white", 3)), cex =0.9, bty = "n")
## 2019 plot
par(mfrow = c(1, 1))
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Challenge then Dulling", xlim = c(0.5, 3.5), ylim = c(-0.1, 25),
xlab = "", ylab = "Average Nestling Mass (g)")
axis(1, c(-1, 1, 2, 3, 10), c("", "Day 6", "Day 12", "Day 15", ""))
axis(2, seq(-30, 30, 5), las = 2)
points(nm2_19$xpos + nm2_19$p_dodge, nm2_19$mass, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(snm2_19)){
lines(rep(snm2_19$xpos[i] + snm2_19$p_dodge[i], 2),
c(snm2_19$mu[i] - snm2_19$se[i], snm2_19$mu[i] + snm2_19$se[i]), lwd = 2)
}
points(snm2_19$xpos + snm2_19$p_dodge, snm2_19$mu, pch = snm2_19$p_shape, cex = 1.4,
bg = as.character(snm2_19$p_col))
abline(v = 1.5, lty = 2)
abline(v = 2.5, lty = 2)
legend("bottomright", c("No Color", "Dull Color", "Control", "Predator", "Handicap"),
pch = c(21, 21, 21, 24, 22), pt.bg = c(col_sham, col_dull, rep("white", 3)), cex =0.9, bty = "n")
## HB and wing
par(mfrow = c(1, 2))
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Challenge then Dulling", xlim = c(1.5, 2.5), ylim = c(15, 60),
xlab = "", ylab = "Length on Day 12 (mm)")
axis(1, c(-1, 1, 2, 10), c("", "Head + Bill", "Wing", ""))
axis(2, seq(-30, 200, 5), las = 2)
points(hbw2_19$xpos + hbw2_19$p_dodge, hbw2_19$length, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(shbw2_19)){
lines(rep(shbw2_19$xpos[i] + shbw2_19$p_dodge[i], 2),
c(shbw2_19$mu[i] - shbw2_19$se[i], shbw2_19$mu[i] + shbw2_19$se[i]), lwd = 2)
}
points(shbw2_19$xpos + shbw2_19$p_dodge, shbw2_19$mu, pch = shbw2_19$p_shape, cex = 1.4,
bg = as.character(shbw2_19$p_col))
plot(1, 1, type = "n", xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n", bty = "n",
main = "Challenge then Dulling", xlim = c(0.5, 1.5), ylim = c(19, 29),
xlab = "", ylab = "Length on Day 12 (mm)")
axis(1, c(-1, 1, 2, 10), c("", "Head + Bill", "Wing", ""))
axis(2, seq(-30, 200, 2), las = 2)
points(hbw2_19$xpos + hbw2_19$p_dodge, hbw2_19$length, pch = 16, col = alpha("black", 0.3))
#abline(h = 0)
for(i in 1:nrow(shbw2_19)){
lines(rep(shbw2_19$xpos[i] + shbw2_19$p_dodge[i], 2),
c(shbw2_19$mu[i] - shbw2_19$se[i], shbw2_19$mu[i] + shbw2_19$se[i]), lwd = 2)
}
points(shbw2_19$xpos + shbw2_19$p_dodge, shbw2_19$mu, pch = shbw2_19$p_shape, cex = 1.4,
bg = as.character(shbw2_19$p_col))
legend("bottomright", c("No Color", "Dull Color", "Control", "Predator"),
pch = c(21, 21, 21, 24), pt.bg = c(col_sham, col_dull, rep("white", 2)), cex =0.9, bty = "n")
There isn't really anything very obvious going on here. Maybe there is a hint of reduced growth in predation treatments. This is just plotting the average value for each nest, but I really probably should plot these as individual nestling points with nest as a random effect for the day 12 & 15 measures. One thing to consider is that because there seems to be higher nestling attrition in the predation treatments, it may just be that the lightest nestlings in those broods are already gone by the time day 12/15 measurements are taken.
There are some other data types that we have on these birds and could look at or incorporate here in various ways (although there is already a lot in here). Those include:
I'm not sure exactly how this all gets put together into a paper. It seems like the main takeaways are:
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] scales_1.1.1 lme4_1.1-23 Matrix_1.2-18 here_0.1
## [5] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.1 purrr_0.3.4
## [9] readr_1.3.1 tidyr_1.1.1 tibble_3.0.3 ggplot2_3.3.2
## [13] tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.5 lubridate_1.7.9 lattice_0.20-41 assertthat_0.2.1
## [5] rprojroot_1.3-2 digest_0.6.25 plyr_1.8.6 R6_2.4.1
## [9] cellranger_1.1.0 backports_1.1.8 reprex_0.3.0 evaluate_0.14
## [13] httr_1.4.2 pillar_1.4.6 rlang_0.4.7 readxl_1.3.1
## [17] rstudioapi_0.11 minqa_1.2.4 nloptr_1.2.2.2 blob_1.2.1
## [21] rmarkdown_2.3 labeling_0.3 splines_3.6.1 statmod_1.4.34
## [25] munsell_0.5.0 broom_0.7.0 compiler_3.6.1 modelr_0.1.8
## [29] xfun_0.16 pkgconfig_2.0.3 mgcv_1.8-31 htmltools_0.5.0
## [33] tidyselect_1.1.0 fansi_0.4.1 crayon_1.3.4 dbplyr_1.4.4
## [37] withr_2.2.0 MASS_7.3-51.6 grid_3.6.1 nlme_3.1-148
## [41] jsonlite_1.7.0 gtable_0.3.0 lifecycle_0.2.0 DBI_1.1.0
## [45] pacman_0.5.1 magrittr_1.5 cli_2.0.2 stringi_1.4.6
## [49] farver_2.0.3 fs_1.5.0 xml2_1.3.2 ellipsis_0.3.1
## [53] generics_0.0.2 vctrs_0.3.2 boot_1.3-25 tools_3.6.1
## [57] glue_1.4.1 hms_0.5.3 yaml_2.2.1 colorspace_1.4-1
## [61] rvest_0.3.6 knitr_1.29 haven_2.3.1