Uncertainty Quantification for ECG Signal Classification

Compute Nonconformity Scores

Use the calibration set of ECG signals to compute conformity scores by using the custom nonconformityScores function, which is defined at the end of this example. The function accepts a deep neural network (netRobust), ECG signals (XCalib), and true class labels for the calibration set (TCalib). For each observation, the function computes the APS nonconformity score, which is the minimum cumulative probability mass that must include the true class when classes are sorted by predicted probability. Larger scores indicate higher nonconformity because a larger probability mass is required to cover the true label.

scores = nonconformityScores(netRobust,XCalib,TCalib);

Compute Conformal Quantile

Use the nonconformity scores and error rate to compute the conformal threshold using the custom conformalQuantile function, which is defined at the end of this example. The function accepts a vector of nonconformity scores (scores) and a desired error rate (errorRate). The function returns the conformal threshold.

The threshold defines the minimum cumulative probability mass that each prediction set must include to ensure that, on average, the set excludes the true label no more than the specified error rate. Therefore, if you include classes until their cumulative probability exceeds this threshold, then your prediction sets will satisfy the error rate guarantee. A lower target error rate requires a higher threshold, resulting in larger prediction sets.

errorRate = 0.75;
threshold = conformalQuantile(scores,errorRate)

threshold = single

0.8314

Plot the distribution of calibration scores and mark the threshold.

figure
histogram(scores, 20)
xline(threshold,"--")
xlabel("Nonconformity Scores")
ylabel("Frequency")
title("Score Distribution with Conformal Threshold")
legend("Scores","Threshold")

Construct Prediction Sets for the Test Set

Apply the conformal threshold to the test set by using the custom conformalPredict function, which is defined at the end of this example. This function accepts a deep neural network (netRobust), a set of ECG signals (XTest), a conformal threshold (threshold), and a vector of class names (classNames).

The function returns three cell arrays. The first array contains the prediction set for each test sample as a categorical vector. The second array contains the softmax probabilities for the labels in each prediction set in the same order. The third array contains all the class scores, in the order specified by the class names.

classNames = categories(TTest);
[YTestSet,YTestScores,allScores] = conformalPredict(netRobust,XTest,threshold,classNames);

For each ECG signal, the prediction set will contain either N, A, or both. If the set contains only one class, then for the specified error rate, the classifier is confident in that label. If the set contains both classes, then the model is unable to confidently distinguish between "N" and "A" for that signal. For example, if the top prediction of a classifier is N, but the prediction set is [A,N], then the model cannot confidently rule out atrial fibrillation. This is especially important in clinical settings, where missing an atrial fibrillation episode (a false negative) can have serious consequences.

Plot the number of signals that have the label set "A", "N", or "A and N" (represented by "AN" in the plot).

