The fraternal birth-order effect as a statistical artefact: convergent evidence from probability calculus, simulated data, and multiverse meta-analysis

Notation

The most widely used effect size for quantifying the number of older brothers among homosexual men in comparison to heterosexual men, the older brothers odds ratio (OBOR), is defined at the level of siblings, as reported by study participants. Hence, much of the observational evidence for the FBOE (including most of the meta-analyses) rests on a rather unusual approach for comparing the number of older brothers between groups. Intuitively, birth-order effects are analyzed by considering the sampled individuals as the units of analysis (e.g., Rohrer, Egloff & Schmukle, 2015). For instance, to quantify the association between an individual’s sexuality and their number of older brothers, one could use either of these variables as the outcome and the other as one of the predictors in a (generalized) linear model. However, the OBOR (but also other recommended statistical approaches) does not build on this intuition and treats the reported siblings as the units of analysis. We thus need to introduce some notation to be able to show how precisely the OBOR fails to address the effect of interest.

If the reported siblings are regarded as the sample, the following events can be defined: A given sibling can be either an older brother, an older sister, a younger brother, or a younger sister of the study participant who reported him or her. These events are denoted by OB, OS, YB and YS, respectively. The complement of any event is denoted by the superscript c (as in “complement”). For instance, OB^c denotes the event that a given sibling is either an older sister, a younger brother, or a younger sister. The event OB^c will mostly be referred to as Other, meaning that the sibling is not an older brother, but any of the other three possible sibling types. There are instances in this work, where Other refers to siblings who are not older sisters, which should then be clear from the given context. Furthermore, we denote the event that a given sibling is an older sibling (i.e., the union of OB and OS) by Older. Its complement Older^c denotes the event that a given sibling is a younger sibling (for the sake of simplicity, twins, triplets, etc., are not accounted for).

A sibling can either be reported by a homosexual or by a heterosexual study participant. These events are denoted by Hom and Het, respectively, and are regarded as complementary. To aid readability, we use Het instead of Hom^c, as the events Hom and Het are extensively referred to in the following.

The probability that an event A occurs, is denoted by P(A). The probability of an intersection of two events A and B is denoted by P(A, B). The probability of an event A conditional on an event B (i.e., the conditional probability of A given B) is denoted by P(A|B).

In order to distinguish the events which pertain to individual siblings from the number of times these events occur in a sample, we use the number sign (#), which in this context should be read as “number of.” Hence, #OB simply denotes the number of all reported older brothers in a sample. The subscripts Het and Hom are used to denote whether the total number of a given sibling type refers to that reported by homosexual or heterosexual participants. For instance, the number of older brothers reported by homosexual participants is denoted by #OB_Hom.

We treat the terms “probability of an event” and “proportion of times an event occurred” interchangeably, with the latter being an estimate of the former.

Blanchard (2014, 2018a) and Blanchard et al. (2020, 2021) warned that the specific (average) difference in the number of older brothers between homosexual and heterosexual individuals due to the FBOE may go undetected (or be obscured), if the mean number of all siblings reported by the group of homosexual men is appreciably smaller than the mean number of all siblings reported by the group of heterosexual men.³

To counter potential mitigation of the relationship between homosexual orientation and the number of older brothers due to the suspected confounding variable “differences in family size between groups”, Blanchard suggested that it is necessary to adjust for either the number of other siblings (#Other; Blanchard, 2014, 2018a, 2020) or the number of all siblings (#All; Blanchard, 2014) in a statistical model (i.e., test) comparing these two groups.⁴

Over the years, three new methods, designed to achieve this goal, have been developed, and increasingly are used in primary publications and meta-analyses on the FBOE. These are the already mentioned OBOR (Apostolou, 2020; Blanchard, 2018a, 2018b; Skorska & Bogaert, 2020), the modified ratio of older brothers (MROB; Blanchard, 2014), and the modified proportion of older brothers (MPOB; Blanchard, 2014; e.g., Skorska & Bogaert, 2020).

For now, we show that the OBOR, MROB, and MPOB quantify the more general average difference (or association) in the number of older siblings (i.e., brothers and sisters) between homosexual and heterosexual men, rather than the more specific difference in the number of older brothers only (or the association between the number of brothers and sexual orientation, only).

The OBOR

As the name implies, the OBOR (see Table 2) is an odds ratio, and it can thus be expressed in terms of conditional probabilities:

(1) $$\begin{array}{c}\begin{array}{c}OBOR={\displaystyle \frac{\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{o}\mathrm{m})/[1-\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{o}\mathrm{m})]}{\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{e}\mathrm{t})/[1-\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{e}\mathrm{t})]}.}\end{array}\end{array}$$

Now suppose we observe an OBOR > 1, which, according to the interpretation of Blanchard (2018a, 2018b, 2018c) and Blanchard et al. (2020, 2021), should be regarded as evidence for the FBOE.

For the OBOR to become greater than unity, the numerator in Eq. (1) must be greater than the denominator. Formally, this relationship can be stated as

(2) $$\begin{array}{c}\begin{array}{ccc}OBOR& >1\iff \mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{o}\mathrm{m})& >\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{e}\mathrm{t}),\end{array}\end{array}$$

Equation (2) is just an algebraic transformation of Eq. (1) (see Supplement S3 for a derivation).

By conditioning on the event Older, the law of total probability (e.g., Blitzstein & Hwang, 2019) allows for the factorisation of both P(OB|Hom) and P(OB|Het), as follows:

$$\begin{array}{cc}\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{o}\mathrm{m})& =\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r},\mathrm{H}\mathrm{o}\mathrm{m})\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}|\mathrm{H}\mathrm{o}\mathrm{m})\\ & +\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}{\mathrm{r}}^{\mathrm{c}},\mathrm{H}\mathrm{o}\mathrm{m})\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}{\mathrm{r}}^{\mathrm{c}}|\mathrm{H}\mathrm{o}\mathrm{m})\mathrm{a}\mathrm{n}\mathrm{d},\\ \mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{e}\mathrm{t})& =\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r},\mathrm{H}\mathrm{e}\mathrm{t})\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}|\mathrm{H}\mathrm{e}\mathrm{t})\\ & +\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}{\mathrm{r}}^{\mathrm{c}},\mathrm{H}\mathrm{e}\mathrm{t})\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}{\mathrm{r}}^{\mathrm{c}}|\mathrm{H}\mathrm{e}\mathrm{t}).\end{array}$$

It is impossible for a sibling to be both an older brother and a younger sibling (i.e., the event Older^c and OB is impossible). Therefore, P(OB|Older^c,Hom) and P(OB|Older^c,Het) are both zero, reducing the above factorisations to

(3) $$\begin{array}{rl}\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{o}\mathrm{m})=& \mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r},\mathrm{H}\mathrm{o}\mathrm{m})\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}|\mathrm{H}\mathrm{o}\mathrm{m}),\mathrm{a}\mathrm{n}\mathrm{d}\\ \mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{e}\mathrm{t})=& \mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r},\mathrm{H}\mathrm{e}\mathrm{t})\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}|\mathrm{H}\mathrm{e}\mathrm{t}).\end{array}$$

