CAMARADES Berlin Meta-Analysis Tutorial

Welcome

This web-based tutorial complements the workshop: Introduction to Systematic Review & Meta-Analysis of Animal Studies

We will use the R programming language to work through a meta-analysis of animal data but you do not need to have used R previously to complete the exercises.

For the tutorial, we will use example data from the publication:

McCann SK, Cramond F, Macleod MR et al. Systematic Review and Meta-Analysis of the Efficacy of Interleukin-1 Receptor Antagonist in Animal Models of Stroke: an Update. Translational Stroke Research 7, 395–406 (2016).

Background information on the data set:
Interleukin-1 receptor antagonist (IL-1RA) is an anti-inflammatory protein used clinically to treat rheumatoid arthritis and is considered a promising candidate therapy for stroke. In this systematic review and meta-analysis, we sought to assess the efficacy of IL-1RA in animal models of ischaemic stroke, the range of conditions under which efficacy has been tested, and whether the data appear to be confounded due to reported study quality and publication bias. The outcomes we will assess in this tutorial are post-stroke infarct volume (brain damage), mortality, and neurobehavioural outcomes e.g. motor deficits, in IL-1RA-treated vs. control animals.

Overview

We will analyse the effect of IL-1RA treatment on each of our three outcomes (infarct volume, mortality, and neurobehaviour) separately.

For infarct volume and mortality, all of the R code is provided and you just need to click on the blue Run Code button on the top right of each code chunk to view and interpret the results.

For neurobehaviour, you will have the opportunity to write your own code to analyse the data set provided.

Within each section, there are multiple-choice questions to test your understanding. If you answer a question incorrectly the first time, you can click the orange Try Again button.

The tutorial comprises five main sections and at the end of the first three sections are data dictionaries describing the variables included in each data set:

1. Meta-analysis of continuous data (infarct volume)

In the first section of the tutorial we will cover the steps to meta-analyse continuous data:
1.1 Calculating effect sizes using normalised mean difference (NMD)
1.2 Performing a random effects meta-analysis
1.3 Exploring heterogeneity using univariable meta-regression
1.4 Data dictionary: infarct volume

2. Meta-analysis of binary data (mortality)

In the second section we will look at the binary outcome, mortality, calculate odds ratio effect sizes and perform a random effects meta-analysis.

2.1 Data dictionary: mortality

3. Test your skills: meta-analysis (neurobehaviour)

In this section we provide you with a data set that describes neurobehavioural outcomes (continuous data) from the studies included in the systematic review, along with prompts to insert code to analyse the data.

You can copy and edit the code provided in section 1 to perform the analyses. If you are unsure about what code to use, click on the Hint(s) buttons. You can also compare your results to those in the paper, McCann et al. 2016.

3.1 Data dictionary: neurobehaviour

4. Small study effects (infarct volume)

In this section we will assess the infarct volume data set for the presence of small study effects using funnel plot-based methods including Egger’s regression and trim and fill.

5. Test your skills: small study effects (neurobehaviour)

Lastly, you will have the opportunity to apply your skills in assessing small study effects using a neurobehaviour data set.

For further information on systematic review and meta-analysis of animal studies, please see the CAMARADES wiki page.

1. Meta-analysis of continuous data

First, let’s load the data we will use for the meta-analysis of infarct volume:

infarctVol <- read.csv("IL-1RA_InfarctVolume_section1.csv", 
                       fileEncoding = "UTF-8-BOM")

It’s good practice to take a look at the data you import to make sure there are no errors. After you run the next section of code, the imported data will open in a table below. To explore the table you can use the arrowhead to the right of the column headings to move through the columns and the blue numbers at the bottom right to move through the rows.

You can find a description of each variable in the data dictionary table in section 1.4.

infarctVol %>% 
  mutate_if(is.numeric, round, digits = 2)

Examine the data set and answer the following:

1.1 Calculating effect sizes

Now we will calculate an effect size for each experiment in the data set. An effect size is a measure of the difference between two experimental cohorts.

Infarct volume is a continuous variable and we will calculate a normalised mean difference (NMD) effect size using the control (vehicle-treated) and treatment (IL-1RA-treated) groups. The effect size is expressed as the proportional or percentage improvement in the treatment group compared with the control group. NMD is not widely used in clinical systematic reviews and there is no NMD effect size calculator included in the meta-analysis package we are using. Therefore, we will calculate NMD effect sizes manually.

