Inspecting_data.Rmd

---
title: "Inspecting data"
output: html_notebook
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, results = 'hide', message = FALSE, warning = FALSE)
```
---

Setting up everything
```{r}
options(repos = "http://nexus.ki.se/repository/cran.r-project.org")
.libPaths("W:/C6_AmorisGuest/Sara Hagg/Karolina Gustavsson/shun_trial/lib")

library(tidyverse)
library(flexsurv)
library(haven)
library(survival)
library(broom)
library(ggplot2)
library(checkmate) 
library(htmlTable)
library(survey)
library(tidyr)
library(dplyr)
library(psych)
```

Read the data
```{r}
amoris <- read_sas("W:/C6_AmorisGuest/Sara Hagg/Karolina Gustavsson/Data/StudyPop.sas7bdat")

#Optional
#load("W:/C6_AmorisGuest/Sara Hagg/Karolina Gustavsson/BioAge-master/BioAge/data/NHANES3.rda")
```

Biomarkers
```{r}
#select our 17 biomarkers
biomarkers <- c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")
```


Distribution pre processesing 
```{r}
library(ggplot2)
library(tidyr)

long_data <- gather(amoris, key = "biomarker", value = "value", biomarkers)

# Plot histograms
ggplot(long_data, aes(x = value)) +
  geom_histogram(bins = 50, fill = "blue", color = "black") +
  facet_wrap(~ biomarker, scales = "free", ncol = 4) +
  theme_minimal() +
  labs(x = "Value", y = "Frequency", title = "Histograms of Biomarkers")

rm(long_data)
```
Density plots before 
```{r}
library(ggplot2)
library(tidyr)

long_data <- gather(amoris, key = "biomarker", value = "value", biomarkers)

# Plot density plots
ggplot(long_data, aes(x = value, fill = biomarker)) +
  geom_density(alpha = 0.7,adjust=4) +
  facet_wrap(~ biomarker, scales = "free", ncol = 4) +
  theme_minimal() +
  labs(x = "Value", y = "Density", title = "Density Plots of Biomarkers") +
  theme(legend.position = "none") # Hides the legend, adjust as needed

