# install.packages("tidyverse") # install.packages("report") # install.packages("broom") # install.packages("survival") # install.packages("ggfortify") # install.packages("survminer") library(tidyverse) library(report) library(broom) library(survival) library(ggfortify) library(survminer) #### ------------- ## I. Importing the dataset for the glue comparison tests (retention) + explanation of terms #### ------------- #### Importing the data for the comparison of retention rates of tags glued with either e2c or 2oc - this dataset contains the full list of observations since it also includes the males. glue_compare_dataset <- read_delim(file = "data_RFIDmethods_glue_retention_comparison.csv", delim = ";") glimpse(glue_compare_dataset) ## Explanation of terms in "glue_compare_dataset": # location - where the observation is from; the two tanks are K4 and K5. For example, each "K4 fish" is one fish sampled from tank K4, while each "K4 filter" is that day's observation of tank K4's outlet filter where detached salmon lice were found. "K5 bucket" are lice found loose in the anesthetic bucket, and "K5 euthanized" is an observation from a fish euthanized early. # treatment - which glue the female salmon lice in the tank were glued with; 2oc = 2-octyl cyanoacrylate and e2c = ethyl 2-cyanoacrylate. # date - the date of the observation. # days_since_tagging - the number of days since the first observation, when all females in a tank were tagged. # female_chip_yes - the number of females observed with a p-Chip attached (all females on day 0). # female_chip_no - the number of females observed with no p-Chip attached. # male - the number of males observed. #### ------------- ## II. Importing the dataset for the glue comparison tests (mortality) + explanation of terms #### ------------- #### Importing the data for the analysis of mortality over time for tagged females glued with either e2c or 2oc and untagged males in the same tank. glue_mortality_dataset <- read_delim(file = "data_RFIDmethods_glue_mortality_comparison.csv", delim = ";") glimpse(glue_mortality_dataset) ## Explanation of terms in "glue_mortality_dataset": # Group - the combined variables of sex + treatment; mainly for labeling purposes. # Treatment - which glue the female salmon lice in the tank were glued with; 2oc = 2-octyl cyanoacrylate and e2c = ethyl 2-cyanoacrylate. # Sex - female or male. All females were tagged; males were not tagged. # date - the date of the observation. # days_since_tagging - the number of days since the first observation, when all females in a tank were tagged. # days_post_infection - the number of days since the salmon in the tank were infected; reflects the total age of the lice minus their time as free-swimming larvae. # n_alive - the cumulative estimate of living individual lice remaining, adjusted by individuals tallied as dead or missing from a population census. In the e2c treatment, the number of expected live lice is also reduced after day 46 when a host was euthanized early, removing 3 females and 2 males. # n_dead - the cumulative tally of individual lice known to be either dead (found in the outlet filter) or missing from a census observation. The 3 females and 2 males in the e2c treatment removed on day 46 were not included in this tally. #### ------------- ## III. Importing the dataset for the analysis of glue retention time + explanation of terms #### ------------- #### Importing the dataset for the analyses of glue retention time. retention_dataset <- read_delim(file = "data_RFIDmethods_retention_time.csv", delim = ";") glimpse(retention_dataset) ## Explanation of terms in "retention_dataset": # female_ID - the ID number of the female salmon louse. # tag_date - the date on which the female was first tagged. # tag_who - the identity of the person who first applied the RFID tag. # lost_date - the date when the female was recorded as lost. # lost_ID_method - the method used to determine that the female was lost. "chip" = p-Chip was scanned, i.e. the date the louse was found in the tank's outlet filter; "absent" = the first week the female was observed to be missing (if not found with a p-Chip in the outlet filter); "photo" = ID'd from photo from the outlet filter; "killed" = killed during sampling (not included in analysis) # lifespan_days_after_tagging - how many additional days the louse lived after first being tagged (i.e. lost_date - tag_date). # lineage_origin - where the louse's lineage was collected from; different lineages are in separate tanks. # host_density - how many fish are in the tank. "high" = 15 (or 16) fish; "low" = 5 fish. # retagged_tf - whether the louse was ever re-tagged, TRUE (yes) or FALSE (no). # first_retag_cause - the reason the original RFID tag was replaced. "loss" = the tag was missing; "malfunction" = the tag stopped being scannable. # days_death_after_1st_retagging - how many additional days the louse lived after the first time it was re-tagged. # retag1_date : retag4_date - the date(s) a louse was observed to no longer have a (working) RFID tag the first, second, third, and fourth time (if applicable). (NB: this is not necessarily the actual date that the louse WAS re-tagged, since they were not always re-tagged immediately.) # nr_retags - the number of times a louse was re-tagged. # interval_retag0_to_1 : interval_retag3_to_4 - the interval of days until a tag needed to be replaced, from the original to the first, from the first to the second, and so on. Only "interval_retag0_to_1" is used here, i.e. the number of days until the original RFID tag needed to be replaced (retag1_date - tag_date). (Also, see the note for the "retag1_date" regarding subsequent intervals.) #### ------------- ## IV. Importing the dataset looking at impacts on egg extrusion + explanation of terms #### ------------- ## Importing the data looking at impacts on egg extrusion for females first glued while they had 0, 1, or 2 egg strings extruded, with the number of egg strings observed to be extruded in subsequent observations. Excludes females that did not live through at least 2 observations, and those that were re-tagged within 14 days of their original tagging. repro_impacts_dataset <- read_delim(file = "data_RFIDmethods_reproductive_effects.csv", delim = ";") glimpse(repro_impacts_dataset) ## Explanation of terms in "repro_impacts_dataset": # female_ID - the ID number of the female salmon louse. # tag_date - the date on which the female was first glued with an RFID tag. # loss_date - the date when the female was recorded as lost. # lifespan_days_after_tagging - how many additional days the louse lived after first being tagged (loss_date - tag_date) # first_obs_n_eggstrings - the number of egg strings a female was observed to have extruded on the first observation (maximum = 2; one per side); these include "empty" egg strings. # second_obs_n_eggstrings - the number of egg strings a female was observed to have extruded on the second observation. # third_obs_n_eggstrings - the number of egg strings a female was observed to have extruded on the third observation. # max_strings_after - the maximum number of egg strings a female was observed to have extruded in the second or third observation (i.e. after being tagged) # followup_obs_n_eggstrings - for those females that produced fewer egg strings than expected after being tagged, the number of egg strings observed after looking at subsequent observations (when possible). # comment - comments for females that produced fewer egg strings than expects, e.g. the date of the followup observation or whether they instead died shortly after, whether blockages were visible, comments on pre-existing infertility, etc. # first_retag_date - if a louse was re-tagged, the date that a new tag was glued to the louse. (NB: while some lice included here were re-tagged as early as day 17, i.e. on the second observation, but by this time females have reliably extruded their second clutch.) # second_retag_date - the date of the second re-tagging, if this occurred. # time_to_1st_retag - the number of days between the original tagging date and the first re-tagging (if applicable; i.e. 1st_retag_date - tag_date). # time_to_2nd_retag - the number of days between the original tagging date and the second re-tagging (if applicable; i.e. 2nd_retag_date - tag_date). #### ------------- ## V. Importing the dataset for the followup mortality analysis + explanation of terms #### ------------- ## Importing the data for the followup analysis on possible effects on mortality from gluing tags to female salmon lice with 2oc, using lice that were re-tagged 7 days after being glued with their original chip (the "second dose"), along with a cohort of females that were never re-tagged and that were also still alive on day 7. cohort_mortality_dataset <- read_delim(file = "data_RFIDmethods_cohort_mortality.csv", delim = ";") glimpse(cohort_mortality_dataset) # female_ID - the ID number of the female salmon louse. # infect_date - the date when the salmon hosts were infected, resulting in this female # tag_date - the date on which the female was first tagged. # lost_date - the date when the female was recorded as lost (either in the outlet filter or absent during a weekly check). # lifespan_days_post_infection - how many days the female has been alive since first introduced to a host as an infective juvenile. # lifespan_days_after_retag_day - how many additional days the louse lived after either being either re-tagged or not on day 7 (i.e. lost_date - tag_date - 7) # retagged_tf - whether the louse was ever re-tagged, TRUE (yes) or FALSE (no). # nr_retags - how many times the louse was re-tagged. # retag1_date - the date when a louse *was* re-glued with a new RFID tag. (NB: This is in contrast to the dataset in section III, which lists the first date when the RFID tag was missing, but not necessarily replaced.) # retag2_date : retag4_date - the date(s) a louse was *observed* to no longer have a (working) RFID tag the second, third, and fourth time (if applicable). (NB: for these, this is not necessarily the actual date that the louse WAS re-tagged, and this is not used in the analysis.) # time_retag1 : time_retag4 - the interval of days until a tag needed to be replaced, from the original to the first, from the first to the second, and so on. Only "time_retag1" is used here, i.e. the number of days until the original RFID tag needed to be replaced (retag1_date - tag_date). (Also, see the note for the "retag2_date : retag4_date" regarding subsequent intervals.) #### --------------- ## 1.1 - 1.2 Glue comparison tests #### --------------- #### 1.1 - This segment is for testing for a difference in RFID tag retention between the glues "2oc" and "e2c", using a GLM fitted with a binomial distribution where the binomial response variable combines the number of tagged female salmon lice scored as successes and that of untagged females scored as failures. ## Making a table with selected columns: glue_retention_compare <- glue_compare_dataset |> select(treatment, days_since_tagging, female_chip_yes, female_chip_no) ## Setting up the model: glm_glue_retention_compare <- glm(cbind(female_chip_yes, female_chip_no) ~ treatment * days_since_tagging, glue_retention_compare, family = "binomial") ## The model output: summary(glm_glue_retention_compare) ## Diagnostic plots: plot(glm_glue_retention_compare) ## This uses the R package "report" to quickly obtain the beta and 95% confidence intervals for each effect: report(glm_glue_retention_compare) #### 1.2 - This segment is for testing for differences in mortality between tagged female salmon lice glued with either "2oc" or "e2c" as well as untagged males in the same tank, using a GLM fitted with a binomial distribution where the binomial response variable combines the number of living salmon lice scored as successes and that of dead salmon lice scored as failures. ## Setting up the model: glm_glue_mortality_compare <- glm(cbind(n_alive, n_dead) ~ Sex + Treatment + days_since_tagging, glue_mortality_dataset, family = "binomial") ## The model output: summary(glm_glue_mortality_compare) ## Diagnostic plots: plot(glm_glue_mortality_compare) ## This uses the R package "report" to quickly obtain the beta and 95% confidence intervals for each effect: report(glm_glue_mortality_compare) #### --------------- ## 1.3 Glue comparison graphs #### --------------- #### This segment is for generating the graphs used in the article, regarding the glue comparison. ## 1.3.1 - Generating the summary dataset from the main dataset, finding the proportion of tagged females on population census days, and adding a standard error based on proportion, and the minimum and maximum SE: glue_chip_proportion <- glue_compare_dataset |> filter(location == "K4 fish" | location == "K5 fish" | location == "K5 bucket") |> group_by(treatment, days_since_tagging) |> summarize(sum_chip_yes = sum(female_chip_yes), sum_chip_no = sum(female_chip_no)) |> mutate(proportion_chipped = (sum_chip_yes / (sum_chip_yes + sum_chip_no))) |> mutate(SE_proportion_chipped = sqrt((proportion_chipped * (1 - proportion_chipped))/(sum_chip_yes + sum_chip_no))) |> mutate(min_chipped_SE = (proportion_chipped - SE_proportion_chipped), max_chipped_SE = (proportion_chipped + SE_proportion_chipped)) ## 1.3.2 - Visualization of the proportion of females retaining their tag when glued with either 2oc or e2c, during census events: ggplot(glue_chip_proportion, aes(x = days_since_tagging, y = proportion_chipped)) + geom_point(aes(shape = treatment), size = 2) + geom_line(aes(group = treatment)) + geom_errorbar(aes(ymin=min_chipped_SE, ymax=max_chipped_SE), width=0.01) + scale_shape_manual(values = c("2oc" = 16, "e2c" = 1)) + scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) + labs(x = "Time after tagging (days)", y = "Proportion tagged") + theme_classic() ## 1.3.3 - Visualization of the cumulative number of females and males recorded as dead over time; females are tagged with either 2oc or e2c glue, and males are untagged but grouped in under the same treatment as the females in the same tank: ggplot(glue_mortality_dataset, aes(x = days_since_tagging, y = n_dead, shape = Group)) + geom_line(aes(group = Group, linetype = Treatment), show.legend=FALSE) + geom_point(size = 2)+ scale_shape_manual(values = c("2oc female" = 16, "2oc male" = 17, "e2c female" = 1, "e2c male" = 2)) + scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) + labs(x = "Time after tagging (days)", y = "Number recorded dead") + theme_classic() #### --------------- ## 2.1 Glue retention over time - dataframe construction (all observations) #### --------------- #### This segment is for the building the dataframe used in the analysis of glue retention over time, using data from lice tagged as part of another, ongoing research project. #### NB: This entire section should be run prior to sections 2.2 - 2.5 ## 2.1 - This section involves building the dataframe to find the proportion of living females still retaining their original chip. ## This makes sure that all dates are read correctly as dates: retention_time <- retention_dataset |> mutate(tag_date = dmy(tag_date), lost_date = dmy(lost_date), retag1_date = dmy(retag1_date), retag2_date = dmy(retag2_date), retag3_date = dmy(retag3_date), retag4_date = dmy(retag4_date)) ## This converts the post-registration lifespan of females and the time when their first RFID tag was observed to need replacement into weeks instead of days: chiploss_bins <- retention_time |> drop_na(lost_date) |> # drops 1 female that was accidentally killed mutate(week_lost = cut(lifespan_days_after_tagging, breaks=c(-1, 7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84, 91, 98, 105, 112, 119, 126, 133, 140, 147, 154, 161, 168, 175, 182, 189, 196, 203, 210, 217, 224, 231, 238, 245, 252, 259, 266))) |> # converts the lifespan of all females from days into weeks mutate(week_lost = recode(week_lost, "(-1,7]"="1", "(7,14]"="2", "(14,21]"="3", "(21,28]"="4", "(28,35]"="5", "(35,42]"="6", "(42,49]"="7", "(49,56]"="8", "(56,63]"="9", "(63,70]"="10", "(70,77]"="11", "(77,84]"="12", "(84,91]"="13", "(91,98]"="14", "(98,105]"="15", "(105,112]"="16", "(112,119]"="17", "(119,126]"="18", "(126,133]"="19", "(133,140]"="20", "(140,147]"="21", "(147,154]"="22", "(154,161]"="23", "(161,168]"="24", "(168,175]"="25", "(175,182]"="26", "(182,189]"="27", "(210,217]"="31", "(217,224]"="32", "(224,231]"="33", "(238,245]"="35", "(252,259]"="37", "(259,266]"="38")) |> # re-writes "weeks_lost" into the number of weeks instead of the range of days mutate(week_retagged = cut(interval_retag0_to_1, breaks=c(-1, 7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84, 91, 98, 105, 112, 119, 126, 133, 140, 147, 154, 161, 168, 175, 182, 189, 196, 203, 210, 217, 224, 231, 238, 245, 252, 259, 266)))|> # converts the time until the original RFID tag was lost from days into weeks mutate(week_retagged = recode(week_retagged, "(-1,7]"="1", "(7,14]"="2", "(14,21]"="3", "(21,28]"="4", "(28,35]"="5", "(35,42]"="6", "(42,49]"="7", "(49,56]"="8", "(56,63]"="9", "(63,70]"="10", "(70,77]"="11", "(84,91]"="13", "(91,98]"="14", "(105,112]"="16", "(112,119]"="17", "(126,133]"="19", "(133,140]"="20", "(147,154]"="22")) # re-writes "week_retagged" into the number of weeks instead of the range of days ## This creates a cumulative total of how many females died per week: chiploss_bins_total <- chiploss_bins |> group_by(week_lost) |> summarise(dead = n()) |> # how many females died each week mutate(week_lost = as.numeric(week_lost))|> add_row(week_lost = 0, dead = 0)|> # adding missing observations (when no females died) add_row(week_lost = 28, dead = 0)|> add_row(week_lost = 29, dead = 0)|> add_row(week_lost = 30, dead = 0)|> add_row(week_lost = 34, dead = 0)|> add_row(week_lost = 36, dead = 0)|> arrange(week_lost)|> # makes sure the weeks are arranged consecutively mutate(cumulative_dead = cumsum(dead)) |> # converts to a cumulative total over time mutate(total_living = 948 - cumulative_dead) # finds the number of living lice per week (based on the starting sample size of 948) ## Creating a table with the number of lice per week that are alive and no longer retaining their original tag, obtained by filtering + summarizing for weeks 0 to 38: week_lost <- 0:38 n_live_retagged <- c(0, 74, 98, 105, 98, 98, 85, 80, 75, 73, 73, 71, 59, 56, 51, 35, 30, 24, 23, 19, 16, 14, 12, 11, 9, 5, 5, 3, 3, 3, 3, 2, 0, 0, 0, 0, 0, 0, 0) chiploss_retagged_week <- tibble(week_lost, n_live_retagged) ## Joining the two tables together so that dataframe has both the number of females that are still alive each week as well as how many of those females have lost their original tag: chiploss_prop_join <- left_join(chiploss_bins_total, chiploss_retagged_week, by = "week_lost") |> rename(week = week_lost, living_lost = n_live_retagged) ## The final working version of the dataset for tag retention, finding the proportion of living female lice still retaining their original RFID tag within each week: chiploss_proportion <- chiploss_prop_join |> mutate(days = week*7) |> # converting the weeks back to days mutate(living_retained = total_living - living_lost) |> # finding the number of living females that retained their original RFID tag mutate(prop_retained = living_retained/total_living) |> # the proportion of living females retaining their original RFID tag #### --------------- ## 2.2 Glue retention over time - bar graph (all observations) #### --------------- #### This section creates a bar graph showing status of living female lice as either retaining their original tag or having lost it, over each week. #### NB: Requires section 2.1 to be run first. ## This pivots the data into a long format, so that the number of females retaining or having lost their tag on a week are unique observations: chiploss_status_long <- chiploss_proportion |> select(days, living_retained, living_lost) |> rename(original = living_retained, retagged = living_lost) |> pivot_longer(cols = original:retagged, names_to = "status", values_to = "count") # pivoting the table; female counts have the status of either "original" or "retagged" ## This generates the bar graph used in the supplementary materials: ggplot(chiploss_status_long, aes(x = days, y = count, fill = status)) + geom_bar(stat="identity") + labs(x = "Time after chipping (days)", y = "Number of females") + scale_fill_grey(start = 0.7, end = 0.2, name = "Original tag status", labels = c("Retained", "Lost")) + scale_y_continuous(breaks = scales::pretty_breaks(n = 10))+ scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+ theme_classic()+ theme(legend.position = c(x = 0.99, y = 0.70), legend.justification = c(x = "right", y = "bottom")) #### --------------- ## 2.3 Glue retention over time - number of replacements and rate of loss (all observations) #### --------------- #### This section is for analyses of glue retention for observations between the dates of 11th January 2021 to 16th December 2022. #### NB: Requires section 2.1 to be run first. #### 2.3.1 - Basic descriptive statistics: ## The range of days until death (after initial tagging) is 0 to 263; the oldest replaced tag was at 150 days. chiploss_bins |> summarise(maxlife = max(lifespan_days_after_tagging, na.rm=TRUE), minlife = min(lifespan_days_after_tagging, na.rm=TRUE), maxchip = max(interval_retag0_to_1, na.rm=TRUE)) ## Finding the sample size: chiploss_bins |> summarise(n = n()) # n = 948 ## The number of re-tagged lice vs. lice that retained their original tags until death: chiploss_bins |> group_by(retagged_tf) |> summarise(n = n()) # n re-tagged = 216; n retained = 732 ## Summarizing how many tags were replaced due to loss vs. malfunction: chiploss_bins |> filter(retagged_tf == "TRUE") |> group_by(first_retag_cause) |> summarise(n = n()) # n lost = 212; n malfunctioned = 4 #### 2.3.2 - This section uses the proportion of living female salmon lice retaining their original RFID tag over time found in section 2.1 to estimate a rate of loss. (Note: it is of course possible to run a GLM with a binomial distribution as in section 1.1, preserving sample size, but the linear regression is better for expressing an easily-understood rate.) ## Filtering the dataset to weeks when there are more than 20 living individuals (below this level, the proportion varies wildly due to small sample size): chiprate_line <- chiploss_proportion |> filter(total_living > 20) ## *** This is the linear regression showing the rate of loss over time for all initial RFID tags glued in the period between 11th January 2021 to 16th December 2022: lm_retention_full <- lm(prop_retained ~ days, data = chiprate_line) summary(lm_retention_full) # model output plot(lm_retention_full) ## This uses the R package "report" to obtain the beta value: report(lm_retention_full) #### 2.3.3 - Graph showing the rate of RFID tag loss over time (as a proportion) with a fitted line, using the linear model from section 2.3.2. ## Creating predicted values: preds_lm_retention_full <- broom::augment(lm_retention_full, interval = "confidence", conf.level = 0.95) ## Scatter plot with the linear regression from 2.3.2: ggplot(chiprate_line, aes(x = days, y = prop_retained)) + geom_point() + geom_line(aes(y = .fitted), data = preds_lm_retention_full) + labs(x = "Time after tagging (days)", y = "Proportion retaining original RFID tag") + scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+ scale_y_continuous(breaks = scales::pretty_breaks(n = 5))+ theme_classic() #### --------------- ## 2.4 Glue retention over time - divided by year (batch) - building the dataframes #### --------------- #### This segment is for the building the dataframe used to compare the rate of tag loss between the years 2021 (two vials of glue) and 2022 (one vial of glue), using data from lice tagged as part of another, ongoing research project. (Note: "2022" is an estimate; the third vial was likely opened at the end of January or the start of February, but this is off by at most 33 lice out of 948.) #### NB: Requires section 2.1 to be run first. #### NB: This entire section should be run prior to section 2.5. ## First building the dataframe for 2021: chiploss2021_bins <- chiploss_bins |> filter(tag_date < dmy("01.01.2022")) |> # filtering for lice originally chipped in 2021 group_by(week_lost) |> summarise(dead_2021 = n()) |> mutate(week_lost = as.numeric(week_lost))|> add_row(week_lost = 0, dead_2021 = 0)|> # adding observations for weeks when no females died add_row(week_lost = 28, dead_2021 = 0)|> add_row(week_lost = 29, dead_2021 = 0)|> add_row(week_lost = 30, dead_2021 = 0)|> add_row(week_lost = 34, dead_2021 = 0)|> add_row(week_lost = 36, dead_2021 = 0)|> arrange(week_lost)|> # makes sure the weeks are arranged consecutively mutate(cumulative_dead_2021 = cumsum(dead_2021)) |> # converts to a cumulative total over time mutate(total_living_2021 = 599 - cumulative_dead_2021) # finds the number of living lice per week (based on the starting sample size of 599) ## Now building the dataframe for 2022: chiploss2022_bins <- chiploss_bins |> filter(tag_date > dmy("31.12.2021")) |> # filtering for lice originally chipped in 2021 group_by(week_lost) |> summarise(dead_2022 = n()) |> mutate(week_lost = as.numeric(week_lost))|> add_row(week_lost = 0, dead_2022 = 0)|> add_row(week_lost = 21:26, dead_2022 = 0)|> add_row(week_lost = 28:38, dead_2022 = 0)|> arrange(week_lost)|> # makes sure the weeks are arranged consecutively mutate(cumulative_dead_2022 = cumsum(dead_2022)) |> mutate(total_living_2022 = 349 - cumulative_dead_2022) # finds the number of living lice per week (based on the starting sample size of 349) ## Creating a table with the number of lice per week that are alive and no longer retaining their original tag, split between those originally tagged in 2021 and those originally tagged in 2022. Obtained by filtering + summarizing for weeks 0 to 38: week_lost <- 0:38 lost_2021 <- c(0, 29, 45, 53, 50, 51, 40, 38, 36, 36, 35, 35, 29, 28, 28, 19, 18, 19, 19, 16, 15, 13, 11, 10, 8, 4, 4, 3, 3, 3, 3, 2, 0, 0, 0, 0, 0, 0, 0) lost_2022 <- c(0, 45, 53, 52, 48, 47, 45, 42, 39, 37, 38, 36, 30, 28, 23, 16, 12, 5, 4, 3, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) chiploss_compare_week <- tibble(week_lost, lost_2021, lost_2022) ## Joining together the tables for the cumulative tallies of living and dead lice for both years: chiploss_2021_2021_total <- left_join(chiploss2021_bins, chiploss2022_bins, by = "week_lost") ## Joining the tables so that the dataframe has both the number of females that are still alive each week as well as how many of those females have lost their original tag: chiploss_year_comparison <- left_join(chiploss_2021_2021_total, chiploss_compare_week, by = "week_lost")|> rename(week = week_lost) |> mutate(days = week*7) |> # converting the weeks back to days mutate(retained_2021 = total_living_2021 - lost_2021, retained_2022 = total_living_2022 - lost_2022) |> # finding the number of living females that retained their original RFID tag mutate(prop_retained_2021 = retained_2021/total_living_2021, prop_retained_2022 = retained_2022/total_living_2022) # the proportion of living females retaining their original RFID tag ## Adding the values from both years combined (from section 2.1) to the same dataframe to directly compare the loss rates for 2021 and 2022 against that of the entire dataset: all_chiploss_comparison <- left_join(chiploss_year_comparison, chiploss_prop_join, by = "week") |> rename(dead_total = dead, cumulative_dead_total = cumulative_dead, total_living_total = total_living, lost_total = living_lost) |> mutate(retained_total = total_living_total - lost_total) |> mutate(prop_retained_total = retained_total/total_living_total) ## Cleaning the dataframes up a little to pivot them to be longer, splitting them into one line per observation by year (2021, 2021, or total (2021 & 2022)), and joining them together with the number of living females next to the proportion of those females still retaining their original p-Chip. chiprate_line_compare_1 <- all_chiploss_comparison |> select(days, prop_retained_2021, prop_retained_2022, prop_retained_total) |> # Selecting values for the proportion of retained tags rename("2021" = prop_retained_2021, "2022" = prop_retained_2022, Total = prop_retained_total) |> # Renaming the columns to variable names pivot_longer(cols = "2021":Total, names_to = "Year", values_to = "proportion") |> # Pivoting drop_na(proportion) chiprate_line_compare_2 <- all_chiploss_comparison |> select(days, total_living_2021, total_living_2022, total_living_total) |> # Selecting values for the number of living lice in a 7-day period rename("2021" = total_living_2021, "2022" = total_living_2022, Total = total_living_total) |> # Renaming the columns to variable names pivot_longer(cols = "2021":Total, names_to = "Year", values_to = "pop_size") # Pivoting all_chiprate_line <- left_join(chiprate_line_compare_1, chiprate_line_compare_2, by = c("days", "Year")) |> # Joining the two pivoted tables mutate(Year = as.factor(Year)) #### --------------- ## 2.5 Glue retention over time - divided by year (batch) - number of replacements and rate of loss #### --------------- #### This section is for analyses of glue retention for observations between the dates of 11th January 2021 to 16th December 2022, comparing the rate of tag loss between the years 2021 (two vials of glue) and 2022 (one vial of glue). This is to test for differences in effectiveness between batches. #### NB: Requires sections 2.1 and 2.4 to be run first. ## Section 2.5.1 - Basic descriptive statistics: ## Finding the sample size by year: chiploss_bins |> mutate(tag_year = ifelse(tag_date > dmy("31.12.2021"), 2022, 2021)) |> group_by(tag_year) |> summarise(n = n()) # n = 599 in 2021; n = 349 in 2022 ## Finding the number of lice that either retained their original RFID tag their entire lives or were re-tagged in each year: chiploss_bins |> mutate(tag_year = ifelse(tag_date > dmy("31.12.2021"), 2022, 2021)) |> group_by(tag_year, retagged_tf) |> summarise(n = n()) # n re-tagged = 117 in 2021; n re-tagged = 99 in 2022 #### 2.5.2 - This section uses the proportion of living female salmon lice retaining their original RFID tag over time found in sections 2.1 and 2.4 to estimate a rate of loss for lice originally tagged in 2021 vs. 2022, as well as the combined total for both years. ## Filtering the data to weeks when there are more than 20 living individuals (below this level, the proportion varies wildly due to small sample size): all_chiprate_line_compare <- all_chiprate_line |> filter(pop_size > 20) ## Filtering the data to exclude the values for both years combined ("Total", i.e. 2021 & 2022), leaving only values for the years 2021 and 2022 separate: year_chiprate_line <- all_chiprate_line_compare |> filter(Year != "Total") ## *** This is the linear regression comparing the rate of RFID tag loss over time between females initially tagged in 2021 (two vials of glue) vs. those initially tagged in 2022 (a third vial of glue), testing for batch variation: lm_year_compare <- lm(proportion ~ days * Year, data = year_chiprate_line) summary(lm_year_compare) # model output plot(lm_year_compare) ## This uses the R package "report" to obtain the beta values: report(lm_year_compare) #### 2.5.3 - Graph showing the rate of RFID tag loss over time (as a proportion) based lice originally tagged in 2021 or 2022, compared to the rate of loss for both years combined. ## Running a linear regression to compare the rate of loss for both years combined to the rate of loss for either 2021 or 2022: lm_lossrate_compare <- lm(proportion ~ days * Year, data = all_chiprate_line_compare) summary(lm_lossrate_compare) ## Creating predicted values: preds_lm_lossrate_compare <- broom::augment(lm_lossrate_compare, interval = "confidence", conf.level = 0.95) ## Visualization of the rate of loss over time ("Total" is renamed here to "2021 + 2022"): ggplot(all_chiprate_line_compare, aes(x = days, y = proportion, shape = Year)) + geom_point() + geom_line(aes(y = .fitted, linetype = Year), data = preds_lm_lossrate_compare) + labs(x = "Time after tagging (days)", y = "Proportion retaining original tag") + scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+ scale_y_continuous(breaks = scales::pretty_breaks(n = 10))+ scale_shape_manual(values = c("2021" = 2, "2022" = 1, "Total" = 15), labels = c("2021", "2022", "2021 + 2022"))+ scale_linetype_manual(values = c("2021" = "dashed", "2022" = "dashed", "Total" = "solid"), labels = c("2021", "2022", "2021 + 2022"))+ theme_classic()+ theme(legend.position = c(x = 0.99, y = 0.70), legend.justification = c(x = "right", y = "bottom")) #### --------------- ## 3.1 Effects of tagging on reproduction #### --------------- #### This section checks for differences in the numbers of egg strings (max = 2) that female salmon lice extruded after being glued with an RFID tag, shortly before or after their first clutch. ## Finding the difference in the number of egg strings extruded after being glued: repro_impacts <- repro_impacts_dataset |> mutate(change_n_eggstrings = (max_strings_after - first_obs_n_eggstrings)) ## Total sample size: repro_impacts |> summarise(n = n()) # n = 233 ## Sample size of each starting value (i.e. the number of egg strings first observed): repro_impacts |> group_by(first_obs_n_eggstrings) |> summarise(n = n()) # n started with 0 = 55; n started with 1 = 6; n started with 2 = 172 ## The number of lice in each starting group that produced 2 egg strings in the observation period (i.e. the ideal number): repro_impacts |> filter(max_strings_after == 2) |> group_by(first_obs_n_eggstrings) |> summarise(n = n()) # 51 of those that began with 0, 4 of those than began with 1, and 164 of those that began with 2 ## The change in number of egg strings in each group (based on starting value), where they did NOT produce 2 egg strings in the observation period: repro_impacts |> filter(max_strings_after != 2) |> group_by(first_obs_n_eggstrings, change_n_eggstrings) |> summarise(n = n()) # 2 that began with 0 extruded 0 (change_n_eggstrings = 0) # 2 that began with 0 extruded 1 (change_n_eggstrings = 1) # 2 that began with 1 extruded 1 (change_n_eggstrings = 0) # 2 that began with 2 extruded 0 (change_n_eggstrings = -2) # 6 that began with 2 extruded 1 (change_n_eggstrings = -1) #### --------------- ## 4.1 Effects of tagging on mortality #### --------------- #### This section checks for increased mortality as an effect of being glued with 2oc by comparing lice that were glued twice within a short period of time (i.e. a "second dose" 7 days after the first) to a manually-selected cohort of lice that were initially tagged on the same day, but which were never re-tagged. ## Sample size: cohort_mortality_dataset|> group_by(retagged_tf) |> summarise(n = n()) # n re-tagged = 49; n not re-tagged = 270; total n = 319 #### 4.1.2 - Survival analysis using a Cox Proportional Hazard model on the remaining lifespan for females either re-tagged 7 days after being glued with their first p-Chip, compared to a cohort that was never re-tagged. ## Selecting the data: cohort_mortality <- cohort_mortality_dataset |> select(lifespan_days_after_retag_day, retagged_tf) |> mutate(status = 2) # to indicate that they are dead for the survival analysis ## The Cox Proportional Hazards model: cox_model_cohort <- coxph(Surv(lifespan_days_after_retag_day, status) ~ retagged_tf, data = cohort_mortality) summary(cox_model_cohort) # Model output; no significant difference (p = 0.693) ## Model diagnostics: ## Testing the proportional hazards assumption of the model (not statistically significant): test.ph = cox.zph(cox_model_cohort) test.ph ## Diagnostic plots (requires the R package "survminer"): ggcoxzph(test.ph) # Schoenfeld visualization; solid line is horizontal ggcoxdiagnostics(cox_model_cohort, type = "dfbeta", linear.predictions = TRUE) # One observation seems to have a high weight ggcoxdiagnostics(cox_model_cohort, type = "deviance", linear.predictions = TRUE) # Deviance residuals appear roughly symmetrical #### 4.1.3 - Fitting the survival data in order to visualize it as a Kaplan-Meier plot: km <- with(cohort_mortality, Surv(lifespan_days_after_retag_day, status)) head(km,80) km_fit <- survfit(Surv(lifespan_days_after_retag_day, status) ~ 1, data=cohort_mortality) summary(km_fit, times = c(1,10*(1:14))) # fitting the data into 14 bins in intervals of 10 days, up to days 140-150 (since the maximum lifespan was 145 days) km_chip_fit <- survfit(Surv(lifespan_days_after_retag_day, status) ~ retagged_tf, data=cohort_mortality) ## Visualization of the model (requires the R package "ggfortify"): autoplot(km_chip_fit) #### 4.1.4 - Chi-square test for the number of females that were either alive or dead after 14 days (i.e. (lost_date - tag_date - 7) > 14 days) according to whether they were re-tagged or not. ## Finding the number of females that were *still alive* after day 14, after either being re-tagged or not (i.e. > 21 days after their initial tagging): cohort_mortality_dataset |> filter(lifespan_days_after_retag_day > 14) |> group_by(retagged_tf) |> summarize(n_alive = n()) # n not re-tagged = 175; n re-tagged = 32 ## Finding the number of females that were *dead* after day 14, after either being re-tagged or not (i.e. > 21 days after their initial tagging): cohort_mortality_dataset |> filter(lifespan_days_after_retag_day < 15) |> group_by(retagged_tf) |> summarize(n_dead = n()) # n not re-tagged = 95; n re-tagged = 17 ## Creating a table of the values found above: retagged_chisq <- as.table(rbind(c(32, 17), c(175, 95))) dimnames(retagged_chisq) <- list("re-tagged" = c("yes", "no"), "status" = c("alive","dead")) print(retagged_chisq) ## Running the chi-square test (highly non-significant): chisq.test(retagged_chisq)