mydata_csv = read.csv("mydata.csv")
mydata_csv

# Read this dataset excluding the first column
mydata_csv = read.csv("mydata.csv")[, -1]
mydata_csv

# Read the dataset from the .txt file
mydata_txt = read.table("mydata.txt", header = TRUE, sep = " ")
mydata_txt

# Read the dataset from the .xlsx file
mydata_xlsx = read.xlsx("mydata.xlsx")
mydata_xlsx

mydata_xlsx = read_xlsx("mydata.xlsx", sheet = "Sheet 1", range = "A1:E20")
mydata_xlsx

# ``, '', "", [ row , column ], $, ~ are used in R.
# ``, '', "" are used for quoting.
# [ row , column ] is used for indexing.
mydata_xlsx[1:2, 2] # First and second row, second column

# $ is used to access the columns of a dataset.
mydata_xlsx$age # Access the age column

mydata_xlsx$family[3:5] # Access the family column from the 3rd to 5th row

mydata_xlsx[-c(2, 6), -2] # Access dataset except 2nd and 6th rows, and 2nd column

# ~ is used in formulas.
# For example, y ~ x means y is modeled as a function of x.
boxplot(age ~ family, data = mydata_xlsx) # Boxplot of age by family

pacman::p_load(car) # qqPlot() function is in the car package, you can load this at the beginning.
qqPlot(age ~ family, data = mydata_xlsx)

mydata_xlsx %>% 
  select(age, family, land) %>% 
  filter(age >40, family %in% c(1, 2)) # Select serial, age, family columns and filter age > 45 and family 1 or 2

# Stack the age and family columns into one column named merged
long = mydata_xlsx %>% 
  pivot_longer(cols = c(age, family) , names_to = 'merged', values_to = 'values') # Combine age and family columns into merged column
long

# Unstack the merged column into age and family columns
short = long %>% 
  pivot_wider(names_from = merged, values_from = values) # Split merged column into age and family columns
short

# Create new variables (income_per_member) using mutate() function
short = short %>% 
  mutate(income_per_member = income/family) # Create new variable income_per_member
head(short)

# Create new variable (status: poor < 1000 income_per_member) using conditional statement
short = short %>% 
  mutate(status = ifelse(income_per_member < 1000, "poor", "rich")) # Create new variable status
head(short)

# Create new variable (status1: poor < 1000 income_per_member) using case_when() function
short = short %>% 
  mutate(status1 = case_when(
    income_per_member < 1000 ~ "poor",
    income_per_member >= 1000 ~ "rich"
  )) # Create new variable status1
head(short)

# Create new variable (status2: poor < 1000, middle 1000-2000, rich > 2000 income_per_member) using ifelse() function and dollar sign
# This is not tidy thinking
# Tidy thinking uses  %>% 
short$status2 = ifelse(short$income_per_member < 1000, "poor", 
    ifelse(short$income_per_member <= 2000, "middle", "rich")) # Create new variable status2
head(short)

# Cut the variable (income_per_member) into intervals (poor < 1000, middle 1000-2000, rich > 2000) using cut() function and save as status3
# right = FALSE means the intervals are left-closed
short$status3 = cut(short$income_per_member, right = FALSE,
                    breaks = c(-Inf, 1000, 1999.9999, Inf), 
                    labels = c("poor", "middle", "rich")) # Create new variable status3
head(short)

# Create a frequency table of status3
short %>% select(status3) %>% table() # tidy way
table(short$status3) # traditional way

# Create a cross-tabulation of status3 and land
table(short$status3, short$land) # Cross-tabulation of status3 and land
# Tidy way
short %>% select(status3, land) %>% table()

# Mean, SD, Median, Min, Max, Range, IQR, Summary functions
# Select a variable using $ sign and apply function
mean(mydata_xlsx$income) # Mean of income_per_member
sd(mydata_xlsx$income) # Standard deviation of income_per_member
median(mydata_xlsx$income) # Median of income_per_member
min(mydata_xlsx$income) # Minimum of income_per_member
max(mydata_xlsx$income) # Maximum of income_per_member
range(mydata_xlsx$income) # Range of income_per_member
IQR(mydata_xlsx$income) # Interquartile range of income_per_member
summary(mydata_xlsx$income) # Summary of income_per_member

# Tidy way
# se and cv: Standard error and coefficient of variation
# se = sd/sqrt(n), cv(%) = sd/mean*100
mydata_xlsx %>% summarise(mean = mean(income),
                         sd = sd(income),
                         iqr = IQR(income),
                         cv = sd(income)*100/mean(income),
                         se = sd(income)/sqrt(n()),
                         frequency = n())

