
% Look at variables selected and estimated coeffs of them in top k
% models and compute the corresponding conditional variable inclusion probabilities
% Also looks at the overall model averaged fit and prediction and splits KM curves for
% exploratory visuals -- samples are partitioned by the median of the fitted 
% regression in th etraining sample, and the resulting two KM curves are plotted for 
% the training and (if any) validation data separately

% -----------
% 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
    % element 3+p:                - the posterior mean of the Weibull index parameter alpha in this model
    % elements (4+p):(4+2p)       - posterior mean of the regression parameter vector beta (includes intercept)
    % elements (5+2p):(8+6p+p*p)  - posterior variance matrix of (alpha,beta) (beta includes intercept) in vectorised form 
%

pmk=pm(1:k)/sum(pm(1:k));                           % condition on only the top k models chosen
X=std_rows(x);
Fit=zeros(n,k);  alpha=zeros(1,k);                  % to save fits and predictions 
pmp=zeros(N,1);                                     % to save marginal inclusive probs
figure(1); clf

for j=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);         	        % save posterior probs on predictors 
          alp=models(j,3+p);                        % posterior mean of Weibull index parameter
          b=models(j,4+p+(0:p))';                   % post mean of regn parameters
          A=[ones(n,1) X(ig,:)'];          	        % design matrix
      else  
          alp=models(j,3); b=models(j,4); A=[ones(n,1)]; 
      end
      Fit(:,j)=A*b;                                 % fitted & predicted linear regn
      alpha(j)=alp;   
      % predictions for hold-out validation cases are already in the above
      mu=exp(-Fit(:,j)); Fit(:,j)=(mu*log(2)).^(1/alp); % median predicted survival times 
      if (p>0) 
          thres=median(Fit(itrain,j)); 
          figure(1)
          km([y(itrain) 1-r(itrain) Fit(itrain,j)<=thres],'br','--')
          legend(char('<50%ile','>50%ile'),'location','SouthWest'); legend boxoff
          xlabel('Time');ylabel('Survival probability');
          title(['Training data - Model ',int2str(j),' : p=',int2str(p),' : Prob=',num2str(pm(j),2)])
          figure(2)
          km([y(ivalid) 1-r(ivalid) Fit(ivalid,j)<=thres],'br','--')
          xlabel('Months');ylabel('Survival probability');
          title(['Test data - Model ',int2str(j),' : p=',int2str(p),' : Prob=',num2str(pm(j),2)])
          legend(char('<50%ile','>50%ile'),'location','SouthWest'); legend boxoff      
      else,  'Null model - no graph' ; end    
      fprintf('Hit a key to continue to next model ....\n') 
pause
end

aveFit=Fit*pmk;
figure(1); 
thres=median(aveFit(itrain));
km([y(itrain) 1-r(itrain) aveFit(itrain)<=thres],'br','--')
xlabel('Months');ylabel('Survival probability');title(['Model average: Training cases'])
legend(char('<50%ile','>50%ile'),'location','SouthWest'); legend boxoff 
figure(2)
km([y(ivalid) 1-r(ivalid) aveFit(ivalid)<=thres],'br','--')
legend(char('<50%ile','>50%ile'),'location','SouthWest');legend boxoff
ylabel('Survival probability');title(['Model average: Test cases'])

 
