library(dplyr)
library(lubridate)
library(stringr)
library(writexl)
library(tidyr)
library(readr)
library(ggplot2)
library(broom)
library(ggalluvial)
library(car)
library(lmtest)
library(pROC)
library(ResourceSelection)
library(pscl)
library(splines)

##################################################
# ---- Andmete sisselugemine ja korrastamine ----
##################################################
cac_alg <- read.csv("/home/merit.matt/Documents/analyys/kaltsiumskoor.csv", stringsAsFactors = FALSE, na.string = c("", "NA"))
mootmised_alg <- read.table(file = "/home/merit.matt/Documents/analyys/measurements.tsv", sep = '\t', header = TRUE)
nmr_alg <- read.table(file = "/home/merit.matt/Documents/analyys/nmr.tsv", sep = '\t', header = TRUE)
ravimid_alg <- read.table(file = "/home/merit.matt/Documents/analyys/diagnosis-drugs.tsv", sep = '\t', header = TRUE)
kmi_suits_alg <- read.table(file = "/home/merit.matt/Documents/analyys/bmi-smoking.tsv", sep = '\t', header = TRUE)
prs <- read.table(file = "/home/merit.matt/Documents/analyys/prs.txt", sep = ' ', header = TRUE)

# --------------------------------------------
# ------- kaltsiumiskoor defineerimine -------
# --------------------------------------------
sum(is.na(cac_alg$value))
typeof(cac_alg$epicrisis_date_start)
cac_alg %>% summarise(sum(epicrisis_date_start == ""))

# ühe CACS ja mõõtmiskuupäeva saamine igale doonorile
f_cac_puhas <- function(cac_alg){
  df_puhas <- cac_alg %>%
    mutate(across(c(procedures_entry_date, date_from_text, epicrisis_date_start, epicrisis_date_end), ~ {
      as.numeric(parse_date_time(., orders = c("ymd HMS", "ymd HM", "ymd"), tz = "UTC", quiet = TRUE))
    })) %>%
    
    mutate(kuup_jrjk = coalesce(procedures_entry_date, date_from_text, epicrisis_date_start)) %>%
    
    group_by(skood) %>% 
    mutate(
      koik_samad = n_distinct(value, na.rm = TRUE) == 1,
      max_v = max(value, na.rm = TRUE)
    ) %>%
    
    mutate(
      final_date = ifelse(
        koik_samad, 
        max(epicrisis_date_start, na.rm = TRUE), 
        max(kuup_jrjk[value == max_v], na.rm = TRUE)
      )
    ) %>%
    
    mutate(sort = if_else(koik_samad, epicrisis_date_start, value)) %>%
    
    arrange(desc(sort), desc(kuup_jrjk), .by_group = TRUE) %>%
    slice(1) %>%
    ungroup() %>%
    
    mutate(across(c(final_date, procedures_entry_date, date_from_text, epicrisis_date_start, epicrisis_date_end), ~ as.POSIXct(., origin="1970-01-01", tz = "UTC"))) %>%
    select(skood, value, final_date)
  
  return(df_puhas)
}

cac_puhas <- f_cac_puhas(cac_alg)


##################################################
# ----------- Riskifaktorite määramine -----------
##################################################
# kuupäevade tüübi muutmine
f_muuda_kuupaevadeks <- function(mootmised_df, diagnoosid_ravimid_df, kmi_suits_df) {
  
  mootmised_df <- mootmised_df %>%
    mutate(across(contains("Date"), as.Date))
  
  diagnoosid_ravimid_df <- diagnoosid_ravimid_df %>%
    mutate(across(contains("Date"), as.Date))
  
  kmi_suits_df <- kmi_suits_df %>%
    mutate(across(contains("Date"), as.Date))
  
  return(list(
    mootmised = mootmised_df,
    diagnoosid = diagnoosid_ravimid_df,
    kmi_suits = kmi_suits_df
  ))
}

teisendatud_kuup <- f_muuda_kuupaevadeks(mootmised_alg, ravimid_alg, kmi_suits_alg)
mootmised_alg <- teisendatud_kuup$mootmised
diagnoosid_alg   <- teisendatud_kuup$diagnoosid
kmi_suits_alg <- teisendatud_kuup$kmi_suits

# --------------------------------------------
# ---- riskifaktorite mõõtmiste kontroll -----
# --------------------------------------------
cac_diagnoos <- cac_puhas %>%
  left_join(ravimid_alg, by = "skood")%>%
  left_join(mootmised_alg, by = "skood") %>%
  filter(Hypertension == 0 & Antihypertensives == 0 & loinc == "8480-6")
sum(is.na(cac_diagnoos$measurementResultNumeric))

cac_diagnoos <- cac_puhas %>%
   left_join(ravimid_alg, by = "skood")%>%
   left_join(mootmised_alg, by = "skood") %>%
   filter(Hypertension == 0 & Antihypertensives == 0 & loinc == "8462-4")
sum(is.na(cac_diagnoos$measurementResultNumeric))

proov <- cac_diagnoos %>% filter(loinc == "22748-8")
unique(proov$measurementUnit)
cac_diagnoos %>% filter(loinc == "22748-8" & measurementUnit == "ng/mL")
proov2 <- cac_diagnoos %>% filter(loinc == "22748-8" & measurementUnit == "mmol")
proov3 <- cac_diagnoos %>% filter(loinc == "22748-8" & measurementUnit == "mmol/L")
sum(is.na(proov2$measurementUnit))
sum(is.na(proov3$measurementUnit))

proova <- cac_diagnoos %>% filter(loinc == "14647-2")
unique(proova$measurementUnit)
proovb <- cac_diagnoos %>% filter(loinc == "4548-4")
unique(proovb$measurementUnit)
proovc <- cac_diagnoos %>% filter(loinc == "8480-6" | loinc == "8462-4")
unique(proovc$measurementUnit)

