A simulations approach for meta-analysis of genetic association studies based on additive genetic model

# A simulations approach for meta-analysis of genetic association studies based on additive genetic model

Majnu John111Corresponding author: Community Drive, Manhasset, NY 11030. e-mail: mjohn5@northwell.edu, Phone: +01 718 470 8221, Fax: +01 718 343 1659, Todd Lencz, Anil K Malhotra, Christoph U Correll, Jian-Ping Zhang Center for Psychiatric Neuroscience,
The Feinstein Institute of Medical Research,
Manhasset, NY.
Psychiatry Research, Zucker Hillside Hospital,
Northwell Health System,
Glen Oaks, NY.
Department of Mathematics,
Hofstra University,
Departments of Psychiatry and of Molecular Medicine,
Hofstra Northwell School of Medicine,
###### Abstract

Genetic association studies are becoming an important component of medical research. To cite one instance, pharmacogenomics which is gaining prominence as a useful tool for personalized medicine is heavily reliant on results from genetic association studies. Meta-analysis of genetic association studies is being increasingly used to assess phenotypic differences between genotype groups. When the underlying genetic model is assumed to be dominant or recessive, assessing the phenotype differences based on summary statistics, reported for individual studies in a meta-analysis, is a valid strategy. However, when the genetic model is additive, a similar strategy based on summary statistics will lead to biased results. This fact about the additive model is one of the things that we establish in this paper, using simulations. The main goal of this paper is to present an alternate strategy for the additive model based on simulating data for the individual studies. We show that the alternate strategy is far superior to the strategy based on summary statistics.

###### keywords:
genetic association studies, meta analysis, additive model, summary statistics, simulations
journal: arXiv

## 1 Introduction

Over the last decade, genetic association studies (-both candidate based designs and genome-wide designs) have become one of the cornerstone techniques in detecting specific genetic variants related to any phenotype of interest. In association studies, genetic variation is typically measured using genotyping based on single nucleotide polymorphisms (SNPs). In order to fix ideas, let us assume that a SNP of interest is named and the corresponding genotypes AA, AB and BB, with A as the major allele and B as the minor allele (the risk allele). The phenotype of interest in our motivating example-studies was weight gain due to antipsychotic treatment. With the proliferation of association studies, it is often seen that the standardized differences in weight gain (“effect sizes”) between any pair of genotypes vary across the studies. For this reason, meta-analysis has become the method that is being increasingly employed to assess the phenotypic differences across the genotypic groups.

Combining the effect sizes between genotype groups in a meta-analysis depends on the underlying genetic model. The most commonly considered genetic models are dominant, recessive and additive models. In a dominant model, allele B increases risk and the number of copies of B doesn’t matter; that is, genotypes AB and BB carry the same risk. In the recessive model, two copies of allele B are required for increased risk; that is, AB is grouped together with AA, and both groups are assumed to have no risk. In an additive model, the increase in risk is proportional to the number of copies of allele B; that is, if the risk is increased -fold for the genotype group AB, then it is increased -fold for BB.

Typically, in the association studies considered for the meta-analysis, summary statistics (mean, median, standard deviation etc.) for weight gain will be reported for each of the three genotype groups. If we consider the underlying model to be dominant, one may take the average of the mean (or median) weight gains reported for the groups AB and BB, and calculate the standardized difference between this average and the corresponding mean (or median) weight gain reported for AA, to get the effect size for each study. It is ideal to weight the average for AB and BB by their respective sample sizes. This is a valid strategy in the sense that if the meta-analyst had access to each study’s data, and if he or she grouped together the weight gains for the individual patients in groups AB and BB, and took the standardized mean difference with the individual data from AA group, then the effect size obtained will be same as the one obtained using summary statistics mentioned above; at least conceptually same - practically there might be some differences, for example, if one did not weight by sample sizes. The same argument applies for the recessive model too. However, for the additive model, the results from crude approach based on combining summary statistics will not correspond with the results using the approach based on the individual data, if the analyst was to have access to the individual data. The crude approach for the additive model would be roughly as follows. Take the difference of the summary statistics (mean or median) for AB and BB, then do similarly for AA and AB, and then take the average of these two standardized pairwise differences to obtain the corresponding effect size. One of the things that we show in this paper, via simulations, is that this crude approach leads to biased results.

Since the approach based on summary statistics leads to biased results for the additive model, the analyst will have to fall back on subject level data from each study to calculate the appropriate effect size. But, it is well known that very rarely does a meta-analyst have access to such data. The solution to this dilemma, that we suggest in this paper, is to generate subject level data via simulations based on the summary level data, in order to calculate the effect size for each study under the additive genetic model. We show via extensive simulations and real data that the simulations-based-approach that we propose is far superior to the crude approach, and gives results very similar to those that might be obtained if individual level data was available. Before causing any confusion, we would like to clarify the word ‘simulations’ used twice in the previous sentence. The new method proposed in this paper is based on simulations, and we compare this new method with the old method (that is, the crude approach based on summary statistics) also using simulations.

In section 2, we explain the methods more clearly using notations. In section 3, we explain the simulations comparing the two methods, and in section 4, we present the results of the simulations. Section 5 illustrates how well the simulations-based-method works for a real data example; that is, for a data set consisting of studies for which we had subject level data. In sections 2 to 5, the phenotype of interest (e.g. weight gain) was considered as a continuous variable. But, phenotypes could be dichotomous as well (e.g. overweight - yes/no). In section 6, we show how our approach can be extended for dichotomous phenotypes as well. Finally, we summarize our findings in the last section.

## 2 Methods

In this section we describe two methods for meta-analysis assuming additive genetic model in genetic association studies. We are interested in the additive effect of the 3 genotypic groups on a phenotype of interest. In our motivating example, the phenotype of interest was the change in body weight due to antipsychotic treatment in patients with schizophrenia. We assume that the means and standard deviations of the phenotype of interest are given for the three groups along with the sample sizes .

Crude Approach. One way to find the additive effect would be to stack the three means and into a column vector and regress it against a column vector of group indicators . The thus obtained is just the average of the pairwise mean differences and . The R codes corresponding to this simple approach would be

M.crude c(m1, m2, m3) #means
G.crude c(1, 2, 3) #group indicator variable
beta.crude summary(lm(M.crude~G.crude))$coeff[2,1] In order to calculate the standard deviation of (that is, the denominator of Cohen’s ), there is no straightforward approach. However, one could take the standard deviation for the mean difference between samples 1 and 2, , and similarly for samples 2 and 3, , and take the average of these two standard deviations, as the standard deviation for . and are within group standard deviations, pooled across the corresponding groups. As mentioned by Borenstein and co-authors, the reason to pool is that even if the underlying population standard deviations are the same, it is unlikely that the sample estimates and (used in the calculation of ) and or (for ) will be equal. The R codes for this part are SD.crude c(sd1, sd2, sd3) SS.crude c(n1, n2, n3) sd12 sqrt((((SS.crude[1]-1)*SD.crude[1]^2) + ((SS.crude[2]-1)*SD.crude[2]^2))/(SS.crude[1] + SS.crude[2] - 2)) sd23 sqrt((((SS.crude[2]-1)*SD.crude[2]^2) + ((SS.crude[3]-1)*SD.crude[3]^2))/(SS.crude[2] + SS.crude[3] - 2)) sd.beta.crude mean(c(sd12, sd23)) We use the following well-known formula to calculate an approximate variance for Cohen’s :  n1+n2n1n2+d22(n1+n2). Note that the above formula can be used only for a pair of samples. So, we use this formula to calculate and separately for the sample pairs 1 and 2, and 2 and 3. Similarly the correction factor  J=1−34(n1+n2−2)−1, for the conversion of to Hedge’s is also used separately for the two pairs to get and . Then and variance of for each pair are calculated as  g12.crude=J12.crude×dandg23.crude=J23.crude×d,  Vg.12.crude=J212.crude×Vd.12.crudeandVg.23.crude=J223.crude×Vd.23.crude. Finally, and its variance are obtained as a weighted average,  gcrude=(g12.crudeVg.12.crude+g23.crudeVg.23.crude)/(1Vg.12.crude+1Vg.23.crude). Note that this weighted average was suggested by Hedges and co-authors for combining effect size estimates. Simulations Approach. Another approach to the same problem would be based on simulations. In this approach, we generate individual data from the normal distribution via simulations using, for example, the rnorm function in R, based on the mean, standard deviation and sample size for each genotypic group, reported for each study. A group indicator variable is created with ’s, ’s and ’s, and the simulated individual data is regressed against the indicator variable to obtain , the numerator of Cohen’s . A one-way ANOVA test comparing the individual data across the 3 genotypic groups is then conducted. The denominator of Cohen’s for this approach is obtained as the square root of the ratio of between-mean-squares and the statistic. This process is repeated 10,000 times and the average of Cohen’s ’s across all the 10,000 iterations is taken as the estimated Cohen’s for this approach. The corresponding R codes are given below. Sim.Num 10000 ## number of iterations beta.sim array(, Sim.Num) sd.sim array(, Sim.Num) for(m in 1:Sim.Num) { sample1 rnorm(n1, m1, sd1)  sample2 rnorm(n2, m2, sd2)  sample3 rnorm(n3, m3, sd3)  M c(sample1, sample2, sample3)  G c(rep(1, n1), rep(2, n2), rep(3, n3))  beta.sim[m] summary(lm(M~G))$coeff[2,1]
 sd.sim[m]  sqrt(anova(lm(M~G))[1,3]/anova(lm(M~G))[1,4]) }

