Expected overlaps

R/tests.R

@@ -116,3 +116,167 @@ upset_test = function(
     result = comparison[order(comparison$fdr), ]
     result
 }
+
+
+test_set_overlaps_of_degree = function(
+    data, intersect, degree, remaining_degrees, diagnosis_rates, test, min_size, max_degree,
+    sets_as_columns=FALSE
+) {
+    if (remaining_degrees == 0) {
+        data_subset = data[diagnosis_rates$name]
+        diagnosed_with_all = sum(apply(data_subset, 1, all), na.rm=TRUE)
+
+        if (is.na(diagnosed_with_all)) {
+            diagnosed_with_all = 0
+        }
+
+        if (diagnosed_with_all < min_size) {
+            return(list())
+        }
+
+        having_data_for_all = sum(apply(!is.na(data_subset), 1, all), na.rm=TRUE)
+
+        expected_overlap = prod(diagnosis_rates$rate)
+
+        if (degree == 1) {
+            # testing overlap does not make sense for degree=1 and could mislead FDR calculations
+            p = list(
+                statistic=NA,
+                p.value=NA
+            )
+        } else {
+            p = test(
+                x=c(
+                    'with_all_overlaping'=diagnosed_with_all,
+                    'without_all_overlapping'=having_data_for_all - diagnosed_with_all
+                ),
+                p=c(
+                    'with_all_overlapping'=expected_overlap,
+                    'without_all_overlapping'=1 - expected_overlap
+                )
+            )
+        }
+        result = list(
+            'expected_overlap'=expected_overlap,
+            'expected_n'=expected_overlap * having_data_for_all,
+            'observed_n'=diagnosed_with_all,
+            'observed_overlap'=diagnosed_with_all / having_data_for_all,
+            'sample_size'=having_data_for_all,
+            'identifier'=paste(sanitize_names(diagnosis_rates$name), collapse='-'),
+            'statistic'=p$statistic,
+            'p_value'=p$p.value
+        )
+
+        if (sets_as_columns) {
+            sets = c(
+                diagnosis_rates$name,
+                rep(NA, max_degree - length(diagnosis_rates$name))
+            )
+            names(sets) = paste0('set_', seq(1, max_degree))
+            result = c(sets, result)
+        }
+
+        return (result)
+    } else {
+        expected_overlaps = list()
+
+        for (i in 1:length(intersect)) {
+            disease = intersect[i]
+
+            if (disease %in% diagnosis_rates$name) {
+                next()
+            }
+
+            data_disease = data[[disease]]
+            diagnoses = sum(data_disease, na.rm=TRUE)
+
+            if (diagnoses < min_size) {
+                next()
+            }
+
+            status_known_for = sum(!is.na(data_disease))
+            diagnosis_rate = diagnoses / status_known_for
+
+            diagnosis_rates_temp = rbind(
+                diagnosis_rates,
+                list(
+                    'name'=disease,
+                    'rate'=diagnosis_rate,
+                    'diagnoses'=diagnoses,
+                    'status_known_for'=status_known_for
+                )
+            )
+
+            expected_overlap = test_set_overlaps_of_degree(
+                data=data,
+                intersect=intersect[i:length(intersect)],
+                diagnosis_rates=diagnosis_rates_temp,
+                degree=degree,
+                remaining_degrees=remaining_degrees - 1,
+                test=test,
+                min_size=min_size,
+                sets_as_columns=sets_as_columns,
+                max_degree=max_degree
+            )
+
+            if (length(expected_overlap)) {
+                expected_overlaps[[length(expected_overlaps) + 1]] = expected_overlap
+            }
+        }
+        return (data.frame(do.call(rbind, expected_overlaps)))
+    }
+}
+
+
+test_set_overlaps = function(
+    data, intersect,
+    min_degre=1, max_degree=3,
+    test=chisq.test,
+    min_size=0,
+    encode=FALSE,
+    fdr_method='fdr',
+    sets_in_columns=FALSE
+) {
+    data = data[intersect]
+
+    # TODO: shared?
+    is_column_logical = sapply(data[, intersect], is.logical)
+    if (any(!is_column_logical)) {
+        non_logical = names(is_column_logical[is_column_logical == FALSE])
+        data[, non_logical] = sapply(data[, non_logical], as.logical)
+    }
+
+    if (encode) {
+        intersect = encode_names(intersect, avoid=c())
+        colnames(data) = intersect
+    }
+
+    all_expected_overlaps = list()
+
+    for (degree in seq(min_degre, max_degree)) {
+
+        expected_overlaps = test_set_overlaps_of_degree(
+            data=data,
+            intersect=intersect,
+            degree=degree,
+            remaining_degrees=degree,
+            diagnosis_rates=data.frame(name=character(), rate=numeric()),
+            test=test,
+            min_size=min_size,
+            max_degree=max_degree,
+            sets_as_columns=sets_in_columns
+        )
+
+        if (length(expected_overlaps)) {
+            expected_overlaps$degree = degree
+            all_expected_overlaps[[length(all_expected_overlaps) + 1]] = expected_overlaps
+        }
+    }
+
+    df = data.frame(do.call(rbind, all_expected_overlaps))
+    df = data.frame(lapply(df, unlist))
+    df$enrichment = df$observed_overlap - df$expected_overlap
+    df$phi = sqrt(df$statistic^2 / df$sample_size)
+    df$fdr = p.adjust(df$p_value, method = fdr_method)
+    df = df[order(-df$phi), ]
+}