Note that green text following a hash symbol (#) in the code chunk is not read by R but simply provides a description of the code to make it easier to follow.

First, specify the columns in the data set that contain the mean, standard deviation, and number of animals in the control groups:

# save the group means to a variable called mCi
mCi <- infarctVol$control_group_mean

# save the group standard deviations to a variable called sdCi 
sdCi <- infarctVol$control_group_sd

# save the number of animals to a variable called nCi
nCi <- infarctVol$control_group_calc_n

…and the columns that contain the mean, standard deviation, and number of animals in the treatment groups:

# save the group means to a variable called mTi
mTi <- infarctVol$treatment_group_mean

# save the group standard deviations to a variable called sdTi 
sdTi <- infarctVol$treatment_group_sd 

# save the number of animals to a variable called nTi
nTi <- infarctVol$treatment_group_n

You have just created six new variables that you will use to calculate effect sizes:

Variable	Description
mCi	mean in the control group
sdCi	standard deviation in the control group
nCi	number of animals in the control group
mTi	mean in the treatment group
sdTi	standard deviation in the treatment group
nTi	number of animals in the treatment group

You can view variables to check that they contain the correct data. For example, running the following section of code will show you the list of values in the mCi variable (data for the mean in the control group). This should match the corresponding column in the data table you saw in the previous section (use the navigation bar on the left to go back and check). You can edit this code to view the other variables if you wish:

print(mCi, digits =4)

Now let’s calculate an NMD effect size for each experiment (saved as the variable “ESi”) and add the effect size data as a new column in the infarctVol data set. The function cbind takes a table of data and binds a new column to the right side. Remember, each variable in this equation is a list of numbers, with each entry corresponding to a single experiment.

Note: we are using a truncated form of the NMD equation here because we know that a sham animal that has not suffered a stroke will have no infarct. This reduces the corresponding parts of the NMD equation to zero; see the CAMARADES wiki page for more information.

# calculate the NMD effect size for each experiment and save the results to the ESi variable
ESi <- (mCi - mTi) / mCi * 100

# add the effect size data as a new column
infarctVol <- cbind(infarctVol, ESi)

Next, we calculate the normalised standard deviation for each group:

# control group
nsdCi <- sdCi / mCi * 100 

# treatment group
nsdTi <- sdTi / mCi * 100

We will now calculate the normalised standard error for each experiment (saved as the variable “SEi”) and add the standard error data as a new column in the infarctVol data set:

# calculate the NMD standard error
SEi <- sqrt((nsdCi^2) / nCi + (nsdTi^2) / nTi)

# add the standard error data as a new column in the table
infarctVol <- cbind(infarctVol, SEi)

Check that your new variables have been added as columns to the right hand side of the data set:

Variable	Description
ESi	effect size for each comparison (NMD)
SEi	standard error for each comparison

infarctVol %>% 
  mutate_if(is.numeric, round, digits = 2)

1.2 Random effects meta-analysis

Building the meta-analysis model

Next we combine the individual effect sizes with meta-analysis using the function metagen from the {meta} package. This gives us an overall summary estimate of the effect of IL-1RA on infarct volume, compared to control. Note that to begin with, we will print both the random effects and fixed effect estimates so that we can compare the results:

infarctVolMA <- metagen(
  # specify the variable that contains the effect size for each experiment
  ESi,   
  
  # specify the variable that contains the standard error for each experiment
  SEi,  
  
  # specify the data set
  data = infarctVol,  
  
  # specify the study labels
  studlab = study_id,  
  
  # specify a random effects model
  comb.random = TRUE, 
  
  # specify a fixed effect model
  comb.fixed = TRUE,
  
  # specify which method is used to estimate the between-study variance, τ^2: restricted maximum-likelihood (REML) estimator
  method.tau = "REML"  
)

View the results of the meta-analysis:

infarctVolMA

Answer the following questions about the results of the meta-analysis:

Forest plot

A forest plot is a way to visualise the effect size and confidence interval for each individual study in a meta-analysis. We can create a forest plot using the function forest:

infarctVolFor <- forest(
  # specify the meta-analysis to plot
  infarctVolMA,  
  
  # sort the data according to effect size
  sortvar = TE,   
  
  # do not plot the fixed effect estimate
  comb.fixed = FALSE,  
  
  # plot the random effects estimate
  comb.random = TRUE,  
  
  # specify the columns to print on the left side of the plot
  leftcols = "studlab",
  
  # specify the labels for the left side plot columns
  leftlabs = "Study",  
  
  # specify the labels for the right side plot columns
  rightlabs = c("NMD", "95% CI", "Weight"),
  
  # specify the label for the heterogeneity statistics
  hetlab = "Het: ",
  
  # specify the x axis label
  xlab = "NMD (% reduction in infarct volume)",  
  
  # specify the graph label on right side of plot
  label.right = "Favours IL-1RA",  
  
  # specify the graph label on left side of plot
  label.left = "Favours Control",  
  
  # specify the x axis limits
  xlim = c(-150,150),  
  
  # specify the size of text (in points)
  fontsize = 10,  
  
  # specify the width of the plotting region
  plotwidth = "10cm",  
  
  # specify minimal number of significant digits for treatment effects
  digits = 1,  
  
  # specify minimal number of significant digits for standard errors
  digits.se = 1  
)

Questions:

1.3 Exploring heterogeneity

Next we will investigate sources of heterogeneity in our data set using meta-regression. A meta-regression model uses one or more explanatory variables to predict differences in the true effect sizes of our outcome variable. These explanatory variables can be continuous or categorical and are sometimes called ‘potential effect modifiers’ or ‘covariates’. We will investigate each variable individually, using univariable models, but multivariable meta-regression models can also be implemented, if there is an adequate number of included studies. For more information on multivariable, or multiple, meta-regression see here.

We will investigate four variables in different subsets of the infarct volume data, as performed in the publication McCann et al. 2016. For this analysis we will use the function metareg from the {meta} package.

Meta-regression: blinding

First, investigate the variable “blinded induction of ischaemia” in the complete data set. This variable addresses whether authors reported blinding of researchers during stroke induction in the animals. Blinding is an important variable when assessing risks of bias.

Categorical variables like this one can be added to a meta-regression model using what are known as “dummy variables”. In this tutorial we will allow R to handle the coding of these dummy variables automatically, but you can also code them manually. For more information see here.

# check the different levels of the variable and number of experiments in each
table(infarctVol$blinded_induction_of_ischaemia)  

infarctVolBlind <- metareg(
  # specify the meta-analysis
  infarctVolMA,  
  
  # specify the experimental variable of interest
  ~ blinded_induction_of_ischaemia   
)  

# view the meta-regression
infarctVolBlind

You can create a bar chart of these results using the code below. Don’t worry about trying to understand the code at this stage, it’s very long!

infarctVolBlind.p <- predict.rma(infarctVolBlind, addx = TRUE)
blind <- as_tibble(infarctVolBlind.p)%>% group_by(X.blinded_induction_of_ischaemiaTRUE)
blind.s <- blind %>% summarise(
  Blind = mean(X.blinded_induction_of_ischaemiaTRUE),
  es = mean(pred),
  se = mean(se),
  ci.lb = mean(ci.lb),
  ci.ub = mean(ci.ub))
blind.s <- blind.s %>% 
  mutate(Blind=replace(Blind,Blind==0, "No"),
  Blind=replace(Blind,Blind==1, "Yes"))
N <- (infarctVol$treatment_group_n + infarctVol$control_group_calc_n)
infarctVol <- cbind(infarctVol, N)
blind.n <- infarctVol %>% 
          group_by(blinded_induction_of_ischaemia) %>%
          summarise(
          sum.N = sum(N))

ggplot(data=blind.s, 
       aes(x=Blind, y=es)) +
  geom_rect(xmin=0.4, 
            xmax=2.6, 
            ymin = infarctVolMA$lower.random, 
            ymax = infarctVolMA$upper.random,
            alpha=0.2) +
  geom_bar(data=blind.s, 
           aes(x=Blind, y=es),
           color = "cadetblue4", 
           fill = "cadetblue", 
           width = c(sqrt(blind.n$sum.N[2])*0.03, sqrt(blind.n$sum.N[1])*0.03),
           stat = "identity") +
  geom_hline(yintercept = 0)  +
  geom_errorbar(data=blind.s, 
                aes(x = Blind, ymin = ci.lb, ymax = ci.ub),
                width=0.1, 
                colour="cadetblue4", 
                alpha=0.9, 
                size=1) +
  scale_y_continuous(name="NMD (% reduction in infarct volume)", 
                     limits=c(0, 50), 
                     breaks = seq(0, 50, 10), expand = c(0,0)) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"),
        axis.title.x = element_text(margin = margin(t = 20, b = 20), size = 16),
        axis.title.y = element_text(margin = margin(r = 20), size = 16),
        axis.text.x= element_text(size = 14),
        axis.text.y= element_text(size = 14))+ 
  labs(x = "Blinded induction of ischaemia", caption = "The horizontal grey bar represents the 95% CIs of the global estimate.
         The width of the bars is proportional to the number of animals in each subgroup.")

Questions:

Meta-regression: route of drug delivery

Next, investigate the variable “route of delivery” in the subset of data where the drug delivery mode is “Protein”.

In this data set, IL-1RA was administered using several different delivery modes (column drug_delivery_mode). We will restrict subsequent analyses to those experiments where IL-1RA was delivered in protein form, due to uncertainty around the timing and effective dose achieved in transgenic and transfection studies (the alternative drug delivery modes).

Route of delivery refers to how IL-1RA was administered e.g. intravenous, subcutaneous. Several different methods were used in these experiments, which may have an influence on effect sizes.

First we will create a new data set including only the protein studies using the function subset:

infarctVol2 <- subset(
  # specify the data to subset
  infarctVol,  
  
  # specify the relevant variable
  !(drug_delivery_mode %in%  
      
      # specify the variable levels to exclude from analysis  
      c("Transgenic", "Vector", "IL-1RA Tg BM cells", "IL-1 RA Tg BM cells"))  
) 

# View the new data subset
infarctVol2 %>% 
  mutate_if(is.numeric, round, digits = 2)

Have a look at the table above and check that the only studies remaining are those where the drug delivery mode is “Protein”.

Next we will run a meta-analysis on this subset of the data:

infarctVolMA2 <- metagen(
  # specify the variable that contains the effect size for each experiment 
  ESi,  
  
  # specify the variable that contains the standard error for each experiment
  SEi,  
  
  # specify the data set
  data = infarctVol2,  
  
  # specify the study labels
  studlab = study_id,  
  
  # specify a random effects model
  comb.random = TRUE,  
  
  # specify that the fixed effect estimate should not be printed
  comb.fixed = FALSE,
  
  # specify which method is used to estimate the between-study variance τ^2, restricted maximum-likelihood (REML) estimator
  method.tau = "REML",  
)

# view the results of the meta-analysis 
infarctVolMA2

Now check the variable data and perform the meta-regression investigating “route of delivery”:

# check the different delivery routes and number of experiments using each
table(infarctVol2$route_of_delivery)
 
infarctVolRoute <- metareg(
  # specify the meta-analysis
  infarctVolMA2,   
  
  # specify the experimental variable of interest
  ~ factor(route_of_delivery)  
)  
     
# view the results of the meta-regression 
infarctVolRoute

We can also graph the results of a meta-regression using a forest plot:

infarctVolRoute.p <- predict.rma(infarctVolRoute, addx = TRUE)

infarctVol2<- cbind(infarctVol2, infarctVolRoute.p$pred) 

route <- as_tibble(infarctVolRoute.p)%>% group_by(pred)

route.s <- route %>% summarise(
  Route = mean(pred),
  es = mean(pred),
  se = mean(se),
  ci.lb = mean(ci.lb),
  ci.ub = mean(ci.ub),
  IP = mean(X.factor.route_of_delivery.IPeritoneal),
  IV = mean(X.factor.route_of_delivery.IVenous),
  Stereo = mean(X.factor.route_of_delivery.Stereotactic),
  SubCut = mean(X.factor.route_of_delivery.SubCut),
  ICV = mean(X.intrcpt)
)

route.s <- route.s %>% 
  mutate(Route=replace(Route,IP==1, "IP"),
  Route=replace(Route,IV==1, "IV"),
  Route=replace(Route,Stereo==1, "Stereo"),
  Route=replace(Route,SubCut==1, "SubCut"),
  Route=replace(Route,IP==0&IV==0&Stereo==0&SubCut==0&ICV==1, "ICV")
)

infarct.for.route <- forest.rma(
  infarctVolRoute,
  xlab = "NMD (% reduction in infarct volume)",
  mlab = "NMD",
  ylim = c(-2.5, 81),
  xlim = c(-300, 300),
  alim = c(-150, 150),
  rows=c(75.5:60.5, 56.5:40.5, 36.5:13.5, 9.5:6.5, 2:-1),
  header = c("Study label", "NMD [95% CI]"),
  order = order(-infarctVol2$`infarctVolRoute.p$pred`, infarctVol2$ESi),
  slab = infarctVol2$study_id,
  addfit = FALSE,
  digits = 1,
  cex = 0.9,
  cex.lab = 1,
  efac = 0.5,
  #top = 1
  )
text(-300, 
     c(77, 58, 38, 11, 3.5), 
     c("Intracerebroventricular", "Intravenous", "Subcutaneous", "Intraperitoneal", "Stereotactic"), 
     adj = c(0,1),
     font = 4,
     cex = 0.9
     )
addpoly(rev(route.s$es), 
        ci.lb = rev(route.s$ci.lb), 
        ci.ub   = rev(route.s$ci.ub),
        rows = c(59, 39, 12, 5, -2.5),
        cex = 0.9,
        col = "cadetblue",
        border = "cadetblue4",
        digits = 1,
        efac = 0.8
)

Questions:

Meta-regression: species

Using the same data set (infarctVol2) and meta-analysis (infarctVolMA2), now investigate whether animal species is a source of heterogeneity:

table(infarctVol2$species)

infarctVolspecies <- metareg(
  
  # specify the meta-analysis
  infarctVolMA2,   
  
  # specify the experimental variable of interest
  ~ factor(species)   
)   

# view the results of the meta-regression
infarctVolspecies

Questions:

Meta-regression: dose of IL-1RA

Lastly, investigate the variable “dose” in the subset of data where the drug delivery mode is “Protein” and the location of administration is “Central”.

Now let’s check if the dose of IL-1RA is associated with effect size differences (i.e. is there a dose response) using only data where the drug was administered centrally (i.e. directly into the brain) and excluding data where the drug was administered peripherally.

First we will create a new data set including only the centrally-administered protein studies, again using the function subset:

infarctVol3 <- subset(
  # specify the data to subset
  infarctVol2,   
  
  # specify the relevant variable
  drug_delivery_location %in%  
    
    # specify the experiments to include in the analysis
    "Central"  
)     

# View the new data subset
infarctVol3 %>% 
  mutate_if(is.numeric, round, digits = 2)

Have a look at the table above and check that the only studies remaining are those where the drug delivery location is “Central”.

Next run a meta-analysis on this new data set:

infarctVolMA3 <-  metagen(
  # specify the variable that contains the effect size for each experiment
  ESi,  
  
  # specify the variable that contains the standard error for each experiment
  SEi,  
  
  # specify the data set
  data = infarctVol3,  
  
  # specify the study labels
  studlab = study_id,  
  
  # specify a random effects model
  comb.random = TRUE,  
  
  # specify that the fixed effect estimate should not be printed
  comb.fixed = FALSE,
  
  # specify which method is used to estimate the between-study variance τ^2, restricted maximum-likelihood (REML) estimator
  method.tau = "REML",  
)

# view the results of the meta-analysis 
infarctVolMA3

Finally, perform the meta-regression investigating “dose”:

infarctVoldose <- metareg(
  # specify the meta-analysis
  infarctVolMA3,   
  
  # specify the experimental variable of interest
  ~ "dose"   
)  
 
# view the results of the meta-regression
infarctVoldose

To graph a meta-regression of a continuous variable we can use a bubble plot:

# specify the meta-regression analysis to plot 
bubble(infarctVoldose, 
       
  # specify the x axis label 
  xlab = "IL-1RA dose (µg)", 
  
  # specify the y axis label
  ylab = "NMD (% reduction in infarct volume)"  
)

Questions:

1.4 Data dictionary: infarct volume

Variable	Description	Levels
study_id	a study identifier for each experimental group within a study	Identifier variable - author, year of publication, experimental group
species	the animal species used in the experiment	Mouse, Rat
strain	the animal strain used in the experiment	Balb/c, C57BL/6, C57BL/6J, C57/SV129, CD-1, Sprague dawley, Unknown, Wistar
sex	the sex of the animals used in the experiment	Male, Unknown
drug	the drug administered in the experiment	IL-1RA
drug_delivery_mode	the mode of the drug delivered to the animals	IL-1 RA Tg BM cells, Protein, Transgenic, Vector
route_of_delivery	the route used to deliver the drug to the animals	ICerebVentricular, IPeritoneal, IVenous, Stereotactic, SubCut
drug_delivery_location	the location of the drug delivery into the animals	Central, Peripheral
dose	the dose of the drug administered to the animals	Numerical - continuous
unit	the unit of the drug dose administered to the animals	cells, mg/kg, ng/g, Tg, µg
time_of_administration	the time of the drug administration, relative to stroke induction, in minutes	Numerical - continuous (mins)
number_of_doses	the number of doses of the drug administered to the animals	Continuous, Multiple, Single
outcome_measure	the outcome measured in the animals	Infarct volume
time_of_assessment	the time the outcome was assessed, relative to stroke induction, in hours	Numerical - continuous (hours)
blinded_assessment_of_infarct_volume	reporting of blinded assessment of outcome	FALSE, TRUE
blinded_induction_of_ischaemia	reporting of blinded induction of ischaemia (stroke)	FALSE, TRUE
random_allocation	reporting of random allocation of animals to treatment or control groups	FALSE, TRUE
sample_size_calculation	reporting of a sample size calculation	FALSE, TRUE
use_of_comorbid_animals	use of animals with a comorbidity	FALSE, TRUE
animal_welfare_compliance	reporting of a statement of compliance with animal welfare regulations	FALSE, TRUE
statement_of_potential_coi	reporting of a statement of potential conflicts of interest	FALSE, TRUE
duration_of_ischaemia	the duration of ischaemic stroke in minutes	Numerical - continuous (mins)
type_of_ischaemia	the type of ischaemic stroke	Permanent, Temporary, Thrombotic
anaesthetic	the anaesthetic used during stroke induction	Chloral Hydrate, Halothane, Isoflurane, Ketamine, Tribromoethanol, Unknown
method_of_mcao	the method used for middle cerebral artery occlusion (mcao) to induce stroke	electrocoagulation, Intraluminal filament/suture, Ligation, thrombin
control_group_mean	the mean value of the outcome measure in the control group	Numerical - continuous
treatment_group_mean	the mean value of the outcome measure in the treatment group	Numerical - continuous
control_group_n	the number of animals in the control group	Numerical - discrete
treatment_group_n	the number of animals in the treatment group	Numerical - discrete
control_group_calc_n	the number of animals in the control group after adjusting for the number of treatment groups it serves	Numerical - continuous
treatment_groups_per_control	the number of treatment groups compared to each control group wihin an experiment	Numerical - discrete
control_group_sd	the standard deviation in the control group	Numerical - continuous
treatment_group_sd	the standard deviation in the treatment group	Numerical - continuous

2. Meta-analysis of binary data

We will now move on to meta-analyse our second outcome of interest, mortality. Mortality is a binary outcome (animals are alive or deceased) so we need to use an effect size appropriate for this data type. There are several such effect sizes but we will use an odds ratio, which represents the odds that our outcome (death) will occur given IL-1RA treatment, compared to the odds of our outcome occurring in the absence of treatment (control).

For this exercise, we will use the function metabin for meta-analysis of binary outcomes. metabin will calculate an odds ratio effect size directly if we provide the number of animals and events (deaths) in the treatment and control groups, so there is no need to manually calculate these values as we did for NMD effect sizes in section 1.

Load and view the data

mortality <- read.csv("IL-1RA_Mortality_section2.csv", 
                      fileEncoding = "UTF-8-BOM")

mortality %>% 
  mutate_if(is.numeric, round, digits = 2)

Take a look at the table above and compare the outcome data to that from the infarct volume (continuous) data in section 1.

You can find a description of each variable in the data dictionary table in section 2.1.

Question:

Meta-analysis

Now perform a meta-analysis using the function metabin:

mortalityMA <- metabin(
  # specify the number of events in experimental group
  event.e = number_affected_in_treatment_group,
  
  # specify the number of observations in the experimental group
  n.e = treatment_group_n, 
  
  # specify the number of events in control group
  event.c = number_affected_in_control_group, 
  
  # specify the number of observations in the control group
  n.c = control_group_n,  
  
  # specify that the meta-analysis should include studies with zero or all events in both groups
  allstudies = TRUE,  
  
  # specify the study labels
  studlab = study_id,   
  
  # specify the data set
  data = mortality,    
  
  # specify which method is used to estimate the between-study variance τ^2, restricted maximum-likelihood (REML) estimator
  method.tau = "REML", 
  
  # specify the summary measure, "OR" = odds ratio
  sm = "OR", 
  
  # specify a random effects model
  comb.random = TRUE,  
 )
 
  # view the results of the meta-analysis
  mortalityMA

Questions:

Forest plot

Next we will create a forest plot of the mortality data using the function forest:

mortalityFor <- forest(
 # specify the meta-analysis to plot
 mortalityMA,  
 
 # sort the data according to effect size
 sortvar = TE,   
 
 # do not plot the fixed effect estimate
 comb.fixed = FALSE,  
 
 # plot the random effects estimate
 comb.random = TRUE,  
 
 # specify the effect size label
 smlab = "Odds ratio",  
 
 # specify the graph label on right side of plot
 label.right = "Favours Control",
 
 # specify the graph label on left side of plot
 label.left = "Favours IL-1RA",  
 
 # specify the size of text (in points)
 fontsize = 11,  
 
 # specify minimal number of significant digits for treatment effects
 digits = 2, 
 
 # specify minimal number of significant digits for standard errors
 digits.se = 2  
)

Questions:

2.1 Data dictionary: mortality

Variable	Description	Levels
study_id	a study identifier for each experimental group within a study	Identifier variable - author, year of publication, experimental group
drug	the drug administered in the experiment	IL-1RA
number_affected_in_control_group	the number of animals affected by the outcome measure in the control group	Numerical - discrete
number_affected_in_treatment_group	the number of animals affected by the outcome measure in the treatment group	Numerical - discrete
control_group_n	the number of animals in the control group	Numerical - discrete
treatment_group_n	the number of animals in the treatment group	Numerical - discrete
outcome_measure	the outcome measured in the animals	Mortality
blinded_assessment_of_outcome	reporting of blinded assessment of outcome	FALSE, TRUE
blinded_induction_of_ischaemia	reporting of blinded induction of ischaemia (stroke)	FALSE, TRUE
method_of_mcao	the method used for middle cerebral artery occlusion (mcao) to induce stroke	electrocoagulation, Intraluminal filament/suture, thrombin
random_allocation	reporting of random allocation of animals to treatment or control groups	FALSE, TRUE
sample_size_calculation	reporting of a sample size calculation	FALSE, TRUE
sex	the sex of the animals used in the experiment	Male
type_of_ischaemia	the type of ischaemic stroke	Permanent, Temporary, Thrombotic
comorbidity	the comorbidity suffered by the animals	Aged, None
dose	the dose of the drug administered to the animals	Numerical - continuous
unit	the unit of the drug dose administered to the animals	mg/kg
time_of_administration	the time of the drug administration, relative to stroke induction, in minutes	Numerical - continuous (mins)
time_of_assessment	the time the outcome was assessed, relative to stroke induction, in hours	Numerical - continuous (hours)

3. Test your skills: meta-analysis

Now you can explore the neurobehavioural outcomes data set by writing your own code. We have only included experiments where IL-1RA was administered in protein form and motor / sensory outcomes were assessed.

In this data set, a normalised mean difference effect size (column ESi) and standard error (column SEi) have already been calculated for each experiment.

You can edit the code we used to analyse the infarct volume data above and compare your results to those in the paper, McCann et al. 2016. If you are unsure what code to use, click on the Hint(s) buttons.

Start by loading and viewing the data set, you can find a description of each variable in the data dictionary table in section 3.1:

# load and view the data set "IL-1RA_Neurobehaviour_nested_section3.csv" using the function read.csv().

# load and view the data set "IL-1RA_Neurobehaviour_nested_section3.csv" using the function read.csv().
neurobeh <- read.csv("IL-1RA_Neurobehaviour_nested_section3.csv", fileEncoding = "UTF-8-BOM")

head(neurobeh)

Then, perform a random effects meta-analysis:

# use the function metagen()

# Perform a random effects meta-analysis using the function metagen().
neurobehMA <- metagen(
  # specify the variable that contains the effect size for each experiment 
  ESi,  
  
  # specify the variable that contains the standard error for each experiment
  SEi,
  
  # specify the data set
  data = neurobeh,  
  
  # specify the study labels
  studlab = study_id,  
  
  # specify a random effects model
  comb.random = TRUE,  
  
  # specify not to print the fixed effect model
  comb.fixed = FALSE,
  
  # specify which method is used to estimate the between-study variance τ^2, restricted maximum-likelihood (REML) estimator
  method.tau = "REML"  
)

 # view the meta-analysis
neurobehMA

Next, visualise the results of the meta-analysis with a forest plot:

# Create a forest plot of the data using the function forest().

# Create a forest plot of the data using the function forest().
neurobehFor <- forest(
  neurobehMA,                                   # specify the meta-analysis to plot
  sortvar = TE,                                 # sort the data according to effect size
  comb.fixed = FALSE,                           # do not plot the fixed effect estimate
  comb.random = TRUE,                           # plot the random effects estimate
  leftlabs = c("Study", "ESi", "SEi"),                           # specify the labels for the left hand plot columns
  xlab = "NMD (% improvement in neurobehavioural outcomes)", # specify the x axis label
  smlab = "NMD",                                # specify the effect size label
  label.right = "Favours IL-1RA",               # specify the graph label on right side of plot
  label.left = "Favours Control",               # specify the graph label on left side of plot
  xlim = c(-150,150),                           # specify the x axis limits
  fontsize = 10,                                # specify the size of text (in points)
  plotwidth = "10cm",                           # specify the width of the plotting region
  digits = 1,                                   # specify minimal number of significant digits for treatment effects
  digits.se = 1                                 # specify minimal number of significant digits for standard errors
)

Question:

Now you can select some variables to investigate as possible sources of heterogeneity using the function metareg e.g. species, sex, dose, time of administration, route of delivery, use of comorbid animals, method of induction of ischaemia, type of ischaemia, anaesthetic used, time of assessment, single vs. multiple administrations.

Remember that the variable name you use in your meta-regression model must match the column name in the data sheet exactly. To see a list of all the column names, you can use the function colnames or check the data dictionary in section 3.1:

# use the metareg() function

# check the names of all the columns in the data set
colnames(neurobeh)

# check the different levels of the variable and number of experiments in each
table(neurobeh$species)

neurobehspecies <- metareg(
  
  # specify the meta-analysis
  neurobehMA,   
  
  # specify the experimental variable of interest
  ~ factor(species)   
)   

# view the results of the meta-regression
neurobehspecies

Questions:

3.1 Data dictionary: neurobehaviour

Variable	Description	Levels
study_id	a study identifier for each experimental group within a study	Identifier variable - author, year of publication, experimental group
species	the animal species used in the experiment	Mouse, Rat
strain	the animal strain used in the experiment	Balb/c, C57BL/6, C57BL/6J, Unknown, Wistar
sex	the sex of the animals used in the experiment	Male, Unknown
drug	the drug administered in the experiment	IL-1RA
drug_delivery_mode	the mode of the drug delivered to the animals	Protein
drug_delivery_location	the location of the drug delivery into the animals	Peripheral
route_of_delivery	the route used to deliver the drug to the animals	IPeritoneal, IVenous, SubCut
dose	the dose of the drug administered to the animals	Numerical - continuous
unit	the unit of the drug dose administered to the animals	mg/kg
number_of_doses	the number of doses of the drug administered to the animals	Multiple, Single
time_of_administration	the time of the drug administration, relative to stroke induction, in minutes	Numerical - continuous (mins)
blinded_induction_of_ischaemia	reporting of blinded induction of ischaemia (stroke)	FALSE, TRUE
blinded_assessment_of_outcome	reporting of blinded assessment of outcome	FALSE, TRUE
random_allocation	reporting of random allocation of animals to treatment or control groups	FALSE, TRUE
sample_size_calculation	reporting of a sample size calculation	FALSE, TRUE
use_of_comorbid_animals	use of animals with a comorbidity	FALSE, TRUE
animal_welfare_compliance	reporting of a statement of compliance with animal welfare regulations	FALSE, TRUE
statement_of_potential_coi	reporting of a statement of potential conflicts of interest	FALSE, TRUE
method_of_mcao	the method used for middle cerebral artery occlusion (mcao) to induce stroke	electrocoagulation, Intraluminal filament/suture, Ligation
type_of_ischaemia	the type of ischaemic stroke	Permanent, Temporary
duration_of_ischaemia	the duration of ischaemic stroke in minutes	Numerical - continuous (mins)
anaesthetic	the anaesthetic used during stroke induction	Halothane, Isoflurane, Ketamine, Tribromoethanol, Unknown
outcome_measure	the outcome measured in the animals	NeurobehaviouralScore
outcome_measure_type	the type of behaviour the outcome assessed	Motor/sensory
time_of_assessment	the time the outcome was assessed, relative to stroke induction, in hours	Numerical - continuous (hours)
ESi	the normalised mean difference effect size	Numerical - continuous
SEi	the normalised mean difference standard error	Numerical - continuous

4. Small study effects

Small study effects in meta-analysis can be assessed using a funnel plot, which is a scatter plot of individual study effect sizes against a measure of each study’s precision. Ideally, the plot should approximately resemble a symmetrical (inverted) funnel. If there is bias due to missing results, for example because smaller studies without statistically significant effects remain unpublished, this will lead to asymmetry in the appearance of the funnel plot.

In this exercise, we will use the complete infarct volume data set and meta-analysis from the first section of this tutorial.

First, draw a funnel plot using the function funnel:

infarctVolFun <- funnel(
 # specify the meta-analysis to plot
 infarctVolMA,  
 
 # plot the random effects estimate
 comb.random = TRUE,  
 
 # specify the x axis label
 xlab = "NMD (% reduction in infarct volume)",  
 
 # specify the measure plotted on the y axis, "invse" = inverse of the standard error
 yaxis = "se"   
)

# infarctVolFun

Examine the funnel plot and answer the following questions:

Funnel plot asymmetry

In addition to assessing the funnel plot visually, we can test for asymmetry statistically. Here we will use Egger’s test, which is a weighted linear regression of the treatment effect on its standard error.

Test for funnel plot asymmetry using the function metabias:

infarctVolBias <- metabias(
 # specify the meta-analysis to analyse
 infarctVolMA,  
 
 # specify the method behind the test statistic
 method.bias = "Egger", 
 
 # plot the results
 plotit = TRUE  
)  
        
# view the results of the funnel plot asymmetry analysis
infarctVolBias

Question

Trim and fill

We can also estimate and adjust for the number and outcomes of theoretically missing studies in a meta-analysis, using a method called trim and fill. We can test for and impute missing studies using the function trimfill:

infarctVolTrim <- trimfill(
  # specify the meta-analysis to plot
  infarctVolMA,
  
  # specify the side of the plot where missing studies are expected
  left = TRUE,  
  
  # plot the random effects estimate
  comb.random = TRUE,  
  
  # specify which method is used to estimate the between-study variance τ^2, restricted maximum-likelihood (REML) estimator
  method.tau = "REML"  
)           

# view the results of the trim-and-fill analysis
infarctVolTrim

You can draw a funnel plot of the filled data using the function funnel:

infarctVolTrimFun <- funnel(infarctVolTrim,
  # plot the random effects estimate
  comb.random = TRUE,  
  
  # specify the x axis label
  xlab = "NMD (% reduction in infarct volume)", 
  
  # specify the measure plotted on the y axis, "invse" = inverse of the standard error
  yaxis = "se"  
)

Questions:

5. Test your skills: small study effects

Now you can assess and visaulise small study effects in the neurobehavioural data set. Create a funnel plot, test for asymmetry, and impute missing studies using the funnel, metabias, and trimfill functions. In section 3. Test your skills: meta-analysis we used only those experiments where IL-1RA was administered in protein form and motor / sensory outcomes were assessed. To assess for small study effects, we will use the entire neurobehavioural data set.

We have provided the first part of the code to load the data and perform the meta-analysis. Start by writing code to create a funnel plot visualising these data:

# first, load the data set "IL-1RA_Neurobehaviour_prenest_section5.csv" using the function read.csv().
neurobeh2 <- read.csv("IL-1RA_Neurobehaviour_prenest_section5.csv", fileEncoding = "UTF-8-BOM")

# perform a random effects meta-analysis using the function metagen().
neurobehMA2 <- metagen(
  ESi,  
  SEi,  
  data = neurobeh2,  
  comb.random = TRUE,
  comb.fixed = FALSE,
  method.tau = "REML"  
)
# view the results of the meta-analysis
neurobehMA2





# use the funnel() function to generate a funnel plot and visually assess the plot for asymmetry

neurobehFun <- funnel(
 # specify the meta-analysis to plot
 neurobehMA2,  
 
 # plot the random effects estimate
 comb.random = TRUE,  
 
 # specify the x axis label
 xlab = "NMD (% improvement in neurobehavioural outcome)",  
 
 # specify the measure plotted on the y axis, "se" = inverse of the standard error
 yaxis = "se"   
)

Next, use Egger’s test to assess funnel plot symmetry statistically:

# use the metabias() function to test for funnel plot asymmetry

neurobehBias <- metabias(
 # specify the meta-analysis to analyse
 neurobehMA2,  
 
 # specify the method behind the test statistic
 method.bias = "Egger", 
 
 # plot the results
 plotit = TRUE  
)  
        
# view the results of the funnel plot asymmetry analysis
neurobehBias

Now use trim and fill to assess whether there are any theoretically missing studies in the data set and then create a funnel plot of the results:

# use the trimfill() function to test for "missing" studies


















# use the funnel() function to plot the filled data set

neurobehTrim <- trimfill(
  # specify the meta-analysis to assess
  neurobehMA2,
  
  # specify the side of the plot where missing studies are expected
  left = TRUE,  
  
  # plot the random effects estimate
  comb.random = TRUE,  
  
  # specify which method is used to estimate the between-study variance τ^2, restricted maximum-likelihood (REML) estimator
  method.tau = "REML"  
)         

# view the results of the trim-and-fill analysis
neurobehTrim

neurobehTrimFun <- funnel(neurobehTrim,
  # plot the random effects estimate
  comb.random = TRUE,  
  
  # specify the x axis label
  xlab = "NMD (% improvement in neurobehavioural outcome)", 
  
  # specify the measure plotted on the y axis, "se" = inverse of the standard error
  yaxis = "se"  
)

Question:

Notes and resources

We have created this tutorial to guide users through a preclinical meta-analysis.

Meta-analysis comes with many researcher degrees of freedom and it’s important that when performing your own meta-analysis, your statistical analysis plan is specified a priori in your study protocol.

In this tutorial we have implemented standard analyses, as performed in a published meta-analysis, but in the functions we have used (and the many others available), multiple details can be tailored to suit different analyses.

If you use R to analyse your data (or any other statistical software package!), you should have a thorough understanding of how the analysis is run and why, and the assumptions and default settings that are built into the functions you use in different software applications.

There are many more things to explore in R, with lots of options to customise graphs and data visualisations for presentations and publications.

For more information on any of the R packages and functions we have used, you can use the help function in R studio, which is described in detail on the Getting Help with R page, along with multiple options on where else to look for information and help.

A useful online resource that describes meta-analysis with R in detail is “Doing Meta-Analysis with R: A Hands-On Guide” .

If you have specific problems developing code in R, there is also a thriving online community with whom to engage. Many platforms offer forums to submit problems and solutions, one of the most popular is Stack Exchange.

If you have questions about your own systematic review of animal data, you can contact CAMARADES Berlin. Alternatively, you can contact the CAMARADES SyRF help desk.

This resource has been created using CC-BY-4.0 sources including: Viechtbauer W. Conducting meta-analyses in R with the metafor package. Journal of Statistical Software 36(3), 1-48 (2010).

This resource was last updated on: 20 June, 2026

About CAMARADES

The materials presented here have been developed by CAMARADES Berlin. CAMARADES (Collaborative Approach to Meta-Analysis and Review of Animal Data from Experimental Studies) provide a framework supporting researchers to perform systematic reviews of data from experimental animal studies. This support includes best practice methodological assistance and access to collaborative tools such as our online review facility, SyRF.

If you have questions about this resource, or would like to seek support for your systematic review, please email us.

CAMARADES Berlin are supported by Charité 3R. For more information about the 3Rs at Charité – Universitätsmedizin Berlin, visit the Charité 3R Toolbox.

CAMARADES Berlin are located at the QUEST Center for Responsible Research, Berlin Institute of Health (BIH) at Charité