Weighted Likelihood Empirical Bayes
Compute empirical Bayes moderated parameter estimators using a weighted likelihood approach.
WLEB(theta, loglik, prior.n = 5, covariate = NULL, trend.method = "locfit", span = NULL,
overall = TRUE, trend = TRUE, individual = TRUE, m0 = NULL, m0.out = FALSE)theta |
numeric vector of values of the parameter at which the log-likelihoods are calculated. |
loglik |
numeric matrix of log-likelihood of all the candidates at those values of parameter. |
prior.n |
numeric scaler, estimate of the prior weight, i.e. the smoothing parameter that indicates the weight to put on the common likelihood compared to the individual's likelihood. |
covariate |
numeric vector of values across which a parameter trend is fitted |
trend.method |
method for estimating the parameter trend. Possible values are |
span |
width of the smoothing window, as a proportion of the data set. |
overall |
logical, should a single value of the parameter which maximizes the sum of all the log-likelihoods be estimated? |
trend |
logical, should a parameter trend (against the covariate) which maximizes the local shared log-likelihoods be estimated? |
individual |
logical, should individual estimates of all the candidates after applying empirical Bayes method along the trend be estimated? |
m0 |
numeric matrix of local shared log-likelihoods. If |
m0.out |
logical, should local shared log-likelihoods be included in the output? |
This function implements a very general empirical Bayes strategy outlined by McCarthy et al (2012). McCarthy et el used the method to estimate genewise negative binomial dispersion parameters, but here the method is generalized to apply to any parameter and any likelihood function. The method gives similar results to parametric empirical Bayes with a conjugate prior, when a conjugate prior exists, but does not require the prior distribution to be specified. The prior distribution is instead inferred from the pooled likelihood, i.e., from the likelihood that would be arise from pooling all the cases or from pooling all the cases with similar covariate values.
The function assumes a series of cases. Each case leads to data set from which a parameter (theta) is to be estimated. For each case, the log-likelihood function has been evaluated over a grid of possible values for theta. The function takes as input the matrix of log-likelihood values where the rows correspond to cases and the columns correspond to putative parameter values.
Each case is associated with a covariate value that might affect theta.
The "overall" parameter estimate is the maximum likelihood estimator of theta that arises if the likelihood is averaged over cases and then maximized over theta.
The "trend" parameter estimates are estimates for theta that arise if each column of loglik is replaced by a smooth trend with respect to the covariate.
The "individual" parameter estimate for each case is a compromise between the maximum likelihood estimate for that case alone and a global parameter estimate computed from all the cases, the latter being either the overall estimate (if trend.method="none" or the trend estimate (otherwise).
A list with the following:
overall |
the parameter estimate that maximizes the sum of all the log-likelihoods. |
trend |
the estimated trended parameters against the covariate. |
individual |
the individual estimates of all the candidates after applying empirical Bayes method along the trend. |
shared.loglik |
the estimated numeric matrix of local shared log-likelihoods |
Yunshun Chen, Gordon Smyth
McCarthy, DJ, Chen, Y, Smyth, GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Research 40, 4288-4297. https://doi.org/10.1093/nar/gks042
locfitByCol, movingAverageByCol and loessByCol implement the local fit, moving average or loess smoothers.
y <- matrix(rpois(100, lambda=10), ncol=4) theta <- 7:14 loglik <- matrix(0,nrow=nrow(y),ncol=length(theta)) for(i in 1:nrow(y)) for(j in 1:length(theta)) loglik[i,j] <- sum(dpois(y[i,], theta[j] ,log=TRUE)) covariate <- log(rowSums(y)) out <- WLEB(theta, loglik, prior.n=3, covariate) out
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.