mootmised_alg %>% 
  group_by(skood) %>%
  filter(
    loinc == "4548-4",
    all(measurementUnit != "%"),
    any(measurementUnit == "mmol/mol")
  ) %>%
  ungroup()

# --------------------------------------------
# ---- mõõtmiste piiritlemine (outlierid) ----
# --------------------------------------------
# füsioloogiliste erindite eraldamine
# glükoosi puhul jäetud välja "14743-9", kuna pole paastuglükoos
f_mootmised_puhas <- function(df) {
  df_puhas <- df %>%
    filter(
      (loinc == "8480-6" & measurementResultNumeric >= 60 & measurementResultNumeric <= 280) |
        (loinc == "8462-4" & measurementResultNumeric >= 30 & measurementResultNumeric <= 160 )|
        ((loinc %in% c("14749-6", "14770-2", "14771-0", "72516-8") & measurementUnit == "mmol/L" | loinc == "14771-0" & is.na(measurementUnit)) & measurementResultNumeric >= 2 & measurementResultNumeric <= 30) |
        (loinc == "22748-8" & measurementResultNumeric >= 0.3 & measurementResultNumeric <= 12) |
        (loinc %in% c("33914-3", "50210-4", "62238-1", "A-4823") & (measurementUnit == "mL/min/1.73m2" | is.na(measurementUnit)) & measurementResultNumeric >= 5 & measurementResultNumeric <= 200) |
        (loinc == "14647-2" & measurementResultNumeric >= 1.5 & measurementResultNumeric <= 15) |
        (loinc == "4548-4" & (measurementUnit == "%" | is.na(measurementUnit))) |
        (loinc == "4548-4" & measurementUnit == "mmol/mol" & measurementResultNumeric >= 29 & measurementResultNumeric <= 42) 
    ) %>%
    mutate(
      measurementUnit = case_when(
        loinc == "4548-4" & is.na(measurementUnit) & measurementResultNumeric < 15 ~ "%",
        loinc == "4548-4" & is.na(measurementUnit) & measurementResultNumeric > 25 ~ "mmol/mol",
        TRUE ~ measurementUnit
      )
    ) %>%
    mutate(
      measurementResultNumeric = ifelse(
        (loinc == "4548-4" & measurementUnit == "mmol/mol"), 
        measurementResultNumeric * 0.09148 + 2.152, 
        measurementResultNumeric
      ),
      measurementUnit = ifelse((loinc == "4548-4" & measurementUnit == "mmol/mol"), "%", measurementUnit)
    )
  
  return(df_puhas)
}

mootmised_puhas <- f_mootmised_puhas(mootmised_alg)

# --------------------------------------------
# ------- kretiniinist eGFR arvutamine -------
# --------------------------------------------
tbl_info <- ravimid_alg %>% 
  transmute(skood, sex,
            ageAtSample = as.numeric(substr(firstSampleDate,1,4)) - birthYear
  )

tbl_nmr <- nmr_alg

# Metaboliidid on määratud erinevatel aegadel (tunnus batch): 2010-2012 üks metoodika ja 2022 teine.
# See ei mõjuta kõiki, aga mõnda metaboliiti. Antud NMR andmetest mõjutab kreateniini 
tbl_nmr %>%
  group_by(batch) %>%
  summarise(
    n = n(),
    mean = mean(Creatinine, na.rm = TRUE)
  )

# Seega tuleb määrata batchide 2010-2012 korral kreatiniin puuduvaks
tbl_nmr <- tbl_nmr %>%
  mutate(Creatinine = ifelse(batch == "b_2022", Creatinine, NA),
         creat_mgdl = Creatinine / 88.4 # vaja ühikud teisendada µmol/L -> mg/dl
  ) 

summary(tbl_nmr$Creatinine)
hist(tbl_nmr$Creatinine)

# CKD-EPI 2021 valem eGFR arvutamiseks (kreatiniin + vanus + sugu → eGFR)
# Viide: https://www.kidney.org/professionals/ckd-epi-creatinine-equation-2021
# Publikatsioon: https://www.nejm.org/doi/full/10.1056/NEJMoa2102953

# eGFR funktsioon
eGFR_CKD_EPI <- function(creat_mgdl, age, sex) {
  # sex: "male" or "female"
  
  kappa <- ifelse(sex == "female", 0.7, 0.9)
  alpha <- ifelse(sex == "female", -0.241, -0.302)
  female_factor <- ifelse(sex == "female", 1.018, 1)
  
  scr_kappa <- creat_mgdl / kappa
  
  eGFR <- 142 *
    pmin(scr_kappa, 1)^alpha *
    pmax(scr_kappa, 1)^(-1.200) *
    0.9938^age *
    female_factor
  
  return(eGFR)
}

tbl_egfr <- left_join(tbl_info, tbl_nmr, by = "skood") %>%
  mutate(
    eGFR = eGFR_CKD_EPI(creat_mgdl, ageAtSample, sex)
  )

# eGFR normväärtused: 15-140
summary(tbl_egfr$eGFR)
hist(tbl_egfr$eGFR)

