# check your R version
R.Version()

# or just
#R.version.string

# Installing (ONLY the 1st time)
utils::install.packages('here')

# OR (same)
install.packages('here')

# Load pckgs for this R session
library(fs)        # file/directory interactions
library(here)      # tools find your project's files, based on working directory
library(janitor)   # tools for examining and cleaning data
library(skimr)     # tools for summary statistics 
library(dplyr)     # {tidyverse} tools for manipulating and summarising tidy data 
library(forcats)   # {tidyverse} tool for handling factors
library(ggplot2)   # {tidyverse} tools for plotting
library(ggridges)  # alternative to plot density functions 

# Opening Help page on package/function
?here

??here

# [NOT REPRODUCIBLE] hard coding your file path  -----------------------

# File path on Mac:
dataset <- readr::read_csv("/Users/testuser/R4biostats/input_data/dataset.csv")
# Same file path on Windows:
dataset <- readr::read_csv("C:\Users\testuser\R4biostats\input_data\dataset.csv")

?read.csv # to learn about function and arguments 

autism_data_url <- read.csv(
  file = "https://raw.githubusercontent.com/Sydney-Informatics-Hub/lessonbmc/gh-pages/_episodes_rmd/data/autism_data.csv", 
  header = TRUE, # 1st line is the name of the variables
  sep = ",", # which is the field separator character.
  na.strings = c("?") # specific values R should interpret as NA
)

# Check my working directory location
# here::here()

# Use `here` in specifying all the subfolders AFTER the working directory 
autism_data_file <- read.csv(
  file = here("practice", "data_input", "01_datasets", "autism_data.csv"), 
  header = TRUE, # 1st line is the name of the variables
  sep = ",", # which is the field separator character.
  na.strings = c("?"),# specific values R should interpret as NA
  row.names = NULL) 

View(autism_data_file)

# What data type is this data?
class(autism_data_file)

# What variables are included in this dataset?
base::colnames(autism_data_file)

autism_data <- janitor::clean_names(autism_data_file, 
                                     case = "none") 
# check change
colnames(autism_data)

dim(autism_data)

autism_data$id
autism_data$A1_Score
autism_data$gender
autism_data$autism

# rename dataframe 
autism_pids <- autism_data

# create a new column 
autism_pids$pids <- paste("PatientID_" , autism_data$id, sep = "")

# check change in df structure
base::colnames(autism_pids)

dim(autism_data)

dim(autism_pids)

# what do I have in the environment? 
ls() 

# remove all EXCEPT for "autism_pids" 
rm("autism_data", "autism_data_file", "autism_data_url" ) 

# With the `$` sign I extract a variable (column name)
head(autism_pids$id) 

head(autism_pids$pids)

head(autism_pids$A1_Score)

head(autism_pids$ethnicity)

# Indexing to pick `[ , #col]`  
head(autism_pids[ ,1] )# empty rows means all 

head(autism_pids[ ,23])

head(autism_pids[ ,2])

head(autism_pids[ ,14])

# Indexing to pick `[#row, ]`  
head(autism_pids[1 , ] ) # empty cols means all 

head(autism_pids[50,])

head(autism_pids[25:26 ,])

# Indexing to pick `[#row, #col]`  
autism_pids[1:3,1]

autism_pids[1:3,23]

autism_pids[1:3,2]

autism_pids[1:3,14]

# What are the data types of the variables? ---------------------------------
str(autism_pids) # integer and character

# What values can the variables take? ---------------------------------
summary(autism_pids$pids)

length(unique(autism_pids$pids)) # N unique values

sum(is.na(autism_pids$pids)) # N missing values

summary(autism_pids$ethnicity)

length(unique(autism_pids$ethnicity)) # N unique values

sum(is.na(autism_pids$ethnicity)) # N missing values

