
library(readr)
library(dplyr)
library(purrr)
library(ggplot2)


violinPlot <- function(data, varx, vary, title, xlab, ylab, alpha=0.5, jitter=F) {
    
    if (jitter) {
        
        p <- ggplot(data=data, aes_string(x=varx, y=vary)) +
            geom_violin(show.legend = T) +
            geom_jitter(width=0.15, alpha=alpha) +
            labs(title=title, xlab=xlab, ylab=ylab) +
            theme(
                plot.title = element_text(size=24),
                axis.title = element_text(size=18),
                axis.text = element_text(size=14)
            )
        
    } else {
        
        p <- ggplot(data=data, aes_string(x=varx, y=vary)) +
            geom_violin(show.legend = T) +
            geom_boxplot(width=0.1, color="darkgrey", alpha=0.5) +
            labs(title=title, xlab=xlab, ylab=ylab) +
            theme(
                plot.title = element_text(size=24),
                axis.title = element_text(size=18),
                axis.text = element_text(size=14)
            )
        
    }
    
    return(p)
}

# Zainspirowane następującym kodem: 
# https://github.com/shakedzy/dython/blob/06aa19f3332de4f80478f5e8bf3ba868f7ddfb63/dython/nominal.py#L194
correlation_ratio <- function(categories, measurements, numeric_replace_value = 0) {
    
    measurements[is.na(measurements)] <- numeric_replace_value
    
    fcat <- as.numeric(categories)
    cat_num <- max(fcat)
    y_avg_array <- rep(0, cat_num)
    n_array <- rep(0, cat_num)
    
    for (i in 1:(cat_num)) {
        cat_measures <- measurements[fcat==i]
        n_array[i] <- length(cat_measures)
        y_avg_array[i] = mean(cat_measures)
    }
    
    y_total_avg <- sum(y_avg_array * n_array) / sum(n_array)
    
    numerator <- sum((y_avg_array - y_total_avg)^2 * n_array)
    
    denominator <- sum((measurements - y_total_avg)^2)
    
    eta <- ifelse(numerator == 0, 0, sqrt(numerator / denominator))
    
    return(eta)
    
}


calc_corr <- function(data, row_name, col_name, numeric_replace_value = 0, theil_uncert=TRUE) {

    row_vec <- data[[row_name]]
    col_vec <- data[[col_name]]
    
    row_data_type <- class(row_vec)
    col_data_type <- class(col_vec)

    corr <- NA
    
    if (row_name == col_name) {
        
        corr <- 1.0
    
    } else if (row_data_type == 'numeric' & col_data_type == 'numeric') {
    
        col_vec[is.na(col_vec)] <- numeric_replace_value
        row_vec[is.na(row_vec)] <- numeric_replace_value
        
        c <- tibble(row_vec, col_vec)
        names(c) <- c(row_name, col_name)
        
        corr <- (c %>% corrr::correlate(method = 'pearson', quiet = T))[[1,3]]
        
    } else if (row_data_type == 'numeric' & (col_data_type == 'character' | col_data_type == 'factor')) {
        
        if (col_data_type == 'character') {
            col_vec <- addNA(as.factor(col_vec))
        }
        
        corr <- correlation_ratio(categories = col_vec, measurements = row_vec,
                                  numeric_replace_value = 0)
        
    } else if ((row_data_type == 'character' | row_data_type == 'factor') & col_data_type == 'numeric') {
        
        if (row_data_type == 'character') {
            row_vec <- addNA(as.factor(row_vec))
        }
        
        corr <- correlation_ratio(categories = row_vec, measurements = col_vec,
                                  numeric_replace_value = 0)
        
    } else if ((row_data_type == 'character' | row_data_type == 'factor') & (col_data_type == 'character' | col_data_type == 'factor')) {
        
        if (row_data_type == 'character') {
            row_vec <- addNA(as.factor(row_vec))
        }
        
        if (col_data_type == 'character') {
            col_vec <- addNA(as.factor(col_vec))
        }
        
        if (theil_uncert) {
            corr <- DescTools::UncertCoef(row_vec, col_vec, direction = 'row')
        } else {
            corr <- rstatix::cramer_v(x=row_vec, y=col_vec)
        }
        
    }
    
    return(corr)

}


# Załaduj zestaw danych o katastrofie Titanica
dataset_url <- 'http://bit.ly/titanic-dataset-csv'

