There are four models I'm looking at and trying to evaluate which will give me the best predictions going forward. So I loop through each model and look at the error rate, the Brier Error and the AIC. In some cases the AIC and Brier Error are not consistent (results for model 3 and 4 below). If they give contrasting indications, which one should be paid attention to? Why? Or am I looking at this the wrong way? For those interested, my R code and results are below:
line<- seq(71,79, by=.5) models<- list("target~total+tot_eit_h_h1+tot_both_h_h1","target~total", "target~total*spread*as.factor(loc)", "target~total+spread*as.factor(loc)") for(j in 1:length(models)){ for(i in 1:length(line)){ data$target <- ifelse(data$result>line[i], 1, ifelse(data$result<line[i], 0, rbinom(dim(data)[1],1,.5))) model<-glm(models[[j]], data=data, family="binomial") data$prediction <- predict(model, data, type='response') brier=with(data, mean((prediction-target)^2)) error_rate <- mean((data$prediction>0.5 & data$target==0) | (data$prediction<0.5 & data$target==1)) aics<-c(aics, model$aic) briers<-c(briers, brier) error_rates<-c(error_rates, error_rate) } print(paste(mean(error_rates), "error rates")) print(paste(mean(aics), "AIC")) print(paste(mean(briers), "BRIER")) cat("n") error_rates<-c() aics<-c() briers<-c() } RESULTS
[1] "0.381418853255588 error rates" [1] "1633.50190616553 AIC" [1] "0.232600213678421 BRIER" [1] "0.38491739552964 error rates" [1] "1641.54699258567 AIC" [1] "0.234847243715002 BRIER" [1] "0.314091350826045 error rates" [1] "1451.51883987259 AIC" [1] "0.200957523053587 BRIER" [1] "0.324509232264334 error rates" [1] "1467.88024134523 AIC" [1] "0.205168360087164 BRIER" Best Answer
You need to divide by $n$ to get the Brier score, and there should be no 100 in its formula. It is just the mean squared error. The "error rate" is not meaningful and will give misleading results, so I would ignore it. AIC is a measure of predictive discrimination whereas the Brier score is a combined measure of discrimination + calibration.
models can be a list of models, e.g.
models <- list(y ~ x1, y ~ x1 + x2, ...) What does line represent and should you really be fitting ordinal regression models?
Please update your code to show us the correct Brier score, along with AIC.