rm(long_data)
```
Correlation with age for entire cohort
```{r}
correlations_spearman <- cor(amoris$Alder, amoris[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "spearman", use = "complete.obs")
correlations_pearson <- cor(amoris$Alder, amoris[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "pearson", use = "complete.obs")

print(correlations_spearman)
print(correlations_pearson)

#Over 0.1 Spearman:
#S_Alb, S_Krea, TC, TG, S_K, S_P, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fs_Gluk

#Over 0.15 Spearman:
#S_Alb, TC, TG, S_P, S_Urea, S_LD, S_Alp, S_Urat, S_FAMN, S_Hapt, fs_Gluk

#Over 0.2 Spearman:
#S_Alb, TC, TG, S_Urea, S_LD, S_Hapt, fs_Gluk

#Over 0.1 Pearson:
#S_Alb, S_Krea, TC, TG, S_K, S_P, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fs_Gluk

#Over 0.15 Pearson:
#S_Alb, TC, S_P, S_Urea, S_LD, S_Urat, S_FAMN, S_Hapt, fs_Gluk

#Over 0.2 Pearson:
#S_Alb, TC, S_Urea, S_LD, fs_Gluk
```

Correlation separately by sex
```{r}
# Filter data based on the 'Kon' column
amoris_men <- amoris[amoris$Kon == 1, ]
amoris_women <- amoris[amoris$Kon == 2, ]


## check pearson correlation
# Calculate Pearson correlations between Age and biomarkers
correlations_pearson_m <- cor(amoris_men$Alder, amoris_men[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "pearson", use = "complete.obs")
correlations_spearman_m <- cor(amoris_men$Alder, amoris_men[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "spearman", use = "complete.obs")
correlations_pearson_f <- cor(amoris_women$Alder, amoris_women[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "pearson", use = "complete.obs")
correlations_spearman_f <- cor(amoris_women$Alder, amoris_women[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "spearman", use = "complete.obs")
print("correlations_pearson_m:")
print(correlations_pearson_m)#Fe_maet, fS_TIBIC, fS_Jaern, S_Alp are under abs(0.1)
print("correlations_spearman_m:")
print(correlations_spearman_m)#Fe_maet, fS_TIBIC, fS_Jaern, S_Alp are under abs(0.1)
print("correlations_pearson_f:")
print(correlations_pearson_f)#S_P,Fe_maet,S_Ca (depends on rounding), fS_Jarn, are under abs(0.1)
print("correlations_spearman_f:")
print(correlations_spearman_f)#S_P,Fe_maet, S_Ca (depends on rounding), fS_Jarn, are under abs(0.1)

#print(abs(correlations_pearson_m) > 0.1)

#"S_Alb", "S_Krea", "TC", "TG", "S_K", "S_Urea", "S_LD", "S_Ca", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk": the ones over 0.1 rounded up, separatley for each sex

#pearson_m > 0.1: S_Alb, S_Krea, TC, TG, S_K, S_P, S_Urea, S_LD, S_Ca, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#pearson_m > 0.15: S_Alb, S_Krea, TC, S_P, S_Urea, S_Ca, S_Hapt, fS_Gluk
#pearson_m > 0.2: S_Alb, TC, S_Urea, S_Ca, fS_Gluk

#spearman_m > 0.1: S_Alb, S_Krea, TC, TG, S_K, S_P, S_Urea, S_LD, S_Ca, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#spearman_m > 0.15: S_Alb, S_Krea, TC, TG, S_P, S_Urea, S_Ca, S_Hapt, fS_Gluk
#spearman_m > 0.2: S_Alb, TC, S_P, S_Urea, S_Ca, S_Hapt

#pearson_f >0.1: S_Alb, S_Krea, TC, TG, S_K, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#pearson_f > 0.15: S_Alb, S_Krea, TC, TG, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#pearson_f > 0.2: S_Alb, TC, TG, S_Urea, S_LD, S_Alp, S_Urat, fS_Gluk

#spearman_f >0.1: S_Alb, S_Krea, TC, TG, S_K, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#spearman_f > 0.15: S_Alb, S_Krea, TC, TG, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#spearman_f > 0.2: S_Alb, TC, TG, S_Urea, S_LD, S_Alp, S_Urat, fS_Gluk

```
```{r}
rm(amoris_men)
rm(amoris_women)
rm(correlations_pearson)
rm(correlations_pearson_m)
rm(correlations_pearson_f)
rm(correlations_spearman)
rm(correlations_spearman_m)
rm(correlations_spearman_f)
```

Summary statistics for cont. variables
```{r}
# Identify continuous (numeric) variables
continuous_vars <- sapply(amoris, is.numeric) & !sapply(amoris, is.factor)

# Custom summary function including standard deviation
custom_summary <- function(x) {
  c(Mean = mean(x, na.rm = TRUE),
    SD = sd(x, na.rm = TRUE),
    Median = median(x, na.rm = TRUE),
    '1st Qu' = quantile(x, probs = 0.25, na.rm = TRUE),
    '3rd Qu' = quantile(x, probs = 0.75, na.rm = TRUE),
    Min = min(x, na.rm = TRUE),
    Max = max(x, na.rm = TRUE))
}

# Apply custom summary function to each continuous variable
summary_stats <- lapply(amoris[continuous_vars], custom_summary)

# Convert the list to a dataframe for a cleaner presentation
summary_df <- do.call(rbind, summary_stats)

# View the summary statistics
print(summary_df)
```


========================================================================================================
PROCESSING
========================================================================================================



Filtering out outliers (5 Standard deviations) ##### This should be changed to sep per sex
```{r}
library(dplyr)

AMORIS <- amoris %>%
  group_by(Kon) %>%
  filter(
    if_else(
      Kon == 1, 
      across(all_of(biomarkers), ~ .x > mean(.x, na.rm = TRUE) - 5 * sd(.x, na.rm = TRUE) & .x < mean(.x, na.rm = TRUE) + 5 * sd(.x, na.rm = TRUE)),
      across(all_of(biomarkers), ~ .x > mean(.x, na.rm = TRUE) - 5 * sd(.x, na.rm = TRUE) & .x < mean(.x, na.rm = TRUE) + 5 * sd(.x, na.rm = TRUE))
    )
  ) %>%
  ungroup()

AMORIS <- AMORIS[complete.cases(AMORIS[, biomarkers]), ]
```

Checking histograms again
```{r}
library(ggplot2)
library(tidyr)

long_data <- gather(AMORIS, key = "biomarker", value = "value", biomarkers)

# Plot histograms
ggplot(long_data, aes(x = value)) +
  geom_histogram(bins = 50, fill = "blue", color = "black") +
  facet_wrap(~ biomarker, scales = "free", ncol = 4) +
  theme_minimal() +
  labs(x = "Value", y = "Frequency", title = "Histograms of Biomarkers")

rm(long_data)
```
Checking density plots again
```{r}
library(ggplot2)
library(tidyr)

long_data <- gather(AMORIS, key = "biomarker", value = "value", biomarkers)

# Plot density plots
ggplot(long_data, aes(x = value, fill = biomarker)) +
  geom_density(alpha = 0.7,adjust=4) +
  facet_wrap(~ biomarker, scales = "free", ncol = 4) +
  theme_minimal() +
  labs(x = "Value", y = "Density", title = "Density Plots of Biomarkers") +
  theme(legend.position = "none") # Hides the legend, adjust as needed

rm(long_data)
```

Correlation with age for entire cohort
```{r}
correlations_spearman <- cor(AMORIS$Alder, AMORIS[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "spearman", use = "complete.obs")
correlations_pearson <- cor(AMORIS$Alder, AMORIS[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "pearson", use = "complete.obs")

print(correlations_spearman)
print(correlations_pearson)

#Over 0.1 Spearman:
#S_Alb, S_Krea, TC, TG, S_K, S_P, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fs_Gluk

#Over 0.15 Spearman:
#S_Alb, TC, TG, S_P, S_Urea, S_LD, S_Alp, S_Urat, S_Hapt, fs_Gluk
# NOT S_FAMN

#Over 0.2 Spearman:
#S_Alb, TC, TG, S_Urea, S_LD, fs_Gluk

#Over 0.1 Pearson:
#S_Alb, S_Krea, TC, TG, S_K, S_P, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fs_Gluk

#Over 0.15 Pearson:
#S_Alb, TC, TG, S_P, S_Urea, S_LD, S_Alp, S_Urat, S_Hapt, fs_Gluk

#Over 0.2 Pearson:
#S_Alb, TC, S_Urea, S_LD, fs_Gluk

```

Correlation separately per sex
```{r}
# Filter data based on the 'Kon' column
amoris_men <- amoris[amoris$Kon == 1, ]
amoris_women <- amoris[amoris$Kon == 2, ]


## check pearson correlation
# Calculate Pearson correlations between Age and biomarkers
correlations_pearson_m <- cor(amoris_men$Alder, amoris_men[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "pearson", use = "complete.obs")
correlations_spearman_m <- cor(amoris_men$Alder, amoris_men[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "spearman", use = "complete.obs")
correlations_pearson_f <- cor(amoris_women$Alder, amoris_women[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "pearson", use = "complete.obs")
correlations_spearman_f <- cor(amoris_women$Alder, amoris_women[, c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_P", "Fe_maet", "S_Urea", "S_LD", "S_Ca", "fS_TIBC", "fS_Jaern", "S_Alp", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")], method = "spearman", use = "complete.obs")
print(correlations_pearson_m)#Fe_maet, fS_TIBIC, fS_Jaern, S_Alp are under abs(0.1)
print(correlations_spearman_m)#Fe_maet, fS_TIBIC, fS_Jaern, S_Alp are under abs(0.1)
print(correlations_pearson_f)#S_P,Fe_maet,S_Ca, fS_Jarn, are under abs(0.1)
print(correlations_spearman_f)#S_P,Fe_maet, S_Ca, fS_Jarn, are under abs(0.1)

#print(abs(correlations_pearson_m) > 0.1)

#pearson_m > 0.1: S_Alb, S_Krea, TC, TG, S_K, S_P, S_Urea, S_LD, S_Ca, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#pearson_m > 0.15: S_Alb, S_Krea, TC, S_P, S_Urea, S_Ca, S_Hapt, fS_Gluk
#pearson_m > 0.2: S_Alb, TC, S_P, S_Urea, S_Ca, fS_Gluk

#spearman_m > 0.1: S_Alb, S_Krea, TC, TG, S_K, S_P, S_Urea, S_LD, S_Ca, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#spearman_m > 0.15: S_Alb, TC, TG, S_P, S_Urea, S_Ca, S_Hapt, fS_Gluk
#spearman_m > 0.2: S_Alb, TC, S_P, S_Urea, S_Ca, S_Hapt, fS_Gluk

#pearson_f >0.1: S_Alb, S_Krea, TC, TG, S_K, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#pearson_f > 0.15: S_Alb, S_Krea, TC, TG, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#pearson_f > 0.2: S_Alb, TC, TG, S_Urea, S_LD, S_Alp, S_Urat, fS_Gluk

#spearman_f >0.1: S_Alb, S_Krea, TC, TG, S_K, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#spearman_f > 0.15: S_Alb, S_Krea, TC, TG, S_Urea, S_LD, fS_TIBC, S_Alp, S_Urat, S_FAMN, S_Hapt, fS_Gluk
#spearman_f > 0.2: S_Alb, TC, TG, S_Urea, S_LD, S_Alp, S_Urat, fS_Gluk
```

```{r}
rm(amoris_men)
rm(amoris_women)
rm(correlations_pearson)
rm(correlations_pearson_m)
rm(correlations_pearson_f)
rm(correlations_spearman)
rm(correlations_spearman_m)
rm(correlations_spearman_f)
```


Now we know which biomarkers to use in the BA analysis
```{r}
#We will now redefine our biomarkers to this combination
#over both 0.10 separately for both men and women, 25-75:
#biomarkers = c("S_Alb", "S_Krea", "TC", "TG", "S_Urea", "S_LD", "S_Ca", "S_FAMN", "S_Hapt", "fS_Gluk")

##new
#biomarkers = c("S_Alb", "S_Krea", "TC", "TG", "S_K", "S_Urea", "S_LD", "S_Ca", "S_Urat", "S_FAMN", "S_Hapt", "fS_Gluk")
```

Now let's check scatter plot against age
```{r}
library(ggplot2)
library(tidyr)

# Convert the data to long format
long_data <- gather(AMORIS, key = "biomarker", value = "value", biomarkers)
long_data$age <- rep(AMORIS$Alder, length(biomarkers))  # Repeat age for each biomarker

# Plot scatter plots
ggplot(long_data, aes(x = Alder, y = value)) +
  geom_point(alpha = 0.1, size = 0.1) +
  facet_wrap(~ biomarker, scales = "free") +
  theme_minimal() +
  labs(x = "Age", y = "Value", title = "Scatter Plots of Biomarkers Against Age")

rm(long_data)
```

Scatter with 1% of the data randomly
```{r}
library(dplyr)
library(tidyr)
library(ggplot2)

# Assuming 'biomarkers' is a vector of the biomarker column names
# For example: biomarkers <- c("biomarker1", "biomarker2", ...)

# Calculate 10% of the number of patients
sample_size <- floor(0.01 * nrow(AMORIS))

# Randomly sample patient rows
sampled_AMORIS <- AMORIS[sample(nrow(AMORIS), size = sample_size), ]

# Select only relevant columns (age and biomarkers) and convert to long format
long_data <- sampled_AMORIS %>%
  select(Alder, all_of(biomarkers)) %>%
  pivot_longer(
    cols = -Alder, # Exclude the age column from reshaping
    names_to = "biomarker",
    values_to = "value"
  )

# Plot scatter plots
ggplot(long_data, aes(x = Alder, y = value)) +
  geom_point(alpha = 0.1, size = 0.1) +
  facet_wrap(~ biomarker, scales = "free") +
  theme_minimal() +
  labs(x = "Age", y = "Value", title = "Scatter Plots of Biomarkers Against Age")

```
One on one view
```{r}
# Assuming 'biomarkers' is a vector of the biomarker column names
# For example: biomarkers <- c("biomarker1", "biomarker2", ...)

# Calculate 10% of the number of patients
sample_size <- floor(0.10 * nrow(AMORIS))

# Randomly sample patient rows
sampled_AMORIS <- AMORIS[sample(nrow(AMORIS), size = sample_size), ]

# Loop through each biomarker and create a scatter plot
for (biomarker in biomarkers) {
  long_data <- sampled_AMORIS %>%
    select(Alder, all_of(biomarker)) %>%
    pivot_longer(
      cols = -Alder, # Exclude the age column from reshaping
      names_to = "biomarker",
      values_to = "value"
    )

  # Create scatter plot for each biomarker
  plot <- ggplot(long_data, aes(x = Alder, y = value)) +
    geom_point(alpha = 0.1, size = 0.1) +
    theme_minimal() +
    labs(x = "Age", y = biomarker, title = paste("Scatter Plot of", biomarker, "Against Age"))
  
  print(plot) # Print the plot
}

rm(sampled_AMORIS)
rm(long_data)
rm(biomarker)
rm(number_of_biomarkers)
rm(sample_size)
```



Preparing AMORIS for BioAge
```{r}
colnames(AMORIS)[colnames(AMORIS) == "Kon"] <- "gender"
colnames(AMORIS)[colnames(AMORIS) == "Alder"] <- "age"
colnames(AMORIS)[colnames(AMORIS) == "Id"] <- "sampleID"
```

```{r}

AMORIS <- AMORIS %>%
  mutate(status = ifelse(is.na(DODSDATn), 0, 1)) %>%
  mutate(time = as.numeric(difftime(lastDate, firstDate, units = "weeks")) / 4.33)

negative_counts <- sum(AMORIS$time < 0)
print(paste("Number of negative survival times:", negative_counts))

# Filtering out negative times
AMORIS <- AMORIS %>% filter(time >= 0)
```
Figuring out the size of our training set
```{r}
### Training on a amoris subset the same size as nhanes 3, in the range 25-75
## Figure out number of nhanes participants and participants in range 25-50
load("W:/C6_AmorisGuest/Sara Hagg/Karolina Gustavsson/BioAge-master/BioAge/data/NHANES3.rda")
# Checking the size of 25-75, and 30-75 NHANES people  
# lets check the size!
count_df <- NHANES3[NHANES3$age >= 30 & NHANES3$age <= 75, ]
# Count the number of rows in the filtered dataframe
num_rows <- nrow(count_df)
num_rows
# okay so number of people in 25-75 in nhanes3 is 14377
# okay so number of people in 30-75 in nhanes3 is 12535, which is what they trained on, we also want this amount
```

```{r}
rm(NHANES3)
rm(num_rows)
rm(count_df)
rm(negative_counts)
```



Considering the distribution of dementia cases
```{r}
# Distribution for unprocessed data
amoris$Dementia_binary <- ifelse(amoris$Dementia == "", "No Dementia", "Dementia")
ggplot(amoris, aes(x = Alder, fill = Dementia_binary)) + 
  geom_histogram(position = "identity", alpha = 0.5, bins = 30) +
  scale_fill_manual(values = c("No Dementia" = "red", "Dementia" = "blue")) +
  theme_minimal() +
  labs(title = "Age Distribution by Dementia Status",
       x = "Age",
       y = "Count")
#okay so this confirms that it's ok to use the cutoff 75 years old
filtered_data_ <- subset(amoris, Alder >= 25 & Alder <= 75)
table_counts <- table(filtered_data_$Dementia_binary)
print(table_counts)
# Dementia 21457 for 25-75
# No dementia 221434 for 25-75

AMORIS$Dementia_binary_AM <- ifelse(AMORIS$Dementia == "", "No Dementia", "Dementia")
ggplot(AMORIS, aes(x = age, fill = Dementia_binary_AM)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 30) +
  scale_fill_manual(values = c("No Dementia" = "red", "Dementia" = "blue")) +
  theme_minimal() +
  labs(title = "Age Distribution by Dementia Status",
       x = "Age",
       y = "Count")
#okay so this confirms that it's ok to use the cutoff 75 years old
filtered_data_ <- subset(AMORIS, age >= 25 & age <= 75)
table_counts <- table(filtered_data_$Dementia_binary_AM)
print(table_counts)
# Dementia 21457 for 25-75
# No dementia 221434 for 25-75
```
```{r}
rm(table_counts)
rm(filtered_data_)
```

Now let's plot per subtype, first we need to categorise
```{r}
alzheimers_codes <- c("304", "290", "290A", "F00", "305", "290B", "G30", "311A", "G308", "G309", "G301", "F009", "F002", "G300", "F001", "F000")
vascular_codes <- c("293", "290E", "F01", "293,1", "F019", "F012", "F018", "F013", "F011", "F010", "F0183")
other_codes <- c("306", "290X", "F02", "290W", "F03", "294B", "G311", "311B", "G318A", "311C", "F051", "311X", "F039", "F020", "F024", "F028", "290", "290,1", "F022", "F021", "F023", "F03-P") 
unclear_type <- c("293,9","293,4", "293,3", "293,5", "F0010", "F0013", "F0180")

amoris$ADRD <- ifelse(amoris$Dementia == "", 0, 1)
amoris$Alzheimers <- ifelse(amoris$Dementia %in% alzheimers_codes, 1, 0)
amoris$Vascular <- ifelse(amoris$Dementia %in% vascular_codes, 1, 0)
amoris$Other <- ifelse(amoris$Dementia %in% other_codes, 1, 0)
amoris$Unclear <- ifelse(amoris$Dementia %in% unclear_type, 1, 0)

sum(amoris$ADRD) #23955
sum(amoris$Alzheimers) #9086
sum(amoris$Vascular) #3630
sum(amoris$Other) #11199
sum(amoris$Unclear) #24

AMORIS$ADRD <- ifelse(AMORIS$Dementia == "", 0, 1)
AMORIS$Alzheimers <- ifelse(AMORIS$Dementia %in% alzheimers_codes, 1, 0)
AMORIS$Vascular <- ifelse(AMORIS$Dementia %in% vascular_codes, 1, 0)
AMORIS$Other <- ifelse(AMORIS$Dementia %in% other_codes, 1, 0)
AMORIS$Unclear <- ifelse(AMORIS$Dementia %in% unclear_type, 1, 0)

sum(AMORIS$ADRD) #22987
sum(AMORIS$Alzheimers) #8889
sum(AMORIS$Vascular) #3506
sum(AMORIS$Other) #10581
sum(AMORIS$Unclear) #11
```

Now we plot
```{r}
library(dplyr)

amoris <- amoris %>%
  mutate(ADRD_Category = case_when(
    Vascular == 1 ~ 'Vascular',
    Alzheimers == 1 ~ 'Alzheimers',
    Other == 1 ~ 'Other',
    Unclear == 1 ~ 'Unknown',
    ADRD == 0 ~ 'No_ADRD',
    TRUE ~ 'Uncategorized'  # Catch-all for any other cases
  ))

# Create a single histogram with overlapping categories
ggplot(amoris, aes(x = Alder, fill = ADRD_Category)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 30) +
  scale_fill_manual(values = c("Vascular" = "red", 
                               "Alzheimers" = "blue",
                               "Other" = "green",
                               "Unknown" = "yellow",
                               "No_ADRD" = "grey",
                               "Uncategorized" = "black")) +
  theme_minimal() +
  labs(title = "Age Distribution by ADRD Category", x = "Age", y = "Count")

AMORIS <- AMORIS %>%
  mutate(ADRD_Category = case_when(
    Vascular == 1 ~ 'Vascular',
    Alzheimers == 1 ~ 'Alzheimers',
    Other == 1 ~ 'Other',
    Unclear == 1 ~ 'Unknown',
    ADRD == 0 ~ 'No_ADRD',
    TRUE ~ 'Uncategorized'  # Catch-all for any other cases
  ))

# Create a single histogram with overlapping categories
ggplot(AMORIS, aes(x = age, fill = ADRD_Category)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 30) +
  scale_fill_manual(values = c("Vascular" = "red", 
                               "Alzheimers" = "blue",
                               "Other" = "green",
                               "Unknown" = "yellow",
                               "No_ADRD" = "grey",
                               "Uncategorized" = "black")) +
  theme_minimal() +
  labs(title = "Age Distribution by ADRD Category", x = "Age", y = "Count")

```


```{r}
library(ggplot2)
library(dplyr)

# Create the ADRD_Category variable
amoris <- amoris %>%
  mutate(ADRD_Category = case_when(
    Vascular == 1 ~ 'Vascular',
    Alzheimers == 1 ~ 'Alzheimers',
    Other == 1 ~ 'Other',
    TRUE ~ NA_character_  # Assign NA for cases not covered above
  ))

# Filter out NA in ADRD_Category to include only specific ADRD subtypes
amoris_adrd_subtypes <- filter(amoris, !is.na(ADRD_Category))

# Create a histogram for specific ADRD subtypes
ggplot(amoris_adrd_subtypes, aes(x = Alder, fill = ADRD_Category)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 30) +
  scale_fill_manual(values = c("Vascular" = "red", 
                               "Alzheimers" = "blue",
                               "Other" = "green")) +
  theme_minimal() +
  labs(title = "Age Distribution by Specific ADRD Subtypes", x = "Age", y = "Count")

AMORIS <- AMORIS %>%
  mutate(ADRD_Category = case_when(
    Vascular == 1 ~ 'Vascular',
    Alzheimers == 1 ~ 'Alzheimers',
    Other == 1 ~ 'Other',
    TRUE ~ NA_character_  # Assign NA for cases not covered above
  ))

# Filter out NA in ADRD_Category to include only specific ADRD subtypes
AMORIS_adrd_subtypes <- filter(AMORIS, !is.na(ADRD_Category))

# Create a histogram for specific ADRD subtypes
ggplot(AMORIS_adrd_subtypes, aes(x = age, fill = ADRD_Category)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 30) +
  scale_fill_manual(values = c("Vascular" = "red", 
                               "Alzheimers" = "blue",
                               "Other" = "green")) +
  theme_minimal() +
  labs(title = "Age Distribution by Specific ADRD Subtypes", x = "Age", y = "Count")


```


Let's look at the other variables:
- Economic variables
- BMI
- Education level 

Let's start with economic variables: DispInk, DispInkFam, DispInkKE, ArbInk (before dementia)
```{r}
library(ggplot2)
library(tidyr)

# Reshaping the data to long format for easy plotting
AMORIS_long <- AMORIS %>%
  select(DispInk, DispInkKE, DispInkFam, ArbInk) %>%
  pivot_longer(cols = everything(), names_to = "IncomeType", values_to = "Income")

# Density plot of raw data
ggplot(AMORIS_long, aes(x = Income, fill = IncomeType)) +
  geom_density(alpha = 0.5,na.rm = TRUE) +
  labs(title = "Density Plot of Various Income Types", x = "Income", y = "Density") +
  theme_minimal()

print("done")
```
Okay so we have some very extreme outliers, we can not use average/STD, we have to use median/MAD
```{r}
library(ggplot2)
library(dplyr)
library(tidyr)

# Assuming these are your income variables
income_vars <- c("DispInk", "DispInkKE", "DispInkFam", "ArbInk")

# Reshape and then filter outliers based on median and MAD
AMORIS_filtered <- AMORIS %>%
  select(all_of(income_vars)) %>%
  pivot_longer(cols = everything(), names_to = "IncomeType", values_to = "Income") %>%
  group_by(IncomeType) %>%
  mutate(
    median_income = median(Income, na.rm = TRUE),
    mad_income = mad(Income, constant = 1, na.rm = TRUE)  # MAD as a robust measure
  ) %>%
  filter(between(Income, median_income - 5 * mad_income, median_income + 5 * mad_income)) %>%
  ungroup() %>%
  select(-median_income, -mad_income)

# Density plot of data with outliers removed
ggplot(AMORIS_filtered, aes(x = Income, fill = IncomeType)) +
  geom_density(alpha = 0.5, na.rm = TRUE) +
  labs(title = "Density Plot of Various Income Types (Outliers Removed)", x = "Income", y = "Density") +
  theme_minimal()

```

BMI
```{r}
ggplot(AMORIS, aes(x = BMI)) +
  geom_density(fill = "blue", alpha = 0.5, na.rm = TRUE) +
  labs(title = "Density Plot of BMI", x = "BMI", y = "Density") +
  theme_minimal()
```
Distribution per age
```{r}
ggplot(AMORIS, aes(x = Alder)) + 
  geom_density(fill = "blue", alpha = 0.5) + 
  labs(title = "Density Plot of Ages", x = "Age", y = "Density") +
  theme_minimal()
```


Education level
```{r}
library(ggplot2)

# Convert the education level to a factor, treating empty strings as a separate level
AMORIS$UtbNiva <- factor(AMORIS$UtbNiva, levels = unique(AMORIS$UtbNiva))

# Create a bar plot for education levels
ggplot(AMORIS, aes(x = UtbNiva)) +
  geom_bar(fill = "blue", alpha = 0.7) +
  labs(title = "Histogram of Education Level", x = "Education Level", y = "Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))  # Adjust the x-axis text angle for better readability

```
Summary statistics 
```{r}
# Identify continuous (numeric) variables
continuous_vars <- sapply(AMORIS, is.numeric) & !sapply(AMORIS, is.factor)

# Custom summary function including standard deviation
custom_summary <- function(x) {
  c(Mean = mean(x, na.rm = TRUE),
    SD = sd(x, na.rm = TRUE),
    Median = median(x, na.rm = TRUE),
    '1st Qu' = quantile(x, probs = 0.25, na.rm = TRUE),
    '3rd Qu' = quantile(x, probs = 0.75, na.rm = TRUE),
    Min = min(x, na.rm = TRUE),
    Max = max(x, na.rm = TRUE))
}

# Apply custom summary function to each continuous variable
summary_stats <- lapply(AMORIS[continuous_vars], custom_summary)

# Convert the list to a dataframe for a cleaner presentation
summary_df <- do.call(rbind, summary_stats)

# View the summary statistics
print(summary_df)
```


```{r}
rm(AMORIS_long)
rm(AMORIS_filtered)
rm(income_vars)
rm(agevar)
rm(agevar_ap)
rm(label)
rm(label_ap)
```

Tables
```{r}
# Convert Dementia into binary and create a new column ADRD
AMORIS$`Total ADRD` <- ifelse(AMORIS$Dementia == "", 0, 1)

# Create binary columns for dementia subtypes
AMORIS$Alzheimers <- ifelse(AMORIS$Dementia %in% alzheimers_codes, 1, 0)
AMORIS$`Vascular Dementia` <- ifelse(AMORIS$Dementia %in% vascular_codes, 1, 0)
AMORIS$`Other Dementias` <- ifelse(AMORIS$Dementia %in% other_codes, 1, 0)
AMORIS$`Unclear Type` <- ifelse(AMORIS$Dementia %in% unclear_type, 1, 0)

# Compute statistics for continuous variables
compute_stats_cont <- function(var) {
  total_count <- sum(!is.na(AMORIS[[var]]))
  men_count <- sum(!is.na(filter(AMORIS, gender == "1")[[var]]))
  women_count <- sum(!is.na(filter(AMORIS, gender == "2")[[var]]))
  
  data.frame(
    Characteristics = var,
    Total = paste0(round(mean(AMORIS[[var]], na.rm = TRUE), 2), " ± ", round(sd(AMORIS[[var]], na.rm = TRUE), 2), " (N=", total_count, ")"),
    Men = paste0(round(mean(filter(AMORIS, gender == "1")[[var]], na.rm = TRUE), 2), " ± ", round(sd(filter(AMORIS, gender == "1")[[var]], na.rm = TRUE), 2), " (N=", men_count, ")"),
    Women = paste0(round(mean(filter(AMORIS, gender == "2")[[var]], na.rm = TRUE), 2), " ± ", round(sd(filter(AMORIS, gender == "2")[[var]], na.rm = TRUE), 2), " (N=", women_count, ")")
  )
}

# Compute statistics for categorical variables
compute_stats_cat <- function(var) {
  
  total_count <- sum(AMORIS[[var]], na.rm = TRUE)
  men_count <- sum(filter(AMORIS, gender == "1")[[var]], na.rm = TRUE)
  women_count <- sum(filter(AMORIS, gender == "2")[[var]], na.rm = TRUE)
  
  data.frame(
    Characteristics = var,
    Total = paste0(round(mean(AMORIS[[var]], na.rm = TRUE) * 100, 2), "% (N=", total_count, ")"),
    Men = paste0(round(mean(filter(AMORIS, gender == "1")[[var]], na.rm = TRUE) * 100, 2), "% (N=", men_count, ")"),
    Women = paste0(round(mean(filter(AMORIS, gender == "2")[[var]], na.rm = TRUE) * 100, 2), "% (N=", women_count, ")")
  )
}

# Compute p-values
compute_p_values <- function(var) {
  if (var %in% c("status", "ADRD", "Alzheimers", "Vascular Dementia", "Other Dementias", "Unclear Type")) {
    tbl <- table(AMORIS$gender, AMORIS[[var]])
    p_val <- tryCatch(chisq.test(tbl)$p.value, error = function(e) NA)
  } else {
    p_val <- tryCatch(t.test(AMORIS[[var]] ~ AMORIS$gender)$p.value, error = function(e) NA)
  }
  return(p_val)
}

# List of variables
vars <- c("age", 
          "S_Alb", "S_Krea", "TG", "TC", "S_Urea", "S_LD", "S_Ca", 
          "S_FAMN", "S_Hapt", "fS_Gluk")

# Compute stats
stats_df <- bind_rows(lapply(vars, function(v) {
  if (v %in% c("status", "Total ADRD", "Alzheimers", "Vascular Dementia", "Other Dementias", "Unclear Type")) {
    compute_stats_cat(v)
  } else {
    compute_stats_cont(v)
  }
}))

# Compute p-values and add to the dataframe
stats_df$`p value` <- sapply(vars, compute_p_values)
stats_df$`p value` <- ifelse(stats_df$`p value` < 0.001, "<0.001", round(stats_df$`p value`, 3))

# Make sure that 'Characteristics' column is a character and not factor to avoid issues
stats_df$Characteristics <- as.character(stats_df$Characteristics)

# Create the gt table
gt_table <- gt(stats_df) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      rows = stats_df$Characteristics %in% c("BA measures", "Biomarkers", "ADRD"),
      columns = c("Characteristics")
    )
  ) %>%
  tab_options(row.striping.include_table_body = FALSE)

print(gt_table)
```

Checking out SveDem data
```{r}
SveDem <- read_sas("W:/C6_AmorisGuest/Sara Hagg/Karolina Gustavsson/Data/SveDem.sas7bdat")
SveDem_format <- read_sas("W:/C6_AmorisGuest/Sara Hagg/Karolina Gustavsson/Data/SveDem_fomats.sas7bdat")
```

Checking out the Svedem unique codes
```{r}
number_of_unique_strings <- length(unique(SveDem$Dementia))
print(number_of_unique_strings) # 26 unique codes, it says 27 but one is the empty chr str
unique_strings <- unique(SveDem$Dementia)
print(unique_strings)

Svedem_ADRD <- c("G308", "F019",  "G301",  "F013",  "F023",  "F039",  "F002",  "F012",  "G300",  "F011",  "311X",  "F001", "G318A", "G309",  "F020",  "F018",  "F009",  "F051",  "F028",  "G311",  "F000",  "F03-P", "F010",  "290B",  "F024", "290A")
Svedem_Alzheimers <- c("304",	"290",	"290A",	"F00", "305",	"290B",	"G30", "311A")
Svedem_Vascular <- c("293",	"290E",	"F01", "293,1")
Svedem_Other <- c("306",	"290X",	"F02", "290W",	"F03", "294B",	"G311", "311B",	"G318A", "311C",	"F051", "311X")

#more categories
```
```{r}
SveDem$S_ADRD <- ifelse(SveDem$Dementia %in% Svedem_ADRD, 1, 0)
SveDem$S_Alzheimers <- ifelse(SveDem$Dementia %in% Svedem_Alzheimers, 1, 0)
SveDem$S_Vascular <- ifelse(SveDem$Dementia %in% Svedem_Vascular, 1, 0)
SveDem$S_Other <- ifelse(SveDem$Dementia %in% Svedem_Other, 1, 0)
```

Do we have entires in SveDem where Dementia diagnosis is not added to amoris?
```{r}
# Step 1: Find Matching IDs
matching_ids <- intersect(amoris$Id, SveDem$Id)

# Step 2: Merge Data based on Matching IDs
merged_data <- merge(amoris[amoris$Id %in% matching_ids, c("Id", "Dementia")], 
                     SveDem[SveDem$Id %in% matching_ids, c("Id", "Dementia")], 
                     by.x = "Id", by.y = "Id")

# Step 3: Compare 'Dementia' Columns
# Assuming the Dementia columns are named Dementia.x and Dementia.y after merging
mismatches <- merged_data$Dementia.x != merged_data$Dementia.y

# Step 4: Count Mismatches
num_mismatches <- sum(mismatches)
print(paste("Number of unmatched Dementia entries:", num_mismatches))

# Okay so we already have them all, SveDem has more accurate diagnosis so we should run separate SveDem Diagnosis 
```
Checking if the SveDem codes exist in the non-matching IDs in AMORIS ie are they unique to Svedem
```{r}
# Step 1: Identify Non-Matching IDs in amoris
non_matching_ids_amoris <- setdiff(amoris$Id, SveDem$Id)

# Step 2: Check if Svedem_ADRD Codes are in Non-Matching IDs
codes_as_ids <- as.numeric(Svedem_ADRD)  # Convert codes to numeric if necessary
matching_codes_in_amoris <- intersect(codes_as_ids, non_matching_ids_amoris)

# Print result
if (length(matching_codes_in_amoris) > 0) {
  print(paste("Matching codes found in the non-matching IDs of amoris:", 
              paste(matching_codes_in_amoris, collapse = ", ")))
} else {
  print("None of the specified codes are found in the non-matching IDs of amoris.")
}

#It seems these codes are specific to SveDem
```


Optional
```{r}
rm(amoris)
rm(AMORIS)
rm(biomarkers)
```