autism_pids %>% 
  skimr::skim(pids, ethnicity) %>%
  dplyr::select(#skim_variable, 
                skim_type, 
                complete_rate,
                n_missing, 
                character.n_unique)

autism_pids %>% 
  skimr::skim() 

#### char 2 factor -------------------------------------------------------------
# Say I want to treat some variables as factors
autism_pids$gender <- as.factor(autism_pids$gender)
autism_pids$ethnicity <- as.factor(autism_pids$ethnicity)
autism_pids$contry_of_res <- as.factor(autism_pids$contry_of_res)
autism_pids$relation <- as.factor(autism_pids$relation)

# check 
class(autism_pids$gender)

class(autism_pids$ethnicity)

class(autism_pids$contry_of_res)

class(autism_pids$relation)

autism_pids_temp <- autism_pids # copy df for test 

to_factor <- c("gender", "ethnicity", "contry_of_res", "relation") # vector of col names 
autism_pids_temp[ ,to_factor] <-  lapply(X =  autism_pids[ ,to_factor], FUN = as.factor)

# check 
class(autism_pids_temp$gender)

class(autism_pids_temp$ethnicity)

class(autism_pids_temp$contry_of_res)

class(autism_pids_temp$relation)

# now I have Variable type: factor

levels(autism_pids$ethnicity)

table(autism_pids$ethnicity,useNA = "ifany")

janitor::tabyl(autism_pids$ethnicity) %>% 
  adorn_totals() %>% 
  adorn_pct_formatting()

which(is.na(autism_pids$ethnicity)) # indices of TRUE elements in vector

which(is.na(autism_pids$relation))  # indices of TRUE elements in vector

# observe a subset of some columns 
autism_subset <- autism_pids [1:5, c("gender","jaundice", "autism","age_desc",
                                     "Class_ASD","pids")]
# View(autism_subset)

# recode "age_desc" as LOGICAL new var "age_desc_log"
autism_pids$age_desc_log <- ifelse(autism_pids$age_desc == "18 and more", TRUE, FALSE )
class(autism_pids$age_desc)

class(autism_pids$age_desc_log)

autism_pids$autism_dummy <- ifelse(autism_pids$autism == 'yes', 1, 0)
class(autism_pids$autism)

class(autism_pids$autism_dummy)

colnames(autism_pids)

head(autism_pids)   #return fist 6 obs
tail(autism_pids)   #return last 6 obs

head(autism_pids, n = 2) #return fist 2 obs

tail(autism_pids, n = 2) #return last 2 obs

# run...
sum(is.na(autism_pids$age)) 

# or 
skimr::n_missing(autism_pids$age)

missing_age_subset <- autism_pids[is.na(autism_pids$age), 
                                  c("pids", "age", "autism_dummy") ]
missing_age_subset

missing_age_subset2 <- autism_pids[which(is.na(autism_pids$age)), 
                                   c("pids", "age", "autism_dummy")] 
missing_age_subset2

# arguments allow me to specify rows and cols 
missing_age_subset3 <- subset(x = autism_pids, 
                              subset = is.na(autism_pids$age), # 1 logical condition
                              select = c("pids", "age", "autism_dummy") # which cols
                              ) 
missing_age_subset3

# Creates a SUBSET based on MORE conditions (`age` and `ethnicity`)
twocond_base_subset <- subset(x = autism_pids, 
                       # 2 logical conditions      
                       subset = age < 25 & contry_of_res == "Brazil", 
                       # pick a few cols 
                       select = c("pids", "age", "contry_of_res",
                                  "autism_dummy")) 

twocond_base_subset

## here the filtering (rows) and selecting (columns) is done in sequence
twocond_dplyr_subset <- autism_pids %>% 
  dplyr::filter(age < 25 & contry_of_res == "Brazil") %>%  # which rows
  dplyr::select (pids, age, contry_of_res, autism_dummy)   # which cols

twocond_dplyr_subset

summary(autism_pids$A1_Score)     # min, max quartiles, mean, median