# --------------------------------------------
# ------- riskifaktorite defineerimine -------
# --------------------------------------------
# peamiste riskifaktorite defineerimine
f_arvuta_riskid <- function(skoorid_df, diagnoosid_ravimid_df, mootmised_df, kmi_suits_df, egfr_df){
  riskid_df <- skoorid_df %>%
    left_join(diagnoosid_ravimid_df, by = "skood") %>%
    left_join(mootmised_df, by = "skood") %>%
    left_join(kmi_suits_df, by = "skood") %>%
    left_join(egfr_df, by = "skood") %>%
    
    mutate(
      synniaeg = ymd(paste0(birthYear, "-07-01")),
      vanus_indeksil = floor(time_length(interval(synniaeg, final_date), "years"))
    ) %>%  
    
    group_by(skood) %>%
    
    mutate(
      # ---- hüpertensioon ----
      hyp_diag = any(Hypertension == 1 & Hypertension_earlistDate <= final_date + years(1), na.rm = TRUE),
      hyp_ravim = any(Antihypertensives == 1 & Antihypertensives_earlistDate <= final_date + years(1), na.rm = TRUE),
      
      on_paaris_mootmine = length(intersect(
        measurementDate[which(loinc == "8480-6" & measurementDate >= final_date - years(5) & measurementDate <= final_date)],
        measurementDate[which(loinc == "8462-4" & measurementDate >= final_date - years(5) & measurementDate <= final_date)]
      )) > 0,
      
      hyp_labor_koos = length(intersect(
        measurementDate[which(loinc == "8480-6" & measurementResultNumeric >= 140 & measurementDate >= final_date - years(5) & measurementDate <= final_date)],
        measurementDate[which(loinc == "8462-4" & measurementResultNumeric >= 90 & measurementDate >= final_date - years(5) & measurementDate <= final_date)]
      )) > 1,
      
      hyp_on_mootmine = any(loinc %in% c("8480-6", "8462-4") & measurementDate >= final_date - years(5) & measurementDate <= final_date, na.rm = TRUE),
      
      hyp_labor_on_analyys = case_when(
        hyp_labor_koos ~ "Kõrge risk (laborist)", 
        on_paaris_mootmine ~ "Paarismõõtmised normis või üksikud",
        hyp_on_mootmine ~ "Vähe mõõtmisi (pole paaris)",
        TRUE ~ "Mõõtmised puuduvad"
      ),
      
      hypertensioon_risk = case_when(
        hyp_diag | hyp_ravim | hyp_labor_koos ~ 1,   
        hyp_on_mootmine ~ 0,                  
        TRUE ~ 2                        
      ),
      
      # ---- hüperlipideemia ----
      lip_diag = any(Hyperlipidemia == 1 & Hyperlipidemia_earlistDate <= final_date + years(1), na.rm = TRUE),
      lip_ravim = any(LLM == 1 & LLM_earlistDate <= final_date + years(1), na.rm = TRUE),
      
      lip_labor_ldl = any(loinc == "22748-8" & measurementResultNumeric >= 3.0 & measurementDate >= final_date - years(5) & measurementDate <= final_date, na.rm = TRUE) |
        any(LDL_C >= 3.0 & firstSampleDate >= final_date - years(5) & firstSampleDate <= final_date, na.rm = TRUE),
        
      lip_labor_koles = any(loinc == "14647-2" & measurementResultNumeric >= 5.0 &  measurementDate >= final_date - years(5) & measurementDate <= final_date, na.rm = TRUE) |
        any(Total_C >= 5.0 & firstSampleDate >= final_date - years(5) & firstSampleDate <= final_date, na.rm = TRUE),
    
      lip_on_mootmine = any(loinc %in% c("14647-2", "22748-8") & measurementDate >= final_date - years(5) & measurementDate <= final_date, na.rm = TRUE) |
        any(!is.na(LDL_C) & firstSampleDate >= final_date - years(5) & firstSampleDate <= final_date, na.rm = TRUE) |
        any(!is.na(Total_C) & firstSampleDate >= final_date - years(5) & firstSampleDate <= final_date, na.rm = TRUE),
      
      hyperlipideemia_risk = case_when(
        lip_diag | lip_ravim | lip_labor_ldl | lip_labor_koles ~ 1,   
        lip_on_mootmine ~ 0,                  
        TRUE ~ 2                        
      ),
      
      # ---- diabeet ----
      diab_diag = any(Diabetes == 1 & Diabetes_earlistDate <= final_date + years(1), na.rm = TRUE),
      diab_ravim = any(diabetes_drugs == 1 & diabetes_drugs_earlistDate <= final_date + years(1), na.rm = TRUE),

      diab_labor_HbA1c = any(loinc == "4548-4" & measurementUnit == "%" & measurementResultNumeric >= 6.5 & measurementDate >= final_date - years(5) & measurementDate <= final_date, na.rm = TRUE),
      
      diab_on_mootmine = any(loinc == "4548-4" & measurementDate >= final_date - years(5) & measurementDate <= final_date, na.rm = TRUE),
      
      diabeet_risk1 = case_when(
        diab_diag | diab_ravim | diab_labor_HbA1c ~ 1,   
        diab_on_mootmine ~ 0,                  
        TRUE ~ 2
      ),
      
      # ---- neerufunktsioon ----
      neerud_labor_egfr = any(eGFR < 60 & firstSampleDate >= final_date - years(5) & firstSampleDate <= final_date, na.rm = TRUE),
      neerud_labor_loinc = any(loinc %in% c("33914-3", "50210-4", "62238-1", "A-4823") & measurementResultNumeric < 60 & measurementDate >= final_date - years(5) & measurementDate <= final_date, na.rm = TRUE),
      
      neerud_on_mootmine = any(loinc %in% c("33914-3", "50210-4", "62238-1", "A-4823") & measurementDate >= final_date - years(5) & measurementDate <= final_date, na.rm = TRUE) |
        any(!is.na(eGFR) & firstSampleDate >= final_date - years(5) & firstSampleDate <= final_date, na.rm = TRUE),
      
      neerud_risk = case_when(
        neerud_labor_egfr | neerud_labor_loinc ~ 1,   
        neerud_on_mootmine ~ 0,                  
        TRUE ~ 2
      ),
    
      # ---- suitsetamine ----
      lahim_suits = case_when(
        any(!is.na(lastSmokingStatus) & smokingDate > final_date, na.rm = TRUE) ~ 
          lastSmokingStatus[which.min(if_else(!is.na(lastSmokingStatus) & smokingDate > final_date, smokingDate, as.Date("2100-01-01")))],
        TRUE ~ NA
      ),
      suits_risk = case_when(
        any(lastSmokingStatus %in% c("Former", "Current") & smokingDate <= final_date, na.rm = TRUE) ~ 1,
        any(!is.na(lastSmokingStatus) & smokingDate <= final_date, na.rm = TRUE) ~ 0,
        
        lahim_suits %in% c("Former", "Current") ~ 1,
        !is.na(lahim_suits) ~ 0,
        
        TRUE ~ 2
      ),
      
      # ---- kmi ----
      viimane_kmi = case_when(
        any(!is.na(lastBmi), na.rm = TRUE) ~ 
          lastBmi[which.min(if_else(!is.na(lastBmi), abs(as.numeric(as.Date(bmiDate)-as.Date(final_date))), Inf))],
        TRUE ~ NA
      ) 
    ) %>%
    
    ungroup()
    
  return(riskid_df)
}

