Solved – Meta Analysis of Cox Regression Coefficients

cox-modelmeta-analysismeta-regressionsurvival

I'm looking for a bit of guidance on this, as my stats skills are leaving me treating this as a bit more black-boxish than I'm comfortable with.

So I have several datasets where I am employing Cox Regression on survival data and some molecular measurements. I'm using the coxph routine in R from the Survival analysis package. From this, I get a Beta coefficient, standard error, a z-score, p-value, etc.

All of that I'm comfortable with.

What I'd like to do is a meta-analysis over these coefficients from the different datasets. I've been using the Metafor package in R to do this. As input, I've been giving it the beta coefficients and associated standard errors output by the coxph routine.

I'm not entirely sure that the standard errors provided by the coxph are the variances I need to do the meta analysis, or whether I should be using some weighting scheme that take the size of each dataset (number of samples) into effect more explicitly. (e.g. I may have 40 samples in one data set and 700 in another)

thanks for any guidance

Best Answer

Just supply the beta coefficients and corresponding standard errors to the rma() function. So, your syntax should be like this:

rma(coef, sei=se, data=dat)

where coef is the name of the variable in dataset dat denoting the coefficients and se is the name of the variable for the corresponding standard errors. The standard errors already include the information about the number of samples (and actually, it's the number of "events", not the sample sizes, that determine the size of the standard errors).

Related Solutions

Solved – Meta-analysis of trials given hazard ratios

Assuming that you have the log hazard ratios (collected in the vector yi) and the standard errors of the log hazard ratios (collected in the vector sei), then you can fit a fixed-effects model to these data with:

res <- rma(yi, sei=sei, data=dat, method="FE")

As with several other outcome measures typically used for meta-analyses (e.g., risk ratios, odds ratios, incidence rate ratios), we do not meta-analyze the hazard ratios directly, but first log-transform them. This has two purposes:

The sampling distribution of the raw hazard ratio is typically very skewed. However, standard meta-analytic models assume that the sampling distribution of the outcome measure is (at least approximately) normal and that assumption is much more appropriate for the log-transformed hazard ratio.
Raw hazard ratios are not symmetric around 1. So, if one were to average a hazard ratio of 0.5 and 2.0 (which are exact opposites of each other), then we would get 1.25, which makes no sense. However, the log-transformed values are -0.6931 and +0.6931, whose average is 0. After back-transformation (exponentiation), we get 1, and the appropriate conclusion that, on average, the hazard ratio is 1.

Also note that analyses using the Cox proportional hazards model (which is often the source of hazard ratios obtained for a meta-analysis) are done on the log-scale. So, when the standard error is reported together with the hazard ratio, that standard error in all likelihood refers to the log-transformed hazard ratio, not the raw one (in fact, given the skewness of the distribution of the raw hazard ratio, the usefulness of a standard error for the raw hazard ratio, if it were indeed reported, is rather questionable).

Similarly, if a confidence interval is reported, then in all likelihood it was constructed on the log-scale and the bounds were exponentiated afterwards. We can then easily back-calculate the SE of the log hazard ratio. See also this question, which deals exactly with that issue.

The same goes for hazard ratios that are reported in conjunction with the results from the log-rank test. Not surprisingly, as the name of the test implies, the underlying analysis is also done on the log scale. So, a standard error that is computed based on such results again refers to log hazard ratio, not the raw one.

A few additional comments:

The sampling variance of the log hazard ratio is obtained by simply squaring the standard error. So, you could just use rma(yi, sei^2, data=dat, method="FE") or first create vi <- dat$sei^2 and then rma(yi, vi, data=dat, method="FE"). Do not multiply/divide by the sample size or square-root thereof.
The standard error of the log hazard ratio is not the same as taking the log of the standard error of the raw hazard ratio. One cannot just apply the same transformation to the standard error (or variance). In some cases, one can use the delta method to approximate the standard error (or variance) of some random variable following a transformation.

Multilevel Meta-Regression – Detailed Insights into Multilevel Multivariate Meta-Regression

As you note, the model that adds random effects for each study and random effects for each outcome is a model that accounts for hierarchical dependence. This model allows the true outcomes/effects within a study to be correlated. This is the Konstantopoulos (2011) example you link to.

But this model still assumes that the sampling errors of the observed outcomes/effects within a study are independent, which is definitely not the case when those outcomes are assessed within the same individuals. So, as in the Berkey et al. (1998) example you link to, ideally you need to construct the whole variance-covariance matrix of the sampling errors (with the sampling variances along the diagonal). The chapter by Gleser and Olkin (2009) from the Handbook of research synthesis and meta-analysis describes how the covariances can be computed for various outcomes measures (including standardized mean differences). The analyses/methods from that chapter are replicated here (you are dealing with the multiple-endpoint case).

And as you note, doing this requires knowing how the actual measurements within studies are correlated. Using your example, you would need to know for study 1 how strong the correlation was between the two measurements for "Phonological loop" (more accurately, there are two correlations, one for the first and one for the second group, but we typically assume that the correlation is the same for the two groups), and how strongly those measurements were correlated with the "Central Executive" measurements. So, three correlations in total.

Obtaining/extracting these correlations is often difficult, if not impossible (as they are often not reported). If you really cannot obtain them (even after contacting study authors in an attempt to obtain the missing information), there are several options:

One can still often make a rough/educated guess how large the correlations are. Then we use those 'guestimates' and conduct sensitivity analyses to ensure that conclusions remain unchanged when the values are varied within a reasonable range.
One could use robust methods -- in essence, we then consider the assumed variance-covariance matrix of the sampling errors to be misspecified (i.e., we assume it is diagonal, when in fact we know it isn't) and then estimate the variance-covariance matrix of the fixed effects (which are typically of primary interest) using consistent methods even under such a model misspecification. This is in essence the approach described by Hedges, Tipton, and Johnson (2010) that you mentioned.
Resampling methods (i.e., bootstrapping and permutation testing) may also work.
There are also some alternative models that try to circumvent the problem by means of some simplification of the model. Specifically, in the model/approach by Riley and colleagues (see, for example: Riley, Abrams, Lambert, Sutton, & Thompson, 2007, Statistics in Medicine, 26, 78-97), we assume that the correlation among the sampling errors is identical to the correlation among the underlying true effects, and then we just estimate that one correlation. This can work, but whether it does depends on how well that simplification matches up with reality.
There is always another option: Avoid any kind of statistical dependence via data reduction (e.g., selecting only one estimate, conducting separate analyses for different outcomes). This is still the most commonly used approach for 'handling' the problem, because it allows practitioners to stick to (relatively simple) models/methods/software they are already familiar with. But this approach can be wasteful and limits inference (e.g., if we conduct two separate meta-analyses for outcomes A and B, we cannot test whether the estimated effect is different for A and B unless we can again properly account for their covariance).

Note: The same issue was discussed on the R-sig-mixed-models mailing list and in essence I am repeating what I already posted there. See here.

For the robust method, you could try the robumeta package. If you want to stick to metafor, you will find these, blog, posts by James Pustejovsky of interest. He is also working on another package, called clubSandwich which adds some additional small-sample corrections. You can also try the development version of metafor (see here) -- it includes a new function called robust() which you can use after you have fitted your model to obtain cluster robust tests and confidence intervals. And you can find some code to get you started with bootstrapping here.

Best Answer

Related Solutions

Solved – Meta-analysis of trials given hazard ratios

Multilevel Meta-Regression – Detailed Insights into Multilevel Multivariate Meta-Regression

Related Question