summary(autism_pids$ethnicity)    # counts of levels' frequency (included NA!)

summary(autism_pids$age_desc_log) # counts of TRUE 

table(autism_pids$gender)

table(autism_pids$age) # automatically drops missing...

table(autism_pids$age, useNA = "ifany") #...unless specified

table(autism_pids$gender, autism_pids$age_inputed)

table(autism_pids$ethnicity, autism_pids$autism_dummy)

# by(data$column, data$grouping_column, mean)
by(data = autism_pids$age_inputed, INDICES = autism_pids$gender, FUN = mean)

autism_pids %>% 
  dplyr::group_by(gender) %>% 
  dplyr::summarise(mean(age_inputed))  # returns a dataframe!

autism_pids %>% 
  dplyr::group_by(gender) %>% 
  dplyr::summarise(mean_age = mean(age_inputed),  
                   N_obs = n(), 
                   N_with_autism = sum(autism_dummy == 1)
  ) 

f_calc_mode  <- function(x) { 
  # `unique` returns a vector of unique values 
  uni_x <- unique(x)  
  # `match` returns the index positions of 1st vector against 2nd vector
  match_x <- match(x, uni_x)
  # `tabulate` count the occurrences of integer values in a vector.
  tab_x  <-  tabulate(match_x) 
  # returns element of uni_x that corresponds to max occurrences
  uni_x[tab_x == max(tab_x)]
}

f_calc_mode(autism_pids$age)

f_calc_mode(autism_pids$age_inputed)

f_calc_mode(autism_pids$ethnicity)

var(autism_pids$age, na.rm = TRUE)

var(autism_pids$age_inputed)

sd(autism_pids$age, na.rm = TRUE)

sd(autism_pids$age_inputed)

two_col_palette <-  c("#9b2339", "#005ca1")

contrast_cols_palette <- c("#E7B800","#239b85", "#85239b", "#9b8523","#23399b",
                "#d8e600", "#0084e6", "#399B23", "#e60066",
                "#00d8e6", "#e68000")

ggplot(data = autism_pids, mapping = aes(x=age_inputed)) + 
  geom_histogram() + 
  theme_bw()

autism_pids %>% 
  ggplot(aes(x = age_inputed )) + 
  geom_histogram() + 
  theme_bw()

autism_pids %>% 
  ggplot(aes(x = age_inputed )) + 
  # specify to avoid warning if we fail to specify the number of bins 
  geom_histogram(bins=40) + 
  theme_bw()

autism_pids %>% 
  ggplot(aes(x = age_inputed )) + 
  geom_histogram(bins=40) + 
  # add mean vertical line
  geom_vline(xintercept = mean(autism_pids$age_inputed),
             na.rm = FALSE,
             lwd=1,
             color="#9b2339") +
  # add annotations with the mean value
  annotate("text",                        
           x = mean(autism_pids$age_inputed) * 1.2, # coordinates for positioning
           y = mean(autism_pids$age_inputed) * 2.5,
           label = paste("Mean =", round(mean(autism_pids$age_inputed), digits = 2)),
           col = "#9b2339",
           size = 4)+
  # add also sd +1 and -1 
  geom_vline(aes(xintercept = mean(autism_pids$age_inputed) + sd(autism_pids$age_inputed)), 
             color = "#000000", size = 1, linetype = "dashed") +
  geom_vline(aes(xintercept = mean(autism_pids$age_inputed) - sd(autism_pids$age_inputed)), 
             color = "#000000", size = 1, linetype = "dashed") +
  theme_bw() 

autism_pids %>% 
  ggplot(aes(x = age_inputed)) +
  geom_density()+
  theme_bw() 

autism_pids %>% 
  ggplot(aes( x=age_inputed)) +
  geom_density(fill="#85239b", color="#4c4c4c", alpha=0.5)+
  theme_bw() 