One can think of P(OB|Older,Hom) and P(OB|Older,Het) as first restricting the sampling space from the set of all reported siblings to the subset which contains only older siblings and then computing the proportion (i.e., probability) of brothers within this subset of older siblings for the homosexual and heterosexual group, respectively.

In combining Eqs. (2 and 3), we obtain the following equivalent inequalities:

(4) $$\begin{array}{c}\hfill \begin{array}{cc}OBOR& >1\\ \iff \mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{o}\mathrm{m})& >\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{H}\mathrm{e}\mathrm{t})\\ \iff \mathrm{P}(\mathrm{O}\mathrm{B}\left|\mathrm{H}\mathrm{o}\mathrm{m},\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r})\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}\right|\mathrm{H}\mathrm{o}\mathrm{m})& >\mathrm{P}(\mathrm{O}\mathrm{B}\left|\mathrm{H}\mathrm{e}\mathrm{t},\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r})\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}\right|\mathrm{H}\mathrm{e}\mathrm{t}).\end{array}\end{array}$$

Notice that Eq. (4) (i.e., an OBOR larger than unity) can hold even if P(OB|Hom,Older) ≤ P(OB|Het,Older). That is, one may observe that the proportion (or probability) of older brothers in the subset of older siblings is smaller in the homosexual group, as opposed to the heterosexual group. However, the OBOR (considering the entire set of reported siblings) may still be greater than unity, due to the proportion (or probability) of older siblings being sufficiently greater in the homosexual group (i.e., P(Older|Hom) > P(Older|Het)).

This contradicts the key necessary observation of the FBOE, as the homosexual group must show a greater proportion (or probability) of older brothers among older siblings. This observation is required by the claim that older brothers increase the odds of homosexual orientation in men, but older sisters do not (or to a lesser extent). The OBOR thus insufficiently accounts for a relevant confounder (here: Older) and could potentially lead to the conclusion that homosexual men have more older brothers relative to older sisters, when in fact the opposite is true. In other words, the OBOR is ambiguous: It is affected by both, the proportions of older brothers among older siblings as well as the proportions of older siblings among all siblings.

It follows that–contrary to previous claims (Blanchard, 2018a, 2020; Blanchard et al., 2021)–for an OBOR of one to be considered to adequately represent the null hypothesis of no group difference in the number of older brothers relative to older sisters, the proportion (or probability) of observing an older sibling must be identical in both groups (this is equivalent to the requirement that the number of older siblings be equal in both groups).

To see this, suppose that P(OB|Older,Het) and P(OB|Older,Hom) are equal to 0.515, which is a widely agreed upon population estimate of the proportion of male births, and hence for the probability of being born male (e.g., Grech & Mamo, 2020). Further assume that homosexual men have more older brothers and sisters compared to heterosexual men; that is, there is no specific difference in the number older brothers that exceeds the difference in the number of older siblings. This can be achieved by requiring that the probabilities of observing an older sibling in the homosexual and heterosexual groups are equal to the median proportions of older siblings of P(Older|Hom) = 0.58 and P(Older|Het) = 0.48 reported across the 45 samples included in the three meta-analyses using the OBOR (Blanchard, 2018a, 2018b; Blanchard et al., 2021). The OBOR under the null hypothesis is then given by

(5) $$\begin{array}{c}\hfill OBO{R}_{{\mathrm{H}}_0}={\displaystyle \frac{\left(0.515\times 0.58\right)/\left(1-0.515\times 0.58\right)}{\left(0.515\times .48\right)/\left(1-0.515\times 0.48\right)}=\mathrm{1.30.}}\end{array}$$

Equation (5) states a more adequate value for a null hypothesis, which would be expected if the probability of observing an older brother among the older siblings of homosexual and heterosexual men were equal to the population probability of being born male (0.515), with the additional assumption that the probability of observing an older sibling (brother or sister) is greater (0.58) in the homosexual as compared to the heterosexual (0.48) group. That is, this null model assumes that the two groups differ only with respect to the proportion of older siblings, but that there is no specific (or additional) difference in the proportion of older brothers in the homosexual group (as the FBOE requires).

The fourth meta-analysis (Blanchard, 2018a) reported an estimated random-effects mean OBOR of 1.47, 95% CI [1.33–1.62], combining a total of 30 (31; see Part III) samples. For a distinct subset of 18 samples denoted “non-feminine/cisgender” men, Blanchard (2018a) reported a mean OBOR of 1.27, 95% CI [1.20–1.35]. Both of these estimates can be regarded as more or less compatible with the value of 1.30 derived in Eq. (5). In a commentary to Blanchard (2018a), Zietsch (2018) cautioned about the evidence in favor of the FBOE contained in this meta-analysis, mainly due to concerns about questionable and arbitrary inclusion criteria used by Blanchard (2018a), which led to the exclusion of all probability samples known at that time.⁵

Blanchard (2018b) replied by putting forth new evidence for the specific difference in the number of older brothers between homosexual and heterosexual men in a fifth meta-analysis, this time amalgamating probability samples only. The estimated mean OBOR over six such probability samples was 1.21, 95% CI [1.13–1.30], and was thus judged to be incompatible with the null hypothesis of OBOR = 1.

However, in light of the above Eq. (5), this estimated OBOR appears compatible with a model of an unspecific group difference in older brothers and sisters. That is, the observed OBOR > 1 may not reflect the claimed specific difference in the number of older brothers. It goes without saying that a more general phenomenon like a larger number of older siblings in homosexual men as opposed to heterosexual men might be seen as a theoretically interesting observation on its own. The claims surrounding the FBOE, however, are highly specific and the necessary observational consequence (conditional on the assumptions underlying the FBOE) of a surplus of older brothers relative to older sisters cannot be quantified using the OBOR.

Addressing group differences with respect to the number of older sisters

It appears that Blanchard and colleagues were aware of the OBOR’s inability to distinguish between a specific difference in the number of older brothers and a more general difference in the number older siblings between homosexual and heterosexual men, as shortly after the fifth meta-analysis, the older sisters odds ratio (OSOR; see Table 2) was introduced (Blanchard, 2018c). Analogous to the OBOR, the OSOR is simply the odds ratio of older sisters among all the siblings reported by homosexual and heterosexual men. Blanchard (2018c) reported that for 29 out of the 36 samples included in the fourth and fifth meta-analyses, the OBOR was numerically greater than the OSOR, concluding that the group difference with respect to the number of older brothers is greater than the difference with respect to the number of older sisters as predicted by the FBOE. However, this result reveals nothing about the size of the difference in the OBOR and the OSOR and might have led to a different conclusion if these differences had been weighted and combined in a meta-analysis. In similar vein, the seventh meta-analysis (Blanchard et al., 2021) considered both the OBOR and the OSOR in a set of 24 samples and estimated a random-effects mean OBOR of 1.28, 95% CI [1.22–1.35] for the entire set (which again is an estimate very compatible with the value derived in Eq. (5)). In addition, Blanchard et al. (2021) reported estimated random-effects mean OSORs of 1.11, 95% CI [1.05–1.17], for what they denoted the “teleiophiles” subgroup and 1.15, 95% CI [0.99–1.34], for the “pedophiles and hebephiles” subgroup. Blanchard et al. (2021) describe pedophiles as men predominantly attracted to children before puberty, hebephiles as men predominantly attracted to children in puberty, and teleiophiles as men attracted to postpubertal, mature individuals.

