Learning Objectives

Use a published dataset to explore some common meta-analysis techniques.

Install the dmetar R package directly from the developer’s github repository.

Perform fixed-effect and random-effects meta-analysis using the meta and metafor packages.

Explain heterogeneity among studies using subgroup analyses.

Plot results of your meta-analysis using forest plots.

Determine whether publication bias is present in your data.

Why Meta-Analysis?

Throughout this course, we have run into the issue of generality of experimental results. Researchers may find a large effect in a particular study, but how are we to know that this result was not a matter of chance, or that this result would be the same if repeated in a different system? Enter systematic review and meta-analysis, which provide a scientific method for performing reviews of the literature.
Systematic review involves reporting the methods for a review the same as we would for any other study. When performing a systematic review, we must report the databases used to search for studies (e.g., Web of Science, Google Scholar, etc.), a list of specific search terms used (e.g. “plant” AND “denitrification”, OR “veg” AND “denit*“, etc.), a list of all studies returned by these searches, and if/why studies were eliminated from further consideration.
Meta-analysis refers to a set of statistical techniques that allows us to synthesize the results of several independent studies to test hypotheses using the body of results available in the literature. The strength of meta-analysis is that it allows us to combine the results of independent studies, weighted by the sample size and precision of those studies, to estimate mean effects and to estimate uncertainty in those effects. Meta-analysis can be used to test our preconceived notions of what we have learned to date, uncover emerging trends, or learn what knowledge is still missing in a field.
Today we will learn to use three meta-analysis packages that have been developed for R.

Modified from xkcd.com/2295

Getting started

If you haven’t already, be sure to open your Environmental Data Analysis project and save a new script for this workshop.

Installing Packages

The meta and metafor packages are installed on the computers in Hudson 032. If using your personal computer, you can install them using the install.packages() function.

The package dmetar is not available in the cran repository, so we cannot install it using the install.packages() function or the Packages window. This is common for packages developed by novice programmers, as compiling a package for cran is difficult, and keeping track of updates is time consuming. However, you can download the package directly from the developer’s github repository, which is a useful thing to know how to do!

First you will need to use install.packages() to install the devtools package. Then you can use install_github() to install the dmetar package from the developer’s repository.

install.packages("devtools")
devtools::install_github("MathiasHarrer/dmetar")

This code will initiate the installation process, which may take a bit of time.
During the process, the installation manager may ask if you want to update existing R packages. The safest answer is to enter 3 for None and then hit Enter.

If the installation manager asks if you wish to “install from sources the package which needs compilation?”, answer n for no.

Loading packages

When the installation completes successfully, use the library() function to load all of the packages we will use. (Note: If the dmetar package fails to load, you will still be able to complete the assignment.)

library(meta)
library(metafor)
library(dmetar)

Formatting Data for Meta-Analysis

To conduct a meta-analysis that will allow you to statistically synthesize and analyze the results of several studies, you must be able to extract each of the following values from each study:

A mean for both the control and treatment group.
A measure of variation for both the control and treatment means (standard deviation, standard error, or variation).
The number of replicate samples used to calculate the mean.

You may need to extract these values from graphs using image-analysis software such as ImageJ. Though reporting has improved over the years, studies will often be unusable due to the failure of the authors to clearly report one or more of these values. The moral of the story is, if you want your experiment cited in a meta-analysis (and you do), be sure to clearly report these values!

All of the best practices for data spreadsheets apply to today’s assignment. Additionally, the following variables (or variable types) must be included in a dataset that is to be used for meta-analysis:

Column	Description
Author	This signifies the column for the study label (i.e., the first author)
Me	The mean of the experimental/intervention group
Se	The standard deviation of the experimental/intervention group
Mc	The Mean of the control group
Sc	The Standard Deviation of the control group
Ne	The number of participants in the experimental/intervention group
Nc	The number of participants in the control group
Subgroup	This is the label for one of your subgroup codes. It is not that important how you name this column, so you can give it a more informative name (e.g.<a0>population). In this column, each study should then be given an subgroup code, which should be exactly the same for each subgroup, including upper/lowercase letters. Of course, you can also include more than one subgroup column with different subgroup codings, but the column name has to be unique

Importing Data for Today’s Exercise

Here we will look at data that I collected for a paper that was eventually published in Ecological Applications (Alldred and Baines 2016). The data include measurements of denitrification rates, or rates of microbial nitrogen removal, that occurred in the rhizospheres of several different plant species. I have simplified the data for our activity today, but a link to the full dataset is available in References.

You can download the data file for today “Alldred2016_MetaAnalysisData.csv” from the course Brightspace page. Save the data file to your project folder, and then import the data using read.csv().

# Import data
madata <- read.csv("Alldred2016_MetaAnalysisData.csv")
#Use str() to view structure of data frame
str(madata)

Data Dictionary for Today’s Dataset

Variable	Definition
Author	Reference for study
Species	Plant species growing in soil/sediments where denitrification was measured
Genus	Genus of the plant species
Family	Family of the plant species
Monocot_Dicot	Whether the plant species was a monocot or a dicot
Growth_form	The growth form of the plant species (i.e., emergent, tree-shrub, forb-emmergent, or submersed)
Wetland_System	Type of wetland system in which the study was performed (i.e., salt marsh, riparian, constructed, microcosm experiments, or mesocosm experiments)
Denitrification_Method	Method used to quantify denitrification rates (i.e., N2:Ar ratios, direct N2 flux measurements using chambers, DEA denitrification enzyme activity, Acetylene inhibition, 15N-NH4 additions, or NO3 removal rates)
Ne	Number of observations for vegetated sediments
Me	Mean denitrification rate for vegetated sediments in g m^(-2) h^(-1)
Se	Standard deviation of measurements in vegetated sediments
Nc	Number of observations for non-vegetated sediments
Mc	Mean denitrification rate for non-vegetated sediments in g m^(-2) h^(-1)
Sc	Standard deviation of measurements in non-vegetated sediments
lnR	Log-response ratio effect size for denitrification rates in vegetated versus non-vegetated sediments
var_lnR	Estimate of variation for log-response ratio effect size for denitrification rates in vegetated versus non-vegetated sediments

Pooling Effect Sizes: Fixed v. Random Effects

In a meta-analysis model, we calculate an effect size to determine the overall effect of a treatment. Effect sizes may be estimated as fixed effects or random effects. Which one is more appropriate is the source of endless debate. The following guidance is adapted from Gurevitch and Hedges (2001) and is generally accepted in the field of ecology:

Typically, you would use a fixed effects model when you are able to assume that your data contain similar studies with similar characteristics that should share a common true effect size. For example, if predators should be expected to decrease prey populations by 20 individuals or 20%, and that effect should be measurable regardless of the method used to measure the prey population, then a fixed effect would be appropriate to estimate the effect of predation.

A random effects model is more appropriate when you expect that studies contain a common true effect but that other differences among the studies also contribute to random variation in the effect size (in addition to sampling variation). Returning to the earlier example, if the studies all measured the effect of predators on prey populations, but some populations were exposed to only one species of predators while others were exposed to several predator species, it may be more appropriate to use a random effects model.

Fig. 1 Example of fixed versus random effects in meta-analysis (Nikolakopoulou et al. 2014).

Choosing an Appropriate Effect Size

The meta package can calculate three types of effect sizes:

Hedge’s d: the mean difference (MD) between the treatment and control group.

Hedge’s g: the standardized mean difference (SMD) between the treatment and control group. (Works better for studies with small sample sizes)

Log-response ratio: the log-transformed ratio of means (ROM) between the treatment group and the control group. (Works well when a proportional response is expected)