cohen.d.sim  mean(beta.sim/sd.sim)

The calculations for variance for Cohen’s , the correction factor for conversion from to , and then for and its variance are all done in a pairwise manner as was done for the crude approach (See above).

## 3 Monte Carlo Simulations

In order to assess the performance of the two methods described above under various scenarios, we conducted Monte Carlo simulations. By “various scenarios” we mean the scenarios obtained by varying the following parameters: number of studies in the meta analysis, the underlying distribution from which the sample for each study was obtained, the sample sizes for each study, and the true underlying effect size for each study. The last parameter was varied by varying the means for the three genotypic groups in each study, as well the within-study standard deviation. Our goals were a) to see whether the newly proposed methods improved with larger number of studies in the meta-analysis, b) whether they improved with larger sample sizes for the genotypic groups in the studies, c) whether the performance depended on either the underlying distribution or the underlying true effect sizes. It was intuitive to hypothesize that the performance of the methods would improve with increasing the number of studies and the sample sizes. So, our questions for parts a) and b) were really how large was good enough.

The key step in the simulations-based method presented above is the generation of (simulated) individual data from normal distributions using the “reported” means, standard deviations and sample sizes of the three genotypic groups from each study. (Note that, for the simulation analysis presented in this section, there was actually no reported means or standard deviations, but these values were generated in such a way that they matched with reported values from one the real data sets (EUFEST) discussed in section 5. See below the discussion related to mean.vec for more details.) Thus, the first question that we considered was whether the performance of the simulations-based method was affected when the original data distribution was not normally distributed. For this purpose, we conducted Monte Carlo simulations in which the original data were sampled from 3 distributions other than normal: strongly right skewed, asymmetric bimodal and heavily kurtotic. For comparison purposes, we assessed the performance of the methods, by sampling the original data from the normal distribution also. The shapes of the probability density functions for the four distributions considered are shown in Figure 1. If denote the density function for a normal distribution with mean and variance , the formulas for the densities used for our simulation study are

 f1:g(0,1)(standard normal% )f2:7∑l=018g(3{(23)l−1},(23)2l)(strongly right % skewed)f3:34g(0,1)+14g(32,(13)2)(asymmetric bimodal)f4:23g(0,1)+13g(0,110)(heavily kurtotic)

Each of the above densities is a normal mixture. These densities correspond to the first, third, eighth and fourth densities considered in Marron and Wand (1992), and have been used in simulation studies previously (for example, in John and Priebe (2007)).

Once we chose a distribution from among , , and , we then selected the number of studies, , to be included in the simulation analysis. We considered scenarios with and . After the distribution and number of studies were fixed, the next task was to select a “mean vector”, mean.vec, from among the following 3 triplets: and , and then to select a within-studies standard deviation to be either or . The approximate effect size corresponding to the various mean.vec and are given in table 1 below.

 Table 1. Effect size table σWS M1 M2 M3 Approximate true value of g 1 4 5.5 7 0.82 4 5.5 9 1.28 4 5.5 11 1.64 5 4 5.5 7 0.30 4 5.5 9 0.48 4 5.5 11 0.65

In the next step, we generated means for the genotypic group () from a normal distribution with mean as the component of the mean.vec and standard deviation equals 2, using, for example, the function in R. Thus the average mean triplet across all studies will be the mean.vec that we selected for this scenario, but we create variation across the studies by generating the mean triplet for each study from a normal distribution. Similarly, we generated (within-study) standard deviations for the genotypic group () from normal distributions using the component of the selected. The codes for these steps with , mean.vec and are given below.

L 10 ## number of studies

mean.vec c(4, 5.5, 9)
 m1 rnorm(L, mean.vec[1], 2); m1[m1 < 0] mean.vec[1] ## L means for the 1st genotypic group
 m2 rnorm(L, mean.vec[2], 2); m2[m2 < 0] mean.vec[1] ## L means for the 2nd genotypic group
 m3 rnorm(L, mean.vec[3], 2); m3[m3 < 0] mean.vec[1] ## L means for the 3rd genotypic group

sd.ws 1
 sd1 rnorm(L, sd, 2); sd1[sd1 < 0] sd.ws ## L within-study SD’s for the 1st genotypic group
 sd2 rnorm(L, sd, 2); sd2[sd2 < 0] sd.ws ## L within-study SD’s for the 2nd genotypic group
 sd3 rnorm(L, sd, 2); sd3[sd3 < 0] sd.ws ## L within-study SD’s for the 3rd genotypic group

With the means and standard deviations selected for the studies, we generate the three samples corresponding to the three genotypic groups for each study. The sample size considered at this step was also varied. We considered 8 sample size triplets (each triplet giving sample sizes for each of the three allelic groups), ranging from small to medium to very large .

Each of the scenarios described above was repeatedly generated 500 times (that is, we used 500 Monte Carlo iterations). At each iteration, the absolute difference between the values from the original sample, and the values obtained from either of the new approaches proposed in this paper was first calculated, and the mean difference across all studies was then calculated. The mean difference thus obtained was considered as the bias for the values corresponding to either of the new methods at each iteration. The average of the above bias across all Monte Carlo iterations was presented as the mean absolute bias for in the appendix tables A1 to A4. In the above case, the bias was calculated for the for each study.

At each iteration in our Monte Carlo analysis, we also conducted a meta analysis based on random effects model to calculate a mean effect size, -WM (that is, a weighted average of the ’s across all studies) from the original sample and a mean effect size based on the ’s obtained via either of the newly proposed methods. The absolute difference between the two mean effect sizes was averaged across all iterations and was labeled as the mean absolute bias for -WM. Thus, across different scenarios, we compared the bias of the ’s for each of the individual studies as well as the bias of the weighted-averaged estimate of obtained by a meta-analysis using the random effects model.

## 4 Simulation Results

Overall, the performance of the simulations based method was better than the crude approach, in all cases that we considered. In general, the performance of the simulations based approach improved with larger number of studies and larger sample sizes for the studies, although the same cannot be said for the crude approach. The simulations based method performed fairly well even when the underlying density was strongly skewed or heavily kurtotic, but when the underlying density substantially deviated from normal (as in the case of asymmetric bimodal density), the performance was slightly worse. The simulations based approach performed consistently well for all underlying effect sizes shown in table 1; however, the performance of the crude approach substantially worsened for larger effect sizes. We provide more details below.