# glükoosimõõtmistest tuleneva riski defineerimine
f_arvuta_glykoos_risk <- function(df) {
  df %>%
    transmute(
      skood, 
      final_date,
      
      loinc, 
      tulemus_loinc = as.numeric(measurementResultNumeric),
      kuup_loinc = measurementDate,
      
      tulemus_nmr = as.numeric(Glucose),
      kuup_nmr = firstSampleDate
    ) %>%
    
    pivot_longer(
      cols = c(tulemus_loinc, tulemus_nmr, kuup_loinc, kuup_nmr),
      names_to = c(".value", "source"),
      names_pattern = "(tulemus|kuup)_(loinc|nmr)"
    ) %>%
    
    mutate(
      on_sobiv = !is.na(tulemus),
      
      sobiv = case_when(
        source == "loinc" ~ loinc %in% c("14771-0", "14749-6", "14770-2", "72516-8") &
          tulemus >= 7.0 &
          kuup >= final_date - years(5) &
          kuup <= final_date,
        source == "nmr" ~ tulemus >= 7.0 & 
          kuup >= final_date - years(5) &
          kuup <= final_date,
        TRUE ~ FALSE
      )
    ) %>%
    
    group_by(skood) %>%
    
    summarise(
      kokku_teste = sum(on_sobiv, na.rm = TRUE),
      korgeid_naitajaid = sum(sobiv, na.rm = TRUE),
      korge_nait_kuupaevi = n_distinct(kuup[sobiv == TRUE], na.rm = TRUE),
      .groups = "drop"
    ) %>%
    
    mutate(
      glykoos_risk = case_when(
        (korgeid_naitajaid > 1 & korge_nait_kuupaevi > 1) ~ 1,
        kokku_teste == 0 ~ 2, 
        TRUE ~ 0
      ),
      
      diab_glyk_labor = (korgeid_naitajaid > 1 & korge_nait_kuupaevi > 1),
      diab_glyk_on_mootmine = (kokku_teste > 0)
    )
}

# lõpliku diabeediriski defineerimine
f_loplik_diabeedirisk <- function(skoor_df, glykoos_df) {
  skoor_df %>% 
    left_join(glykoos_df, by = "skood") %>%
    
    mutate(
      glykoos_risk = replace_na(glykoos_risk, 2),
      diab_glyk_labor = replace_na(diab_glyk_labor, FALSE),
      diab_glyk_on_mootmine = replace_na(diab_glyk_on_mootmine, FALSE)
    ) %>%
    
    group_by(skood) %>%
    mutate(diabeet_risk = case_when(
      any(diabeet_risk1 == 1 | glykoos_risk == 1) ~ 1,
      all(diabeet_risk1 == 2 & glykoos_risk == 2) ~ 2,
      TRUE ~ 0
    ),
    
    diab_lopp_labor = any(diab_labor_HbA1c, na.rm = TRUE) | any(diab_glyk_labor, na.rm = TRUE),
    diab_lopp_mootmine = any(diab_on_mootmine, na.rm = TRUE) | any(diab_glyk_on_mootmine, na.rm = TRUE)
    ) %>%
    
    ungroup()
}

# lõplikud riskifaktorid
f_loplikud_riskifaktorid <- function(df) {
  df %>%
    group_by(skood) %>%
    slice(1) %>%
    ungroup()
}

# --------------------------------------------
# --------- riskifaktorite kokkuvõte ---------
# --------------------------------------------
# millise kriteeriumi alusel määrati riskigrupp
f_kriteeriumite_koondtbl <- function(df_riskid){
  kriteeriumid <- list(
    "Hüpertensioon" = c("hyp_diag", "hyp_ravim", "hyp_labor_koos"),
    "Hüperlipideemia" = c("lip_diag", "lip_ravim", "lip_labor_ldl", "lip_labor_koles"),
    "Diabeet" = c("diab_diag", "diab_ravim", "diab_labor_HbA1c", "diab_glyk_labor"),
    "Neerufunktsioon" = c("neerud_labor_egfr", "neerud_labor_loinc")
  )
  
  koond <- lapply(names(kriteeriumid), function(haigus) {
    risk_veerg <- case_when(
      haigus == "Hüpertensioon" ~ "hypertensioon_risk",
      haigus == "Hüperlipideemia" ~ "hyperlipideemia_risk",
      haigus == "Diabeet" ~ "diabeet_risk",
      haigus == "Neerufunktsioon" ~ "neerud_risk"
    )
  
  df_riskid %>%
     group_by(sex.x, Riskigrupp = .data[[risk_veerg]]) %>%
     summarise(
       Patsiente_kokku = n(),
       across(all_of(kriteeriumid[[haigus]]), ~sum(.x, na.rm = TRUE)),
       .groups = "drop"
     ) %>%
     
     pivot_longer(
       cols = all_of(kriteeriumid[[haigus]]),
       names_to = "Kriteerium",
       values_to = "N"
     ) %>%
    mutate(Haigus = haigus)
  }) %>%
    bind_rows()
  
  koond <- koond %>%
    select(Haigus, sex.x, Riskigrupp, Patsiente_kokku, Kriteerium, N) %>%
    arrange(Haigus, Riskigrupp)
  
  return(koond)
}

