##
##  Moving plot that shows how an ROC curve traces out as the point of integration in the classifier changes.
##  Example simulated data is transitional logit.
##
##  jamesHonaker  2/3/07
##

rm(list=ls())

t<-1:5000
p<-runif(length(t))
q<-runif(length(t))

y<-seq(from=1,to=1,length=length(t))
t<-seq(from=1,to=1,length=length(t))
all.pi<-matrix(1,length(t),1)

for(i in 2:length(t)){
  xb<- -1 -p[i] + y[i-1]*(1 + q[i]) 
  current.pi<-1/(1+exp(-xb))
  y[i]<-runif(1)<current.pi
  all.pi[i,1]<-current.pi
  if(y[i]==y[i-1]){
    t[i]<- t[i-1]+1
  }
}  

lagy<-c(-9,y[1:(length(y)-1)])
lagt<-c(-9,t[1:(length(y)-1)])


dataset<-cbind(y,t,p,q,lagy,lagt)
dataset<-as.data.frame(dataset)

# Calculate survival times
flag<-t==1                                               # Flag denotes new regimes beginning
survival.times<-t[1:(length(t)-1)][flag[2:length(t)]]    # Notice we can further subset[] a subset[]


output.logit<-glm(y ~ (p + q),family=binomial,data=dataset[2:nrow(dataset),] )
output.trans<-glm(y ~ (p + q)*(lagy),family=binomial,data=dataset[2:nrow(dataset),] )

print(summary(output.logit))
print(summary(output.trans))

yhat.logit<-predict(output.logit,type="response")
yhat.trans<-predict(output.trans,type="response")

ytrim<-y[2:length(y)]


par(ask=FALSE)

n<-1000
history<-matrix(0,n+1,3)

for(i in 0:n){
  set.pi<-i/n

  history[i+1,1]<-i/n
  history[i+1,2]<-sum(yhat.trans>i/n & ytrim==1)/sum(ytrim==1)
  history[i+1,3]<-sum(yhat.trans>i/n & ytrim==0)/sum(ytrim==0) 

  par(mfrow=c(2,1))

  plot(density(yhat.trans[ytrim==0]),xlim=c(0,1),ylim=c(0,10),col="red" ,main="",xlab="density")
  rect(set.pi,0,1,10,col="grey")
  par(new=TRUE)  
  plot(density(yhat.trans[ytrim==0]),xlim=c(0,1),ylim=c(0,10),col="red" ,main="",xlab="density")
  par(new=TRUE)
  plot(density(yhat.trans[ytrim==1]),xlim=c(0,1),ylim=c(0,10),col="blue",main="",xlab="density")
  
  abline(v=set.pi)

  plot(history[1:(i+1),3],history[1:(i+1),2],type="l",ylim=c(0,1),xlim=c(0,1),ylab="True Positive Rate",xlab="False Positive Rate",main="ROC Curve",lwd=2,xaxs="i",yaxs="i")
  abseq<-(1:4)/5
  abline(v=abseq,lty=3,col="grey")
  abline(h=abseq,lty=3,col="grey")
  abline(a=0,b=1)
}