autism_pids %>% 
  ggplot(aes( x=age_inputed)) +
  geom_density(fill="#85239b", color="#4c4c4c", alpha=0.5)+
  theme_bw() + 
  # increase number of x axis ticks 
  scale_x_continuous(breaks = seq(10, 100,5 ), limits = c(16, 86))

# trying to improve readability 
autism_pids %>% 
  ggplot(mapping = aes(x = age_inputed, fill = gender )) + 
  # bars next to each other with `position = 'dodge'`
  geom_histogram(bins=40, position = 'dodge')  + 
  scale_fill_manual(values = c("#e07689","#57b7ff")) +
  scale_color_manual(values = c("#9b2339","#005ca1")) +
  theme_bw()  

autism_pids %>% 
  ggplot(aes(x = age_inputed, fill = gender )) + 
  geom_histogram(color="#e9ecef", alpha=0.8, position = 'dodge') + 
  theme_bw() + 
  # splitting the gender groups, specifying `ncol` to see one above the other
  facet_wrap(~gender, ncol = 1)  + 
  scale_fill_cyclical(values = c("#9b2339","#005ca1"))

group_stats <- autism_pids %>% 
  dplyr::group_by(gender) %>% 
  dplyr::summarize(mean_age = mean(age_inputed),
                   median_age = median (age_inputed)) 

group_stats

autism_pids %>% 
  ggplot(aes(x = age_inputed, fill = gender)) + 
  # geom_histogram from dataframe 1
  geom_histogram(bins=30,color="#e9ecef", alpha=0.8, position = 'dodge') + 
  facet_wrap(~gender, ncol = 1) + 
  scale_fill_manual(values = c("#9b2339","#005ca1"))  +
  # geom_vline from dataframe 2
  geom_vline(data = group_stats_long, 
             mapping = aes(xintercept = Value, color = Stat),
             lwd=1,
             linetype=1) + 
  scale_color_manual(values = c( "#f0a441" , "#d8cf71")) +
  theme_bw()  

hist_plot <- autism_pids %>% 
  ggplot(aes(x = age_inputed, fill = gender)) + 
  # geom_histogram from dataframe 1
  geom_histogram(bins=30,color="#e9ecef", alpha=0.8, position = 'dodge') + 
  facet_wrap(~gender, ncol = 1) + 
  scale_fill_manual(values = c("#9b2339","#005ca1"))  +
  # geom_vline from dataframe 2
  geom_vline(data = group_stats_long, 
             mapping = aes(xintercept = Value, color = Stat),
             lwd=1.5,
             linetype=6) + 
  scale_color_manual(values = c( "#e68000", "#d8cf71")) +
  # increase number of x axis ticks 
  scale_x_continuous(breaks = seq(10, 100,10 ), limits = c(10,70)) +
  # Additional theme details 
  labs(x = "age brackets", y = "n of individuals",
       color = "Stats",
       title = "Distribution of observations by gender",
       subtitle = "",
       caption = "Source: Thabtah,Fadi (2017) https://doi.org/10.24432/C5F019.") +
  theme(legend.position = "right",
        plot.title = element_text(face = "bold")) 
hist_plot

autism_pids %>% 
  # this takes also `y` = group
  ggplot(aes(x=age_inputed, y = gender, fill = gender)) +
  ggridges::geom_density_ridges() +
  # I can add quantile lines (2 is the median)
  stat_density_ridges(quantile_lines = TRUE, quantiles = c(0.5), alpha = 0.75)+  
  # increase number of x axis ticks 
  scale_x_continuous(breaks = seq(10, 100,10 ), limits = c(16, 86)) + 
  scale_fill_cyclical(values = c("#9b2339","#005ca1")) + 
  theme_bw() 

# Let's take a variable that we recoded as `factor`
class(autism_pids$ethnicity)

#### ... no formatting ---------------------------------- 
autism_pids %>% 
  ggplot(aes(x = ethnicity )) + 
  geom_bar() +   
  theme_bw() 