Here we will use the log-response ratio ROM because I expect to observe a proportional change in denitrification rates when plants are present, i.e., a X% increase, as opposed to a X difference in g m^-2 h^-1. If you do expect a fixed difference, Hedge’s d or Hedge’s g may be more appropriate as a measure of effect size. Also note that log-response ratios do not work for data that are less than or equal to 0. You would need to remove any non-positive data rows using a filter() function prior to analysis to avoid an error.

Fixed Effects Model

To calculate an overall effect, we average the effect size of independent studies, but we give a greater weight to studies whose estimates have greater precision and/or include a greater number of replicates. To do this, we calculate an inverse of the sampling variance for each study and apply this as a weight when calculating the effect sizes.

Let’s use our data to calculate a fixed effect size for the effect of plants on denitrification rates using a fixed-effect model. Because our data include the raw means, standard deviations, and sample sizes, we will use the metacont() function in the meta package. The format of the function is:

Parameter	Function
Ne	The number of participants (N) in the intervention group
Me	The mean (M) of the intervention group
Se	The standard deviation (SD) of the intervention group
Nc	The number of participants (N) in the control group
Mc	The mean (M) of the control group
Sc	The standard deviation (SD) of the control group
data=	After <91>=<92>, paste the name of your dataset here
studlab=paste()	This tells the function were the labels for each study are stored. If you named the spreadsheet columns as advised, this should be studlab=paste(Author)
comb.fixed=	Whether to use a fixed-effects-model
comb.random	Whether to use a random-effects-model
prediction=	Whether to print a prediction interval for the effect of future studies based on present evidence
sm=	The summary measure we want to calculate. We can either calculate the mean difference (MD), Hedges<92> g (SMD), or a log-response ratio (ROM)
backtransf=	For an “ROM” summary measure, whether to backtransform results in printouts and plots
method.tau=	For a random-effects model, the method used to measure between-study variance

Therefore, we specify our fixed-effect model as follows:

#Fixed effect model, using log ratio as a measure of effect size
fe_mod <- metacont(Ne,
                   Me,
                   Se,
                   Nc,
                   Mc,
                   Sc,
                   data = madata,
                   studlab = paste(Author),
                   comb.fixed = T,
                   comb.random = F,
                   prediction = T,
                   sm = "ROM",
                   backtransf = T)
#View results of fixed-effect model
fe_mod_results <- summary(fe_mod)
fe_mod_results

The output of our model will include:

The individual effect sizes for each study, a confidence interval around the effect size, and the weight of the study in the overall calculation of effect size.
The total number of included studies (k).
The overall effect (in our case, ROM = 1.3875) and its confidence interval and p-value.
Measures of between-study heterogeneity, such as tau² or I², and a Q-test of heterogeneity.

You can save the results of the meta-analysis to your working directory as a .txt file using this command:

sink("fe_mod_results.txt")
print(fe_mod_results)
sink()

Challenge 1

Refer to the results of our fixed-effects model to answer the following questions:
1. Was the overall effect of plants on denitrification significant at alpha = 0.05?
2. Was the mean effect of plants on denitrification positive or negative?
3. How many studies were included in this analysis?
4. Which study was given the greatest weight in the standardized mean effect size?
5. Was there significant heterogeneity in effect among studies?

Random Effects Model

The fixed-effect model assumes that your studies all contain the same true effect, or that the studies all come from the same population of studies. In practice, you may not be able to make that assumption. A random-effects model assumes that there is a distribution of true effect sizes, and that all studies contain sampling error as well as additional sources of random variation.
Specifying a random-effects model in R is very similar to the code we used to specify a fixed-effect model. The main difference is that we will need to specify a method for estimating the variance (or distribution) of true effect sizes. Different disciplines have different conventions, but in ecology we often use maximum-likelihood estimators, which we will specify with the argument method.tau = "ML". You should always consult the literature to determine which estimator will work best for your data.

re_mod <- metacont(Ne,
                   Me,
                   Se,
                   Nc,
                   Mc,
                   Sc,
                   data = madata,
                   studlab = paste(Author),
                   comb.fixed = F,
                   comb.random = T,
                   method.tau = "ML",
                   prediction = TRUE,
                   sm = "ROM")
re_mod_results <- summary(re_mod)
re_mod_results

Save the results of the meta-analysis to your working directory as a .txt file using this command:

sink("re_mod_results.txt")
print(re_mod_results)
sink()

Challenge 2

Refer to the results of our random-effects model to answer the following questions:
1. How does the p-value associated with the overall effect of plants compare to the p-value in the fixed-effects model?
2. How does the effect size estimated in this model compare to the effect size estimated in the fixed-effects model?
3. How does the size of the confidence interval for the effect size estimated in this model compare to the confidence interval estimated in the fixed-effects model?
4. How do the weights given to individual studies differ from the fixed-effects model?
5. Was there significant heterogeneity in effect among studies?

Subgroup Analysis

In both of our models, we found significant heterogeneity in effect among studies. We can attempt to explain some of this heterogeneity in effect using subgroup analyses. For example, when I collected these data, I also recorded information about the studies including the type of system in which the study was conducted, the plant species examined in the study, and higher-level taxonomic information about the plant species. I can group studies by these variables to explain heterogeneity in the effect of plants on denitrification rates.

Fixed v. Random Effects in Subgroup Analysis

If the subgroups are pre-chosen and represent specific effects of interest, then a fixed effect model is appropriate. Typically, the only time you would use a fixed-effect model is for limited, mutually exclusive groups (female v. male OR marine v. freshwater v. terrestrial) or fixed treatment effects (fertilized v. unfertilized OR predator v. no predator).

If the subgroups are randomly chosen or randomly available categories, then a random effect model is more appropriate. One example could include populations from different regions (North America v. Europe). In this case, samples could have been collected or categorized from multiple regions, but we chose only two, or had only two available. This random sample introduces new sampling error, which we must control for using a random-effects model.

We will first use the subgroup.analysis.mixed.effects() function in the dmetar package to conduct subgroup analyses. This model assumes a fixed effect for variance among subgroups, and a random effect for variance within subgroups.

If you were unable to install the dmetar package, you can copy and paste the subgroup.analysis.mixed.effects() function source code below into your console and run the code to create the function. If you were able to install dmetar, you can skip this step.

View source code for subgroup.analysis.mixed.effects() function