tbl <- read_csv(dataset_url)
tbl


# Przekształć zmienne kategoryczne na tekstowe kolumny
tbl <- tbl %>% 
    mutate( across(c('Survived', 'Pclass'), as.factor) )

tbl

# Obliczenie współczynnika V Craméra dla kolumn Survived i Pclass
rstatix::cramer_v(x=tbl$Survived, y=tbl$Pclass)

# Weryfikujemy, czy funkcja V Craméra jest symetryczna
rstatix::cramer_v(x=tbl$Survived, y=tbl$Pclass) == rstatix::cramer_v(x=tbl$Pclass, y=tbl$Survived)


# Możesz także narysować dla tych zmiennych wykres mozaikowy
vcd::mosaic(~ Survived + Pclass, data = tbl,
            main = "Survived VS Passenger class", shade = TRUE)

# W ten sposób możesz uzyskać interesującą tabelę krzyżową
ct <- tbl %>%
    sjPlot::sjtab(fun = "xtab", var.labels=c("Survived", "Pclass"),
                  show.row.prc=T, show.col.prc=T, show.summary=T, show.exp=T, show.legend=T)

ct



# Skorzystaj z asymetrii funkcji U Theila poprzez obliczenie jej dla tych samych zmiennych.
# To jest wartość U(Survived|Pcalss), tzn. "U dla wartości Survived oraz danej wartości Pclass"
DescTools::UncertCoef(tbl$Survived, tbl$Pclass, direction = 'row')

# Just check that the opposite direction gives you a different result
DescTools::UncertCoef(tbl$Pclass, tbl$Survived, direction = 'row')

# Możesz także wyświetlić rozkład trzeciego wymiaru (w tym przypadku Survived)
# poprzez pokolorowanie połowy wykresu skrzypcowego
violinPlot(data = tbl, varx = 'Pclass', vary = 'Age',
           title = 'Passenger Age VS Passenger Class',
           xlab = 'Pclass', ylab = 'Age', alpha = 0.4)

# Obliczenie współczynnika korelacji pomiędzy kategoryczną zmienną Pclass, a 
# liczbową zmienną Age
correlation_ratio( categories = tbl$Pclass, measurements = tbl$Age, numeric_replace_value = 0)

# Pokażmy, w jaki sposób zmieni się współczynnik korelacji po zmianie rozproszenia obserwacji dla każdej kategorii 
t1 <- tibble(
    topic=c('Algebra','Algebra','Algebra','Algebra','Algebra','Geometry','Geometry','Geometry','Geometry','Statistics','Statistics','Statistics','Statistics','Statistics','Statistics'),
    score=c(45,70,29,15,21,40,20,30,42,65,95,80,70,85,73)
)


t1 %>%
    violinPlot(varx = 'topic', vary = 'score',
               title = 'η=0.84', xlab = '', ylab = '',
               alpha = 0.8, jitter = T)

t2 <- tibble(
    topic=c('Algebra','Algebra','Algebra','Algebra','Algebra','Geometry','Geometry','Geometry','Geometry','Statistics','Statistics','Statistics','Statistics','Statistics','Statistics'),
    score=c(36,36,36,36,36,33,33,33,33,78,78,78,78,78,78)
)

t2 %>%
    violinPlot(varx = 'topic', vary = 'score',
               title = 'η=0.84', xlab = '', ylab = '',
               alpha = 0.8, jitter = T)


# Utwórz dwie ramki danych z jedną kolumną zawierającą 
# wartości w postaci nazw kolumn obiektu tibble
row <- data.frame(row=names(tbl))
col <- data.frame(col=names(tbl))

# Utwórz z poprzednich dwóch ramek danych ramkę danych złączenia krzyżowego 
ass <- tidyr::crossing(row, col)

# Dodaj kolumnę corr zawierającą wartości korelacji 
corr_tbl <- ass %>% 
    mutate( corr = map2_dbl(row, col, ~ calc_corr(data = tbl, row_name = .x, col_name = .y, theil_uncert = T)) )

corr_tbl

# Wykreślenie mapy cieplnej z wykorzystaniem obiektu tibble korelacji 
corr_tbl %>% 
    ggplot( aes(x=row, y=col, fill=corr) ) +
    geom_tile() +
    geom_text(aes(row, col, label = round(corr, 2)), color = "white", size = 4)

