Functionality/option to combine estimates

[prabasaj]

Suggesting an option to combine multiple estimates of the same parameter set (shape and scale of gamma distribution fitted to incubation period for Ebola, for example, across 5 studies). Here's some basic code to start with:

Weighted combination of multivariate estimates and population covariance
User supplied weights may be some combination: (weight) = (effective sample size) X (quality weight) X (relevance weight)
Neglects sample size corrections and is based on definition
Includes option of weighting by inverse radius of covariance matrix
estimate set is a list of estimate-lists: list(est1 = list(mean, covariance, weight), ...)

weighted_estimates = function(estimate_set, inverse_variance_weighting=FALSE) {

if (inverse_variance_weighting) {
wts = proportions(sapply(estimate_set, function(ll) 1/max(eigen(ll$covariance)$values) ))
} else {
wts = proportions(sapply(estimate_set, function(ll) ll$weight))
}

mu = c(sapply(estimate_set, function(ll) ll$mean) %% wts)
sigma = matrix(sapply(estimate_set, function(ll) (ll$mean - mu) %% t(ll$mean - mu)) %*% wts, nrow=length(mu)) +
matrix(rowSums(mapply(*, wts, lapply(estimate_set, function(ll) ll$covariance))), nrow=length(mu))
list(mean = mu, covariance = sigma)
}

Example: combining bivariate parameter estimates (shape and scale of gamma distributions mean=c(shape, scale), for example)

estimate_set = list(
list(mean=c(1, 2), covariance=cbind(c(0.25, 0.1), c(0.1, 0.36)), weight=2),
list(mean=c(1.2, 2.5), covariance=cbind(c(0.25, 0.1), c(0.1, 0.36)), weight=3)
)
weighted_estimates(estimate_set, inverse_variance_weighting=TRUE)
weighted_estimates(estimate_set)