The implication of the above seems to be that, in order to ensure that an observed group difference in the number of older brothers via the OBOR is not an artefact of a more general group difference in the number of older siblings (i.e., older brothers and sisters), one should interpret the OBOR and the OSOR in tandem.

However, having to interpret two separate, but correlated, statistics introduces an unnecessary degree of complexity. Moreover, the apparent need to report both the OBOR and the OSOR, as conveyed by Blanchard et al. (2021), illustrates that it is the number (or equivalently, the proportion) of older siblings that needs to be adjusted for–not total family size or number of “other” siblings, as has repeatedly been asserted (Blanchard, 2014, 2018a, 2018b).

If researchers want to use an odds ratio in order to quantify the difference in the number of older brothers relative to older sisters between homosexual and heterosexual men, then the OBOR can be fixed by simply omitting the younger siblings from the denominator of the OBOR. Then, only the relevant subset of older siblings is considered. The odds of observing an older brother in either group are now given by the ratio of older brothers to older sisters, and consequently the odds ratio (OR) is given by

(6) $$\begin{array}{c}\hfill OR={\displaystyle \frac{\mathrm{\#}\mathrm{O}{\mathrm{B}}_{\mathrm{H}\mathrm{o}\mathrm{m}}/\mathrm{\#}\mathrm{O}{\mathrm{S}}_{\mathrm{H}\mathrm{o}\mathrm{m}}}{\mathrm{\#}\mathrm{O}{\mathrm{B}}_{\mathrm{H}\mathrm{e}\mathrm{t}}/\mathrm{\#}\mathrm{O}{\mathrm{S}}_{\mathrm{H}\mathrm{e}\mathrm{t}}}.}\end{array}$$

This OR adjust for the confounding effect of any differences in the number of older siblings in either group. Alternatively, one could also use the difference in the proportion of older brothers among older siblings as an effect size (see supplemental materials to Part III, below).

The numerator and denominator of Eq. (6) may be conceived as the male-to-female sex ratio within the subset of older siblings among homosexual and heterosexual participants, respectively. Under the assumption of the FBOE and the implied greater number of older brothers among older siblings of homosexual men, one would expect this ratio of sex ratios (i.e., an odds ratio) to be greater than unity.

We note again that the use of ORs (the OBOR and the OR in Eq. (6)) implies a shift from the level of participants to the siblings reported by these participants. For instance, suppose participant i reports to have three siblings, one older brother and two older sisters. Use of the OBOR or OR implies that the number of older brothers is compared to the number of other sibling types. In the example just given, each of the three siblings provides one observation of the binary (or Bernoulli) variable “older brother” (1 = sibling is older brother, 0 = sibling is not older brother). It follows that in this example, a single participant provides multiple (three) data points and thus naturally the units of analysis (i.e., siblings) are nested, or grouped, within participants.

The modified ratio and proportion of older brothers

The MROB (see Table 2 for a definition) can be regarded as an estimator of the odds of observing an older brother among all the siblings of the $i$th participant, whereas the MPOB (see Table 2 for a definition) can be regarded as an estimator of the probability of this event. Hence, given the one-to-one relation between odds and probabilities, the MROB and the MPOB are related as follows:⁶

$$MROB\approx {\displaystyle \frac{MPOB}{1-MPO{B}^{\mathrm{\prime }}}\mathrm{a}\mathrm{n}\mathrm{d}MPOB\approx {\displaystyle \frac{MROB}{1+MROB}.}}$$

Keeping in mind that we use proportions and probabilities interchangeably, the MPOB of the $i$th participant is defined as the probability of observing an older brother within the sibship of the $i$th participant, which can be written as

(7) $$\begin{array}{c}\hfill MPO{B}_i\approx \mathrm{P}(\mathrm{O}\mathrm{B}|i)=\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r},i)\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}|i).\end{array}$$

Applying once more the law of total probability, it follows that the MROB for the ith participant is given by:

(8) $$\begin{array}{c}\hfill MRO{B}_i\approx {\displaystyle \frac{\mathrm{P}(\mathrm{O}\mathrm{B}|i)}{1-\mathrm{P}(\mathrm{O}\mathrm{B}|i)}={\displaystyle \frac{\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r},i)\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}|i)}{1-\mathrm{P}(\mathrm{O}\mathrm{B}|\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r},i)\mathrm{P}(\mathrm{O}\mathrm{l}\mathrm{d}\mathrm{e}\mathrm{r}|i)}.}}\end{array}$$

Blanchard (2014, 2020) recommended either using the MROB or the MPOB as predictors in a logistic regression model of participants’ sexual orientation, as opposed to using the raw number of older brothers. Equations (7 and 8) suggest that it would not be difficult to come up with a scenario in which, on average, homosexual participants report just as many (or fewer) older brothers among their older siblings as heterosexual participants do, yet, due to more general group differences in the proportion or number of older siblings, the odds (or the probability) of observing an older brother still are greater for homosexual participants. Consequently, any positive association between the probability of homosexual orientation and the odds (or probability) of observing an older brother among all of the siblings of a participant may be compatible with both models: a specific group difference in the number of older brothers (conditional on the number of older siblings), but also with a more general difference in the number of older brothers and sisters.

Blanchard (2014) defined two complementary indices for older sisters, the modified proportion of older sisters (MPOS) and the modified ratio of older sisters (MROS; see Table 2 for definitions), and interpreted the presence of a statistically significant coefficient for the MROB (MPOB) and the simultaneous lack of a statistically significant coefficient for the MROS (MPOS) as evidence for the hypothesis that the number of older brothers, but not the number of older sisters, are related to the probability of homosexual orientation in men. Declaring a difference between two effects (that of the MROB/MPOB and that of the MROS/MPOS) significant, based on the observation that one effect statistically is significantly different from the null hypothesis, whereas the other is not, is a well-known fallacy in statistical significance testing (Gelman & Stern, 2006). A pattern of statistically significant vs statistically nonsignificant coefficients does not inform about whether the MROB’s (MPOB’s) association with sexual orientation can be taken to be greater than that of the MROS (MPOS).

As discussed next, the rationale for using the MPOB and MROS seems to be based on a common misconception about ratios, namely, that these would adjust for the variable in their denominator (Sollberger & Ehlert, 2016), which was the stated purpose for introducing these ratios into the FBOE literature in the first place (Blanchard, 2014).

Ratios do not adjust for confounding variables