Performance based on the number of studies. As hypothesized, the performance of both the methods improved, in general, when larger number of studies were included in the simulation analyses. For example, when the original data distribution was asymmetric bimodal, with a within-study standard deviation of 1, and with the means, and the corresponding sample sizes of the three genotypic groups being and , respectively, the bias for -WM (the weighted average of effect sizes obtained using a meta-analysis based on random effects model) was , when the simulations based method was used and the number of studies was . For the same setting, when the number of studies was increased to , the corresponding bias reduced to (about reduction), and when the number of studies was further increased to , the bias was further reduced to (about reduction). A similar trend was observed for the crude approach also: the bias values corresponding to study-sizes of , and were respectively , and , but the percent reductions, which amounted to and were much lower. Although the above pattern was observed in most of the scenarios there were certain cases, especially when the sample sizes were low, where the crude approach did worse with increasing number of studies. For example, when the original data distribution was strongly right skewed, with a within-study standard deviation of , and with the means, and the corresponding sample sizes of the three genotypic groups being and , respectively, the bias for -WM went up from to and then changed to as the number of studies went up from to , and then to , when the crude approach was used. On the other hand, for the same setting, the corresponding bias values for the simulations based approach decreased from to ( reduction), and further to ( reduction) as the number of studies increased from to , and then to .

Performance based on the sample sizes within studies. When all other parameters were kept the same, the bias values for both the crude approach and the simulations based approach went down, in general, as the sample size values within each study was increased. For example, when the original data was generated from the strongly right skewed density, with the three genotypic group means as , with within-study standard deviation as and number of studies equals , the mean absolute bias for -WM based on the simulations approach decreased from to ( reduction), and then to ( reduction), as the within-study sample size triplets increased from (small) to (medium), and then to (large); the corresponding bias values based on the crude approach were , ( reduction) and ( reduction), respectively (See appendix table A2b). However there were a few scenarios for which the bias values based on the crude approach decreased initially as the sample sizes increased from small to medium, but then increased as the sample sizes increased from medium to large, even as the bias based on the simulations based approach reduced steadily as the sample sizes increased. The following example based on appendix table A1c illustrates the above point. With the underlying data normally distributed, with the three genotypic group means as , with within-study standard deviation as and number of studies equal to , the bias values based on the crude approach based on the sample size triplets , and were respectively, , and (that is, an initial reduction of , but then an increase of ). Within the same setting and the same sample sizes, the bias based on the simulations approach decreased steadily from to , and then to .

Performance based on the underlying densities. As expected, when the underlying data distribution is normal, both the methods work better than when the data distribution is not normal. When the underlying distribution was asymmetric bimodal, the bias values of both the methods were somewhat larger compared to other distributions. When the underlying distribution was either strongly right skewed or heavily kurtotic, the bias values were slightly larger in general, yet comparable to those from the normal data distribution. For example, with a within-study standard deviation of , with the means for the three genotypic groups as , and the sample sizes as , the mean absolute bias for -WM, using the crude approach, for the four densities (normal), (strongly right skewed), (asymmetric bimodal) and (heavily kurtotic) are respectively and ; for the simulations based approach the corresponding four bias values were respectively and . As mentioned above, the values under and were similar to those for , but values for were much larger. The pattern remained the same even for larger sample sizes. With all the parameters exactly same as in the above example, but with sample sizes , the four bias values corresponding to the densities to for the crude approach were and , and for the simulations based approach were and .

Performance based on effect sizes. For small effect sizes, although the simulations based method edged out slightly over the crude approach, the performance of both the methods were more or less the same. For example, when the underlying density was strongly right skewed and number of studies was (appendix table A2b), with a within-study standard deviation of and means (that is, an effect size approximately ; see table 1), when the sample size triplet was , the bias value for the crude approach was and for the simulations based approach was . As the sample size increased for the same underlying effect size, the performance of both the methods improved substantially, even though the simulations based method still edged out a little bit - for the same parameters as above but with the sample size triplet as , the corresponding bias values were and . As the underlying effect size increased, the performance of the crude approach got markedly worse than that of the simulations based approach, as seen by the examples that follow. When the underlying density was strongly right skewed and number of studies was , with a within-study standard deviation of and means (that is, an effect size approximately ), when the sample size triplet was , the bias value for the crude approach was and for the simulations based approach was . With same underlying density, same number of studies and same sample size triplet as in the above example, but with the within-study standard deviation as and the mean triplet as (effect size approximately ), the corresponding bias values for the two methods were and .

## 5 Real Data Example

In order to further assess the performance of the methods proposed in this paper, we applied them to a real data example for which the subject-level data was also available for the studies included in the meta analysis. This example is part of the analyses reported by Zhang and co-authors. Here we include only pertinent details sufficient to illustrate the methods proposed in the present paper and assess their performances. More details can be found in Zhang et al’s paper. The goal of meta-analysis was to determine specific genetic variants associated with weight gain related to antipsychotic drugs. Primary outcome was change in body mass index (BMI). Literature search with inclusion/exclusion criteria narrowed down reports from independent samples. For among these independent cohorts, patient-level data were available. The three cohorts were the Second-Generation Antipsychotic Treatment Indications, Effectiveness and Tolerablity in Youth (SATIETY) study, European First Episode Schizophrenia Trial (EUFEST), and Zucker Hillside Hospital First Episode Schizophrenia Trial (ZHH-FE). Since patient-level data were available for these three cohorts, we could directly assess the additive effect of the alleles and compare it with the additive effects obtained from the crude and simulations based approaches (proposed in this paper) which utilizes reported summary means, standard deviations and sample sizes only. In the analysis reported in Zhang et al, SNP in the ADRA2A (adrenoceptor alpha 2A) gene was significantly associated with antipsychotic drug induced weight gain across 6 studies, including the 3 studies mentioned above, with the G allele increasing the risk. The SNP is located in the upstream of ADRA2A, and may be a binding site for transcription factors. The reported summary data for the three studies mentioned above, for the allelic groups for the SNP are presented in the table 2 below; the subscripts 1, 2 and 3 correspond to the genotypes , and .

 Table 2. Reported Summary Statistics for the 3 cohorts in the real data example Study/Cohort M1 M2 M3 SD1 SD2 SD3 SS1 SS2 SS3 SATIETY 11.45 12.16 14.73 8.29 8.38 9.63 63 63 42 EUFEST 4.04 5.35 4.67 5.11 5.88 6.44 74 40 9 ZHH-FE 3.24 2.44 3.64 2.11 1.23 2.42 25 24 21

Table 3 below shows the ’s (numerator of Cohen’s ), standard deviation of the ’s (denominator of Cohen’s ) and the Cohen’s ’s based on the patient level data, the crude approach and the simulations based approach. The values obtained using the simulations approach match very closely with the values calculated using patient level data, while as the values based on crude approach differed substantially. All the ’s and all the Cohen’s ’s from the simulations approach differed from those from the patient-level data only by or less (percent error ranging from to ). On the other hand, the percent error for crude approach ranged from to .

Table 3. Effect Sizes based on patient-level data and the two methods proposed in the paper
Patient Level Data Crude Approach Simulations Approach
Study/Cohort

SATIETY
1.552 8.661 0.179 1.625 8.675 0.187 1.563 8.680 0.180
EUFEST 0.746 5.460 0.137 0.313 5.470 0.057 0.742 5.474 0.136
ZHH-FE 0.168 2.009 0.084 0.199 1.965 0.101 0.171 2.009 0.085

## 6 The case when the phenotype of interest is a binary variable.

The methods that we have discussed so far are applicable only when the phenotype of interest is a continuous variable, like for example, change in body weight. Sometimes, the phenotype of interest is a binary variable, and in such cases odds ratios or relative risks comparing the genotypic groups are reported. That is, if we denote the genotypic groups as, say, AA, AB and BB, then usually an odds ratio for the dichotomous phenotype of interest comparing AB and AA, and another comparing BB and AB are reported, and typically these two odds ratios differ from each other. Now the question is, if we do a meta-analysis with additive genetic model as the underlying model, how do we combine the pair of odds ratios reported within a study to get a single odds ratio for the additive model within that study. In this section, we present methods that addresses this question. Although we focus on odds ratios to explain our method, a very similar method can be worked out for relative risk as well. The method that we present here is based on a simple adaptation of the reconstruction of four-fold cell frequencies for meta-analysis by Di Pietrantonj. This reconstruction was further studied in Veroniki and co-authors8.

An overview of our adaptation is easy to describe:

Step 1): Within each study, reconstruct the table for each odds ratio using Di Pietrantonj’s method. Thus, after reconstruction, with our above notation, there will be a table comparing the dichotomous phenotype of interest between AB and AA, and another one comparing BB and AB. Note that Di Pietrantonj’s method can possibly give rise to two tables for each odds ratio. Di Pietrantonj, and Veroniki (and co-authors) select the correct one among these two, for each odds ratio, using a cut-off based on the event rate in the ”treatment group or the control group” (- in our case, this would be event rate in either of the genotype group). Since, in genetic association studies, this event rate is rarely reported, we use a different approach to select the correct table for each odds ratio. More details are provided further down.

