function [ out, roc ] = stat_calc_struct( pred, target, lessStr )
%stat_calc_struct Unleashes many statistical tests to determine the efficacy
% of a given set of predictions, pred, on a set of dichotomous targets,
% target. Predictions should range between 0-1, and represent the
% probability estimate that the target is 1.
% [ out ] = stat_calc_struct(pred,target) calculates the following
% statistical tests and stores them in the corresponding structure fields
% listed:
%
% AUROC - Area under the receiver operator (characteristic) curve
% HL - Hosmer-Lemeshow Statistic
% HLD - Hosmer-Lemeshow Statistic normalized to variable range
% acc - Mean accuracy of the prediction using a threshold of 0.5
% R - Shapiro's R, the geometric mean of positive outcomes
% B - Brier's score, the mean square error
% SMR - Standardized mortality ratio. Mean observed outcomes
% divided by mean predicted outcomes.
% coxA - Alpha coefficient in cox calibration.
% coxB - Beta coefficient in cox calibration.
% Cox calibration involves a linear regression of the
% predictions onto the targets. coxA is the constant
% offset, while coxB is the slope.
% TP - True positives.
% FP - False positives.
% TN - True negatives.
% FN - False negatives.
% sens - Sensitivity. Calculated as: TP/(TP+FN)
% spec - Specificity. Calculated as: TN/(TN+FP)
% ppv - Positive predictive value. Calculated as: TP/(TP+FP)
% npv - Negative predictive value. Calculated as: TN/(TN+FN)
% ll - Log likelihood
%
% [ out,roc ] = stat_calc_struct(pred,target) additionally calculates the
% values pairs which form the ROC curve.
% roc.x - Stores 1-specificity, plotted on the x-axis
% roc.y - Stores Sensitivity, plotted on the y-axis
% Copyright 2011 Alistair Johnson
% $LastChangedBy: alistair $
% $LastChangedDate: 2012-06-01 17:44:43 +0100 (Fri, 01 Jun 2012) $
% $Revision: 25 $
% Originally written on GLNXA64 by Alistair Johnson
% Contact: alistairewj@gmail.com
out=[]; roc=[];
%=== Save computation time by skipping cox
if nargin>2 && strcmp(lessStr,'less')
lessFlag = true;
else
%=== Default, calculate cox calibration curves
lessFlag = false;
end
%=== Add in some error checks
if isempty(pred)
out=defaultStruct();
return;
elseif sum(isnan(pred))==length(pred)
% Pred is all NaNs
out = defaultStruct();
return;
end
if max(pred>1) || max(pred<0)
warning('stat_calc_struct:IllConditioned',...
['\nThis function will return potentially erroneous values \n' ...
'when predictions are not in the range [0,1].\n']);
end
% c-stat
out.AUROC=cstat(pred,target);
% Shapiro's R (geometric mean of positive targetcomes)
% transform into naural logarithm, take arithmetic mean, transform back
out.R=exp(sum(log(pred(target==1)))/sum(target));
% Brier's Score, B (mean square error)
out.B=mean((target-pred).^2);
% Accuracy
out.acc=1-(sum(abs(round(target-pred)),1)/length(target));
% hosmer-lemeshow
out.HL=lemeshow(pred,target);
out.HLmine=hlCStat(pred,target);
out.HLD=lemeshowNormalized(pred,target);
% SMR
out.SMR=sum(target)/sum(pred);
idx1 = target==1;
idx0 = target==0;
out.TP=sum(pred(idx1) >= 0.5);
out.FN=sum(pred(idx1) < 0.5);
out.FP=sum(pred(idx0) >= 0.5);
out.TN=sum(pred(idx0) < 0.5);
out.sens=out.TP/(out.TP+out.FN);
out.spec=out.TN/(out.TN+out.FP);
out.ppv=out.TP/(out.TP+out.FP);
out.npv=out.TN/(out.TN+out.FN);
out.ll = sum(-target.*log(pred)-(1-target).*log(1-pred));
% cox linear regression testing
if ~lessFlag
b=glmfit(pred,target,'binomial','link','identity');
out.coxA=b(1); out.coxB=b(2);
end
if nargout>1
if any(idx1) && any(idx0)
[roc.x,roc.y]=perfcurve(target,pred,1);
else
roc.x=[]; roc.y=[];
end
end
end
function [ H ] = lemeshow( pred, outcome )
data = [outcome,pred];
D=10; % Use D-iles (i.e., deciles as in Lemeshow's original paper)
N=length(data(:,1));
M=floor(N/D); % M should always be 400 for the challenge
data=sortrows(data,2); % rearrange the matrix by estimated risk
CM=zeros(D,3)+NaN;
centroid=zeros(D-1,1)+NaN;
for i=1:D
%Set number of deaths proportional to Y true value
ind2=i*M; % highest bin in decile D
ind1=ind2-M+1; % lowest bin
Obs=sum([data(ind1:ind2,1)]); % number of observed deaths within decile D
centroid(i)=mean(data(ind1:ind2,2)); % mean estimated risk within decile
Exp=M*centroid(i); % expected number of deaths within decile
Htemp=(Obs-Exp)^2/(M*centroid(i)*(1-centroid(i)) + 0.001);
CM(i,:)=[Obs Exp Htemp];
end
% Hosmer-Lemeshow H statistic, normalized by range of decile risk
H=sum(CM(:,3));
end
function [ H ] = lemeshowNormalized( pred, outcome )
data = [outcome,pred];
D=10; % Use D-iles (i.e., deciles as in Lemeshow's original paper)
N=length(data(:,1));
M=floor(N/D); % M should always be 400 for the challenge
data=sortrows(data,2); % rearrange the matrix by estimated risk
CM=zeros(D,3)+NaN;
centroid=zeros(D-1,1)+NaN;
for i=1:D
%Set number of deaths proportional to Y true value
ind2=i*M; % highest bin in decile D
ind1=ind2-M+1; % lowest bin
Obs=sum([data(ind1:ind2,1)]); % number of observed deaths within decile D
centroid(i)=mean(data(ind1:ind2,2)); % mean estimated risk within decile
Exp=M*centroid(i); % expected number of deaths within decile
Htemp=(Obs-Exp)^2/(M*centroid(i)*(1-centroid(i)) + 0.001);
CM(i,:)=[Obs Exp Htemp];
end
% Hosmer-Lemeshow H statistic, normalized by range of decile risk
H=sum(CM(:,3))/(centroid(10)-centroid(1));
end
function [ G ] = hlCStat( pred, outcome )
%hltest hosmer lemeshow test
% [ G ] = hltest(pred, outcome) calculates the hosmer-lemeshow C statistic
% for a given set of predictions and outcomes. Steps include:
% 1) Sort data from lowest pred to highest
% 2) Split the data into 10 bins with equal size
% 3) Calculate the mean square error for each bin
% 4) Divide the each bin by the number of samples in each bin
% 5) Divide each bin by the mean prediction in that bin
% 6) Divide each bin by (1-mean prediction) in that bin
% 7) Sum across all bins
%
bins=30;
% sort predictions from lowest to highest
[pred,ind]=sort(pred,1,'ascend');
outcome=outcome(ind);
G=zeros(10,1);
bin = floor(1+length(pred)*(0:bins)/bins);
for q = 1:bins
int = bin(q):(bin(q+1)-1); %bin indexes
N = length(int); % number of patients in bin
E=sum((pred(int))); % expected
O=sum(round(outcome(int))); % observed
Eprob = mean(pred(int)) ; % expected probability
G(q) = (E-O)^2 / (Eprob*N*(1-Eprob));
end
G = sum(G);
end
function [ c ] = cstat( pred, target )
%cstat Calculates the c-statistic
% [ c ] = cstat (pred, out)
% pred - vector containing the predicted targets
% target - vector containing the true targets
%
% c is equivalent to the AUROC, a measure of model discrimination
% Mathematically: Pr(pred|out==1 > pred|out==0)
%=== Arrange predictions
[pred,idx] = sort(pred,1,'ascend');
target=target(idx);
[N,P] = size(pred);
%=== Find location of negative targets
c=zeros(1,P); % 1xP where P is # of AUROCs to calculate
negative=false(N,P);
for n=1:P
negative(:,n) = target(:,n)==0;
%=== Count the number of negative targets below each element
negativeCS = cumsum(negative(:,n),1);
%=== Only get positive targets
pos = negativeCS(~negative(:,n));
c(n) = sum(pos);
end
count=sum(negative,1); %=== count number who are negative
count = count .* (N-count); %=== multiply by positives
c=c./count;
end
function [defStruct] = defaultStruct()
defStruct = struct('AUROC',[],...
'HL',[],...
'acc',[],...
'R',[],...
'B',[],...
'SMR',[],...
'coxA',[],...
'coxB',[],...
'TP',[],...
'FP',[],...
'TN',[],...
'FN',[],...
'sens',[],...
'spec',[],...
'ppv',[],...
'npv',[]);
end