As the sole predictor in a linear model, the MROB (and/or the MROS) is equivalent to including only the interaction term between the number of older brothers (+0.33) and the reciprocal of the number of all other siblings (+1) into the model (Kronmal, 1993; see the corresponding equation in Table 2). However, the constituent variables of the interaction, (#OB + 0.33) and 1/(#OS + #YB + #YS + 1), are omitted from the model, which implies that the regression coefficients for #OB + 0.33 and 1/(#OS + #YB + #YS + 1) are both set to zero. That is, when using the MROB as a predictor, it is not clear whether the statistical effect is driven by the number of older brothers (numerator), or the reciprocal of the number of other siblings (denominator), or some interaction between these two. Similar considerations apply to the use of the MPOB and MPOS.

Ratios (often referred to as “indices”) are ubiquitous in many areas in the social and behavioural sciences, and often their substantive interpretation appears straightforward and meaningful. The statistical analysis of ratios, however, is not at all straightforward (Wiseman, 2009; see also Kronmal, 1993; Sollberger & Ehlert, 2016). Most importantly, it is not the case that ratios adjust for the variable in the denominator (e.g., Kronmal, 1993; Wiseman, 2009). Instead, in order to adjust for the influence of an assumed confounding variable, one could simply add the constituents of a ratio as a predictor variable of its own to the statistical model.

For instance, a researcher might posit that the relationship between homosexual orientation and the number of older brothers should be positive, but that the number of other siblings (#Other) could attenuate the estimate of this relationship (as suggested by Blanchard (2014)). Regressing sexual orientation (via an appropriate link function) on #OB and #Other would adjust for the confounding effect of #Other. This regression model is a simple representation of the assumed theoretical model. In using the MROB (MPOB) for modelling this relationship, the associated regression coefficient does not correspond to the regression coefficient of #OB in the more adequate model, and it is generally not clear how to interpret this regression coefficient. At worst, a spurious relationship between the predictor and outcome variables is introduced (Kronmal, 1993; Wiseman, 2009).

Nevertheless, in trying to interpret the regression coefficient for #OB, it becomes clear that adjusting for the number of other siblings does not rule out the possibility of an older-sibling effect (i.e., an effect of older brothers and sisters). In a logistic regression (as recommended by Blanchard (2014)), a positive regression coefficient for #OB would indicate that–while holding the number of other siblings constant–the logit of the probability of homosexual orientation increases as a function of the number of older brothers. The problem with this model lies in holding #Other constant, as there are numerous combinations of its constituent variables #OS, #YB, and #YS, which all could add up to one and the same value for #Other. That is, if the regression coefficient for #OB were in fact driven by a greater number of older brothers and sisters, participants who have more older brothers would also have more older sisters, whereas fewer younger brothers and younger sisters. In adjusting for #Other, this information is lost. By analogy, the variable #All (the total number of siblings) is equally unfit for ruling out an older sibling effect, since numerous combinations of #OB, #OS, #YB and #YS would sum up to identical values of #All.

Thus, there is no reasonable justification for adjusting for #Other and #All (see also Zietsch, 2018). However, for the goal of quantifying the specific association between the number of older brothers and homosexual orientation, one could adjust for the confounding effect of the number of older siblings, #Older (or, equivalently, the proportion of older siblings; Frisch & Hviid, 2006; Gelman & Stern, 2006; Zietsch, 2018).

To this end, Gelman & Stern (2006, p. 330) suggested the difference between #OB and #OS as one predictor and #Older as a second predictor (in a logistic regression model). If solely the number of older brothers, but not the number of older sisters (or the number of older sisters to a lesser extent than the number of older brothers), were associated with homosexual orientation, then a positive regression coefficient for the difference between #OB and #OS should be observed (the number of older siblings). Alternatively, and to obtain an estimate of the increase in the odds of homosexual orientation, a model with the predictors #OB and #Older could be fitted to the data.

Part I conclusions

Our findings in Part I boil down to two overarching themes. First, assuming the claims made by Blanchard et al. about which relationships between variables should be observable as given, the OBOR, the MROB, and the MPOB are all intended to adjust for the confounding effect of total family (or sibship) size in statistical analysis. Yet, it is the number of older siblings that should be adjusted for instead, as has already been pointed out before (Frisch & Hviid, 2006; Gelman & Stern, 2006; Zietsch, 2018). Second, ratios do not adjust for the variable in the denominator, which is a common misconception surrounding the use of ratio variables (e.g., Wiseman, 2009). Using basic probability calculus, the statistical clarification provided here has thus shown that ratios better should not be used in the way they are used in extant research on the FBOE, and that analyses need to adjust for the number of older siblings, but not for total family (or sibship) size.

Methods

A detailed description of the simulation study is provided in Supplement S5. We investigated scenarios of (1) a 33% increase in the odds of homosexual orientation per older brother (θ = 0.33), as is reported throughout the literature (Blanchard, 2001; Blanchard et al., 1996); (2) no distinct association between the number of older brothers and sexual orientation (θ = 0); and (3) a 33% decrease in the odds of homosexual orientation per older brother (θ = −0.33). Each older brother increased these odds by a factor of (1 + θ). For the homosexual sample, the study employed three different values for the proportion of older siblings in a given sample, π, namely (1) 0.5, (2) 0.6, and (3) 0.7, while for the heterosexual sample π was fixed at 0.5. That is, we simulated differences in the number (or proportion) of older siblings, assuming that the proportion of older siblings is equal to the proportion of younger siblings among homosexual individuals, but increasing in the homosexual group.

The median of the mean numbers of all siblings for homosexual participants, μ, across the 45 samples in Blanchard (2018a, 2018b), and Blanchard et al. (2021) was 2.45. In the equal condition (1) of the simulation study, this value served as the mean number of all siblings for both the homosexual and heterosexual group. In the unequal condition (2), the mean number of siblings for each group were taken from the Mismatch 2 sample in Blanchard (2014), where the mean number of siblings in the homosexual group was 2.19, and the mean number of siblings in the heterosexual group was 3.31. Blanchard used the Mismatch 2 sample to demonstrate the inability of tests for mean differences and logistic regression to detect an older brother effect and to promote the use of the MROB and MPOB.

Thus, there were 3 × 3 × 2 = 18 possible combinations of conditions in this simulation study. For each combination, we fitted 10 different models (see Table 3).

Models and evaluation

Table 3 lists the equations of the models we fitted to the simulated data. For each model and each combination of μ, π, and θ, 1,000 replications with a sample size of 700 participants per group each were carried out (i.e., 18,000 replications in total).

The performance of the models was primarily assessed with respect to inferring the state of θ with null-hypothesis significance testing (NHST), as it is the case throughout the FBOE literature. Only Models 1, 6, 7, and 9 in Table 3 could be interpreted as providing an estimate of θ (by exponentiating the estimates and subtracting (1)). However, in the FBOE literature, studies rarely report estimates of $\mathrm{\theta }$. Instead, all of the currently advocated practices involving the OBOR, the MROB, and MPOB infer about θ ≠ 0 given a test of a parameter β, which does not correspond to θ. β represents any of the regression coefficients of interest in Table 3. It is evident that the hypothesis of β = 0 vs. β ≠ 0 need not correspond to a hypothesis test about θ = 0 vs. θ ≠ 0. Inferences about θ are nevertheless made in the FBOE literature based on inferences about β, using $\mathrm{\alpha }$ = 0.05 (two-sided, with only a few exceptions). The simulation aims to show that relying on NHST of the wrong model inevitably leads to biased conclusions about the true state of θ.

Summary of Parts I and II

The results of Parts I and II demonstrate that currently recommended effect sizes, variable transformations, and their implementation in statistical models in the FBOE literature neither are needed nor adequate for quantifying and distinguishing the specific association between the number of older brothers and homosexual orientation (i.e., the theoretical estimand) from a more general association of older siblings and sexual orientation. Furthermore, models adjusting for #Other or #All (Models 6 and 7) can be regarded as misspecifications and thus are prone to lead to false or biased conclusions.

We also demonstrated the necessity of adjusting for the number of older siblings, as group differences in the proportion of older siblings are bound to lead to the false decision that solely the number of older brothers is associated with sexual orientation when using the OBOR, MROB, or MPOB. It is quite plausible to suppose such group differences, since the difference in the median proportions of older siblings across the 45 homosexual and heterosexual samples in Blanchard (2018a, 2018b), and Blanchard et al. (2021) amounted to a difference of 10 percentage points. We note that such differences in the number of older siblings are not independent of differences in the number of older brothers between groups, and it would therefore be too simplistic to conclude that such differences solely are due to the associations between the number of older sisters and homosexual orientation and the number of older brothers and homosexual orientation being of equal magnitude.

The OBOR, MROB, and MPOB fail to quantify the very thing, which they are advertised as quantifying. These indices and their recommended usage are not rigorous enough in exposing the proclaimed specific association between the number of older brothers and sexual orientation to conditions under which observations could indicate the absence of the specificity of such an association.

Given that the most extensive meta-analyses to date (totaling 45 samples of homo- and heterosexual males; Blanchard, 2018a, 2018b; Blanchard et al., 2021) all employed the OBOR, it follows that these meta-analyses actually provided no conclusive evidence for the key observation of homosexual men having more older brothers relative to older sisters than heterosexual men. Our findings, of course, do not preclude the possibility that there may nevertheless be a meaningful stronger association between the number of older brothers and sexual orientation than between the number of older sisters and sexual orientation in men as required by the FBOE. This is one of the reasons for why in Part III we provide a set of new meta-analyses, using more adequate effect-size metrics, along with an advanced meta-analytic framework.

Methods

Data analysis

In choosing an effect size for meta-analysis, we followed the approach of Blanchard et al. (2020) in modelling reported siblings as independent observations (Blanchard, 2018a, 2018b). Note that no further information on the reported siblings was available than the number of times each sibling type occurred in a sample.

We quantified the difference in the number of older brothers between the homosexual and heterosexual group through the natural logarithm of the odds ratio (lnOR) of older brothers vs. older sisters, as described in Eq. (6) above and Model 10 in Table 3. As shown in Parts I and II, this lnOR adjusts for the number (or the proportion) of older siblings and remains unaffected by group differences in the total number of siblings. In order to compute this lnOR, the total numbers of older brothers and sisters of the homosexual and heterosexual participants are needed, instead of the mean numbers listed in Tables 5 and 6. These totals can be obtained by multiplying the Tables 5 and 6 mean values with the corresponding sample size (column N) and rounding to the nearest integer.¹⁰ Alternatively, one could use the difference in the proportion of older brothers among older siblings between homosexual and heterosexual men and women as an effect size. Corresponding meta-analyses using the difference in proportions as an effect size can be found in Supplement S6 (Subsection 1).

We fitted random-effects models to five specific study sets comprising the 31 male samples in Blanchard (2018a), the six male probability samples in Blanchard (2018b), the 24 male samples in Blanchard et al. (2021), the entire set of 64 male samples, and the 17 female samples (results of fixed-effect analyses are provided in Supplement 6, Subsection 2). For all random-effects models, the between-sample variance ${\tau }^2$ was estimated using the restricted maximum-likelihood (REML) estimator as implemented in the R package metafor (Viechtbauer, 2010).

To further assess the extent to which the effect might differ between men and women, we also computed a combined meta-analysis over the entire set of 81 (male and female) samples, with participant sex (male vs. female) serving as a predictor variable. Because each female sample was embedded in a study alongside a male sample (to our knowledge, there are no all-female studies in this literature), there may be some degree of dependency between effect-size estimates from the same study. This study-level dependency can be modelled using a three-level-meta-analysis (e.g., Van den Noortgate et al., 2013), which requires the estimation of a second, study-level, variance component ( ${\tau }_2^2$) in addition to the first, sample-level, variance component ( ${\tau }_1^2$).

Specification-curve and multiverse meta-analysis

Invariably, the results of any meta-analysis are subject to a sequence of researcher-dependent decisions (or degrees of freedom) pertaining to the inclusion of studies, the assumed meta-analytic model and the estimation of its parameters (Voracek, Kossmeier & Tran, 2019). For instance, Blanchard (2018a) defined the proximal assessment of sexual orientation via same-sex marriage as sufficiently invalid to justify the exclusion of one of the largest samples (Frisch & Hviid, 2006) available. In contrast, samples comprised of individuals whose sexual orientation was classified using phallometric testing (Blanchard et al., 2006, 2000) were considered worthy of inclusion; as were studies which categorized children as young as 2 years of age as “prehomosexual” (Blanchard et al., 1995; VanderLaan et al., 2014), based on the “femininity” rating of their behaviours; or studies which drew on expert ratings of criminal or clinical records (Blanchard & Bogaert, 1998; Blanchard & Sheridan, 1992; Bogaert et al., 1997) to categorize individuals. While the inclusion criteria of Blanchard (2018a) are based on subject-matter considerations, other researchers could disagree on substantive grounds and argue that same-sex marriage constitutes a less error-prone indicator of homosexual orientation than changes in penile blood volume upon exposure to certain stimuli (as is the case in phallometric testing) and that the latter, not the former study, should be noneligible for meta-analysis. Consequently, other researchers could opt to analyse a different set of samples than Blanchard (2018a) and therefore obtain different results. The perceived validity of differing approaches for classifying sexual orientation is just one out of various crucial steps in the sequence of researchers’ decisions, which lead to a particular meta-analysis.

With regards to the decision of which meta-analytic model should be fitted to the data, Blanchard (2018a, 2018b) and Blanchard et al. (2021) solely reported the results of random-effects meta-analyses. Even though some researchers may argue that in the presence of effect-size heterogeneity random-effects models need to be fitted to the data, other researchers may disagree and choose a fixed-effect model instead (see Hedges & Vevea, 1998; Rice, Higgins & Lumley, 2018).

There are many more such decision steps or researcher degrees of freedom in specifying a meta-analysis, thus demanding researchers to choose one out of several possible alternatives and, for each combination of such alternatives, specify a different meta-analytic model. Thus, for each reported meta-analysis, there is a sizeable collection of unreported, alternatively specified, meta-analyses, which may lead to considerably different conclusions.

Instances of disagreement among researchers about how a given set of data should be analyzed are not confined to the realm of meta-analysis. The numerous degrees of freedom researchers have at their hands in primary data analysis are well documented (Silberzahn et al., 2018; Simmons, Nelson & Simonsohn, 2011; Simonsohn, Simmons & Nelson, 2020; Steegen et al., 2016). Pointing out the issue of researcher degrees of freedom in a given analysis may occasionally be taken up as an accusation of questionable research practices (e.g., “p-hacking,” “fishing for statistical significance”; Gelman & Loken, 2013). Yet, even in the absence of questionable research practices, researchers may still disagree considerably with respect to the goals of a given study and how the data should be analysed, and thus how the research question should be optimally answered (e.g., Auspurg & Brüderl, 2021; Silberzahn et al., 2018).

Given the obvious structural similarities between primary data analysis and meta-analysis (effect sizes correspond to individual participant data), specification-curve and multiverse meta-analysis (Voracek, Kossmeier & Tran, 2019) builds on two almost identical approaches for tackling the issue of researcher degrees of freedom, taken from the realm of primary data analysis, namely specification-curve analysis (Simonsohn, Simmons & Nelson, 2020) and multiverse analysis (Steegen et al., 2016). These approaches identify all reasonably specified analysis plans and report all of the corresponding results. The set of results obtained then serves as the foundation for an extra data-analytic step of inferring the state of the substantive research question. For the identification of reasonable analysis plans, Simonsohn, Simmons & Nelson (2020) and Voracek, Kossmeier & Tran (2019) suggested to survey potential decision steps in the process of devising an analysis and to determine all possible combinations of the alternatives offered at each decision step. All nonredundant data-analytic specifications retrieved in this manner are then considered as reasonable specifications. For an application of specification-curve analysis in the closely related field of birth-order effects on personality traits, see Rohrer, Egloff & Schmukle (2017).

We follow the terminology introduced by Voracek, Kossmeier & Tran (2019) and refer to the decision steps in specifying a meta-analysis as Which and How factors. The Which factors (or, the data multiverse) refer to decisions leading up to the set of samples (or studies) included in a meta-analysis (i.e., which primary studies are included in a meta-analysis) and are mostly determined by features of the samples themselves. The How factors (or, the model multiverse) subsume the decisions leading up to the technical implementation of a meta-analysis over a given set of samples, such as the selected meta-analytic model and the effect-size metric (i.e., how the primary studies are meta-analysed). Table 6 lists the 10 Which factors we derived by scrutinizing the inclusion criteria of prior related meta-analyses (Blanchard, 2018a, 2018b; Blanchard & VanderLaan, 2015; Blanchard et al., 2021), from one critical commentary by Zietsch (2018), and from our own concerns.

1. The factor In Previous Meta-Analysis indicated whether a sample had been included in any of the seven previous meta-analyses (Blanchard, 2004, 2018a, 2018b, 2020; Blanchard & VanderLaan, 2015; Jones & Blanchard, 1998). Given that all of these were conducted by the same group of researchers (spearheaded by Blanchard), some samples may well have been judged more favourable for meta-analytic inclusion (possibly based on prior convictions that the FBOE is real).

2. The factor Lab indicated whether Blanchard, Bogaert, Zucker, or VanderLaan were listed as co-authors of the publications the samples appeared in. These researchers frequently collaborated on FBOE publications.

3. The factor Archives indicated whether a study was published in the journal Archives of Sexual Behavior. Thirty-two out of the 64 male samples appeared in articles published in this journal. As of mid-2022, Zucker served as Editor and Blanchard, Bogaert and VanderLaan as Editorial Board members of this journal (https://www.springer.com/journal/10508/editors).

4. The factor Feminine indicated whether a given sample had been declared a “feminine” sample by Blanchard (2018a), who theorized that the association between older brothers and homosexual orientation is stronger in “feminine” samples, as opposed to “non-feminine/cisgender” samples.

5. The factor GDY indicated whether samples comprised individuals diagnosed with GDY. With exception of Vasey & VanderLaan (2007) and the replication sample in VanderLaan & Vasey (2011), all “feminine” samples comprised individuals diagnosed with GDY. It might be argued that GDY samples should be excluded from meta-analysis.

6. The factor Clinical indicated whether participants were recruited in a clinical setting. This issue was pointed out by Zietsch (2018), who questioned the representativeness of such samples.

7. The factor Classification distinguished between four different methods for classifying sexual orientation, as used in the FBOE literature. We assigned each sample’s classification scheme to one of the three broad categories (i.e., factor levels), which reflect different levels of reliability of classification. The factor level Single-Item referred to studies which determined sexual orientation via a single piece of information (such as asking participants which category they would ascribe themselves to, or whether participants were in a same-sex relationship; Frisch & Hviid, 2006), as opposed to a summary score of multiple items or expert ratings. The factor level Various subsumed samples for which sexual orientation had been determined by computing summary scores over multiple items (such as the Kinsey scale), or which combined participants whose sexual orientation had been assessed via a mixture of methods, such as the “Blanchard” subsample in Blanchard et al. (2006). The factor level Indirect was assigned to samples wherein participants were categorized by either employing behavioural measures (e.g., phallometric testing, feminine behaviours in children as predictors of homosexual orientation), criminal offence records, or expert ratings to determine sexual orientation. Also, two studies (Blanchard & Bogaert, 1996a, Study 2; Schagen et al., 2012) did not classify sexual orientation per se (but used GDY diagnosis as a proxy). These studies are indicated by the factor level none.

8. The factor Probability indicated whether a sample was declared a “probability sample”, i.e., samples comprised of participants who were randomly selected with respect to sexual orientation Zietsch (2018), or samples that had been obtained through some kind of random sampling procedure (Skorska & Bogaert, 2020; Xu, Norton & Rahman, 2019. All other samples were regarded as non-probabilistic convenience samples.

9. The factor Stopping Rule indicated whether the authors of a primary study or Blanchard (2018a) assumed the presence of a stopping rule, such as one-child policies (e.g., Xu & Zheng, 2014), or parents ceasing to reproduce after having at least one male and one female offspring (Blanchard & Lippa, 2007). Blanchard (2018a) noted that the proposed excess of older brothers may go undetected in the presence of a stopping rule.

10. Finally, the factor Pedo-/Hebephiles indicated whether a sample comprised individuals who were classified as pedophile or hebephile by Blanchard et al. (2021), who investigated this factor as a moderator.

We considered only a single factor in terms of how to analyse these subsets of samples from the corpus of primary studies (i.e., the How factor), namely, fitting a random-effects model, as in Blanchard (2018a, 2018b) and Blanchard et al. (2021), or a fixed-effect model.

Counting all the factor-level combinations (considering also their union, where applicable) of the 10 Which factors and the single How factor resulted in a total of 3 (In Previous Meta-Analysis) × 3 (Lab) × 3 (Archives) × 3 (Feminine) × 3 (GDY) × 3 (Clinical) × 5 (Classification) × 3 (Probability) × 3 (Stopping Rule) × 3 (Pedo-/Hebephiles) × 2 (fixed effect vs. random effects) = 196,830 possible meta-analytic specifications for the 64 male samples. This set of specifications comprised 1,628 unique, non-redundant combinations of at least two from the 64 male samples in Table 4 (i.e., the set of reasonable specifications; see analysis code provided online: https://osf.io/3wnhu/). Owing to the lack of previous meta-analyses of an older brother effect in female samples, the factor In Previous Meta-Analysis was not applicable when determining the subset of possible meta-analyses for female samples, as was the case for the factors Feminine and Pedo-/Hebephile. For the remaining applicable factors, there were 3 (Lab) × 3 (Archives) × 3 (GDY) × 3 (Clinical) × 5 (Classification) × 3 (Probability) × 3 (Stopping Rule) × 2 (fixed effect vs. random effects) = 7,290 possible meta-analytic specifications. Of these, 212 turned out to be unique and non-redundant specifications comprising at least two from the 17 available female samples.

All analyses were carried out in R (version 4.0.3; R Core Team, 2021). Meta-analyses were computed using the metafor package (version 3.0.2; Viechtbauer, 2010), graphical displays of results were generated using metaviz (Version 0.3.1; Kossmeier, Tran & Voracek, 2021), and ggplot2 (Wickham et al., 2019; see online material for code; https://osf.io/3wnhu/).

We follow recommendations (Morey et al., 2016) of stating how the reported interval estimates were obtained. Confidence intervals for summary lnORs were based on the normal approximation to the sampling distribution of the point estimator of the summary effect (i.e., the lnORs; these intervals being the default option in metafor; Viechtbauer, 2010). Confidence intervals for the between-sample standard deviation, $\tau$, were based on the profile likelihood of ${\tau }^2$ (Hardy & Thompson, 1996), which is the default when fitting meta-analyses using the rma.mv() function in the metafor package (Viechtbauer, 2010). Prediction intervals were obtained using the quantiles of the standard normal distribution. This is also the default option in metafor.

Further, in reporting the results of the meta-analytic specification-curves for male and female samples, we also made use of the S value (Rafi & Greenland, 2020), a mapping of the p value onto an additive scale without an upper bound, given by S = log₂(1/p), where log₂ denotes the base two logarithm. Like the p value, the S value indicates the degree of compatibility between the observed data and the statistical null model, which not only comprises the null hypothesis, but also a host of statistical background assumptions, pertaining to the assumed data generating process, the absence of selection bias, and the sufficiently correct specification of the model, among other possible assumptions (Greenland et al., 2022; Rafi & Greenland, 2020). The S value is interpretable as bits of information against the null model. For instance, a (just significant) p value of 0.05 corresponds to an S value of 4.32, i.e., conveys 4.32 bits of information against the null model. The unit of bits can further be interpreted as a surprisal value for observing a given S value. An S value of 4.32 would be as surprising as observing four heads out of four tosses of a fair coin–which is not that much surprising and certainly would not already justify the conclusion that the coin is unfair.

Meta-analyses of the five specific study sets

Table 7 provides a summary of the meta-analyses over the five specific subsets (Men included in Blanchard, 2018a; Men included in Blanchard, 2018b; Men included in Blanchard et al., 2021; Men full set; Women full set) of male and female samples using lnOR as an effect size (the fixed-effect analysis results are provided in Supplement 6, Subsection 2). There was little disagreement between the estimated summary lnORs across all four specific male study sets. Depending on which specific study set is interpreted, the odds for observing an older brother among the set of all older siblings reported by homosexual participants (male or female) were between 7% (for the Women full set) and 17% (for the 31 samples included in Blanchard, 2018a) greater than those same odds for the heterosexual participants. However, the 95% CIs suggest that these estimates were compatible with a 6% decrease as well as with a 35% increase (i.e., the respective lower and upper bounds of the 95% CI of the summary estimate for the six probability samples included in Blanchard, 2018a) for these odds. While the precision of these summary lnORs varied considerably across the specific study sets (the shortest 95% CI spanned 0.11 lnORs, the widest 95% CI spanned 0.36 lnORs), all of the point estimates were directionally consistent with a greater number of older brothers among the older siblings of homosexual men as predicted by the FBOE. This greater number of older brothers was not specific to the male samples, however, as the summary lnORs for the female samples were also consistent with a greater number of older brothers in homosexual women, with an estimated 7% greater odds of observing an older brother among the older siblings of homosexual vs heterosexual women.

As stated above, we also combined the male and female samples into a single data set and fitted random-effects meta-anlyses using sex (male vs female) as a predictor variable to compare lnORs between these two groups of samples. The regression coefficient for sex can be interpreted as an estimate of the difference in lnOR between male and female samples.

This difference amounted to 0.07, 95% CI [0.08–0.22], in the random-effects model. The three-level meta-analysis returned results almost identical to this model, with a between-group difference of 0.06, 95% CI [−0.08 to 0.20], ${\hat{\tau }}_1$ = 0.12, 95% CI [0.00–0.22], ${\hat{\tau }}_2$ = 0.11, 95% CI [0.00–0.22]. These results suggest that conditional on the assumed statistical model(s), the available data do not convey enough information to confidently declare that the male and female samples differ with respect to the greater number of older brothers among the older siblings of homosexual vs. heterosexual individuals. To put things in perspective: this inconclusive finding is based on the data of 30,000 homosexual individuals (17,134 older brothers and 15,286 older sisters reported). Moreover, it calls into question previous confident claims about the excess of older brothers in homosexual individuals being specific to men (e.g., Bogaert et al., 1997; Bogaert & Skorska, 2011).

With regards to effect-size heterogeneity, the $\hat{\tau }$s for the random-effects models were all of the same magnitude as the corresponding summary lnORs themselves. Population effects thus appeared highly variable across samples. We computed 95% prediction intervals (PIs) to illustrate effect-size heterogeneity. A 95% PI can be interpreted as providing an estimate of where 95% of the population effects of future studies might fall (IntHout et al., 2016). The 95% PI suggest a wide range of effects, which are consistent with both a greater as well as a lower number of older brothers relative to older sisters among homosexual vs. heterosexual men, depending on the population.

To facilitate the comparison of the magnitude of the OR for older brothers vs older sisters to that of the OBOR, Table 8 lists the results of the meta-analyses using the lnOBOR as an effect size. It is evident that if the OBOR is falsely interpreted as quantifying the effect of interest, one would exaggerate the magnitude of the effect in male samples considerably.

Blanchard (2018a; fourth meta-analysis) found that some of the heterogeneity of the OBOR could be explained by two distinct subgroups of samples, namely “feminine” samples and “non-feminine/cisgender” samples. The re-analysis of Blanchard (2018a) using the lnOR returned random-effects summary estimates of lnOR = 0.09, 95% CI [0.03–0.17], and lnOR = 0.29, 95% CI [0.17–0.41]. Across all 47 samples which were not denoted “feminine” by Blanchard (2018a), we obtained random-effects estimates of lnOR = 0.08, 95% CI [0.02–0.13]. This subgroup analysis suggests a sizeable difference in the effect between “feminine” and “non-feminine/cisgender” individuals.

Figure 3 displays a raindrop forest plot (Kossmeier, Tran & Voracek, 2021; Schild & Voracek, 2015; see Kossmeier, Tran & Voracek, 2020, for the advantages of this graphical device) accompanying the random-effects meta-analysis of all male-sample lnOR effect sizes (the respective plot of the fixed-effect meta-analysis is provided in Supplement S6, Subsection 2). Therein, the observed effect-size estimates are arranged (from top to bottom) in ascending order of their meta-analytic weight. Two patterns are visible in Fig. 3. First, only seven out of 64 (11%) confidence intervals’ lower bounds were greater than zero. That is, only seven samples returned a statistically significant effect in the right direction when using the (adequate) lnOR effect size metric. This contrasts the “consistent” and “reliable” (e.g., Blanchard, 2004; Blanchard & VanderLaan, 2015) greater number of older brothers relative to older sisters among homosexual vs. heterosexual men previously reported in all these primary studies.

Second, with increasing study weight, the effect estimates seemed to approach zero. This is indicative of small-study effects (i.e., smaller, less precise, studies report larger effects than larger, more precise, studies). The funnel plot in Fig. 4 further illustrates this observation. Assuming the absence of publication bias, one would expect the observed effect sizes to be distributed symmetrically around the summary effect estimate (Vevea, Coburn & Sutton, 2019). This was visually not the case in Fig. 4 and the observed distribution of effect sizes also appeared to be incompatible with a model assuming symmetry of effect sizes, as was indicated by Egger et al. (1997)’s regression: Regressing the observed effect sizes on their standard errors, the regression coefficient for the standard error was 0.80, 95% CI [0.28–1.31], in the random-effects model. This indicates that (when rescaling these coefficients by dividing them by 10), on average, two effect sizes differing by 0.1 standard errors were expected to differ by 0.08 lnORs. This difference was sizeable in comparison to all summary estimates listed in Table 7.

Specification-curve and multiverse meta-analyses

To demonstrate which other sample-specific factors could potentially influence the estimates of the meta-analytic summary estimates and thus the qualitative conclusion about the presence of a greater number of older brothers relative to older sisters among homosexual vs heterosexual men, we now consider the results of the specification curve meta-analyses described above.

The results of the 1,628 and 212 non-redundant meta-analytic specifications for the set of male and female samples are summarized in Figs. 5 and 6. From top to bottom, the panels encode each specification’s estimated summary effect and 95% confidence interval, number of samples included (density plot), and combination of factor levels (tile plot). The point estimates are arranged in ascending order and are connected by a black solid line, representing the specification curve. In addition, each specification is color-coded, using hues of six distinct colors (red, orange, yellow, green, blue, violet) on the spectrum from red to violet. A specification’s color and hue are indicative of the number of samples included in the meta-analysis, with the number of samples increasing in the following order: Red, orange, yellow, green, blue, violet.

In terms of magnitude, the estimated mean lnORs of the various subsets of male samples ranged from −0.29 (OR = 0.74) to 0.93 (OR = 2.53), with an interquartile range of 0.05 to 0.14 (OR = 1.05 to 1.15). In total, 96% of the estimated means were greater than 0, and 47.4% of these had 95% CIs that did not include 0 (i.e., estimated summary effects greater than 0, with p < 0.05, two-tailed).

The overall pattern of the specification curve seems to indicate that the specifications with many samples and/or narrow confidence intervals were closer to zero, as opposed to specifications with only a few samples and/or wider confidence intervals.

The Which factors GDY and Feminine systematically influenced the magnitude of the summary effect. Specifications that included samples containing gender-dysphoric individuals, or “feminine” (Blanchard, 2018b) samples only, produced some of the largest effects. When including GDY as a predictor, the random-effects summary estimates for the 52 non-GDY and the 12 GDY samples were lnOR = 0.09, 95% CI [0.03–0.15], and lnOR = 0.31, 95% CI [0.16–0.46].

Figure 5 also indicates that as the meta-analytic summary estimates increased, so did the frequency of specifications encompassing only those samples, which had already been included in at least one of the previous six meta-analyses (Blanchard, 2004, 2018a, 2018b; Blanchard et al., 2020, 2021; Blanchard & VanderLaan, 2015; Jones & Blanchard, 1998). Including the Which factor In Previous MA as a moderator in the meta-analysis consisting of all 64 samples, the random-effects summary estimates were lnOR = 0.00, 95% CI [−0.13 to 0.13], and lnOR = 0.15, 95% CI [0.09–0.21] (a difference of 0.15, 95% CI [−0.01 to 0.29], in favour of studies included in previous meta-analyses).

The meta-analytic specifications across the female samples returned similar results. The summary effects ranged from −0.22 to 0.54, with an interquartile range of 0.05 to 0.16. A total of 84.9% of estimated means were greater than 0 of which 32.2% did not include zero. Again, tighter confidence intervals were observed for summary effects close to zero. No systematic pattern of Which/How factor level combinations was apparent.

Voracek, Kossmeier & Tran (2019) recommended assessing the overall evidence from a specification-curve and multiverse meta-analysis against the respective null model with parametric bootstrapping (inferential specification-curve plots) and an additional histogram of p values. Here, we extend the idea of the histogram of p values to the S value. A derivation of the distribution of S values under the null model (exponential with rate parameter λ = ln(2)), with which the observed distribution can be compared, is provided in Supplement S7. The inferential specification-curve plots are presented in the Supplement S8.

Figure 7 (left) displays histograms and kernel density estimates of the 1,628 observed S values, as well as the probability density function of the exponential distribution with rate parameter λ = ln(2). Clearly, the distribution of the observed S values deviates from both the expected and simulated distribution. Thus, conditional on all statistical background assumptions being met¹¹ (Greenland et al., 2022), a true mean lnOR of 0 (i.e., the null hypothesis) appears incompatible with the results of the specification curve meta-analysis. We can get an idea of the degree of incompatibility (or distance) between the data and test model by thinking of the distribution of these S value in terms of surprise. The interquartile range of observed S values was (2.18 to 6.22), the median S value was 4.02. That is, under the null model, the central 50% of results would be approximately as surprising as observing all heads in 2, 3, 4, 5 or 6 tosses of a fair coin.

The interquartile range of the distribution of S values for the female samples (Fig. 7, right) was 0.73 to 4.47, with a median of 1.67.

A specification curve analysis for the difference between the summary estimates in men and women can be found in Supplement 9.