combineClasses = strings(size(YTestSet));
for i = 1:numel(YTestSet)
    x = YTestSet{i};
    sortedChars = sort(char(x));
    combineClasses(i) = string(sortedChars');
end

[uniqueElements,~,idx] = unique(combineClasses);
counts = accumarray(idx(:),1);

figure
bar(counts)
xticks(1:numel(uniqueElements))
xticklabels(uniqueElements)

xlabel("Class Set")
ylabel("Count")
title("Frequency of Each Class Set")

To better understand the prediction set, visualize the cumulative probability mass for an individual test signal by using the explainPredictionSet function, which is defined at the end of this example. This function takes as input the class scores, the class names, and the threshold.

The function returns a plot showing the score for each class and a horizontal line showing the threshold. Conformal prediction adds the classes in order of descending probability, and all classes up to and including the first one that causes the cumulative probability to meet or exceed the threshold are included in the prediction set. In the plot, blue bars represent included classes, and red bars represent excluded classes. If both the bars are blue, then the label set for that signal contains both classes.

Explain the prediction set of idx observation.

idx = 711;
explainPredictionSet(allScores{idx},categories(TTest),threshold);

Verify Empirical Coverage

For the test set, check how often the true class label is not included in the prediction set. The proportion of such cases is the empirical error rate.

numTest = numel(TTest);
containsTrueLabel = false(numTest,1);

for i = 1:numTest
    containsTrueLabel(i) = any(YTestSet{i} == TTest(i));
end

empiricalErrorRate = 1- mean(containsTrueLabel);
disp("Target error rate: " + errorRate)

Target error rate: 0.75

disp("Empirical error rate: " + empiricalErrorRate)

Empirical error rate: 0.081368

If the empirical error rate is less than or equal to the target error rate, then that confirms that the conformal prediction method provides valid uncertainty quantification for your ECG classifier.

In the next step of the workflow, you will use the gradient-weighted class activation mapping (Grad-CAM) technique to explain the most likely prediction.

Supporting Functions

Compute Nonconformity Scores

The nonconformityScores function accepts a deep neural network (net), ECG signals (XCalib) and true class labels for the calibration set (TCalib). For each observation, the function computes the APS nonconformity score, which is the minimum cumulative probability mass needed to include the true class when classes are sorted by predicted probability. Larger scores indicate higher nonconformity because the true label requires more probability mass to be covered.

The nonconformityScores function performs these steps:

function scores = nonconformityScores(net,XCalib,TCalib)
YCalib = minibatchpredict(net,XCalib,InputDataFormats="CTB");
[sortedPredictions,sortedIndex] = sort(YCalib,2,"descend"); 
[~,trueClassRank] = max(sortedIndex == double(TCalib),[],2);
sortedSum = cumsum(sortedPredictions,2);
trueClassIndex = sub2ind(size(sortedSum),(1:length(sortedIndex))',trueClassRank);
scores = sortedSum(trueClassIndex);
end

Compute Conformal Quality

The conformalQuantile function accepts a vector of nonconformity scores (calibrationScores) and a desired error rate (errorRate). This function returns the conformal threshold. The threshold defines the minimum cumulative probability mass that each prediction set must include to ensure that, on average, the set excludes the true label at a rate no higher than the specified error rate. Therefore, if you include classes until their cumulative probability exceeds this threshold, your prediction sets will satisfy the error rate guarantee. A lower target error rate requires a higher threshold, resulting larger prediction sets.

The conformalQuantile function performs these steps:

function threshold = conformalQuantile(calibrationScores,errorRate)
sortedScores = sort(calibrationScores,"ascend");
numObservations = numel(calibrationScores);
coverageIndex = ceil((numObservations+1)*(1-errorRate));
threshold = sortedScores(coverageIndex);
end

Apply Conformal Prediction Threshold

The conformalPredict function accepts a deep neural network (net), a set of ECG signals (X), a conformal threshold (threshold), and a vector of class names (classNames). The function returns three cell arrays. The first array contains the prediction set for each test sample as a categorical vector. The second array contains the softmax probabilities for the labels in each prediction set in the same order. The third array contains all the class scores, in the order specified by the class names.

This function performs these steps:

function [predictionSets,predictionScores,allScores] = conformalPredict(net,X,threshold,classNames)
Y = minibatchpredict(net,X,InputDataFormats="CTB");

[sortedPredictions,sortedIndex] = sort(Y,2,"descend");
sortedSum = cumsum(sortedPredictions,2);

meetsThreshold = sortedSum >= threshold;
[~,cutoffIndex] = max(meetsThreshold,[],2);

numSamples = size(Y,1);
predictionSets = cell(1,numSamples);
predictionScores = cell(1,numSamples);
allScores = cell(1,numel(classNames));

for i = 1:numSamples
    includedIndex = sortedIndex(i,1:cutoffIndex(i));
    includedLabels = classNames(includedIndex);
    predictionSets{i} = categorical(includedLabels,classNames);
    predictionScores{i} = sortedPredictions(i, 1:cutoffIndex(i));
    allScores{i} = Y(i,:);
end
end

Explain Prediction Set

The explainPredictionSet function takes as input the class scores for a single signal (signalScores), class names (classes), and the threshold (threshold). The function returns a plot that shows the score for each class, with a horizontal line marking the threshold. Conformal prediction adds the classes in order of descending probability, and all classes up to and including the first one that causes the cumulative probability to meet or exceed the threshold are included in the prediction set. In the plot, blue bars represent included classes, and red bars represent excluded classes. If both the bars are blue, then the label set for that signal contains both classes.

The function performs these steps:

function explainPredictionSet(signalScores,classes,threshold)

% Sort probabilities
[sortedProbs,sortedIdx] = sort(signalScores,"descend");
sortedNames = classes(sortedIdx);
cumProb = cumsum(sortedProbs);

% Find smallest top-k set meeting threshold
cutoff = find(cumProb >= threshold,1,"first");

if isempty(cutoff)
    cutoff = numel(sortedProbs);
end

% Plot results
RGB = orderedcolors("gem");
colorBar = repmat(RGB(2,:),numel(cumProb),1);  
colorBar(1:cutoff,:) = repmat(RGB(1,:),cutoff,1);

figure
labels = ["Included","Excluded"];

hold on
if cutoff ~= numel(classes)
    for i = 1:numel(sortedIdx)
        bb(i) = bar(i,sortedProbs(i),FaceColor=colorBar(i,:));
    end
    legend(bb,labels{cutoff:end},Location="eastoutside");
else
    bb = bar(1:cutoff,sortedProbs,FaceColor=colorBar(1,:));
    legend(bb,labels{1},Location="eastoutside");
end

yline(threshold,"-",DisplayName="Threshold");
xticks(1:numel(sortedNames))
xticklabels(sortedNames)
ylabel("Score")
xlabel("Class")
ylim([0,1])
title("Classes Included in Label Set")
hold off
end

References

[1] Romano, Yaniv, et al. “Classification with Valid and Adaptive Coverage.” arXiv:2006.02544, arXiv, 3 Jun. 2020. arXiv.org, https://doi.org/10.48550/arXiv.2006.02544.