# riskifaktorite esinemine
f_kokkuvote <- function(x, type = c("pidev", "bin")) {
  type <- match.arg(type)
  
  if (type == "pidev") {
    m <- mean(x, na.rm = TRUE)
    s <- sd(x, na.rm = TRUE)
    return(sprintf("%.1f (%.1f)", m, s))
  } else {
    n <- sum(x == 1, na.rm = TRUE)
    kokku <- sum(!is.na(x))
    protsent <- 100*n/kokku
    return(sprintf("%d (%.1f%%)", n, protsent))
  }
}

# riskifaktorite esinemine tabelina
f_kokkuvotte_tbl <- function(risk_df) {
  df_naised <- risk_df[risk_df$sex.x == "female", ]
  df_mehed <- risk_df[risk_df$sex.x == "male", ]
  
  arvuta <- function(df) {
    c(
      f_kokkuvote(df$sex.x == "female", "bin"),
      f_kokkuvote(df$vanus_indeksil, "pidev"),
      f_kokkuvote(df$hypertensioon_risk, "bin"),
      f_kokkuvote(df$hyperlipideemia_risk, "bin"),
      f_kokkuvote(df$diabeet_risk, "bin"),
      f_kokkuvote(df$suits_risk, "bin"),
      f_kokkuvote(df$viimane_kmi, "pidev"),
      f_kokkuvote(df$neerud_risk, "bin")
    )
  }
  
  data.frame(
    muutuja = c("Sugu (naine)", "Vanus", "Hüpertensioon", "Hüperlipideemia", "Diabeet", "Suitsetamine", "KMI", "Neerufunktsioon"),
    kokku = arvuta(risk_df),
    naised = arvuta(df_naised),
    mehed = arvuta(df_mehed)
  )
}

# --------------------------------------------
# ----- riskifaktorid KS-ga doonoritele ------
# --------------------------------------------
cac_riskid_alg <- f_arvuta_riskid(
  skoorid_df = cac_puhas,
  diagnoosid_ravimid_df = diagnoosid_alg,
  mootmised_df = mootmised_puhas,
  kmi_suits_df = kmi_suits_alg,
  egfr_df = tbl_egfr
)

mootmised_glykoos <- f_arvuta_glykoos_risk(cac_riskid_alg)

cac_riskid <- f_loplik_diabeedirisk(cac_riskid_alg, mootmised_glykoos)

# lõplikud riskifaktorid
cac_riskifaktorid <- f_loplikud_riskifaktorid(cac_riskid)

# prs-de lisamine
prs_cac <- inner_join(cac_riskifaktorid, prs, by = "skood")
sum(is.na(prs_cac$raw_PGS003726)) # 27 inimesel pole kasutatavat PRS-i

prs_cac <- prs_cac %>%
  filter(!is.na(raw_PGS003726))

sum(is.na(prs_cac$viimane_kmi))
prs_cac <- prs_cac %>%
  filter(!is.na(viimane_kmi))

# ---- valimi suuruse kontroll ----
krit_tbl_cac <- f_kriteeriumite_koondtbl(prs_cac)
table(prs_cac$hyp_labor_on_analyys)

sapply(prs_cac[, c("hypertensioon_risk", "hyperlipideemia_risk", "diabeet_risk", "suits_risk", "neerud_risk")], table)
by(prs_cac[, c("hypertensioon_risk", "hyperlipideemia_risk", "diabeet_risk", "suits_risk", "neerud_risk")], prs_cac$sex.x, function(x) {
  sapply(x, table)
})

cac_kokkuvote_tbl <- f_kokkuvotte_tbl(prs_cac)
cac_kokkuvote_tbl

# kõigi riskifaktorite olemasolu/ eksisteerimine
prs_cac %>% filter(hypertensioon_risk == 1, hyperlipideemia_risk == 1, diabeet_risk == 1, neerud_risk == 1, suits_risk == 1, is.na(lastBmi))
prs_cac %>% filter(hypertensioon_risk != 2, hyperlipideemia_risk != 2, diabeet_risk != 2, neerud_risk != 2, suits_risk != 2, is.na(lastBmi))

# kmi ajaline erinevus (kmi mõõtmine ja indeks)
cac_kmi_aeg <- prs_cac %>%
  mutate(erinevus_p = as.numeric(difftime(bmiDate, final_date, units = "days")))

ggplot(cac_kmi_aeg, aes(x = erinevus_p)) +
  geom_histogram(na.rm = TRUE) +
  theme_minimal() +
  labs(x = "Erinevus päevades", y = "Sagedus")


##################################################
# ------------ Kirjeldav statistika -------------
##################################################
nrow(prs_cac)

# üldine
prs_cac %>%
  mutate(cac_riskigrupp = case_when(
    value >= 1000 ~ "Väga kõrge",
    value >= 300 ~ "Kõrge",
    value >= 100 ~ "Keskmine",
    TRUE ~ "Madal"
  )) %>%
  count(cac_riskigrupp)

sum(prs_cac$value == 0)

# naistele
prs_cac %>%
  filter(sex.x == "female") %>%
  mutate(cac_riskigrupp = case_when(
    value >= 1000 ~ "Väga kõrge",
    value >= 300 ~ "Kõrge",
    value >= 100 ~ "Keskmine",
    TRUE ~ "Madal"
  )) %>%
  count(cac_riskigrupp)

sum(prs_cac$value == 0 & prs_cac$sex.x == "female")

# meestele
prs_cac %>%
  filter(sex.x == "male") %>%
  mutate(cac_riskigrupp = case_when(
    value >= 1000 ~ "Väga kõrge",
    value >= 300 ~ "Kõrge",
    value >= 100 ~ "Keskmine",
    TRUE ~ "Madal"
  )) %>%
  count(cac_riskigrupp)

