# packages ####

#ochaRd
# 
# install.packages("devtools")
# install.packages("tidyverse")
# #install.packages("metafor")
# install.packages("patchwork")
# install.packages("R.rsp")
# 
# devtools::install_github("daniel1noble/orchaRd", force = TRUE, build_vignettes = TRUE)
# remotes::install_github("rvlenth/emmeans", dependencies = TRUE, build_opts = "") 
# 
# #emmeans
# remotes::install_github("rvlenth/emmeans", dependencies = TRUE, build_opts = "")
# # metafor
# install.packages("remotes")
# remotes::install_github("wviechtb/metafor")

# loading
pacman::p_load(tidyverse,
               metafor,
               pander,
               stringr,
               ape,
               kableExtra,
               patchwork,
               lme4,
               readxl,
               #emmeans,
               rotl,
               orchaRd,
               clubSandwich,
               MuMIn,
               png,
               grid,
               here,
               formatR,
               naniar,
               GoodmanKruskal,
               ggalluvial
)

# you may need to run this (we will plan to fix this)
library(groundhog)
groundhog.library("emmeans", '2022-04-01')

# need for metafor to understand MuMin 
eval(metafor:::.MuMIn)

# custom functions

# function for getting lnRR for proportional data (mortality)

lnrrp <- function(m1, m2, n1, n2) {
    # arcsine transforamtion
    asin_trans <- function(p) {
        asin(sqrt(p))
    }
    # SD for arcsine distribution (see Wiki -
    # https://en.wikipedia.org/wiki/Arcsine_distribution)
    var1 <- 1/8
    var2 <- 1/8
    # lnRR - with 2nd order correction
    lnrr <- log(asin_trans(m1)/asin_trans(m2)) + 0.5 * ((var1/(n1 * asin_trans(m1)^2)) -
        (var2/(n2 * asin_trans(m2)^2)))

    var <- var1/(n1 * asin_trans(m1)^2) + var1^2/(2 * n1^2 * asin_trans(m1)^4) +
        var2/(n2 * asin_trans(m2)^2) + var2^2/(2 * n2^2 * asin_trans(m2)^4)

    invisible(data.frame(yi = lnrr, vi = var))
}


# function to get to lnRR for longevity data (CV required) The method proposed
# in Nakagawa et al (2022) - missing SD method

lnrrm <- function(m1, m2, n1, n2, cv21, cv22) {
    # lnRR - with 2nd order correction
    lnrr <- log(m1/m2) + 0.5 * ((cv21/n1) - (cv22/n2))

    var <- (cv21/n1) + ((cv21^2)/(2 * n1^2)) + (cv22/n2) + ((cv22^2)/(2 * n2^2))

    invisible(data.frame(yi = lnrr, vi = var))
}

# for folded normal distribution see:
# https://en.wikipedia.org/wiki/Folded_normal_distribution

# folded mean
folded_mu <- function(mean, variance) {
    mu <- mean
    sigma <- sqrt(variance)
    fold_mu <- sigma * sqrt(2/pi) * exp((-mu^2)/(2 * sigma^2)) + mu * (1 - 2 * pnorm(-mu/sigma))
    fold_mu
}

# folded variance
folded_v <- function(mean, variance) {
    mu <- mean
    sigma <- sqrt(variance)
    fold_mu <- sigma * sqrt(2/pi) * exp((-mu^2)/(2 * sigma^2)) + mu * (1 - 2 * pnorm(-mu/sigma))
    fold_se <- sqrt(mu^2 + sigma^2 - fold_mu^2)
    # adding se to make bigger mean
    fold_v <- fold_se^2
    fold_v
}

# dat_full <- read_csv(here('data', 'dat_07072021.csv'), na = c('', 'NA'))
# dat_full <- read_csv(here('data', 'data_15022022.csv'), na = c('', 'NA'))
dat_full <- read_csv(here("data", "data_05052022.csv"), na = c("", "NA"))
# glimpse(dat_full)

# loading data ####

dat_full %>%
    filter(is.na(Treatment_lifespan_variable) == FALSE) %>%
    filter(Type_of_sterilization != "Vasectomy") %>%
    mutate_if(is.character, as.factor) -> dat


dim(dat)

dim(dat_full)

# separating two kinds

effect_type <- ifelse(str_detect(dat$Lifespan_parameter, "Me"), "longevity", "mortality")

# fix a typo in species name
dat$Species_Latin <- gsub("Macaca Fascicularis", "Macaca fascicularis", dat$Species_Latin)
dat$Species_Latin <- gsub("Equus caballus", "Equus ferus", dat$Species_Latin)

# creating the phylo column
dat$Phylogeny <- sub(" ", "_", dat$Species_Latin)
dat$Effect_type <- effect_type
dat$Effect_ID <- 1:nrow(dat)
# key variables
names(dat)

unique(dat$Species_Latin)

