###################################### # Bayesian Model Averaging Program # ###################################### # This version: 2008-11-05 # Martin Feldkircher # martin.feldkircher@gzpace.net, http:feldkircher.gzpace.net # Stefan Zeugner # ##################### # The main code starts at line 371 with the function "fls=function(....)" and is written # by Martin Feldkircher and Stefan Zeugner as part of their work at the Institute for Advanced Studies (IHS), # Vienna in 2006/2007. Descriptions of the algorithms and priors used can be found in # Gary Koop ("Bayesian Econometrics", Wiley & Sons), Fernandez, C., E. Ley and M.F.J. Steel (2001b) # “Model Uncertainty in CrossCountry Growth Regressions,” Journal of Applied Econometrics and # and Fernandez, C., E. Ley and M.F.J. Steel (2001a) “Benchmark Priors for Bayesian Model Averaging,” # Journal of Econometrics, 100: 381–427. #################### # USAGE # ################################################################################################### # fls(X.data,burn=100,iter=1000,g="bric",start.value=41,theta="random", prior.msize=7,mcmc="bd",nmodel=100, # logfile=FALSE,logstep=100000,beta.save=TRUE,exact=FALSE,printRes=TRUE,int=FALSE,ask.set=FALSE) # # #################### # PARAMETERS # # # ################################################################################################### #X.data submit a data frame or a matrix, where the first column corresponds to the dependent variable # followed by the covariates, including a constant term is not necessary since y and X is # demeaned automatically. #burn is the number of burn-in draws #iter is the number of posterior draws #g is the hyperparameter on Zellner's g-prior for the regression coefficients. You can specify # g="bric" corresponding to the benchmark prior suggestion from FLS (2001), i.e g=max(1/N, 1/K^2) # with K denoting the total number of regressors and N the number of observations #start.value specifies the starting model. You can choose either a specific model by the corresponding # column indices (e.g. starting.model=numeric(ncol(X)) starts from the null model including # solely a constant term) or you set a number (e.g. starting.model=20). In the latter case # randomly 20 covariates are chosen and the starting model is identified by those regressors # with a t-statistics>0.2. # start.value=0 or start.value=NULL starts from the null model # start.value=NA sets start.value=min(N-1,K) #theta regards the prior on model size. It can be either "random" or "fixed" and is based on the # working paper "On the Effect of Prior Assumptions in Bayesian Model Averaging with Applications # to Growth Regression", by Ley and Steel (2008). Theta denotes the a priori inclusion probability # of a regressor. Their suggestion is to use a binomial-beta hyperprior on theta (i.e. theta=random) # in order to be noninformative on model size. You can use theta=fix if you have strong prior # information on model size. #prior.msize corresponds to the expected value of the model size prior. For theta=random there is little # impact on results by varying prior.msize. For fixed theta (i.e. informative) prior this is a # sensible choice. #mcmc default is "bd" which corresponds to a birth / death MCMC alogrithm. You can choose # also "revjump" where we have implemented a true reversible jump algorithm where we have # added a "move" step to the birth / death steps from "bd". #nmodel is the number of best models for which information is stored. Convergence analysis of # the sampler by means of the correlation between analytic posterior model probabilities # (pmp's) and hat of the MCMC sampler is based on the number you have set in nmodels. Also # if you want to save the regression coefficients (beta.save=T), they are taken from the # nmodel best models. Setting nmodel500 slows down the MCMC sampler. Note that posterior # inclusion probabilites, and mean calculations are based on the MCMC frequencies as opposed # to exact analytical calculations (as in fls). # Set nmodel=0 if you do not want to save this information. #beta.save can be either TRUE or FALSE, and specifies whether you want to save the regression coefficients # for the "nmodel" best models. #logfile setting logfile=TRUE produces a logfile named "fls.log" in your current working directory, # in order to keep track of the sampling procedure. # setting logfile equal to some filepath (like "subfolder/bla.txt") puts the logfile # into that specified position. #logstep specifies at which number of posterior draws information is written to the log file #printRes print out results to command line after ending the routine #ask.set TRUE or FALSE whether R should pause between the graphs drawn at the end of the routine (only matters in interactive mode) ############################################################################################## #Output # # # ################################################################################################### #estimates a matrix....PIP denotint Posterior Inclusion Probabilities based on the relative frequencies # of the MCMC sampler, post mean = posterior mean, post SD = posterior standard deviation # and Idx is the column index of the specific regressor in the submitted data matrix #info contains information about the "Mean nr. of Regressors" (not counting the constant term), # "Draws"=posterior draws, "Burnins"=nr. of burnins taken, "Time" denotes total time elapsed # since calling the "fls" function, "Nr. of models visited" counts each time a model is accepted. # Note that we do not take into account the case of revisiting models by the sampler. # Modelspace is simply indicating the whole model space (2^K) and percentage visited is # the nr. of models visited as a percentage of 2^K. "Corr PMP" is the correlation between # analytic and MCMC posterior model probabilites, where a correlation of 0.99 indicates # excellent convergence. For "nmodel=100" the best 100 models are considered for the correlation # analysis. Finally, Nr. of Observations is given in the output as well. # # -----The following objects are saved in your output object------- #topmodels matrix of dimension K x nmodel and consists of the nmodel best models where zero stands for # exclusion of a regressor and 1 for inclusion #k.vec is a vector of length K+1 and indicates how many model werde drawn with a size of k.vec[i] regressors, # where the first entry in k.vec corresponds to the nullmodel (hence K+1 entries) #start.pos is the chosen starting model #beta.draws matrix of stored regression coefficients if specified via beta.save=TRUE #topmodprob analytical posterior model probabilities of the best "nmodel" models #pmp.10 Matrix where the first column represents analytical and the second MCMC model probabilities # they are normalized such that each column sums up to 1. #topmodobject function to save pmp and betas #theta theta as specified by the user #K nr. of regressors #prior.msize prior model size as specified by the user # The plot of the model size shows the prior as compared to the posterior model size and is done automatically after the # routine finishes. To redo the plot use modelsize.plot(X=result,line.width=2,sub="This is my plot"), # where comes from result=fls(......................) #fls.call the function call as it was entered into the command line ###### FLS MAIN FUNCTION ######################### fls <-function(X.data,burn=100,iter=1000,g="bric",start.value=41,theta="random", prior.msize=7, mcmc="bd",nmodel=100,logfile=FALSE,logstep=100000,beta.save=TRUE,printRes=T,ask.set=F){ # AUXILIARY FUNCTIONS TO FUNCTION "fls()"############## #START: SUBFUNCTIONS ################ # function calculates the likelihood lprob=function(positions,g0,g2,yty,k,N,null.lik,K=K,XtX.big=XtX.big,Xty.big=Xty.big,...){ lg02=log(g0*g2); Nlg2=(N-1)*log(g2) # if nullmodel if(sum(k)==0){ lprob=null.lik;b1new=rep(0,K);stdevnew=rep(0,K); Xty=0;XtX=0;XtXinv=0;bhat=0;ymy=yty;positions=1; } else { XtX<-XtX.big[positions,positions] Xty<-Xty.big[positions] #do cholesky split: A=XtX=LL' #get lower triangular matrix from cholesky split L=chol(XtX) XtXinv<-chol2inv(L) bhat<-crossprod(XtXinv,Xty) ymy=yty-as.vector(crossprod(Xty,bhat)) lprob = .5*(k*lg02-(N-1)*log(ymy + g0*yty)-Nlg2)#.5*(k*log(g0*g2)-(N-1)*log((ymy + g0*yty)*g2))# b1new = g2*bhat bcov=((g0*yty+ymy)*(g2^2)/(N-2))*XtXinv #vs2=(g0*yty+ymy)*g2, stdevnew = diag(bcov)+b1new^2 } return(list(lprob=lprob,b1new=b1new,stdevnew=stdevnew, child.lprob = function(addix=0,dropix=0,...) { if (!any(as.logical(c(addix,dropix)))) {return(lprob)} if (all(as.logical(c(addix,dropix)))) { #swap jhere=(1:k)[positions==dropix]; poshere=positions[-jhere];Xj=XtXinv[,jhere];Xtxi=XtX.big[poshere,addix] bxlessj=crossprod(XtXinv,XtX.big[positions,addix])-Xj*XtX.big[addix,dropix]; bhatx=bxlessj[-jhere]-Xj[-jhere]*bxlessj[jhere]/Xj[jhere] child.ymy = as.vector(ymy+bhat[jhere]^2/Xj[jhere]-(Xty.big[addix]-crossprod(Xty.big[poshere],bhatx))^2/(XtX.big[addix,addix]-crossprod(bhatx,Xtxi))) return(.5*(k*lg02-(N-1)*log(child.ymy + g0*yty)-Nlg2)) } else { if (addix==0) { #drop jhere=(1:k)[positions==dropix] child.ymy=ymy+bhat[jhere]^2/XtXinv[jhere,jhere] return(.5*((k-1)*lg02-(N-1)*log(child.ymy + g0*yty)-Nlg2)) } else { #add Xtxi=XtX.big[positions,addix] bhatx=as.vector(crossprod(XtXinv,Xtxi)) child.ymy = as.vector(ymy - (Xty.big[addix]-crossprod(bhatx,Xty))^2 /(XtX.big[addix,addix]-crossprod(bhatx,Xtxi))) return(.5*((k+1)*lg02-(N-1)*log(child.ymy + g0*yty)-Nlg2)) } } })) } # function calculates the posterior odds pmp.ran= function(ki=knew,K=K,m=prior.msize, lprobnew=lprobnew){ post.odds1=lprobnew+lgamma(1+ki)+lgamma((K-m)/m+K-ki) return(post.odds1) } pmp.fix= function(ki=knew,K=K,m=prior.msize, lprobnew=lprobnew){ post.odds1=lprobnew+ki*log((m/K))+(K-ki)*log(1-m/K) return(post.odds1) } pmp.unif= function(ki=knew,K=K,m=prior.msize, lprobnew=lprobnew){ return(lprobnew) } #"starter" draws random start vector of size start.value, and keeps those regressors with a t-stat>0.2 starter=function(K,start.value,y,N=N,XtX.big=XtX.big,Xty.big=Xty.big,X=X){ # in case user submitted a single number as start.value # we randomly draw a model with start.value regressors # finally only regressors with t-stats>0.2 are kept and build the starting model #some input checks if (length(start.value)==0) {start.value=numeric(K)} if (length(start.value)==1) {if (start.value==0) {start.value=numeric(K)}} if (is.na(start.value[1])) {start.value=min((N-2),K)} if(length(start.value)==1){ if (start.value>min((N-2),K)){ cat("Submitted Start value is too large, used\n min(N-2,K) as starting model size instead\n\n") start.value=min((N-2),K) } # draw randomly sorter=rnorm(K) sequ<-seq(1:K) start.position=order(sorter,sequ)[1:start.value] # calculate bhats and t-stats XtX.start<-XtX.big[start.position,start.position] XtXinv.start<-chol2inv(chol(XtX.start)) bhat=XtXinv.start%*%Xty.big[start.position] e=y-X[,start.position]%*%bhat sse=crossprod(e) s2=as.numeric(sse/(N-length(start.position))) bcov=s2*XtXinv.start bt=bhat/sqrt(diag(bcov)) # choose only regressors with t-stat>0.2 molddraw=rep(0,K) goodguy=as.numeric(abs(bt)>.2) molddraw[start.position]=goodguy return(list(molddraw=molddraw,bhat=bhat,start.position=start.position)) } # else we start with the user specified starting model # in this case we do not test whether the t-stat's are greater than 0.2 # but start with exactly the selected model no matter whether it is a # good starting model or not # in case you want to start with the null model if(length(start.value)>1 && sum(start.value)==0){ return(list(molddraw=rep(0,K),bhat=rep(0,K),start.position=integer(0))) } if(length(start.value)>1 && sum(start.value)>0){ if(length(start.value)!=K){ # in case user specified a too big model stop("Starting Model contains unequal to K regressors,please respecify") } start.position=which(start.value==1) XtX.start<-XtX.big[start.position,start.position] XtXinv.start<-chol2inv(chol(XtX.start)) bhat=XtXinv.start%*%Xty.big[start.position] molddraw=rep(0,K); molddraw[start.position]=1 return(list(molddraw=molddraw,bhat=bhat,start.position=start.position)) } } ########################################################################################################################### #Sample Functions ############################################################################################################################ #First, we have implemented the original FLS Sample Function; here a variable is drawn from the set of # K regressors, then #conditional on whether it is included in the current model it can be discarded or added. fls.samp=function(molddraw=molddraw,K=K,mMinus=NA,mPlus=NA){ indch=ceiling(runif(1)*(K)) #rounding to the smallest integer part by floor, uniform distr. [0,1] if (molddraw[indch]==1){ #dropping addvar=0;dropvar=indch; mnewdraw=molddraw mnewdraw[indch]=0 positionnew=(1:K)[mnewdraw==1] } else { #adding addvar=indch;dropvar=0; mnewdraw=molddraw mnewdraw[indch]=1 positionnew=(1:K)[mnewdraw==1] } return(list(mnewdraw=mnewdraw,positionnew=positionnew,addi=addvar,dropi=dropvar,indch=indch)) } ############################################################################################################################# #Second, we have implemented a reversible jump algorithm, where we #have added a move step. See below ############################################################################################################################## #Reversible Jump Algorithm rev.jump=function(molddraw=molddraw,K=K,mMinus=NA,mPlus=NA){ rev.idx=ceiling(runif(1,0,2)) #rev.idx is a flag that indicates the three possible steps of #the reversible jump algorithm, 1=birth or death and 2=move. # Perform Death, Birth or Move Step # if rev.idx is 1, do the same as in fls sampler (i.e. increase or decrease depending # on variables already included in the model if(rev.idx==1){ birth.death=fls.samp(molddraw=molddraw,K=K) mnewdraw=birth.death$mnewdraw positionnew=birth.death$positionnew addvar=birth.death$addi; dropvar=birth.death$dropi } #move step if(rev.idx==2){ var.in=(1:K)[as.logical(molddraw)] #positions of the variables that are currently in the model var.out=(1:K)[!as.logical(molddraw)] #positions of the variables that are currently out of the model var.in.rand=ceiling(length(var.in)*runif(1)); addvar=var.out[ceiling(length(var.out)*runif(1))] dropvar=var.in[var.in.rand] positionnew=c(var.in[-var.in.rand],addvar) mnewdraw=molddraw; mnewdraw[addvar]=1; mnewdraw[dropvar]=0; #mnewdraw[positionnew]=1 } return(list(mnewdraw=mnewdraw,positionnew=positionnew,addi=addvar,dropi=dropvar,indch=indch)) } ########################################################################################################################### #Sample Functions ############################################################################################################################ #First, we have implemented the original FLS Sample Function; here a variable is drawn from the set of # K regressors, then #conditional on whether it is included in the current model it can be discarded or added. fls.samp.int=function(molddraw=molddraw,K=K,mPlus=mPlus,mMinus=mMinus){ indch=ceiling(runif(1)*(K)) #rounding to the smallest integer part by floor, uniform distr. [0,1] # have to make sure that we delete all the regressors if (molddraw[indch]==1){ #dropping mnewdraw=as.numeric(molddraw>mMinus[indch,]) #mnewdraw=molddraw-mMinus[indch,] #mnewdraw[which(mnewdraw<0)]=0 dropvar = (1:K)[xor(molddraw,mnewdraw)]; addvar=0; } else{ #adding mnewdraw=as.numeric(molddraw|mPlus[indch,]) #mnewdraw=molddraw+mPlus[indch,] #mnewdraw[which(mnewdraw>1)]=1 addvar = (1:K)[xor(molddraw,mnewdraw)]; dropvar=0; } positionnew=which(mnewdraw==1) return(list(mnewdraw=mnewdraw,positionnew=positionnew,addi=addvar,dropi=dropvar,indch=indch)) } ############################################################################################################################# #Second, we have implemented a reversible jump algorithm, where we #have added a move step. See below ############################################################################################################################## #Reversible Jump Algorithm rev.jump.int=function(molddraw=molddraw,K=K,mPlus=mPlus,mMinus=mMinus){ rev.idx=ceiling(runif(1,0,2)) #rev.idx is a flag that indicates the three possible steps of #the reversible jump algorithm, 1=birth or death and 2=move. # Perform Death, Birth or Move Step # if rev.idx is 1, do the same as in fls sampler (i.e. increase or decrease depending # on variables already included in the model if(rev.idx==1){ birth.death=fls.samp.int(molddraw=molddraw,K=K,mPlus=mPlus,mMinus=mMinus) mnewdraw=birth.death$mnewdraw positionnew=birth.death$positionnew addvar=birth.death$addi; dropvar=birth.death$dropi } #move step if(rev.idx==2){ var.in=(1:K)[as.logical(molddraw)] #positions of the variables that are currently in the model var.out=(1:K)[!as.logical(molddraw)] #positions of the variables that are currently out of the model mnewdraw=(molddraw>mMinus[var.in[ceiling(length(var.in)*runif(1))],]) #mnewdraw=molddraw-mMinus[var.in[var.in.rand],] #mnewdraw[which(mnewdraw<0)]=0 mnewdraw=mnewdraw|mPlus[var.out[ceiling(length(var.out)*runif(1))],] #mnewdraw=mnewdraw+mPlus[var.out[ceiling(length(var.out)*runif(1))],] #mnewdraw[which(mnewdraw>1)]=1 positionnew=(1:K)[mnewdraw] addvar = (1:K)[molddrawmnewdraw]; } return(list(mnewdraw=as.numeric(mnewdraw),positionnew=positionnew,addi=addvar,dropi=dropvar,rev.idx=rev.idx)) } ############################################################################################## #Function that does some posterior calculations ############################################################################################## post.calc=function(inccount,iter,msize,b1mo,b2mo,y.mean,X.mean,k.vec,null.count,models.visited,nmodel=nmodel,spatFilt=F, beta.save,topmods,X.data,burn,timed,theta=theta,m=prior.msize,exact=exact,SAR=F,rho.vec=NULL){ K=ncol(X.data)-1 N=nrow(X.data) beta.draws=NA # posterior inclusion probabilities pip=inccount/iter # mean number of regressors in models' nr.reg=msize/iter b1mo=b1mo/iter b2mo=b2mo/iter bsd=sqrt(b2mo-b1mo^2) # retrieve constant term cons=y.mean-b1mo%*%X.mean # adjust k.vec to include the nullmodel k.vec=c(null.count,k.vec) # Calculating Posterior Model Probabilities model.space=2^K fraction.model=models.visited/model.space*100 fraction.topmodel=sum(topmods$ncount())/iter*100 # this is the number of visits the topmodels accounted for relative to if (fraction.topmodel==0) {fraction.topmodel=NA} # the number of iterations topmod = topmods$bool_binary() if(is.null(colnames(X.data)[-1]) || colnames(X.data)[-1]==""){ reg.names<-paste("beta",1:K) colnames(X.data)=c("y",reg.names) } else{ #reg.names=unlist(dimnames(X.data[,-1])[2]) reg.names=colnames(X) } colnames(topmod)<-topmods$bool() if(SAR){ rownames(topmod)=c(reg.names,"rho") # Prob of top "nmodel" models, analytical (after deleting rest of models)' #LeSage Integration likvec=topmods$fixed_vector() lt1=trap(lmarginal=likvec,N=detlength,gridsize=gridsize)$probs } if(spatFilt){ topmod=topmod[-nrow(topmod),] if(length(reg.names)==nrow(topmod)){ rownames(topmod)=reg.names } colnames(topmod)<-topmods$bool() #c(paste("Model",1:nmodel)) lt1=topmods$lik() - max(topmods$lik()) # do this to get only positive probabilities lt1=exp(lt1)/sum(exp(lt1)) } if(!spatFilt & !SAR){ colnames(topmod)<-topmods$bool() #c(paste("Model",1:nmodel)) rownames(topmod)=reg.names # Prob of top "nmodel" models, analytical (after deleting rest of models)' lt1=topmods$lik() - max(topmods$lik()) # do this to get only positive probabilities lt1=exp(lt1)/sum(exp(lt1)) } # if (topmods$nbmodels>0) { # if(colSums(topmod[1:K,])[ncol(topmod)]==0){ # topmod=topmod[,-ncol(topmod)] # } # } #MCMC Prob of top "nmodel" models, numerical (after deleting rest of models)' lt2=topmods$ncount()/sum(topmods$ncount()) cpoint=min(c(ncol(topmod),length(lt1),length(lt2))) # do this because last value is initialized one #rbind the probs to the topmodmatrix topmod=rbind(topmod[,1:cpoint],lt1[1:cpoint], lt2[1:cpoint]) rownames(topmod)[(nrow(topmod)-1):nrow(topmod)]=c("PMP (Exact)","PMP (MCMC)") #do correlation test pmp.10<-cbind(lt1,lt2) corr.pmb<-try(cor(lt1,lt2),silent=T) if(class(corr.pmb)=="try-error"){ corr.pmb=NA } colnames(pmp.10)<-c("PMP Analytical", "PMP MCMC") prior=paste(theta, "/", m) # Info Object info<-as.character(c(format(round(nr.reg,4),nsmall=4),format(iter,nsmall=0),format(burn,nsmall=0), format(timed,nsmall=4),models.visited,format(model.space,digits=2), format(fraction.model,digits=2),format(fraction.topmodel,digits=2),format(round(corr.pmb,4),nsmall=4), format(N,nsmall=4),prior)) names(info)<-c("Mean nr. of Regressors", "Draws","Burnins", "Time", "Nr. of models visited", "Modelspace 2^K", "Percentage visited","Percentage Topmodels","Corr PMP","Nr.of Obs.", "Model Prior") # store information for the correlation test for different chains models=cbind(topmods$bool(),topmods$ncount()) if(exact){ exact.pmp=pmp.10[,1] pip=c() for(i in 1:(nrow(topmod)-2)){ pip[i]=sum(exact.pmp[which(topmod[i,]==1)]) } betas=topmods$betas() betas.bma=rowSums(betas*matrix(exact.pmp,nrow=K,ncol=ncol(betas),byrow=T)) #multiply betas with corr. pmp's and sum up stdev.bma=sqrt(rowSums(betas^2*matrix(exact.pmp,nrow=K,ncol=ncol(betas),byrow=T))- rowSums(betas*matrix(exact.pmp,nrow=K,ncol=ncol(betas),byrow=T))^2) #set in main script, if(exact==T){beta.save=T,nmodel=min(100,nmodel)} idx=2:(length(pip)+1) post.mean<-cbind(pip,betas.bma,stdev.bma,idx) # rownames(post.mean)<-reg.names colnames(post.mean)<-c("PIP","Post Mean", "Post SD","Idx") post.mean<-post.mean[order(-post.mean[,1]),] #order the results table according to PIP } else{ # Posterior mean and stand dev of each coefficient' idx=2:(length(pip)+1) post.mean<-cbind(pip,b1mo,bsd,idx) # rownames(post.mean)<-reg.names colnames(post.mean)<-c("PIP","Post Mean", "Post SD","Idx") post.mean<-post.mean[order(-post.mean[,1]),] #order the results table according to PIP } # have to do some reordering cause of the final integration with the saved liks if(SAR){ topmod=topmod[,order(-topmod["PMP (Exact)",])] } # assign names to the saved betas if(beta.save & !SAR){ beta.draws=as.matrix(topmods$betas()) #kick out the initialized zero column if nr. of models0){ colnames(beta.draws)=beta.names[-c(which(beta.names=="0"))] } } return(list(post.mean=post.mean,info=info,topmod=topmod,models=models, pmp.10=pmp.10,cons=cons,k.vec=k.vec,reg.names=reg.names,beta.draws=beta.draws)) } combine_chains <- function(...) { #to combine outputs of the the function fls #EXAMPLE: #bma1<-fls(X.data=t5.within[,1:20],burn=100,iter=1000,g=TRUE,nmodel=10,logfile=TRUE,beta.save=FALSE,start.value=41,step=1000) #bma2<-fls(X.data=X.data,burn=1000,iter=10000,g=TRUE,nmodel=10,logfile=TRUE,beta.save=FALSE,start.value=41,step=1000) #out=combine_chains(bma1,bma2) combine_2chains <- function(flsobject1,flsobject2) { if (nrow(flsobject1$estimates)!=nrow(flsobject2$estimates)) { stop("Both chains need to have the same number (and order) of regressors!") } if ("beta.draws"%in%names(flsobject1) & "beta.draws"%in%names(flsobject2)) { beta.save = TRUE } else { beta.save = FALSE } topmodobj=combine_topmods(flsobject1$topmodobject,flsobject2$topmodobject) iter1=flsobject1$info[2] iter2=flsobject2$info[2] msize1=flsobject1$info[1]*iter1 msize2=flsobject2$info[1]*iter2 inccount1=flsobject1$estimates[,1]*iter1 inccount2=flsobject2$estimates[,1]*iter2 b1mo1=flsobject1$estimates[,2]*iter1 b1mo2=flsobject2$estimates[,2]*iter2 b2mo1=(flsobject1$estimates[,3]^2+flsobject1$estimates[,2]^2)*iter1 b2mo2=(flsobject2$estimates[,3]^2+flsobject2$estimates[,2]^2)*iter2 model.space=flsobject1$info[5] models.visited1=flsobject1$info[4] models.visited2=flsobject2$info[4] iter=iter1+iter2 msize=msize1+msize2 inccount=inccount1+inccount2 b1mo=b1mo1+b1mo2 b2mo=b2mo1+b2mo2 models.visited=models.visited1+models.visited2 timed=flsobject1$info[3]+flsobject2$info[3] #do it as der_fliega does pip=inccount/iter nr.reg=msize/iter b1mo=b1mo/iter b2mo=b2mo/iter bsd=sqrt(b2mo-b1mo^2) #Calculating Posterior Model Probabilities fraction.model=models.visited/model.space*100 top10mod = topmodobj$bool_binary() colnames(top10mod)<-topmodobj$bool() rownames(top10mod)<-rownames(b1mo1) #Prob of top 10 models, analytical (after deleting rest of models)' lt1=topmodobj$lik() - max(topmodobj$lik()) lt1=exp(lt1)/sum(exp(lt1)) #Prob of top 10 models, numerical (after deleting rest of models)' lt2=topmodobj$ncount()/sum(topmodobj$ncount()) pmp.10<-cbind(lt1,lt2) corr.pmb<-cor(lt1,lt2) colnames(pmp.10)<-c("PMP Analytical", "PMP MCMC") #Prob of top 10 models out of total number of models, numerical' topmodprob=topmodobj$ncount()/iter #Info Object info<-c(nr.reg, iter,timed,models.visited,model.space,fraction.model,corr.pmb) names(info)<-c("Mean nr. of Regressors", "Draws", "Time", "Nr. of models visited","Modelspace 2^K", "Percentage visited","Corr PMP") #store information for the correlation test for different chains models=cbind(topmodobj$bool(),topmodobj$ncount()) #Posterior mean and stand dev of each coefficient' post.mean<-cbind(pip,b1mo,bsd) #maybe order? rownames(post.mean)<-rownames(flsobject1$estimates) colnames(post.mean)<-c("PIP","Post Mean", "Post SD") stpos1=as.matrix(flsobject1$start.pos);stpos2=as.matrix(flsobject2$start.pos) start.position=cbind(rbind(stpos1,matrix(0,max(0,nrow(stpos2)-nrow(stpos1)),ncol(stpos1))),rbind(stpos2,matrix(0,max(0,nrow(stpos1)-nrow(stpos2)),ncol(stpos2)))) if(beta.save){ return(list(estimates=post.mean,time=timed,topmodels=top10mod,info=info,models=models,start.pos=start.position,beta.draws=topmodobj$betas(),topmodprob=topmodprob,pmp.10=pmp.10,print(post.mean),print(info),topmodobject=topmodobj)) } return(list(estimates=post.mean,time=timed,topmodels=top10mod,info=info,models=models,start.pos=start.position,topmodprob=topmodprob,pmp.10=pmp.10,print(post.mean),print(info),topmodobject=topmodobj)) } ############################################################################################################# #this is the rest of the combine function; the combine function is iteratively used to combine as many chains #as are specified by (...) ############################################################################################################ arglist=list(...) combined_output <- combine_2chains(arglist[[1]],arglist[[2]]) if (nargs()>2) { for (inarg in 3:nargs()) { combined_output <- combine_2chains(arglist[[inarg]],combined_output) } } ############################################################################################################ return(combined_output) } combine_topmods <- function(topmodobj1,topmodobj2) { #to combine top10 models objects of function fls() #e.g. ppp=combine_topmods(test1$topmodobject,test2$topmodobject) #retrieve the necessary properties nregs1=topmodobj1$nregs nregs2=topmodobj2$nregs if (nregs1!=nregs2) {stop("The number of regressors in both BMA chains has to be the same!")} k1=length(topmodobj1$ncount()) k2=length(topmodobj2$ncount()) nbmodels1=topmodobj1$nbmodels nbmodels2=topmodobj2$nbmodels ncount1=topmodobj1$ncount() ncount2=topmodobj2$ncount() lik1=topmodobj1$lik() lik2=topmodobj2$lik() bool1=topmodobj1$bool() bool2=topmodobj2$bool() betas1=topmodobj1$betas() betas2=topmodobj2$betas() if (all(betas1==0)|all(betas2==0)) {dobetas=F} else {dobetas=T} #first look which models of 1 are already in 2 and #for these just update the ncounts (note: this is quite easy, since this subset in 1 has the same order as in 2) idxin2_boolof1in2=match(bool1,bool2) idxin1_boolof1in2=which(!is.na(idxin2_boolof1in2)) idxin2_boolof1in2=idxin2_boolof1in2[!is.na(idxin2_boolof1in2)] ncount2[idxin2_boolof1in2]=ncount2[idxin2_boolof1in2]+ncount1[idxin1_boolof1in2] #strip 1 of all the models that were already in 2 ncount1=ncount1[-idxin1_boolof1in2] lik1=lik1[-idxin1_boolof1in2] bool1=bool1[-idxin1_boolof1in2] if (dobetas) {betas1=betas1[,-idxin1_boolof1in2]} #now do A u (B\(AnB)) lika=c(lik2,lik1) orderlika=order(lika,decreasing=TRUE) lika=lika[orderlika] ncounta=c(ncount2,ncount1)[orderlika] boola=c(bool2,bool1)[orderlika] if (dobetas) { betasa=cbind(betas2,betas1)[,orderlika] betasa_not0=betasa!=0 vecka=colSums(betasa_not0) vbetaa=as.vector(betasa); vbetaa=vbetaa[betasa_not0] } else { vecka=0;vbetaa=0 } return(top10(nregs1,nbmodels1+nbmodels2,TRUE,inivec_lik=lika,inivec_bool=boola,inivec_count=ncounta,inivec_vbeta=vbetaa,inivec_veck=vecka)) } hexcode.binvec.convert <- function(length.of.binvec) { #function to initialise conversion betwwen logical vector (such as c(1,0,0,0)) and character hexcode vector (such as "f") #length.of.binvec is the desired length of the inserted and resulting logical vectors; #the initialisation will fit some leading zeros to make it convertible into hexcode (length of bin. vector as a multiple of 4) if (length(length.of.binvec)>1) length.of.binvec=length(length.of.binvec) addpositions=4-length.of.binvec%%4; positionsby4=(length.of.binvec+addpositions)/4; hexvec=c(0:9,"a","b","c","d","e","f"); #lookup list for converting from binary to hexadecimal hexcodelist=list("0"=numeric(4),"1"=c(0,0,0,1),"2"=c(0,0,1,0),"3"=c(0,0,1,1),"4"=c(0,1,0,0),"5"=c(0,1,0,1),"6"=c(0,1,1,0), "7"=c(0,1,1,1),"8"=c(1,0,0,0),"9"=c(1,0,0,1),"a"=c(1,0,1,0),"b"=c(1,0,1,1),"c"=c(1,1,0,0),"d"=c(1,1,0,1),"e"=c(1,1,1,0),"f"=c(1,1,1,1)); #lookup list for converting from hexadecimal to binary return(list( as.hexcode = function(binvec) { #convert logical vector to hexcode character incl=c(numeric(addpositions),binvec);dim(incl)=c(4,positionsby4); #split into elements of four positions return(paste(hexvec[crossprod(incl,2L^(3:0))+1],collapse="")); }, as.binvec = function(hexcode) { #convert hexcode character to numeric vector (e.g. "a" to c(1,0,1,0)) return(unlist(hexcodelist[unlist(strsplit(hexcode,"",fixed=T),recursive=F,use.names=F)],recursive=F,use.names=F)[-(1:addpositions)]) })) } hex2bin<-function(hexcode) { #user-friendly function to convert some hexcode character to numeric vector (e.g. "a" to c(1,0,1,0)) if (!is.character(hexcode)) stop("please input a character like '0af34c'"); hexobj<-hexcode.binvec.convert(length(hexcode)*16L) return(hexobj$as.binvec(hexcode)) } bin2hex<-function(binvec) { #user-friendly function to convert some logical vector to hexcode character(e.g. c(1,0,1,0) or c(T,F,T,F) to "a") if (!is.logical(binvec)) {if (is.numeric(binvec)) {binvec=as.logical(binvec)} else {stop("need to supply a logical vector like c(T,F) or c(1,0)")} } hexobj<-hexcode.binvec.convert(length(binvec)) return(hexobj$as.hexcode(binvec)) } top10=function(nmaxregressors=10,nbmodels=10,bbeta=FALSE,lengthfixedvec=0,..., inivec_lik=-10e10,inivec_bool=NA,inivec_count=0,inivec_vbeta=0,inivec_veck=0){ #object used by fls to save the best models #set up the variables to be augmented in the process #use top10(...) to initalise this object hexobject<-hexcode.binvec.convert(nmaxregressors) #initialize object for hexcode to logical vector conversion if (is.na(inivec_bool[1])) {inivec_bool=numeric(nmaxregressors)} #initialise varianles top10_lik=inivec_lik;top10_bool=hexobject$as.hexcode(inivec_bool);top10_count=inivec_count;mylik=100 top10_vbeta=inivec_vbeta; top10_veck=inivec_veck; min.top10_lik=inivec_lik;nbmodel=length(top10_count)+1;top10_fixvec=numeric(0); if (nbmodels==0) {min.top10_lik=Inf}; lastvec01=inivec_bool; modidx=length(top10_lik); #tcalls=0;inlik=0;added=0;maintained=0; list( addmodel=function(mylik,vec01,vbeta=numeric(0),fixedvec=numeric(0)) { #mylik: scalar likelihood, vec01: numeric vector of 0s and 1s, vbeta: small vector of betas (if bbeta was set to TRUE) that does NOT contain restricted zeros #use this function to add a model: #if its already among the best "nbmodels" models, its counter will be incremented by one #if it is not already in the best "nbmodels" models though its likelihood justifies that, the model will be added to the list #tcalls<<-tcalls+1 if (mylik>=min.top10_lik|length(top10_count) add model to list #added<<-added+1 nbmodel<<-min(nbmodels,length(top10_count)+1) modidx1=sum(top10_lik>mylik) top10_lik <<- append(top10_lik[-nbmodels],mylik,after=modidx1) top10_bool <<- append(top10_bool[-nbmodels],mybool,after=modidx1) top10_count <<- append(top10_count[-nbmodels],1,after=modidx1) min.top10_lik <<- top10_lik[nbmodel] if (bbeta) { #adding the model's beta vector kidx=sum(top10_veck[0:modidx1]) top10_veck <<- append(top10_veck,length(vbeta),after=modidx1)[1:nbmodel] top10_vbeta <<- append(top10_vbeta,vbeta,after=kidx)[1:sum(top10_veck[1:nbmodel])] } if (lengthfixedvec>0) { top10_fixvec <<- append(top10_fixvec,fixedvec,after=(modidx1*lengthfixedvec))[1:(nbmodel*lengthfixedvec)] } modidx<<-modidx1+1 } else { #the model is already contained in the bestof list -> just adjust counter top10_count[modidx]<<-top10_count[modidx]+1 } } } }, lik = function(){return(top10_lik)}, #return a vector of the best "nbmodels" likelihoods bool = function(){return(top10_bool)}, #return a vector of the best "nbmodels" codes as hexadecimal (e.g. model c(0,1,0,0) as "4") ncount = function(){return(top10_count)}, #return a vector of how each of the best models was chosen #counters = function(){return(c(tcalls,inlik,added,maintained))}, # for programming checks nbmodels = nbmodels, # return the maximum number of best mdoels saved in this object nregs = nmaxregressors, # return K, the maximum number of regressors overall betas_raw = function(){return(top10_vbeta)}, # return the vector of beta coefs. of the best models in one line without the zeros kvec_raw = function(){return(top10_veck)}, #return a vector that details how many coefs. each of the best models has bool_binary = function(){return(sapply(as.list(top10_bool),hexobject$as.binvec))}, #return a matrix: each column: one of the best models; rows: logical of coeficcient inclusion betas = function() { # return a matrix: rows: betas per model (including zeros); columns: model bins=(sapply(as.list(top10_bool),hexobject$as.binvec)) nbmodels_new=length(top10_lik) betamat=matrix(0,nmaxregressors,nbmodels_new) #betamat[which(bins==1)]=top10_vbeta betamat[which(bins==1)]=top10_vbeta[which(top10_vbeta!=0)] #previously =top10_vbeta return(betamat) }, fixed_vector = function(){return(matrix(top10_fixvec,lengthfixedvec))} ) } # lprob.fix returns the likelihood for the fixed vector sampling routine, hence it takes as arguments # lists for XtX.big and Xty.big, K.filt is a vector indicating the nr. of fixed regressors for each setting lprob.fix=function(positions,g0,g2,yty,k,N,null.lik,K=K,XtX.big=XtX.big,Xty.big=Xty.big,K.filt=K.filt,...){ lg02=log(g0*g2); Nlg2=(N-1)*log(g2) positions=c(1:K.filt,positions+K.filt) # if nullmodel if(sum(k)==0){ lprob=null.lik;b1new=rep(0,K);stdevnew=rep(0,K); Xty=0;XtX=0;XtXinv=0;bhat=0;ymy=yty;positions=1; } else{ k=k+K.filt XtX<-XtX.big[positions,positions] Xty<-Xty.big[positions] #do cholesky split: A=XtX=LL' #get lower triangular matrix from cholesky split L=chol(XtX) XtXinv<-chol2inv(L) bhat<-crossprod(XtXinv,Xty) ymy=yty-as.vector(crossprod(Xty,bhat)) lprob = .5*(k*lg02-(N-1)*log(ymy + g0*yty)-Nlg2)#.5*(k*log(g0*g2)-(N-1)*log((ymy + g0*yty)*g2))# b1new = g2*bhat bcov=((g0*yty+ymy)*(g2^2)/(N-2))*XtXinv #vs2=(g0*yty+ymy)*g2, stdevnew = diag(bcov)+b1new^2 } return(list(lprob=lprob,b1new=b1new,stdevnew=stdevnew, ymyOld=ymy, child.lprob = function(addix=0,dropix=0) { if (!any(as.logical(c(addix,dropix)))) {return(lprob)} if (all(as.logical(c(addix,dropix)))) { #swap jhere=(1:k)[positions==dropix]; poshere=positions[-jhere];Xj=XtXinv[,jhere];Xtxi=XtX.big[poshere,addix] bxlessj=crossprod(XtXinv,XtX.big[positions,addix])-Xj*XtX.big[addix,dropix]; bhatx=bxlessj[-jhere]-Xj[-jhere]*bxlessj[jhere]/Xj[jhere] child.ymy = as.vector(ymy+bhat[jhere]^2/Xj[jhere]-(Xty.big[addix]-crossprod(Xty.big[poshere],bhatx))^2/(XtX.big[addix,addix]-crossprod(bhatx,Xtxi))) return(.5*(k*lg02-(N-1)*log(child.ymy + g0*yty)-Nlg2)) } else { if (addix==0) { #drop jhere=(1:k)[positions==dropix] child.ymy=ymy+bhat[jhere]^2/XtXinv[jhere,jhere] return(.5*((k-1)*lg02-(N-1)*log(child.ymy + g0*yty)-Nlg2)) } else { #add Xtxi=XtX.big[positions,addix] bhatx=as.vector(crossprod(XtXinv,Xtxi)) child.ymy = as.vector(ymy - (Xty.big[addix]-crossprod(bhatx,Xty))^2 /(XtX.big[addix,addix]-crossprod(bhatx,Xtxi))) return(.5*((k+1)*lg02-(N-1)*log(child.ymy + g0*yty)-Nlg2)) } } })) } # this function returns a vector of null model likelihood values for wNr models, # each null model is the model including only the fixed terms. nullLik=function(g0,g2,yty,N,K=K,XtX.big=XtX.big,Xty.big=Xty.big,K.filt=K.filt){ lg02=log(g0*g2); Nlg2=(N-1)*log(g2) wNr=length(K.filt) nullVec=c() for(i in 1:wNr){ positions=c(1:K.filt[i]) k=K.filt[i] XtX<-XtX.big[[i]][positions,positions] Xty<-Xty.big[[i]][positions] #do cholesky split: A=XtX=LL' #get lower triangular matrix from cholesky split L=chol(XtX) XtXinv<-chol2inv(L) bhat<-crossprod(XtXinv,Xty) ymy=yty-as.vector(crossprod(Xty,bhat)) lprob = .5*(k*lg02-(N-1)*log(ymy + g0*yty)-Nlg2)#.5*(k*log(g0*g2)-(N-1)*log((ymy + g0*yty)*g2))# b1new = g2*bhat bcov=((g0*yty+ymy)*(g2^2)/(N-2))*XtXinv #vs2=(g0*yty+ymy)*g2, stdevnew = diag(bcov)+b1new^2 nullVec[i]=lprob } return(nullVec) } lik.fix=function(positions,g0,g2,yty,k,N,null.lik,XtX.big=XtX.big,Xty.big=Xty.big,K.filt=K.filt,wSampleOld=wSampleOld, K.filtOld=K.filtOld,ymyOld=ymyOld,wSampleNew=wSampleNew,...){ lg02=log(g0*g2); Nlg2=(N-1)*log(g2) positions=c(1:K.filt,positions+K.filt) # if nullmodel if(sum(k)==0){ lprob=null.lik[wSampleNew];b1new=rep(0,K);stdevnew=rep(0,K); Xty=0;XtX=0;XtXinv=0;bhat=0;ymy=yty;positions=1; lprobOld=null.lik[wSampleOld] } else{ k=k+K.filt #nr. of regressors XtX<-XtX.big[positions,positions] Xty<-Xty.big[positions] #do cholesky split: A=XtX=LL' #get lower triangular matrix from cholesky split L=chol(XtX) XtXinv<-chol2inv(L) bhat<-crossprod(XtXinv,Xty) ymy=yty-as.vector(crossprod(Xty,bhat)) lprob = .5*(k*lg02-(N-1)*log(ymy + g0*yty)-Nlg2)#.5*(k*log(g0*g2)-(N-1)*log((ymy + g0*yty)*g2))# lprobOld=.5*(K.filtOld*lg02-(N-1)*log(ymyOld + g0*yty)-Nlg2) b1new = g2*bhat bcov=((g0*yty+ymy)*(g2^2)/(N-2))*XtXinv #vs2=(g0*yty+ymy)*g2, stdevnew = diag(bcov)+b1new^2 } return(list(lprob=lprob,b1new=b1new,stdevnew=stdevnew,lprobold=lprobOld)) } #END SUBFUNCTIONS ################ #SOME PLOTS################ # plots the posterior model size modelsize.plot<<-function(X=fls.model.fix,line.width=1.5,Subset=NULL){ #library(RColorBrewer) #Colours=brewer.pal(8,"Dark2")[c(2,7)] Colours=c("#D95F02", "#A6761D") # get required information K=X$K theta=X$theta m=X$prior.msize if(X$exact){ modelSmean=sum(apply(X$topmodels,2,function(x) length(which(x==1)))*X$pmp.10[,1]) modelS.var=sum(apply(X$topmodels,2,function(x) length(which(x==1)))^2*X$pmp.10[,1])-modelSmean^2 x = apply(X$topmodels,2,function(x) length(which(x==1))) y = X$pmp.10[,1] result = c() for( i in sort(unique(x)) ) result = c(result, sum(y[which(x==i)])) names(result) = sort(unique(x)) kvec=rep(0,(K+1)) kvec[(as.numeric(names(result))+1)]=result } else{ k.vec=X$k.vec #calculate expected value summi=sum(k.vec) #modelSmean -1 because we have to deduct the constant from the mean model size modelSmean=sum((1:length(k.vec))*(k.vec/summi))-1 kvec=k.vec/sum(k.vec) modelSmean.sq=sum(((1:length(k.vec))^2)*(k.vec/summi)) modelS.var=modelSmean.sq-modelSmean^2 #var(x)=E(X^2)-(E(X))^2 } upper=min(ceiling(modelSmean+5*modelS.var), K)+1 #because of constant lower=max(floor(modelSmean-5*modelS.var),0) if(theta=="random"){ beta.bin=function(a=1,b=(K-m)/m,K=K,w=0:K){ part1=lgamma(a+b)-(lgamma(a)+lgamma(b)+lgamma(a+b+K)) part2=log(choose(K,w))+lgamma(a+w)+lgamma(b+K-w) bb=part1+part2 return(bb) } prior=beta.bin(a=1,b=(K-m)/m,K=K,w=0:K) prior=exp(prior) } if(theta!="random"){ prior=dbinom(x=0:K,size=K,prob=m/K,log=F) } mat=cbind(kvec,prior) #do the plot upper.ylim=max(kvec,prior) par(mar=rep(2,4)) #par(mar=c(1.7,2,1.7,1.7)) layout(matrix(1:2,ncol=1,nrow=2),height=c(0.85,0.15)) if(is.null(Subset)){ Subset=lower:upper } matplot(mat[Subset,],type="l",xaxt="n",col=c("steelblue3","tomato"), main=paste("Posterior Model Size Distribution","\n","Mean:", round(modelSmean,4)),cex.main=0.8,ylab="", ylim=c(0,1.1*upper.ylim),lwd=line.width) grid() points(kvec[lower:upper],cex=0.8,pch=4) if(lower==0){ axis(1, las=1, at =1:length(lower:upper), label = c(0:K)[lower:(upper+1)],cex.axis=0.7) } else{ axis(1, las=1, at =1:length(lower:upper), label = c(0:K)[lower:upper],cex.axis=0.7) } #box(lwd=2) plot(1,type="n",xlab="",ylab="",axes=FALSE) par(mar=rep(0,4)) legend(x="center",lty=c(1,2),legend=c("Posterior","Prior"),col=c("steelblue3","tomato"),ncol=2,bty="n",cex=1,lwd=line.width) # plot(k.vec[range],type="l",main=paste("Posterior Model Size\n","Mean(",round(modelSmean,digits=4),")",sep=""), # cex.main=0.8,ylab="",xlab="Nr. of Regressors",xaxt="n",col="white",yaxs="i", # ylim=c(0,max(k.vec)*1.10)) # polygon(c(-0.1,k.vec[range],-0.1),col="lavender",border="gray",lwd=2) layout(matrix(1)) return(list(mean=modelSmean,var=modelS.var)) } #modelsize.plot(model1,line.width=2,Subset=c(0:25)) # plots the density if you have saved betas in your bma routine # reg specifies the regressor name you want to plot, if null # all regressors are plotted subsequentially, # Range specifies the range of the density plot if you want to # cusomize your plot for publication output for example plotDensity<<-function(X,lwd=2,ask.set=T,reg=NULL,Range=NULL ,PIP=NULL,Title=""){ if(length(grep("bma",class(X)))>0){ # if class X is bma (e.g. bma sar) betas=X$beta.draws if(sum(betas)==0){ stop("Run estimation again and save betas via setting beta.save=TRUE") } if(class(X)=="bma.sar"){ X$estimates= rbind(X$estimates,1) rownames(X$estimates)[nrow(X$estimates)]="rho" #include rho in the estimation table in order to get a "pip" for it } reg.names=X$reg.names rownames(betas)[1:length(reg.names)]=reg.names beta.nr=nrow(as.matrix(betas)) PIP=X$estimates[reg,1] Betas=betas[reg,,drop=F] } if(is.vector(X)){ reg=NA # have to set it to something that is not null cause if X is a submitted matrix, then cannot plot it for more regs in a loop Betas=X } densPlot=function(betas=betas,pip=pip,Range=NULL,Title=""){ if(!is.null(pip)){ par(mar=c(0,2,0,2)) layout(matrix(1:2,ncol=1,nrow=2),height=c(0.1,0.87)) plot(1,type="n",xlim=c(0,1),ylim=c(0,0.5),xaxt="n",ylab="",xlab="",bty="n",yaxt="n") rect(0,0,1,0.25,col="lavender") rect(0,0,pip,0.25,col="darkred") } par(mar=c(3,2,2,2)) #default is c(5, 4, 4, 2) c(bottom, left, top, right) dens=density(betas) x=dens$x;y=dens$y if(Title=="" & !is.null(rownames(Betas))){ Title=paste("Unconditional Posterior Distribution of", rownames(betas), "(", ncol(betas), "Obs.)") } plot(dens,lwd=lwd,col="steelblue4", xlab="",xlim=Range, main=Title) grid() up=which(x>mean(betas))[1] lines(rep(mean(betas),2),c(0,mean(y[(up-1):up])),col="darkred",lwd=2) up=which(x>median(betas))[1] lines(rep(median(betas),2),c(0,mean(y[(up-1):up])),col="darkseagreen",lwd=2,lty=2) legend(x="topleft",lty=c(1,2),lwd=lwd,col=c("darkred","darkseagreen"),legend=c("Mean", "Median")) # set back the options for standard plotting par(mar=c(5.1,4.1,4.1,2.1)) layout(1) } if(is.null(reg)){ par(ask=ask.set) for(i in 1:beta.nr){ densPlot(betas=betas[i,,drop=F],pip=X$estimates[reg.names[i],1],Range=Range,Title=paste("Unconditional Posterior Distribution of", rownames(betas)[i], "(", ncol(betas), "Obs.)")) } } else{ par(ask=F) densPlot(betas=Betas,pip=PIP,Range=Range,Title=Title) } } #plotDensity(modelEBL$topmodobject$fixed_vector()[1,]) #plotDensity(modelEBL,reg="LabForce",Range=c(0.001,0.002)) #compare MCMC and exact PMP's plotConv<<-function(X,lwd=2){ if(class(X)!="bma" & class(X)!="bma.sar" & class(X)!="bma.fcast"){ stop("submit an object of class bma") } # now get PMP analytical and MCMC mat=t(X$topmodels[c("PMP (Exact)","PMP (MCMC)"),]) op=par(mar=c(0,2,4,2),bg="white",cex=0.7) layout(matrix(1:2,ncol=1,nrow=2),height=c(0.85,0.15)) matplot(mat,type="l",lty=1,col=c("tomato","steelblue3"),lwd=lwd, main=paste("Posterior Model Probabilities\n(Corr:",X$info["Corr PMP"],")",sep=""), xlab="Index of Models",ylab="") grid() par(mar=c(0,2,0,2)) plot(1,type="n",xlab="",ylab="",axes=FALSE,xlim=c(-1,1),ylim=c(-1,1)) legend("center",lty=1,legend=c("PMP (Exact)", "PMP (MCMC)"),col=c("tomato","steelblue3"),ncol=2,bty="n",cex=1,lwd=2) layout(matrix(1)) par(mar=c(5, 4, 4, 2) + 0.1) } # End plots ################################################################################################# # # This is the end of the functions we have sourced (auxiliary and plots) # ############################################################################################################## # Here the bma program starts with initializing and calculating some necessary quantities # exact=F;int=F y<-as.matrix(X.data[,1]) X<-as.matrix(X.data[,2:ncol(X.data)]) N<-nrow(X) K=ncol(X) # User Checks: if(prior.msize>=K){ cat("Submitted prior model size is >= than the nr. of regressors\n, used K/2 instead\n\n") prior.msize=K/2 } if(exact & !beta.save){ cat("For exact posterior inference betas have to be saved\n, have set beta.save=TRUE instead\n\n") beta.save=TRUE } if (nmodel[1]<0|is.na(nmodel[1])) {nmodel=0} # subtract mean from all regressors as in FLS y.mean=mean(y) y<-y-matrix(y.mean,N,1,byrow=T) X.mean=colMeans(X) X<-X-matrix(X.mean,N,K,byrow=T) # multiply the whole matrix stuff out before going into the simulation loops XtX.big=crossprod(X) Xty.big=crossprod(X,y) yty = as.vector(crossprod(y)) ###################################################################################################################################### #The function Starter selects randomly a start matrix and runs a #regression. From this regression, the #start Design matrix is that for which the t-stats are >0.2. So we #can circumvent starting from a very bad start point. start.list=starter(K,start.value,y,N=N,XtX.big=XtX.big,Xty.big=Xty.big,X=X) molddraw=start.list$molddraw start.position=start.list$start.position kold=sum(molddraw) position=(1:K)[molddraw==1] ###################################################################################################################################### # for the case that X contains interaction terms if(int){ if(length(grep("[.x.]",colnames(X.data),extended=T))==0){ stop("Please separate column names of interaction terms by .x. (e.g. A.x.B)") } # this file constructs the molddraw matrix source(paste(sPath,"interaction_Sampler.r",sep=""),local=T) # overwrite matrix multiplication and start model X=cbind(X.base,X.int) XtX.big=crossprod(X) Xty.big=crossprod(X,y) # we start with a model drawn from the pool of base models (without the interaction terms) start.list=starter(K=k.base,start.value,y,N=N,XtX.big=XtX.big,Xty.big=Xty.big,X=X) molddraw=c(start.list$molddraw,rep(0,k.int)) start.position=start.list$start.position kold=sum(molddraw) position=which(molddraw==1) } else{ mMinus<-mPlus<-NA } ######################################################################################################################################## #specify g0 for the g-prior # use benchmark ric criterion as standard setting if(g=="bric" | !is.numeric(g)){ if (N<=(K^2)){ g0=1/(K^2) } else{ g0=1/N } } # if user specifies a number for g, use this number instead of bric if(is.numeric(g)){ g0=g } g2=1/(g0+1) ################################################################################################ #Initializing ############################################################################################### # calculate Likelihood for NullModel null.lik=((1-N)/2)*log(yty) null.count=0 models.visited=0 # initialize top 10 function lik.list=lprob(positions=(1:K)[molddraw==1],g0=g0,g2=g2,yty=yty,k=kold,N=N,null.lik=null.lik, K=K,XtX.big=XtX.big,Xty=Xty.big) lprobold=lik.list$lprob b1=lik.list$b1new stdev=lik.list$stdevnew # specify model.odds function 1) "fixed theta", 2) or "uniform prior" or 3) for all other specifications use the random theta specification pmp = switch(theta,"fix"=pmp.fix,"uniform"={prior.msize<<-K/2; pmp.unif},pmp.ran) # calculate the posterior model probability for the first model pmp2=pmp(ki=kold,K=K,m=prior.msize, lprobnew=lprobold) # initialize top 10 function topmods=top10(nmaxregressors=K,nbmodel=nmodel,bbeta=beta.save, inivec_lik=pmp(ki=0,K=K,m=prior.msize,lprobnew=null.lik)) topmods$addmodel(mylik=pmp2,vec01=molddraw,vbeta=b1) # Initialize the rest mstart=molddraw lprobstart=pmp2 inccount=rep(0,K) msize=0 k.vec=rep(0,K) b1mo=numeric(K) #calculate first and second moment of all coefficients ab=numeric(K) #Initialize them here b2mo=numeric(K) bb=numeric(K) mnewdraw=numeric(K) int.terms.changed=FALSE ####################################################################### ############################################################################################# nrep=burn+iter #Burnins + Draws set.seed(as.numeric(Sys.time())) #Set Seed randomly for number generator t1<-Sys.time() #Save time before going into the loop # generate logfile if desired if(logfile){ if (is.character(logfile)) { sfilename=logfile} else { sfilename="fls.log" } file.create(sfilename) cat(as.character(Sys.time()),": starting loop ... \n",append=TRUE, file=sfilename) #write one line fact=max(floor(nrep/100),logstep) } if(mcmc!="bd" & !int){ sampling=rev.jump } if(mcmc!="bd" & int){ sampling=rev.jump.int } if(mcmc=="bd" & int){ sampling=fls.samp.int } if(mcmc=="bd" & !int){ sampling=fls.samp } ########################################################################################### #START MAIN LOOP ########################################################################################### for (i in 1:nrep){ if(logfile){ if (i %% fact==0) { cat(as.character(Sys.time()),":",i,"current draw \n", "kold:",kold, "\n",append=TRUE, file=sfilename) #write one line } } ########################################################################################## #Start sampling program ########################################################################################### a=sampling(molddraw=molddraw,K=K,mMinus=mMinus,mPlus=mPlus) mnewdraw=a$mnewdraw positionnew=a$positionnew knew=length(positionnew) if (int) {if (length(a$addi)>1|length(a$dropi)>1) {int.terms.changed=TRUE} else {int.terms.changed=FALSE}} #int.terms.changed = TRUE if there were multiple regs dropped or added due to interaction terms if (int.terms.changed) { lik.list.int=lprob(positions=positionnew,g0=g0,g2=g2,yty=yty,k=knew,N=N,null.lik=null.lik, K=K,XtX.big=XtX.big,Xty=Xty.big) #in case of changing interaction terms, draw the big likelihood lprobnew = lik.list.int$lprob } else { lprobnew=lik.list$child.lprob(a$addi,a$dropi) #if standard sampling, use Frisch-Waugh to get the new lprob (faster) } pmp1=pmp(ki=knew,K=K,m=prior.msize, lprobnew=lprobnew) lratio=pmp1-pmp2 #if (int.terms.changed) browser() #Now decide whether to accept candidate draw if(log(runif(1,min=0,max=1))< lratio){ if (!int.terms.changed) { # in case one has used Frisch-Waugh and the new model got accepted, # calculate the 'real' inverse in order not to make copying mistakes lik.list=lprob(positions=positionnew,g0=g0,g2=g2,yty=yty,k=knew,N=N,null.lik=null.lik, K=K,XtX.big=XtX.big,Xty=Xty.big) } else { lik.list=lik.list.int } position = positionnew lprobold=lik.list$lprob pmp2=pmp1 # get posterior odds for new model if accepted molddraw=mnewdraw kold=knew models.visited=models.visited+1 #does not account for revisiting models } # Collect Posterior Draws ######################################################################## if (i>burn){ stdev=lik.list$stdevnew; b1=lik.list$b1new inccount = inccount + molddraw # compute average size of models msize=msize + kold # collect nr. of regressors in vector k.vec[kold]=k.vec[kold]+1 if(!any(kold)){ null.count=null.count+1 } # add log(lik)*p(M) to topmodels topmods$addmodel(mylik=pmp2,vec01=molddraw,vbeta=b1) #calculating posterior properties of coefficients means b1mo[position]=b1mo[position]+b1; b2mo[position]=b2mo[position]+stdev; } } ########################################################################################### #END MAIN LOOP ########################################################################################### ########################################################################################### t2<-Sys.time() timed<-difftime(t2,t1) post.inf=post.calc(inccount=inccount,iter=iter,msize=msize,b1mo=b1mo,b2mo=b2mo,y.mean=y.mean,X.mean=X.mean, k.vec=k.vec,null.count=null.count,models.visited=models.visited,exact=exact,beta.save=beta.save,SAR=F, topmods=topmods,X.data=X.data,burn=burn,timed=timed,theta=theta,m=prior.msize) ########################################################################################### # print results to console if(printRes){ print(post.inf$post.mean) cat("\n") print(post.inf$info) cat("\n") print(timed) } result=list(estimates=post.inf$post.mean,time=timed,topmodels=post.inf$topmod,info=post.inf$info, k.vec=post.inf$k.vec,start.pos=sort(start.position),beta.draws=post.inf$beta.draws, topmodprob=post.inf$topmodprob,pmp.10=post.inf$pmp.10,topmodobject=topmods,theta=theta, K=K, prior.msize=prior.msize,cons=post.inf$cons,reg.names=post.inf$reg.names,exact=exact,fls.call=sys.call(0)) class(result)=c("bma") # do modelsize plot op=par(ask=ask.set) try(modelsize.plot(X=result,line.width=2),silent=T) try(plotConv(X=result,lwd=2)) par(op) return(result) } ########################################################################################################################################## # Example: # read in the data (see link "fls_data") #put the fls_data.txt file in your working directory dataM=as.matrix(read.table("fls_data.txt",header=T)) # we don't need the first four columns K=ncol(dataM)-1 # nr. of regressors # this setting corresponds to a uniform prior on the model space (prior.msize=K/2 and theta="fix") # and the ric specification since K^2> N (with N the nr. of observations) as suggested by fls model.ric=fls(X.data=dataM,burn=1e05,iter=1e06,g=(1/K^2),nmodel=100,theta="fix", prior.msize=K/2,logfile=T,mcmc="bd",start.value=rep(0,K),beta.save=T) # make a posterior density plot for a specific variable plotDensity(model.ric,reg="LifeExp") # look at the objects saved in the bma object names(model.ric) # access the objects of interest model.ric$topmodels # this is the matrix of the "nmodel" best models with (normalized) posterior model probabilities given at the very bottom # End Example #########################################################################################################################################