M is for Meta-Analysis

I've blogged before about meta-analysis (some examples here, here, and here), but haven't really gone into detail about what exactly it is. It actually straddles the line between method and analysis. Meta-analysis is a set of procedures and analyses that allow you to take multiple studies on the same topic, and aggregate their results (using different statistical techniques to combine results).

Meta-analysis draws upon many of the different concepts I've covered so far this month. Aggregating across studies increases your sample size, maximizing power and providing a better estimate of the true effect (or set of effects). It's an incredibly time-intensive process, but it is incredibly rewarding and the results are very valuable for helping to understand (and come to a consensus on) an area of research and guide future research on the topic.

First of all, you gather every study you can find on a topic, including studies you ultimately might not include. And when I say every study, I mean every study. Not just journal articles but conference presentations, doctoral dissertations, unpublished studies, etc. Some of it you can find in article databases, but some of it you have to find by reaching out to people who are knowledgeable about an area or who have research published on that topic. You'd be surprised how many of them have another study on a topic they've been unable to publish (what we call the "file drawer problem" and relatedly, "publication bias"). The search then weeding through is a pretty intensive process. It helps to have a really clear idea of what you're looking for, and what aspects of a study might result in it being dropped from the meta-analysis.

Next, you would "code" the studies on different characteristics you think might be important. That is, even if you have very narrow criteria for including a study in your meta-analysis, there are going to be differences in how the study was conducted. Maybe the intervention used was slightly different across studies. Maybe the samples were drawn from college freshmen for some studies and community-dwelling adults in others. You decide which of these characteristics are important to examine, then create a coding scheme to pull that information from the articles. To make sure your coding scheme is clear, you'd want to have another person code independently with the same scheme and see if you get the same results. (Yes, this is one of the times I used Cohen's kappa in my research.)

You would use the results of the study (the means/standard deviation, statistical analyses, etc.) to generate an effect size (or effect sizes) for the study. I'll talk more about this later, but basically an effect size allows you to take the results of the study and convert it to a standard metric. Even if the different studies you included in the meta-analysis examined the data in different ways, you can find a common metric so you can compare across studies. At this point, you might average these effect sizes together (using a weighted average - so studies with more people have more impact on the average than studies with fewer people), or you might use some of the characteristics you coded for to see if they have any impact on the effect size.