subgroup.analysis.mixed.effects = sgame = function(x, subgroups, exclude = "none") {

    # Define variables
    m = x
    subgroups = subgroups
    exclude = exclude

    # Levels of subgroup
    subgroups = as.factor(subgroups)
    k = as.vector(summary(subgroups))
    levels = levels(subgroups)
    k.level.df = data.frame(level = levels, k = k)

    # Out Loop for wrong input
    if (length(subgroups) != length(m$studlab)) {
        stop("Subgroup variable does not contain the same number of cases as the 'meta' object. You need to define a variable which provides a subgroup value for each effect size included in your 'meta' results object.")
    }

    # get 'Exclude' Subgroup level names
    if (exclude[1] != "none") {
        levels = levels[!levels %in% exclude]
        k = k.level.df[(k.level.df$level %in% levels), ]$k
    }

    # Create Loop for subgroups
    list = list()
    for (x in levels) {
        list[[x]] = which(subgroups %in% c(paste(x)))
    }

    # Loop over list to generate subgroup results
    sg.results = list()
    for (x in 1:length(list)) {
        sg.results[[x]] = update.meta(m, subset = list[[x]])
    }

    # Loop over sg.results to get effect size estimates
    ES = vector()
    SE = vector()
    Qsg = vector()
    I2sg = vector()
    I2sg.lower = vector()
    I2sg.upper = vector()
    for (x in 1:length(sg.results)) {
        ES[x] = sg.results[[x]]$TE.random
        SE[x] = sg.results[[x]]$seTE.random
        Qsg[x] = sg.results[[x]]$Q
        I2sg[x] = sg.results[[x]]$I2
        I2sg.lower[x] = sg.results[[x]]$lower.I2
        I2sg.upper[x] = sg.results[[x]]$upper.I2
    }

    me.data = data.frame(Subgroup = levels, TE = ES, seTE = SE)

    # Fixed Meta-Analysis betweens subgroups
    meta = metagen(TE, seTE, data = me.data,
                   comb.fixed = TRUE, comb.random = FALSE,
                   byvar = Subgroup, hakn = FALSE)

    # Create full output dataset

    if (m$sm %in% c("RR", "RD", "OR", "ASD", "HR", "IRR")){

      subgroup.results = data.frame(Subgroup = me.data$Subgroup,
                                    k = k,
                                    TE = me.data$TE,
                                    seTE = me.data$seTE,
                                    sm = exp(me.data$TE),
                                    LLCI = round(exp(meta$lower), 3),
                                    ULCI = round(exp(meta$upper), 3),
                                    p = meta$pval,
                                    Q = Qsg,
                                    I2 = round(I2sg, 2),
                                    I2.lower = round(I2sg.lower, 2),
                                    I2.upper = round(I2sg.upper, 2))
      colnames(subgroup.results)[5] = m$sm

    } else {

      subgroup.results = data.frame(Subgroup = me.data$Subgroup,
                                    k = k,
                                    TE = me.data$TE,
                                    seTE = me.data$seTE,
                                    LLCI = round(meta$lower, 3),
                                    ULCI = round(meta$upper, 3),
                                    p = meta$pval,
                                    Q = Qsg,
                                    I2 = round(I2sg, 2),
                                    I2.lower = round(I2sg.lower, 2),
                                    I2.upper = round(I2sg.upper, 2))
      if (m$sm != ""){
        colnames(subgroup.results)[3] = m$sm
        colnames(subgroup.results)[4] = "SE"
      }

    }

    mixedeffects.results = data.frame(Q = meta$Q,
                                      df = meta$df.Q,
                                      p = meta$pval.Q,
                                      row.names = "Between groups")

    # Create Forest plot data
    forest.m = update.meta(m, byvar=subgroups, comb.fixed = FALSE)

    if (exclude[1] != "none"){

      exclude.level = levels(subgroups)[levels(subgroups) %in% exclude]
      exclude.nums = which(subgroups %in% exclude.level)
      not.exclude.nums = 1:m$k
      not.exclude.nums = not.exclude.nums[!(not.exclude.nums %in% exclude.nums)]
      forest.m = update.meta(m, byvar=subgroups, comb.fixed = FALSE, subset = not.exclude.nums)

    }

    forest.m$TE.random = meta$TE.fixed
    forest.m$lower.random = meta$lower.fixed
    forest.m$upper.random = meta$upper.fixed
    forest.m$Q = meta$Q
    forest.m$df.Q = meta$df.Q
    forest.m$pval.Q = meta$pval.Q

    rownames(subgroup.results) = subgroup.results$Subgroup
    subgroup.results$Subgroup = NULL

    reslist  = list(within.subgroup.results = subgroup.results,
               subgroup.analysis.results = mixedeffects.results,
               m.random = forest.m,
               method.tau = m$method.tau,
               k = sum(k))

    class(reslist) = "subgroup.analysis.mixed.effects"

    invisible(reslist)

    reslist

}

For the subgroup.analysis.mixed.effects() function, the following parameters have to be set:

Code	Description
x	An object of class meta, generated by the metabin, metagen, metacont, metacor, metainc, or metaprop function.
subgroups	A character vector of the same length as the number of studies within the meta-analysis, with a unique code for the subgroup each study belongs to. Must have the same order as the studies in the meta object.
exclude	Single string or concatenated array of strings. The name(s) of the subgroup levels to be excluded from the subgroup analysis. If <91>none<92> (default), all subgroup levels are used for the analysis.
plot	Logical. Should a forest plot for the mixed-effect subgroup analysis be generated? Calls the forest.meta function internally. TRUE by default.

Let’s try a subgroup analysis in which we group studies by whether the plant species was a monocot or dicot.

s_modi <- subgroup.analysis.mixed.effects(x = re_mod,
                                subgroups = madata$Monocot_Dicot)
s_modi

Save the results of the meta-analysis to your working directory as a .txt file using this command:

sink("s_modi_results.txt")
print(s_modi)
sink()

Challenge 3

Refer to the results of our mixed-effects subgroup analysis to answer the following questions:
1. Did the effect of plants on denitrification differ significantly between studies conducted on dicots and studies conducted on monocot plants at alpha = 0.05?
2. Did the effect of plants in either group differ significantly from zero at alpha = 0.05?

Monocot v. dicot is an example of a mutually exclusive category for which a fixed-effect model is appropriate. However, we may also want to examine groupings that were randomly chosen or randomly available in the literature, such as the type of wetland system where the study was performed. We can also run a random-effects subgroup analysis using the update.meta() function in the meta package.

For this function, we have to set two parameters:

Code	Description
byvar	Here, we specify the variable in which the subgroup of each study is stored
comb.random	Weather we want to use a random-effects-model for between-subgroup-differences. In this case, we have to set comb.random = TRUE

Let’s try a subgroup analysis to determine whether the effect of plants differs significantly among wetland systems.

s_wsys <- update.meta(re_mod,
                       byvar = Wetland_System,
                       comb.random = T,
                       comb.fixed = F)
s_wsys

Save the results of the meta-analysis to your working directory as a .txt file using this command:

sink("s_wsys_results.txt")
print(s_wsys)
sink()

Challenge 4

Refer to the results of our random-effects subgroup analysis to answer the following questions:
1. Did the effect of plants on denitrification differ significantly among wetland systems at alpha = 0.05?
2. Did the effect of plants in any wetland system differ significantly from zero at alpha = 0.05?

Forest Plots

Now that we created our meta-analysis models, we can visualize the data in a more digestible way using forest plots, which are conventional tools for visualizing results in meta-analysis publications.

Let’s begin by calling the default function forest() to visualize the results of our random-effects model.

forest(re_mod)

Notice that we have so many studies that it is difficult to visualize them all in one plot! We also don’t necessarily care about the results of individual studies as much as we care about heterogeneity explained in our subgroup analyses. We can add the argument study.results = F to the forest() function to exclude the results of individual studies and show only the differences among subgroups.
The following code will generate a forest plot for our monocot v dicot mixed-effects subgroup analysis:

meta::forest(s_modi,
             study.results = F)

The following code will generate a forest plot for our wetland system mixed-effects subgroup analysis:

meta::forest(s_wsys,
             study.results = F)

Try calling the help file for the forest function ?meta::forest() to see what other aspects of the forest plot you can change.

Publication Bias

A common problem in research is that studies with low or non-significant effect sizes are less likely to be selected for publication. A related problem is that studies with large and significant effect sizes are more likely to be cited, so authors will tend to put more effort into the considerable ordeal of publishing these studies, and let lower-impact results languish in the file drawer (or computer folder). This issue is often referred to as publication bias or file-drawer problem.

Publication bias presents a considerable challenge for meta-analysis. It may lead us to overestimate effect sizes relative to true effect sizes, due to the many missing low-effect or non-significant study results that were never published.

Small Sample Bias: Funnel Plots

Small sample bias methods focus on whether small studies of small effect are underrepresented in the meta-analysis dataset. The reasoning behind these methods is:

Large studies are likely to be published whether or not the results are significant.
With larger sample sizes, results are more likely to be significant, whether or not effect sizes are large.
Small studies are at the greatest risk of being non-significant, so significant small studies are likely to be overrepresented relative to non-significant small studies.

The best way to visualize whether small studies with small effect sizes are missing is to construct a funnel plot.
We can generate a funnel plot for our re_mod meta-analysis output using the funnel() function in meta.

funnel(re_mod, xlab = "Ratio of Difference")

The funnel plot consists of a funnel and two axes, with the effect size on the x axis, and the standard error on the y-axis. We would expect large studies with small standard errors to fall out toward the top of this plot, while small studies with large standard errors would fall out toward the bottom. The ideal would be a symmetrical scatter of studies about the mean effect (vertical dotted line).
Small studies of small effect would show up in the bottom left corner. If studies are missing in the bottom left corner, relative to the bottom right, this would indicate publication bias. This is likely the case for our data, which includes one small study of large effect.

I can add labels to the study points by adding the argument studlab = T. The labels will allow me to determine which study is having a large effect relative to its size (i.e., the outlier in the bottom right corner).

funnel(re_mod, xlab = "Ratio of Difference", studlab = T)

The label indicates that the Ottosen et al. 1999 study has a large effect for its size. I can revisit my model to see how much weight this study was given in my analysis. Small-study bias is a bigger problem for random-effects models, in which they are given greater relative weight than fixed-effects models.

P-Curve

P-curve analysis addresses a different source of publication bias. Rather than assume that researchers fail to publish non-significant results, this set of methods tests for whether researchers “play” with their data (e.g., selectively remove outliers, choose different outcomes, control for different variables, inflate sample size, etc.) until results become significant. This practice is called p-hacking and is very common in certain fields of science.

We will use pcurve() function in the dmetar package. If you were unable to install the dmetar package, you can copy and paste the following source code into the console and run it. The function requires the stringr and poibin packages to work correctly. Don’t panic if you can’t get this to work. I can send you the p-curve output if need be.

View source code for pcurve() function

