######################## CALCULATING CLIFF'S DELTA ########################

# Supplementary material for:
# Meissel, K., & Yao, E. S. (2024). Using Cliff's Delta as a non-parametric 
# effect size measure: An accessible web app  and R tutorial. Practical 
# Assessment, Research, & Evaluation. 

# Step 1: Install and load the effsize package ----------------------------

# install.packages("effsize") # run if not already installed
library(effsize) # load package
help(effsize) # view package documentation

# Step 2: Load prepared dataset -------------------------------------------

# Dataset should be a CSV file with two columns of data, and a header row.

df <- read.csv("C:/Users/UserName/Documents/example-data.csv", # select CSV file
               fileEncoding = "UTF-8-BOM") # prevents "ï.." added in front of Column 1 name

head(df) # view first 6 rows of data

# Step 3: Calculate Cliff's delta -----------------------------------------

cliff.delta(df$School2, df$School1)

# If School 2 is specified as the first argument (i.e., the comparison group): 
# - Positive result means the distribution of School 2 is *larger* than School 1.
# - Negative result means the distribution of School 2 is *smaller* than School 1.
# If School 1 is specified as the first argument, the signs will be reversed.

# Step 4: Calculate additional statistics (optional) ----------------------

# Number of non-missing cases 
length(na.omit(df$School1)) # n for School 1
length(na.omit(df$School2)) # n for School 2

# Probability of superiority
result <- cliff.delta(df$School2, df$School1, return.dm = TRUE) # calculate dominance matrix
sum(result$dm <  0)/length(result$dm) # probability School 1 > School 2
sum(result$dm == 0)/length(result$dm) # probability School 1 = School 2
sum(result$dm >  0)/length(result$dm) # probability School 1 < School 2

# Cohen's d converted from Cliff's delta
result <- cliff.delta(df$School2, df$School1) # save Cliff's delta results
2*qnorm(-1/(result$estimate-2)) # Cohen's d converted from Cliff's delta
2*qnorm(-1/(result$conf.int[[1]]-2)) # Lower 95% CI for Cohen's d
2*qnorm(-1/(result$conf.int[[2]]-2)) # Upper 95% CI for Cohen's d

# Cohen's d calculated from raw data (with Hedges' correction)
cohen.d(df$School2, df$School1, pooled = TRUE, hedges = TRUE, na.rm = TRUE)
?cohen.d # view options for cohen.d function