Nothing
#require(multicore)
#
#Compute marginal for rho using marginal likelihhod x prior
#Compute the marginals for all the other parameters using BMA
#
#Version of splinefun to return 0 outside x-range
#
mysplinefun<-function(x, y = NULL,
method = c("fmm", "periodic", "natural", "monoH.FC")[1],
ties = mean)
{
xmin<-min(x)
xmax<-max(x)
ff<-splinefun(x=x, y=y, method=method, ties=ties)
fff<-function(x, deriv)
{
sapply(x, function(X){
if(X<xmin |X>xmax)
return(0)
else
return(ff(X))
})
}
}
#Compute and re-scale marginal using splines
#x: x-valuies
#logy: log(y)-values (posterior LOG-density)
#logp: log-prior for x-values
#usenormal: Use Normal approximation??
#
fitmarg<-function(x, logy, logp=0, usenormal=FALSE)
{
if(!usenormal)
{
#func = splinefun(x, logy-max(logy))
#post = exp(logp+func(x))
logpost<-logy-max(logy)+logp-max(logp)
post<-exp(logpost)
post.func = splinefun(x, post)
## normalize
xx = seq(min(x), max(x), len = 1000)
#z = sum(post.func(xx)) * diff(xx)[1]
z=integrate(post.func, min(x), max(x))$value
post.func = mysplinefun(x=xx, y=(post.func(xx) / z) )
}
else
{
xx = seq(min(x), max(x), len = 1000)
meany=sum(x*exp(logp+logy))/sum(exp(logp+logy))
sdy=sqrt( sum(((x-meany)^2)*exp(logp+logy))/sum(exp(logp+logy)) )
post.func = function(x){dnorm(x, mean=meany, sd=sdy)}
}
return(post.func)
}
#
#Fits a marginal using BMA
#margs: list of matrices with the marginals are in INLA (i.e., 2-col matrix)
#ws: weights
fitmargBMA<-function(margs, ws, len=100)
{
ws<-ws/sum(ws)
xmin <- quantile((unlist(lapply(margs, function(X){min(X[,1])}))), 0.25)
xmax <- quantile(unlist(lapply(margs, function(X){max(X[,1])})), 0.75)
xx<-seq(xmin, xmax, len=len)
margsws<-lapply(1:length(margs), function(i){
func<-fitmarg(margs[[i]][,1], log(margs[[i]][,2]))
ws[i]*func(xx)
})
margsws<-do.call(cbind, margsws)
d<-data.frame(x=xx, y=apply(margsws, 1, sum))
names(d)<-c("x", "y")
return(d)
}
#
#Weighted sum of summary results in matrices
#
fitmatrixBMA<-function(models, ws, item)
{
lmatrix<-lapply(models, function(X){X[[item]]})
auxbma<-ws[1]*lmatrix[[1]]
for(i in 2:length(lmatrix)){auxbma<-auxbma+ws[i]*lmatrix[[i]]}
return(auxbma)
}
#
#Weighted sum of summary results in lists
#
fitlistBMA<-function(models, ws, item)
{
nlist<-names(models[[1]][[item]])
auxlist<-as.list(rep(NA, length(nlist)))
names(auxlist)<-nlist
for(ele in nlist)
{
lmatrix<-lapply(models, function(X){X[[item]][[ele]]})
auxbma<-ws[1]*lmatrix[[1]]
for(i in 2:length(lmatrix)){auxbma<-auxbma+ws[i]*lmatrix[[i]]}
auxlist[[ele]]<-auxbma
}
return(auxlist)
}
#
#Weighted sum of marginals
#
#The elements to handle are lists with different marginals (in matrix form)
#
fitmargBMA2<-function(models, ws, item)
{
if(is.null(models[[1]][[item]]))
return(NULL)
nlist<-names(models[[1]][[item]])
auxlist<-as.list(rep(NA, length(nlist)))
names(auxlist)<-nlist
for(ele in nlist)
{
lmatrix<-lapply(models, function(X){X[[item]][[ele]]})
#xxr<-range(unlist(lapply(lmatrix, function(X){X[,1]})))
#Compute range using a weighted sum of min and max
xxr<-c(NA, NA)
xxr[1]<-sum(ws*unlist(lapply(lmatrix,function(X){min(X[,1])})))
xxr[2]<-sum(ws*unlist(lapply(lmatrix,function(X){max(X[,1])})))
xx<-seq(xxr[1], xxr[2], length.out=81)
auxbma<-rep(0, length(xx))
for(i in 1:length(lmatrix)){
auxspl<- mysplinefun(x=lmatrix[[i]][,1], y=lmatrix[[i]][,2])
auxbma<-auxbma+ws[i]*auxspl(xx)
}
auxlist[[ele]]<-cbind(xx, auxbma)
}
return(auxlist)
}
#Returns a summary of the fitted.values using BMA
#
#models: list of INLA models
#rrho: Vector of values of rho (same as number of models)
#logrhoprior:vector of logdensities of prior for rho
BMArho<-function(models, rho, logrhoprior=rep(1, length(rho)) )
{
mlik<-unlist(lapply(models, function(X){X$mlik[1]}))
post.func <- fitmarg(rho, mlik, logrhoprior)
# mlik.func = splinefun(rho, mlik - max(mlik))
#
# post = exp(logrhoprior) * exp(mlik.func(rho))
# post.func = splinefun(rho, post)
#
# ## normalize
# rrho = seq(min(rho), max(rho), len = 1000)
# z = sum(post.func(rrho)) * diff(rrho)[1]
# post.func = splinefun(rrho, post.func(rrho) / z)
#Weights for BMA
ws<-(post.func(rho))
ws<-ws/sum(ws)
#Fitted values with BMA
#Compute BMA
fvalues<-mclapply(1:length(models), function(X){ws[X]*models[[X]]$summary.fitted.values[,1]})
fvalues<-data.frame(fvalues)
fvalues<-apply(fvalues, 1, sum)
return(fvalues)
}
#Elements worth BMA
#
# "summary.fixed" "marginals.fixed"
# "summary.lincomb" "marginals.lincomb"
# "summary.random" "marginals.random"
# "summary.linear.predictor" "marginals.linear.predictor"
# "summary.fitted.values" "marginals.fitted.values"
#Performs BMA on a number of things
#
#
#
INLABMA<-function(models, rho, logrhoprior=rep(1, length(rho)), impacts=FALSE, usenormal=FALSE )
{
#require(INLA)
mlik<-unlist(lapply(models, function(X){X$mlik[1]}))
post.func <- fitmarg(rho, mlik, logrhoprior, usenormal)
#Weights for BMA
ws<-(post.func(rho))
ws<-ws/sum(ws)
mfit<-list(rho=list())
mfit$rho$marginal<-data.frame(x=seq(min(rho), max(rho), len=100))
mfit$rho$marginal$y<-post.func(mfit$rho$marginal$x)
mfit$rho$marginal<-as.matrix(mfit$rho$marginal)
# mfit$rho$mean<-sum(mfit$rho$marginal$x*mfit$rho$marginal$y)*diff(mfit$rho$marginal$x)[1]
# mfit$rho$sd<-sqrt(sum(((mfit$rho$marginal$x-mfit$rho$mean))^2*mfit$rho$marginal$y)*diff(mfit$rho$marginal$x)[1])
#Compute some quantiles
#fquant<-function(x, ff, qtile){integrate(ff, 0, x)$value-qtile}
#qvect<-c(.025, .5, .975)
#mfit$rho$quantiles<-sapply(qvect,
# function(qtile){
# uniroot(fquant, interval=c(min(rho),max(rho)), ff=post.func, qtile=qtile)$root
# }
#)
#names(mfit$rho$quantiles)<-as.character(qvect)
#Summary statistics
margsum <- INLA::inla.zmarginal(mfit$rho$marginal, TRUE)
mfit$rho$mean<-margsum$mean
mfit$rho$sd<-margsum$sd
mfit$rho$quantiles<-unlist(margsum[-c(1:2)])
#Results which are stored in matrices
mateff<-c("summary.fixed", "summary.lincomb",
#"summary.random",
"summary.linear.predictor", "summary.fitted.values")
#summary.fixed
#mfit$summary.fixed<-fitmatrixBMA(lapply(models, function(X){X$summary.fixed}), ws)
lmat<-mclapply(mateff, function(X){fitmatrixBMA(models, ws, X)})
names(lmat)<-mateff
mfit<-c(mfit, lmat)
#Results stored in lists
listeff<-c("dic", "cpo")
leff<-mclapply(listeff, function(X){fitlistBMA(models, ws, X)})
names(leff)<-listeff
mfit<-c(mfit, leff)
#Results are marginals
listmarg<-c("marginals.fixed", "marginals.lincomb",
"marginals.lincomb.derived", #"marginals.random",
"marginals.linear.predictor", #"marginals.fitted.values",
"marginals.hyperpar", "marginals.spde2.blc")
margeff<-mclapply(listmarg, function(X){fitmargBMA2(models, ws, X)})
names(margeff)<-listmarg
mfit<-c(mfit, margeff)
#Impacts
mfit$impacts<-FALSE
if(impacts)
{
mfit$impacts<-TRUE
summimp<-c("summary.total.impacts", "summary.direct.impacts",
"summary.indirect.impacts")
matsummimp<-mclapply(summimp, function(X){fitmatrixBMA(models, ws, X)})
names(matsummimp)<-summimp
mfit<-c(mfit, matsummimp)
margimp<-c("marginals.total.impacts","marginals.direct.impacts",
"marginals.indirect.impacts")
lmargimp<-mclapply(margimp, function(X){fitmargBMA2(models, ws, X)})
names(lmargimp)<-margimp
mfit<-c(mfit, lmargimp)
#Recompute impacts summaries
mfit<-recompute.impacts(mfit)
}
return(mfit)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.