pcurve = function(x, effect.estimation = FALSE, N, dmin = 0, dmax = 1){

  # Rename x to metaobject, remove x
  metaobject = x
  rm(x)

  # Stop if metaobject is not meta or does not contain TE or seTE column
  if (!(class(metaobject)[1] %in% c("metagen", "metabin", "metacont", "metacor", "metainc", "meta", "metaprop"))){

    for (i in 1:length(colnames(metaobject))){

      te.exists = FALSE

      if (colnames(metaobject)[i]=="TE"){
        te.exists = TRUE
        break
      } else {

      }

    }

    for (i in 1:length(colnames(metaobject))){

      sete.exists = FALSE

      if (colnames(metaobject)[i]=="seTE"){
        sete.exists = TRUE
        break
      } else {

      }

    }

    for (i in 1:length(colnames(metaobject))){

      studlab.exists = FALSE

      if (colnames(metaobject)[i]=="studlab"){
        studlab.exists = TRUE
        break
      } else {

      }

    }

    if(te.exists == FALSE | sete.exists ==FALSE | studlab.exists ==FALSE){
      stop("x must be a meta-analysis object generated by meta functions or a data.frame with columns labeled studlab, TE, and seTE.")
    }

  }


  #Disable scientific notation
  options(scipen=999)

  # Calculate Z
  zvalues.input = abs(metaobject$TE/metaobject$seTE)


  ##############################################
  # 1. Functions ###############################
  ##############################################

  getncp.f =function(df1,df2, power)   {
    error = function(ncp_est, power, x, df1,df2) pf(x, df1 = df1, df2=df2, ncp = ncp_est) - (1-power)
    xc=qf(p=.95, df1=df1,df2=df2)
    return(uniroot(error, c(0, 1000), x = xc, df1 = df1,df2=df2, power=power)$root)  }

  getncp.c =function(df, power)   {
    xc=qchisq(p=.95, df=df)
    error = function(ncp_est, power, x, df)      pchisq(x, df = df, ncp = ncp_est) - (1-power)
    return(uniroot(error, c(0, 1000), x = xc, df = df, power=power)$root)   }

  getncp=function(family,df1,df2,power) {
    if (family=="f") ncp=getncp.f(df1=df1,df2=df2,power=power)
    if (family=="c") ncp=getncp.c(df=df1,power=power)
    return(ncp)  }

  percent <- function(x, digits = 0, format = "f", ...)   {
    paste(formatC(100 * x, format = format, digits = digits, ...), "%", sep = "")
  }

  pbound=function(p) pmin(pmax(p,2.2e-16),1-2.2e-16)

  prop33=function(pc)
  {

    prop=ifelse(family=="f" & p<.05,1-pf(qf(1-pc,df1=df1, df2=df2),df1=df1, df2=df2, ncp=ncp33),NA)
    prop=ifelse(family=="c" & p<.05,1-pchisq(qchisq(1-pc,df=df1),  df=df1, ncp=ncp33),prop)
    prop
  }

  stouffer=function(pp) sum(qnorm(pp),na.rm=TRUE)/sqrt(sum(!is.na(pp)))



  ###############################################################################
  # 2.  Process data ############################################################
  ###############################################################################

  # Note: due to reliance on the pcurve-app function, z-scores are pasted into characters first
  # and then screened to generate variables necessary for further computation

  zvalues.input =  paste("z=", zvalues.input, sep="")
  filek = "input"

  raw = zvalues.input
  raw=tolower(raw)
  ktot=length(raw)

  k=seq(from=1,to=length(raw))

  stat=substring(raw,1,1)
  test=ifelse(stat=="r","t",stat)

  # Create family
  family=test
  family=ifelse(test=="t","f",family)
  family=ifelse(test=="z","c",family)
  #family: f,c        converting t-->f and z-->c

  # Find comma,parentheses,equal sign
  par1 =str_locate(raw,"\\(")[,1]
  par2 =str_locate(raw,"\\)")[,1]
  comma=str_locate(raw,",")[,1]
  eq   =str_locate(raw,"=")[,1]

  # DF for t-tests
  df=as.numeric(ifelse(test=="t",substring(raw,par1+1,par2 -1),NA))

  # DF1
  df1=as.numeric(ifelse(test=="f",substring(raw,par1+1,comma-1),NA))
  df1=as.numeric(ifelse(test=="z",1,df1))
  df1=as.numeric(ifelse(test=="t",1,df1))
  df1=as.numeric(ifelse(test=="c",substring(raw,par1+1,par2 -1),df1))

  # DF2
  df2=as.numeric(ifelse(test=="f",substring(raw,comma+1,par2-1),NA))
  df2=as.numeric(ifelse(test=="t",df,df2))

  equal=abs(as.numeric(substring(raw,eq+1)))

  value=ifelse((stat=="f" | stat=="c"),equal,NA)
  value=ifelse(stat=="r", (equal/(sqrt((1-equal**2)/df2)))**2,value)
  value=ifelse(stat=="t", equal**2 ,value)
  value=ifelse(stat=="z", equal**2 ,value)

  p=ifelse(family=="f",1-pf(value,df1=df1,df2=df2),NA)
  p=ifelse(family=="c",1-pchisq(value,df=df1),p)
  p=pbound(p)  #Bound it to level of precision, see function 3 above

  ksig= sum(p<.05,na.rm=TRUE)     #significant studies
  khalf=sum(p<.025,na.rm=TRUE)    #half p-curve studies

  if (ksig <= 2){
    stop("Two or less effect sizes were detected, so p-curve analysis cannot be conducted.")
  }



  ##############################################################################
  # 3. PP-values ###############################################################
  ##############################################################################

  # Right Skew, Full p-curve
  ppr=as.numeric(ifelse(p<.05,20*p,NA))
  ppr=pbound(ppr)


  # Right Skew, half p-curve
  ppr.half=as.numeric(ifelse(p<.025,40*p,NA))
  ppr.half=pbound(ppr.half)

  # Power of 33%
  ncp33=mapply(getncp,df1=df1,df2=df2,power=1/3,family=family)

  # Full-p-curve
  pp33=ifelse(family=="f" & p<.05,3*(pf(value, df1=df1, df2=df2, ncp=ncp33)-2/3),NA)
  pp33=ifelse(family=="c" & p<.05,3*(pchisq(value, df=df1, ncp=ncp33)-2/3),pp33)
  pp33=pbound(pp33)

  # half p-curve
  prop25=3*prop33(.025)
  prop25.sig=prop25[p<.05]

  #Compute pp-values for the half
  pp33.half=ifelse(family=="f" & p<.025, (1/prop25)*(pf(value,df1=df1,df2=df2,ncp=ncp33)-(1-prop25)),NA)
  pp33.half=ifelse(family=="c" & p<.025, (1/prop25)*(pchisq(value,df=df1, ncp=ncp33)-(1-prop25)),pp33.half)
  pp33.half=pbound(pp33.half)


  ##############################################################################
  # 4. Stouffer & Binomial test ################################################
  ##############################################################################

  # Convert pp-values to Z scores, using Stouffer function above
  Zppr = stouffer(ppr)
  Zpp33 = stouffer(pp33)
  Zppr.half = stouffer(ppr.half)
  Zpp33.half = stouffer(pp33.half)

  # Overall p-values from Stouffer test
  p.Zppr = pnorm(Zppr)
  p.Zpp33 = pnorm(Zpp33)
  p.Zppr.half = pnorm(Zppr.half)
  p.Zpp33.half = pnorm(Zpp33.half)

  # Save results to file
  main.results=as.numeric(c(ktot, ksig, khalf, Zppr,
                            p.Zppr, Zpp33, p.Zpp33, Zppr.half,
                            p.Zppr.half, Zpp33.half, p.Zpp33.half))

  # BINOMIAL
  # Observed share of p<.025
  prop25.obs=sum(p<.025)/sum(p<.05)
  # Flat null
  binom.r=1-pbinom(q=prop25.obs*ksig- 1, prob=.5, size=ksig)
  # Power of 33% null
  binom.33=ppoibin(kk=prop25.obs*ksig,pp=prop25[p<.05])
  # Save binomial results
  binomial=c(mean(prop25.sig), prop25.obs, binom.r, binom.33)

  # Beautifyier Function
  cleanp=function(p)
  {
    p.clean=round(p,4)           #Round it
    p.clean=substr(p.clean,2,6)  #Drop the 0
    p.clean=paste0("= ",p.clean)
    if (p < .0001) p.clean= " < .0001"
    if (p > .9999) p.clean= " > .9999"
    return(p.clean)
  }

  #If there are zero p<.025, change Stouffer values for half-p-curve tests for "N/A" messages
  if (khalf==0) {
    Zppr.half ="N/A"
    p.Zppr.half ="=N/A"
    Zpp33.half ="N/A"
    p.Zpp33.half ="=N/A"
  }


  #If there are more than 1 p<.025, round the Z and beutify the p-values
  if (khalf>0) {
    Zppr.half =round(Zppr.half,2)
    Zpp33.half =round(Zpp33.half,2)
    p.Zppr.half=cleanp(p.Zppr.half)
    p.Zpp33.half=cleanp(p.Zpp33.half)
  }

  #Clean  results for full test
  Zppr=round(Zppr,2)
  Zpp33=round(Zpp33,2)
  p.Zppr=cleanp(p.Zppr)
  p.Zpp33=cleanp(p.Zpp33)
  binom.r=cleanp(binom.r)
  binom.33=cleanp(binom.33)

  ################################################
  # 5. Power  ####################################
  ################################################

  powerfit=function(power_est)
  {
    ncp_est=mapply(getncp,df1=df1,df2=df2,power=power_est,family=family)
    pp_est=ifelse(family=="f" & p<.05,(pf(value,df1=df1,df2=df2,ncp=ncp_est)-(1-power_est))/power_est,NA)
    pp_est=ifelse(family=="c" & p<.05,(pchisq(value,df=df1,ncp=ncp_est)-(1-power_est))/power_est,pp_est)
    pp_est=pbound(pp_est)

    return(stouffer(pp_est))
  }


  fit=c()
  fit=abs(powerfit(.051))
  for (i in 6:99)   fit=c(fit,abs(powerfit(i/100)))

  mini=match(min(fit,na.rm=TRUE),fit)

  hat=(mini+4)/100

  x.power=seq(from=5,to=99)/100

  get.power_pct =function(pct)   {
    #Function that finds power that gives p-value=pct for the Stouffer test
    #for example, get.power_pct(.5) returns the level of power that leads to p=.5  for the stouffer test.
    #half the time we would see p-curves more right skewed than the one we see, and half the time
    #less right-skewed, if the true power were that get.power_pct(.5). So it is the median estimate of power
    #similarliy, get.power_pct(.1) gives the 10th percentile estimate of power...
    #Obtain the normalized equivalent of pct, e.g., for 5% it is -1.64, for 95% it is 1.64
    z=qnorm(pct)  #convert to z because powerfit() outputs a z-score.
    #Quantify gap between computed p-value and desired pct
    error = function(power_est, z)  powerfit(power_est) - z
    #Find the value of power that makes that gap zero, (root)
    return(uniroot(error, c(.0501, .99),z)$root)   }

  # Boundary conditions
  p.power.05=pnorm(powerfit(.051)) #Proability p-curve would be at least at right-skewed if power=.051
  p.power.99=pnorm(powerfit(.99))  #Proability p-curve would be at least at right-skewed if power=.99

  # Lower end of ci
  if (p.power.05<=.95) power.ci.lb=.05
  if (p.power.99>=.95) power.ci.lb=.99
  if (p.power.05>.95 && p.power.99<.95)  power.ci.lb=get.power_pct(.95)

  # Higher end of CI
  if (p.power.05<=.05) power.ci.ub=.05
  if (p.power.99>=.05) power.ci.ub=.99
  if (p.power.05>.05 && p.power.99<.05) power.ci.ub=get.power_pct(.05)

  # Save power fit
  power_results=c(power.ci.lb,hat,power.ci.ub)


  ##############################################################################
  # 6. Plot  ###################################################################
  ##############################################################################

  # Green line (Expected p-curve for 33% power)
  gcdf1=prop33(.01)
  gcdf2=prop33(.02)
  gcdf3=prop33(.03)
  gcdf4=prop33(.04)

  green1=mean(gcdf1,na.rm=TRUE)*3
  green2=mean(gcdf2-gcdf1,na.rm=TRUE)*3
  green3=mean(gcdf3-gcdf2,na.rm=TRUE)*3
  green4=mean(gcdf4-gcdf3,na.rm=TRUE)*3
  green5=mean(1/3-gcdf4,na.rm=TRUE)*3
  green=100*c(green1,green2,green3,green4,green5)


  # Blue line (observed p-curve)
  ps=ceiling(p[p<.05]*100)/100
  blue=c()
  for (i in c(.01,.02,.03,.04,.05)) blue=c(blue,sum(ps==i,na.rm=TRUE)/ksig*100)


  # Red line
  red=c(20,20,20,20,20)

  # Make the graph
  x = c(.01,.02,.03,.04,.05)
  par(mar=c(6,5.5,1.5,3))
  moveup=max(max(blue[2:5])-66,0)
  ylim=c(0,105+moveup)
  legend.top=100+moveup
  plot(x,blue,   type='l', col='dodgerblue2',  main="",
       lwd=2, xlab="", ylab="", xaxt="n",yaxt="n", xlim=c(0.01,0.051),
       ylim=ylim, bty='L', las=1,axes=F)
  x_=c(".01",".02",".03",".04",".05")
  axis(1,at=x,labels=x_)
  y_=c("0%","25%","50%","75%","100%")
  y=c(0,25,50,75,100)
  axis(2,at=y,labels=y_,las=1,cex.axis=1.2)
  mtext("Percentage of test results",font=2,side=2,line=3.85,cex=1.25)
  mtext("p            ",font=4,side=1,line=2.3,cex=1.25)
  mtext(" -value",      font=2,side=1,line=2.3,cex=1.25)
  points(x,blue,type="p",pch=20,bg="dodgerblue2",col="dodgerblue2")
  text(x+.00075,blue+3.5,percent(round(blue)/100),col='black', cex=.75)
  lines(x,red,   type='l', col='firebrick2',    lwd=1.5, lty=3)
  lines(x,green, type='l', col='springgreen4',  lwd=1.5, lty=5)
  tab1=.017          #Labels for line at p=.023 in x-axis
  tab2=tab1+.0015    #Test results and power esimates at tab1+.0015
  gap1=9             #between labels
  gap2=4             #between lable and respective test (e.g., "OBserved p-curve" and "power estimate")
  font.col='gray44'
  text.blue=paste0("Power estimate: ",percent(hat),", CI(",
                   percent(power.ci.lb),",",
                   percent(power.ci.ub),")")
  text(tab1,legend.top,     adj=0,cex=.85,bquote("Observed "*italic(p)*"-curve"))
  text(tab2,legend.top-gap2,adj=0,cex=.68,text.blue,col=font.col)
  text.red=bquote("Tests for right-skewness: "*italic(p)*""[Full]~.(p.Zppr)*",  "*italic(p)*""[Half]~.(p.Zppr.half))
  #note: .() within bquote prints the value rather than the variable name
  text(tab1,legend.top-gap1,    adj=0,cex=.85, "Null of no effect" )
  text(tab2,legend.top-gap1-gap2,  adj=0,cex=.68, text.red, col=font.col )
  text.green=bquote("Tests for flatness: "*italic(p)*""[Full]~.(p.Zpp33)*",  "*italic(p)*""[half]~.(p.Zpp33.half)*",  "*italic(p)*""[Binomial]~.(binom.33))
  text(tab1,legend.top-2*gap1,    adj=0,cex=.85,"Null of 33% power")
  text(tab2,legend.top-2*gap1-gap2,  adj=0,cex=.68,text.green,col=font.col)
  segments(x0=tab1-.005,x1=tab1-.001,y0=legend.top,y1=legend.top,      col='dodgerblue2',lty=1,lwd=1.5)
  segments(x0=tab1-.005,x1=tab1-.001,y0=legend.top-gap1,  y1=legend.top-gap1,col='firebrick2',lty=3,lwd=1.5)
  segments(x0=tab1-.005,x1=tab1-.001,y0=legend.top-2*gap1,y1=legend.top-2*gap1,col='springgreen4',lty=2,lwd=1.5)
  rect(tab1-.0065,legend.top-2*gap1-gap2-3,tab1+.032,legend.top+3,border='gray85')
  msgx=bquote("Note: The observed "*italic(p)*"-curve includes "*.(ksig)*
                " statistically significant ("*italic(p)*" < .05) results, of which "*.(khalf)*
                " are "*italic(p)*" < .025.")
  mtext(msgx,side=1,line=4,cex=.65,adj=0)
  kns=ktot-ksig
  if (kns==0) ns_msg="There were no non-significant results entered."
  if (kns==1) ns_msg=bquote("There was one additional result entered but excluded from "*italic(p)*"-curve because it was "*italic(p)*" > .05.")
  if (kns>1)  ns_msg=bquote("There were "*.(kns)*" additional results entered but excluded from "*italic(p)*"-curve because they were "*italic(p)*" > .05.")
  mtext(ns_msg,side=1,line=4.75,cex=.65,adj=0)



  ##############################################################################
  # 7 Save Calculations  #######################################################
  ##############################################################################

  # table_calc
  table_calc=data.frame(raw, p, ppr, ppr.half, pp33, pp33.half,
                        qnorm(ppr),  qnorm(ppr.half), qnorm(pp33), qnorm(pp33.half))
  headers1=c("Entered statistic","p-value", "ppr", "ppr half", "pp33%","pp33 half",
             "Z-R","Z-R half","Z-33","z-33 half")
  table_calc=setNames(table_calc,headers1)

  # table_figure
  headers2=c("p-value","Observed (blue)","Power 33% (Green)", "Flat (Red)")
  table_figure=setNames(data.frame(x,blue,green,red),headers2)


  ################################################
  # 8. Cumulative p-curves (Deprecated) ##########
  ################################################


  #7.1 FUNCTION THAT RECOMPUTES OVERALL STOUFFER TEST WITHOUT (K) MOST EXTREME VALUES, ADJUSTING THE UNIFORM TO ONLY INCLUDE RANGE THAT REMAINS
  dropk=function(pp,k,droplow)
  {
    #Syntax:
    #pp: set of pp-values to analyze sensitivity to most extremes
    #k:  # of most extreme values to exclude
    #dropsmall: 1 to drop smallest, 0 to drop largest

    pp=pp[!is.na(pp)]                             #Drop missing values
    n=length(pp)                                  #See how many studies are left
    pp=sort(pp)                                   #Sort the pp-value from small to large
    if (k==0) ppk=pp                              #If k=0 do nothing for nothing is being dropped
    #If we are dropping low values
    if (droplow==1 & k>0)
    {
      #Eliminate lowest k from the vector of pp-values
      ppk=(pp[(1+k):n])
      ppmin=min(pp[k],k/(n+1))   #Boundary used to define possible range of values after exclusion
      ppk=(ppk-ppmin)/(1-ppmin)  #Take the k+1 smallest pp-value up to the highest, subtract from each the boundary value, divide by the range, ~U(0,1) under the null
      #This is explained in Supplement 1 of Simonsohn, Simmons Nelson, JEPG 2016 "Better p-curves" paper. See https://osf.io/mbw5g/

    }

    #If we are dropping high values
    if (droplow==0 & k>0)
    {
      #Eliminate lowest k from the vector of pp-values
      ppk=pp[1:(n-k)]
      ppmax=max(pp[n-k+1],(n-k)/(n+1))  #Find new boundary of range
      ppk=ppk/ppmax                      #Redefine range to make U(0,1)
    }
    #In case of a tie with two identical values we would have the ppk be 0 or 1, let's replace that with almost 0 and almost 1
    ppk=pmax(ppk,.00001)                       #Adds small constant to the smallest redefined p-value, avoids problem if dropped p-value is "equal" to next highest, then that pp-value becomes 0
    ppk=pmin(ppk,.99999)                       #Subtract small constant to the largest  redefined pp-value, same reason
    Z=sum(qnorm(ppk))/sqrt(n-k)
    return(pnorm(Z))
  } #End function dropk


  #7.2 Apply function, in loop with increasing number of exclusions, to full p-curve
  #Empty vectors for results
  droplow.r=droplow.33=drophigh.r=drophigh.33=c()

  #Loop over full p-curves
  for (i in 0:(round(ksig/2)-1))
  {
    #Drop the lowest k studies in terms of respective overall test
    #Right skew
    droplow.r= c(droplow.r,   dropk(pp=ppr,k=i,droplow=1))
    drophigh.r=c(drophigh.r,  dropk(pp=ppr,k=i,droplow=0))
    #Power of 33%
    droplow.33=c(droplow.33,  dropk(pp=pp33,k=i,droplow=1))
    drophigh.33=c(drophigh.33, dropk(pp=pp33,k=i,droplow=0))
  }

  #Half p-curves

  if (khalf>0)
  {
    droplow.halfr=drophigh.halfr=c()
    for (i in 0:(round(khalf/2)-1))
    {
      #Drop the lowest k studies in terms of respective overall test
      droplow.halfr= c(droplow.halfr,   dropk(pp=ppr.half,k=i,droplow=1))
      drophigh.halfr=c(drophigh.halfr,  dropk(pp=ppr.half,k=i,droplow=0))
    } #End loop
  }#End if that runs calculations only if khalf>0

  #7.3 FUNCTION THAT DOES THE PLOT OF RESULTS
  plotdrop=function(var,col)
  {
    k=length(var)
    #Plot the dots
    plot(0:(k-1),var,xlab="",ylab="",type="b",yaxt="n",xaxt="n",main="",
         cex.main=1.15,ylim=c(0,1),col=col)
    #Add marker in results with 0 drops
    points(0,var[1],pch=19,cex=1.6)
    #Red line at p=.05
    abline(h=.05,col="red")
    #Y-axis value labels
    axis(2,c(.05,2:9/10),labels=c('.05','.2','.3','.4','.5','6','7','.8','.9'),las=1,cex.axis=1.5)
    axis(1,c(0:(k-1)),las=1,cex.axis=1.4)
  }


  ######################################################################################
  # 9. Effect Estimation ###############################################################
  ######################################################################################

  if (effect.estimation == TRUE){

    # Define ci.to.t function
    ci.to.t = function(TE, lower, upper, n){

      z.to.d = function(z, n){
        d = (2*z)/sqrt(n)
        return(abs(d))
      }

      ci.to.p = function(est, lower, upper){
        SE = (upper-lower)/(2*1.96)
        z = abs(est/SE)
        p = exp(-0.717*z - 0.416*z^2)
        return(p)
      }

      d.to.t = function(d, n){
        df = n-2
        t = (d*sqrt(df))/2
        return(t)
      }

      p = ci.to.p(TE, lower, upper)
      z = abs(qnorm(p/2))
      d = z.to.d(z, n)
      t = d.to.t(d, n)

      return(t)

    }

    #Function 13 - loss function
    loss=function(t_obs,df_obs,d_est) {

      #1.Convert all ts to the same sign (for justification see Supplement 5)
      t_obs=abs(t_obs)

      #2 Compute p-values
      p_obs=2*(1-pt(t_obs,df=df_obs))

      #3 Keep significant t-values and corresponding df.
      t.sig=subset(t_obs,p_obs<.05)
      df.sig=subset(df_obs,p_obs<.05)

      #4.Compute non-centrality parameter implied by d_est and df_obs
      #df+2 is total N.
      #Becuase the noncentrality parameter for the student distribution is ncp=sqrt(n/2)*d,
      #we add 2 to d.f. to get N,  divide by 2 to get n, and by 2 again for ncp, so -->df+2/4
      ncp_est=sqrt((df.sig+2)/4)*d_est

      #5.Find critical t-value for p=.05 (two-sided)
      #this is used below to compute power, it is a vector as different tests have different dfs
      #and hence different critical values
      tc=qt(.975,df.sig)

      #4.Find power for ncp given tc, again, this is a vector of implied power, for ncp_est,  for each test
      power_est=1-pt(tc,df.sig,ncp_est)

      #5.Compute pp-values
      #5.1 First get the overall probability of a t>tobs, given ncp
      p_larger=pt(t.sig,df=df.sig,ncp=ncp_est)

      #5.2 Now, condition on p<.05
      ppr=(p_larger-(1-power_est))/power_est  #this is the pp-value for right-skew

      #6. Compute the gap between the distribution of observed pp-values and a uniform distribution 0,1
      KSD=ks.test(ppr,punif)$statistic        #this is the D statistic outputted by the KS test against uniform
      return(KSD)
    }

    if(missing(N)){
      stop("If 'effect.estimation=TRUE', argument 'N' must be provided.")
    }

    if (length(N) != length(metaobject$TE)){
      stop("N must be of same length as the number of studies contained in x.")
    }

    lower = metaobject$TE - (metaobject$seTE*1.96)
    upper = metaobject$TE + (metaobject$seTE*1.96)
    t_obs = ci.to.t(metaobject$TE, lower, upper, N)
    df_obs = N-2

    #Results will be stored in these vectors, create them first
    loss.all=c()
    di=c()

    #Compute loss for effect sizes between d=c(dmin,dmax) in steps of .01
    for (i in 0:((dmax-dmin)*100))
    {
      d=dmin+i/100                   #effect size being considered
      di=c(di,d)                     #add it to the vector (kind of silly, but kept for symmetry)
      options(warn=-1)               #turn off warning becuase R does not like its own pt() function!
      loss.all=c(loss.all,loss(df_obs=df_obs,t_obs=t_obs,d_est=d))
      #apply loss function so that effect size, store result
      options(warn=0)                #turn warnings back on
    }

    #find the effect leading to smallest loss in that set, that becomes the starting point in the optimize command
    imin=match(min(loss.all),loss.all)       #which i tested effect size lead to the overall minimum?
    dstart=dmin+imin/100                     #convert that i into a d.

    #optimize around the global minimum
    dhat=optimize(loss,c(dstart-.1,dstart+.1), df_obs=df_obs,t_obs=t_obs)
    options(warn=-0)

    #Plot results
    plot(di,loss.all,xlab="Effect size\nCohen-d", ylab="Loss (D stat in KS test)",ylim=c(0,1), main="How well does each effect size fit? (lower is better)")
    points(dhat$minimum,dhat$objective,pch=19,col="red",cex=2)
    text(dhat$minimum,dhat$objective-.08,paste0("p-curve's estimate of effect size:\nd=",round(dhat$minimum,3)),col="red")


  }

  ######################################################################################
  # 10. Prepare Results for Return #####################################################
  ######################################################################################

  # Get results
  main.results = round(main.results, 3)
  ktotal = round(main.results[1]) # Get the total number of inserted TEs
  k.sign = round(main.results[2]) # Get the total number of significant TEs
  k.025 =  round(main.results[3]) # Get the number of p<0.25 TEs
  skew.full.z = main.results[4] # Get the Z-score for the full curve skewness test
  skew.full.p = main.results[5] # Get the p-value for the full curve skewness test
  flat.full.z = main.results[6] # Get the Z-score for the full curve flatness test
  flat.full.p = main.results[7] # Get the p-value for the full curve flatness test
  skew.half.z = main.results[8] # Get the Z-score for the half curve skewness test
  skew.half.p = main.results[9] # Get the p-value for the half curve skewness test
  flat.half.z = main.results[10] # Get the Z-score for the half curve flatness test
  flat.half.p = main.results[11] # Get the p-value for the half curve flatness test
  skew.binomial.p = round(binomial[3], 3) # Get the skewness binomial p-value
  flat.binomial.p = round(binomial[4], 3) # Get the flatness binomial p-value

  # Make data.frame
  skewness = c(skew.binomial.p, skew.full.z, skew.full.p, skew.half.z, skew.half.p)
  flatness = c(flat.binomial.p, flat.full.z, flat.full.p, flat.half.z, flat.half.p)
  colnames.df = c("pBinomial", "zFull", "pFull", "zHalf", "pHalf")
  rownames.df = c("Right-skewness test", "Flatness test")
  pcurveResults = rbind(skewness, flatness)
  colnames(pcurveResults) = colnames.df
  rownames(pcurveResults) = rownames.df

  # Power results
  power_results = round(power_results, 3)
  powerEstimate = power_results[2]
  powerLower = power_results[1]
  powerUpper = power_results[3]
  Power = as.data.frame(cbind(powerEstimate, powerLower, powerUpper))
  rownames(Power) = ""

  # Presence and absence of evidential value
  # - If the half p-curve test is right-skewed with p<.05 or both the half and full test
  #   are right-skewed with p<.1, then p-curve analysis indicates the presence of evidential value
  # - Evidential value is inadequate or absent if the 33% power test is p<.05 for the full p-curve
  #   or both the half p-curve and binomial 33% power test are p<.1

  if (skew.half.p < 0.05 | (skew.half.p < 0.1 & skew.full.p < 0.1)){
    presence.ev = "yes"
  } else {
    presence.ev = "no"
  }

  if (flat.full.p < 0.05 | (flat.half.p < 0.1 & flat.binomial.p < 0.1)){
    absence.ev = "yes"
  } else {
    absence.ev = "no"
  }

  # Plot Data
  PlotData = round(table_figure, 3)

  # Input Data
  table_calc[,1] = NULL
  colnames(table_calc) = c("p", "ppSkewFull", "ppSkewHalf", "ppFlatFull", "ppFlatHalf", "zSkewFull", "zSkewHalf",
                           "zFlatFull", "zFlatHalf")
  Input = cbind(metaobject$TE, round(table_calc,3))
  rownames(Input) = paste(1:length(metaobject$TE), metaobject$studlab)
  colnames(Input)[1] = "TE"


  if (effect.estimation==TRUE){

    dEstimate = round(dhat$minimum, 3)
    return.list = list("pcurveResults" = pcurveResults,
                   "Power" = Power,
                   "PlotData" = PlotData,
                   "Input" = Input,
                   "EvidencePresent" = presence.ev,
                   "EvidenceAbsent" = absence.ev,
                   "kInput" = ktot,
                   "kAnalyzed" = k.sign,
                   "kp0.25" = k.025,
                   "dEstimate" = dEstimate,
                   "I2" = metaobject$I2,
                   "class.meta.object" = class(metaobject)[1])

    class(return.list) = c("pcurve", "effect.estimation")

  } else {

    return.list = list("pcurveResults" = pcurveResults,
                   "Power" = Power,
                   "PlotData" = PlotData,
                   "Input" = Input,
                   "EvidencePresent" = presence.ev,
                   "EvidenceAbsent" = absence.ev,
                   "kInput" = ktot,
                   "kAnalyzed" = k.sign,
                   "kp0.25" = k.025,
                   "I2" = metaobject$I2,
                   "class.meta.object" = class(metaobject)[1])

    class(return.list) = c("pcurve", "no.effect.estimation")
  }

  cat("  ", "\n")
  invisible(return.list)
  return.list

}