sum(prs_cac$value == 0 & prs_cac$sex.x == "male")

# KS jaotus
ks_jaotus <- ggplot(prs_cac, aes(x = value)) +
  geom_histogram(bins = 50) +
  theme_minimal() +
  theme(text = element_text(size = 10)) +
  labs(x = "Kaltsiumskoor", y = "Inimeste arv")

ks_jaotus
ggsave("ks_jaotus.png", plot = ks_jaotus, width = 6, height = 4, dpi = 300)


##################################################
# ----------- Mudelid riskiskooridele ------------
##################################################
# defineeritud riskifaktorite eeltöötlus enne mudelit
f_tbl_mudeliks <- function(riskiskoorid_df) {
   tbl <- riskiskoorid_df %>%
    mutate(
      hypertensioon_risk = factor(hypertensioon_risk, levels = c(0, 1, 2)),
      hyperlipideemia_risk = factor(hyperlipideemia_risk, levels = c(0, 1, 2)),
      diabeet_risk = factor(diabeet_risk, levels = c(0, 1, 2)),
      suits_risk = factor(suits_risk, levels = c(0, 1, 2)),
      neerud_risk = factor(neerud_risk, levels = c(0, 1, 2))
    )
   return(tbl)
}

# --------------------------------------------
# --------- mudel kaltsiumiskoorile ----------
# --------------------------------------------
cac_mudel_andmed <- f_tbl_mudeliks(prs_cac)
cac_mudel_andmed$cacs_bin <- ifelse(cac_mudel_andmed$value > 0, 1, 0)

# ---- kaheastmelised mudelid ----
# vanus x KS
ggplot(cac_mudel_andmed, aes(x = vanus_indeksil, y = value)) +
  geom_point() +
  theme_minimal()

# vanus x log(KS)
j_vanus <- ggplot(cac_riskifaktorid, aes(x = vanus_indeksil, y = log(value+1))) +
  geom_point() +
  theme_minimal() + 
  theme(text = element_text(size = 10)) +
  labs(x = "Vanus",
       y = "Kaltsiumiskoor (logaritmitud)")

j_vanus
ggsave("j_vanus.png", plot = j_vanus, width = 6, height = 4, dpi = 300)

ggplot(cac_riskifaktorid, aes(x = vanus_indeksil, y = log(value+1), color = sex.x)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~ sex.x) +
  theme_minimal() +
  theme(text = element_text(size = 10)) +
  labs(x = "Vanus",
       y = "Kaltsiumiskoor (logaritmitud)")

# ---- logistiline mudel
m_bin <- glm(cacs_bin ~  sex.x + vanus_indeksil + neerud_risk + viimane_kmi + suits_risk + hypertensioon_risk + diabeet_risk + hyperlipideemia_risk + sex.x * vanus_indeksil, family = binomial, data = cac_mudel_andmed, na.action = na.exclude)
summary(m_bin)
exp(coef(m_bin))
exp(confint(m_bin))

# eelduste kontroll
vif(m_bin, type = "terms") # multikollineaarsuse kontroll
roc_kover <- roc(m_bin$y, m_bin$fitted.values)
auc(roc_kover)
plot(roc_kover)
hoslem.test(m_bin$y, m_bin$fitted.values, g = 10)
pR2(m_bin)

# ---- lineaarsed mudelid
m_pos <- lm(log(value) ~ sex.x + vanus_indeksil + neerud_risk + viimane_kmi + suits_risk + hypertensioon_risk + diabeet_risk + hyperlipideemia_risk + sex.x * vanus_indeksil, data = cac_mudel_andmed[cac_mudel_andmed$value > 0,], na.action = na.exclude)
summary(m_pos)
AIC(m_pos)

# parima splanidega mudeli leidmine
m_pos_3 <- lm(log(value) ~ sex.x * ns(vanus_indeksil, df = 3) + sex.x + ns(vanus_indeksil, df = 3)  + neerud_risk + viimane_kmi + suits_risk + hypertensioon_risk + diabeet_risk + hyperlipideemia_risk, data = cac_mudel_andmed[cac_mudel_andmed$value > 0,], na.action = na.exclude)
m_pos_4 <- lm(log(value) ~ sex.x * ns(vanus_indeksil, df = 4) + sex.x + ns(vanus_indeksil, df = 4) + neerud_risk + viimane_kmi + suits_risk + hypertensioon_risk + diabeet_risk + hyperlipideemia_risk, data = cac_mudel_andmed[cac_mudel_andmed$value > 0,], na.action = na.exclude)
m_pos_5 <- lm(log(value) ~ sex.x * ns(vanus_indeksil, df = 5) + sex.x + ns(vanus_indeksil, df = 5) + neerud_risk + viimane_kmi + suits_risk + hypertensioon_risk + diabeet_risk + hyperlipideemia_risk, data = cac_mudel_andmed[cac_mudel_andmed$value > 0,], na.action = na.exclude)
summary(m_pos_3)
summary(m_pos_4)
summary(m_pos_5)
AIC(m_pos_3, m_pos_4, m_pos_5)
summary(m_pos_3)$adj.r.squared
summary(m_pos_4)$adj.r.squared
summary(m_pos_5)$adj.r.squared
# vähim AIC on df = 3 korral

# vanuse ruudu ja kuubi testimine
m_pos_ruut <- lm(log(value) ~ sex.x + vanus_indeksil + I(vanus_indeksil^2) + neerud_risk + viimane_kmi + suits_risk + hypertensioon_risk + diabeet_risk + hyperlipideemia_risk + sex.x * vanus_indeksil, data = cac_mudel_andmed[cac_mudel_andmed$value > 0,], na.action = na.exclude)
m_pos_kuup <- lm(log(value) ~ sex.x + vanus_indeksil + I(vanus_indeksil^2) + I(vanus_indeksil**3) + neerud_risk + viimane_kmi + suits_risk + hypertensioon_risk + diabeet_risk + hyperlipideemia_risk + sex.x * vanus_indeksil, data = cac_mudel_andmed[cac_mudel_andmed$value > 0,], na.action = na.exclude)
summary(m_pos_ruut)
summary(m_pos_kuup)
AIC(m_pos_ruut, m_pos_kuup)
summary(m_pos_ruut)$adj.r.squared
summary(m_pos_kuup)$adj.r.squared

