#### Commands to select a 'best' model for logistic regression
###

#Say the output is stored in 'logitmodel'

## Fit the model to all the predictors in x

library(MASS)

x = as.data.frame(x)

logitmodel = glm(y ~ ., data=x, family=binomial(link="logit"))

bestmodel = stepAIC(logitmodel, trace=FALSE)	#trace=TRUE shows outputs on the screen as it
						# looks through various combinations

### You can explore other arguments in stepAIC

bestmodel$anova		#gives you the initial and final model and AIC in each step



## you can now use the bestmod to look other other quantities, R2, etc. and 
## predict at new values of x

## You can use predict.glm commands to predict at new points - similar to 
### thhe commands used in predict.lm

### RPSS  
ypred = predict(bestmodel,type="response")
acc1=y
N=length(acc1)
p1=length(acc1[acc1 == 0])/N
p2 = 1-p1

climo <- c(p1,p2)
rps1 <- rpsc <- 0
obs <- cbind(1-acc1,acc1)
pred = cbind(1-ypred,ypred)

for(i in 1:N){
   # Prediction model score
   rps1 <- rps1 + sum((cumsum(obs[i,])-cumsum(pred[i,]))^2)
   # Null model (climatology) score
   rpsc <- rpsc + sum((cumsum(obs[i,])-cumsum(climo))^2)
}

RPSS <- 1 - (rps1/rpsc)

### using verification package..
library(verification)
acctemp=acc1+1 # Add 1 to get categories 1 and 2
rps(acctemp, pred, baseline=climo )$rpss