We can compare P1-P3 to the manual annotation using external spatial
+metrics.
+
Dataset level
+
+
Let’s first calculate two dataset-level metrics, SpatialARI and
+SpatialAccuracy:
+
+res3 <- getSpatialExternalMetrics(true=data$reference, pred=data$P3,
+ location=data[,c("pxl_col_in_fullres","pxl_row_in_fullres")], level="dataset",
+ metrics=c("SpatialARI","SpatialAccuracy"), k=6,
+ fuzzy_true = TRUE, fuzzy_pred = FALSE)
+
+res2 <- getSpatialExternalMetrics(true=data$reference, pred=data$P2,
+ location=data[,c("pxl_col_in_fullres","pxl_row_in_fullres")], level="dataset",
+ metrics=c("SpatialARI","SpatialAccuracy"), k=6,
+ fuzzy_true = TRUE, fuzzy_pred = FALSE)
+
+res1 <- getSpatialExternalMetrics(true=data$reference, pred=data$P1,
+ location=data[,c("pxl_col_in_fullres","pxl_row_in_fullres")], level="dataset",
+ metrics=c("SpatialARI","SpatialAccuracy"), k=6,
+ fuzzy_true = TRUE, fuzzy_pred = FALSE)
+
+cbind(bind_rows(list(res1, res2, res3), .id="P")) %>%
+ pivot_longer(cols=c("SpatialARI", "SpatialAccuracy"),
+ names_to="metric", values_to="value") %>%
+ ggplot(aes(x=P, y=value, group=metric)) +
+ geom_point(size=3, aes(color=P)) +
+ facet_wrap(~metric, scales = "free") +
+ theme_bw() + labs(x="Prediction")
+
+
+
Class/cluster level
+
+
We can further calculate the class/cluster-level metrics, SpatialAWH
+and SpatialAWC, to get more insights about the errors our predictions
+make:
+
+res3 <- getSpatialExternalMetrics(true=data$reference, pred=data$P3,
+ location=data[,c("pxl_col_in_fullres","pxl_row_in_fullres")], level="class",
+ metrics=c("SpatialAWH","SpatialAWC"), k=6,
+ fuzzy_true = TRUE, fuzzy_pred = FALSE)
+
+res2 <- getSpatialExternalMetrics(true=data$reference, pred=data$P2,
+ location=data[,c("pxl_col_in_fullres","pxl_row_in_fullres")], level="class",
+ metrics=c("SpatialAWH","SpatialAWC"), k=6,
+ fuzzy_true = TRUE, fuzzy_pred = FALSE)
+
+res1 <- getSpatialExternalMetrics(true=data$reference, pred=data$P1,
+ location=data[,c("pxl_col_in_fullres","pxl_row_in_fullres")], level="class",
+ metrics=c("SpatialAWH","SpatialAWC"), k=6,
+ fuzzy_true = TRUE, fuzzy_pred = FALSE)
+res1
+## SpatialAWH SpatialAWC class cluster
+## 1 NA 0.7963689 1 NA
+## 2 NA 0.8145576 2 NA
+## 3 NA 0.7561886 3 NA
+## 4 NA 0.8184988 4 NA
+## 5 NA 0.7949846 5 NA
+## 6 NA 0.8179805 6 NA
+## 7 NA 0.7884382 7 NA
+## 8 0.7704762 NA NA 1
+## 9 0.8727733 NA NA 2
+## 10 0.8274801 NA NA 3
+## 11 0.5853277 NA NA 4
+## 12 0.8989078 NA NA 5
+## 13 0.6277691 NA NA 6
+## 14 0.6240081 NA NA 7
+
Note that the indices in columns “class” and “cluster” correspond to
+the levels of original factors passed to true and
+pred. We align them back to the previous factor values, and
+then plot them in heatmap:
+
+awh1 <- na.omit(res1[,c("SpatialAWH", "cluster")]) %>% mutate(cluster = levels(data$P1)[cluster])
+awh2 <- na.omit(res2[,c("SpatialAWH", "cluster")]) %>% mutate(cluster = levels(data$P2)[cluster])
+awh3 <- na.omit(res3[,c("SpatialAWH", "cluster")]) %>% mutate(cluster = levels(data$P3)[cluster])
+
+awh <- cbind(bind_rows(list(awh1, awh2, awh3), .id="P")) %>%
+ pivot_wider(names_from = cluster, values_from = SpatialAWH) %>%
+ subset(select = -c(P))
+awh <- as.matrix(awh)
+rownames(awh) <- c("P1", "P2", "P3")
+awh <- awh[,c("1", "2", "3", "4", "5", "6", "7", "8")]
+
+column_ha = HeatmapAnnotation(class=c("1", "2", "3", "4", "5", "6", "7", "8"),
+ col = list(class = my_cols),
+ show_annotation_name = FALSE,
+ show_legend = FALSE)
+
+p4 <- awh %>%
+ Heatmap(
+ name = "AWH",
+ cluster_columns = FALSE,
+ cluster_rows = FALSE,
+ column_order = c("1","2","3","4","5","6","7","8"),
+ bottom_annotation = column_ha,
+ width = ncol(awh)*unit(6, "mm"),
+ height = nrow(awh)*unit(7, "mm"))
+
+p4 <- grid.grabExpr({
+ draw(p4, merge_legend = TRUE, padding = unit(c(12, 2, 2, 2), "mm"))
+})
+
+awc1 <- na.omit(res1[,c("SpatialAWC", "class")]) %>% mutate(class = levels(data$reference)[class])
+awc2 <- na.omit(res2[,c("SpatialAWC", "class")]) %>% mutate(class = levels(data$reference)[class])
+awc3 <- na.omit(res3[,c("SpatialAWC", "class")]) %>% mutate(class = levels(data$reference)[class])
+
+awc <- cbind(bind_rows(list(awc1, awc2, awc3), .id="P")) %>%
+ pivot_wider(names_from = class, values_from = SpatialAWC) %>%
+ subset(select = -c(P))
+awc <- as.matrix(awc)
+rownames(awc) <- c("P1", "P2", "P3")
+
+col_ls <- my_cols[1:7]
+names(col_ls) <- c("WM", "Layer6", "Layer5", "Layer4", "Layer3", "Layer2", "Layer1")
+column_ha = HeatmapAnnotation(class=c("WM", "Layer6", "Layer5", "Layer4", "Layer3", "Layer2", "Layer1"),
+ col = list(class = col_ls),
+ show_annotation_name = FALSE,
+ show_legend = FALSE)
+p5 <- awc %>%
+ Heatmap(
+ name = "AWC",
+ cluster_columns = FALSE,
+ cluster_rows = FALSE,
+ column_order = c("WM", "Layer6", "Layer5", "Layer4", "Layer3", "Layer2", "Layer1"),
+ bottom_annotation = column_ha,
+ width = ncol(awc)*unit(6, "mm"),
+ height = nrow(awc)*unit(7, "mm"))
+
+p5 <- grid.grabExpr({
+ draw(p5, merge_legend = TRUE)
+})
+
+
+
The class-level AWC highlights that in P2, Layer3 has low
+completeness. This align with our simulation that Layer3 is splitted
+into 3 clusters in P2. Similarly, the cluster-level AWH highlights that
+in P3, cluster 3 has low homogeneity, consistent with the merging of
+layer 4 and 5.
+
+
Element level
+
+
One can also calculate element-level metric, SPC, for
+visualization.
+
+res1 <- cbind(getSpatialExternalMetrics(true=data$reference, pred=data$P1,
+ location=data[,c("pxl_col_in_fullres","pxl_row_in_fullres")], level="element",
+ metrics=c("SpatialSPC"), k=6,
+ fuzzy_true = TRUE, fuzzy_pred = FALSE, useNegatives = FALSE),
+ data[,c("pxl_col_in_fullres", "pxl_row_in_fullres")])
+
+res2 <- cbind(getSpatialExternalMetrics(true=data$reference, pred=data$P2,
+ location=data[,c("pxl_col_in_fullres","pxl_row_in_fullres")], level="element",
+ metrics=c("SpatialSPC"), k=6,
+ fuzzy_true = TRUE, fuzzy_pred = FALSE, useNegatives = FALSE),
+ data[,c("pxl_col_in_fullres", "pxl_row_in_fullres")])
+
+res3 <- cbind(getSpatialExternalMetrics(true=data$reference, pred=data$P3,
+ location=data[,c("pxl_col_in_fullres","pxl_row_in_fullres")], level="element",
+ metrics=c("SpatialSPC"), k=6,
+ fuzzy_true = TRUE, fuzzy_pred = FALSE, useNegatives = FALSE),
+ data[,c("pxl_col_in_fullres", "pxl_row_in_fullres")])
+
+
+cbind(bind_rows(list(res1, res2, res3), .id="P")) %>%
+ pivot_longer(cols=c("SpatialSPC"),
+ names_to="metric", values_to="value") %>%
+ ggplot(aes(x = pxl_col_in_fullres, y = - pxl_row_in_fullres, color = value)) +
+ scale_colour_gradient(high="white", low ="deeppink4") +
+ geom_point(size=0.3) +
+ facet_wrap(~P, scales = "free") +
+ theme_bw() + labs(x="Prediction", y="", color="SpatialSPC")
+
+This clear highlights the low concordance regions in each prediction as
+expected.
+