Step 2): Once the tables, one for each odds ratio within a study, are selected, we merge them to get a table. That is, three rows, one for each genotypic group, and two columns, one for each category of the dichotomous variable.

Step 3): Generate (via simulations) two variables - 1) the binary variable for the phenotype of interest and 2) a categorical indicator variable for the three genotype groups - that reflect exactly the cell frequencies in the table obtained in step 2 above.

Step 4): Run a logistic regression with (the logit of) the simulated binary outcome variable as the dependent variable and the simulated 3-genotype-groups indicator as the independent variable. Exponentiate the co-efficient for the independent variable to get the combined odds ratio for the genetic additive model. Standard errors for the can be utilized to calculated the confidence intervals for the combined odds ratio.

We explain and illustrate the details of the above approach using an example. Before we get to the details of our approach, we describe the example first. We generated sample data for each of the three genotype groups (AA, AB, BB) with sample size 30 in each group, that eventually yielded the following tables:

Table 4a. AB phenotype present phenotype absent 18 12 10 20
Table. 4b. BB phenotype present phenotype absent 18 12 18 12

Details of how we generated the data and how we obtained the above tables from the data are somewhat irrelevant to our discussion, but a reader interested in these can find the details in the appendix B.

Thus, in our example, the odds ratio comparing AB with AA is and the odds ratio comparing BB with AB is . If we use the formula

 √(1a+1b+1c+1d)

for the standard error of the log-odds-ratio (here, and denote the cell frequencies of a generic table), then we get the corresponding confidence intervals for the two odds ratios as and . Typically, the above odds ratios and confidence intervals are the only information reported when the results from the study’s analysis are published (and hence the only information available for the meta-analyst). The tables shown in Table 4 above are typically not reported.

However, someone who has the original data for this study (that is, the data that we simulated for this example and that we pretend as from an original study) could use a logistic regression with the binary phenotype (present/absent) variable as the dependent variable and the 3-genotype-groups indicator as the independent variable and then exponentiate the corresponding to obtain the odds ratio with confidence interval for the additive genetic model as (or with more digits). If this (combined) odds ratio is available for the meta-analyst, then s/he doesn’t have to look further (that is, s/he won’t need the method described in this section). However, the combined odds ratio is seldom reported, but the two separate odds ratios (and their confidence intervals) for pairs of genotypic group comparisons are the ones that are typically reported. The method described in this section would help the meta-analyst to recover/estimate the unknown (-unknown to the meta-analyst, but not perhaps to the original data analyst-), combined odds ratio (, in our example) and the confidence interval, from the two reported odds ratios (3.00 for AB vs. AA and 1.00 for BB vs. AB) and their corresponding confidence intervals.

For a generic table given below,

Generic phenotype phenotype genotype group
Table present absent totals
Genotype group 2 a b
Genotype group 1 c d

Di Pietrantonj’s formula for recovering is

 a=−λ±√λ2−4αγ2α,

where

 α=[(1−OR)2+ORm2^SE2ln(OR)]
 λ=ORm1[2(1−OR)−m2^SE2ln(OR)]
 γ=[ORm1(ORm1+m2)].

OR in these formulas denote the odds ratio obtained from the generic table and is the estimate of the standard error of obtained from the upper limit, UL, and lower limit, LL of the confidence interval, as

 ln(UL)−ln(LL)2×1.96.

Based on the estimated, we can obtain and as

 b=m1−a;c=am2ORm1+a(1−OR);d=ORm2(m1−a)ORm1+a(1−OR).

Note that there are two possible solutions for , leading to two possible values for and , and hence two possible estimates for the unknown table. Let us denote the cell entries in the two possible tables as and . As mentioned above, Di Pietrantonj and later Veroniki and co-authors devised methods based on event rates for one of the marginal groups that would help to decide on which among these possible tables to be picked. Such event rates are sometimes reported for epidemiological studies or clinical trials, but rarely for genetic association studies. We devised a simple scheme, more suitable to our context, that’d help us recover the correct table.

Recall that there are two odds ratios that we are considering - one for AB vs. AA and the other for BB vs. AB, denoted from now on as and , respectively. For each of these odds ratios, Di Pietrantonj’s formulas give two tables. Thus in total there are four tables. Our scheme for selecting the correct tables is best explained by working it out for the example that we are considering at the moment. For this example, the two possible tables for are

Table 5a. AB phenotype present phenotype absent 18 12 10 20
Table. 5b. AB phenotype present phenotype absent 20 10 12 18

and for are

Table 6a. BB phenotype present phenotype absent 12 18 12 18
Table 6b. BB phenotype present phenotype absent 18 12 18 12

There are four different ways to match the tables:

Table 5a Table 6a,    Table 5a Table 6b,    Table 5b Table 6a,    Table 5b Table 6b.

A careful look at all the tables will immediately reveal that among the four combinations that match a table corresponding to with one of the tables corresponding to , the only one that makes sense (that is, plausible) is Table 5a Table 6b, because this is the only pair where the rows corresponding to AB match. In all the other three table-matching-combinations, the estimates for the row corresponding to the AB genotype group doesn’t match.

In this particular example, there was one combination (Table 5a Table 6b) where the rows for AB matched exactly. But, there could be other examples, where the AB-rows do not match exactly for any of the combinations. In such cases, we can calculate the “distance” for the corresponding AB-rows for all four table match combinations listed above (e.g. the “distance” could be the Euclidean distance, if one may consider the pair of values for the AB-row as a 2-dimensional vector), and select the combination where this distance is the minimum. Once we decide on the table-combination based on this criteria, then we can take the column-average of the corresponding AB-rows to get the AB-row for the final table. The distances calculated for our example is given below:

Distance table. Distance between the corresponding AB rows
Table 5a Table 6a
Table 5a Table 6b
Table 5b Table 6a
Table 5b Table 6b

Again, based on this minimum distance criterion, the combination that we will select is Table 5a Table 6b. In this particular example, the AB-rows for Table 5a and Table 6b are the same, and . Hence the average is also , so that by merging Table 5a and Table 6b, we get the final table as

Table 7. BB phenotype present phenotype absent 18 12 18 12 10 20

This completes step 2. In the next step we generate a phenotype (present/absent) variable and a 3-genotype-group indicator variable using the following R codes:

al.gp.ind c(rep(1, 30), rep(2, 30), rep(3, 30))
 # 30 is the sample size for each genotype group
ev.ind c(sample(c(rep(1, 10), rep(0, 20))),
 sample(c(rep(1, 18), rep(0, 12))),
 sample(c(rep(1, 18), rep(0, 12))))

In the R codes above, ev.ind is the phenotype (present/absent) variable with 1 = present and 0 = absent, and al.gp.ind is the genotype-group-indicator variable with 1 = AA, 2 = AB and 3 = BB. This completes step 3. In the final step, we run a logistic regression with ev.ind as the dependent binary variable and al.gp.ind as the independent variable, and then exponentiate the for al.gp.ind to get the combined odds ratio for the additive genetic model as . This is remarkably very close to the unknown combined odds ratio from the original study in our example, . The difference is less than . Using the standard error for the , we calculate/recover the confidence interval as , which is also remarkably close to the original interval . Thus our slight adaptation of Di Pietrantonj’s method works very well for this example. We did not conduct extensive simulations for the method presented in this section, since Di Pietrantonj and Veroniki and co-authors have studied elaborately the core method.

## 7 Conclusion

Meta analysis is being increasingly used to estimate the phenotype differences between genotype groups from genetic association studies. When the underlying genetic model is dominant or recessive, summary statistics from individual studies can be combined to get the pooled estimate in the meta-analysis. However, we show via simulations that when the underlying model is additive, the pooled estimate based on summary statistics leads to biased results. Since, data from individual studies is rarely available for the meta-analyst, we recommend using simulated data based on the summary statistics. We show in this paper that the method based on such simulations leads to much improved results.

## 8 Acknowledgements