For the pcurve() function, the following parameters need to be specified:

Parameter	Function
x	The meta-analysis results object generated by meta functions.
effect.estimation	Logical. Should the true effect size underlying the p-curve be estimated? If set to TRUE, a vector containing the total sample size for each study must be provided for N. FALSE by default.
N	A numeric vector of same length as the number of effect sizes included in x specifiying the total sample size N corresponding to each effect. Only needed if effect.estimation = TRUE.
dmin	If effect.estimation = TRUE: lower limit for the effect size (d) space in which the true effect size should be searched. Must be greater or equal to 0. Default is 0.
dmax	If effect.estimation = TRUE: upper limit for the effect size (d) space in which the true effect size should be searched. Must be greater than 0. Default is 1.

We plug the re_mod meta-analysis object into the pcurve() function to generate the p-curve.

pcurve(re_mod)

## P-curve analysis 
##  ----------------------- 
## - Total number of provided studies: k = 64 
## - Total number of p<0.05 studies included into the analysis: k = 41 (64.06%) 
## - Total number of studies with p<0.025: k = 36 (56.25%) 
##    
## Results 
##  ----------------------- 
##                     pBinomial   zFull pFull   zHalf pHalf
## Right-skewness test     0.000 -20.773     0 -20.925     0
## Flatness test           0.997  16.469     1  20.696     1
## Note: p-values of 0 or 1 correspond to p<0.001 and p>0.999, respectively.   
## Power Estimate: 99% (99%-99%)
##    
## Evidential value 
##  ----------------------- 
## - Evidential value present: yes 
## - Evidential value absent/inadequate: no