This is just an overview, of course. I could probably teach a full semester course on meta-analysis. (In fact, that's something I would love to do, since meta-analysis is one of my areas of expertise.) They're a lot of work, but also lots of fun: you get to read and code studies (don't ask me why but this is something I really enjoy doing), and you end up with tons of data to analyze (ditto). If you're interested in learning more about meta-analysis, I recommend starting with this incredible book:

It's a really straightforward, step-by-step approach to conducting a meta-analysis (giving attention to the statistical aspect but mostly focusing on the methods). For a more thorough introduction to the different statistical analyses you can conduct for meta-analysis, I highly recommend the work of Michael Borenstein.

How Statisticians Solve Disagreements

I'm currently taking an online course on meta-analysis, which is a set of statistical and methodological techniques that allow you to combine multiple studies on a topic and generate an estimate (or set of estimates) about the true effect. It's almost like crowd-sourcing data - you're taking advantage of all the work others have done and capitalizing on the strength of having an increased number of participants, difference treatment methods, and so on. I did a candidacy exam in grad school on meta-analysis, and have conducted one before (on pretrial publicity effects), so I know a bit about the topic. This course is devoted to using the R Statistical Package, an open-source program with powerful analysis and graphing capabilities, to conduct a meta-analysis.

For the first week, we were assigned to read up on the R packages we'll be using, as well as an article from the creator of meta-analysis, Gene Glass. I've read some of Glass's work before, but for some reason, didn't encounter this article until now, which tells the reason meta-analysis was created. In addition to wanting to contribute to the field, and have a good topic to introduce in his Presidential Address to the American Educational Research Association, it was really developed to solve a disagreement.

Glass, like many grad students, left grad school with a brand new PhD and a case of depression. He found his way into psychotherapy and was so pleased with his progress, he began studying clinical psychology and became psychotherapy's biggest fan. However, another researcher, Hans Eysenck, became psychotherapy's biggest critic, constantly arguing that any effects were merely placebo:

I found this conclusion personally threatening—it called into question not only the preoccupation of about a decade of my life but my scholarly judgment (and the wisdom of having dropped a fair chunk of change) as well. I read Eysenck's literature reviews and was impressed primarily with their arbitrariness, idiosyncrasy and high-handed dismissiveness. I wanted to take on Eysenck and show that he was wrong: psychotherapy does change lives and make them better.

Glass goes through the decisions Eysenck made in conducting his literature review on the subject, and it's easy to see why, based on these decisions, Eysenck concluded psychotherapy was ineffective - or rather, it easy to see that because Eysenck strongly believed going in that psychotherapy was ineffective, he looked for evidence that supported and ignored evidence that refuted his conclusion. First, he refused to include any research that was not published in a peer reviewed journal, even studies that have to undergo another form of peer review, such as dissertations, theses, or conference presentations. But there is much reason to believe that peer reviewed articles could be biased.

Next, he eliminated any study that didn't have a control group (a group that received no treatment). So if a study compared two forms of therapy, it was tossed out. This left only 11 studies. He then did a vote count, which involves tallying up the number of studies finding a significant difference and the number finding no significant difference. "All that Eysenck considered worth noting about an experiment was whether the differences reached significance at the .05 level. If it reached significance at only the .07 level, Eysenck classified it as showing 'no effect for psychotherapy.'"

And finally, here's the real gem: if he didn't like the outcome they used (that is, he considered it subjective), he discounted the finding, and if a study found differences for one outcome but not for a second one, he also discounted it, calling it "inconsistent." This was the case even if one of the outcomes was something that might be only show a small change due to therapy, such as GPA, versus an outcome that would show a big difference, such as a measure of symptom severity. Eysenck's review didn't even take into account effect sizes: what outcomes would show big differences after psychotherapy and what would show small difference.

And that's where meta-analysis comes in:

Looking back on it, I can almost credit Eysenck with the invention of meta-analysis by anti-thesis. By doing everything in the opposite way that he did, one would have been led straight to meta-analysis. Adopt an a posteriori attitude toward including studies in a synthesis, replace statistical significance by measures of strength of relationship or effect, and view the entire task of integration as a problem in data analysis where "studies" are quantified and the resulting data-base subjected to statistical analysis, and meta-analysis assumes its first formulation. (Thank you, Professor Eysenck.)

So the TL;DR is, how to statisticians solve disagreements? They create new statistics, and then publish pithy articles where they thank the person they disagreed with. Love. It.

E is for Effect Sizes

Today is the first in a three-part series this month on how to conduct meta-analysis using R, plus a fourth Statistics Sunday post tying it all together. As I mentioned in yesterday's post, we'll be using the metafor package. If you didn't install the package then, you'll want to install it now.

install.packages("metafor")

## Installing package into '\\marge/users$/slocatelli/My Documents/R/win-library/3.4'
## (as 'lib' is unspecified)

library(metafor)

## Loading required package: Matrix

## Loading 'metafor' package (version 2.0-0). For an overview 
## and introduction to the package please type: help(metafor).

For this post, we'll focus on computing effect sizes from the summary data. If this sounds like Greek to you, or if you need a refresher, you'll want to first review this post on meta-analysis and this post introducing effect sizes. But to briefly summarize, to conduct a meta-analysis, you're taking information from individual studies to convert into an effect size. You collect these from multiple studies on the same topic to try to estimate the true effect size - the value that all of these individual studies, using different methods, samples, and (yes) flaws are trying to estimate. Later on, we'll go into how you can aggregate those individual study effect sizes.

Let's assume, since all studies are looking at the same topic, that they report the same kind of summary data that we convert into an effect size. In practice, this isn't always true, which is when you have to get into converting between effect sizes. But that's an advanced topic.

Meta-analyses frequently use one of three types of effect sizes: correlation, standardized mean difference (which requires mean and standard deviations for the 2 groups or time points being compared), or odds ratio (or similar metric used for binary data). In some medical meta-analysis, you might also look at event counts over time, such as number of strokes occurring in a set period of time. If you're meta-analyzing correlations, studies often give you your effect size directly, though there are corrections you may want to apply. But if you're meta-analyzing mean differences, ratios, or counts, you have to pull in the summary statistics and convert those to your effect size.

The metafor package will do these calculations easily, using the escalc (effect size calculation) function. You'll state in the function what kind of measure you want, and this will determine what data is needed.

Standardized Mean Difference

If you're computing a standardized mean difference, abbreviated as SMD in metafor, you'll need your study data file to include: mean for each group, standard deviation for each group, and sample size per group. In grad school, I did a meta-analysis on pretrial publicity. Though most studies used guilty/not guilty verdicts as their study outcome, a handful used guilt ratings. Just for fun, I pulled some of those studies out of my study dataset. Here's the data you need to create a data frame we can analyze in metafor:

smd_meta<-data.frame(
  id = c("005","005","029","031","038","041","041","058","058","067","067"),
  study = c(1,2,3,1,1,1,2,1,2,1,2),
  author_year = c("Ruva 2007","Ruva 2007","Chrzanowski 2006","Studebaker 2000",
                  "Ruva 2008","Bradshaw 2007","Bradshaw 2007","Wilson 1998",
                  "Wilson 1998","Locatelli 2011","Locatelli 2011"),
  n1 = c(138,140,144,21,54,78,92,31,29,90,181),
  n2 = c(138,142,234,21,52,20,18,15,13,29,53),
  m1 = c(5.29,5.05,1.97,5.95,5.07,6.22,5.47,6.13,5.69,4.81,4.83),
  m2 = c(4.08,3.89,2.45,3.67,3.96,5.75,4.89,3.80,3.61,4.61,4.51),
  sd1 = c(1.65,1.50,1.08,1.02,1.65,2.53,2.31,2.51,2.51,1.20,1.19),
  sd2 = c(1.67,1.61,1.22,1.20,1.76,2.17,2.59,2.68,2.78,1.39,1.34)
)

ID is a number I assigned to keep track of every source I examined for the meta-analysis, including ones I didn't end up using. Study refers to the study number within the source, since some sources had multiple studies. I used the study number the authors gave, so I could easily refer back to the source if necessary. Not all studies qualified for the meta-analysis, so you'll notice some numbers are skipped.

You may also notice that study ID 067 is mine. Since the meta-analysis, originally done in 2009, was on the same topic as my dissertation, which I completed in 2011, I reran the meta-analysis for one of my dissertation chapters, adding in the data from my own dissertation. In the data above, the group 1s are the treatment group, who saw pretrial publicity, and group 2s are the control group, who did not. This is important to remember when it comes time to interpret the direction of the effect sizes. A positive value means the treatment group gave a higher guilt rating. While each study may use different scales for the guilt rating, since we're standardizing, these differences don't matter. The mean differences are in standard deviation units.

We now apply the escalc function for standardized mean difference, referencing the smd_meta data frame so it appends effect sizes (as a column called "yi") to the data:

smd_meta <- escalc(measure="SMD", m1i=m1, m2i=m2, sd1i=sd1, sd2i=sd2, n1i=n1, n2i=n2,
                   data=smd_meta)

Odds Ratio

Next, we can conduct the same kind of effect size calculation for odds ratio, which is what I used for most of the pretrial publicity studies. As I said, most used verdicts, so each study had a 2x2 table of results:

Once again, here's a data frame pulled from my original meta-analysis dataset:

or_meta<-data.frame(
  id = c("001","003","005","005","011","016","025","025","035","039","045","064","064"),
  study = c(1,5,1,2,1,1,1,2,1,1,1,1,2),
  author_year = c("Bruschke 1999","Finkelstein 1995","Ruva 2007","Ruva 2007",
                  "Freedman 1996","Keelen 1979","Davis 1986","Davis 1986",
                  "Padawer-Singer 1974","Eimermann 1971","Jacquin 2001",
                  "Ruva 2006","Ruva 2006"),
  tg = c(58,26,67,90,36,37,17,17,47,15,133,68,53),
  cg = c(49,39,22,50,12,33,19,17,33,11,207,29,44),
  tn = c(72,60,138,140,99,120,60,55,60,40,136,87,74),
  cn = c(62,90,138,142,54,120,52,57,60,44,228,83,73)
)

I provided guilty counts by group (tg = treatment guilty verdicts, cg = control guilty verdicts), as well as sample sizes per group, which I can use to get not guilty verdicts. We can now request odds ratios (note: metafor gives log odds ratios) with the escalc function:

or_meta <- escalc(measure="OR", ai=tg, bi=(tn-tg), ci=cg, di=(cn-cg), data=or_meta)

But I could request other types of binary effect sizes, such as risk ratios (once again, log-transformed automatically), or risk difference. The metafor package gives more information on the different measures you can request, how they're interpreted, and some sample datasets using the different measures.

In part 2, we'll talk about variances, and in part 3, weights, with a Statistics Sunday post talking about using the results from escalc to generate an aggregate effect size, one of the main goals of meta-analysis. Check back later this month for those posts! And tomorrow, check back for a post on conducting confirmatory factor analysis in R, once again using the Facebook dataset.

W is for (Meta-Analysis) Weights

Yesterday, I talked about the variance calculation for meta-analysis, which is based on sample size - the exact calculation varies by the type of effect size you're using. It's a good idea to inspect those variance calculations. There are many places where your numbers for the meta-analysis can be incorrect.

One place is when extracting that information from the studies themselves, which is why having more than one person extract information and calculating some sort of interrater reliability is a good idea. Another place is at data entry. You should always double-check those numbers, either with a visual examination of the dataset or, my preference especially for large meta-analyses, through data visualization, as I did yesterday with ggplot. Values that don't show an inverse relationship between sample size and variance should be checked. This was the case with study 045 in the or_meta file. A visual check, and knowledge of the formula for variance, confirmed that this variance was correct; it was large because there were a lot more people in the control group than treatment group.

Once we're satisfied that our numbers are correct, we can move forward with the meta-analysis. Part of that calculation involves study weights. Weights and variance have an inverse relationship; in fact, one value can be used to compute the other:

w_i = 1/v_i
v_i = 1/w_i

When you take the inverse of a small number, such as variance, you get a large number, weight. And when you take the inverse of a large number, you get a small number. metafor calculates those weights for you once you conduct the actual meta-analysis, using the rma function. So let's do that. First, I'll load my data frames, and calculate effect sizes and variances.

smd_meta<-data.frame(
  id = c("005","005","029","031","038","041","041","058","058","067","067"),
  study = c(1,2,3,1,1,1,2,1,2,1,2),
  author_year = c("Ruva 2007","Ruva 2007","Chrzanowski 2006","Studebaker 2000",
                  "Ruva 2008","Bradshaw 2007","Bradshaw 2007","Wilson 1998",
                  "Wilson 1998","Locatelli 2011","Locatelli 2011"),
  n1 = c(138,140,144,21,54,78,92,31,29,90,181),
  n2 = c(138,142,234,21,52,20,18,15,13,29,53),
  m1 = c(5.29,5.05,1.97,5.95,5.07,6.22,5.47,6.13,5.69,4.81,4.83),
  m2 = c(4.08,3.89,2.45,3.67,3.96,5.75,4.89,3.80,3.61,4.61,4.51),
  sd1 = c(1.65,1.50,1.08,1.02,1.65,2.53,2.31,2.51,2.51,1.20,1.19),
  sd2 = c(1.67,1.61,1.22,1.20,1.76,2.17,2.59,2.68,2.78,1.39,1.34)
)

or_meta<-data.frame(
  id = c("001","003","005","005","011","016","025","025","035","039","045","064","064"),
  study = c(1,5,1,2,1,1,1,2,1,1,1,1,2),
  author_year = c("Bruschke 1999","Finkelstein 1995","Ruva 2007","Ruva 2007",
                  "Freedman 1996","Keelen 1979","Davis 1986","Davis 1986",
                  "Padawer-Singer 1974","Eimermann 1971","Jacquin 2001",
                  "Ruva 2006","Ruva 2006"),
  tg = c(58,26,67,90,36,37,17,17,47,15,133,68,53),
  cg = c(49,39,22,50,12,33,19,17,33,11,207,29,44),
  tn = c(72,60,138,140,99,120,60,55,60,40,136,87,74),
  cn = c(62,90,138,142,54,120,52,57,60,44,228,83,73)
)

library(metafor)

## Loading required package: Matrix

## Loading 'metafor' package (version 2.0-0). For an overview 
## and introduction to the package please type: help(metafor).

smd_meta <- escalc(measure="SMD", m1i=m1, m2i=m2, sd1i=sd1, sd2i=sd2, n1i=n1,
                   n2i=n2,data=smd_meta)
or_meta <- escalc(measure="OR", ai=tg, bi=(tn-tg), ci=cg, di=(cn-cg),
                  data=or_meta)

These two datasets use different effect size metrics, so I can't combine them unless I convert the effect sizes from one dataset. For now, we'll meta-analyze them separately. First up will be the dataset using standardized mean difference (Cohen's d).

There are two overall methods of meta-analysis - fixed effects and random effects. Fixed effects means that you suspect there is only one underlying effect size that your studies are trying to estimate. You can access this by setting method="FE" (for fixed effects) in the rma function.

smd.rma<-rma(yi,vi,method="FE",data=smd_meta)
summary(smd.rma)

## 
## Fixed-Effects Model (k = 11)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -41.4601   97.3274   84.9202   85.3181   85.3647  
## 
## Test for Heterogeneity: 
## Q(df = 10) = 97.3274, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub     
##   0.3492  0.0526  6.6388  <.0001  0.2461  0.4523  ***
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

I'll go into the various parts of this output on Sunday. For now, we'll look at the Model Results, which includes the estimate (the overall effect size).

smd.rma$beta

##              [,1]
## intrcpt 0.3491923

The overall Cohen's d for this set of studies is 0.349, a moderate effect. In plain English, the average guilt rating among people who saw pretrial publicity is about 0.349 standard deviations higher than people in the control group. There's an associated Z-value and p-value, which tests whether that metric is significantly different from 0.

options(scipen=999)
smd.rma$zval

## [1] 6.638784

smd.rma$pval

## [1] 0.00000000003162821

The Z-test confirms that this value is significantly different from 0. I'll go into conducting a random effects model (and how to tell if you should) on Sunday. Now, let's run the meta-analysis for the log odds ratio dataset.

or.rma<-rma(yi,vi,method="FE",data=or_meta)
summary(or.rma)

## 
## Fixed-Effects Model (k = 13)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -23.0842   47.8107   48.1684   48.7334   48.5321  
## 
## Test for Heterogeneity: 
## Q(df = 12) = 47.8107, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub     
##   0.7456  0.0987  7.5566  <.0001  0.5522  0.9389  ***
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

or.rma$beta

##              [,1]
## intrcpt 0.7455659

or.rma$zval

## [1] 7.556647

or.rma$pval

## [1] 0.00000000000004135935

Our estimated log odds ratio for this dataset is 0.7456, which is significantly different from 0. This number is difficult to interpret in isolation, so we want to convert it back to odds ratio, which is easier to interpret.

exp(or.rma$beta)

##             [,1]
## intrcpt 2.107634

Overall, people who saw pretrial publicity were over 2 times as likely to convict as people in the control group.

metafor will also let you supply user-defined weights, using the weights argument in the rma function. If you're applying a correction to your effect sizes - such as converting your Pearson's r coefficients to Fisher's Z - you'll probably want to supply custom weights, because the formula is different. If you leave that argument out, as I did, rma will create weights as the inverse of variance by default.

I should note that I only picked a handful of studies to include for this example, so these don't match the results I found in my actual meta-analysis, which included many more studies. (If you're conducting research on pretrial publicity and would like the actual results of my meta-analysis, please contact me. You can find contact info on my about page. I'm happy to share my actual results. These numbers are for demonstration purposes only. I still hope to publish this work someday, which is why I'm not using the full dataset or actual results here, but at this point, I need to collect more data on the new pretrial publicity studies that have been conducted since 2011.)

More on meta-analysis Sunday! And tune in tomorrow for the letter X!

Statistics Sunday: Converting Between Effect Sizes for Meta-Analysis

I'm currently working on my promised video on mixed effects meta-analysis, and was planning on covering this particular topic in that video - converting between effect sizes. But I decided to do this as a separate post that I can reference in the video, which I hope to post next week.

As a brief refresher, meta-analysis is aimed at estimating the true effect (or effects) in an area of study by combining findings from multiple studies on that topic. Effect sizes, the most frequently used being Cohen's d, Pearson's r, and log odds ratio, are estimated from information in study reports and presentations. There's a lot of variation in how clearly reports and documents describe the findings and the information given to estimate the study's overall effect. But when you conduct the meta-analysis, whether using fixed, random, or mixed effects analysis, you need to use only one type of effect size. That means that, sometimes, studies will give you a different type of effect size than you plan to use. Fortunately, there are ways to convert between effect sizes and use different types of statistical information to generate your estimates.

First up, converting between those key effect sizes. In the meta-analysis I performed in grad school, I examined the effect of pretrial publicity on guilt. There are two ways guilt was frequently operationalized in the studies: as a guilty/not guilty verdict or as a continuous guilt rating. For those outcomes, we would likely use, respectively, log odds ratio and Cohen's d. The escalc function in the metafor package can compute log odds ratio for guilty/not guilty counts, and Cohen's d for mean and standard deviation of the guilt ratings. But studies may use different types of information when presenting their results, so you may not be able to simply compute those effect sizes.

For instance, a study using verdict may present a chi-square and one of its effect sizes, Cramer's V, which is very similar to a correlation coefficient. How can I convert that into log odds ratio?

To convert from one effect size to the other, you need to follow a prescribed path, which can be seen in the diagram below. What this diagram tells you is which effect sizes you can convert between directly: you can directly convert between log odds ratio and Cohen's d, and between Cohen's d and Pearson's r. If you wanted to convert between Pearson's r and log odds ratio, you'll first need to convert to Cohen's d. You'll need to do the same thing for variance - compute it for the native effect size metric, then convert that to the new effect size metric.

Let's start by setting up functions that will convert between our effect sizes for us, beginning with Cohen's d and log odds ratio. Then we'll demonstrate with some real data.

#Convert log odds ratio to d
ltod <- function(lor) {
  d = lor * (sqrt(3)/pi)
  return(d)
}
vltovd <- function(vl) {
  vd = vl * (3/pi^2)
  return(vd)
}

#Convert d to log odds ratio
dtol <- function(d) {
  lor = d*(pi/sqrt(3))
  return(lor)
}
vdtovl <- function(vd) {
  vl = vd*(pi^2/3)
  return(vl)
}

You'll notice a mathematical symmetry in these equations - the numerators and denominators switch between the equations. Now let's set up equations to r and d. These equations are slightly more complex and will require a few additional arguments. For instance, converting the variance of r to variance of d requires both the variance of r and r itself. Converting from d to r requires group sample sizes, referred to as n1 and n2.

#Convert r to d
rtod <- function(r) {
  d = (2*r)/(sqrt(1-r^2))
  return(d)
}
vrtovd <- function(vr,r) {
  vd = (4*vr)/(1-r^2)^3
  return(vd)
}

#Convert d to r
dtor <- function(n1,n2,d) {
  a = (n1+n2)^2/(n1*n2)
  r = d/(sqrt(d^2+a))
  return(r)
}
vdtovr <- function(n1,n2,vd,d) {
  a = (n1+n2)^2/(n1*n2)
  vr = a^2*vd/(d^2+a)^3
  return(vr)
}

Remember that the metafor package can compute effect sizes and variances for you, so you might want to run the escalc on the native effect sizes so that you have the estimates and variances you need to run these functions. But if you ever find yourself having to compute those variances by hand, here are the equations, which we'll use in the next step.

vard <- function(n1,n2,d) {
  vd = ((n1+n2)/(n1*n2)) + (d^2/(2*(n1+n2)))
  return(vd)
}

varr <- function(r,n) {
  vr = (1-r^2)^2/(n-1)
  return(vr)
}

varlor <- function(a,b,c,d) {
  vl = (1/a)+(1/b)+(1/c)+(1/d)
  return(vl)
}

One of the studies I included in my meta-analysis gave Cramer's V. It had a sample size of 42, with 21 people in each group. I'd like to convert that effect size to log odds ratio. Here's how I could do it.

cramerv <- 0.67
studyd <- rtod(cramerv)
studyvr <- varr(0.67,42)
studyvd <- vrtovd(studyvr,cramerv)
dtol(studyd)

## [1] 3.274001

vdtovl(studyvd)

## [1] 0.5824038

I can now include this study in my meta-analysis of log odds ratios.

What if my study gives different information? For instance, it might have given me a chi-square or a t-value. This online effect size calculator, created by David Wilson, coauthor of Practical Meta-Analysis, can compute effect sizes for you from many different types of information. In fact, spoiler alert: I used an earlier version of this calculator extensively for my meta-analysis. Note that this calculator returns odds ratios, so you'll need to convert those values into a log odds ratio.

Free Tools for Meta-Analysis

My boss is attending a two-day course on meta-analysis, and shared these tools with me, available through Brown School of Health:

• The Systematic Review Data Repository - as the name suggests, this is a repository of systematic review data, so you pull out data relevant to your own systematic review as well as contribute your own data for others to use. Systematic reviews are a lot of work, so a tool that lets you build off of the work of others can help systematic reviews be performed (and their findings disseminated and used to make data-driven decisions) much more quickly
• Abstrackr - a free, open-source tool for the citation screening process. Conducting a systematic review or meta-analysis involves an exhaustive literature review, and those citations then have to be inspected to see if they qualify to be included in the study. It isn't unusual to review 100s of studies only to include a couple dozen (or fewer). This tool lets you upload abstracts, and invite reviewers to examine abstracts for inclusion. This tool is still in beta, but they're incorporating machine learning to automate some of the screening process in the future. Plus, they use "automagically" in the description, which is one of my favorite portmanteaus.
• Open Meta-Analyst - another free, open-source tool for conducting meta-analysis. You can work with different types of data (binary, continuous, diagnostic), conduct fixed- or random-effects models, and even use different estimation methods, like maximum likelihood or Bayesian. 
• Open MEE - a free, open-source tool based on Open Meta-Analyst, with extra tools for ecological and evolutionary meta-analysis. This might be the tool to use in general, because it has the ability to conduct meta-regression with multiple covariates. 

I think of all of these, I'm looking forward to trying out Abstrackr the most.

And of course, there are many great meta-analysis packages for R. I'm currently working on a methods article describing how to conduct a mini meta-analysis to inform a power analysis using R tools - something I did for my dissertation, but not something everyone knows how to do. (By working, I mean I have an outline and a few paragraphs written. But I'm hoping to have more time to dedicate to it in the near future. I'm toying with the idea of spending NaNoWriMo this year on scholarly pursuits, rather than a novel.)

BTW, if you like free stuff, check out these free data science and statistics resources (and let me know if you know of any not on the list).

Meta-Analysis and Reproducibility: More from Sara's Week in Psychological Science

Yesterday, I attended a symposium on meta-analysis, which is a set of methods and statistical techniques to combine results of multiple studies. The purpose of a meta-analysis is to determine the true effect of some phenomenon, and also to understand what characteristics of samples, studies, and so on impact study results. I did one of my grad school candidacy exams on meta-analysis, and have a meta-analysis I've been working on and trying to get published.

During yesterday's session, one of the presenters talked about an Open Science Framework project, in which the researchers are trying to reproduce findings from multiple meta-analyses in psychology. The issue is that, as is true with many psychological studies, people often have difficulty reproducing findings from past meta-analyses, in part because of the strong impact of subjective decisions when conducting a meta-analysis. That is, decisions like what types of studies to include, how to divide up data within studies, and the specific analysis technique can be highly subjective and left up to the discretion of the researcher.

This project will involve meta-analyzing 2000 studies. Obviously, they are going to need help, and if you're so inclined, you can volunteer to be part of the project.

V is for (Meta-Analysis) Variance

For the letter E, I introduced the metafor package to compute effect sizes. That is, you provide a data frame with the study information and the data needed to compute the effect size(s), and metafor does that for you. But what I didn't mention is that metafor also gives you additional information that will become very important when conducting an actual meta-analysis - and that information will be the focus of today's and tomorrow's posts.

Today, I'll talk about the concept of variance in meta-analysis and discuss how metafor gives you that information. But first, we should set up our data frames with study information and data. I'll be reusing the two data frames I created for the E post - one to compute standardized mean difference (Cohen's d) and one to compute log odds ratio.

smd_meta<-data.frame(
  id = c("005","005","029","031","038","041","041","058","058","067","067"),
  study = c(1,2,3,1,1,1,2,1,2,1,2),
  author_year = c("Ruva 2007","Ruva 2007","Chrzanowski 2006","Studebaker 2000",
                  "Ruva 2008","Bradshaw 2007","Bradshaw 2007","Wilson 1998",
                  "Wilson 1998","Locatelli 2011","Locatelli 2011"),
  n1 = c(138,140,144,21,54,78,92,31,29,90,181),
  n2 = c(138,142,234,21,52,20,18,15,13,29,53),
  m1 = c(5.29,5.05,1.97,5.95,5.07,6.22,5.47,6.13,5.69,4.81,4.83),
  m2 = c(4.08,3.89,2.45,3.67,3.96,5.75,4.89,3.80,3.61,4.61,4.51),
  sd1 = c(1.65,1.50,1.08,1.02,1.65,2.53,2.31,2.51,2.51,1.20,1.19),
  sd2 = c(1.67,1.61,1.22,1.20,1.76,2.17,2.59,2.68,2.78,1.39,1.34)
)

or_meta<-data.frame(
  id = c("001","003","005","005","011","016","025","025","035","039","045","064","064"),
  study = c(1,5,1,2,1,1,1,2,1,1,1,1,2),
  author_year = c("Bruschke 1999","Finkelstein 1995","Ruva 2007","Ruva 2007",
                  "Freedman 1996","Keelen 1979","Davis 1986","Davis 1986",
                  "Padawer-Singer 1974","Eimermann 1971","Jacquin 2001",
                  "Ruva 2006","Ruva 2006"),
  tg = c(58,26,67,90,36,37,17,17,47,15,133,68,53),
  cg = c(49,39,22,50,12,33,19,17,33,11,207,29,44),
  tn = c(72,60,138,140,99,120,60,55,60,40,136,87,74),
  cn = c(62,90,138,142,54,120,52,57,60,44,228,83,73)
)

Now we'll rerun the code from that previous post that we used to generate effect sizes. We reference the original data frame when we run this code if we want the computed variables to be appended to the existing data frame.

library(metafor)

## Loading required package: Matrix

## Loading 'metafor' package (version 2.0-0). For an overview 
## and introduction to the package please type: help(metafor).

smd_meta <- escalc(measure="SMD", m1i=m1, m2i=m2, sd1i=sd1, sd2i=sd2, n1i=n1,
                   n2i=n2,data=smd_meta)
or_meta <- escalc(measure="OR", ai=tg, bi=(tn-tg), ci=cg, di=(cn-cg),
                  data=or_meta)

Let's take a quick look at the updated data frame, starting with smd_meta.

smd_meta

##     id study      author_year  n1  n2   m1   m2  sd1  sd2      yi     vi
## 1  005     1        Ruva 2007 138 138 5.29 4.08 1.65 1.67  0.7269 0.0154
## 2  005     2        Ruva 2007 140 142 5.05 3.89 1.50 1.61  0.7433 0.0152
## 3  029     3 Chrzanowski 2006 144 234 1.97 2.45 1.08 1.22 -0.4099 0.0114
## 4  031     1  Studebaker 2000  21  21 5.95 3.67 1.02 1.20  2.0087 0.1433
## 5  038     1        Ruva 2008  54  52 5.07 3.96 1.65 1.76  0.6464 0.0397
## 6  041     1    Bradshaw 2007  78  20 6.22 5.75 2.53 2.17  0.1893 0.0630
## 7  041     2    Bradshaw 2007  92  18 5.47 4.89 2.31 2.59  0.2444 0.0667
## 8  058     1      Wilson 1998  31  15 6.13 3.80 2.51 2.68  0.8927 0.1076
## 9  058     2      Wilson 1998  29  13 5.69 3.61 2.51 2.78  0.7867 0.1188
## 10 067     1   Locatelli 2011  90  29 4.81 4.61 1.20 1.39  0.1592 0.0457
## 11 067     2   Locatelli 2011 181  53 4.83 4.51 1.19 1.34  0.2603 0.0245

You'll notice that, in addition to the information from the data frame creation code above, I have two new variables: yi and vi. yi refers to the effect size for that study. vi is the variance associated with the effect size for that study. This information is important for computing the study weight, which I'll talk about tomorrow, as well as generating confidence intervals around the overall effect size, which I'll talk about Sunday. But you might be wondering where this variance comes from. For that, stand back - I'm about to use math.

Meta-analysis uses weights based on sample size. The reason for that is because meta-analysis is attempting to take information from multiple studies on a topic and estimate the underlying effect size that these studies are trying to estimate. We weight on sample size because studies done on more people are expected to be better able to estimate the true population value, so larger studies get more weight. On the flip side, larger studies are expected to have more precision in their estimate, so larger studies have smaller variance. In meta-analysis, the variance is estimated based on sample size. For standardized mean difference, the variance is computed as:

v_i = ((n₁+n₂)/(n₁*n₂)) + (y_i²/(2*(n₁+n₂)))

So for the first effect size, which is 0.7269, we can recreate the computed variance (within some rounding error):

((138+138)/(138*138)) + (0.7269/(2*(138+138)))

## [1] 0.0158096

That calculation is done for each effect size, and fortunately, metafor will do it for you. And since sample size is part of the calculation, you'll note that variances are larger for smaller studies and smaller for larger studies.

library(ggplot2)
ggplot(smd_meta, aes(x=vi, y=(n1+n2))) + geom_point() +
  geom_text(aes(label=id),hjust=0, vjust=0) +
  scale_x_continuous(breaks=seq(0,0.15,0.01)) +
  scale_y_continuous(breaks=seq(0,400,25)) +
  labs(x = "Variance", y = "Total Sample Size") + theme_bw()

I did the same calculations for the or_meta data frame. Variance for odds ratio is calculated based on the sample size in each of the four cells (treatment-outcome1 (a), treatment-outcome2 (b), control-outcome1 (c), and control-outcome2 (d)):

v_i = (1/n_a) + (1/n_b) + (1/n_c) + (1/n_d)

So again, larger studies will have smaller variances, though imbalance among those four cells (for instance, having a lot more people in one group than another) will inflate variance somewhat. We can see that for study 045, which has a lot more people in the control group than treatment group.

or_meta

##     id study         author_year  tg  cg  tn  cn      yi     vi
## 1  001     1       Bruschke 1999  58  49  72  62  0.0945 0.1860
## 2  003     5    Finkelstein 1995  26  39  60  90  0.0000 0.1131
## 3  005     1           Ruva 2007  67  22 138 138  1.6046 0.0831
## 4  005     2           Ruva 2007  90  50 140 142  1.1976 0.0620
## 5  011     1       Freedman 1996  36  12  99  54  0.6931 0.1508
## 6  016     1         Keelen 1979  37  33 120 120  0.1615 0.0809
## 7  025     1          Davis 1986  17  19  60  52 -0.3759 0.1650
## 8  025     2          Davis 1986  17  17  55  57  0.0513 0.1690
## 9  035     1 Padawer-Singer 1974  47  33  60  60  1.0845 0.1655
## 10 039     1      Eimermann 1971  15  11  40  44  0.5878 0.2279
## 11 045     1        Jacquin 2001 133 207 136 228  1.5035 0.3933
## 12 064     1           Ruva 2006  68  29  87  83  1.8968 0.1203
## 13 064     2           Ruva 2006  53  44  74  73  0.5089 0.1237

ggplot(or_meta, aes(x=vi, y=(tn+cn))) + geom_point() +
  geom_text(aes(label=id),hjust=0, vjust=0) +
  scale_x_continuous(breaks=seq(0,0.4,0.05)) +
  scale_y_continuous(breaks=seq(0,375,25)) +
  labs(x = "Variance", y = "Total Sample Size") + theme_bw()

Tune in tomorrow, when I'll talk about weights for meta-analysis, which is based on this variance calculation!