Majnu John's work was supported in part by grants from the National Institute of Mental Health for an Advanced Center for Intervention and Services Research (P30 MH090590) and a Center for Intervention Development and Applied Research (P50 MH080173).

## References

• (1) Lewis CM. Genetic association studies: design, analysis and interpretation. Brief Bioinform. 2002 Jun;3(2):146-53. Review.
• (2) Zhang JP, Lencz T, Zhang RX, Nitta M, Maayan L, John M, Robinson DG, Fleischhacker WW, Kahn RS, Ophoff RA, Kane JM, Malhotra AK, Correll CU. Pharmacogenetic Associations of Antipsychotic Drug-Related Weight Gain: A Systematic Review and Meta-analysis. Schizophr Bull. 2016 Nov;42(6):1418-1437.
• (3) Borenstein M, Hedges LV, Higgins JPT, Rothstein HR. Introduction to Meta-analysis. Chichester (UK): John Wiley Sons, 2009.
• (4) Hedges LV, Shymansky KA, Woodworth G. Practical guide to modern methods of Meta-analysis. National Science Teachers Association, 1989.
• (5) Marron S, and Wand M. Exact Mean Integrated Squared Error. Annals of Statistcs. 1992 20: 712-736.
• (6) John M, Priebe CE. A data-adaptive methodology for finding an optimal weighted generalized Mann-Whitney-Wilcoxon statistic. Computational Statistics Data Analysis. 2007 51: 4337-4353
• (7) Di Pietrantonj C. Four-fold table cell frequencies imputation in meta analysis. Statistics in Medicine. 2006 25: 2299-2322.
• (8) Veroniki AA, Pavlides M, Patsopoulos NA, Salanti G. Reconstructing contingency tables from odds ratios using the Di Pietrantonj method: difficulties, constraints and impact in meta-analysis results. Res Synth Methods. 2013 Mar;4(1):78-94.

