# Evaluation of performance

## The scaled Brier score

Gary Weissman posted an excellent illustration of how we may calculate the maximum Brier score, which we need for the scaled Brier score. Three approaches were found to be equivalent, while the intuition for the maximum score may be best for the first formulation, as discussed at |Twitter.

### Code for scaled Brier score calculation

```brier_score <- function(obs, pred) { mean((obs - pred)^2) } # obs: 0/1 outcome y; pred: predicted probability p̂
scaled_brier_score_1 <- function(obs, pred) {
1 - (brier_score(obs, pred) / brier_score(obs, mean(obs))) } # mean(obs): ȳ
scaled_brier_score_2 <- function(obs, pred) {
1 - (brier_score(obs, pred) / (mean(obs) * (1 - mean(obs)))) }
scaled_brier_score_3 <- function(obs, pred) {
1 - (brier_score(obs, pred) / (mean(obs) * (1 - mean(obs))^2 + (1 - mean(obs)) * mean(obs)^2)) }
```

## Illustration of the Hosmer-Lemeshow goodness of fit test

### Code for HL test calculation

```# Jan 07, Ewout Steyerberg
# function to perform Hosmer-Lemeshow test for external validation
# instead of chi square and corresponding p-value, this function provides the number of subjects per group,
# and the mean values of p and y per group.
##########################
# p	: predicted probability
# Y	: outcome variable
# g	: number of groups to calculate H-L (10 is default)
#
# NB: the library Hmisc need to be attached in order to be able to run hl.ext2
##########################

hl.ext2<-function(p,y,g=10, df=g-1)
{
matres	<-matrix(NA,nrow=g,ncol=5)
sor	<-order(p)
p	<-p[sor]
y	<-y[sor]
groep	<-cut2(p,g=g)				#g more or less equal sized groups
len	<-tapply(y,groep,length)		#n per group
sump	<-tapply(p,groep,sum)			#expected per group
sumy	<-tapply(y,groep,sum)			#observed per group
meanp	<-tapply(p,groep,mean)			#mean probability per group
meany	<-tapply(y,groep,mean)			#mean observed per group
matres	<-cbind(len,meanp,meany, sump, sumy)	#matrix for results
contr<-((sumy-sump)^2)/(len*meanp*(1-meanp))	#contribution per group to chi square
chisqr<-sum(contr)				#chi square total
pval<-1-pchisq(chisqr,df)	#p-value corresponding to chi square with df degrees of freedom
cat("\nChi-square",chisqr," p-value", pval,"\n")
dimnames(matres)	<-list(c(1:g),Cs(n,avg(p),avg(y), Nexp, Nobs))
result  <- list(table(groep), matres,chisqr,pval)
}
```
additional/chapter15.txt · Last modified: 2019/05/03 08:50 by ewsteyerberg 