
# 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

# -----------
# 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):(2+2p)       - posterior mean of the regression parameter vector beta
    # elements (3+2p):(2+2p+p*p)  - posterior variance matrix of beta in vectorised form
    # final two elements:         - (s,d) - the residual SD estimate, and the posterior degrees of freedom 
#

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

my <- mean(y[itrain,]); sy <- sqrt(var(y[itrain,])); # again only stanardize based on training data
Y <- as.matrix((y - my)/sy);

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

tau <- 1; delta <- 3;
aa <- c(min(y),max(y))
eend <- dim(models)[2]
                     
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,2+p+(1:p)]),ncol=1,nrow=p)
                                                 # post mean of regn parameters
                                                    
       A <- as.matrix(X[,ig])                    # design matrix
   }
   if(p==0)
   {
       b <- 0; A <- 0; 
   }
   sig[1,j] <- models[j, eend-1]; df[1,j] <- models[j, eend];   
   Fit[,j] <- my + (A%*%b)*rep(sy,n);                   # fitted & predicted linear regn
      
   # predictions for hold-out validation cases are already in the above
   bb <- c(min(aa,Fit[,j]),max(aa,Fit[,j]))
   plot(bb, bb, ylab="data",xlab="Fitted/Predicted Values",type="n",main=paste("Model",j,": p =",p," : Prob =",round(pmk[j],4)),xlim=c(bb[1],bb[2]+0.05*(bb[2]-bb[1])))
   legend("bottomright", c("obs < mean","obs > mean","pred < mean","pred > mean"),pch=c(16,16,16,16),col=c("blue","red","cyan","magenta"),bg="gray90")
      
   if(p>0)
   {
      scattertv(Fit[,j], y, Y>0, itrain, ivalid)
   }
   if(p==0)
   {
      ii <- rep(1,n); ii[ivalid] <- 0; show(t(y),ii,1);
   }
   abline(0,1,lty=2);                         
   par(ask=T)
}
   
aveFit <- Fit%*%pmk
bb <- c(min(aa,aveFit),max(aa,aveFit))
plot(bb, bb, xlab="Averaged Fit/Prediction",ylab="data", xlim=c(bb[1],bb[2]+0.05*(bb[2]-bb[1])),type="n",main=paste("Model Average"))
legend("bottomright", c("obs < mean","obs > mean","pred < mean","pred > mean"),pch=c(16,16,16,16),col=c("blue","red","cyan","magenta"),bg="gray90")
scattertv(aveFit,y,Y>0,itrain,ivalid)
abline(0,1,lty=2)     