autism_pids %>% 
  ggplot(aes(x = ethnicity )) + 
  geom_bar(fill = "steelblue") +
  # reference line  
  geom_hline(yintercept=100, color = "#9b2339", size=0.5, ) +
  # labels, title, etc 
  labs(x = "ethnicity", y = "n of individuals",
       color = "Stats",
       title = "Distribution of observations by ethnicity",
       subtitle = "",
       caption = "Autism study")  +
  theme_bw() +
  # specification son axis labels
  theme(axis.text.x = element_text(angle=50, vjust=0.75), 
        axis.text.y = element_text(size=10,face="bold"))

autism_pids %>% 
    # we modify our x like so 
    ggplot(aes(x = forcats::fct_infreq(ethnicity ))) + 
    geom_bar(fill = "steelblue") +
    geom_hline(yintercept=100, color = "#9b2339", size=0.5, ) +
    labs(x = "ethnicity", y = "n of individuals",
         color = "Stats",
         title = "Distribution of observations by ethnicity",
         subtitle = "",
         caption = "Autism study")  +
    # --- wrap long x labels (flipped ) !!!
    #  scale_x_discrete(labels = function(x) stringr::str_wrap(x, width = 10)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle=50, vjust=0.75), 
          axis.text.y = element_text(size=10,face="bold"))  

autism_pids %>%
  ## --- prep the dataframe 
  dplyr::mutate(# Add a factor variable with two levels
    highlight = forcats::fct_other(ethnicity, 
                                   keep = "NA", 
                                   other_level = "All Groups")) %>% 
  ## --- now plot 
  # In `aes mapping` we map color to a variable (`fill = highlight`)
  ggplot(aes(x = forcats::fct_infreq(ethnicity), fill = highlight)) + 
  geom_bar()+
  # Use custom color palettes
  scale_fill_manual(values=c("#0084e6")) +
  # Add a line at a significant level 
  geom_hline(yintercept=100, color = "#9b2339", size=0.5, ) +
  theme_bw() +
  # make some more theme specifications  
  labs(x = "ethnicity", y = "n of individuals",
       color = "Stats",
       title = "Distribution of observations by ethnicity",
       subtitle = "",
       caption = "Autism study")  +
  theme(axis.text.x = element_text(angle=50, vjust=0.75), 
        axis.text.y = element_text(size=10,face="bold"))  +
  theme(legend.position = "none") 

autism_pids %>% 
  ggplot(aes(x = result )) +
  geom_boxplot(alpha=0.5)+
  # switch to vertical orientation
  coord_flip() +
  theme_bw()    

autism_pids %>% 
  ggplot(aes(x = ethnicity,  y= result, fill = ethnicity)) +
  geom_boxplot(alpha=0.5)+
  scale_fill_manual(values =  contrast_cols_palette)   +
  theme_bw()+
  # make x axis labes readable
  theme(axis.text.x = element_text(angle=50, vjust=0.75)) +
  # drop legend and Y-axis title
  theme(legend.position = "none") 

autism_pids %>% 
  ggplot(mapping = aes(y = age_inputed, x = ethnicity, fill = ethnicity)) +
  geom_violin(alpha=0.5) +
  geom_point(position = position_jitter(width = 0.1), size = 0.5)+ 
  scale_fill_manual(values =  contrast_cols_palette)  +
  # make x axis labes readable
  theme(axis.text.x = element_text(angle=50, vjust=0.75)) +
  # drop legend and Y-axis title
  theme(legend.position = "none") 

ggsave (hist_plot, 
        filename = here::here("practice",  "data_output", "hist_plot.png"))

saveRDS (object = autism_pids, 
         file =  here::here("practice",  "data_output", "autism_pids_v2.Rds"))

# to load it later I will use 
readRDS(here::here("practice",  "data_output", "autism_pids_v2.Rds"))