# Tidy way: summary by another factor (land)
# se and cv: Standard error and coefficient of variation
# se = sd/sqrt(n), cv(%) = sd/mean*100
mydata_xlsx %>% group_by(land) %>% summarise(mean = mean(income),
                         sd = sd(income),
                         iqr = IQR(income),
                         cv = sd(income)*100/mean(income),
                         se = sd(income)/sqrt(n()),
                         frequency = n())

# Multiple variable, multiple functions
# transpose using t()
# split column names and sort according to variable
# rename() the columns
# separate() into multiple or unite() into one column
# format using kble() from kableExtra package
mydata_xlsx %>% 
    select(age, income, family) %>% 
    summarise_all(list(mean = mean, max = max)) %>% 
    t() %>% 
    as.data.frame() %>% 
    round(2) %>% 
    rownames_to_column('variable') %>% 
    separate(variable, c('variable', 'statistics'), sep = '_') %>% 
    arrange(variable) %>% 
    pivot_wider(names_from = statistics, values_from = V1)

# Custom function: Create a function to calculate the standard error
# Function name is followed by parentheses and curly brackets.
se = function(x) {
  return(sd(x)/sqrt(length(x)))
}

se(mydata_xlsx$income) %>% round(2) # Standard error of income_per_member

# IQR=Q3−Q1, Q2 = Median
# Boxplot: Q1, Q2, Q3, Whiskers, Outliers
# The lower whisker extends from Q1 to the minimum value within 1.5 × IQR below Q1.
# The upper whisker extends from Q3 to the maximum value within 1.5 × IQR above Q3.
# Outliers are the values beyond the whiskers.
boxplot(mydata_xlsx$income) # Boxplot of income_per_member

head(mydata_xlsx)

# We want to change the variable types and save (overwrite) in the original name of the dataset
mydata_xlsx = mydata_xlsx %>% mutate_if(is.character, as.factor)
head(mydata_xlsx)

mydata_xlsx = mydata_xlsx %>% mutate_at(c('serial', 'family'), as.integer)
head(mydata_xlsx)

mydata_xlsx = mydata_xlsx %>% mutate_at('family', as.numeric) %>% mutate_at(c('age', 'income'), round, 2)
head(mydata_xlsx)

set.seed(123)
n = 30 # Sample size
B = 1000 # Number of samples
sample_means = replicate(B, mean(sample(mydata_xlsx$income, n, replace = TRUE))) # Generate B sample means

# Calculate the mean and standard deviation of the sample means
(mean_sample_means = mean(sample_means)) # Mean of the sample means
(sd_sample_means = sd(sample_means)) # Standard deviation of the sample means
(se_sample_means = sd_sample_means) # Standard error of the sample means

# Plot the sampling distribution of the sample mean
hist(sample_means, breaks = 30, col = "skyblue", main = "Sampling Distribution of the Sample Mean", xlab = "Sample Mean")

# use second bracket to combine multiline codes
{ 
  plot(density(sample_means), col = "blue", 
       main = "Sampling Distribution of the Sample Mean", 
       xlab = "x-axis => Sample Mean or z-values \n Probability = Area under the curve.", 
       ylab = "Density or relative frequency", lty = 1, lwd = 2)
  abline(v = mean_sample_means, col = "red", lty = 2, lwd = 2) # Add a vertical line for the mean of the sample means
  abline(v = mean_sample_means + c(-1, 1)*se_sample_means*1.96, col = "black", lty = 1, lwd = 2) # Add vertical lines for the mean ± 1.96 SE
  text(3000, 0.0005, paste('mean =', round(mean_sample_means)), pos = 4, col = "red", srt = 90) # Add a text for the mean of the sample means
  text(mean_sample_means, 0, '0', pos = 3, col = "black") 
  text(3400, 0.0005, 'mean+1.96se', pos = 4, col = "black", srt=90)
  text(3200, 0.0000, '+1.96', pos = 4, col = "black", srt=0)
  text(2400, 0.0005, 'mean-1.96se', pos = 4, col = "black", srt=90)
  text(2400, 0.0000, '-1.96', pos = 4, col = "black", srt=0)
  text(2900, 0.0020, '95% CI [2557, 3341]', pos = 1, col = "black", srt=0)
}

# Now sort the mydata_xlsx$income dataset by income in ascending order
income_ascending = mydata_xlsx %>% select(income) %>% arrange()
income_ascending

# Now sort the mydata_xlsx$income dataset by income in ascending order
income_ascending = sort(mydata_xlsx$income, decreasing = FALSE)
# Or tidy thinking
income_ascending = mydata_xlsx %>% arrange(income, decreasing = FALSE)
head(income_ascending)

