Meta-analysis#
meta_analysis() combines, for every gene, the effects of
several studies and their standard errors into an estimate of a common effect, as
metafor’s rma.uni does without moderators: the fixed-effect (inverse-variance)
estimate, or a random-effects estimate with the between-study variance by
DerSimonian-Laird, REML or Paule-Mandel, with Cochran’s Q, I-squared, H-squared and
z-tests or the Knapp-Hartung adjustment. A gene measured in some studies only is
analysed over those.
differential_expression_meta_analysis() does it for the
results of limma_voom() or
limma_trend() on several cubes, such as the
same cell types in several cohorts: it matches groups, contrasts and genes by their ids,
reads each table’s log2 fold changes and moderated standard errors, and returns a table
per cell type and contrast with p-values adjusted over its genes, the studies used and
left out with the reason, and the provenance of every study.
Meta-analysis of the effects of many features across studies, as metafor computes it.
meta_analysis() combines, for every feature (gene), the estimates of one effect
in several studies, each with its standard error, into an estimate of a common effect,
with the between-study variance and the heterogeneity statistics, as metafor’s
rma.uni does for one feature without moderators (Viechtbauer 2010).
differential_expression_meta_analysis() does it for the differential-expression
results of several cubes, such as the same cell types in several cohorts, tissues or
platforms, matched by group, contrast and feature id.
- The model:
Study i estimates the effect of a feature by y_i, with standard error s_i, and its sampling variance v_i = s_i^2 is taken as known. In the random-effects model the true effects of the studies vary around a mean mu with the between-study variance tau^2: y_i = mu + u_i + e_i, with u_i ~ N(0, tau^2) and e_i ~ N(0, v_i) (DerSimonian and Laird 1986). With the weights w_i = 1 / (v_i + tau^2), mu is estimated by mu_hat = sum w_i y_i / sum w_i, with standard error 1 / sqrt(sum w_i). The fixed-effect model, which assumes one true effect (metafor’s equal-effects model,
method = "EE"), is the case tau^2 = 0: the inverse-variance weighted mean.- The between-study variance (
method), over the k studies that measured a feature: "reml", metafor’s default: the restricted maximum likelihood estimate (Viechtbauer 2005), the tau^2 >= 0 that maximizes -1/2 sum log(v_i + tau^2) - 1/2 log(sum w_i) - 1/2 sum w_i (y_i - mu_hat)^2. This function can have more than one maximum, as when an imprecise study lies far from precise ones that agree. Its maxima are zero, where the score is not positive there, and the roots where the score, sum w_i^2 (y_i - mu_hat)^2 - tr(P), with P = W - w w’ / sum w_i for the diagonal matrix W of the weights, falls from positive to negative. The score is negative beyond sum (y_i - y_bar)^2 / (k - 1) + 2 max(v_i). The roots are bracketed by the signs of the score at 8 points per decade of tau^2, from 1e-3 min(v_i) to twice that bound, so that rounding cannot hide one near the bound; each is found by Newton’s method with the observed information, kept in its bracket by bisection, and the estimate is the maximum with the largest restricted likelihood."dl": the DerSimonian and Laird (1986) estimate, the method of moments with the weights 1 / v_i: max(0, (Q - (k - 1)) / (sum w_i - sum w_i^2 / sum w_i)), with Q below."pm": the Paule and Mandel (1982) estimate, at which the generalized Q statistic sum w_i (y_i - mu_hat)^2 equals its expectation k - 1, or 0 where Q is at most k - 1. The statistic falls with tau^2 and is convex, and the root lies between sum (y_i - y_bar)^2 / (k - 1) minus the largest and minus the smallest v_i, so Newton’s method from the lower end converges to it from below."fixed": tau^2 = 0.
- Tests and confidence intervals (
test): "z": z = mu_hat / SE against the standard normal distribution, and the interval mu_hat +/- z_(1 - alpha / 2) SE, as metafor’s defaulttest = "z"."knapp_hartung": the standard error is multiplied by sqrt(q), with q = sum w_i (y_i - mu_hat)^2 / (k - 1), and the statistic is referred to a t distribution with k - 1 degrees of freedom (Hartung 1999; Sidik and Jonkman 2002; Knapp and Hartung 2003; metafor’stest = "knha"), which metafor’s documentation recommends for few studies. Where q < 1, as when the studies agree more closely than their standard errors imply, the adjusted standard error is smaller than the unadjusted one; with the fixed-effect model it is the standard error scaled by sqrt(Q / (k - 1)), which metafor computes with a warning that the method is not meant for that model.
- Heterogeneity:
Cochran’s (1954) Q = sum (y_i - mu_fixed)^2 / v_i, with the fixed-effect estimate mu_fixed, and its p-value from the chi-squared distribution with k - 1 degrees of freedom. I^2 and H^2 are those metafor reports: I^2 = tau^2 / (tau^2 + s^2) and H^2 = (tau^2 + s^2) / s^2, where s^2 = (k - 1) sum w_i / ((sum w_i)^2 - sum w_i^2), with w_i = 1 / v_i, is the typical within-study variance (Higgins and Thompson 2002, equation 9), so that they follow the estimate of tau^2. For the fixed-effect model they are I^2 = max(0, (Q - (k - 1)) / Q) and H^2 = Q / (k - 1), the conventional definitions, which equal the others for DerSimonian-Laird where Q > k - 1; where Q <= k - 1, I^2 is 0 for both, and H^2 is Q / (k - 1) for the fixed effect and 1 for DerSimonian-Laird. I^2 is a proportion here; metafor reports it in percent.
- Studies that did not measure a feature:
NaN marks them, in the effect or the standard error. Each feature is analysed over the k studies that measured it, as metafor leaves out studies with missing values. A feature measured in one study has that study’s estimate and standard error and its z-test, also with
"knapp_hartung", which needs two studies; tau^2 = 0, Q = 0 with a p-value of 1, I^2 = 0 and H^2 = 1: metafor’s results for one study (found by probing). A feature measured in no study has NaN.- Accuracy:
Every estimate is computed for the deviations of the effects from their mean and the standard errors divided, exactly, by a power of two near the largest of them, so that the results do not depend on the scale of the data, and no square overflows or underflows however much larger the spread of the effects is than their standard errors. The sums sum w_i - sum w_i^2 / sum w_i and tr(P^2), which cancel when one study is far more precise than the others, are formed from the sum of the weights of the other studies for each study, which does not cancel.
- Deliberate differences from metafor (
validation/meta_analysis.pyreports them): The REML estimate is the largest maximum of the restricted likelihood. metafor finds a maximum by Fisher scoring from the Hedges estimate, max(0, sum (y_i - y_bar)^2 / (k - 1) - mean(v_i)), and replaces it by zero, with a warning that Fisher scoring “may have gotten stuck at a local maximum”, where the restricted likelihood is larger at zero (found by probing). Where the likelihood has another, larger maximum, metafor’s estimate is the smaller one or zero. Fisher scoring also oscillates without converging where the observed information exceeds twice the expected one: metafor then stops with an error, or, at its default threshold, within the oscillation, and its help page’s remedy
control = list(stepadj = 0.5)converges to the estimate here.The iterative estimates converge to rounding. metafor stops REML when tau^2 changes by less than 1e-5 from one step to the next, and Paule-Mandel at the tolerance of R’s
uniroot, about 1.2e-4 on tau^2;control = list(threshold = 1e-13)andcontrol = list(tol = 1e-14)make it converge to the values here.Where the restricted likelihood is largest at zero although the Hedges estimate is positive, the REML estimate is zero; metafor’s step halving approaches zero and stops near its threshold, at about 1e-6 by default (found by probing).
metafor’s Paule-Mandel search stops with an error when the root lies above max(100, 10 mad(y)^2); here it has no upper limit.
Standard errors must be positive and finite where they are given. For a zero standard error metafor stops (fixed effect, DerSimonian-Laird, Paule-Mandel) or fits without Q, I^2 and H^2 (REML); for an infinite one it stops.
References
Benjamini Y, Hochberg Y (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B 57, 289-300. doi:10.1111/j.2517-6161.1995.tb02031.x
Cochran WG (1954). The combination of estimates from different experiments. Biometrics 10, 101-129. doi:10.2307/3001666
DerSimonian R, Laird N (1986). Meta-analysis in clinical trials. Controlled Clinical Trials 7, 177-188. doi:10.1016/0197-2456(86)90046-2
Hartung J (1999). An alternative method for meta-analysis. Biometrical Journal 41, 901-916. doi:10.1002/(SICI)1521-4036(199912)41:8<901::AID-BIMJ901>3.0.CO;2-W
Higgins JPT, Thompson SG (2002). Quantifying heterogeneity in a meta-analysis. Statistics in Medicine 21, 1539-1558. doi:10.1002/sim.1186
Knapp G, Hartung J (2003). Improved tests for a random effects meta-regression with a single covariate. Statistics in Medicine 22, 2693-2710. doi:10.1002/sim.1482
Paule RC, Mandel J (1982). Consensus values and weighting factors. Journal of Research of the National Bureau of Standards 87, 377-385. doi:10.6028/jres.087.022
Sidik K, Jonkman JN (2002). A simple confidence interval for meta-analysis. Statistics in Medicine 21, 3153-3159. doi:10.1002/sim.1262
Viechtbauer W (2005). Bias and efficiency of meta-analytic variance estimators in the random-effects model. Journal of Educational and Behavioral Statistics 30, 261-293. doi:10.3102/10769986030003261
Viechtbauer W (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software 36(3), 1-48. doi:10.18637/jss.v036.i03
- class cytepool.meta_analysis.MetaAnalysisFit(estimate, standard_error, statistic, p_value, ci_lower, ci_upper, tau_squared, q_statistic, q_p_value, i_squared, h_squared, n_studies, converged, method, test, level)#
Meta-analyses of many features, from
meta_analysis().The per-feature arrays have the shape of the effects without their first (study) axis: one value per feature, or 0-d arrays for the effects of one feature given as a vector. They are read-only.
- Variables:
estimate – The estimate of the common effect mu (metafor’s
beta): the weighted mean of the studies’ effects, with the weights 1 / (v_i + tau^2).standard_error – Its standard error (
se), adjusted by Knapp and Hartung withtest="knapp_hartung".statistic –
estimate / standard_error(zval): a z-statistic, or withtest="knapp_hartung"a t-statistic onn_studies - 1degrees of freedom.p_value – Its two-sided p-value (
pval).ci_lower – The lower limit of the confidence interval of the estimate at
level(ci.lb).ci_upper – The upper limit (
ci.ub).tau_squared – The estimate of the between-study variance tau^2 (
tau2); 0 for the fixed-effect model.q_statistic – Cochran’s Q, the heterogeneity statistic of the fixed-effect fit (
QE).q_p_value – Its p-value, from the chi-squared distribution with
n_studies - 1degrees of freedom (QEp).i_squared – I^2, as a proportion in [0, 1] (metafor’s
I2is in percent).h_squared – H^2 (
H2).n_studies – The number of studies that measured each feature (
k).converged – Whether the iterations of REML or Paule-Mandel met their convergence criterion; True for the other methods.
method –
"reml","dl","pm"or"fixed".test –
"z"or"knapp_hartung".level – The confidence level of the intervals.
Explained in: Meta-analysis
- class cytepool.meta_analysis.MetaAnalysisResult(summary, studies, tables, provenance)#
Meta-analyses of the differential expression of several studies.
Built by
differential_expression_meta_analysis(). Tables hold, for every group and contrast meta-analysed, one row per feature that at leastmin_studiesstudies tested, in the order in which the studies, taken in the order given, first list them:feature_idand, when a study’s tables have one,symbol: the first symbol that a study gives for the feature, taking the studies in the order given;n_studies: the number of studies that tested the feature;log2_fold_changeandstandard_error: the estimate of the common log2 fold change and its standard error;statisticandp_value: the z-statistic, or the Knapp-Hartung t-statistic onn_studies - 1degrees of freedom, and its two-sided p-value. A feature that one study tested, kept withmin_studies=1, has the z-test of that study’s log2 fold change against its standard error, also with"knapp_hartung";adjusted_p_value: Benjamini-Hochberg over the features of the table (Benjamini and Hochberg 1995);ci_lowerandci_upper: the confidence interval of the log2 fold change;tau_squared: the between-study variance;i_squaredandh_squared: I^2, as a proportion, and H^2;q_statisticandq_p_value: Cochran’s Q and its p-value.
The summary has one row per group and contrast:
group_id,contrast,status("analysed";"too_few_studies"when fewer thanmin_studiesstudies tested it;"no_features"when no feature was tested by that many),reason(empty when analysed),n_studies(the studies that tested it),n_features(in its table) andn_features_left_out(tested by fewer thanmin_studiesstudies). The table of studies has one row per group, contrast and study:group_id,contrast,study,tested(whether the study’s results have a table of the group and contrast),reason(why not; empty when tested) andn_features(in the study’s table; 0 when not tested).Immutable.
- Raises:
TypeError – If
summaryorstudiesis not aTable,tablesnot a mapping of (group id, contrast name) toTable, orprovenancenot a mapping that JSON can write.ValueError – If
provenanceholds a NaN or an infinity.
Explained in: Meta-analysis
- property summary#
The status, studies and features of every group and contrast.
- property studies#
Every study of every group and contrast, whether it tested them, and why not.
- property tables#
The table of every group and contrast meta-analysed, by (group, contrast).
- property settings#
The settings of the meta-analysis, as a new dict (also in the provenance).
- property provenance#
A JSON-serialisable record of what produced every number, as a new dict.
It holds the CytePool, numpy, numba and Python versions; the settings; the provenance of every study’s differential-expression result, by study name, which identifies its cube by its fingerprint and gives its design and settings and the samples of every group’s fit; and, for every group and contrast, its status, the studies that tested it and those that did not with the reason, the number of features in its table and left out of it, and the features whose estimate of tau^2 did not converge (
not_converged).
- table(group, contrast)#
Return the table of one group and contrast.
- Parameters:
group – A group id.
contrast – A contrast name.
- Raises:
KeyError – If the group or contrast was not analysed, fewer than
min_studiesstudies tested it, or none of its features; the message says why.
Explained in: Meta-analysis
- cytepool.meta_analysis.differential_expression_meta_analysis(results, *, groups=None, contrasts=None, method='reml', test='z', level=0.95, min_studies=2)#
Meta-analyse the differential expression of several studies, feature by feature.
Each study is the result of
limma_voom()orlimma_trend()on its own cube, such as one cohort, tissue or platform. For every group (cell type) and contrast, the log2 fold changes of the studies that tested it, with their standard errors (the tables’log2_fold_changeandstandard_error), are matched by feature id and meta-analysed bymeta_analysis(). Groups are matched by group id and contrasts by name, so the studies must name them alike, and a contrast should mean the same comparison in every study. A feature that some studies did not test, because their expression filter left it out or their cube does not have it, is analysed over the studies that tested it.The standard errors are the moderated ones of limma’s empirical Bayes moderation, with which the tables’ t-statistics are formed; the meta-analysis takes them as known, as it does any sampling variance, and so does not use their degrees of freedom.
- Parameters:
results – The differential-expression result of every study, by study name; at least two.
groups – A group id or a list of them; by default every group of any study, in the order in which the studies first list them.
contrasts – A contrast name or a list of them; by default every contrast of any study, in the order in which the studies first list them.
method – The estimator of the between-study variance, as for
meta_analysis():"reml","dl","pm"or"fixed".test –
"z"or"knapp_hartung", as formeta_analysis().level – The confidence level of the intervals, in (0, 1).
min_studies – The fewest studies that must test a feature for it to be meta-analysed, from 1 to the number of studies; 2 by default, since one study gives no information about the between-study variance. A group and contrast that fewer studies test is not analysed. With 1, a feature that one study tested has the z-test of that study’s log2 fold change against its standard error, taken as known, as metafor gives it for one study: its p-value is smaller than that of the study’s own moderated t-test, whose degrees of freedom are finite.
- Returns:
The table of every group and contrast meta-analysed, a summary of every group and contrast, the studies that tested them, and a provenance record.
- Warns:
RuntimeWarning – If the estimate of tau^2 of a feature did not converge; the message names the group, the contrast and the features, which the record lists too.
- Raises:
TypeError – If
resultsis not a mapping of study names toDifferentialExpressionResult, an option has the wrong type, orgroups,contrastsormin_studiesis or holds a numpy masked array.ValueError – If there are fewer than two studies, a group or contrast is listed twice, an option is invalid,
min_studiesexceeds the number of studies, or a study’s table has nostandard_errorcolumn, a log2 fold change that is infinite, or a standard error that is not positive and finite (NaN marks a feature the study did not test).KeyError – If a group or contrast is in no study’s results.
Explained in: Choosing a method, Meta-analysis
- cytepool.meta_analysis.meta_analysis(effects, standard_errors, *, method='reml', test='z', level=0.95)#
Meta-analyse the effects of every feature across studies, as metafor’s rma.uni.
For every feature, the effects of the studies that measured it are combined into an estimate of a common effect by the fixed-effect or a random-effects model, with the between-study variance tau^2, Cochran’s Q and its p-value, I^2 and H^2, as
rma.uni(yi, sei = ..., method = ..., test = ...)gives them; the module description gives the model, the estimators of tau^2 and the tests.- Parameters:
effects –
(n_studies, n_features)estimates of the effect of every feature in every study, such as log fold changes, or(n_studies,)for one feature. Finite, or NaN where a study did not measure a feature.standard_errors – Their standard errors, of the same shape: positive and finite, or NaN where a study did not measure a feature. A study is left out of a feature’s analysis where its effect or its standard error is NaN.
method – The estimator of the between-study variance:
"reml"(restricted maximum likelihood, metafor’s default),"dl"(DerSimonian-Laird),"pm"(Paule-Mandel), or"fixed"for the fixed-effect model.test –
"z"for z-tests and normal confidence intervals, or"knapp_hartung"for the Knapp-Hartung adjustment, with t-tests onn_studies - 1degrees of freedom.level – The confidence level of the intervals, in (0, 1). metafor takes it in percent (
level = 95).
- Returns:
The estimate, standard error, statistic, p-value and confidence interval of every feature, with tau^2, Q and its p-value, I^2, H^2 and the number of studies.
- Raises:
TypeError – If
effectsorstandard_errorsdoes not hold real numbers or is a numpy masked array, or an option has the wrong type.ValueError – If the shapes differ or are not 1-D or 2-D, there is no study, an effect or a standard error is infinite, a standard error is not positive,
methodortestis unknown, orlevelis not in (0, 1).
Explained in: Choosing a method, Meta-analysis