## Appendix A Tables of simulation results

 Table A1a. Original Data from standard normal density Number of Studies = 5 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.3829 0.2299 0.1995 0.0966 15 20 10 0.3510 0.2203 0.1603 0.0794 15 20 30 0.2972 0.1701 0.1272 0.0687 15 45 30 0.3706 0.2579 0.1253 0.0623 35 45 30 0.3290 0.2498 0.1055 0.0497 75 100 60 0.3066 0.2401 0.0700 0.0298 150 200 120 0.2764 0.2240 0.0472 0.0239 300 400 240 0.2902 0.2508 0.0315 0.0164 4 5.5 9 10 15 5 0.6053 0.4559 0.2209 0.1099 15 20 10 0.4947 0.4000 0.1819 0.0881 15 20 30 0.4005 0.2902 0.1394 0.0754 15 45 30 0.4482 0.3438 0.1345 0.0670 35 45 30 0.4907 0.4231 0.1156 0.0515 75 100 60 0.4722 0.4286 0.0787 0.0315 150 200 120 0.4660 0.4344 0.0529 0.0257 300 400 240 0.4899 0.4528 0.0376 0.0196 4 5.5 11 10 15 5 0.9037 0.7497 0.2483 0.1199 15 20 10 0.7837 0.6901 0.1916 0.0997 15 20 30 0.6235 0.5160 0.1568 0.0840 15 45 30 0.6453 0.5342 0.1563 0.0845 35 45 30 0.7449 0.6751 0.1282 0.0589 75 100 60 0.7291 0.6827 0.0867 0.0318 150 200 120 0.7380 0.7120 0.0564 0.0290 300 400 240 0.7850 0.7387 0.0419 0.0212 5 4 5.5 7 10 15 5 0.2085 0.0961 0.1988 0.0935 15 20 10 0.1679 0.0864 0.1610 0.0773 15 20 30 0.1417 0.0714 0.1303 0.0656 15 45 30 0.1522 0.0729 0.1327 0.0590 35 45 30 0.1026 0.0444 0.0950 0.0395 75 100 60 0.0843 0.0434 0.0681 0.0292 150 200 120 0.0576 0.0266 0.0464 0.0203 300 400 240 0.0516 0.0268 0.0320 0.0154 4 5.5 9 10 15 5 0.2224 0.1033 0.2013 0.0920 15 20 10 0.1785 0.0922 0.1621 0.0772 15 20 30 0.1490 0.0766 0.1336 0.0682 15 45 30 0.1542 0.0680 0.1326 0.0604 35 45 30 0.1142 0.0532 0.0976 0.0401 75 100 60 0.0981 0.0590 0.0699 0.0294 150 200 120 0.0667 0.0458 0.0328 0.0156 300 400 240 0.0726 0.0428 0.0479 0.0209 4 5.5 11 10 15 5 0.2444 0.1262 0.2050 0.0902 15 20 10 0.1966 0.1068 0.1632 0.0768 15 20 30 0.1571 0.0813 0.1383 0.0704 15 45 30 0.1608 0.0703 0.1335 0.0625 35 45 30 0.1326 0.0714 0.1014 0.0417 75 100 60 0.1191 0.0838 0.0721 0.0296 150 200 120 0.0967 0.0719 0.0493 0.0214 300 400 240 0.0935 0.0731 0.0337 0.0158
 Table A1b. Original Data from standard normal density Number of Studies = 10 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.4157 0.1933 0.2326 0.0764 15 20 10 0.3377 0.1787 0.1777 0.0552 15 20 30 0.3083 0.1707 0.1342 0.0474 15 45 30 0.3551 0.2128 0.1297 0.0462 35 45 30 0.3108 0.2279 0.0998 0.0338 75 100 60 0.2816 0.2125 0.0705 0.0222 150 200 120 0.2853 0.2324 0.0494 0.0145 300 400 240 0.2712 0.2183 0.0371 0.0130 4 5.5 9 10 15 5 0.6007 0.4092 0.2432 0.0787 15 20 10 0.4995 0.3650 0.1938 0.0586 15 20 30 0.3896 0.2476 0.1633 0.0583 15 45 30 0.4670 0.3326 0.1475 0.0616 35 45 30 0.4782 0.3970 0.1169 0.0420 75 100 60 0.4540 0.4019 0.0786 0.0263 150 200 120 0.4658 0.4253 0.0554 0.0148 300 400 240 0.4688 0.4372 0.0399 0.0136 4 5.5 11 10 15 5 0.9001 0.7244 0.2559 0.0768 15 20 10 0.7728 0.6490 0.1991 0.0620 15 20 30 0.5594 0.4099 0.1855 0.0729 15 45 30 0.6814 0.5550 0.1595 0.0610 35 45 30 0.7433 0.6610 0.1279 0.0456 75 100 60 0.7045 0.6577 0.0850 0.0300 150 200 120 0.7453 0.7031 0.0617 0.0176 300 400 240 0.7683 0.7344 0.0433 0.0138 5 4 5.5 7 10 15 5 0.2234 0.0685 0.2186 0.0695 15 20 10 0.1732 0.0540 0.1672 0.0482 15 20 30 0.1353 0.0487 0.1227 0.0467 15 45 30 0.1370 0.0438 0.1216 0.0405 35 45 30 0.1007 0.0330 0.0939 0.0283 75 100 60 0.0812 0.0315 0.0699 0.0235 150 200 120 0.0606 0.0194 0.0486 0.0145 300 400 240 0.0520 0.0221 0.0343 0.0107 4 5.5 9 10 15 5 0.2335 0.0722 0.2220 0.0695 15 20 10 0.1818 0.0632 0.1705 0.0492 15 20 30 0.1400 0.0507 0.1256 0.0479 15 45 30 0.1396 0.0466 0.1249 0.0434 35 45 30 0.1102 0.0410 0.0969 0.0277 75 100 60 0.0944 0.0408 0.0706 0.0239 150 200 120 0.0749 0.0381 0.0498 0.0143 300 400 240 0.0665 0.0419 0.0352 0.0108 4 5.5 11 10 15 5 0.2530 0.0935 0.2241 0.0692 15 20 10 0.1985 0.0842 0.1727 0.0504 15 20 30 0.1464 0.0555 0.1308 0.0501 15 45 30 0.1492 0.0528 0.1289 0.0461 35 45 30 0.1278 0.0591 0.0999 0.0280 75 100 60 0.1164 0.0736 0.0714 0.0244 150 200 120 0.0979 0.0667 0.0510 0.0141 300 400 240 0.0927 0.0714 0.0361 0.0111
 Table A1c. Original Data from standard normal density Number of Studies = 15 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.4404 0.1842 0.2143 0.0738 15 20 10 0.3529 0.1937 0.1619 0.0369 15 20 30 0.3204 0.1837 0.1329 0.0329 15 45 30 0.3502 0.2064 0.1380 0.0337 35 45 30 0.3003 0.2071 0.1000 0.0269 75 100 60 0.2791 0.2014 0.0718 0.0213 150 200 120 0.3066 0.2366 0.0494 0.0121 300 400 240 0.2952 0.2396 0.0368 0.0110 4 5.5 9 10 15 5 0.6261 0.3907 0.2375 0.0803 15 20 10 0.5208 0.3784 0.1799 0.0401 15 20 30 0.4280 0.2952 0.1510 0.0407 15 45 30 0.4476 0.3060 0.1558 0.0403 35 45 30 0.4666 0.3795 0.1124 0.0298 75 100 60 0.4278 0.3717 0.0797 0.0261 150 200 120 0.5009 0.4542 0.0560 0.0155 300 400 240 0.4818 0.4419 0.0422 0.0124 4 5.5 11 10 15 5 0.9440 0.7006 0.2601 0.0829 15 20 10 0.7859 0.6528 0.2014 0.0450 15 20 30 0.6324 0.4832 0.1736 0.0523 15 45 30 0.6350 0.4979 0.1701 0.0446 35 45 30 0.7311 0.6449 0.1214 0.0371 75 100 60 0.6668 0.6208 0.0873 0.0283 150 200 120 0.7919 0.7504 0.0614 0.0183 300 400 240 0.7577 0.7195 0.0460 0.0135 5 4 5.5 7 10 15 5 0.2224 0.0697 0.2110 0.0645 15 20 10 0.1698 0.0394 0.1589 0.0317 15 20 30 0.1350 0.0409 0.1244 0.0368 15 45 30 0.1439 0.0381 0.1263 0.0324 35 45 30 0.1032 0.0283 0.0960 0.0273 75 100 60 0.0805 0.0262 0.0712 0.0213 150 200 120 0.0607 0.0193 0.0485 0.0134 300 400 240 0.0506 0.0196 0.0347 0.0095 4 5.5 9 10 15 5 0.2355 0.0777 0.2137 0.0654 15 20 10 0.1778 0.0482 0.1604 0.0322 15 20 30 0.1399 0.0434 0.1286 0.0372 15 45 30 0.1483 0.0409 0.1295 0.0333 35 45 30 0.1124 0.0337 0.0984 0.0280 75 100 60 0.0917 0.0402 0.0723 0.0213 150 200 120 0.0748 0.0381 0.0495 0.0136 300 400 240 0.0659 0.0378 0.0356 0.0099 4 5.5 11 10 15 5 0.2580 0.0997 0.2151 0.0656 15 20 10 0.1939 0.0678 0.1623 0.0327 15 20 30 0.1464 0.0444 0.1336 0.0383 15 45 30 0.1577 0.0480 0.1331 0.0341 35 45 30 0.1277 0.0475 0.1010 0.0289 75 100 60 0.1107 0.0651 0.0736 0.0215 150 200 120 0.0991 0.0662 0.0505 0.0137 300 400 240 0.0907 0.0664 0.0366 0.0103
 Table A2a. Original Data from strongly right skewed density Number of Studies = 5 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.3808 0.2275 0.2511 0.1179 15 20 10 0.3704 0.2077 0.1987 0.0935 15 20 30 0.2933 0.1860 0.1419 0.0744 15 45 30 0.4627 0.3399 0.1297 0.0640 35 45 30 0.3384 0.2615 0.1136 0.0598 75 100 60 0.2580 0.1966 0.0835 0.0393 150 200 120 0.2727 0.2134 0.0556 0.0268 300 400 240 0.3476 0.2986 0.0416 0.0203 4 5.5 9 10 15 5 0.5351 0.3683 0.2822 0.1350 15 20 10 0.5413 0.3808 0.2308 0.1152 15 20 30 0.4209 0.2793 0.1828 0.0914 15 45 30 0.5838 0.4453 0.1520 0.0755 35 45 30 0.4889 0.4182 0.1377 0.0746 75 100 60 0.4375 0.3800 0.0954 0.0446 150 200 120 0.4442 0.4113 0.0629 0.0279 300 400 240 0.5385 0.4977 0.0515 0.0267 4 5.5 11 10 15 5 0.7905 0.6018 0.3165 0.1513 15 20 10 0.8277 0.6335 0.2590 0.1273 15 20 30 0.6337 0.4587 0.2109 0.1026 15 45 30 0.7701 0.6212 0.1811 0.1001 35 45 30 0.7286 0.6331 0.1609 0.0867 75 100 60 0.7081 0.6646 0.1038 0.0490 150 200 120 0.7416 0.7093 0.0722 0.0328 300 400 240 0.8530 0.8260 0.0550 0.0298 5 4 5.5 7 10 15 5 0.2323 0.1054 0.2224 0.0915 15 20 10 0.1687 0.0760 0.1632 0.0774 15 20 30 0.1344 0.0666 0.1236 0.0622 15 45 30 0.1443 0.0655 0.1279 0.0534 35 45 30 0.1138 0.0527 0.1058 0.0486 75 100 60 0.0888 0.0396 0.0778 0.0340 150 200 120 0.0676 0.0391 0.0517 0.0240 300 400 240 0.0532 0.0311 0.0347 0.0157 4 5.5 9 10 15 5 0.2498 0.1139 0.2300 0.0957 15 20 10 0.1822 0.0842 0.1724 0.0824 15 20 30 0.1416 0.0746 0.1326 0.0684 15 45 30 0.1528 0.0655 0.1324 0.0560 35 45 30 0.1251 0.0616 0.1104 0.0537 75 100 60 0.1021 0.0510 0.0803 0.0344 150 200 120 0.0835 0.0591 0.0533 0.0242 300 400 240 0.0721 0.0510 0.0372 0.0172 4 5.5 11 10 15 5 0.2725 0.1320 0.2397 0.1021 15 20 10 0.2000 0.0995 0.1802 0.0854 15 20 30 0.1519 0.0813 0.1428 0.0732 15 45 30 0.1613 0.0724 0.1392 0.0589 35 45 30 0.1433 0.0741 0.1153 0.0587 75 100 60 0.1232 0.0742 0.0834 0.0349 150 200 120 0.1085 0.0892 0.0557 0.0248 300 400 240 0.1006 0.0808 0.0396 0.0187
 Table A2b. Original Data from strongly right skewed density Number of Studies = 10 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.4718 0.2340 0.2556 0.0832 15 20 10 0.3490 0.2008 0.1856 0.0631 15 20 30 0.3304 0.1959 0.1436 0.0450 15 45 30 0.4147 0.2617 0.1481 0.0520 35 45 30 0.3467 0.2429 0.1172 0.0432 75 100 60 0.2840 0.2164 0.0859 0.0274 150 200 120 0.2817 0.1962 0.0555 0.0183 300 400 240 0.2703 0.2227 0.0388 0.0139 4 5.5 9 10 15 5 0.6651 0.4418 0.2747 0.0926 15 20 10 0.5253 0.3661 0.2243 0.0765 15 20 30 0.4522 0.2946 0.1760 0.0545 15 45 30 0.5396 0.3827 0.1763 0.0689 35 45 30 0.5122 0.4274 0.1362 0.0495 75 100 60 0.4682 0.4082 0.0970 0.0286 150 200 120 0.4416 0.3939 0.0620 0.0210 300 400 240 0.4437 0.4107 0.0449 0.0168 4 5.5 11 10 15 5 0.9679 0.7447 0.2990 0.1016 15 20 10 0.7988 0.6482 0.2401 0.0855 15 20 30 0.6595 0.4927 0.1987 0.0601 15 45 30 0.7407 0.5690 0.2088 0.0803 35 45 30 0.7808 0.7001 0.1531 0.0512 75 100 60 0.7543 0.6962 0.1040 0.0308 150 200 120 0.7155 0.6840 0.0668 0.0261 300 400 240 0.7349 0.7039 0.0495 0.0186 5 4 5.5 7 10 15 5 0.2418 0.0723 0.2241 0.0684 15 20 10 0.1705 0.0549 0.1639 0.0511 15 20 30 0.1395 0.0458 0.1283 0.0407 15 45 30 0.1456 0.0473 0.1214 0.0426 35 45 30 0.1094 0.0431 0.1010 0.0358 75 100 60 0.0836 0.0311 0.0720 0.0212 150 200 120 0.0616 0.0253 0.0485 0.0163 300 400 240 0.0482 0.0225 0.0337 0.0117 4 5.5 9 10 15 5 0.2530 0.0812 0.2277 0.0693 15 20 10 0.1818 0.0631 0.1701 0.0539 15 20 30 0.1481 0.0448 0.1338 0.0423 15 45 30 0.1513 0.0492 0.1273 0.0453 35 45 30 0.1198 0.0535 0.1028 0.0362 75 100 60 0.0969 0.0468 0.0744 0.0217 150 200 120 0.0742 0.0419 0.0498 0.0177 300 400 240 0.0628 0.0404 0.0352 0.0122 4 5.5 11 10 15 5 0.2736 0.1043 0.2323 0.0700 15 20 10 0.2003 0.0754 0.1760 0.0558 15 20 30 0.1593 0.0452 0.1411 0.0439 15 45 30 0.1610 0.0569 0.1348 0.0490 35 45 30 0.1395 0.0727 0.1056 0.0357 75 100 60 0.1194 0.0744 0.0771 0.0224 150 200 120 0.0967 0.0700 0.0514 0.0190 300 400 240 0.0861 0.0679 0.0368 0.0125
 Table A2c. Original Data from strongly right skewed density Number of Studies = 15 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.4125 0.2003 0.2456 0.0814 15 20 10 0.3673 0.1902 0.1871 0.0441 15 20 30 0.3131 0.1634 0.1491 0.0477 15 45 30 0.4248 0.2665 0.1403 0.0443 35 45 30 0.2919 0.1947 0.1167 0.0340 75 100 60 0.3022 0.2353 0.0713 0.0243 150 200 120 0.2887 0.2166 0.0556 0.0132 300 400 240 0.3033 0.2501 0.0362 0.0126 4 5.5 9 10 15 5 0.6115 0.4281 0.2695 0.0799 15 20 10 0.5609 0.3867 0.2168 0.0608 15 20 30 0.4175 0.2576 0.1756 0.0563 15 45 30 0.5572 0.4020 0.1652 0.0509 35 45 30 0.4557 0.3743 0.1445 0.0393 75 100 60 0.4956 0.4433 0.0816 0.0283 150 200 120 0.4721 0.4225 0.0649 0.0178 300 400 240 0.4903 0.4553 0.0431 0.0159 4 5.5 11 10 15 5 0.9083 0.7306 0.2865 0.0821 15 20 10 0.8513 0.6581 0.2434 0.0763 15 20 30 0.6007 0.4358 0.2027 0.0642 15 45 30 0.7771 0.6164 0.1854 0.0520 35 45 30 0.7512 0.6585 0.1601 0.0440 75 100 60 0.7873 0.7333 0.0896 0.0321 150 200 120 0.7584 0.7153 0.0726 0.0239 300 400 240 0.7852 0.7526 0.0491 0.0176 5 4 5.5 7 10 15 5 0.2427 0.0703 0.2258 0.0631 15 20 10 0.1779 0.0474 0.1732 0.0432 15 20 30 0.1389 0.0390 0.1320 0.0383 15 45 30 0.1490 0.0362 0.1260 0.0315 35 45 30 0.1097 0.0348 0.1044 0.0289 75 100 60 0.0758 0.0268 0.0668 0.0189 150 200 120 0.0604 0.0182 0.0480 0.0110 300 400 240 0.0521 0.0225 0.0332 0.0104 4 5.5 9 10 15 5 0.2603 0.0918 0.2316 0.0636 15 20 10 0.1905 0.0557 0.1798 0.0460 15 20 30 0.1471 0.0424 0.1386 0.0396 15 45 30 0.1533 0.0404 0.1320 0.0348 35 45 30 0.1223 0.0476 0.1106 0.0309 75 100 60 0.0885 0.0450 0.0689 0.0199 150 200 120 0.0744 0.0382 0.0499 0.0121 300 400 240 0.0697 0.0443 0.0348 0.0113 4 5.5 11 10 15 5 0.2867 0.1237 0.2392 0.0651 15 20 10 0.2109 0.0728 0.1878 0.0492 15 20 30 0.1588 0.0460 0.1467 0.0411 15 45 30 0.1612 0.0479 0.1393 0.0382 35 45 30 0.1423 0.0673 0.1170 0.0322 75 100 60 0.1116 0.0722 0.0713 0.0210 150 200 120 0.0977 0.0677 0.0522 0.0137 300 400 240 0.0967 0.0748 0.0365 0.0121
 Table A3a. Original Data from asymmetric bimodal density Number of Studies = 5 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.4474 0.2842 0.2214 0.1132 15 20 10 0.3959 0.2646 0.1746 0.0903 15 20 30 0.3257 0.2465 0.1462 0.0800 15 45 30 0.4514 0.3609 0.1366 0.0744 35 45 30 0.3268 0.2500 0.1193 0.0661 75 100 60 0.3091 0.2501 0.0866 0.0514 150 200 120 0.2938 0.2506 0.0731 0.0511 300 400 240 0.3109 0.2684 0.0694 0.0549 4 5.5 9 10 15 5 0.6536 0.5111 0.2486 0.1185 15 20 10 0.5830 0.4879 0.1939 0.0983 15 20 30 0.4565 0.3746 0.1738 0.1017 15 45 30 0.6238 0.5616 0.1634 0.0938 35 45 30 0.5173 0.4797 0.1441 0.0951 75 100 60 0.5008 0.4781 0.1010 0.0768 150 200 120 0.4750 0.4525 0.0944 0.0796 300 400 240 0.4774 0.4572 0.0953 0.0875 4 5.5 11 10 15 5 0.9596 0.8128 0.2660 0.1255 15 20 10 0.8741 0.7817 0.2143 0.1132 15 20 30 0.7028 0.6035 0.1945 0.1149 15 45 30 0.9031 0.8340 0.1903 0.1163 35 45 30 0.8144 0.7876 0.1634 0.1228 75 100 60 0.7995 0.7845 0.1186 0.0952 150 200 120 0.7363 0.7202 0.1108 0.1001 300 400 240 0.7428 0.7342 0.1136 0.1086 5 4 5.5 7 10 15 5 0.2297 0.0972 0.2157 0.0952 15 20 10 0.1742 0.0882 0.1671 0.0857 15 20 30 0.1329 0.0660 0.1292 0.0640 15 45 30 0.1582 0.0877 0.1302 0.0712 35 45 30 0.1113 0.0571 0.1005 0.0525 75 100 60 0.0861 0.0450 0.0747 0.0372 150 200 120 0.0649 0.0411 0.0507 0.0287 300 400 240 0.0598 0.0441 0.0442 0.0288 4 5.5 9 10 15 5 0.2458 0.1125 0.2209 0.0975 15 20 10 0.1848 0.1027 0.1676 0.0879 15 20 30 0.1368 0.0664 0.1339 0.0662 15 45 30 0.1602 0.0865 0.1356 0.0751 35 45 30 0.1270 0.0800 0.1053 0.0610 75 100 60 0.1050 0.0767 0.0792 0.0456 150 200 120 0.0876 0.0711 0.0588 0.0402 300 400 240 0.0862 0.0785 0.0532 0.0427 4 5.5 11 10 15 5 0.2720 0.1443 0.2248 0.0993 15 20 10 0.2094 0.1302 0.1686 0.0918 15 20 30 0.1450 0.0711 0.1401 0.0712 15 45 30 0.1647 0.0860 0.1422 0.0795 35 45 30 0.1547 0.1193 0.1119 0.0708 75 100 60 0.1347 0.1167 0.0841 0.0557 150 200 120 0.1192 0.1089 0.0673 0.0524 300 400 240 0.1237 0.1200 0.0631 0.0562
 Table A3b. Original Data from asymmetric bimodal density Number of Studies = 10 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.4523 0.2457 0.2135 0.0889 15 20 10 0.3894 0.2341 0.1819 0.0701 15 20 30 0.3141 0.2013 0.1441 0.0752 15 45 30 0.3725 0.2546 0.1363 0.0581 35 45 30 0.3137 0.2492 0.1181 0.0602 75 100 60 0.3209 0.2641 0.0887 0.0524 150 200 120 0.3472 0.2926 0.0761 0.0480 300 400 240 0.3588 0.3041 0.0625 0.0461 4 5.5 9 10 15 5 0.6790 0.5392 0.2332 0.0945 15 20 10 0.6293 0.5597 0.1983 0.1115 15 20 30 0.4620 0.3713 0.1701 0.0968 15 45 30 0.4931 0.4135 0.1536 0.0906 35 45 30 0.4798 0.4423 0.1391 0.0873 75 100 60 0.5533 0.5240 0.1074 0.0793 150 200 120 0.5766 0.5568 0.1055 0.0946 300 400 240 0.5523 0.5378 0.0906 0.0840 4 5.5 11 10 15 5 0.9939 0.8866 0.2901 0.1455 15 20 10 0.8933 0.8201 0.2084 0.1137 15 20 30 0.6706 0.5891 0.2006 0.1131 15 45 30 0.7514 0.6904 0.1825 0.1214 35 45 30 0.7834 0.7431 0.1494 0.0990 75 100 60 0.8285 0.8051 0.1326 0.1098 150 200 120 0.8009 0.7893 0.1221 0.1127 300 400 240 0.7994 0.7925 0.1072 0.1034 5 4 5.5 7 10 15 5 0.2262 0.0786 0.2213 0.0727 15 20 10 0.1741 0.0648 0.1725 0.0567 15 20 30 0.1405 0.0463 0.1325 0.0427 15 45 30 0.1364 0.0457 0.1197 0.0417 35 45 30 0.1106 0.0456 0.1061 0.0375 75 100 60 0.0866 0.0451 0.0746 0.0316 150 200 120 0.0727 0.0435 0.0545 0.0261 300 400 240 0.0570 0.0346 0.0398 0.0214 4 5.5 9 10 15 5 0.2435 0.0941 0.2192 0.0693 15 20 10 0.1846 0.0810 0.1724 0.0585 15 20 30 0.1517 0.0502 0.1434 0.0529 15 45 30 0.1505 0.0584 0.1330 0.0559 35 45 30 0.1244 0.0645 0.1075 0.0377 75 100 60 0.1083 0.0785 0.0764 0.0408 150 200 120 0.0928 0.0733 0.0609 0.0369 300 400 240 0.0904 0.0782 0.0534 0.0397 4 5.5 11 10 15 5 0.2924 0.1567 0.2362 0.0953 15 20 10 0.1970 0.1134 0.1654 0.0608 15 20 30 0.1398 0.0558 0.1345 0.0660 15 45 30 0.1454 0.0481 0.1298 0.0607 35 45 30 0.1495 0.1042 0.1143 0.0540 75 100 60 0.1403 0.1176 0.0873 0.0526 150 200 120 0.1205 0.1085 0.0660 0.0501 300 400 240 0.1305 0.1271 0.0619 0.0541
 Table A3c. Original Data from asymmetric bimodal density Number of Studies = 15 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.4665 0.2749 0.2312 0.0802 15 20 10 0.3695 0.2537 0.1786 0.0683 15 20 30 0.3311 0.2384 0.1471 0.0652 15 45 30 0.4150 0.3044 0.1397 0.0659 35 45 30 0.3250 0.2479 0.1128 0.0479 75 100 60 0.3188 0.2597 0.0931 0.0529 150 200 120 0.3334 0.2703 0.0786 0.0519 300 400 240 0.3371 0.2937 0.0655 0.0511 4 5.5 9 10 15 5 0.6454 0.5243 0.2491 0.0996 15 20 10 0.6392 0.5347 0.2050 0.0894 15 20 30 0.4335 0.3251 0.1659 0.0753 15 45 30 0.4990 0.4051 0.1632 0.0882 35 45 30 0.5162 0.4633 0.1374 0.0770 75 100 60 0.5730 0.5502 0.1028 0.0743 150 200 120 0.5270 0.5108 0.0907 0.0776 300 400 240 0.5399 0.5279 0.0882 0.0825 4 5.5 11 10 15 5 0.9689 0.8433 0.2737 0.0884 15 20 10 0.8685 0.7752 0.2066 0.0927 15 20 30 0.6852 0.5711 0.1976 0.1001 15 45 30 0.7605 0.6873 0.1766 0.1064 35 45 30 0.8375 0.7917 0.1507 0.0964 75 100 60 0.8397 0.8227 0.1215 0.0986 150 200 120 0.9333 0.8284 0.1132 0.1017 300 400 240 0.8358 0.8286 0.0993 0.0948 5 4 5.5 7 10 15 5 0.2233 0.0789 0.2126 0.0679 15 20 10 0.1793 0.0642 0.1698 0.0556 15 20 30 0.1320 0.0369 0.1242 0.0374 15 45 30 0.1399 0.0564 0.1248 0.0429 35 45 30 0.1032 0.0358 0.0977 0.0299 75 100 60 0.0841 0.0391 0.0747 0.0265 150 200 120 0.0712 0.0432 0.0566 0.0277 300 400 240 0.0620 0.0381 0.0432 0.0222 4 5.5 9 10 15 5 0.2547 0.1045 0.2292 0.0663 15 20 10 0.1826 0.0826 0.1655 0.0562 15 20 30 0.1376 0.0466 0.1304 0.0484 15 45 30 0.1502 0.0478 0.1320 0.0495 35 45 30 0.1213 0.0678 0.1074 0.0432 75 100 60 0.1103 0.0745 0.0795 0.0368 150 200 120 0.0952 0.0761 0.0630 0.0372 300 400 240 0.0829 0.0721 0.0479 0.0359 4 5.5 11 10 15 5 0.2713 0.1487 0.2230 0.0762 15 20 10 0.2018 0.1072 0.1702 0.0547 15 20 30 0.1419 0.0502 0.1403 0.0650 15 45 30 0.1453 0.0447 0.1300 0.0574 35 45 30 0.1575 0.1141 0.1153 0.0562 75 100 60 0.1440 0.1258 0.0828 0.0550 150 200 120 0.1205 0.1083 0.0637 0.0428 300 400 240 0.1207 0.1153 0.0602 0.0526
 Table A4a. Original Data from heavily kurtotic density Number of Studies = 5 Crude Approach Simulations Approach σWS M1 M2 M3 N1 N2 N3 Mean absolute Mean absolute Mean absolute Mean absolute bias for g bias for g-WM bias for g bias for g-WM 1 4 5.5 7 10 15 5 0.4617 0.2848 0.2477 0.1173 15 20 10 0.3794 0.2681 0.1833 0.0823 15 20 30 0.3278 0.1952 0.1527 0.0698 15 45 30 0.3809 0.2349 0.1434 0.0559 35 45 30 0.3478 0.2534 0.1208 0.0598 75 100 60 0.3805 0.3192 0.0825 0.0404 150 200 120 0.3348 0.2887 0.0536 0.0259 300 400 240 0.2578 0.1993 0.0388 0.0160 4 5.5 9 10 15 5 0.7155 0.5305 0.2708 0.1407 15 20 10 0.5752 0.4717 0.2082 0.0961 15 20 30 0.4638 0.2877 0.1840 0.0862 15 45 30 0.5120 0.3831 0.1587 0.0684 35 45 30 0.5492 0.4414 0.1414 0.0712 75 100 60 0.6003 0.5469 0.0996 0.0540 150 200 120 0.5332 0.4889 0.0643 0.0322 300 400 240 0.3970 0.3512 0.0465 0.0186 4 5.5 11 10 15 5 1.0573 0.8815 0.2714 0.1471 15 20 10 0.8831 0.7722 0.2264 0.0996 15 20 30 0.6700 0.4892 0.2044 0.1016 15 45 30 0.7480 0.6143 0.1681 0.0753 35 45 30 0.8052 0.7034 0.1511 0.0784