kable(dat, "html") %>%
    kable_styling("striped", position = "left") %>%
    scroll_box(width = "100%", height = "250px")

# let's get CVs

dat %>%
    group_by(Study) %>%
    summarise(cv2_cont = mean((Error_control_SD/Control_lifespan_variable)^2, na.rm = T),
        cv2_trt = mean((Error_experimental_SD/Treatment_lifespan_variable)^2, na.rm = T),
        cv2_opst = mean((Error_opposite_sex_SD/Opposite_sex_lifespan_variable)^2,
            na.rm = T), n_cont = mean(Sample_size_control, na.rm = T), n_trt = mean(Sample_size_sterilization,
            na.rm = T), n_opst = mean(Sample_size_opposite_sex, na.rm = T)) %>%
    ungroup() %>%
    summarise(cv2_cont = weighted.mean(cv2_cont, n_cont, na.rm = T), cv2_trt = weighted.mean(cv2_trt,
        n_trt, na.rm = T), cv2_opst = weighted.mean(cv2_opst, n_opst, na.rm = T)) ->
    cvs

# lnRR using CV
dat$yi <- ifelse(effect_type == "longevity", lnrrm(dat$Treatment_lifespan_variable,
    dat$Control_lifespan_variable, dat$Sample_size_sterilization, dat$Sample_size_control,
    cvs[["cv2_trt"]], cvs[["cv2_cont"]])[[1]], lnrrp(dat$Treatment_lifespan_variable,
    dat$Control_lifespan_variable, dat$Sample_size_sterilization, dat$Sample_size_control)[[1]])

dat$vi <- ifelse(effect_type == "longevity", lnrrm(dat$Treatment_lifespan_variable,
    dat$Control_lifespan_variable, dat$Sample_size_sterilization, dat$Sample_size_control,
    cvs[["cv2_trt"]], cvs[["cv2_cont"]])[[2]], lnrrp(dat$Treatment_lifespan_variable,
    dat$Control_lifespan_variable, dat$Sample_size_sterilization, dat$Sample_size_control)[[2]])

# getting effect size for the long format we create a longer data format

dat1 <- dat
dat2 <- dat

dat1$yi <- ifelse(effect_type == "longevity", lnrrm(dat$Control_lifespan_variable,
    dat$Opposite_sex_lifespan_variable, dat$Sample_size_control, dat$Sample_size_opposite_sex,
    cvs[["cv2_cont"]], cvs[["cv2_opst"]])[[1]], lnrrp(dat$Control_lifespan_variable,
    dat$Opposite_sex_lifespan_variable, dat$Sample_size_control, dat$Sample_size_opposite_sex)[[1]])

dat1$vi <- ifelse(effect_type == "longevity", lnrrm(dat$Control_lifespan_variable,
    dat$Opposite_sex_lifespan_variable, dat$Sample_size_control, dat$Sample_size_opposite_sex,
    cvs[["cv2_cont"]], cvs[["cv2_opst"]])[[2]], lnrrp(dat$Control_lifespan_variable,
    dat$Opposite_sex_lifespan_variable, dat$Sample_size_control, dat$Sample_size_opposite_sex)[[2]])

# here we create CM/F or CF/M
dat2$yi <- ifelse(effect_type == "longevity", lnrrm(dat$Treatment_lifespan_variable,
    dat$Opposite_sex_lifespan_variable, dat$Sample_size_sterilization, dat$Sample_size_opposite_sex,
    cvs[["cv2_trt"]], cvs[["cv2_opst"]])[[1]], lnrrp(dat$Treatment_lifespan_variable,
    dat$Opposite_sex_lifespan_variable, dat$Sample_size_sterilization, dat$Sample_size_opposite_sex)[[1]])

dat2$vi <- ifelse(effect_type == "longevity", lnrrm(dat$Treatment_lifespan_variable,
    dat$Opposite_sex_lifespan_variable, dat$Sample_size_sterilization, dat$Sample_size_opposite_sex,
    cvs[["cv2_trt"]], cvs[["cv2_opst"]])[[2]], lnrrp(dat$Treatment_lifespan_variable,
    dat$Opposite_sex_lifespan_variable, dat$Sample_size_sterilization, dat$Sample_size_opposite_sex)[[2]])




# putting two data frames
dat_long <- rbind(dat1, dat2)

# putt 2 new column

dat_long$Obs <- factor(1:dim(dat_long)[[1]])
dat_long$Comp_type <- as.factor(rep(c("both_normal", "one_castrated"), each = dim(dat_long)[[1]]/2))

dat_long$Comp_type_Sex <- paste(dat_long$Comp_type, dat_long$Sex, sep = "_")

dat_long$yi <- ifelse(dat_long$Comp_type_Sex == "one_castrated_Female" | dat_long$Comp_type_Sex ==
    "both_normal_Female", -1 * dat_long$yi, dat_long$yi)

dat_long %>%
    filter(!is.na(yi), !is.na(vi)) -> dat_long

dim(dat_long)