How to interpret pcurve results:
Among all studies, if there is no true effect, you would expect a uniform distribution of p-values (shown on the plot with a red dotted line). That is, if researchers contrived p-values using one of the dishonest methods, all you would expect to see is a bunch of p-values under 0.05, with no true pattern.
If there is a true effect, you would expect a right-skewed distribution of p-values (shown on the plot with a blue line). That is, if p-values are not contrived, and there is a true effect, you should expect the majority of studies to have low p-values, and increasingly fewer studies with higher p-values.
The plot also shows the distribution of effects you would expect to see if all studies were low-power studies (shown on the plot with a green dashed line).
If authors of our data were guilty of data-mining or p-hacking, wherein they add observations until they attain a “significant” result, then we would expect to observe a left-skewed distribution. This assumes that researchers would stop p-hacking when they achieve a “significant” p-value of 0.05, leaving us with an excess of p-values close to 0.05.

Challenge 5

Refer to the results of our funnel plot and pcurve to answer the following questions:
1. Based on the results of our funnel plot, is small-study publication bias present in our meta-analysis?
2. Based on the results of our pcurve, is academic dishonesty in the form of p-hacking present in our meta-analysis?

That was a short introduction to meta-analysis. If you wish to learn more, please consult the resources listed below. This exercise drew heavily from Harrer et al. (2019).

