



data(dataBoston)
data(WL.boston)
#now do an MC3 sampling over the boston housing data 'dataBoston' with 1000 burn-ins,
#9000 iterations (ex burn-ins),
#and retaining the best 100 models (besides overall MC3 frequencies)
invisible(readline("hit <Return> to estimate a short BMA MC3 sampling chain."))
 mboston =spatFilt.bms(X.data=as.matrix(dataBoston),WList=WL.boston,burn=1000,iter=9000,nmodel=100,user.int=T)
#The console printout shows some information on the resulting bma object 'mboston'.
#The chart details the prior and posterior model size (above), and the likelihoodos and MCMC frequencies of the best 100 models (below).
#Setting user.int=FALSE suppresses his kind of output.\n\n ")

# some results:
 info.bma(mboston)
# info.bma() shows basic aggregate results from MC3 sampling

#### a function on the coefficients
 estimates.bma(mboston)
 #based on MCMC frequencies, this shows
 # the posterior inclusion probabilities (column 1),
 # the unconditional expecteded value of coefficients (column 2),
 # their standard deviations (column 3),
 # the percentage of times the coefficents had a positive sign
 # (conditional on inclusion, column 4)

#### a different: version
 estimates.bma(mboston,exact=T,std.coefs=T)
 #this is similar to estimates.bma(),
 # however here the numbers are based on the exact marginal
 # likelihoods of the best (100) models drawn.
 # Moreover the coefficents are shown in standardized form.

invisible(readline("Hit <Return> for other low-level commands."))
 pmp.bma(mboston)
  #post. model probs. for top 100 models based on MCMC freqs and likelihoods

 pmpW.bma(mboston)
 #post. inclusion probs. of weight matrices for top 100 models based on MCMC freqs and likelihoods
 
beta.draws.bma(mboston[1:3])
 # show the estimates for the best 3 models
 # the column names are the inclusion vectors in hex-code
 # (e.g. 101 for inclusion of variables 1 and 3 is "5" in hexciode )

mboston[3]$topmod
 #shows the third-best model


 
 mboston$wTopModels[,1:3]
 #shows disaggregated results in matrix form of topmodels along with the according post. model probabilities; W-Index refers to the weight matrices stored in WList
 # here results are shown for best 3 models
 
 topmodels.bma(mboston)[,1:3]
 # shows a matrix of topmodels along with the according post. model probabilities (exact and frequencies) after integrating out uncertainty w.r.t. the weight matrices
 # ie if the first two models would be the same in terms of expl. variables but differ wrt to the weight matrix used, the post. mass of the two models is aggregated
 # here results are shown for best 3 models
 
 
 
   library(spdep)
   data(boston.soi)
   mTest=moranTest.bma(mboston,variants="double",W=nb2listw(boston.soi),nmodel=3)
  # test the best 3 models in terms of post. model probabilities for remaining spatial autocorrelation in the residuals based on weight matrix boston soi


  pvalueMat=cbind(sapply(mTest$moran,function(x) x$p.value),sapply(mTest$moranEV,function(x) x$p.value))
  #extract p-values
 
 
 
invisible(readline("Hit <Return> for some plots "))
oldask=par()$ask;par(ask=TRUE)
  density(mboston,reg="CRIM")
#plot density for regressor "BlMktPm"

 image(mboston[1:20],FALSE,cex=0.6)
 #plot signs (pos=blue, neg=red, not included=white) for best 20 models


 plotModelsize(mboston,exact=TRUE)
 #plot prior and posterioro model size based on exact likelhoods
 # of best (100) models



 plot(mboston)
 #a combined plot

 #compare boxplots
  Colours=c("#1B9E77", "#D95F02", "#7570B3", "#E7298A", "#66A61E" ,"#E6AB02", "#A6761D", "#666666")
  par(mfrow=c(1,2))
  boxplot(pvalueMat[,1],main="Box plot of p-values \n OLS-regression",xaxt="n",col=Colours[3],outline=F,ylim=c(0,1),cex.main=0.8)
  grid()
  boxplot(pvalueMat[,2],main="Box plot of p-values \n Eigenvector augmented regression",xaxt="n",col=Colours[3],outline=F,ylim=c(0,1),cex.main=0.8)
  grid()
  par(mfrow=c(1,1))



par(ask=oldask); rm(oldask)



