
# Look at variables selected and estimated coeffs of them in top k models, with
# predictions from each and the model average over top k models.  Also
# compute the corresponding conditional variable inclusion probabilities
# -- the latter in the vector pmp 

# -----------
# models array has following elements:  each row is one model of the 1:nm "top models"
# ordered in decreasing order of posterior probability. The number of columns is defined 
# by the largest model, and entries are NA/NaN for models smaller than the largest
# In each model/row the entries are: 
    # element  1:                 - dimension of the model = mumber of predictors p for this model
    # element  2:                 - log posterior probability of this model (the "score") 
    # elements 3:(2+p)            - the indices of the p variables in this model
    # elements (3+p):(3+2p)       - posterior mode of the regression parameter vector beta (includes intercept)
    # elements (4+2p):(4+4p+p*p)  - estimated posterior variance matrix of beta (includes intercept)
    #                               in vectorised form 
#

pmk <- pm[1:k]/sum(pm[1:k]);                        # condition on only the top k models chosen

train.m <- apply(x[,itrain],1,mean)                 # must only standardize based on the observations
train.sd <- sqrt(apply(x[,itrain],1,var))           # used to fit the models
X <- (t(x) - matrix(train.m,ncol=N,nrow=n,byrow=T)) # subtract the mean
X <- X/matrix(train.sd,ncol=N,nrow=n,byrow=T)       # divide the the sd                    

Fit <- matrix(0,nrow=n,ncol=k);                     # to save fits and predictions
pFit <- matrix(0,nrow=n,ncol=k);                     
pmp <- matrix(0,nrow=N,ncol=1);                     # to save marginal inclusive probs

for(j in 1:k)
{
      p <- models[j,1]                              # dim of this model
      if (p>0)
      {
          ig <- models[j,2+(1:p)]                   # predictors in model    
          pmp[ig] <- pmp[ig] + pmk[j];     	    # posterior probs on predictors 
          
          b <- matrix(as.numeric(models[j,(p+3):(2*p+3)]),ncol=1,nrow=p+1)
                                                    # post mode of regn parameters
                                                    
          A <- matrix(0,nrow=n,ncol=(p+1)); A[,1] <- 1;
          A[,-1] <- X[,ig]                          # design matrix
      }
      if(p==0)
      {
          b <- models[j,3]; A <- matrix(1,ncol=n,nrow=1); 
      }   
      Fit[,j] <- A%*%b;                       # fitted & predicted linear regn
      pFit[,j] <- 1 / (1+exp(-Fit[,j]));      #    ... and plug-in probabilities
      
      # predictions for hold-out validation cases are already in the above
      plot(1:n, pFit[,j], xlab="case id",ylab="Fitted Probabilities",ylim=c(0,1),type="n",main=paste("Model",j,": p =",p," : Prob =",round(pmk[j],4)),xlim=c(1,n + n*0.2),axes=F)
      box(); axis(2); tt <- axTicks(1); tt <- tt[tt<=n]; axis(1,at=tt)
      legend("bottomright", c("obs 0","obs 1","pred 0","pred 1","base"),pch=c(16,16,16,16,-1),lty=c(-1,-1,-1,-1,2),col=c("blue","red","cyan","magenta","black"),merge=T,bg="gray90")
      
      showtv(t(pFit), y, itrain, ivalid, j);
      abline(h=sum(y)/n,lty=2);                         
      par(ask=T)
}
   
avepFit <- pFit%*%pmk
plot(1:n, pFit[,j], xlab="case id",ylab="", ylim=c(0,1),type="n",main=paste("Model Average"),xlim=c(1,n + n*0.2),axes=F)
box(); axis(2); tt <- axTicks(1); tt <- tt[tt<=n]; axis(1,at=tt)
legend("bottomright", c("obs 0","obs 1","pred 0","pred 1","base"),pch=c(16,16,16,16,-1),lty=c(-1,-1,-1,-1,2),col=c("blue","red","cyan","magenta","black"),merge=T,bg="gray90")
showtv(t(avepFit),y,itrain,ivalid,1)
abline(h=sum(y)/n,lty=2)     

