diff --git a/bootstrap/initial/data/config.JSON b/bootstrap/initial/data/config.JSON index 24d763b..e741855 100644 --- a/bootstrap/initial/data/config.JSON +++ b/bootstrap/initial/data/config.JSON @@ -7,6 +7,6 @@ "report_name": "SmartDots_Report_Event_398", "report_title": "SmartDots Report for event 398", "report_tokens": "tokens goes here", - "mode_definition": "standard", + "mode_definition": "multistage", "strata": ["strata"] } diff --git a/bootstrap/smartdots_db.R b/bootstrap/smartdots_db.R index a0d6714..126083e 100644 --- a/bootstrap/smartdots_db.R +++ b/bootstrap/smartdots_db.R @@ -36,9 +36,9 @@ unlink(zipfile) ad <- ad[,names(ad) != "Comment"] # set experise of event organiser to 1 -if (all(is.na(ad$expertise[ad$TypeAnnotation == "eventOrganizer"]))) { +if (all(is.na(ad$expertise[ad$TypeAnnotation == "Organizer" |ad$TypeAnnotation == "eventDelegate" |ad$TypeAnnotation == "eventOrganizer"|ad$TypeAnnotation == "Delegate"]))) { message("setting eventOrganiser expertise to 'expert' as none was provided.") - ad$expertise[ad$TypeAnnotation == "eventOrganizer"] <- 1 + ad$expertise[ad$TypeAnnotation == "Organizer"|ad$TypeAnnotation == "eventDelegate" |ad$TypeAnnotation == "eventOrganizer"|ad$TypeAnnotation == "Delegate"] <- 1 } # Change the Strata variable in the ad database if needed diff --git a/data_checker.R b/data_checker.R index 3af8563..3af7b1d 100644 --- a/data_checker.R +++ b/data_checker.R @@ -1,97 +1,102 @@ -## perform some sanity checks on the data - -library(icesTAF) -library(jsonlite) -library(tidyr) - -# # load configuration -config <- read_json("bootstrap/data/config.json", simplifyVector = TRUE) - -# get data from bootstrap folder ------------------------------- -ad <- read.taf("bootstrap/data/smartdots_db/ad.csv") - -# tag some feilds as missing? - -# some messages to the user ------ -frmt_vector <- function(x) { - namesx <- names(x) - namesx[namesx == ""] <- "" - paste(paste0(namesx, ": ", x), collapse = ", ") -} - -check_ad <- function(ad, what = "ad") { - checks <- - list( - c("Summary of ", what), - c("number of annotations: ", nrow(ad)), - c("samples with missing area: ", sum(ad$ices_area == "")), - c("samples with missing stock: ", sum(is.na(ad$stock) | ad$stock == "")), - c("samples with missing prep_method: ", sum(is.na(ad$prep_method) | ad$prep_method == "")), - c("prep_method names: ", frmt_vector(table(ad$prep_method))), - c("Advanced reader annotations: ", sum(ad$expertise)), - c("Samples with missing strata: ", sum(is.na(ad$strata))) - ) - - check_text <- paste(sapply(checks, paste, collapse = ""), collapse = "\n * ") - - # other checks - multiple_annotations <- - ad %>% - dplyr::group_by(EventID, event_name, ices_area, FishID, reader) %>% - dplyr::count() %>% - dplyr::filter(n > 1) %>% - dplyr::rename(annotations = n) - - if (nrow(multiple_annotations) > 0) { - txt <- paste(capture.output(print(multiple_annotations)), collapse = "\n") - image_urls <- - sprintf( - "https://smartdots.ices.dk/manage/viewDetailsImage?tblEventID=%i&SmartImageID=%i", - multiple_annotations$EventID, - multiple_annotations$FishID) - - check_text <- - paste0(check_text, - "\n\n*****************\n", - "**** Warning ****\n", - "*****************\n\n", - "Some readers have multiple annotations:\n\n", - txt, - "\n\nSee annotated images here:\n\t", - paste(image_urls, collapse = "\n\t") - ) - - } - - - if (sum(ad$expertise) == 0) { - check_text <- - paste0(check_text, - "\n\n*****************\n", - "**** Warning ****\n", - "*****************\n\n", - "** There are no advanced readers! **\n", - "** the report scripts require there to be advanced readers. **" - ) - - } - - - msg(check_text, "\n") -} - - - -if (config$onlyApproved == FALSE) { - # check all data - msg("Checking ALL data for Event: ", config$event_id) - - check_ad(ad, "ALL (approved and unapproved) annotations (sets of dots)") -} - -msg("Checking approved data for Event: ", config$event_id) - -check_ad(ad, "approved annotations (sets of dots)") - - -# done +## perform some sanity checks on the data + + +library(icesTAF) +library(jsonlite) +library(tidyr) + +# # load configuration +#config <- read_json("bootstrap/data/config.json", simplifyVector = TRUE) +config <- read_json("bootstrap/initial/data/config.json", simplifyVector = TRUE) + +# get data from bootstrap folder ------------------------------- +#ad <- read.taf("bootstrap/smartdots_db/ad.csv") +ad <- read.taf("bootstrap/data/smartdots_db/ad.csv") +#ad <- read.taf("bootstrap/ad.csv") + +# tag some feilds as missing? + +# some messages to the user ------ +frmt_vector <- function(x) { + namesx <- names(x) + namesx[namesx == ""] <- "" + paste(paste0(namesx, ": ", x), collapse = ", ") +} + +check_ad <- function(ad, what = "ad") { + checks <- + list( + c("Summary of ", what), + c("number of annotations: ", nrow(ad)), + c("samples with missing area: ", sum(ad$ices_area == "")), + c("samples with missing stock: ", sum(is.na(ad$stock) | ad$stock == "")), + c("samples with missing prep_method: ", sum(is.na(ad$prep_method) | ad$prep_method == "")), + c("prep_method names: ", frmt_vector(table(ad$prep_method))), + c("Advanced reader annotations: ", sum(ad$expertise)), + c("Samples with missing strata: ", sum(is.na(ad$strata))) + ) + + check_text <- paste(sapply(checks, paste, collapse = ""), collapse = "\n * ") + + # other checks + multiple_annotations <- + ad %>% + dplyr::group_by(EventID, event_name, ices_area, FishID, reader) %>% + dplyr::count() %>% + dplyr::filter(n > 1) %>% + dplyr::rename(annotations = n) + + if (nrow(multiple_annotations) > 0) { + txt <- paste(capture.output(print(multiple_annotations)), collapse = "\n") + image_urls <- + sprintf( + "https://smartdots.ices.dk/manage/viewDetailsImage?tblEventID=%i&SmartImageID=%i", + multiple_annotations$EventID, + multiple_annotations$FishID) + + check_text <- + paste0(check_text, + "\n\n*****************\n", + "**** Warning ****\n", + "*****************\n\n", + "Some readers have multiple annotations:\n\n", + txt, + "\n\nSee annotated images here:\n\t", + paste(image_urls, collapse = "\n\t") + ) + + } + + + if (sum(ad$expertise) == 0) { + check_text <- + paste0(check_text, + "\n\n*****************\n", + "**** Warning ****\n", + "*****************\n\n", + "** There are no advanced readers! **\n", + "** the report scripts require there to be advanced readers. **" + ) + + } + + + msg(check_text, "\n") +} + + + +if (config$onlyApproved == FALSE) { + # check all data + msg("Checking ALL data for Event: ", config$event_id) + + check_ad(ad, "ALL (approved and unapproved) annotations (sets of dots)") +} + +msg("Checking approved data for Event: ", config$event_id) + +check_ad(ad, "approved annotations (sets of dots)") + + +# done + diff --git a/data_processing.R b/data_processing.R index 02a87dd..76c6676 100644 --- a/data_processing.R +++ b/data_processing.R @@ -1,157 +1,174 @@ -## Preprocess data, write TAF data tables - -## Before: -## After: report_template.docx, dist.csv, -## ad_long.csv, ad_long_ex.csv, -## ad_wide.csv, ad_wide_ex.csv - -library(icesTAF) -library(jsonlite) -library(tidyr) -library(lubridate) -library(plyr) -library(dplyr) -library(tidyverse) -taf.library(ragree) - -# create data directory -mkdir("data") - -# get utility functions -source("utilities.R") -source("utilities_data.R") - -# load configuration -config <- read_json("bootstrap/data/config.json", simplifyVector = TRUE) - -# get data from bootstrap folder ------------------------------- -ad <- read.taf("bootstrap/data/smartdots_db/ad.csv") - -# prepare data ------------------------------- - -# keep only approved annotations -if (config$onlyApproved) { - ad <- ad[ad$IsApproved == 1, ] -} - -# add date columns -ad <- - within(ad, { - year <- year(parse_date_time(catch_date, "%d/%m/%Y %H:%M:%S")) - qtr <- quarter(parse_date_time(catch_date, "%d/%m/%Y %H:%M:%S")) - month <- month(parse_date_time(catch_date, "%d/%m/%Y %H:%M:%S")) - }) - -# if variables are missing add "-" -ad$ices_area[is.na(ad$ices_area) | ad$ices_area == ""] <- "-" -ad$stock[is.na(ad$stock) | ad$stock == ""] <- "-" -ad$prep_method[is.na(ad$prep_method) | ad$prep_method == ""] <- "-" -ad$Sex[is.na(ad$Sex) | ad$Sex==""]<-"NI" -ad<-ad[ad$Sex!="NI",] - -# if no advanced readers! make them all advanced -if (all(ad$expertise == 0)) { - msg("NOTE: all readers were Basic - all have been converted to Advanced") - ad$expertise[] <- 1 -} - -# convert reader expertise -ad$expertise <- c("Basic", "Advanced")[ad$expertise + 1] - -# Assign weight to the readers based in their ranking-experience -# Add different weight columns to data -weight <- length(sort(unique(ad$reader_number))):1 -reader_number <- sort(unique(ad$reader_number)) -reader <- data.frame(reader_number = reader_number, weight_I = weight, weight_II = 1 / (1 + log(sort(weight, decreasing = FALSE) + 0.0000000001))) -ad <- merge(ad, reader, by.x = "reader_number", by.y = "reader_number", all.x = TRUE) - -ad4webgr <- ad - -# Before calculating the mode, give to the sampleID of readings by eventOrganizer (histological sample) the same name as the samples analyzed by the other readers, so the maturity defined for the histological samples is assigned as mode to the other samples of the same FishID -fishid <- sort(unique(ad$FishID)) -for (i in 1:length(fishid)) -{ - nohist <- ad[ad$FishID == fishid[i] & ad$TypeAnnotation != "eventOrganizer", ] - sampleid_nohist <- unique(nohist$SampleID) - yeshist <- ad[ad$FishID == fishid[i] & ad$TypeAnnotation == "eventOrganizer", ] - if (dim(yeshist)[1] > 0) { - yeshist <- yeshist[rep(row.names(yeshist), length(sampleid_nohist)), ] - yeshist$SampleID <- sampleid_nohist - } - - temp <- rbind(nohist, yeshist) - - if (i == 1) { - result <- temp - } else { - result <- rbind(result, temp) - } -} - -ad <- result - -# Calculate modal maturity stage and coefficient of unalikability of maturity stage -ad_long <- ad %>% - add_modal_trad(varmod = "Maturity", config$ma_method) %>% - add_modal_linearweight(varmod = "Maturity", config$ma_method) %>% - add_modal_negexpweight(varmod = "Maturity", config$ma_method) - -ad_long_ex <- ad[ad$expertise == "Advanced", ] %>% - add_modal_trad(varmod = "Maturity", config$ma_method) %>% - add_modal_linearweight(varmod = "Maturity", config$ma_method) %>% - add_modal_negexpweight(varmod = "Maturity", config$ma_method) - - -# Calculate modal maturity stage and coefficient of unalikability of sex category -ad_long <- ad_long %>% - add_modal_trad(varmod = "Sex", config$ma_method) %>% - add_modal_linearweight(varmod = "Sex", config$ma_method) %>% - add_modal_negexpweight(varmod = "Sex", config$ma_method) - -ad_long_ex <- ad_long_ex %>% - add_modal_trad(varmod = "Sex", config$ma_method) %>% - add_modal_linearweight(varmod = "Sex", config$ma_method) %>% - add_modal_negexpweight(varmod = "Sex", config$ma_method) - -# Choose the final mode (traditional, readers linear weight or negative exponential linear weight) based in the existence of histological samples or not, and, in case there are no histological samples, depending if there is multimodality or not. -ad_long <- select_mode(ad_long, config$ma_method, config$mode_definition) -ad_long_ex <- select_mode(ad_long_ex, config$ma_method, config$mode_definition) - -# prepare data in wbgr output format -# IMAGE,1,2,3,4,5,6,7,8,9,10,11,12,13 -# Expertise level,-,-,-,-,-,-,-,-,-,-,-,-,- -# Stock assessment,no,no,no,no,no,no,no,no,no,no,no,no,no -# 6698256.jpg,1,1,1,1,1,-,2,1,-,2,-,1,- -# 6698257.jpg,3,3,3,3,2,1,3,3,-,3,-,3,- -readers <- sort(unique(ad$reader)) - -webgr_maturity <- spread(ad4webgr[c("FishID", "reader", "Maturity")], key = reader, value = Maturity) -webgr_maturity[] <- paste(unlist(webgr_maturity)) -webgr_maturity[webgr_maturity == "NA"] <- "-" -webgr_maturity <- - rbind( - c("Expertise level", rep("-", length(readers))), - c("Stock assessment", rep("no", length(readers))), - webgr_maturity - ) -names(webgr_maturity) <- c("IMAGE", 1:length(readers)) -head(webgr_maturity) - -webgr_sex <- spread(ad4webgr[c("FishID", "reader", "Sex")], key = reader, value = Sex) -webgr_sex[] <- paste(unlist(webgr_sex)) -webgr_sex[webgr_sex == "NA"] <- "-" -webgr_sex <- - rbind( - c("Expertise level", rep("-", length(readers))), - c("Stock assessment", rep("no", length(readers))), - webgr_sex - ) -names(webgr_sex) <- c("IMAGE", 1:length(readers)) -head(webgr_sex) - -# write out input data tables for use later -write.taf(ad4webgr, "data/data.csv", quote = TRUE) -write.taf(ad_long, "data/ad_long.csv", quote = TRUE) -write.taf(ad_long_ex, "data/ad_long_ex.csv", quote = TRUE) -write.taf(webgr_maturity, "data/WebGR_maturity_ages_all.csv", quote = TRUE) -write.taf(webgr_sex, "data/WebGR_sex_ages_all.csv", quote = TRUE) +## Preprocess data, write TAF data tables + +## Before: +## After: report_template.docx, dist.csv, +## ad_long.csv, ad_long_ex.csv, +## ad_wide.csv, ad_wide_ex.csv + +library(icesTAF) +library(jsonlite) +library(tidyr) +library(lubridate) +library(plyr) +library(dplyr) +library(tidyverse) +taf.library(ragree) +#devtools::install_github("raredd/ragree") +#library(ragree) + +# create data directory +#mkdir("data") + +# get utility functions +source("utilities.R") +source("utilities_data.R") + +# load configuration +#config <- read_json("bootstrap/data/config.json", simplifyVector = TRUE) +config <- read_json("bootstrap/initial/data/config.json", simplifyVector = TRUE) + + +# get data from bootstrap folder ------------------------------- +ad <- read.taf("bootstrap/data/smartdots_db/ad.csv") +#ad <- read.taf("bootstrap/smartdots_db/ad.csv") +#ad <- read.taf("bootstrap/ad.csv") + +# prepare data ------------------------------- + +# keep only approved annotations +if (config$onlyApproved) { + ad <- ad[ad$IsApproved == 1, ] +} + +# add date columns + ad <- + within(ad, { + year <- year(parse_date_time(catch_date, "%d/%m/%Y %H:%M:%S")) + qtr <- quarter(parse_date_time(catch_date, "%d/%m/%Y %H:%M:%S")) + month <- month(parse_date_time(catch_date, "%d/%m/%Y %H:%M:%S")) + }) + + +# if variables are missing add "missing" +ad$ices_area[is.na(ad$ices_area) | ad$ices_area == ""] <- "missing" +ad$stock[is.na(ad$stock) | ad$stock == ""] <- "missing" +ad$prep_method[is.na(ad$prep_method) | ad$prep_method == ""] <- "missing" + +# if no advanced readers! make them all advanced +if (all(ad$expertise == 0)) { + msg("NOTE: all readers were Basic - all have been converted to Advanced") + ad$expertise[] <- 1 +} + +# convert reader expertise +ad$expertise <- c("Basic", "Advanced")[ad$expertise + 1] + +# Assign weight to the readers based in their ranking-experience +# Add different weight columns to data +weight <- length(sort(unique(ad$reader_number))):1 +reader_number <- sort(unique(ad$reader_number)) +reader <- data.frame(reader_number = reader_number, weight_I = weight, weight_II = 1 / (1 + log(sort(weight, decreasing = FALSE) + 0.0000000001))) +ad <- merge(ad, reader, by.x = "reader_number", by.y = "reader_number", all.x = TRUE) + +ad4webgr <- ad + +# Before calculating the mode, give to the sampleID of readings by eventOrganizer (histological sample) the same name as the samples analyzed by the other readers, so the maturity defined for the histological samples is assigned as mode to the other samples of the same FishID +fishid <- sort(unique(ad$FishID)) +for (i in 1:length(fishid)) +{ + nohist <- ad[ad$FishID == fishid[i] & ad$TypeAnnotation != "eventOrganizer", ] + sampleid_nohist <- unique(nohist$SampleID) + yeshist <- ad[ad$FishID == fishid[i] & ad$TypeAnnotation == "eventOrganizer", ] + if (dim(yeshist)[1] > 0) { + yeshist <- yeshist[rep(row.names(yeshist), length(sampleid_nohist)), ] + yeshist$SampleID <- sampleid_nohist + } + + temp <- rbind(nohist, yeshist) + + if (i == 1) { + result <- temp + } else { + result <- rbind(result, temp) + } +} + +ad <- result +ad$TypeAnnotation[ad$TypeAnnotation=="Reader"]<-"reader" +ad$TypeAnnotation[ad$TypeAnnotation=="Delegate"]<-"eventOrganizer" +ad$TypeAnnotation[ad$TypeAnnotation=="eventDelegate"]<-"eventOrganizer" +ad$TypeAnnotation[ad$TypeAnnotation=="Organizer"]<-"eventOrganizer" +ad$reader[ad$reader==""]<-"eventOrganizer" ## to add a name to the Reader column from the Event Organizer +ad$reader[ad$TypeAnnotation=="eventOrganizer"]<-"eventOrganizer" +ad$Sex[is.na(ad$Sex) | ad$Sex==""]<-"NI" +ad<-ad[ad$Sex!="NI",] +ad$Maturity[is.na(ad$Maturity) | ad$Maturity==""]<-"NI" +ad<-ad[ad$Maturity!="NI",] + + +# Calculate modal maturity stage and coefficient of unalikability of maturity stage +ad_long <- ad %>% + add_modal_trad(varmod = "Maturity", config$ma_method) %>% + add_modal_linearweight(varmod = "Maturity", config$ma_method) %>% + add_modal_negexpweight(varmod = "Maturity", config$ma_method) + +ad_long_adv <- ad[ad$expertise == "Advanced", ] %>% + add_modal_trad(varmod = "Maturity", config$ma_method) %>% + add_modal_linearweight(varmod = "Maturity", config$ma_method) %>% + add_modal_negexpweight(varmod = "Maturity", config$ma_method) + + +# Calculate modal maturity stage and coefficient of unalikability of sex category +ad_long <- ad_long %>% + add_modal_trad(varmod = "Sex", config$ma_method) %>% + add_modal_linearweight(varmod = "Sex", config$ma_method) %>% + add_modal_negexpweight(varmod = "Sex", config$ma_method) + +ad_long_adv <- ad_long_adv %>% + add_modal_trad(varmod = "Sex", config$ma_method) %>% + add_modal_linearweight(varmod = "Sex", config$ma_method) %>% + add_modal_negexpweight(varmod = "Sex", config$ma_method) + +# Choose the final mode (traditional, readers linear weight or negative exponential linear weight) based in the existence of histological samples or not, and, in case there are no histological samples, depending if there is multimodality or not. +ad_long <- select_mode(ad_long, config$ma_method, config$mode_definition) +ad_long_adv <- select_mode(ad_long_adv, config$ma_method, config$mode_definition) + +# prepare data in wbgr output format +# IMAGE,1,2,3,4,5,6,7,8,9,10,11,12,13 +# Expertise level,-,-,-,-,-,-,-,-,-,-,-,-,- +# Stock assessment,no,no,no,no,no,no,no,no,no,no,no,no,no +# 6698256.jpg,1,1,1,1,1,-,2,1,-,2,-,1,- +# 6698257.jpg,3,3,3,3,2,1,3,3,-,3,-,3,- +readers <- sort(unique(ad$reader)) + +webgr_maturity <- spread(ad4webgr[c("FishID", "reader", "Maturity")], key = reader, value = Maturity) +webgr_maturity[] <- paste(unlist(webgr_maturity)) +webgr_maturity[webgr_maturity == "NA"] <- "-" +webgr_maturity <- + rbind( + c("Expertise level", rep("-", length(readers))), + c("Stock assessment", rep("no", length(readers))), + webgr_maturity + ) +names(webgr_maturity) <- c("IMAGE", 1:length(readers)) +head(webgr_maturity) + +webgr_sex <- spread(ad4webgr[c("FishID", "reader", "Sex")], key = reader, value = Sex) +webgr_sex[] <- paste(unlist(webgr_sex)) +webgr_sex[webgr_sex == "NA"] <- "-" +webgr_sex <- + rbind( + c("Expertise level", rep("-", length(readers))), + c("Stock assessment", rep("no", length(readers))), + webgr_sex + ) +names(webgr_sex) <- c("IMAGE", 1:length(readers)) +head(webgr_sex) + + +# write out input data tables for use later +write.taf(ad, file = "data.csv", dir = "data", quote = TRUE) +write.taf(ad_long, dir = "data", quote = TRUE) +write.taf(ad_long_adv, dir = "data", quote = TRUE) +write.taf(webgr_maturity, file = "WebGR_maturity_all.csv", dir = "data", quote = TRUE) +write.taf(webgr_sex, file = "WebGR_sex_all.csv", dir = "data", quote = TRUE) diff --git a/model.R b/model.R index 9ee6de8..2b7007b 100644 --- a/model.R +++ b/model.R @@ -1,351 +1,434 @@ -## Run analysis, write model results - -## Before: data/ad_long.csv, data/ad_long_ex.csv, -## data/ad_wide.csv, data/ad_wide_ex.csv -## After: - -library(icesTAF) -library(jsonlite) -#unloadNamespace("dplyr") -library(plyr) # age error matrix -library(dplyr) -library(tidyr) -library(tibble) # bias_test - -library(ggplot2) -library(scales) # rescale_none - -# library ragree contains the function unalike, that is used to estimate the coefficient of unalikeability -taf.library(ragree) -# - -# make model directory -mkdir("model") - -# # load configuration -config <- read_json("bootstrap/data/config.json", simplifyVector = TRUE) - -# load utilities -source("utilities.R") -source("utilities_model.R") - -# read input data -ad_long_all <- read.taf("data/ad_long.csv") -ad_long_ex <- read.taf("data/ad_long_ex.csv") - - -# model maturity range -modal_matur_unique_all <- sort(unique(ad_long_all$modal_maturity)) -modal_matur_unique_ex <- sort(unique(ad_long_ex$modal_maturity)) - -# model sex range -modal_sex_unique_all <- sort(unique(ad_long_all$modal_sex)) -modal_sex_unique_ex <- sort(unique(ad_long_ex$modal_sex)) - -# set strata to NULL if all are NA -# if(length(setdiff("strata", names(config)))==0) {if(all(is.na(ad_long_all[["strata"]]))) config$strata <- NULL} - -# Overview of samples and stagers ############################################## - -# Sample overview -sample_data_overview <- sample_data_overview_table(ad_long_all, strata=config$strata) -write.taf(sample_data_overview, dir = "model") - -# Participants table -stager_data <- reader_data_table(ad_long_all, strata=config$strata) -write.taf(stager_data, dir = "model") - -# Results ############################################## - -# repeat for all and for experts only - -for (group in c("all", "ex")) { - #group <- "all" - #group <- "ex" - - # get the appropriate dataset - ad_long <- get(vname("ad_long")) - modal_matur_unique <- get(vname("modal_matur_unique")) - modal_sex_unique <- get(vname("modal_sex_unique")) - - ad_long$modal_maturity <- factor(ad_long$modal_maturity, levels = modal_matur_unique) - ad_long$modal_sex <- factor(ad_long$modal_sex, levels = modal_sex_unique) - ad_long$reader <- factor(ad_long$reader) - - ################################################################## - # number read table - assign( - vname("num_read_tab_modal_matur_by_reader"), - num_read_table_modal_maturity(ad_long, by = "reader") - ) - write.taf(vname("num_read_tab_modal_matur_by_reader"), dir = "model") - - assign( - vname("num_read_tab_modal_sex_by_reader"), - num_read_table_modal_sex(ad_long, by = "reader") - ) - write.taf(vname("num_read_tab_modal_sex_by_reader"), dir = "model") - - ################################################################################################################################################################ - # Calculate the number of cases with multiple modes when the traditional method is used, the linear weight or the negative exponential weighting. - - # First the Maturity stage - assign( - vname("multimode_cases_tab_traditional_Maturity"), - multimode_cases_table_traditional(ad_long, "Maturity") - ) - write.taf(vname("multimode_cases_tab_traditional_Maturity"), dir = "model") - - assign( - vname("multimode_cases_tab_linear_Maturity"), - multimode_cases_table_linear(ad_long, "Maturity") - ) - write.taf(vname("multimode_cases_tab_linear_Maturity"), dir = "model") - - assign( - vname("multimode_cases_tab_negexp_Maturity"), - multimode_cases_table_negexp(ad_long, "Maturity") - ) - write.taf(vname("multimode_cases_tab_negexp_Maturity"), dir = "model") - - assign( - vname("multimode_cases_tab_multistage_Maturity"), - multimode_cases_table_multistage(ad_long, "Maturity") - ) - write.taf(vname("multimode_cases_tab_multistage_Maturity"), dir = "model") - - # Next the sex category - assign( - vname("multimode_cases_tab_traditional_Sex"), - multimode_cases_table_traditional(ad_long, "Sex") - ) - write.taf(vname("multimode_cases_tab_traditional_Sex"), dir = "model") - - assign( - vname("multimode_cases_tab_linear_Sex"), - multimode_cases_table_linear(ad_long, "Sex") - ) - write.taf(vname("multimode_cases_tab_linear_Sex"), dir = "model") - - assign( - vname("multimode_cases_tab_negexp_Sex"), - multimode_cases_table_negexp(ad_long, "Sex") - ) - write.taf(vname("multimode_cases_tab_negexp_Sex"), dir = "model") - - assign( - vname("multimode_cases_tab_multistage_Sex"), - multimode_cases_table_multistage(ad_long, "Sex") - ) - write.taf(vname("multimode_cases_tab_multistage_Sex"), dir = "model") - - - ############################################################################################################################## - # CU table (coefficient of unalikeability) - assign( - vname("cu_tab_maturity"), - cu_table(ad_long, "Maturity", by = "reader") - ) - write.taf(vname("cu_tab_maturity"), dir = "model") - - assign( - vname("cu_tab_sex"), - cu_table(ad_long, "Sex", by = "reader") - ) - write.taf(vname("cu_tab_sex"), dir = "model") - - - ############################################################################################################################## - # Percent agreement between maturity stagings and modal maturity stage. - assign( - vname("pa_tab_maturity"), - pa_table(ad_long, "Maturity", by = "reader") - ) - write.taf(vname("pa_tab_maturity"), dir = "model") - - assign( - vname("pa_tab_sex"), - pa_table(ad_long, "Sex", by = "reader") - ) - write.taf(vname("pa_tab_sex"), dir = "model") - - ################################################################################################################################################## - # Frequency table (Number for each maturity stage per modal_maturity for each ). This is the equivalent to the relative bias table in the ageing. - assign( - vname("stager_bias_freq_tab_maturity"), - reader_freq_table(ad_long, "Maturity") - ) - write.taf(vname("stager_bias_freq_tab_maturity"), dir = "model") - - assign( - vname("stager_bias_freq_tab_sex"), - reader_freq_table(ad_long, "Sex") - ) - write.taf(vname("stager_bias_freq_tab_sex"), dir = "model") - - ######################################################################################################################################################### - # general Frequency table (Number for each maturity stage per modal_maturity for all s). This is the equivalent to the relative bias table in the ageing. - assign( - vname("general_bias_freq_tab_maturity"), - general_freq_table(ad_long, "Maturity") - ) - write.taf(vname("general_bias_freq_tab_maturity"), dir = "model") - - assign( - vname("general_bias_freq_tab_sex"), - general_freq_table(ad_long, "Sex") - ) - write.taf(vname("general_bias_freq_tab_sex"), dir = "model") - - - - ################################################################## - ## Annex tables################################################### - ################################################################## - - ################# - # data overview - assign( - vname("data_overview_tab_maturity"), - data_overview_table(ad_long, "Maturity", config$report_tokens) - ) - write.taf(vname("data_overview_tab_maturity"), dir = "model") - - - assign( - vname("data_overview_tab_sex"), - data_overview_table(ad_long, "Sex", config$report_tokens) - ) - write.taf(vname("data_overview_tab_sex"), dir = "model") - - ######################## - # maturity composition - assign( - vname("maturity_composition_tab"), - maturity_composition_table(ad_long, by = "reader") - ) - write.taf(vname("maturity_composition_tab"), dir = "model") - - #################### - # sex composition - assign( - vname("sex_composition_tab"), - sex_composition_table(ad_long, by = "reader") - ) - write.taf(vname("sex_composition_tab"), dir = "model") - - ############################################################# - # maturity staging error matrix (MSEM) only for advanced s - assign( - vname("msem"), - mat_stage_error_matrix(ad_long, by = config$strata) - ) - - saveRDS(get(vname("msem")), file = file.path("model", paste0(vname("msem"), ".rds"))) - - ######################################################### - # sex category error matrix (SSEM) only for advanced s - assign( - vname("ssem"), - sex_stage_error_matrix(ad_long, by = config$strata) - ) - - saveRDS(get(vname("ssem")), file = file.path("model", paste0(vname("ssem"), ".rds"))) - - - - ################################################################## - # Results by strata ############################################## - ################################################################## - - ################################################################################################################################################## - # Frequency table (Number for each maturity stage per modal_maturity for each ). This is the equivalent to the relative bias table in the ageing. - if (is.null(config$strata)) config$strata <- numeric() - for (stratum in config$strata) { - assign( - vsname("stager_bias_freq_tab_sex_by"), - reader_freq_table(ad_long, varmod="Sex", by = stratum) - ) - write.taf(vsname("stager_bias_freq_tab_sex_by"), dir = "model") - - - ################################################################## - # number read table - assign( - vsname("num_read_tab_modal_matur_by"), - num_read_table_modal_maturity(ad_long, by = stratum) - ) - write.taf(vsname("num_read_tab_modal_matur_by"), dir = "model") - - assign( - vsname("num_read_tab_modal_sex_by"), - num_read_table_modal_sex(ad_long, by = stratum) - ) - write.taf(vsname("num_read_tab_modal_sex_by"), dir = "model") - - - ############################################################################################################################## - # CU table (coefficient of unalikeability) - assign( - vsname("cu_tab_maturity_by"), - cu_table(ad_long, "Maturity", by = stratum) - ) - write.taf(vsname("cu_tab_maturity_by"), dir = "model") - - assign( - vsname("cu_tab_sex_by"), - cu_table(ad_long, "Sex", by = stratum) - ) - write.taf(vsname("cu_tab_sex_by"), dir = "model") - - - ############################################################################################################################## - # Percent agreement between maturity stagings and modal maturity stage. - assign( - vsname("pa_tab_maturity_by"), - pa_table(ad_long, "Maturity", by = stratum) - ) - write.taf(vsname("pa_tab_maturity_by"), dir = "model") - - assign( - vsname("pa_tab_sex_by"), - pa_table(ad_long, "Sex", by = stratum) - ) - write.taf(vsname("pa_tab_sex_by"), dir = "model") - - ################################## - ## Annex tables ################# - ################################## - - ######################## - # maturity composition - # # assign( - # - # # vname("maturity_composition_tab"), - # # maturity_composition_table(ad_long, by=stratum) - # # ) - # # write.taf(vname("maturity_composition_tab"), dir = "model") - - ############################################################ - # maturity staging error matrix (MSEM) only for advanced s - assign( - vname("msem"), - mat_stage_error_matrix(ad_long, by = stratum) - ) - saveRDS(get(vname("msem")), file = file.path("model", paste0(vname("msem"), ".rds"))) - - ######################################################### - # sex category error matrix (MSEM) only for advanced s - assign( - vname("ssem"), - sex_stage_error_matrix(ad_long, by = stratum) - ) - saveRDS(get(vname("ssem")), file = file.path("model", paste0(vname("ssem"), ".rds"))) - - } - -} - - - -#config$strata=select_strata +## Run analysis, write model results + +## Before: data/ad_long.csv, data/ad_long_ex.csv, +## data/ad_wide.csv, data/ad_wide_ex.csv +## After: + +library(icesTAF) +library(jsonlite) +#unloadNamespace("dplyr") +library(plyr) # age error matrix +library(dplyr) +library(tidyr) +library(tibble) # bias_test + +library(ggplot2) +library(scales) # rescale_none + +# library ragree contains the function unalike, that is used to estimate the coefficient of unalikeability +#taf.library(ragree) +# + +# make model directory +mkdir("model") + +# # load configuration +config <- read_json("bootstrap/data/config.json", simplifyVector = TRUE) +#config <- read_json("bootstrap/initial/data/config.json", simplifyVector = TRUE) + +# load utilities +source("utilities.R") +source("utilities_model.R") + +# read input data +ad_long_all <- read.taf("data/ad_long.csv") +ad_long_adv <- read.taf("data/ad_long_adv.csv") + + +# model maturity range +modal_matur_unique_all <- sort(unique(ad_long_all$modal_maturity)) +modal_matur_unique_adv <- sort(unique(ad_long_adv$modal_maturity)) + +# model sex range +modal_sex_unique_all <- sort(unique(ad_long_all$modal_sex)) +modal_sex_unique_adv <- sort(unique(ad_long_adv$modal_sex)) + +# set strata to NULL if all are NA +# if(length(setdiff("strata", names(config)))==0) {if(all(is.na(ad_long_all[["strata"]]))) config$strata <- NULL} + +# Overview of samples and stagers ############################################## + +# Sample overview +sample_data_overview <- sample_data_overview_table(ad_long_all, strata=config$strata) +write.taf(sample_data_overview, dir = "model") + +# Participants table +stager_data <- reader_data_table(ad_long_all, strata=config$strata) +stager_data <- slice(stager_data, 1:(n() - 1)) +write.taf(stager_data, dir = "model") + +# Results ############################################## + +# repeat for all and for experts only + +for (group in c("all", "adv")) { + #group <- "all" + #group <- "adv" + + # get the appropriate dataset + ad_long <- get(vname("ad_long")) + modal_matur_unique <- get(vname("modal_matur_unique")) + modal_sex_unique <- get(vname("modal_sex_unique")) + + ad_long$modal_maturity <- factor(ad_long$modal_maturity, levels = modal_matur_unique) + ad_long$modal_sex <- factor(ad_long$modal_sex, levels = modal_sex_unique) + ad_long$reader <- factor(ad_long$reader) + + ################################################################## + # number read table + assign( + vname("num_read_tab_modal_matur_by_reader"), + num_read_table_modal_maturity(ad_long, by = "reader") + ) + write.taf(vname("num_read_tab_modal_matur_by_reader"), dir = "model") + + assign( + vname("num_read_tab_modal_sex_by_reader"), + num_read_table_modal_sex(ad_long, by = "reader") + ) + write.taf(vname("num_read_tab_modal_sex_by_reader"), dir = "model") + + ################################################################################################################################################################ + # Calculate the number of cases with multiple modes when the traditional method is used, the linear weight or the negative exponential weighting. + + # First the Maturity stage + assign( + vname("multimode_cases_tab_traditional_Maturity"), + multimode_cases_table_traditional(ad_long, "Maturity") + ) + write.taf(vname("multimode_cases_tab_traditional_Maturity"), dir = "model") + + assign( + vname("multimode_cases_tab_linear_Maturity"), + multimode_cases_table_linear(ad_long, "Maturity") + ) + write.taf(vname("multimode_cases_tab_linear_Maturity"), dir = "model") + + assign( + vname("multimode_cases_tab_negexp_Maturity"), + multimode_cases_table_negexp(ad_long, "Maturity") + ) + write.taf(vname("multimode_cases_tab_negexp_Maturity"), dir = "model") + + assign( + vname("multimode_cases_tab_multistage_Maturity"), + multimode_cases_table_multistage(ad_long, "Maturity") + ) + write.taf(vname("multimode_cases_tab_multistage_Maturity"), dir = "model") + + # Next the sex category + assign( + vname("multimode_cases_tab_traditional_Sex"), + multimode_cases_table_traditional(ad_long, "Sex") + ) + write.taf(vname("multimode_cases_tab_traditional_Sex"), dir = "model") + + assign( + vname("multimode_cases_tab_linear_Sex"), + multimode_cases_table_linear(ad_long, "Sex") + ) + write.taf(vname("multimode_cases_tab_linear_Sex"), dir = "model") + + assign( + vname("multimode_cases_tab_negexp_Sex"), + multimode_cases_table_negexp(ad_long, "Sex") + ) + write.taf(vname("multimode_cases_tab_negexp_Sex"), dir = "model") + + assign( + vname("multimode_cases_tab_multistage_Sex"), + multimode_cases_table_multistage(ad_long, "Sex") + ) + write.taf(vname("multimode_cases_tab_multistage_Sex"), dir = "model") + + + ############################################################################################################################## + # CU table (coefficient of unalikeability) + assign( + vname("cu_tab_maturity"), + cu_table(ad_long, "Maturity", by = "reader") + ) + write.taf(vname("cu_tab_maturity"), dir = "model") + + # CU table (coefficient of unalikeability) - Females + assign( + vname("cu_tab_maturity_females"), + cu_table(ad_long[ad_long$Sex=="F",], "Maturity", by = "reader") + ) + write.taf(vname("cu_tab_maturity_females"), dir = "model") + + + # CU table (coefficient of unalikeability) - Males + assign( + vname("cu_tab_maturity_males"), + cu_table(ad_long[ad_long$Sex=="M",], "Maturity", by = "reader") + ) + write.taf(vname("cu_tab_maturity_males"), dir = "model") + + + # CU table (coefficient of unalikeability) - by Sex + assign( + vname("cu_tab_sex"), + cu_table(ad_long, "Sex", by = "reader") + ) + write.taf(vname("cu_tab_sex"), dir = "model") + + + ############################################################################################################################## + # Percent agreement between maturity stagings and modal maturity stage. + assign( + vname("pa_tab_maturity"), + pa_table(ad_long, "Maturity", by = "reader") + ) + write.taf(vname("pa_tab_maturity"), dir = "model") + + assign( + vname("pa_tab_sex"), + pa_table(ad_long, "Sex", by = "reader") + ) + write.taf(vname("pa_tab_sex"), dir = "model") + + ##females + assign( + vname("pa_tab_maturity_females"), + pa_table(ad_long[ad_long$Sex=="F",], "Maturity", by = "reader") + ) + write.taf(vname("pa_tab_maturity"), dir = "model") + + + ##males + assign( + vname("pa_tab_maturity_males"), + pa_table(ad_long[ad_long$Sex=="M",], "Maturity", by = "reader") + ) + write.taf(vname("pa_tab_maturity_males"), dir = "model") + + + ################################################################################################################################################## + # Frequency table (Number for each maturity stage per modal_maturity for each ). This is the equivalent to the relative bias table in the ageing. + assign( + vname("stager_bias_freq_tab_maturity"), + reader_freq_table(ad_long, "Maturity") + ) + write.taf(vname("stager_bias_freq_tab_maturity"), dir = "model") + + assign( + vname("stager_bias_freq_tab_sex"), + reader_freq_table(ad_long, "Sex") + ) + write.taf(vname("stager_bias_freq_tab_sex"), dir = "model") + + ######################################################################################################################################################### + # general Frequency table (Number for each maturity stage per modal_maturity for all s). This is the equivalent to the relative bias table in the ageing. + assign( + vname("general_bias_freq_tab_maturity"), + general_freq_table(ad_long, "Maturity") + ) + write.taf(vname("general_bias_freq_tab_maturity"), dir = "model") + + ##females + assign( + vname("general_bias_freq_tab_maturity_females"), + general_freq_table(ad_long[ad_long$Sex=="F",], "Maturity") + ) + write.taf(vname("general_bias_freq_tab_maturity_females"), dir = "model") + + + ##males + assign( + vname("general_bias_freq_tab_maturity_males"), + general_freq_table(ad_long[ad_long$Sex=="M",], "Maturity") + ) + write.taf(vname("general_bias_freq_tab_maturity_males"), dir = "model") + + + ### Sex + assign( + vname("general_bias_freq_tab_sex"), + general_freq_table(ad_long, "Sex") + ) + write.taf(vname("general_bias_freq_tab_sex"), dir = "model") + + + + ################################################################## + ## Annex tables################################################### + ################################################################## + + ################# + # data overview + assign( + vname("data_overview_tab_maturity"), + data_overview_table(ad_long, "Maturity", config$report_tokens) + ) + write.taf(vname("data_overview_tab_maturity"), dir = "model") + + + assign( + vname("data_overview_tab_sex"), + data_overview_table(ad_long, "Sex", config$report_tokens) + ) + write.taf(vname("data_overview_tab_sex"), dir = "model") + + ######################## + # maturity composition + assign( + vname("maturity_composition_tab"), + maturity_composition_table(ad_long, by = "reader") + ) + write.taf(vname("maturity_composition_tab"), dir = "model") + + + + + #################### + # sex composition + assign( + vname("sex_composition_tab"), + sex_composition_table(ad_long, by = "reader") + ) + write.taf(vname("sex_composition_tab"), dir = "model") + + ############################################################# + # maturity staging error matrix (MSEM) only for advanced s + assign( + vname("msem"), + mat_stage_error_matrix(ad_long, by = config$strata) + ) + + saveRDS(get(vname("msem")), file = file.path("model", paste0(vname("msem"), ".rds"))) + + + ######################################################### + # sex category error matrix (SSEM) only for advanced s + assign( + vname("ssem"), + sex_stage_error_matrix(ad_long, by = config$strata) + ) + + saveRDS(get(vname("ssem")), file = file.path("model", paste0(vname("ssem"), ".rds"))) + + + + ################################################################## + # Results by strata ############################################## + ################################################################## + + ################################################################################################################################################## + # Frequency table (Number for each maturity stage per modal_maturity for each ). This is the equivalent to the relative bias table in the ageing. + if (is.null(config$strata)) config$strata <- numeric() + for (stratum in config$strata) { + assign( + vsname("stager_bias_freq_tab_sex_by"), + reader_freq_table(ad_long, varmod="Sex", by = stratum) + ) + write.taf(vsname("stager_bias_freq_tab_sex_by"), dir = "model") + + + ################################################################## + # number read table + assign( + vsname("num_read_tab_modal_matur_by"), + num_read_table_modal_maturity(ad_long, by = stratum) + ) + write.taf(vsname("num_read_tab_modal_matur_by"), dir = "model") + + assign( + vsname("num_read_tab_modal_sex_by"), + num_read_table_modal_sex(ad_long, by = stratum) + ) + write.taf(vsname("num_read_tab_modal_sex_by"), dir = "model") + + + ############################################################################################################################## + # CU table (coefficient of unalikeability) + assign( + vsname("cu_tab_maturity_by"), + cu_table(ad_long, "Maturity", by = stratum) + ) + write.taf(vsname("cu_tab_maturity_by"), dir = "model") + + assign( + vsname("cu_tab_sex_by"), + cu_table(ad_long, "Sex", by = stratum) + ) + write.taf(vsname("cu_tab_sex_by"), dir = "model") + + + # CU table (coefficient of unalikeability) - Males + assign( + vname("cu_tab_maturity_males"), + cu_table(ad_long[ad_long$Sex=="M",], "Maturity", by = "reader") + ) + write.taf(vname("cu_tab_maturity_males"), dir = "model") + + + # CU table (coefficient of unalikeability) - by Sex + assign( + vname("cu_tab_sex"), + cu_table(ad_long, "Sex", by = "reader") + ) + write.taf(vname("cu_tab_sex"), dir = "model") + + ############################################################################################################################## + # Percent agreement between maturity stagings and modal maturity stage. + assign( + vsname("pa_tab_maturity_by"), + pa_table(ad_long, "Maturity", by = stratum) + ) + write.taf(vsname("pa_tab_maturity_by"), dir = "model") + + assign( + vsname("pa_tab_sex_by"), + pa_table(ad_long, "Sex", by = stratum) + ) + write.taf(vsname("pa_tab_sex_by"), dir = "model") + + ##females + assign( + vname("pa_tab_maturity_females"), + pa_table(ad_long[ad_long$Sex=="F",], "Maturity", by = "reader") + ) + write.taf(vname("pa_tab_maturity"), dir = "model") + + + ##males + assign( + vname("pa_tab_maturity_males"), + pa_table(ad_long[ad_long$Sex=="M",], "Maturity", by = "reader") + ) + write.taf(vname("pa_tab_maturity_males"), dir = "model") + + ################################## + ## Annex tables ################# + ################################## + + ######################## + # maturity composition + # # assign( + # + # # vname("maturity_composition_tab"), + # # maturity_composition_table(ad_long, by=stratum) + # # ) + # # write.taf(vname("maturity_composition_tab"), dir = "model") + + ############################################################ + # maturity staging error matrix (MSEM) only for advanced s + assign( + vname("msem"), + mat_stage_error_matrix(ad_long, by = stratum) + ) + saveRDS(get(vname("msem")), file = file.path("model", paste0(vname("msem"), ".rds"))) + + + ######################################################### + # sex category error matrix (MSEM) only for advanced s + assign( + vname("ssem"), + sex_stage_error_matrix(ad_long, by = stratum) + ) + saveRDS(get(vname("ssem")), file = file.path("model", paste0(vname("ssem"), ".rds"))) + + } + +} diff --git a/report.R b/report.R index dddbd4e..d147cbb 100644 --- a/report.R +++ b/report.R @@ -8,7 +8,7 @@ library(icesTAF) library(rmarkdown) library(jsonlite) library(knitr) -# +# library(pander) library(ggplot2) library(scales) @@ -24,11 +24,12 @@ source("utilities.R") source("utilities_report.R") # load configuration data -config <- read_json("bootstrap/data/config.json", simplifyVector = TRUE) +#config <- read_json("bootstrap/data/config.json", simplifyVector = TRUE) +config <- read_json("bootstrap/initial/data/config.json", simplifyVector = TRUE) # load data for report ad_long_all <- read.taf("data/ad_long.csv") -ad_long_ex <- read.taf("data/ad_long_ex.csv") +ad_long_adv <- read.taf("data/ad_long_adv.csv") # # set strata to NULL is all are NA # if (all(is.na(ad_long_all[[config$strata]]))) config$strata <- NULL diff --git a/report_full.Rmd b/report_full.Rmd index 50943aa..b347663 100644 --- a/report_full.Rmd +++ b/report_full.Rmd @@ -1,2123 +1,2413 @@ ---- -output: - word_document: - fig_caption: true - fig_height: 10 - fig_width: 10 - reference_docx: bootstrap/data/reportTemplate.docx - toc: true - keep_md: false -params: - report_title: "" - strata: NULL ---- - -```{r introduction, include = FALSE} - -# INTRODUCTION ################################################################ - -# This markdown documents and integrated scripts analyse biological -# maturity readings. -# The output is a .docx template that includes -# the results of the analysis and should be used as a standard for -# reporting of maturity staging comparisons. - -``` - -```{r chunk_setup, include=FALSE} -# CHUNK SETUPS ################################################################# - -knitr::opts_chunk$set(echo = FALSE, warning = FALSE, - message=FALSE, results = 'asis', dpi=400) -``` - -```{r pander_settings, include = FALSE} - -# PANDER OPTIONS ############################################################## - -panderOptions('table.split.table', Inf) -panderOptions('keep.trailing.zeros', TRUE) -panderOptions('table.alignment.default', "center") - -``` - ---- -title: `r config$report_title` ---- - -# Executive summary - -# Terms of reference - -# Agenda and participant list - -The agenda can be found in Annex 1 and the list of participants in Annex 2. - - -# Introduction - -This part should include a background to the species, the workshop/exchange -and what to expect to read about in the report. - - -# Methods - -This report contains statistical analyses and comparisons of sex categorization and maturity stagings -in the form of tables and graphical plots. - -First, an overview of participating maturity stagers and the samples are presented. - -Before each table or plot there is a short explanation of it. This text is -thought as a help to understand the tables/plot and can just be deleted in -the final output report. The document can be edited just like any other -.docx file. New text can be added, additional pictures can be included and -the tables edited. If some tables which are presently in the annexes need to -be moved to the body of the report this is also possible. Only the plots -cannot be changed. - -In the first part of analysis some of the tables and plots from the -Guus Eltink Excel sheet 'Age Reading Comparisons' **(Eltink, A.T.G.W. 2000)** are presented, that can be used for "maturity staging comparisons". Since the sex categories and maturity stages are categorical in comparison to ages that are quantitative, some of the statistics presented in the report of the age exchange events cannot be calculated here, like the average percentage error or the coefficient of variation. The later is being replaced by the coefficient of unlikeability (see bellow). - - -**Pecentage Agreement** - -The percentage agreement tells how large is the part of sex categorizations/maturity stagings that are equal to the modal sex/maturity. The percentage agreement is estimated by modal sex or maturity and stager as the proportion (as percentage) of times that the lectures of that stager agreed with the resulting modal sex or maturity.This percentage is estimated as the number of times that a stager agreed with the modal sex or maturity divided by the total number of gonads analyzed by a stager for that modal maturity stage. - -$$PA = { \frac{ {number \,of \,readings \,that \,agree \,with \,modal \,maturity}} {total \,number \,of \,readings \,by \,modal \,maturity} } \cdot {100 \%} $$ - - -**Coefficient of unalikeability (CU)** - -The concept of unalikeability (Kader and Perry, 2007) focuses on how often observations differ within a group. Specifically, for the sex/maturity staging events, the CU provides a measure of how alike, for each modal maturity stage, the stages decided by each stager were (or all stagers at once). The cu ranges between 0 and 1. The higher the cu value, the more unalike the data are. - -The table presents the CU per modal maturity stage and stager. For the case of a finite number of observations (n), a finite number of categories (m) and a finite number of objects, k_i, within category i, will allow expression of the coefficient of unalikeablity as: - -$$ u = {1 - ∑ p_i ^ 2} $$ - -where - -$$ p_i = {k_i / n} $$ - -To the table is also added the CU of all stagers combined per modal maturity and -a weighted mean of the CU per stager. - - -**Maturity stage error matrix (MSEM)** - -Maturity stage error matrices (AEM) were produced following the same procedures outlined - by WKSABCAL (2014) for calculating the "Age error matrix". MSEM shows the proportion of each modal -maturity stage mis-staged as other stages. The sum of each row is 1, which equals 100%. -The maturity data was analysed twice, the first time all stagers were included -and the second time only the “advanced” stagers were included. If a stager -is “advanced” then they are considered well trained and they provide maturity stages -for stock assessment or similar purposes. When the MSEM is compiled for -assessment purposes it uses only those stagers who provide maturity data for -the stock assessment in that specific area. - - - -# Analysis of maturity staging calibration exercise (ToR?) - - -## Overview of samples and stagers - - -\br - -```{r sample_overview} - -# Table caption -set.caption( - paste( - '**Table X:** Overview of samples used for the exchange event number ', - config$event_id, - '. The number of samples (num) are shown by year, ices area, season (qtr). The range of maturity stages (mat_stages) and sex (sex_cat) assigned by the group or stagers participating in the event is also shown, as well as the range of modal maturity stages (modal_mat_stages) and modal sex (modal_sex_cat) decided for all samples.', - sep="" - ) -) - -# Output table -pander(style_table0(sample_data_overview), style = "simple") - -``` - -\br - -```{r participants_overview} - -# PARTICIPATANTS OVERVIEW ##################################################### - -# Table caption -set.caption('**Table X:** Overview of stagers participating in the event, with their overal expertise (Expertise), and their ranking position based in their experience (only valid if the experience weighting protocol has been applied).') - -# Output table -pander(style_table0(filter(stager_data))) - -``` - - -## Results - -Text? - -### All stagers - -Text? - -#### All samples together - -Text? - -**Multimodal cases** - -Those writing the report put TEXT here describing the results. - -If there were available histological samples for all the fish individuals in the exchange event, no multiple modes are expected, however, if there were no histology samples for some of them, multiple maturity stage modes could be found. In the next table this information is presented for the sex category determination. - - - -```{r summary_multiple_modes_sex_all} - -histN=ad_long_all %>% subset(TypeAnnotation=="eventOrganizer") %>% select(FishID) %>% unique() %>% dim() -allN=ad_long_all %>% select(FishID) %>% unique() %>% dim() -perc_not_Hist=100-round(100*histN[1]/allN[1], 0) - -summary_multiple_modes_sex_all <- - c(NSample = length(unique(ad_long_all$FishID)), - Nhist=histN[1], - Perc_not_Hist=perc_not_Hist, - PercMM_traditional=unique(ifelse(multimode_cases_tab_traditional_Sex_all$NModes_trad=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_traditional_Sex_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), - PercMM_linear_weight=unique(ifelse(multimode_cases_tab_linear_Sex_all$NModes_linear=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_linear_Sex_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), - PercMM_negexp_weight=unique(ifelse(multimode_cases_tab_negexp_Sex_all$NModes_negexp=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_negexp_Sex_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), - PercMM_multistage=unique(ifelse(multimode_cases_tab_multistage_Sex_all$NModes_multistage=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_multistage_Sex_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%")))) - -# Table caption -set.caption('**Table X:** Summary of statistics for sex staging; Total number of fish individuals studied (NSample), number of fish individuals with histological samples (Nhist), percentage of fish individuals without histology (Perc_not_Hist). The percentage of cases (fish samples) with multiple modes depending on the approach to weight the experience of the stager which will be considered when defining the fish sex stage mode. PercMM_traditional shows the percentage of the total samples for which multiple modes are obtained when all the stagers are equally weighted. PercMM_linear_weight shows the percentage of the total samples for which multiple modes are obtained when the weight assigned to the different stagers decreases linearly with the experience, while in the PercMM_negexp the weight applied decreases with a negative exponential shape with the experience. The PercMM_multistage shows the percentage of multiple mode cases when a combination of the different methodologies is used, as explained in the material and methods section') - -# Output table -pander(style_table0(summary_multiple_modes_sex_all), missing="") - -``` - - -If there were available histological samples for all the fish individuals in the exchange event, no multiple modes are expected, however, if there were no histology samples for some of them, multiple maturity stage modes could be found. In the next table this information is presented for the maturity staging. - - - -```{r summary_multiple_modes_maturity_all} - -histN=ad_long_all %>% subset(TypeAnnotation=="eventOrganizer") %>% select(FishID) %>% unique() %>% dim() -allN=ad_long_all %>% select(FishID) %>% unique() %>% dim() -perc_not_Hist=100-round(100*histN[1]/allN[1], 0) - -summary_multiple_modes_maturity_all <- - c(NSample = length(unique(ad_long_all$FishID)), - Nhist=histN[1], - Perc_not_Hist=perc_not_Hist, - PercMM_traditional=unique(ifelse(multimode_cases_tab_traditional_Maturity_all$NModes_trad=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_traditional_Maturity_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), - PercMM_linear_weight=unique(ifelse(multimode_cases_tab_linear_Maturity_all$NModes_linear=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_linear_Maturity_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), - PercMM_negexp_weight=unique(ifelse(multimode_cases_tab_negexp_Maturity_all$NModes_negexp=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_negexp_Maturity_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), - PercMM_multistage=unique(ifelse(multimode_cases_tab_multistage_Maturity_all$NModes_multistage=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_multistage_Maturity_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%")))) - -# Table caption -set.caption('**Table X:** Summary of statistics for maturity staging; Total number of fish individuals studied (NSample), number of fish individuals with histological samples (Nhist), percentage of fish individuals without histology (Perc_not_Hist). The percentage of cases (fish samples) with multiple modes depending on the approach to weight the experience of the stager which will be considered when defining the fish maturity stage mode. PercMM_traditional shows the percentage of the total samples for which multiple modes are obtained when all the stagers are equally weighted. PercMM_linear_weight shows the percentage of the total samples for which multiple modes are obtained when the weight assigned to the different stagers decreases linearly with the experience, while in the PercMM_negexp the weight applied decreases with a negative exponential shape with the experience. The PercMM_multistage shows the percentage of multiple mode cases when a combination of the different methodologies is used, as explained in the material and methods section') - -# Output table -pander(style_table0(summary_multiple_modes_maturity_all), missing="") - - - -``` - - - -```{r summary_statistics_sex_all} -summary_statistics_sex_all <- - c(NSample = length(unique(ad_long_all$FishID)), - CU = cu_tab_sex_all[nrow(cu_tab_sex_all), "Total"], - PA = pa_tab_sex_all[nrow(pa_tab_sex_all), "Total"]) - -summary_statistics_maturity_all <- - c(NSample = length(unique(ad_long_all$FishID)), - CU = cu_tab_maturity_all[nrow(cu_tab_maturity_all), "Total"], - PA = pa_tab_maturity_all[nrow(pa_tab_maturity_all), "Total"]) -``` - -The average percentage agreement by modal sex category for all stagers was `r summary_statistics_sex_all[3]`, with a weighted average CU of `r summary_statistics_sex_all[2]`. Regarding the maturation staging, the percentage agreement by modal maturation stage was `r summary_statistics_maturity_all[3]`, and the weight average CU was `r summary_statistics_maturity_all[2]`. - - -**Sex categorization table** - -Text? - -```{r num_read_sex_all} - -# Sex categorization TABLE - All stagers ################################################# -data=num_read_tab_modal_sex_by_reader_all -#data=sex_composition_tab_all -nstagers=length(unique(stager_data$'Reader code')) -ncols=dim(data[,!colnames(data) %in% c("Total")])[2] -basicols=ncols-nstagers - -z=8 -nfig=round((dim(data)[2]-basicols)/z) -N=min(z,ncols) -for(i in 1:nfig) -{ - if(i==1){ - # Table caption - set.caption('**Table X:** Sex categorization table: presents the - number of categorizations made per expert for each modal sex category.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] - pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt - } else { if(i -The percentage agreement per stager and modal maturity stage represent the proportion of the total number of stagings that are equal to the modal maturity stage. The weighted mean including at the bottom of the table is weighted according to number of maturity stagings. - -```{r percentage_agreement_sex_all} - -# PERCENTAGE AGREEMENT TABLE by modal sex category - All stagers #################################### - -data=pa_tab_sex_all -nstagers=length(unique(stager_data$'Reader code')) -ncols=dim(data[,!colnames(data) %in% c("Total")])[2] -basicols=ncols-nstagers - -z=8 -nfig=round((dim(data)[2]-basicols)/z) -N=min(z,ncols) -for(i in 1:nfig) -{ - if(i==1){ - # Table caption - set.caption('**Table X:** Percentage agreement (PA) table: shows the PA per modal sex category and stager, the PA for all stagers combined and the weighted mean of the PA per stager.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] - pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt - } else { if(i -The frequency bias is shown by sex category for each modal sex category per stager. As for the previous tables, a combined frequency distribution for all stagers are calculated. This frequency distribution presented in the form of a table or a buble plot may be used to explore the bias in the sex determination by the stagers participating in this event. - -```{r freq_sex_all} - -# frequency bias TABLE by sex category- per stagers ########################################### - -data=stager_bias_freq_tab_sex_all -nstagers=length(unique(stager_data$'Reader code')) -ncols=dim(data[,!colnames(data) %in% c("Total")])[2] -basicols=ncols-nstagers - -z=8 -nfig=round((dim(data)[2]-basicols)/z) -N=min(z,ncols) -for(i in 1:nfig) -{ - if(i==1){ - # Table caption - set.caption('**Table X:** Frequency bias table: represents the frequency per - modal sex category per stager and the frequency distribution of all stagers combined - per modal sex category.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] - pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt - } else { if(i -The frequency bias is shown by maturity stages for each modal maturity stage per stager. As for the previous tables, a combined frequency distribution for all stagers are calculated. This frequency distribution presented in the form of a table or a buble plot may be used to explore the bias in the maturity determination by the stagers participating in this event. - -```{r freq_maturity_all} - -# frequency bias TABLE by maturity stage - per stagers ########################################### - -data=stager_bias_freq_tab_maturity_all -nstagers=length(unique(stager_data$'Reader code')) -ncols=dim(data[,!colnames(data) %in% c("Total")])[2] -basicols=ncols-nstagers - -z=8 -nfig=round((dim(data)[2]-basicols)/z) -N=min(z,ncols) -for(i in 1:nfig) -{ - if(i==1){ - # Table caption - set.caption('**Table X:** Frequency bias table: represents the frequency per - modal maturity stage per stager and the frequency distribution of all stagers combined - per modal maturity stage.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] - pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt - } else { if(i= istrata -``` - - - -```{r strata_plus_one} -# second strata -istrata <- istrata + 1 -``` - -```{r strata_title, eval = print_strata()} -stratum <- params$strata[istrata] -# Section title -asis_output(paste0("**Results by ", stratum, "**")) -``` - - -```{r ar_sex_title, eval = print_strata()} -# title -asis_output(paste0("**Sex categorization by ", stratum, "**")) -``` - -```{r ar_sex_by, eval = print_strata()} -# NUMBER OF SEX CATEGORIZATION PER Strata - Advanced stagers ######################### - -# Table caption -set.caption(paste0('**Table X:** Number of sex categorization per ', stratum, ' for all the stagers.')) - -# Output table -pander(style_table1(get(vsname("num_read_tab_modal_sex_by"))), missing = "-") -``` - - -```{r ar_matur_title, eval = print_strata()} -# title -asis_output(paste0("**Maturity staging by ", stratum, "**")) -``` - -```{r ar_matur_by, eval = print_strata()} -# NUMBER OF MATURITY STAGING PER strata - Advanced stagers ######################### - -# Table caption -set.caption(paste0('**Table X:** Number of maturity stagings per ', stratum, ' for all the stagers.')) - -# Output table -pander(style_table1(get(vsname("num_read_tab_modal_matur_by"))), missing = "-") -``` - - - -```{r cu_sex_title, eval = print_strata()} -# title -asis_output(paste0("**Coefficient of Unalikeability by modal sex category per", stratum, "**")) -``` - - -```{r cu_sex_by, eval = print_strata()} - -# Table caption -set.caption(paste0('**Table X:** CU by modal sex category per ', stratum, '.')) - -# Output table -pander(style_table2(get(vsname("cu_tab_sex_by"))), missing = "-") - -``` - - -```{r cu_matur_title, eval = print_strata()} -# title -asis_output(paste0("**Coefficient of Unalikeability by modal maturity stage per", stratum, "**")) -``` - - -```{r cu_matur_by, eval = print_strata()} - -# Table caption -set.caption(paste0('**Table X:** CU by modal maturity stage per ', stratum, '.')) - -# Output table -pander(style_table2(get(vsname("cu_tab_maturity_by"))), missing = "-") - -``` - - -```{r pa_sex_title, eval = print_strata()} -# title -asis_output(paste0("**Percentage Agreement by sex category per ", stratum, "**")) -``` - -```{r pa_sex_by, eval = print_strata()} -# Section title -#asis_("## PA by sex category - Advanced stagers") # # ############################################# - -# Table caption -set.caption(paste0('**Table X:** Percentage Agreement by sex category per ', stratum, '.')) - -# Output table -pander(style_table2(get(vsname("pa_tab_sex_by"))), missing = "-") - -``` - - -```{r pa_matur_title, eval = print_strata()} -# title -asis_output(paste0("**Percentage Agreement by maturation stage per ", stratum, "**")) -``` - -```{r pa_matur_by, eval = print_strata()} -# Section title -#asis_("## PA by maturation stage - Advanced stagers") # # ############################################# - -# Table caption -set.caption(paste0('**Table X:** Percentage Agreement by maturation stage per ', stratum, '.')) - -# Output table -pander(style_table2(get(vsname("pa_tab_maturity_by"))), missing = "-") - -``` - - - - -```{r ref.label='strata_plus_one'} -``` - -```{r ref.label='strata_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_by', eval = print_strata()} -``` - -```{r ref.label='ar_matur_title', eval = print_strata()} -``` - -```{r ref.label='ar_matur_by', eval = print_strata()} -``` - -```{r ref.label='cu_sex_title', eval = print_strata()} -``` - -```{r ref.label='cu_sex_by', eval = print_strata()} -``` - -```{r ref.label='cu_matur_title', eval = print_strata()} -``` - -```{r ref.label='cu_matur_by', eval = print_strata()} -``` - -```{r ref.label='pa_sex_title', eval = print_strata()} -``` - -```{r ref.label='pa_sex_by', eval = print_strata()} -``` - -```{r ref.label='pa_matur_title', eval = print_strata()} -``` - -```{r ref.label='pa_matur_by', eval = print_strata()} -``` - - - - -```{r ref.label='strata_plus_one'} -``` - -```{r ref.label='strata_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_by', eval = print_strata()} -``` - -```{r ref.label='ar_matur_title', eval = print_strata()} -``` - -```{r ref.label='ar_matur_by', eval = print_strata()} -``` - -```{r ref.label='cu_sex_title', eval = print_strata()} -``` - -```{r ref.label='cu_sex_by', eval = print_strata()} -``` - -```{r ref.label='cu_matur_title', eval = print_strata()} -``` - -```{r ref.label='cu_matur_by', eval = print_strata()} -``` - -```{r ref.label='pa_sex_title', eval = print_strata()} -``` - -```{r ref.label='pa_sex_by', eval = print_strata()} -``` - -```{r ref.label='pa_matur_title', eval = print_strata()} -``` - -```{r ref.label='pa_matur_by', eval = print_strata()} -``` - - - - -```{r ref.label='strata_plus_one'} -``` - -```{r ref.label='strata_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_by', eval = print_strata()} -``` - -```{r ref.label='ar_matur_title', eval = print_strata()} -``` - -```{r ref.label='ar_matur_by', eval = print_strata()} -``` - -```{r ref.label='cu_sex_title', eval = print_strata()} -``` - -```{r ref.label='cu_sex_by', eval = print_strata()} -``` - -```{r ref.label='cu_matur_title', eval = print_strata()} -``` - -```{r ref.label='cu_matur_by', eval = print_strata()} -``` - -```{r ref.label='pa_sex_title', eval = print_strata()} -``` - -```{r ref.label='pa_sex_by', eval = print_strata()} -``` - -```{r ref.label='pa_matur_title', eval = print_strata()} -``` - -```{r ref.label='pa_matur_by', eval = print_strata()} -``` - - - - - -## Advanced stagers - -Those writing the report put TEXT here describing the results. - -#### All samples together - -***Multimodal cases*** - -Those writing the report put TEXT here describing the results. - -If there were available histological samples for all the fish individuals in the exchange event, no multiple modes are expected, however, if there were no hitology samples for some of them, multiple maturity stage modes could be found. In the next table this information is presented for the sex category determination. - - - -```{r summary_multiple_modes_sex_ex} - -histN=ad_long_ex %>% subset(TypeAnnotation=="eventOrganizer") %>% select(FishID) %>% unique() %>% dim() -allN=ad_long_ex %>% select(FishID) %>% unique() %>% dim() -perc_not_Hist=100-round(100*histN[1]/allN[1], 0) - -summary_multiple_modes_sex_ex <- - c(NSample = length(unique(ad_long_ex$FishID)), - Nhist=histN[1], - Perc_not_Hist=perc_not_Hist, - PercMM_traditional=unique(ifelse(multimode_cases_tab_traditional_Sex_ex$NModes_trad=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_traditional_Sex_ex)/length(unique(ad_long_ex$SampleID)))*100, digits=0),"%"))), - PercMM_linear_weight=unique(ifelse(multimode_cases_tab_linear_Sex_ex$NModes_linear=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_linear_Sex_ex)/length(unique(ad_long_ex$SampleID)))*100, digits=0),"%"))), - PercMM_negexp_weight=unique(ifelse(multimode_cases_tab_negexp_Sex_ex$NModes_negexp=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_negexp_Sex_ex)/length(unique(ad_long_ex$SampleID)))*100, digits=0),"%"))), - PercMM_multistage=unique(ifelse(multimode_cases_tab_multistage_Sex_ex$NModes_multistage=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_multistage_Sex_ex)/length(unique(ad_long_ex$SampleID)))*100, digits=0),"%")))) - -# Table caption -set.caption('**Table X:** Summary of statistics for sex staging; Total number of fish individuals studied (NSample), number of fish individuals with histological samples (Nhist), percentage of fish individuals without histology (Perc_not_Hist). The percentage of cases (fish samples) with multiple modes depending on the approach to weight the experience of the stager which will be considered when defining the fish sex stage mode. PercMM_traditional shows the percentage of the total samples for which multiple modes are obtained when all the stagers are equally weighted. PercMM_linear_weight shows the percentage of the total samples for which multiple modes are obtained when the weight assigned to the different stagers decreases linearly with the experience, while in the PercMM_negexp the weight applied decreases with a negative exponential shape with the experience. The PercMM_multistage shows the percentage of multiple mode cases when a combination of the different methodologies is used, as explained in the material and methods section') - -# Output table -pander(style_table0(summary_multiple_modes_sex_ex), missing="") - -``` - - -If there were available histological samples for all the fish individuals in the exchange event, no multiple modes are expected, however, if there were no hitology samples for some of them, multiple maturity stage modes could be found. In the next table this information is presented for the maturity staging. - - -```{r summary_multiple_modes_maturity_ex} - -histN=ad_long_ex %>% subset(TypeAnnotation=="eventOrganizer") %>% select(FishID) %>% unique() %>% dim() -allN=ad_long_ex %>% select(FishID) %>% unique() %>% dim() -perc_not_Hist=100-round(100*histN[1]/allN[1], 0) - -summary_multiple_modes_maturity_ex <- - c(NSample = length(unique(ad_long_ex$FishID)), - Nhist=histN[1], - Perc_not_Hist=perc_not_Hist, - PercMM_traditional=unique(ifelse(multimode_cases_tab_traditional_Maturity_ex$NModes_trad=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_traditional_Maturity_ex)/length(unique(ad_long_ex$SampleID)))*100, digits=0),"%"))), - PercMM_linear_weight=unique(ifelse(multimode_cases_tab_linear_Maturity_ex$NModes_linear=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_linear_Maturity_ex)/length(unique(ad_long_ex$SampleID)))*100, digits=0),"%"))), - PercMM_negexp_weight=unique(ifelse(multimode_cases_tab_negexp_Maturity_ex$NModes_negexp=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_negexp_Maturity_ex)/length(unique(ad_long_ex$SampleID)))*100, digits=0),"%"))), - PercMM_multistage=unique(ifelse(multimode_cases_tab_multistage_Maturity_ex$NModes_multistage=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_multistage_Maturity_ex)/length(unique(ad_long_ex$SampleID)))*100, digits=0),"%")))) - -# Table caption -set.caption('**Table X:** Summary of statistics for maturity staging; Total number of fish individuals studied (NSample), number of fish individuals with histological samples (Nhist), percentage of fish individuals without histology (Perc_not_Hist). The percentage of cases (fish samples) with multiple modes depending on the approach to weight the experience of the stager which will be considered when defining the fish maturity stage mode. PercMM_traditional shows the percentage of the total samples for which multiple modes are obtained when all the stagers are equally weighted. PercMM_linear_weight shows the percentage of the total samples for which multiple modes are obtained when the weight assigned to the different stagers decreases linearly with the experience, while in the PercMM_negexp the weight applied decreases with a negative exponential shape with the experience. The PercMM_multistage shows the percentage of multiple mode cases when a combination of the different methodologies is used, as explained in the material and methods section') - -# Output table -pander(style_table0(summary_multiple_modes_maturity_ex), missing="") - - - -``` - - - -```{r summary_statistics_sex_ex} -summary_statistics_sex_ex <- - c(NSample = length(unique(ad_long_ex$SampleID)), - CU = cu_tab_sex_ex[nrow(cu_tab_sex_ex), "Total"], - PA = pa_tab_sex_ex[nrow(pa_tab_sex_ex), "Total"]) - -summary_statistics_maturity_ex <- - c(NSample = length(unique(ad_long_ex$SampleID)), - CU = cu_tab_maturity_ex[nrow(cu_tab_maturity_ex), "Total"], - PA = pa_tab_maturity_ex[nrow(pa_tab_maturity_ex), "Total"]) -``` - -The average percentage agreement by modal sex category for all stagers was `r summary_statistics_sex_ex[3]`, with a weighted average CU of `r summary_statistics_sex_ex[2]`. Regarding the maturation staging, the percentage agreement by modal maturation stage was `r summary_statistics_maturity_ex[3]`, and the weight average CU was `r summary_statistics_maturity_ex[2]` - - -***Sex categorization table*** - -```{r num_read_sex_ex} - -# Sex categorization TABLE - Expert stagers ################################################# - -data=num_read_tab_modal_sex_by_reader_ex -#data=sex_composition_tab_ex -stagers=stager_data[stager_data$Expertise=="Advanced",] -nstagers=length(unique(stagers$'Reader code')) -ncols=dim(data[,!colnames(data) %in% c("Total")])[2] -basicols=ncols-nstagers - -z=8 -nfig=ceiling((dim(data)[2]-basicols)/z) -N=min(z,ncols) -if(nfig==0){ -# Table caption - set.caption('**Table X:** Sex categorization table: presents the - number of categorizations made per expert for each modal sex category.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):dim(data)[2])] - pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt -} else { - - for(i in 1:nfig) -{ - if(i==1){ -# Table caption - set.caption('**Table X:** Sex categorization table: presents the - number of categorizations made per expert for each modal sex category.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] - pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt - } else { if(i -The percentage agreement per stager and modal maturity stage represent the proportion of the total number of stagings that are equal to the modal maturity stage. The weighted mean including at the bottom of the table is weighted according to number of maturity stagings. - -```{r percentage_agreement_sex_ex} - -# PERCENTAGE AGREEMENT TABLE by modal sex category - Expert stagers #################################### - -data=pa_tab_sex_ex -stagers=stager_data[stager_data$Expertise=="Advanced",] -nstagers=length(unique(stagers$'Reader code')) -ncols=dim(data[,!colnames(data) %in% c("Total")])[2] -basicols=ncols-nstagers - -z=8 -nfig=round((dim(data)[2]-basicols)/z, digits=0) -N=min(z,ncols) - -if(nfig==0){ - # Table caption - set.caption('**Table X:** Percentage agreement (PA) table: shows the PA per modal sex category and stager, the PA for all stagers combined and the weighted mean of the PA per stager.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):dim(data)[2])] - pander(style_table1(selec), missing = "-") -} else { - -for(i in 1:nfig) -{ - if(i==1){ - # Table caption - set.caption('**Table X:** Percentage agreement (PA) table: shows the PA per modal sex category and stager, the PA for all stagers combined and the weighted mean of the PA per stager.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] - pander(style_table1(selec), missing = "-") - } else { if(i -The frequency bias is shown by sex category for each modal sex category per stager. As for the previous tables, a combined frequency distribution for all stagers are calculated. This frequency distribution presented in the form of a table or a buble plot may be used to explore the bias in the sex determination by the stagers participating in this event. - -```{r freq_sex_ex} - -# frequency bias TABLE by sex category- Expert stagers ########################################### - -data=stager_bias_freq_tab_sex_ex -stagers=stager_data[stager_data$Expertise=="Advanced",] -nstagers=length(unique(stagers$'Reader code')) -ncols=dim(data[,!colnames(data) %in% c("Total")])[2] -basicols=ncols-nstagers - -z=ncol(data)-2 -nfig=round((dim(data)[2]-basicols)/z, digits=0) -N=min(z,ncols) - -if(nfig==0){ - # Table caption - set.caption('**Table X:** Frequency bias table: represents the frequency per - modal sex category per stager and the frequency distribution of all stagers combined - per modal sex category.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):dim(data)[2])] - pander(style_table0(selec), missing = "-") -} else { - -for(i in 1:nfig) -{ - if(i==1){ - # Table caption - set.caption('**Table X:** Frequency bias table: represents the frequency per - modal sex category per stager and the frequency distribution of all stagers combined - per modal sex category.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] - pander(style_table1(selec), missing = "-") - } else { if(i -The frequency bias is shown by maturity stages for each modal maturity stage per stager. As for the previous tables, a combined frequency distribution for all stagers are calculated. This frequency distribution presented in the form of a table or a buble plot may be used to explore the bias in the maturity determination by the stagers participating in this event. - -```{r freq_maturity_ex} - -# frequency bias TABLE by maturity stage - Expert stagers ########################################### - -data=stager_bias_freq_tab_maturity_ex -stagers=stager_data[stager_data$Expertise=="Advanced",] -nstagers=length(unique(stagers$'Reader code')) -ncols=dim(data[,!colnames(data) %in% c("Total")])[2] -basicols=ncols-nstagers - -z=ncol(data)-2 -nfig=round((dim(data)[2]-basicols)/z, digits=0) -N=min(z,ncols) - -if(nfig==0){ - # Table caption - set.caption('**Table X:** Frequency bias table: represents the frequency per - modal maturity stage per stager and the frequency distribution of all stagers combined - per modal maturity stage.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):dim(data)[2])] - pander(style_table0(selec), missing = "-") -} else { - -for(i in 1:nfig) -{ - if(i==1){ - # Table caption - set.caption('**Table X:** Frequency bias table: represents the frequency per - modal maturity stage per stager and the frequency distribution of all stagers combined - per modal maturity stage.') - selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] - pander(style_table1(selec), missing = "-") - } else { if(i -The Sex Category Error Matrix is calculated per area, and only based on the maturity stagings of the advanced stagers. - -```{r sex_category_error_matrix_ex, results='asis'} - -# SSEM - Advanced stagers ################################################# - -# Loop through each area and output MSEM for that area -for (i in seq_along(ssem_ex)) { - # name - strata_name <- names(ssem_ex)[i] - - # Table caption - cap_i <- paste0("**Table X:** Sex Category Error Matrix (SSEM) for ", strata_name ,". - The SSEM shows the proportional distribution - of age category for each modal age category. Stage column should sum to - one but due to rounding there might be small deviations in - some cases. Only advanced stagers are used for calculating the SSEM.") - - set.caption(cap_i) - pander(style_table3(ssem_ex[[i]]), missing = "-") - -} -``` - - -***Maturity Stage Error Matrix (MSME)*** - - -The Maturity Stage Error Matrix is calculated per area, and only based on the maturity stagings of the advanced stagers. - -```{r maturity_stage_error_matrix_ex, results='asis'} - -# MSEM - Advanced stagers ################################################# - -# Loop through each area and output MSEM for that area -for (i in seq_along(msem_ex)) { - # name - strata_name <- names(msem_ex)[i] - - # Table caption - cap_i <- paste0("**Table X:** Maturity stage error matrix (MSEM) for ", strata_name ,". - The MSEM shows the proportional distribution - of maturity stagings for each modal maturity stage. Stage column should sum to - one but due to rounding there might be small deviations in - some cases. Only advanced stagers are used for calculating the MSEM.") - - set.caption(cap_i) - pander(style_table3(msem_ex[[i]]), missing = "-") - -} -``` - - - - - -#### With samples split by strata - -```{r set_strata_ex} -# initialise strata loop -istrata <- 0 -group <- "ex" -print_strata <- function() length(params$strata) >= istrata -``` - - - -```{r ref.label='strata_plus_one'} -``` - -```{r ref.label='strata_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_by', eval = print_strata()} -``` - -```{r ref.label='ar_matur_title', eval = print_strata()} -``` - -```{r ref.label='ar_matur_by', eval = print_strata()} -``` - -```{r ref.label='cu_sex_title', eval = print_strata()} -``` - -```{r ref.label='cu_sex_by', eval = print_strata()} -``` - -```{r ref.label='cu_matur_title', eval = print_strata()} -``` - -```{r ref.label='cu_matur_by', eval = print_strata()} -``` - -```{r ref.label='pa_sex_title', eval = print_strata()} -``` - -```{r ref.label='pa_sex_by', eval = print_strata()} -``` - -```{r ref.label='pa_matur_title', eval = print_strata()} -``` - -```{r ref.label='pa_matur_by', eval = print_strata()} -``` - - - - -```{r ref.label='strata_plus_one'} -``` - -```{r ref.label='strata_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_by', eval = print_strata()} -``` - -```{r ref.label='ar_matur_title', eval = print_strata()} -``` - -```{r ref.label='ar_matur_by', eval = print_strata()} -``` - -```{r ref.label='cu_sex_title', eval = print_strata()} -``` - -```{r ref.label='cu_sex_by', eval = print_strata()} -``` - -```{r ref.label='cu_matur_title', eval = print_strata()} -``` - -```{r ref.label='cu_matur_by', eval = print_strata()} -``` - -```{r ref.label='pa_sex_title', eval = print_strata()} -``` - -```{r ref.label='pa_sex_by', eval = print_strata()} -``` - -```{r ref.label='pa_matur_title', eval = print_strata()} -``` - -```{r ref.label='pa_matur_by', eval = print_strata()} -``` - - - - -```{r ref.label='strata_plus_one'} -``` - -```{r ref.label='strata_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_by', eval = print_strata()} -``` - -```{r ref.label='ar_matur_title', eval = print_strata()} -``` - -```{r ref.label='ar_matur_by', eval = print_strata()} -``` - -```{r ref.label='cu_sex_title', eval = print_strata()} -``` - -```{r ref.label='cu_sex_by', eval = print_strata()} -``` - -```{r ref.label='cu_matur_title', eval = print_strata()} -``` - -```{r ref.label='cu_matur_by', eval = print_strata()} -``` - -```{r ref.label='pa_sex_title', eval = print_strata()} -``` - -```{r ref.label='pa_sex_by', eval = print_strata()} -``` - -```{r ref.label='pa_matur_title', eval = print_strata()} -``` - -```{r ref.label='pa_matur_by', eval = print_strata()} -``` - - - - -```{r ref.label='strata_plus_one'} -``` - -```{r ref.label='strata_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_title', eval = print_strata()} -``` - -```{r ref.label='ar_sex_by', eval = print_strata()} -``` - -```{r ref.label='ar_matur_title', eval = print_strata()} -``` - -```{r ref.label='ar_matur_by', eval = print_strata()} -``` - -```{r ref.label='cu_sex_title', eval = print_strata()} -``` - -```{r ref.label='cu_sex_by', eval = print_strata()} -``` - -```{r ref.label='cu_matur_title', eval = print_strata()} -``` - -```{r ref.label='cu_matur_by', eval = print_strata()} -``` - -```{r ref.label='pa_sex_title', eval = print_strata()} -``` - -```{r ref.label='pa_sex_by', eval = print_strata()} -``` - -```{r ref.label='pa_matur_title', eval = print_strata()} -``` - -```{r ref.label='pa_matur_by', eval = print_strata()} -``` - - - - -## Discussion - - - - - - -## Conclusion - - - - - - - - - -# Other ToRs - - - - - - - - -# References - -Gary D. Kader & Mike Perry (2007) Variability for Categorical Variables, -Journal of Statistics Education, 15:2, , DOI: 10.1080/10691898.2007.11889465 - -ICES (2014) Report of the Workshop on Statistical Analysis of Biological Calibration Studies (WKSABCAL). ICES CM 2014/ACOM: 35 - - - - -# Annex 1. Agenda - - -# Annex 2. List of participants - -```{r participants} - -# COMPLETE PARTICIPANTS LIST ################################################## - -# Table caption -set.caption('**Table X:** Participants list.') - -# Output table -pander(style_table0(stager_data), missing = "-") - -``` - - - -# Annex 3. Additional results - -***Data Summary Statistics*** - -```{r tab1_sum_stat_sex} - -# DATA SUMMARY STATISTICS for the sex category annotations- All stagers ####################################### - -# Table caption -set.caption('**Table X:** Summary of statistics; Number of samples, CU and PA (%) for the sex category annotations') - -# Output table -pander(style_table0(summary_statistics_sex_all), missing="") - -``` - -```{r tab1_sum_stat_maturity} - -# DATA SUMMARY STATISTICS for the maturity stage annotations- All stagers ####################################### - -# Table caption -set.caption('**Table X:** Summary of statistics; Number of samples, CU and PA (%) for the maturity stage annotations') - -# Output table -pander(style_table0(summary_statistics_maturity_all), missing="") - -``` - - -## Results all stagers - -**All samples included** - -```{r set_annex_all} -group <- "all" -``` - -***Data overview*** - -```{r data_overview_maturity} -# DATA OVERVIEW - All stagers ################################################# - -# Table caption -set.caption('**Table X:** Data overview including modal - maturity and statistics per sample.') - -# Output table -pander(style_table0(get(vname("data_overview_tab_maturity"))), missing="-", style = "simple") -``` - - -```{r data_overview_sex} -# DATA OVERVIEW - All stagers ################################################# - -# Table caption -set.caption('**Table X:** Data overview including modal - sex and statistics per sample.') - -# Output table -pander(style_table0(get(vname("data_overview_tab_sex"))), missing="-", style = "simple") -``` - - - -***Individual frequency plots by sex categorization*** - -```{r freq_bias_plots_sex, fig.width = 5, fig.height = 4, fig.cap = cap_in} -# INDIVIDUAL Frequency PLOTS - All stagers ######################################### - -# Figure caption -cap_in <- '**Figure X:** Frequency plots per stager. The modal sex category is plotted against the different sex categories. The frequency of categies corresponging to each modal sex category is represented together with the 1:1 equilibrium line (solid line).' - -# Output figures -plot_bias_sex(get(vname("ad_long"))) -``` - - -***Individual frequency plots by maturity stage*** - -```{r freq_bias_plots_maturity, fig.width = 5, fig.height = 4, fig.cap = cap_in} -# INDIVIDUAL Frequency PLOTS - All stagers ######################################### - -# Figure caption -cap_in <- '**Figure X:** Frequency plots per stager. The modal maturity stage is plotted against the different maturity stages. The frequency of stagings corresponging to each modal maturity/maturity combinations is represented together with the 1:1 equilibrium line (solid line).' - -# Output figures -plot_bias_matur(get(vname("ad_long"))) -``` - - -***Plot with the CU AND PA statistichs for maturity stagings*** - -```{r stat_plot_maturity, fig.width = 7, fig.height = 4, fig.cap = cap_in} -# PLOT STATISTICS, CU AND PA - All stagers ############################# - -# Figure caption -cap_in <- '**Figure X:** CU and PA are plotted against modal maturity' - -# The coefficient of variation (CV%), percent agreement and -# the standard deviation (STDEV) are plotted against modal maturity - -# Output figure -plot_stat_matur(get(vname("ad_long"))) -``` - - -***Plot with the CU AND PA statistichs for sex categorizations*** - -```{r stat_plot_sex, fig.width = 7, fig.height = 4, fig.cap = cap_in} -# PLOT STATISTICS, CU AND PA - All stagers ############################# - -# Figure caption -cap_in <- '**Figure X:** CU and PA are plotted against modal sex' - -# The coefficient of variation (CV%), percent agreement and -# the standard deviation (STDEV) are plotted against modal sex - -# Output figure -plot_stat_sex(get(vname("ad_long"))) -``` - - - -## Results Advanced stagers - -```{r set_annex_ex} -group <- "ex" -``` - -**All samples included** - -***Data overview*** - -```{r ref.label='data_overview_sex'} -``` - -```{r ref.label='data_overview_maturity'} -``` - -***Individual frequency plots by sex categorization*** - -```{r ref.label='freq_bias_plots_sex', fig.width = 5, fig.height = 4, fig.cap = cap_in} -``` - -***Individual frequency plots by maturity stage*** - -```{r ref.label='freq_bias_plots_maturity', fig.width = 5, fig.height = 4, fig.cap = cap_in} -``` - -***Plot with the CU AND PA statistichs for maturity stagings*** - -```{r ref.label='stat_plot_maturity'} -``` - -***Plot with the CU AND PA statistichs for sex categorizations*** - -```{r ref.label='stat_plot_sex_ex'} -``` - - -# Annex 4. Table of the cases with multiple modes - -## All stagers - -***Multiple mode in the maturity stage*** - -```{r multiple mode cases maturity all stagers} -datos=ad_long_all[ad_long_all$NModes_trad_Maturity>1, c(1, 2, 4, 5, 8, 9, 10, 11, 12, 13, 16:26, 35, 36)] - -# Output table -pander(style_table0(datos), missing="-", style = "simple") - - -``` - -***Multiple mode in the sex category*** - -```{r multiple mode cases sex all stagers} -datos=ad_long_all[ad_long_all$NModes_trad_Sex>1, c(1, 2, 4, 5, 8, 9, 10, 11, 12, 13, 16:26, 56, 57)] - -# Output table -pander(style_table0(datos), missing="-", style = "simple") - - -``` - -## Advanced stagers - -***Multiple mode in the maturity stage*** - -```{r multiple mode cases maturity advanced stagers} -datos=ad_long_ex[ad_long_ex$NModes_trad_Maturity>1, c(1, 2, 4, 5, 8, 9, 10, 11, 12, 13, 16:26, 35, 36)] - -# Output table -pander(style_table0(datos), missing="-", style = "simple") - - -``` - -***Multiple mode in the sex category*** - -```{r multiple mode cases sex advanced stagers} -datos=ad_long_ex[ad_long_ex$NModes_trad_Sex>1, c(1, 2, 4, 5, 8, 9, 10, 11, 12, 13, 16:26, 56, 57)] - -# Output table -pander(style_table0(datos), missing="-", style = "simple") - - -``` - - -# Annex 5. ToRs for next meeting - - - -# Annex 6. Recommendations - - - -# Annex 7. Report specific annexes +--- +output: + word_document: + fig_caption: true + fig_height: 10 + fig_width: 10 + reference_docx: bootstrap/initial/data/reportTemplate.docx + toc: true + keep_md: false +params: + report_title: "" + strata: NULL +--- + +```{r introduction, include = FALSE} + +# INTRODUCTION ################################################################ + +# This markdown documents and integrated scripts analyse biological +# maturity readings. +# The output is a .docx template that includes +# the results of the analysis and should be used as a standard for +# reporting of maturity staging comparisons. + +``` + +```{r chunk_setup, include=FALSE} +# CHUNK SETUPS ################################################################# + +knitr::opts_chunk$set(echo = FALSE, warning = FALSE, + message=FALSE, results = 'asis', dpi=400) +``` + +```{r pander_settings, include = FALSE} + +# PANDER OPTIONS ############################################################## + +panderOptions('table.split.table', Inf) +panderOptions('keep.trailing.zeros', TRUE) +panderOptions('table.alignment.default', "center") + +``` + +--- +title: `r config$report_title` +--- + +# Executive summary + +# Terms of reference + +# Agenda and participant list + +The agenda can be found in Annex 1 and the list of participants in Annex 2. + + +# Introduction + +This part should include a background to the species, the workshop/exchange +and what to expect to read about in the report. + + +# Methods + +This report contains statistical analyses and comparisons of sex categorization and maturity staging’s +in the form of tables and graphical plots. + +First, an overview of participating maturity stagers and the samples are presented. + +Before each table or plot there is a short explanation of it. This text is +thought as a help to understand the tables/plot and can just be deleted in +the final output report. The document can be edited just like any other +.docx file. New text can be added, additional pictures can be included and +the tables edited. If some tables which are presently in the annexes need to +be moved to the body of the report this is also possible. Only the plots +cannot be changed. + +In the first part of analysis some of the tables and plots from the +Guus Eltink Excel sheet 'Age Reading Comparisons' **(Eltink, A.T.G.W. 2000)** are presented, that can be used for "maturity staging comparisons". Since the sex categories and maturity stages are categorical in comparison to ages that are quantitative, some of the statistics presented in the report of the age exchange events cannot be calculated here, like the average percentage error or the coefficient of variation. The latter is being replaced by the coefficient of unlikeability (see below). + + +**Percentage Agreement** + +The percentage agreement tells how large is the part of sex categorizations/maturity staging’s that are equal to the modal sex/maturity. The percentage agreement is estimated by modal sex or maturity and stager as the proportion (as percentage) of times that the lectures of that stager agreed with the resulting modal sex or maturity. This percentage is estimated as the number of times that a stager agreed with the modal sex or maturity divided by the total number of gonads analyzed by a stager for that modal maturity stage. + +$$PA = { \frac{ {number \,of \,readings \,that \,agree \,with \,modal \,maturity}} {total \,number \,of \,readings \,by \,modal \,maturity} } \cdot {100 \%} $$ + + +**Coefficient of unalikeability (CU)** + +The concept of unalikeability (Kader and Perry, 2007) focuses on how often observations differ within a group. Specifically, for the sex/maturity staging events, the CU provides a measure of how alike, for each modal maturity stage, the stages decided by each stager were (or all stagers at once). The cu ranges between 0 and 1. The higher the cu value, the more unalike the data are. + +The table presents the CU per modal maturity stage and stager. For the case of a finite number of observations (n), a finite number of categories (m) and a finite number of objects, k_i, within category i, will allow expression of the coefficient of unalikeablity as: + +$$ u = {1 - ∑ p_i ^ 2} $$ + +where + +$$ p_i = {k_i / n} $$ + +To the table is also added the CU of all stagers combined per modal maturity and +a weighted mean of the CU per stager. + + +**Maturity stage error matrix (MSEM)** + +Maturity stage error matrices (AEM) were produced following the same procedures outlined + by WKSABCAL (2014) for calculating the "Age error matrix". MSEM shows the proportion of each modal +maturity stage mis-staged as other stages. The sum of each row is 1, which equals 100%. +The maturity data was analysed twice, the first time all stagers were included +and the second time only the “advanced” stagers were included. If a stager +is “advanced” then they are considered well trained and they provide maturity stages +for stock assessment or similar purposes. When the MSEM is compiled for +assessment purposes it uses only those stagers who provide maturity data for +the stock assessment in that specific area. +In this report the MSEM presented represents is in the transpose format. +The frequency bias table represents the MSEM in the format proposed by WKSABCAL (2014). + + + +# Analysis of maturity staging calibration exercise (ToR?) + + +## Overview of samples and stagers + + +\br + +```{r sample_overview} + +# Table caption +set.caption( + paste( + '**Table X:** Overview of samples used for the exchange event number ', + config$event_id, + '. The number of samples (num) are shown by year, ices area, season (qtr). The range of maturity stages (mat_stages) and sex (sex_cat) assigned by the group or stagers participating in the event is also shown, as well as the range of modal maturity stages (modal_mat_stages) and modal sex (modal_sex_cat) decided for all samples.', + sep="" + ) +) + +# Output table +pander(style_table0(sample_data_overview), style = "simple") + +``` + +\br + +```{r participants_overview} + +# PARTICIPATANTS OVERVIEW ##################################################### + +# Table caption +set.caption('**Table X:** Overview of stagers participating in the event, with their overall expertise (Expertise), and their ranking position based in their experience (only valid if the experience weighting protocol has been applied).') + +# Output table +pander(style_table0(filter(stager_data))) + +``` + + +## Results + +Text? + +### All stagers + +Text? + +#### All samples together + +Text? + +**Multimodal cases** + +Those writing the report put TEXT here describing the results. + +If there were available histological samples for all the fish individuals in the exchange event, no multiple modes are expected, however, if there were no histology samples for some of them, multiple maturity stage modes could be found. In the next table this information is presented for the sex category determination. + + + +```{r summary_multiple_modes_sex_all} + +#histN=ad_long_all %>% subset(TypeAnnotation=="eventOrganizer") %>% select(FishID) %>% unique() %>% dim() +#allN=ad_long_all %>% select(FishID) %>% unique() %>% dim() +#perc_not_Hist=round(100*histN[1]/allN[1], 0) + +histN=ad %>% subset(DoesSampleHaveHistologyImage=="Yes") %>% select(FishID) %>% unique() %>% dim() +allN=ad_long_all %>% select(FishID) %>% unique() %>% dim() +perc_not_Hist=100-round(100*histN[1]/allN[1], 0) + +summary_multiple_modes_sex_all <- + c(NSample = length(unique(ad_long_all$FishID)), + Nhist=histN[1], + Perc_not_Hist=perc_not_Hist, + PercMM_traditional=unique(ifelse(multimode_cases_tab_traditional_Sex_all$NModes_trad=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_traditional_Sex_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), + PercMM_linear_weight=unique(ifelse(multimode_cases_tab_linear_Sex_all$NModes_linear=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_linear_Sex_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), + PercMM_negexp_weight=unique(ifelse(multimode_cases_tab_negexp_Sex_all$NModes_negexp=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_negexp_Sex_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), + PercMM_multistage=unique(ifelse(multimode_cases_tab_multistage_Sex_all$NModes_multistage=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_multistage_Sex_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%")))) + +# Table caption +set.caption('**Table X:** Summary of statistics for sex staging; Total number of fish individuals studied (NSample), number of fish individuals with histological samples (Nhist), percentage of fish individuals without histology (Perc_not_Hist). The percentage of cases (fish samples) with multiple modes depending on the approach to weight the experience of the stager which will be considered when defining the fish sex stage mode. PercMM_traditional shows the percentage of the total samples for which multiple modes are obtained when all the stagers are equally weighted. PercMM_linear_weight shows the percentage of the total samples for which multiple modes are obtained when the weight assigned to the different stagers decreases linearly with the experience, while in the PercMM_negexp the weight applied decreases with a negative exponential shape with the experience. The PercMM_multistage shows the percentage of multiple mode cases when a combination of the different methodologies is used, as explained in the material and methods section') + +# Output table +pander(style_table0(summary_multiple_modes_sex_all), missing="") + +``` + + +If there were available histological samples for all the fish individuals in the exchange event, no multiple modes are expected, however, if there were no histology samples for some of them, multiple maturity stage modes could be found. In the next table this information is presented for the maturity staging. + + + +```{r summary_multiple_modes_maturity_all} + +histN=ad %>% subset(DoesSampleHaveHistologyImage=="Yes") %>% select(FishID) %>% unique() %>% dim() +allN=ad_long_all %>% select(FishID) %>% unique() %>% dim() +perc_not_Hist=round(100*histN[1]/allN[1], 0) + +summary_multiple_modes_maturity_all <- + c(NSample = length(unique(ad_long_all$FishID)), + Nhist=histN[1], + Perc_not_Hist=perc_not_Hist, + PercMM_traditional=unique(ifelse(multimode_cases_tab_traditional_Maturity_all$NModes_trad=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_traditional_Maturity_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), + PercMM_linear_weight=unique(ifelse(multimode_cases_tab_linear_Maturity_all$NModes_linear=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_linear_Maturity_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), + PercMM_negexp_weight=unique(ifelse(multimode_cases_tab_negexp_Maturity_all$NModes_negexp=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_negexp_Maturity_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%"))), + PercMM_multistage=unique(ifelse(multimode_cases_tab_multistage_Maturity_all$NModes_multistage=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_multistage_Maturity_all)/length(unique(ad_long_all$SampleID)))*100, digits=0),"%")))) + +# Table caption +set.caption('**Table X:** Summary of statistics for maturity staging; Total number of fish individuals studied (NSample), number of fish individuals with histological samples (Nhist), percentage of fish individuals without histology (Perc_not_Hist). The percentage of cases (fish samples) with multiple modes depending on the approach to weight the experience of the stager which will be considered when defining the fish maturity stage mode. PercMM_traditional shows the percentage of the total samples for which multiple modes are obtained when all the stagers are equally weighted. PercMM_linear_weight shows the percentage of the total samples for which multiple modes are obtained when the weight assigned to the different stagers decreases linearly with the experience, while in the PercMM_negexp the weight applied decreases with a negative exponential shape with the experience. The PercMM_multistage shows the percentage of multiple mode cases when a combination of the different methodologies is used, as explained in the material and methods section') + +# Output table +pander(style_table0(summary_multiple_modes_maturity_all), missing="") + + + +``` + + + +```{r summary_statistics_sex_all} +summary_statistics_sex_all <- + c(NSample = length(unique(ad_long_all$FishID)), + CU = cu_tab_sex_all[nrow(cu_tab_sex_all), "Total"], + PA = pa_tab_sex_all[nrow(pa_tab_sex_all), "Total"]) + +summary_statistics_maturity_all <- + c(NSample = length(unique(ad_long_all$FishID)), + CU = cu_tab_maturity_all[nrow(cu_tab_maturity_all), "Total"], + PA = pa_tab_maturity_all[nrow(pa_tab_maturity_all), "Total"]) +``` + +The average percentage agreement by modal sex category for all stagers was `r summary_statistics_sex_all[3]`, with a weighted average CU of `r summary_statistics_sex_all[2]`. Regarding the maturation staging, the percentage agreement by modal maturation stage was `r summary_statistics_maturity_all[3]`, and the weight average CU was `r summary_statistics_maturity_all[2]`. + + +**Sex categorization table** + +Text? + +```{r sex_composition_tab_all} + +# Sex categorization TABLE - All stagers ################################################# +data=sex_composition_tab_all +nstagers=length(unique(stager_data$'Reader code')) +ncols=dim(data[,!colnames(data) %in% c("Total")])[2] +basicols=ncols-nstagers + +z=8 +nfig=round((dim(data)[2]-basicols)/z) +N=min(z,ncols) +for(i in 1:nfig) +{ + if(i==1){ + # Table caption + set.caption('**Table X:** Sex categorization table: presents the + number of categorizations made per expert for each modal sex category.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] + pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt + } else { if(i +The percentage agreement per stager and modal maturity stage represent the proportion of the total number of staging’s that are equal to the modal maturity stage. The weighted mean including at the bottom of the table is weighted according to number of maturity staging’s. + +```{r percentage_agreement_sex_all} + +# PERCENTAGE AGREEMENT TABLE by modal sex category - All stagers #################################### + +data=pa_tab_sex_all +nstagers=length(unique(stager_data$'Reader code')) +ncols=dim(data[,!colnames(data) %in% c("Total")])[2] +basicols=ncols-nstagers + +z=8 +nfig=round((dim(data)[2]-basicols)/z) +N=min(z,ncols) +for(i in 1:nfig) +{ + if(i==1){ + # Table caption + set.caption('**Table X:** Percentage agreement (PA) table: shows the PA per modal sex category and stager, the PA for all stagers combined and the weighted mean of the PA per stager.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] + pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt + } else { if(i +The frequency bias is shown by sex category for each modal sex category per stager. As for the previous tables, a combined frequency distribution for all stagers are calculated. This frequency distribution presented in the form of a table or a bubble plot may be used to explore the bias in the sex determination by the stagers participating in this event. + +```{r freq_sex_all} + +# frequency bias TABLE by sex category- per stagers ########################################### + +data=stager_bias_freq_tab_sex_all +nstagers=length(unique(stager_data$'Reader code')) +ncols=dim(data[,!colnames(data) %in% c("Total")])[2] +basicols=ncols-nstagers + +z=8 +nfig=round((dim(data)[2]-basicols)/z) +N=min(z,ncols) +for(i in 1:nfig) +{ + if(i==1){ + # Table caption + set.caption('**Table X:** Frequency bias table: represents the frequency per + modal sex category per stager and the frequency distribution of all stagers combined + per modal sex category.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] + pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt + } else { if(i +The frequency bias is shown by maturity stages for each modal maturity stage per stager. As for the previous tables, a combined frequency distribution for all stagers are calculated. This frequency distribution presented in the form of a table or a bubble plot may be used to explore the bias in the maturity determination by the stagers participating in this event. + +```{r freq_maturity_all} + +# frequency bias TABLE by maturity stage - per stagers ########################################### + +data=stager_bias_freq_tab_maturity_all +nstagers=length(unique(stager_data$'Reader code')) +ncols=dim(data[,!colnames(data) %in% c("Total")])[2] +basicols=ncols-nstagers + +z=8 +nfig=round((dim(data)[2]-basicols)/z) +N=min(z,ncols) +for(i in 1:nfig) +{ + if(i==1){ + # Table caption + set.caption('**Table X:** Frequency bias table: represents the frequency per + modal maturity stage per stager and the frequency distribution of all stagers combined + per modal maturity stage.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] + pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt + } else { if(i= istrata +``` + + + +```{r strata_plus_one} +# second strata +istrata <- istrata + 1 +``` + +```{r strata_title, eval = print_strata()} +stratum <- params$strata[istrata] +# Section title +asis_output(paste0("**Results by ", stratum, "**")) +``` + + +```{r ar_sex_title, eval = print_strata()} +# title +asis_output(paste0("**Sex categorization by ", stratum, "**")) +``` + +```{r ar_sex_by, eval = print_strata()} +# NUMBER OF SEX CATEGORIZATION PER Strata - Advanced stagers ######################### + +# Table caption +set.caption(paste0('**Table X:** Number of sex categorization per ', stratum, ' for all the stagers.')) + +# Output table +pander(style_table1(get(vsname("num_read_tab_modal_sex_by"))), missing = "-") +``` + + +```{r ar_matur_title, eval = print_strata()} +# title +asis_output(paste0("**Maturity staging by ", stratum, "**")) +``` + +```{r ar_matur_by, eval = print_strata()} +# NUMBER OF MATURITY STAGING PER strata - Advanced stagers ######################### + +# Table caption +set.caption(paste0('**Table X:** Number of maturity staging’s per ', stratum, ' for all the stagers.')) + +# Output table +pander(style_table1(get(vsname("num_read_tab_modal_matur_by"))), missing = "-") +``` + + + +```{r cu_sex_title, eval = print_strata()} +# title +asis_output(paste0("**Coefficient of Unalikeability by modal sex category per", stratum, "**")) +``` + + +```{r cu_sex_by, eval = print_strata()} + +# Table caption +set.caption(paste0('**Table X:** CU by modal sex category per ', stratum, '.')) + +# Output table +pander(style_table2(get(vsname("cu_tab_sex_by"))), missing = "-") + +``` + + +```{r cu_matur_title, eval = print_strata()} +# title +asis_output(paste0("**Coefficient of Unalikeability by modal maturity stage per", stratum, "**")) +``` + + +```{r cu_matur_by, eval = print_strata()} + +# Table caption +set.caption(paste0('**Table X:** CU by modal maturity stage per ', stratum, '.')) + +# Output table +pander(style_table2(get(vsname("cu_tab_maturity_by"))), missing = "-") + +``` + + +```{r pa_sex_title, eval = print_strata()} +# title +asis_output(paste0("**Percentage Agreement by sex category per ", stratum, "**")) +``` + +```{r pa_sex_by, eval = print_strata()} +# Section title +#asis_("## PA by sex category - Advanced stagers") # # ############################################# + +# Table caption +set.caption(paste0('**Table X:** Percentage Agreement by sex category per ', stratum, '.')) + +# Output table +pander(style_table2(get(vsname("pa_tab_sex_by"))), missing = "-") + +``` + + +```{r pa_matur_title, eval = print_strata()} +# title +asis_output(paste0("**Percentage Agreement by maturation stage per ", stratum, "**")) +``` + +```{r pa_matur_by, eval = print_strata()} +# Section title +#asis_("## PA by maturation stage - Advanced stagers") # # ############################################# + +# Table caption +set.caption(paste0('**Table X:** Percentage Agreement by maturation stage per ', stratum, '.')) + +# Output table +pander(style_table2(get(vsname("pa_tab_maturity_by"))), missing = "-") + +``` + + + + +```{r ref.label='strata_plus_one'} +``` + +```{r ref.label='strata_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_by', eval = print_strata()} +``` + +```{r ref.label='ar_matur_title', eval = print_strata()} +``` + +```{r ref.label='ar_matur_by', eval = print_strata()} +``` + +```{r ref.label='cu_sex_title', eval = print_strata()} +``` + +```{r ref.label='cu_sex_by', eval = print_strata()} +``` + +```{r ref.label='cu_matur_title', eval = print_strata()} +``` + +```{r ref.label='cu_matur_by', eval = print_strata()} +``` + +```{r ref.label='pa_sex_title', eval = print_strata()} +``` + +```{r ref.label='pa_sex_by', eval = print_strata()} +``` + +```{r ref.label='pa_matur_title', eval = print_strata()} +``` + +```{r ref.label='pa_matur_by', eval = print_strata()} +``` + + + + +```{r ref.label='strata_plus_one'} +``` + +```{r ref.label='strata_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_by', eval = print_strata()} +``` + +```{r ref.label='ar_matur_title', eval = print_strata()} +``` + +```{r ref.label='ar_matur_by', eval = print_strata()} +``` + +```{r ref.label='cu_sex_title', eval = print_strata()} +``` + +```{r ref.label='cu_sex_by', eval = print_strata()} +``` + +```{r ref.label='cu_matur_title', eval = print_strata()} +``` + +```{r ref.label='cu_matur_by', eval = print_strata()} +``` + +```{r ref.label='pa_sex_title', eval = print_strata()} +``` + +```{r ref.label='pa_sex_by', eval = print_strata()} +``` + +```{r ref.label='pa_matur_title', eval = print_strata()} +``` + +```{r ref.label='pa_matur_by', eval = print_strata()} +``` + + + + +```{r ref.label='strata_plus_one'} +``` + +```{r ref.label='strata_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_by', eval = print_strata()} +``` + +```{r ref.label='ar_matur_title', eval = print_strata()} +``` + +```{r ref.label='ar_matur_by', eval = print_strata()} +``` + +```{r ref.label='cu_sex_title', eval = print_strata()} +``` + +```{r ref.label='cu_sex_by', eval = print_strata()} +``` + +```{r ref.label='cu_matur_title', eval = print_strata()} +``` + +```{r ref.label='cu_matur_by', eval = print_strata()} +``` + +```{r ref.label='pa_sex_title', eval = print_strata()} +``` + +```{r ref.label='pa_sex_by', eval = print_strata()} +``` + +```{r ref.label='pa_matur_title', eval = print_strata()} +``` + +```{r ref.label='pa_matur_by', eval = print_strata()} +``` + + + + + +## Advanced stagers + +Those writing the report put TEXT here describing the results. + +#### All samples together + +***Multimodal cases*** + +Those writing the report put TEXT here describing the results. + +If there were available histological samples for all the fish individuals in the exchange event, no multiple modes are expected, however, if there were no histology samples for some of them, multiple maturity stage modes could be found. In the next table this information is presented for the sex category determination. + + + +```{r summary_multiple_modes_sex_adv} +histN=ad %>% subset(DoesSampleHaveHistologyImage=="Yes") %>% select(FishID) %>% unique() %>% dim() +allN=ad_long_all %>% select(FishID) %>% unique() %>% dim() +perc_not_Hist=100-round(100*histN[1]/allN[1], 0) + +summary_multiple_modes_sex_adv <- + c(NSample = length(unique(ad_long_adv$FishID)), + Nhist=histN[1], + Perc_not_Hist=perc_not_Hist, + PercMM_traditional=unique(ifelse(multimode_cases_tab_traditional_Sex_adv$NModes_trad=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_traditional_Sex_adv)/length(unique(ad_long_adv$SampleID)))*100, digits=0),"%"))), + PercMM_linear_weight=unique(ifelse(multimode_cases_tab_linear_Sex_adv$NModes_linear=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_linear_Sex_adv)/length(unique(ad_long_adv$SampleID)))*100, digits=0),"%"))), + PercMM_negexp_weight=unique(ifelse(multimode_cases_tab_negexp_Sex_adv$NModes_negexp=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_negexp_Sex_adv)/length(unique(ad_long_adv$SampleID)))*100, digits=0),"%"))), + PercMM_multistage=unique(ifelse(multimode_cases_tab_multistage_Sex_adv$NModes_multistage=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_multistage_Sex_adv)/length(unique(ad_long_adv$SampleID)))*100, digits=0),"%")))) + +# Table caption +set.caption('**Table X:** Summary of statistics for sex staging; Total number of fish individuals studied (NSample), number of fish individuals with histological samples (Nhist), percentage of fish individuals without histology (Perc_not_Hist). The percentage of cases (fish samples) with multiple modes depending on the approach to weight the experience of the stager which will be considered when defining the fish sex stage mode. PercMM_traditional shows the percentage of the total samples for which multiple modes are obtained when all the stagers are equally weighted. PercMM_linear_weight shows the percentage of the total samples for which multiple modes are obtained when the weight assigned to the different stagers decreases linearly with the experience, while in the PercMM_negexp the weight applied decreases with a negative exponential shape with the experience. The PercMM_multistage shows the percentage of multiple mode cases when a combination of the different methodologies is used, as explained in the material and methods section') + +# Output table +pander(style_table0(summary_multiple_modes_sex_adv), missing="") + +``` + + +If there were available histological samples for all the fish individuals in the exchange event, no multiple modes are expected, however, if there were no histology samples for some of them, multiple maturity stage modes could be found. In the next table this information is presented for the maturity staging. + + +```{r summary_multiple_modes_maturity_adv} + +histN=ad %>% subset(DoesSampleHaveHistologyImage=="Yes") %>% select(FishID) %>% unique() %>% dim() +allN=ad_long_all %>% select(FishID) %>% unique() %>% dim() +perc_not_Hist=100-round(100*histN[1]/allN[1], 0) + +summary_multiple_modes_maturity_adv <- + c(NSample = length(unique(ad_long_adv$FishID)), + Nhist=histN[1], + Perc_not_Hist=perc_not_Hist, + PercMM_traditional=unique(ifelse(multimode_cases_tab_traditional_Maturity_adv$NModes_trad=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_traditional_Maturity_adv)/length(unique(ad_long_adv$SampleID)))*100, digits=0),"%"))), + PercMM_linear_weight=unique(ifelse(multimode_cases_tab_linear_Maturity_adv$NModes_linear=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_linear_Maturity_adv)/length(unique(ad_long_adv$SampleID)))*100, digits=0),"%"))), + PercMM_negexp_weight=unique(ifelse(multimode_cases_tab_negexp_Maturity_adv$NModes_negexp=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_negexp_Maturity_adv)/length(unique(ad_long_adv$SampleID)))*100, digits=0),"%"))), + PercMM_multistage=unique(ifelse(multimode_cases_tab_multistage_Maturity_adv$NModes_multistage=="zero", paste(0,"%"), paste(round((nrow(multimode_cases_tab_multistage_Maturity_adv)/length(unique(ad_long_adv$SampleID)))*100, digits=0),"%")))) + +# Table caption +set.caption('**Table X:** Summary of statistics for maturity staging; Total number of fish individuals studied (NSample), number of fish individuals with histological samples (Nhist), percentage of fish individuals without histology (Perc_not_Hist). The percentage of cases (fish samples) with multiple modes depending on the approach to weight the experience of the stager which will be considered when defining the fish maturity stage mode. PercMM_traditional shows the percentage of the total samples for which multiple modes are obtained when all the stagers are equally weighted. PercMM_linear_weight shows the percentage of the total samples for which multiple modes are obtained when the weight assigned to the different stagers decreases linearly with the experience, while in the PercMM_negexp the weight applied decreases with a negative exponential shape with the experience. The PercMM_multistage shows the percentage of multiple mode cases when a combination of the different methodologies is used, as explained in the material and methods section') + +# Output table +pander(style_table0(summary_multiple_modes_maturity_adv), missing="") + + + +``` + + + +```{r summary_statistics_sex_adv} +summary_statistics_sex_adv <- + c(NSample = length(unique(ad_long_adv$SampleID)), + CU = cu_tab_sex_adv[nrow(cu_tab_sex_adv), "Total"], + PA = pa_tab_sex_adv[nrow(pa_tab_sex_adv), "Total"]) + +summary_statistics_maturity_adv <- + c(NSample = length(unique(ad_long_adv$SampleID)), + CU = cu_tab_maturity_adv[nrow(cu_tab_maturity_adv), "Total"], + PA = pa_tab_maturity_adv[nrow(pa_tab_maturity_adv), "Total"]) +``` + +The average percentage agreement by modal sex category for all stagers was `r summary_statistics_sex_adv[3]`, with a weighted average CU of `r summary_statistics_sex_adv[2]`. Regarding the maturation staging, the percentage agreement by modal maturation stage was `r summary_statistics_maturity_adv[3]`, and the weight average CU was `r summary_statistics_maturity_adv[2]` + + +***Sex categorization table*** + +```{r sex_composition_tab_adv} + +# Sex categorization TABLE - Advanced stagers ################################################# + +data=sex_composition_tab_adv +stagers=stager_data[stager_data$Expertise=="Advanced",] +nstagers=length(unique(stagers$'Reader code')) +ncols=dim(data[,!colnames(data) %in% c("Total")])[2] +basicols=ncols-nstagers + +#z=8 +z=ncol(data)-2 +nfig=ceiling((dim(data)[2]-basicols)/z) +N=min(z,ncols) +if(nfig==0){ +# Table caption + set.caption('**Table X:** Sex categorization table: presents the + number of categorizations made per expert for each modal sex category.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):dim(data)[2])] + pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt +} else { + + for(i in 1:nfig) +{ + if(i==1){ +# Table caption + set.caption('**Table X:** Sex categorization table: presents the + number of categorizations made per expert for each modal sex category.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] + pander(style_table0(selec), missing = "-", style = "simple") #changed here to not appear the last column in bolt + } else { if(i +The percentage agreement per stager and modal maturity stage represent the proportion of the total number of staging’s that are equal to the modal maturity stage. The weighted mean including at the bottom of the table is weighted according to number of maturity staging’s. + +```{r percentage_agreement_sex_adv} + +# PERCENTAGE AGREEMENT TABLE by modal sex category - Advanced stagers #################################### + +data=pa_tab_sex_adv +stagers=stager_data[stager_data$Expertise=="Advanced",] +nstagers=length(unique(stagers$'Reader code')) +ncols=dim(data[,!colnames(data) %in% c("Total")])[2] +basicols=ncols-nstagers + +z=8 +nfig=round((dim(data)[2]-basicols)/z, digits=0) +N=min(z,ncols) + +if(nfig==0){ + # Table caption + set.caption('**Table X:** Percentage agreement (PA) table: shows the PA per modal sex category and stager, the PA for all stagers combined and the weighted mean of the PA per stager.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):dim(data)[2])] + pander(style_table1(selec), missing = "-") +} else { + +for(i in 1:nfig) +{ + if(i==1){ + # Table caption + set.caption('**Table X:** Percentage agreement (PA) table: shows the PA per modal sex category and stager, the PA for all stagers combined and the weighted mean of the PA per stager.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] + pander(style_table1(selec), missing = "-") + } else { if(i +The frequency bias is shown by sex category for each modal sex category per stager. As for the previous tables, a combined frequency distribution for all stagers are calculated. This frequency distribution presented in the form of a table or a bubble plot may be used to explore the bias in the sex determination by the stagers participating in this event. + +```{r freq_sex_adv} + +# frequency bias TABLE by sex category- Expert stagers ########################################### + +data=stager_bias_freq_tab_sex_adv +stagers=stager_data[stager_data$Expertise=="Advanced",] +nstagers=length(unique(stagers$'Reader code')) +ncols=dim(data[,!colnames(data) %in% c("Total")])[2] +basicols=ncols-nstagers + +#z=ncol(data)-2 +z=8 +nfig=round((dim(data)[2]-basicols)/z, digits=0) +N=min(z,ncols) + +if(nfig==0){ + # Table caption + set.caption('**Table X:** Frequency bias table: represents the frequency per + modal sex category per stager and the frequency distribution of all stagers combined + per modal sex category.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):dim(data)[2])] + pander(style_table1(selec), missing = "-") +} else { + +for(i in 1:nfig) +{ + if(i==1){ + # Table caption + set.caption('**Table X:** Frequency bias table: represents the frequency per + modal sex category per stager and the frequency distribution of all stagers combined + per modal sex category.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] + pander(style_table1(selec), missing = "-") + } else { if(i +The frequency bias is shown by maturity stages for each modal maturity stage per stager. As for the previous tables, a combined frequency distribution for all stagers are calculated. This frequency distribution presented in the form of a table or a bubble plot may be used to explore the bias in the maturity determination by the stagers participating in this event. + +```{r freq_maturity_adv} + +# frequency bias TABLE by maturity stage - Expert stagers ########################################### + +data=stager_bias_freq_tab_maturity_adv +stagers=stager_data[stager_data$Expertise=="Advanced",] +nstagers=length(unique(stagers$'Reader code')) +ncols=dim(data[,!colnames(data) %in% c("Total")])[2] +basicols=ncols-nstagers + +#z=ncol(data)-2 +z=8 +nfig=round((dim(data)[2]-basicols)/z, digits=0) +N=min(z,ncols) + +if(nfig==0){ + # Table caption + set.caption('**Table X:** Frequency bias table: represents the frequency per + modal maturity stage per stager and the frequency distribution of all stagers combined + per modal maturity stage.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):dim(data)[2])] + pander(style_table1(selec), missing = "-") +} else { + +for(i in 1:nfig) +{ + if(i==1){ + # Table caption + set.caption('**Table X:** Frequency bias table: represents the frequency per + modal maturity stage per stager and the frequency distribution of all stagers combined + per modal maturity stage.') + selec=data[,c(1:basicols, (i+basicols+((N-1)*(i-1))):((i+basicols+((N-1)*(i-1)))+(N-1)))] + pander(style_table1(selec), missing = "-") + } else { if(i +The Sex Category Error Matrix is calculated per area, and only based on the maturity staging’s of the advanced stagers. + +```{r sex_category_error_matrix_ex, results='asis'} + +# SSEM - Advanced stagers ################################################# + +# Loop through each area and output MSEM for that area +for (i in seq_along(ssem_adv)) { + # name + strata_name <- names(ssem_adv)[i] + + # Table caption + cap_i <- paste0("**Table X:** Sex Category Error Matrix (SSEM) for ", strata_name ,". + The SSEM shows the proportional distribution + of age category for each modal age category. Stage column should sum to + one but due to rounding there might be small deviations in + some cases. Only advanced stagers are used for calculating the SSEM.") + + set.caption(cap_i) + pander(style_table3(ssem_adv[[i]]), missing = "-") + +} +``` + + +***Maturity Stage Error Matrix (MSME)*** + + +The Maturity Stage Error Matrix is calculated per area, and only based on the maturity staging’s of the advanced stagers. + +```{r maturity_stage_error_matrix_ex, results='asis'} + +# MSEM - Advanced stagers ################################################# + +# Loop through each area and output MSEM for that area +for (i in seq_along(msem_adv)) { + # name + strata_name <- names(msem_adv)[i] + + # Table caption + cap_i <- paste0("**Table X:** Maturity stage error matrix (MSEM) for ", strata_name ,". + The MSEM shows the proportional distribution + of maturity staging’s for each modal maturity stage. Stage column should sum to + one but due to rounding there might be small deviations in + some cases. Only advanced stagers are used for calculating the MSEM.") + + set.caption(cap_i) + pander(style_table3(msem_adv[[i]]), missing = "-") + +} +``` + + + + + +#### With samples split by strata + +```{r set_strata_adv} +# initialise strata loop +istrata <- 0 +group <- "adv" +print_strata <- function() length(params$strata) >= istrata +``` + + + +```{r ref.label='strata_plus_one'} +``` + +```{r ref.label='strata_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_by', eval = print_strata()} +``` + +```{r ref.label='ar_matur_title', eval = print_strata()} +``` + +```{r ref.label='ar_matur_by', eval = print_strata()} +``` + +```{r ref.label='cu_sex_title', eval = print_strata()} +``` + +```{r ref.label='cu_sex_by', eval = print_strata()} +``` + +```{r ref.label='cu_matur_title', eval = print_strata()} +``` + +```{r ref.label='cu_matur_by', eval = print_strata()} +``` + +```{r ref.label='pa_sex_title', eval = print_strata()} +``` + +```{r ref.label='pa_sex_by', eval = print_strata()} +``` + +```{r ref.label='pa_matur_title', eval = print_strata()} +``` + +```{r ref.label='pa_matur_by', eval = print_strata()} +``` + + + + +```{r ref.label='strata_plus_one'} +``` + +```{r ref.label='strata_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_by', eval = print_strata()} +``` + +```{r ref.label='ar_matur_title', eval = print_strata()} +``` + +```{r ref.label='ar_matur_by', eval = print_strata()} +``` + +```{r ref.label='cu_sex_title', eval = print_strata()} +``` + +```{r ref.label='cu_sex_by', eval = print_strata()} +``` + +```{r ref.label='cu_matur_title', eval = print_strata()} +``` + +```{r ref.label='cu_matur_by', eval = print_strata()} +``` + +```{r ref.label='pa_sex_title', eval = print_strata()} +``` + +```{r ref.label='pa_sex_by', eval = print_strata()} +``` + +```{r ref.label='pa_matur_title', eval = print_strata()} +``` + +```{r ref.label='pa_matur_by', eval = print_strata()} +``` + + + + +```{r ref.label='strata_plus_one'} +``` + +```{r ref.label='strata_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_by', eval = print_strata()} +``` + +```{r ref.label='ar_matur_title', eval = print_strata()} +``` + +```{r ref.label='ar_matur_by', eval = print_strata()} +``` + +```{r ref.label='cu_sex_title', eval = print_strata()} +``` + +```{r ref.label='cu_sex_by', eval = print_strata()} +``` + +```{r ref.label='cu_matur_title', eval = print_strata()} +``` + +```{r ref.label='cu_matur_by', eval = print_strata()} +``` + +```{r ref.label='pa_sex_title', eval = print_strata()} +``` + +```{r ref.label='pa_sex_by', eval = print_strata()} +``` + +```{r ref.label='pa_matur_title', eval = print_strata()} +``` + +```{r ref.label='pa_matur_by', eval = print_strata()} +``` + + + + +```{r ref.label='strata_plus_one'} +``` + +```{r ref.label='strata_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_title', eval = print_strata()} +``` + +```{r ref.label='ar_sex_by', eval = print_strata()} +``` + +```{r ref.label='ar_matur_title', eval = print_strata()} +``` + +```{r ref.label='ar_matur_by', eval = print_strata()} +``` + +```{r ref.label='cu_sex_title', eval = print_strata()} +``` + +```{r ref.label='cu_sex_by', eval = print_strata()} +``` + +```{r ref.label='cu_matur_title', eval = print_strata()} +``` + +```{r ref.label='cu_matur_by', eval = print_strata()} +``` + +```{r ref.label='pa_sex_title', eval = print_strata()} +``` + +```{r ref.label='pa_sex_by', eval = print_strata()} +``` + +```{r ref.label='pa_matur_title', eval = print_strata()} +``` + +```{r ref.label='pa_matur_by', eval = print_strata()} +``` + + + + +## Discussion + + + + + + +## Conclusion + + + + + + + + + +# Other ToRs + + + + + + + + +# References + +Gary D. Kader & Mike Perry (2007) Variability for Categorical Variables, +Journal of Statistics Education, 15:2, , DOI: 10.1080/10691898.2007.11889465 + +ICES (2014) Report of the Workshop on Statistical Analysis of Biological Calibration Studies (WKSABCAL). ICES CM 2014/ACOM: 35 + + + + +# Annex 1. Agenda + + +# Annex 2. List of participants + +```{r participants} + +# COMPLETE PARTICIPANTS LIST ################################################## + +# Table caption +set.caption('**Table X:** Participants list.') + +# Output table +pander(style_table0(stager_data), missing = "-") + +``` + + + +# Annex 3. Additional results + +***Data Summary Statistics*** + +```{r tab1_sum_stat_sex} + +# DATA SUMMARY STATISTICS for the sex category annotations- All stagers ####################################### + +# Table caption +set.caption('**Table X:** Summary of statistics; Number of samples, CU and PA (%) for the sex category annotations') + +# Output table +pander(style_table0(summary_statistics_sex_all), missing="") + +``` + +```{r tab1_sum_stat_maturity} + +# DATA SUMMARY STATISTICS for the maturity stage annotations- All stagers ####################################### + +# Table caption +set.caption('**Table X:** Summary of statistics; Number of samples, CU and PA (%) for the maturity stage annotations') + +# Output table +pander(style_table0(summary_statistics_maturity_all), missing="") + +``` + + +## Results all stagers + +**All samples included** + +```{r set_annex_all} +group <- "all" +``` + +***Data overview*** + +```{r data_overview_maturity} +# DATA OVERVIEW - All stagers ################################################# + +# Table caption +set.caption('**Table X:** Data overview including modal + maturity and statistics per sample.') + +# Output table +pander(style_table0(get(vname("data_overview_tab_maturity"))), missing="-", style = "simple") +``` + + +```{r data_overview_sex} +# DATA OVERVIEW - All stagers ################################################# + +# Table caption +set.caption('**Table X:** Data overview including modal + sex and statistics per sample.') + +# Output table +pander(style_table0(get(vname("data_overview_tab_sex"))), missing="-", style = "simple") +``` + + + +***Individual frequency plots by sex categorization*** + +```{r freq_bias_plots_sex, fig.width = 5, fig.height = 4, fig.cap = cap_in} +# INDIVIDUAL Frequency PLOTS - All stagers ######################################### + +# Figure caption +cap_in <- '**Figure X:** Frequency plots per stager. The modal sex category is plotted against the different sex categories. The frequency of categies corresponging to each modal sex category is represented together with the 1:1 equilibrium line (solid line).' + +# Output figures +plot_bias_sex(get(vname("ad_long"))) +``` + + +***Individual frequency plots by maturity stage*** + +```{r freq_bias_plots_maturity, fig.width = 5, fig.height = 4, fig.cap = cap_in} +# INDIVIDUAL Frequency PLOTS - All stagers ######################################### + +# Figure caption +cap_in <- '**Figure X:** Frequency plots per stager. The modal maturity stage is plotted against the different maturity stages. The frequency of staging’s corresponging to each modal maturity/maturity combinations is represented together with the 1:1 equilibrium line (solid line).' + +# Output figures +plot_bias_matur(get(vname("ad_long"))) +``` + + +***Plot with the CU AND PA statistichs for maturity staging’s*** + +```{r stat_plot_maturity, fig.width = 7, fig.height = 4, fig.cap = cap_in} +# PLOT STATISTICS, CU AND PA - All stagers ############################# + +# Figure caption +cap_in <- '**Figure X:** CU and PA are plotted against modal maturity' + +# The coefficient of variation (CV%), percent agreement and +# the standard deviation (STDEV) are plotted against modal maturity + +# Output figure +plot_stat_matur(get(vname("ad_long"))) +``` + + +***Plot with the CU AND PA statistichs for sex categorizations*** + +```{r stat_plot_sex, fig.width = 7, fig.height = 4, fig.cap = cap_in} +# PLOT STATISTICS, CU AND PA - All stagers ############################# + +# Figure caption +cap_in <- '**Figure X:** CU and PA are plotted against modal sex' + +# The coefficient of variation (CV%), percent agreement and +# the standard deviation (STDEV) are plotted against modal sex + +# Output figure +plot_stat_sex(get(vname("ad_long"))) +``` + + + +## Results Advanced stagers + +```{r set_annex_adv} +group <- "adv" +``` + +**All samples included** + +***Data overview*** + +```{r ref.label='data_overview_sex'} +``` + +```{r ref.label='data_overview_maturity'} +``` + +***Individual frequency plots by sex categorization*** + +```{r ref.label='freq_bias_plots_sex', fig.width = 5, fig.height = 4, fig.cap = cap_in} +``` + +***Individual frequency plots by maturity stage*** + +```{r ref.label='freq_bias_plots_maturity', fig.width = 5, fig.height = 4, fig.cap = cap_in} +``` + +***Plot with the CU AND PA statistichs for maturity staging’s*** + +```{r ref.label='stat_plot_maturity'} +``` + +***Plot with the CU AND PA statistichs for sex categorizations*** + +```{r ref.label='stat_plot_sex_adv'} +``` + + +# Annex 4. Table of the cases with multiple modes + +## All stagers + +***Multiple mode in the maturity stage*** + +```{r multiple mode cases maturity all stagers} +datos=ad_long_all[ad_long_all$NModes_trad_Maturity>1, c(1, 2, 4, 5, 8, 9, 10, 11, 12, 13, 16:26, 35, 36)] + +# Output table +pander(style_table0(datos), missing="-", style = "simple") + + +``` + +***Multiple mode in the sex category*** + +```{r multiple mode cases sex all stagers} +datos=ad_long_all[ad_long_all$NModes_trad_Sex>1, c(1, 2, 4, 5, 8, 9, 10, 11, 12, 13, 16:26, 56, 57)] + +# Output table +pander(style_table0(datos), missing="-", style = "simple") + + +``` + +## Advanced stagers + +***Multiple mode in the maturity stage*** + +```{r multiple mode cases maturity advanced stagers} +datos=ad_long_adv[ad_long_adv$NModes_trad_Maturity>1, c(1, 2, 4, 5, 8, 9, 10, 11, 12, 13, 16:26, 35, 36)] + +# Output table +pander(style_table0(datos), missing="-", style = "simple") + + +``` + +***Multiple mode in the sex category*** + +```{r multiple mode cases sex advanced stagers} +datos=ad_long_adv[ad_long_adv$NModes_trad_Sex>1, c(1, 2, 4, 5, 8, 9, 10, 11, 12, 13, 16:26, 56, 57)] + +# Output table +pander(style_table0(datos), missing="-", style = "simple") + + +``` + + +# Annex 5. ToRs for next meeting + + + +# Annex 6. Recommendations + + + +# Annex 7. Report specific annexes diff --git a/report_summary.Rmd b/report_summary.Rmd index 6f3746e..9c32d71 100644 --- a/report_summary.Rmd +++ b/report_summary.Rmd @@ -4,7 +4,7 @@ output: fig_caption: true fig_height: 10 fig_width: 10 - reference_docx: bootstrap/data/summaryTemplate.docx + reference_docx: bootstrap/initial/data/summaryTemplate.docx toc: false keep_md: false params: @@ -42,7 +42,7 @@ panderOptions('keep.trailing.zeros', TRUE) panderOptions('table.alignment.default', "center") panderOptions('knitr.auto.asis', FALSE) -group <- "ex" +group <- "adv" # from extrafont # font_import(pattern=c("Cali")) @@ -110,7 +110,7 @@ If there were available histological samples for all the fish individuals in the ```{r summary_multiple_modes_sex} -histN=ad_long_all %>% subset(TypeAnnotation=="eventOrganizer") %>% select(FishID) %>% unique() %>% dim() +histN=ad %>% subset(DoesSampleHaveHistologyImage=="Yes") %>% select(FishID) %>% unique() %>% dim() allN=ad_long_all %>% select(FishID) %>% unique() %>% dim() perc_not_Hist=100-round(100*histN[1]/allN[1], 0) @@ -140,7 +140,7 @@ If there were available histological samples for all the fish individuals in the if(!isTRUE(config$maturity_hist)){ -histN=ad_long_all %>% subset(TypeAnnotation=="eventOrganizer") %>% select(FishID) %>% unique() %>% dim() +histN=ad %>% subset(DoesSampleHaveHistologyImage=="Yes") %>% select(FishID) %>% unique() %>% dim() allN=ad_long_all %>% select(FishID) %>% unique() %>% dim() perc_not_Hist=100-round(100*histN[1]/allN[1], 0) @@ -166,7 +166,7 @@ pander(style_table0(summary_multiple_modes_maturity_all), missing="") **Sex categorization table** -```{r num_read_sex_ex} +```{r sex_composition_tab_adv} # Sex categorization TABLE - Advanced stagers ################################################# @@ -175,14 +175,14 @@ set.caption('**Table X:** Sex categorization table: presents the number of categorizations made per expert for each modal sex category.') # Output table -pander(style_table0(sex_composition_tab_ex), missing = "-", style = "simple") +pander(style_table0(sex_composition_tab_adv), missing = "-", style = "simple") ``` **Maturity staging table** -```{r num_read_maturity_ex} +```{r num_read_maturity_adv} # maturity stagings TABLE - Advanced stagers ################################################# @@ -191,14 +191,14 @@ set.caption('**Table X:** Maturity staging table: presents the number of stagings made per expert stager for each modal maturity stage.') # Output table -pander(style_table0(maturity_composition_tab_ex), missing = "-", style = "simple") +pander(style_table0(maturity_composition_tab_adv), missing = "-", style = "simple") ``` **Coefficient of Unalikeability (CU) table by modal sex category** -```{r cv_sex_ex} +```{r cv_sex_adv} # CU TABLE by modal sex category - Advanced stagers ################################################# @@ -209,14 +209,14 @@ set.caption('**Table X:** Coefficient of unlikeability (CU) table by modal sex c per stager.') # Output table -pander(style_table1(cu_tab_sex_ex), missing = "-") +pander(style_table1(cu_tab_sex_adv), missing = "-") ``` **Coefficient of Unalikeability (CU) table by modal maturity stage** -```{r cv_maturity_ex} +```{r cv_maturity_adv} # CU TABLE by modal maturity stage - Advanced stagers ################################################# @@ -227,14 +227,14 @@ set.caption('**Table X:** Coefficient of unlikeability (CU) table by modal matur per stager.') # Output table -pander(style_table1(cu_tab_maturity_ex), missing = "-") +pander(style_table1(cu_tab_maturity_adv), missing = "-") ``` **PA table by modal sex category** -```{r percentage_agreement_sex_ex} +```{r percentage_agreement_sex_adv} # PERCENTAGE AGREEMENT TABLE by modal sex category - Advanced stagers ############################### @@ -244,14 +244,14 @@ set.caption('**Table X:** Percentage agreement (PA) table by modal sex category: combined and a weighted mean of the PA per stager.') # Output table -pander(style_table1(pa_tab_sex_ex), missing = "-") +pander(style_table1(pa_tab_sex_adv), missing = "-") ``` **PA table by modal maturity stage** -```{r percentage_agreement_maturity_ex} +```{r percentage_agreement_maturity_adv} # PERCENTAGE AGREEMENT TABLE by modal maturity stage - Advanced stagers ############################### @@ -261,7 +261,7 @@ set.caption('**Table X:** Percentage agreement (PA) table: presents the PA per combined and a weighted mean of the PA per stager.') # Output table -pander(style_table1(pa_tab_maturity_ex), missing = "-") +pander(style_table1(pa_tab_maturity_adv), missing = "-") ``` @@ -269,7 +269,7 @@ pander(style_table1(pa_tab_maturity_ex), missing = "-") **General frequency bias table by modal sex category ** -```{r rb_sex_ex,results='asis'} +```{r rb_sex_adv,results='asis'} # Frequency bias table by modal sex category - Advanced stagers ###################################### @@ -278,14 +278,14 @@ set.caption('**Table X:** Frequency bias table by modal sex category: presents t modal maturity for all the advanced stager together.') # Output table -pander(style_table0(general_bias_freq_tab_sex_ex), missing = "-") +pander(style_table0(general_bias_freq_tab_sex_adv), missing = "-") ``` **General frequency bias table by modal maturity stage** -```{r rb_maturity_ex,results='asis'} +```{r rb_maturity_adv,results='asis'} # Frequency bias table by modal maturity stage - Advanced stagers ###################################### @@ -294,7 +294,7 @@ set.caption('**Table X:** Frequency bias table by modal maturity stage: presents modal maturity stage for all the advanced stager together.') # Output table -pander(style_table0(general_bias_freq_tab_maturity_ex), missing = "-") +pander(style_table0(general_bias_freq_tab_maturity_adv), missing = "-") ``` @@ -310,7 +310,7 @@ cap_in <- '**Figure X:** Sex categorization bias plot by modal sex category for modal sex category and sex category for all the advanced stager together. ' # Output figure -plot_general_freq_sex(ad_long_ex, strata=NULL) +plot_general_freq_sex(ad_long_adv, strata=NULL) ``` @@ -327,7 +327,7 @@ cap_in <- '**Figure X:** Maturity staging bias plot by modal maturity stage for modal maturity and maturity stage for all the advanced stager together. ' # Output figure -plot_general_freq_matur(ad_long_ex, strata=NULL) +plot_general_freq_matur(ad_long_adv, strata=NULL) ``` diff --git a/utilities_model.R b/utilities_model.R index 0f36990..c81b29f 100644 --- a/utilities_model.R +++ b/utilities_model.R @@ -1,4 +1,3 @@ - format_table_modal_maturity <- function(tab, fmt = "%i", extra_rows = "Total", matur_unique = modal_matur_unique) { # tab[] <- lapply(tab, function(x) ifelse(is.na(x), "-", sprintf(fmt, x))) cbind(`Modal maturity` = c(matur_unique, extra_rows), tab) @@ -22,7 +21,6 @@ format_table_sex_stage <- function(tab, fmt = "%i", extra_rows = "Total", sex_un sample_data_overview_table <- function(dat, strata) { - #debug #dat<-ad_long_all #strata<-config$strata @@ -236,7 +234,7 @@ pa_table <- function(ad_long, varmod, by = "reader") { reader_freq_table <- function(ad_long, varmod, by=NULL) { if(is.null(by)) { - reader_freq_tab <- + reader_freq_tab <- ad_long %>% ddply(.(eval(parse(text=paste0("modal_", tolower(varmod)))), get(all_of(varmod)), reader), summarise, count=length(get(all_of(varmod)))) colnames(reader_freq_tab) <- c(paste0("modal_", tolower(varmod)), tolower(varmod), "reader","count") @@ -252,13 +250,13 @@ reader_freq_table <- function(ad_long, varmod, by=NULL) { reader_freq_tab <- reader_freq_tab %>% select(-count.x, -count.y) - reader_freq_tab <- + reader_freq_tab <- reader_freq_tab %>% spread(key=reader, value=freqrelat) #temp <- # expand.grid( - # unique(ad_long[, paste0("modal_", tolower(varmod))]), + # unique(ad_long[, paste0("modal_", tolower(varmod))]), # unique(ad_long[, varmod]) # ) temp <- data.frame(expand.grid(unique(eval(parse(text=paste0("ad_long$modal_", tolower(all_of(varmod)))))), unique(eval(parse(text=paste0("ad_long$", all_of(varmod))))))) @@ -271,10 +269,10 @@ reader_freq_table <- function(ad_long, varmod, by=NULL) { reader_freq_tab <- rbind(reader_freq_tab, c("all", "all", as.numeric(colMeans(reader_freq_tab[,3:dim(reader_freq_tab)[2]], na.rm=TRUE)))) reader_freq_tab[is.na(reader_freq_tab)] <- 0 + reader_freq_tab$Total <- round(rowMeans(as.data.frame(lapply(reader_freq_tab[,3:dim(reader_freq_tab)[2]],as.numeric)), na.rm=TRUE), digits=3) + reader_freq_tab[is.na(reader_freq_tab)] <- "-" - reader_freq_tab2 <- slice(reader_freq_tab, 1:(n() - 1)) - } else { reader_freq_tab <- @@ -306,15 +304,14 @@ reader_freq_table <- function(ad_long, varmod, by=NULL) { reader_freq_tab=rbind(reader_freq_tab, c("all", "all", as.numeric(colMeans(reader_freq_tab[,4:dim(reader_freq_tab)[2]], na.rm=TRUE)))) - reader_freq_tab[is.na(reader_freq_tab)]<-0 + reader_freq_tab[is.na(reader_freq_tab)]=0 reader_freq_tab$Total=round(rowMeans(as.data.frame(lapply(reader_freq_tab[,4:dim(reader_freq_tab)[2]],as.numeric)), na.rm=TRUE), digits=3) - - reader_freq_tab[is.na(reader_freq_tab)]<-"-" - reader_freq_tab2 <- slice(reader_freq_tab, 1:(n() - 1)) + + reader_freq_tab[is.na(reader_freq_tab)]="-" } - return(reader_freq_tab2) + return(reader_freq_tab) } @@ -338,8 +335,7 @@ general_freq_table <- function(ad_long, varmod, by=NULL) { general_freq_tab$modes=as.character(general_freq_tab$modes) general_freq_tab=general_freq_tab[with(general_freq_tab, order(modes)),] - general_freq_tab[is.na(general_freq_tab)]<-"-" - general_freq_tab2 <- slice(general_freq_tab, 1:(n() - 1)) + general_freq_tab[is.na(general_freq_tab)]="-" colnames(general_freq_tab)=c(paste0("modal_", tolower(varmod), "/", tolower(varmod), "_stages"), colnames(general_freq_tab)[-1]) } else { @@ -350,20 +346,21 @@ general_freq_table <- function(ad_long, varmod, by=NULL) { general_freq_tab <- general_freq_tab %>% spread(key=get(all_of(varmod)), value=frequency) - general_freq_tab<-general_freq_tab[with(general_freq_tab, order(get(all_of(by)))),] - general_freq_tab[is.na(general_freq_tab)]<-"-" - general_freq_tab2 <- slice(general_freq_tab, 1:(n() - 1)) + general_freq_tab=general_freq_tab[with(general_freq_tab, order(get(all_of(by)))),] + general_freq_tab[is.na(general_freq_tab)]="-" } - return(general_freq_tab2) + return(general_freq_tab) } data_overview_table <- function(dat, varmod, report_token) { - if (any(dat$TypeAnnotation == "eventOrganizer")) { + #if (any(dat$TypeAnnotation == "eventOrganizer")) { + if (any(dat$DoesSampleHaveHistologyImage == "Yes")) { hist <- - dat[dat$TypeAnnotation == "eventOrganizer", ] %>% + #dat[dat$TypeAnnotation == "eventOrganizer", ] %>% + dat[dat$DoesSampleHaveHistologyImage == "Yes", ] %>% ddply(.(SampleID), summarise, Histology = "yes") %>% select(., c("SampleID", "Histology")) dat <- merge(dat, hist, by.x = "SampleID", by.y = "SampleID", all.x = TRUE) @@ -382,7 +379,7 @@ data_overview_table <- function(dat, varmod, report_token) { readings <- ad_wide %>% select(matches("R[0-9][0-9]*")) - + complete <- complete.cases(readings) ad_wide[c("Mode", "PA %", "CU %")] <- NA @@ -436,6 +433,7 @@ maturity_composition <- function(dat, by = "reader") { as.data.frame } + sex_composition <- function(dat, by = "reader") { # Number of gonads staged per reader and maturity stage dat %>% diff --git a/utilities_report.R b/utilities_report.R index 5e98301..daccb53 100644 --- a/utilities_report.R +++ b/utilities_report.R @@ -1,399 +1,401 @@ - - -# Style output tables ######################################################### - -# These four functions are used to change the style of the output tables. -# Depending on the form of the table different styels are chiosen. - -# Style 0 -style_table0 <- function(tab) { - - # Capitalize first letter of column and make header boldface - names(tab) <- pandoc.strong.return(names(tab)) - - return(tab) -} - -# Style 1 -style_table1 <- function(tab) { - - # Capitalize first letter of column, make header, last column and second - # last row in boldface and make last row italic - names(tab) <- pandoc.strong.return(names(tab)) - emphasize.strong.cols(ncol(tab)) - emphasize.strong.rows(nrow(tab)) - - return(tab) -} - -# Style 2 -style_table2 <- function(tab) { - - # Capitalize first letter of column, make header, last column and - # last row in boldface - names(tab) <- pandoc.strong.return(names(tab)) - emphasize.strong.cols(ncol(tab)) - emphasize.strong.rows(nrow(tab)) - - return(tab) -} - -# Style 3 -style_table3 <- function(tab) { - - # Capitalize first letter of column, make header and first column boldface - names(tab) <- pandoc.strong.return(names(tab)) - emphasize.strong.cols(1) - - return(tab) -} - - - - -# Here the maturity bias are plotted for all readers combined. - -plot_general_freq_matur <- function(ad_long, strata=NULL){ - - # There are conflicts between plyr and ddply, and, it seems that it is necesary unload plyr and load it again for the function ddply and other working properly - unloadNamespace("plyr") - library("plyr") - - freq_tab <- - ad_long %>% - ddply(.(modal_maturity, Maturity), summarize, count=length(Maturity)) - - total_tab <- - ad_long %>% - ddply(.(modal_maturity), summarize, count=length(Maturity)) - - freq_tab=merge(freq_tab, total_tab, by.x="modal_maturity", by.y="modal_maturity") - freq_tab$count=round(freq_tab$count.x/freq_tab$count.y, digits=3) - freq_tab=freq_tab[,-c(3,4)] - - diffmodal=setdiff(sort(unique(ad_long$Maturity)), sort(unique(freq_tab$modal_maturity))) - mat=sort(unique(ad_long$Maturity)) - completedat=as.data.frame(expand.grid(diffmodal,mat)) - completedat$count=rep(0,dim(completedat)[1]) - colnames(completedat)=c("Maturity", "modal_maturity", "count") - - freq_tab=rbind(freq_tab, completedat) - freq_tab=as.data.frame(ddply(freq_tab,.(modal_maturity, Maturity), summarize, count=if(sum(count)==0){NA} else {sum(count, na.rm=TRUE)})) - - ggplot(freq_tab, aes(x = modal_maturity, y = Maturity, size = count, label = count)) + - geom_point(shape = 21, fill = "white") + - scale_size(range = c(2, 15), guide = "legend") + - geom_text(aes(label=as.character(count),hjust=0.5,vjust=-1), size=3) + - geom_abline(mapping = NULL, data = freq_tab, slope=1, 0, na.rm = FALSE, show.legend = NA) + - ggtitle("Relative frequency plot") + - xlab("Modal maturity") + ylab("Maturity stage") + - theme( - plot.title = element_text(color="black", size=20, face="bold", hjust=0.5), - axis.title.x = element_text(color="black", size=15, face="bold"), - axis.title.y = element_text(color="black", size=15, face="bold"), - axis.text.x = element_text(size = 15), - axis.text.y = element_text(size = 15) - ) + - theme(legend.position = "none") + - theme(panel.background = element_rect(fill = "#BFD5E3", colour = "#6D9EC1",size = 2, linetype = "solid")) - - } - - -# Here the sex categorization bias are plotted for all readers combined. - -plot_general_freq_sex <- function(ad_long, strata=NULL){ - - # There are conflicts between plyr and ddply, and, it seems that it is necesary unload plyr and load it again for the function ddply and other working properly - unloadNamespace("plyr") - library("plyr") - - freq_tab <- - ad_long %>% - ddply(.(modal_sex, Sex), summarize, count=length(Sex)) - - total_tab <- - ad_long %>% - ddply(.(modal_sex), summarize, count=length(Sex)) - - freq_tab=merge(freq_tab, total_tab, by.x="modal_sex", by.y="modal_sex") - freq_tab$count=round(freq_tab$count.x/freq_tab$count.y, digits=3) - freq_tab=freq_tab[,-c(3,4)] - - diffmodal=setdiff(sort(unique(ad_long$Sex)), sort(unique(freq_tab$modal_sex))) - sex=sort(unique(ad_long$Sex)) - completedat=as.data.frame(expand.grid(diffmodal,sex)) - completedat$count=rep(0,dim(completedat)[1]) - colnames(completedat)=c("modal_sex", "Sex", "count") - - freq_tab=rbind(freq_tab, completedat) - freq_tab=as.data.frame(ddply(freq_tab,.(modal_sex, Sex), summarize, count=if(sum(count)==0){NA} else {sum(count, na.rm=TRUE)})) - - ggplot(freq_tab, aes(x = modal_sex, y = Sex, size = count, label = count)) + - geom_point(shape = 21, fill = "white") + - scale_size(range = c(2, 25), guide = "legend") + - geom_text(aes(label=as.character(count),hjust=0.5,vjust=-1), size=3) + - geom_abline(mapping = NULL, data = freq_tab, slope=1, 0, na.rm = FALSE, show.legend = NA) + - ggtitle("Relative frequency plot") + - xlab("Modal sex") + ylab("Sex category") + - theme( - plot.title = element_text(color="black", size=20, face="bold", hjust=0.5), - axis.title.x = element_text(color="black", size=15, face="bold"), - axis.title.y = element_text(color="black", size=15, face="bold"), - axis.text.x = element_text(size = 15), - axis.text.y = element_text(size = 15) - ) + - theme(legend.position = "none") + - theme(panel.background = element_rect(fill = "#BFD5E3", colour = "#6D9EC1",size = 2, linetype = "solid")) - - -} - - - -# Sex and maturity bias plots ############################################################# - - -# The age bias plots are made per reader. -# For each reader is the mean age and 2 times standard deviation -# (as error bars) plotted against the modal age. - -plot_bias_matur <- function(ad_long) { - - # There are conflicts between plyr and ddply, and, it seems that it is necesary unload plyr and load it again for the function ddply and other working properly - unloadNamespace("plyr") - library("plyr") - - ad_long$reader=as.character(ad_long$reader) - # return list of plots - lapply(sort(unique(ad_long$reader)), function(ireader) { - # Plot data and make settings/design - freq_tab <- - ad_long %>% - filter(reader == ireader) %>% - ddply(.(modal_maturity, Maturity, reader), summarize, count=length(Maturity)) - - diffmodal=setdiff(sort(unique(ad_long$Maturity)), sort(unique(freq_tab$modal_maturity))) - mat=sort(unique(ad_long$Maturity)) - completedat=as.data.frame(expand.grid(diffmodal, mat)) - completedat$reader=rep(unique(freq_tab$reader),dim(completedat)[1]) - completedat$count=rep(NA,dim(completedat)[1]) - colnames(completedat)=c( "modal_maturity", "Maturity", "reader", "count") - - freq_tab=rbind(freq_tab, completedat) - freq_tab=as.data.frame(ddply(freq_tab,.(modal_maturity, Maturity, reader), summarize, count=sum(count))) - - total_tab <- - ad_long %>% - filter(reader == ireader) %>% - ddply(.(modal_maturity), summarize, count=length(Maturity)) - - freq_tab=merge(freq_tab, total_tab, by.x="modal_maturity", by.y="modal_maturity", all.x=TRUE) - freq_tab$count=round(freq_tab$count.x/freq_tab$count.y, digits=3) - freq_tab=freq_tab[,-c(4,5)] - - ggplot(freq_tab, aes(x = modal_maturity, y = Maturity, size = count, label = count, group = reader)) + - geom_point(shape = 21, fill = "white") + - geom_text(aes(label=as.character(count),hjust=0.5,vjust=-1), size=3) + - scale_size(range = c(2, 15), guide = "legend") + - geom_abline(mapping = NULL, data = freq_tab, slope=1, 0, na.rm = FALSE, show.legend = NA) + - facet_wrap(~ reader, ncol = 2) + - theme(strip.text.x=element_text(size=25, hjust=0.5, vjust=0.5)) + - xlab("Modal maturity") + ylab("Maturity stage") + - theme( - plot.title = element_text(color="black", size=20, face="bold", hjust=0.5), - axis.title.x = element_text(color="black", size=15, face="bold"), - axis.title.y = element_text(color="black", size=15, face="bold"), - axis.text.x = element_text(size = 15), - axis.text.y = element_text(size = 15) - ) + - theme(legend.position = "none") + - theme(panel.background = element_rect(fill = "#BFD5E3", colour = "#6D9EC1",size = 2, linetype = "solid")) - - }) -} - - -plot_bias_sex <- function(ad_long) { - - # There are conflicts between plyr and ddply, and, it seems that it is necesary unload plyr and load it again for the function ddply and other working properly - unloadNamespace("plyr") - library("plyr") - - ad_long$reader=as.character(ad_long$reader) - # return list of plots - lapply(sort(unique(ad_long$reader)), function(ireader) { - # Plot data and make settings/design - freq_tab <- - ad_long %>% - filter(reader == ireader) %>% - ddply(.(modal_sex, Sex, reader), summarize, count=length(Sex)) - - diffmodal=setdiff(sort(unique(ad_long$Sex)), sort(unique(freq_tab$modal_sex))) - mat=sort(unique(ad_long$Sex)) - completedat=as.data.frame(expand.grid(diffmodal, mat)) - completedat$reader=rep(unique(freq_tab$reader),dim(completedat)[1]) - completedat$count=rep(NA,dim(completedat)[1]) - colnames(completedat)=c( "modal_sex", "Sex", "reader", "count") - - freq_tab=rbind(freq_tab, completedat) - freq_tab=as.data.frame(ddply(freq_tab,.(modal_sex, Sex, reader), summarize, count=sum(count))) - - total_tab <- - ad_long %>% - filter(reader == ireader) %>% - ddply(.(modal_sex), summarize, count=length(Sex)) - - freq_tab=merge(freq_tab, total_tab, by.x="modal_sex", by.y="modal_sex", all.x=TRUE) - freq_tab$count=round(freq_tab$count.x/freq_tab$count.y, digits=3) - freq_tab=freq_tab[,-c(4,5)] - - ggplot(freq_tab, aes(x = modal_sex, y = Sex, size = count, label = count, group = reader)) + - geom_point(shape = 21, fill = "white") + - geom_text(aes(label=as.character(count),hjust=0.5,vjust=-1), size=3) + - scale_size(range = c(2, 15), guide = "legend") + - geom_abline(mapping = NULL, data = freq_tab, slope=1, 0, na.rm = FALSE, show.legend = NA) + - facet_wrap(~ reader, ncol = 2) + - theme(strip.text.x=element_text(size=25, hjust=0.5, vjust=0.5)) + - xlab("Modal sex") + ylab("Sex stage") + - theme( - plot.title = element_text(color="black", size=20, face="bold", hjust=0.5), - axis.title.x = element_text(color="black", size=15, face="bold"), - axis.title.y = element_text(color="black", size=15, face="bold"), - axis.text.x = element_text(size = 15), - axis.text.y = element_text(size = 15) - ) + - theme(legend.position = "none") + - theme(panel.background = element_rect(fill = "#BFD5E3", colour = "#6D9EC1",size = 2, linetype = "solid")) - - }) -} - - - -# Plot std, ca and pa ######################################################### - -# Plot overall std, CA and PA per modal maturity stage in same plot -plot_stat_matur <- function(ad_long) { - - - # Again, beacuse of conflicts between plyr and dplyr it is necesary to unload the library plyr - unloadNamespace("plyr") - - # Combine the three data sets - dat_in <- - ad_long %>% - group_by(modal_maturity) %>% - summarise( - cu = cu_II(Maturity), - pa = mean(Maturity == modal_maturity, na.rm = TRUE) * 100) - - # Limit to use for axis - cu_lim <- ceiling(max(dat_in$cu, na.rm = T)) - - p <- - dat_in %>% - filter(!is.na(cu)) %>% - ggplot() + - theme_bw() + - # Coefficient of unalikeability - geom_line(aes(x = modal_maturity, y = cu, colour = "cu")) + - geom_point(aes(x = modal_maturity, y = cu, colour = "cu", - shape = "cu"), size = 3) + - # PA - geom_line(aes(x = modal_maturity, y = pa*cu_lim/100, colour = "pa")) + - geom_point(aes(x = modal_maturity, y = pa*cu_lim/100, colour = "pa", - shape = "pa"), size = 3) + - # Make left side y-axis - scale_y_continuous(name = expression("Coefficient of unalikability CU (years)"), - limits = c(0, cu_lim)) + - # Make right side y-axis - scale_y_continuous(name = expression("Coefficient of unalikability (years)"), - sec.axis = sec_axis(~ . * 100/cu_lim, - name = "Percentage of Agreement PA (%)"), - limits = c(0, cu_lim)) - print(p) - - # # Colors and labels - # p + - # theme(axis.text.y = element_text(color = "#80B1D3"), - # axis.text.y.right = element_text(color = "#FB8072") - # #text = element_text(family="Calibri") - # ) + - # scale_colour_manual(name = "Measure", - # values = c("#80B1D3", "#FB8072"), - # labels = c("CU", "PA")) + - # scale_shape_manual(name = "Measure", values = c(16, 8, 17), - # labels = c("CU", "PA")) + - # labs(x = "Modal maturity", colour = "") - # - # print(p) - # Now, I load again plyr and dplyr in this order: - library(plyr); library(dplyr) - -} - - - - - - -plot_stat_sex <- function(ad_long) { - - - # Again, beacuse of conflicts between plyr and dplyr it is necesary to unload the library plyr - unloadNamespace("plyr") - - # Combine the three data sets - dat_in <- - ad_long %>% - group_by(modal_sex) %>% - summarise( - cu = cu_II(Sex), - pa = mean(Sex == modal_sex, na.rm = TRUE) * 100) - - # Limit to use for axis - cu_lim <- ceiling(max(dat_in$cu, na.rm = T)) - - p <- - dat_in %>% - filter(!is.na(cu)) %>% - ggplot() + - theme_bw() + - # Coefficient of unalikeability - geom_line(aes(x = modal_sex, y = cu, colour = "cu")) + - geom_point(aes(x = modal_sex, y = cu, colour = "cu", - shape = "cu"), size = 3) + - # PA - geom_line(aes(x = modal_sex, y = pa*cu_lim/100, colour = "pa")) + - geom_point(aes(x = modal_sex, y = pa*cu_lim/100, colour = "pa", - shape = "pa"), size = 3) + - # Make left side y-axis - scale_y_continuous(name = expression("Coefficient of unalikability CU (years)"), - limits = c(0, cu_lim)) + - # Make right side y-axis - scale_y_continuous(name = expression("Coefficient of unalikability (years)"), - sec.axis = sec_axis(~ . * 100/cu_lim, - name = "Percentage of Agreement PA (%)"), - limits = c(0, cu_lim)) - print(p) - - # # Colors and labels - # p + - # theme(axis.text.y = element_text(color = "#80B1D3"), - # axis.text.y.right = element_text(color = "#FB8072") - # #text = element_text(family="Calibri") - # ) + - # scale_colour_manual(name = "Measure", - # values = c("#80B1D3", "#FB8072"), - # labels = c("CU", "PA")) + - # scale_shape_manual(name = "Measure", values = c(16, 8, 17), - # labels = c("CU", "PA")) + - # labs(x = "Modal sex", colour = "") - # - # print(p) - # Now, I load again plyr and dplyr in this order: - library(plyr); library(dplyr) - -} - + + +# Style output tables ######################################################### + +# These four functions are used to change the style of the output tables. +# Depending on the form of the table different styels are chiosen. + +# Style 0 +style_table0 <- function(tab) { + + # Capitalize first letter of column and make header boldface + names(tab) <- pandoc.strong.return(names(tab)) + + return(tab) +} + +# Style 1 +style_table1 <- function(tab) { + + # Capitalize first letter of column, make header, last column and second + # last row in boldface and make last row italic + names(tab) <- pandoc.strong.return(names(tab)) + emphasize.strong.cols(ncol(tab)) + emphasize.strong.rows(nrow(tab)) + + return(tab) +} + +# Style 2 +style_table2 <- function(tab) { + + # Capitalize first letter of column, make header, last column and + # last row in boldface + names(tab) <- pandoc.strong.return(names(tab)) + emphasize.strong.cols(ncol(tab)) + emphasize.strong.rows(nrow(tab)) + + return(tab) +} + +# Style 3 +style_table3 <- function(tab) { + + # Capitalize first letter of column, make header and first column boldface + names(tab) <- pandoc.strong.return(names(tab)) + emphasize.strong.cols(1) + + return(tab) +} + + + + +# Here the maturity bias are plotted for all readers combined. + +plot_general_freq_matur <- function(ad_long, strata=NULL){ + + # There are conflicts between plyr and ddply, and, it seems that it is necesary unload plyr and load it again for the function ddply and other working properly + unloadNamespace("plyr") + library("plyr") + + freq_tab <- + ad_long %>% + ddply(.(modal_maturity, Maturity), summarize, count=length(Maturity)) + + total_tab <- + ad_long %>% + ddply(.(modal_maturity), summarize, count=length(Maturity)) + + freq_tab=merge(freq_tab, total_tab, by.x="modal_maturity", by.y="modal_maturity") + freq_tab$count=round(freq_tab$count.x/freq_tab$count.y, digits=3) + freq_tab=freq_tab[,-c(3,4)] + + diffmodal=setdiff(sort(unique(ad_long$Maturity)), sort(unique(freq_tab$modal_maturity))) + mat=sort(unique(ad_long$Maturity)) + completedat=as.data.frame(expand.grid(diffmodal,mat)) + completedat$count=rep(0,dim(completedat)[1]) + colnames(completedat)=c("Maturity", "modal_maturity", "count") + + freq_tab=rbind(freq_tab, completedat) + freq_tab=as.data.frame(ddply(freq_tab,.(modal_maturity, Maturity), summarize, count=if(sum(count)==0){NA} else {sum(count, na.rm=TRUE)})) + + ggplot(freq_tab, aes(x = modal_maturity, y = Maturity, size = count, label = count)) + + geom_point(shape = 21, fill = "white") + + scale_size(range = c(2, 15), guide = "legend") + + geom_text(aes(label=as.character(count),hjust=0.5,vjust=-1), size=3) + + geom_abline(mapping = NULL, data = freq_tab, slope=1, 0, na.rm = FALSE, show.legend = NA) + + ggtitle("Relative frequency plot") + + xlab("Modal maturity") + ylab("Maturity stage") + + theme( + plot.title = element_text(color="black", size=20, face="bold", hjust=0.5), + axis.title.x = element_text(color="black", size=15, face="bold"), + axis.title.y = element_text(color="black", size=15, face="bold"), + axis.text.x = element_text(size = 15), + axis.text.y = element_text(size = 15) + ) + + theme(legend.position = "none") + + theme(panel.background = element_rect(fill = "#BFD5E3", colour = "#6D9EC1",size = 2, linetype = "solid")) + + } + + + + +# Here the sex categorization bias are plotted for all readers combined. + +plot_general_freq_sex <- function(ad_long, strata=NULL){ + + # There are conflicts between plyr and ddply, and, it seems that it is necesary unload plyr and load it again for the function ddply and other working properly + unloadNamespace("plyr") + library("plyr") + + freq_tab <- + ad_long %>% + ddply(.(modal_sex, Sex), summarize, count=length(Sex)) + + total_tab <- + ad_long %>% + ddply(.(modal_sex), summarize, count=length(Sex)) + + freq_tab=merge(freq_tab, total_tab, by.x="modal_sex", by.y="modal_sex") + freq_tab$count=round(freq_tab$count.x/freq_tab$count.y, digits=3) + freq_tab=freq_tab[,-c(3,4)] + + diffmodal=setdiff(sort(unique(ad_long$Sex)), sort(unique(freq_tab$modal_sex))) + sex=sort(unique(ad_long$Sex)) + completedat=as.data.frame(expand.grid(diffmodal,sex)) + completedat$count=rep(0,dim(completedat)[1]) + colnames(completedat)=c("modal_sex", "Sex", "count") + + freq_tab=rbind(freq_tab, completedat) + freq_tab=as.data.frame(ddply(freq_tab,.(modal_sex, Sex), summarize, count=if(sum(count)==0){NA} else {sum(count, na.rm=TRUE)})) + + ggplot(freq_tab, aes(x = modal_sex, y = Sex, size = count, label = count)) + + geom_point(shape = 21, fill = "white") + + scale_size(range = c(2, 25), guide = "legend") + + geom_text(aes(label=as.character(count),hjust=0.5,vjust=-1), size=3) + + geom_abline(mapping = NULL, data = freq_tab, slope=1, 0, na.rm = FALSE, show.legend = NA) + + ggtitle("Relative frequency plot") + + xlab("Modal sex") + ylab("Sex category") + + theme( + plot.title = element_text(color="black", size=20, face="bold", hjust=0.5), + axis.title.x = element_text(color="black", size=15, face="bold"), + axis.title.y = element_text(color="black", size=15, face="bold"), + axis.text.x = element_text(size = 15), + axis.text.y = element_text(size = 15) + ) + + theme(legend.position = "none") + + theme(panel.background = element_rect(fill = "#BFD5E3", colour = "#6D9EC1",size = 2, linetype = "solid")) + + +} + + + +# Sex and maturity bias plots ############################################################# + + +# The maturity bias plots are made per reader. +# For each reader is the mean age and 2 times standard deviation +# (as error bars) plotted against the modal age. + +plot_bias_matur <- function(ad_long) { + + # There are conflicts between plyr and ddply, and, it seems that it is necesary unload plyr and load it again for the function ddply and other working properly + unloadNamespace("plyr") + library("plyr") + + ad_long$reader=as.character(ad_long$reader) + # return list of plots + lapply(sort(unique(ad_long$reader)), function(ireader) { + # Plot data and make settings/design + freq_tab <- + ad_long %>% + filter(reader == ireader) %>% + ddply(.(modal_maturity, Maturity, reader), summarize, count=length(Maturity)) + + diffmodal=setdiff(sort(unique(ad_long$Maturity)), sort(unique(freq_tab$modal_maturity))) + mat=sort(unique(ad_long$Maturity)) + completedat=as.data.frame(expand.grid(diffmodal, mat)) + completedat$reader=rep(unique(freq_tab$reader),dim(completedat)[1]) + completedat$count=rep(NA,dim(completedat)[1]) + colnames(completedat)=c( "modal_maturity", "Maturity", "reader", "count") + + freq_tab=rbind(freq_tab, completedat) + freq_tab=as.data.frame(ddply(freq_tab,.(modal_maturity, Maturity, reader), summarize, count=sum(count))) + + total_tab <- + ad_long %>% + filter(reader == ireader) %>% + ddply(.(modal_maturity), summarize, count=length(Maturity)) + + freq_tab=merge(freq_tab, total_tab, by.x="modal_maturity", by.y="modal_maturity", all.x=TRUE) + freq_tab$count=round(freq_tab$count.x/freq_tab$count.y, digits=3) + freq_tab=freq_tab[,-c(4,5)] + + ggplot(freq_tab, aes(x = modal_maturity, y = Maturity, size = count, label = count, group = reader)) + + geom_point(shape = 21, fill = "white") + + geom_text(aes(label=as.character(count),hjust=0.5,vjust=-1), size=3) + + scale_size(range = c(2, 15), guide = "legend") + + geom_abline(mapping = NULL, data = freq_tab, slope=1, 0, na.rm = FALSE, show.legend = NA) + + facet_wrap(~ reader, ncol = 2) + + theme(strip.text.x=element_text(size=25, hjust=0.5, vjust=0.5)) + + xlab("Modal maturity") + ylab("Maturity stage") + + theme( + plot.title = element_text(color="black", size=20, face="bold", hjust=0.5), + axis.title.x = element_text(color="black", size=15, face="bold"), + axis.title.y = element_text(color="black", size=15, face="bold"), + axis.text.x = element_text(size = 15), + axis.text.y = element_text(size = 15) + ) + + theme(legend.position = "none") + + theme(panel.background = element_rect(fill = "#BFD5E3", colour = "#6D9EC1",size = 2, linetype = "solid")) + + }) +} + + +plot_bias_sex <- function(ad_long) { + + # There are conflicts between plyr and ddply, and, it seems that it is necesary unload plyr and load it again for the function ddply and other working properly + unloadNamespace("plyr") + library("plyr") + + ad_long$reader=as.character(ad_long$reader) + # return list of plots + lapply(sort(unique(ad_long$reader)), function(ireader) { + # Plot data and make settings/design + freq_tab <- + ad_long %>% + filter(reader == ireader) %>% + ddply(.(modal_sex, Sex, reader), summarize, count=length(Sex)) + + diffmodal=setdiff(sort(unique(ad_long$Sex)), sort(unique(freq_tab$modal_sex))) + mat=sort(unique(ad_long$Sex)) + completedat=as.data.frame(expand.grid(diffmodal, mat)) + completedat$reader=rep(unique(freq_tab$reader),dim(completedat)[1]) + completedat$count=rep(NA,dim(completedat)[1]) + colnames(completedat)=c( "modal_sex", "Sex", "reader", "count") + + freq_tab=rbind(freq_tab, completedat) + freq_tab=as.data.frame(ddply(freq_tab,.(modal_sex, Sex, reader), summarize, count=sum(count))) + + total_tab <- + ad_long %>% + filter(reader == ireader) %>% + ddply(.(modal_sex), summarize, count=length(Sex)) + + freq_tab=merge(freq_tab, total_tab, by.x="modal_sex", by.y="modal_sex", all.x=TRUE) + freq_tab$count=round(freq_tab$count.x/freq_tab$count.y, digits=3) + freq_tab=freq_tab[,-c(4,5)] + + ggplot(freq_tab, aes(x = modal_sex, y = Sex, size = count, label = count, group = reader)) + + geom_point(shape = 21, fill = "white") + + geom_text(aes(label=as.character(count),hjust=0.5,vjust=-1), size=3) + + scale_size(range = c(2, 15), guide = "legend") + + geom_abline(mapping = NULL, data = freq_tab, slope=1, 0, na.rm = FALSE, show.legend = NA) + + facet_wrap(~ reader, ncol = 2) + + theme(strip.text.x=element_text(size=25, hjust=0.5, vjust=0.5)) + + xlab("Modal sex") + ylab("Sex stage") + + theme( + plot.title = element_text(color="black", size=20, face="bold", hjust=0.5), + axis.title.x = element_text(color="black", size=15, face="bold"), + axis.title.y = element_text(color="black", size=15, face="bold"), + axis.text.x = element_text(size = 15), + axis.text.y = element_text(size = 15) + ) + + theme(legend.position = "none") + + theme(panel.background = element_rect(fill = "#BFD5E3", colour = "#6D9EC1",size = 2, linetype = "solid")) + + }) +} + + + +# Plot std, ca and pa ######################################################### + +# Plot overall std, CA and PA per modal maturity stage in same plot +plot_stat_matur <- function(ad_long) { + + + # Again, beacuse of conflicts between plyr and dplyr it is necesary to unload the library plyr + unloadNamespace("plyr") + + # Combine the three data sets + dat_in <- + ad_long %>% + group_by(modal_maturity) %>% + summarise( + cu = cu_II(Maturity), + pa = mean(Maturity == modal_maturity, na.rm = TRUE) * 100) + + # Limit to use for axis + cu_lim <- ceiling(max(dat_in$cu, na.rm = T)) + + p <- + dat_in %>% + filter(!is.na(cu)) %>% + ggplot() + + theme_bw() + + # Coefficient of unalikeability + geom_line(aes(x = modal_maturity, y = cu, colour = "cu")) + + geom_point(aes(x = modal_maturity, y = cu, colour = "cu", + shape = "cu"), size = 3) + + # PA + geom_line(aes(x = modal_maturity, y = pa*cu_lim/100, colour = "pa")) + + geom_point(aes(x = modal_maturity, y = pa*cu_lim/100, colour = "pa", + shape = "pa"), size = 3) + + # Make left side y-axis + scale_y_continuous(name = expression("Coefficient of unalikability CU (years)"), + limits = c(0, cu_lim)) + + # Make right side y-axis + scale_y_continuous(name = expression("Coefficient of unalikability (years)"), + sec.axis = sec_axis(~ . * 100/cu_lim, + name = "Percentage of Agreement PA (%)"), + limits = c(0, cu_lim)) + print(p) + + # # Colors and labels + # p + + # theme(axis.text.y = element_text(color = "#80B1D3"), + # axis.text.y.right = element_text(color = "#FB8072") + # #text = element_text(family="Calibri") + # ) + + # scale_colour_manual(name = "Measure", + # values = c("#80B1D3", "#FB8072"), + # labels = c("CU", "PA")) + + # scale_shape_manual(name = "Measure", values = c(16, 8, 17), + # labels = c("CU", "PA")) + + # labs(x = "Modal maturity", colour = "") + # + # print(p) + # Now, I load again plyr and dplyr in this order: + library(plyr); library(dplyr) + +} + + + + + + +plot_stat_sex <- function(ad_long) { + + + # Again, beacuse of conflicts between plyr and dplyr it is necesary to unload the library plyr + unloadNamespace("plyr") + + # Combine the three data sets + dat_in <- + ad_long %>% + group_by(modal_sex) %>% + summarise( + cu = cu_II(Sex), + pa = mean(Sex == modal_sex, na.rm = TRUE) * 100) + + # Limit to use for axis + cu_lim <- ceiling(max(dat_in$cu, na.rm = T)) + + p <- + dat_in %>% + filter(!is.na(cu)) %>% + ggplot() + + theme_bw() + + # Coefficient of unalikeability + geom_line(aes(x = modal_sex, y = cu, colour = "cu")) + + geom_point(aes(x = modal_sex, y = cu, colour = "cu", + shape = "cu"), size = 3) + + # PA + geom_line(aes(x = modal_sex, y = pa*cu_lim/100, colour = "pa")) + + geom_point(aes(x = modal_sex, y = pa*cu_lim/100, colour = "pa", + shape = "pa"), size = 3) + + # Make left side y-axis + scale_y_continuous(name = expression("Coefficient of unalikability CU (years)"), + limits = c(0, cu_lim)) + + # Make right side y-axis + scale_y_continuous(name = expression("Coefficient of unalikability (years)"), + sec.axis = sec_axis(~ . * 100/cu_lim, + name = "Percentage of Agreement PA (%)"), + limits = c(0, cu_lim)) + print(p) + + # # Colors and labels + # p + + # theme(axis.text.y = element_text(color = "#80B1D3"), + # axis.text.y.right = element_text(color = "#FB8072") + # #text = element_text(family="Calibri") + # ) + + # scale_colour_manual(name = "Measure", + # values = c("#80B1D3", "#FB8072"), + # labels = c("CU", "PA")) + + # scale_shape_manual(name = "Measure", values = c(16, 8, 17), + # labels = c("CU", "PA")) + + # labs(x = "Modal sex", colour = "") + # + # print(p) + # Now, I load again plyr and dplyr in this order: + library(plyr); library(dplyr) + +} +