# Find the 2.5th and 97.5th percentiles of the income dataset
income_ascending %>% select(income) %>% summarise_all(quantile, c(0.025, 0.975)) %>% round(2)
# 95% CI [1186.531, 4774.397]

c((mean_income - 1.96*se_income),(mean_income + 1.96*se_income)) %>% round(2)
# 95% CI [2441.083, 3440.422]

# Conditional statement
# if, else if, else
# if (condition) {statement} else {statement}
# if (condition) {statement} else if (condition) {statement} else {statement}

if (mean(mydata_xlsx$income) > 3000) {
  print("The mean income is greater than 3000.")
} else {
  print("The mean income is less than or equal to 3000.")
}

if (mean(mydata_xlsx$income) > 4000) {
  print("The mean income is greater than 4000.")
} else if (mean(mydata_xlsx$income) > 3000) {
  print("The mean income is greater than 3000 but less than or equal to 4000.")
} else {
  print("The mean income is less than or equal to 3000.")
}

# Loop function
# for loop, while loop
# for (variable in sequence) {statement}
for (i in 1:5) {
  print(i)
}

# Iterate the mean function over three columns of the dataset and same the means in means object
means = c() # blank vector for storing calculated means
for (column in c("age", "family", "income")) {
  means = c(means, mean(mydata_xlsx[[column]]))
}
means %>% round(2)

# same output using apply function
apply(mydata_xlsx[c("age", "family", "income")], 2, mean) %>% round(2)

# same output in tidy thinking
mydata_xlsx %>% select(age, family, income) %>% summarise_all(mean) %>% round(2)

# Look the format of the column indexing.
# mydata_xlsx[[column]] is used for indexing the columns of the dataset.
# Here, column = 1, in the first iteration to calculate the mean of age
# column = 2, in the second iteration to calculate the mean of family, and so on.
# Just for your checking
mydata_xlsx[[2]] # First column of the dataset (value only)
mydata_xlsx['age'] # This contains variable name and values
mydata_xlsx[['age']] # This contains only values

# ggplot2 is a powerful package for creating graphics in R.
# It is based on the grammar of graphics, which is a structured way (layer by layer) of creating graphics.
# There are different ways to put the dataset name and aes-thetics as well as geom-etric shapes
ggplot(mtcars) + 
    geom_violin(aes(cyl, mpg, group = cyl))

ggplot(mtcars, aes(cyl, mpg, group = cyl)) + 
    geom_violin()

ggplot(mtcars) + 
    aes(cyl, mpg, group = cyl) + 
    geom_violin()

mtcars %>% ggplot() + 
    aes(cyl, mpg, group = cyl) + 
    geom_violin()

mtcars %>% ggplot(aes(cyl, mpg, group = cyl) ) + 
    geom_violin()

p = mtcars %>% ggplot(aes(cyl, mpg, group = cyl))
p = p + geom_violin()
p = p + theme_bw()
p

# Scattered plot of x=age and y=income
pacman::p_load(jtools) # for theme_apa()
ggplot(mydata_xlsx, aes(x = age, y = income)) + 
  geom_point() + # Scatter plot of age and income
  labs(title = "Scatter plot of Age and Income", x = "Age", y = "Income") + # Add title and axis labels
  theme_bw()

# Boxplot of income by land
ggplot(mydata_xlsx, aes(x = land, y = income)) + 
  geom_boxplot() + # Boxplot of income by land
  labs(title = "Boxplot of Income by Land", x = "Land", y = "Income") + # Add title and axis labels
  theme_apa()

# Histogram of income
ggplot(mydata_xlsx, aes(x = income)) + 
  geom_histogram(binwidth = 500, fill = "skyblue", color = "white") + # Histogram of income
  labs(title = "Histogram of Income", x = "Income", y = "Frequency") + # Add title and axis labels
  theme_test()

# Density plot of income by land
ggplot(mydata_xlsx, aes(x = income, fill = land)) + 
  geom_density(alpha = 0.5) + 
  labs(title = "Fig. 1.1 Density Plot of Income by Land", x = "Income", y = "Density") +

  theme_apa()+
  theme(plot.title = element_text(hjust = 0.5, vjust = -150),
        plot.margin = margin(10, 10, 50, 10),
        legend.position = 'top')

# Bar plot of family by land
ggplot(mydata_xlsx, aes(x = family, fill = land)) + 
  geom_bar(position = "dodge") + # Bar plot of family by land
  labs(title = "Bar Plot of Family by Land", x = "Family", y = "Count") + 
  theme_apa()

