DREAM-derived BMA weights

[milonbrri]

How can calculate Marginal posterior PDFs of the DREAM-derived BMA weights of monthly rainfall ?? How can I get the script? Please help

[florianhartig]

1. fit your model with the DREAM algorithm
2. calculate weights with the marginalLikelihood functions

The usual caveats on ML calculations noted in the help of marginalLikelihood apply!

[SalamBRRI]

Dear Florianhartg Thank you for your nice comment and instruction. So far I understand, MilonBrri want to calculate ensemble model from different climate model of GCM using DREAM algorithm of bayesianTool packages for weight BMA model. To do so, he is searching R script. As example there are 15 model data for rainfall, he want to make ensemble model from 15 model. DREAM is popular algorithm. Is marginalLikelihood used in order to estimate parameter of each model such as: alpha+beta*model_1 separately. I mean 15 times for 15 model needs to run marginalLikelihood function. Can you help him a bit detail.

[florianhartig]

Hi there,

I'm sorry for the slow reply, I was overwhelmed with online teaching this spring. @milonbrri - have you checked the help of ?marginalLikelihood? Doesn't this answer your question?

Best, F

[florianhartig]

For completeness, this is the part of the help I'm referring to:

##############################################################
# Comparison of ML for two regression models

# Creating test data with quadratic relationship
sampleSize = 30
x <- (-(sampleSize-1)/2):((sampleSize-1)/2)
y <-  1 * x + 1*x^2 + rnorm(n=sampleSize,mean=0,sd=10)
# plot(x,y, main="Test Data")

# likelihoods for linear and quadratic model 
likelihood1 <- function(param){
  pred = param[1] + param[2]*x + param[3] * x^2
  singlelikelihoods = dnorm(y, mean = pred, sd = 1/(param[4]^2), log = TRUE)
  return(sum(singlelikelihoods))  
}
likelihood2 <- function(param){
  pred = param[1] + param[2]*x 
  singlelikelihoods = dnorm(y, mean = pred, sd = 1/(param[3]^2), log = TRUE)
  return(sum(singlelikelihoods))  
}

setUp1 <- createBayesianSetup(likelihood1, 
                              lower = c(-5,-5,-5,0.01), 
                              upper = c(5,5,5,30))
setUp2 <- createBayesianSetup(likelihood2, 
                              lower = c(-5,-5,0.01), 
                              upper = c(5,5,30))

out1 <- runMCMC(bayesianSetup = setUp1)
M1 = marginalLikelihood(out1, start = 1000)

out2 <- runMCMC(bayesianSetup = setUp2)
M2 = marginalLikelihood(out2, start = 1000)

### Calculating Bayes factor

exp(M1$ln.ML - M2$ln.ML)

# BF > 1 means the evidence is in favor of M1. See Kass, R. E. & Raftery, A. E. 
# (1995) Bayes Factors. J. Am. Stat. Assoc., Amer Statist Assn, 90, 773-795.

### Calculating Posterior weights

exp(M1$ln.ML) / ( exp(M1$ln.ML) + exp(M2$ln.ML))

# If models have different model priors, multiply with the prior probabilities of each model. 

As I said, the usual caveats of the ML apply, if you go on in the help it shows how to calculate a fractional BF, which I would recommend for calculating posterior weights.