# lineaarsete mudelite võrdlus
anova(m_pos, m_pos_ruut) # ruuduga on parem
anova(m_pos, m_pos_kuup) # kuubiga mudel on parem
anova(m_pos, m_pos_3) # splainiga mudel ei ole parem

exp(coef(m_pos_3))
exp(confint(m_pos_3))
# mudeli eelduste kontroll
vif(m_pos_3)
plot(m_pos_3, which = 2)
plot(m_pos_3, which = 1)

# splainiga argumentide olulisuse kontroll
Anova(m_pos_3, type = 3)

# ---- diskordantsus (jäägid) ----
res_bin <- residuals(m_bin, type = "response")
res_pos <- residuals(m_pos_3)

cac_mudel_andmed$disk_cac <- NA
cac_mudel_andmed$disk_cac[cac_mudel_andmed$value == 0] <- res_bin[cac_mudel_andmed$value == 0]
cac_mudel_andmed$disk_cac[cac_mudel_andmed$value > 0] <- scale(res_pos)

cac_mudel_andmed$p_mittenull <- predict(m_bin, newdata = cac_mudel_andmed, type = "response")
cac_mudel_andmed$mu_log <- predict(m_pos_3, newdata = cac_mudel_andmed)
cac_mudel_andmed$oodatud <- cac_mudel_andmed$p_mittenull * exp(cac_mudel_andmed$mu_log)
cac_mudel_andmed$jaagid <- log(cac_mudel_andmed$value + 1) - log(cac_mudel_andmed$oodatud + 1)

ggplot(cac_mudel_andmed, aes(x = jaagid)) +
  geom_histogram() + 
  theme_minimal() +
  labs(y = "arv")

# kuidas inimesed grupeeruvad diskordantsuse järgi
cac_mudel_andmed$varvid <- ifelse(
  cac_mudel_andmed$oodatud >= cac_mudel_andmed$value, 
  "Resistentne (oodatud > tegelik)", 
  "Vastuvõtlik (oodatud < tegelik)"
)

table(cac_mudel_andmed$varvid)
table(cut(cac_mudel_andmed$jaagid, breaks = c(-Inf, -2, 2, Inf)))