# Scatter plot of income by age
ggplot(mydata_xlsx, aes(x = age, y = income, group = 1)) + 
  geom_point(color = "tomato") + 
  labs(title = "Line Plot of Income by Age", x = "Age", y = "Income") + 
  theme_apa() +
  geom_smooth(method = "lm", se = TRUE, formula = y ~ x)

# Create the pie chart
ggplot(pie_data, aes(x = "", y = count, fill = factor(family))) +
  geom_bar(stat = "identity", width = 1) + 
  coord_polar("y", start = 0) +  
  geom_text(aes(label = paste0(round(percentage, 1), '%')   ), 
            position = position_stack(vjust = 0.5), size = 6) + 
  labs(title = "Pie Chart of Family", fill = "Family") + 
  scale_fill_discrete('Family size') +
  theme_void() +  
  theme(legend.position = "right",
        plot.title = element_text(hjust = 0.5, size = rel(2)))

# Using base pie()
x = pie_data$count
percent = paste0(round(pie_data$percentage, 2), '%')
categories = c('One member', 'Two member', 'Three member', 'Four member', 'Five member')
labels = paste0(categories, '\n [', percent, ']')

pie(x, labels, col = rainbow(5))
title(main = "Fig. Pie Chart of Family", adj = 0.5, line = -27)

# Read specific range (B7:N507) from a sheet (wrangling) in excel file (DataSets.xlsx)
df = read_excel('DataSets.xlsx', sheet = 'wrangling', range = 'B7:N507')
head(df)

# We will exclude sr. column
df %>% select(-sr.) %>% head()

# We will calculate means, se, upper and lower limits from this dataset and plot in ggplot
# Define the lower and upper function
# Reorder factors levels: # plot_df$Months = factor(plot_df$Months, 
# levels = c('jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec'), 
# labels = c('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'))

se = function(x){sd(x)/sqrt(length(x))}
lower95 = function(x){mean(x) - 1.96*se(x)}
upper95 = function(x){mean(x) + 1.96*se(x)}

plot_df = df %>% 
    select(-sr.) %>% 
    summarise_all(list(mean = mean, lower = lower95, upper = upper95)) %>% 
    t() %>% as.data.frame() %>% 
    rownames_to_column('variable') %>% 
    separate(variable, c('Months', 'Statistics')) %>% 
    pivot_wider(names_from = Statistics, values_from = V1) %>% 
    mutate(Months = factor(Months, levels = c('jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec'),
                          labels = c('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec')))
plot_df

ggplot(plot_df, aes(x = Months, y = mean, group = 1, fill = Months)) +
    geom_col() +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.25) +
    geom_text(aes(y = upper + 0.01, label = formatC(mean, format = 'f', digits = 2))) +
    theme_bw() +
    theme(legend.position = 'none') +
    labs(x = '', y = 'Average proportion of water scarcity')

# Create a dataset for the spider plot
set.seed(1)
spider_data = as.data.frame(
  matrix(sample(1:10 , 96 , replace = TRUE),
         ncol=8, byrow = TRUE))
spider_data

# Scale the data between 0 and 1 and keep two decimal places
scaled_data <- round(apply(spider_data, 2, scales::rescale), 2)
plot_data <- as.data.frame(scaled_data)
plot_data

# Create the radar plot for the first row of the dataset
ggradar(plot_data[1, ], 
        values.radar = c('0%', '50%', '100%'),
        grid.min = 0, grid.mid = 0.5, grid.max = 1,
        legend.position = "top")

# Create the radar plot for two rows of the dataset 
ggradar(plot_data[2:3,], 
        values.radar = c('0%', '50%', '100%'),
        grid.min = 0, grid.mid = 0.5, grid.max = 1,
        legend.position = "top")

# Create the radar plot for all rows of the dataset using a loop
# Create an empty list to store the plots
options(repr.plot.height=18)
{radar_plots = list()
  
  for (i in 1:nrow(plot_data)) {
    radar_plot = ggradar(plot_data[i,], 
                         values.radar = c('0%', '50%', '100%'),
                         grid.min = 0, grid.mid = 0.5, grid.max = 1,
                         plot.title = row.names(plot_data)[i],
                         legend.position = "top",
                         group.point.size = 3,
                         group.line.width = 1,
                         axis.label.size = 4,
                         grid.label.size = 4,
                         base.size = 15) +
      theme(plot.title = element_text(size = 15, hjust = 0.5))
    
    radar_plots[[i]] = radar_plot # Store the radar plot in the list
  }
  
  wrap_plots(radar_plots, ncol = 3) # Combine the radar plots into a single plot
}