References

Alldred, Mary 2016. Denitrification measurements in various plant communities. figshare. Dataset. https://doi.org/10.6084/m9.figshare.1541106.v1.

Alldred, M. and S.B. Baines. 2016. Effects of wetland plants on denitrification rates: a meta-analysis. Ecological Applications 26(3): 676-685. https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/14-1525.

Gurevitch, J. and L.V. Hedges. 2001. Meta-Analysis: Combining the results of independent experiments. in S.M. Scheiner and J. Gurevitch, editors. Design and Analysis of Ecological Experiments. Oxford University Press. New York, NY, USA.

Harrer, M., P. Cuijpers, T.A. Furukawa, and D.D. Ebert. 2019. Doing Meta-Analysis in R: A Hands-on Guide. https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/. DOI 10.5281/zenodo.2551803.

Nikolakopoulou, A., D. Mavridis, and G. Salanti. 2014. Demystifying fixed and random effects meta-analysis. Evidence-Based Mental Health 17(2): 53-57. http://dx.doi.org/10.1136/eb-2014-101795.

Assignment 14: Meta-Analysis in R

Mary Alldred

Last updated 2024/01/25

Learning Objectives

Why Meta-Analysis?

Getting started

Installing Packages

Loading packages

Formatting Data for Meta-Analysis

Importing Data for Today’s Exercise

Data Dictionary for Today’s Dataset

Pooling Effect Sizes: Fixed v. Random Effects

Choosing an Appropriate Effect Size

Fixed Effects Model

Challenge 1

Random Effects Model

Challenge 2

Subgroup Analysis

Challenge 3

Challenge 4

Forest Plots

Publication Bias

Small Sample Bias: Funnel Plots

P-Curve

Challenge 5

References