j_disk <- ggplot(cac_mudel_andmed, aes(x = log(oodatud+1), y = log(value+1))) +
  geom_point(aes(color = varvid), alpha = 0.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_abline(slope = 1, intercept = -2, linetype = "dotted") +
  geom_abline(slope = 1, intercept = 2, linetype = "dotted") +
  scale_color_manual(values = c(
    "Resistentne (oodatud > tegelik)" = "lightblue",
    "Vastuvõtlik (oodatud < tegelik)" = "lightpink"
  )) +
  theme_minimal() +
  theme(text = element_text(size = 10), legend.position = "bottom") +
  guides(color = guide_legend(ncol = 2)) +
  labs(color = "Diskordantsus", 
       x = "Oodatud haiguskoormus (logaritmitud)",
       y = "Kaltsiumiskoor (logaritmitud)")

j_disk
ggsave("j_disk.png", plot = j_disk, width = 7, height = 4, dpi = 300)

##################################################
# --------------- Mudelid jääkidele --------------
##################################################
# --------------------------------------------
#  --------- detailne profileerimine ---------
# --------------------------------------------
m_prs <- lm(jaagid ~ raw_PGS003726, data = cac_mudel_andmed)
summary(m_prs)

# kvintiilide leidmine ja referentsgrupi määramine
cac_mudel_andmed <- cac_mudel_andmed %>%
  mutate(
    kvint = ntile(percent_rank_PGS003726, 5),
    kvint = factor(kvint),
    kvint = relevel(kvint, ref = "3")
  )

m_prs_kvint <- lm(jaagid ~ kvint, data = cac_mudel_andmed)
summary(m_prs_kvint)

puhas_kvint <- tidy(m_prs_kvint, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(term = gsub("kvint", "kvintiil ", term))

# kui on suurem prs, siis on suurem diskordantsus
j_kvint <- ggplot(puhas_kvint, aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(size = 0.4) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(data = data.frame(term = "kvintiil 3", estimate = 0),
             aes(x = term, y = estimate),
             inherit.aes = FALSE,
             shape = 18,
             size = 4) +
  theme_minimal() +
  theme(text = element_text(size = 10)) +
  labs(x = "PRS-i kvintiil", y = "Parameeter")

j_kvint
ggsave("j_kvint.png", plot = j_kvint, width = 6, height = 4, dpi = 300)

# --------------------------------------------
#  ---------- reklassifitseerimine -----------
# --------------------------------------------
# ---- kaheastmeline mudel prs-ga ----
# ---- logistiline mudel
m_bin_prs <- glm(cacs_bin ~  sex.x + vanus_indeksil + neerud_risk + viimane_kmi + suits_risk + hypertensioon_risk + diabeet_risk + hyperlipideemia_risk + sex.x * vanus_indeksil + raw_PGS003726, family = binomial, data = cac_mudel_andmed, na.action = na.exclude)
summary(m_bin_prs)
exp(coef(m_bin_prs))
exp(confint(m_bin_prs))

# eelduste kontroll
vif(m_bin_prs) # multikollineaarsuse kontroll
roc_kover <- roc(m_bin_prs$y, m_bin_prs$fitted.values)
auc(roc_kover)
plot(roc_kover)
pR2(m_bin_prs)
hoslem.test(m_bin_prs$y, m_bin_prs$fitted.values, g = 10)

# ---- lineaarne mudel
m_pos_prs <- lm(log(value) ~ sex.x * ns(vanus_indeksil, df = 3) + sex.x + ns(vanus_indeksil, df = 3)  + neerud_risk + viimane_kmi + suits_risk + hypertensioon_risk + diabeet_risk + hyperlipideemia_risk + raw_PGS003726, data = cac_mudel_andmed[cac_mudel_andmed$value > 0,], na.action = na.exclude)
summary(m_pos_prs)
exp(coef(m_pos_prs))
exp(confint(m_pos_prs))

# mudeli eelduste kontroll
vif(m_pos_prs)
plot(m_pos_prs, which = 2)
plot(m_pos_prs, which = 1)

# splainiga argumentide olulisuse kontroll
Anova(m_pos_prs, type = 3)

# ---- võrdlus esialgse kaheastmelise mudeliga
anova(m_bin, m_bin_prs)
anova(m_pos_3, m_pos_prs)

# ---- diskordantsus (jäägid) ----
res_bin_prs <- residuals(m_bin_prs, type = "response")
res_pos_prs <- residuals(m_pos_prs)

cac_mudel_andmed$disk_cac_prs <- NA
cac_mudel_andmed$disk_cac_prs[cac_mudel_andmed$value == 0] <- res_bin_prs[cac_mudel_andmed$value == 0]
cac_mudel_andmed$disk_cac_prs[cac_mudel_andmed$value > 0] <- scale(res_pos_prs)

cac_mudel_andmed$p_mittenull_prs <- predict(m_bin_prs, newdata = cac_mudel_andmed, type = "response")
cac_mudel_andmed$mu_log_prs <- predict(m_pos_prs, newdata = cac_mudel_andmed)
cac_mudel_andmed$oodatud_prs <- cac_mudel_andmed$p_mittenull_prs * exp(cac_mudel_andmed$mu_log_prs)
cac_mudel_andmed$jaagid_prs <- log(cac_mudel_andmed$value + 1) - log(cac_mudel_andmed$oodatud_prs + 1)

ggplot(cac_mudel_andmed, aes(x = jaagid_prs)) +
  geom_histogram() + 
  theme_minimal() +
  labs(y = "arv")

# kuidas inimesed grupeeruvad diskordantsuse järgi
cac_mudel_andmed$varvid <- ifelse(
  cac_mudel_andmed$oodatud_prs >= cac_mudel_andmed$value, 
  "Resistentne (oodatud > tegelik)", 
  "Vastuvõtlik (oodatud < tegelik)"
)

table(cac_mudel_andmed$varvid)
table(cut(cac_mudel_andmed$jaagid_prs, breaks = c(-Inf, -2, 2, Inf)))

j_prs_disk <- ggplot(cac_mudel_andmed, aes(x = log(oodatud_prs+1), y = log(value+1))) +
  geom_point(aes(color = varvid), alpha = 0.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_abline(slope = 1, intercept = -2, linetype = "dotted") +
  geom_abline(slope = 1, intercept = 2, linetype = "dotted") +
  scale_color_manual(values = c(
    "Resistentne (oodatud > tegelik)" = "lightblue",
    "Vastuvõtlik (oodatud < tegelik)" = "lightpink"
  )) +
  theme_minimal() +
  theme(text = element_text(size = 10), legend.position = "bottom") +
  guides(color = guide_legend(ncol = 2)) +
  labs(color = "Diskordantsus", 
       x = "Oodatud haiguskoormus (logaritmitud)",
       y = "Kaltsiumiskoor (logaritmitud)")
j_prs_disk
ggsave("j_prs_disk.png", plot = j_prs_disk, width = 7, height = 4, dpi = 300)

# ---- doonorite reklassifitseerimine prs-i tõttu
reklass_andmed <- cac_mudel_andmed %>%
  mutate(
    vana_grupp = case_when(
      jaagid <= -2 ~ "z < -2",
      jaagid > 2 ~ "z > 2",
      TRUE ~ "-2 < z < 2"
    ),
    uus_grupp = case_when(
      jaagid_prs <= -2 ~ "z < -2",
      jaagid_prs > 2 ~ "z > 2",
      TRUE ~ "-2 < z < 2"
    )
  ) %>%
  count(vana_grupp, uus_grupp, name = "inimesi") %>%
  mutate(
    suund = case_when(
      vana_grupp == uus_grupp ~ "Muutusteta",
      vana_grupp == "z < -2" & uus_grupp %in% c("-2 < z < 2", "z > 2") ~ "Risk väheneb",
      vana_grupp == "-2 < z < 2" & uus_grupp == "z > 2" ~ "Risk väheneb",
      TRUE ~ "Risk suureneb"
    ),
    vana_grupp = factor(vana_grupp, levels = c("z < -2", "-2 < z < 2","z > 2")),
    uus_grupp = factor(uus_grupp, levels = c("z < -2", "-2 < z < 2","z > 2"))
  )
reklass_andmed

j_reklass <- ggplot(reklass_andmed, aes(axis1 = vana_grupp, axis2 = uus_grupp, y = inimesi)) +
  geom_alluvium(aes(fill = suund), width = 1/8, alpha = 0.75, color = "white", linewidth = 0.5) +
  geom_stratum(width = 1/8, color = "black", fill = "gray") +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), size = 4) +
  scale_x_discrete(limits = c("Baasmudel", "PRS-ga mudel"), expand = c(0.15, 0.15)) +
  scale_fill_manual(values = c("Risk suureneb" = "pink", "Risk väheneb" = "lightgreen", "Muutusteta" = "gray")) +
  theme_minimal() +
  theme(text = element_text(size = 10), panel.grid.major.x = element_blank(), panel.grid.major.y = element_blank()) +
  labs(y = "Patsiente", fill = "Riski muutus")

j_reklass
ggsave("j_reklass.png", plot = j_reklass, width = 6, height = 4, dpi = 300)