Shiv has come out of retirement with an excellent post on Meghan Murphy. Go read it. I’m a fan, but I do have a critique: I don’t think it goes far enough.

Customer service: is run by John and Stacy

Shiv has come out of retirement with an excellent post on Meghan Murphy. Go read it. I’m a fan, but I do have a critique: I don’t think it goes far enough.

I’ve seen some confusion about “fourth wave feminism” and what it means. As it stands, the term has two separate meanings. [Read more…]

I thought I knew how this post would play out. EssenceOfThought has gotten some flack for declaring Stephen Woodford to be a “violent transphobe,” which I didn’t think they deserved. They gave a good defense in one of their videos, starting off with a definition of violence.

You see, violence is defined as the following by the World Health Organization. Quote; “the intentional use of physical force

or power, threatened or actual, against oneself, another person, or against a group or community, that either results in, or has a high likelihood of resulting in injury, death, psychological harm, maldevelopment or deprivation.”

EoT points out that controlling someone’s behaviour or social networks by using their finances as leverage can be considered economic violence. They also point out that using legislation to control access to abortion can be considered legislative violence, as it deprives a person of their right to bodily autonomy. And thus, as EoT explains,

When you exclude trans women from women’s sports you’re not simply violating numerous human rights. You’re designating them as not real women, as an invasive force coming to take what doesn’t belong to them. You are cultivating future transphobic violence.

Note the air gap: “cultivating violence” and “violence” are not the same thing, and the definition EoT quoted above places intent front-and-centre. EoT bridges the gap by pointing out they gave Rationality Rules several months to demonstrate he promoted violent policies out of ignorance, rather than with intent. When “he [doubled] down on his violent transphobia,” EoT had sufficient evidence of intent to justify calling him a “violent transphobe.”

At this point I’d shore up their one citation with a few more. This decoupling of physical force and violence is not a new argument in the philosophy and social sciences literature.

Violence often involves physical force, and the association of force with violence is very close: in many contexts the words become synonyms. An obvious instance is the reference to a violent storm, a storm of great force. But in human affairs violence and force, cannot be equated. Force without violence is often used on a person’s body. If a person is in the throes of drowning, the standard Red Cross life-saving techniques specify force which is certainly not violence. To equate an act of rescue with an act of violence would be to lose sight entirely of the significance of the concept. Similarly, surgeons and dentists use force without doing violence.

Violence in human affairs is much more closely connected with the idea of violation than with the idea of force. What is fundamental about violence is that a person is violated. And if one immediately senses the truth of that statement, it must be because a person has certain rights which are undeniably, indissolubly, connected with being a person. One of these is a right to one’s body, to determine what one’s body does and what is done to one’s body — inalienable because without one’s body one would cease to be a person. Apart from a body, what is essential to one’s being a person is dignity. The real dignity of a person does not consist in remaining “dignified”, but rather in the ability to make decisions.

Garver, Newton. “What violence is.” The Nation 209.24 (1968): 819-822.As a point of departure, let us say that

violence is present when human beings are being influenced so that their actual somatic and mental realizations are below their potential realizations. […]The first distinction to be made is between physical and psychological violence. The distinction is trite but important mainly because the narrow concept of violence mentioned above concentrates on physical violence only. […] It is useful to distinguish further between ’biological violence’, […] and ’physical violence as such’, which increases the constraint on human movements – as when a person is imprisoned or put in chains, but also when access to transportation is very unevenly distributed, keeping large segments of a population at the same place with mobility a monopoly of the selected few. But that distinction is less important than the basic distinction between violence that works on the body, and violence that works on the soul; where the latter would include lies, brainwashing, indoctrination of various kinds, threats, etc. that serve to decrease mental potentialities. […]

We shall refer to the type of violence where there is an actor that commits the violence as

personalordirect, and to violence where there is no such actor asstructuralorindirect. In both cases individuals maybe killed or mutilated, hit or hurt in both senses of these words, and manipulated by means of stick or carrot strategies. But whereas in the first case these consequences can be traced back to concrete persons as actors, in the second case this is no longer meaningful. There may not be any person who directly harms another person in the structure. The violence is built into the structure and shows up as unequal power and consequently as unequal life chances.

Galtung, Johan. “Violence, peace, and peace research.” Journal of peace research 6.3 (1969): 167-191.

This expansive definition of “violence” has been influential, Galtung’s fifty-year-old paper from above has been cited from over 6,000 times according to Google Scholar. “Influential” is not a synonym for “consensus,” however.

Nearly all inquiries concerning the phenomenon of violence demonstrate that violence not only takes on many forms and possesses very different characteristics, but also that the current range of definitions is considerable and creates ample controversies concerning the question what violence is and how it ought to be defined (…). Since there are so many different kinds of violence (…) and since violence is studied from different actor perspectives (i.e. perpetrator, victim, third party, neutral observer), existing literature displays a wide variety of definitions based on different theoretical and, sometimes even incommensurable domain assumptions (e.g. about human nature, social order and history). In short, the concept of ‘violence’ is notoriously difficult to define because as a phenomenon it is multifaceted, socially constructed and highly ambivalent. […]

Violence is socially constructed because who and what is considered as violent varies according to specific socio-cultural and historical conditions. While legal scholars may require narrow definitions for punishable acts, the phenomenon of violence is invariably more complex in social reality. Not only do views about violence differ, but feelings regarding physical violence also change under the influence of social and cultural developments. The meanings that participants in a violent episode give to their own and other’s actions and experiences vary and can be crucial for deciding what is and what is not considered as violence since there is no simple relationship between the apparent severity of an attack and the impact that it has upon the victim. For example, in some cases, verbal aggression may prove to be more debilitating than physical attack.

De Haan, Willem. “Violence as an essentially contested concept.” Violence in Europe. Springer, New York, NY, 2008. 27-40.

A major objection to this inclusive definition of violence is that it makes everything violence, creating confusion instead of clarity. One example:

If violence is violating a person or a person’s rights, then every social wrong is a violent one, every crime against another a violent crime, every sin against one’s neighbor an act of violence. If violence is whatever violates a person and his rights of body, dignity, or autonomy, then lying to or about another, embezzling, locking one out of his house, insulting, and gossiping are all violent acts.

Betz, Joseph. “Violence: Garver’s definition and a Deweyan correction.” Ethics 87.4 (1977): 339-351.

The problem with this objection is that it assumes violence is binary: things are either violent, or they are not. Almost nothing in life falls in a binary, sex included, so a much more plausible model for violence is a continuum. I’m convinced that even the people who buy into a violence binary also accept that violence falls on a continuum, as I have yet to hear anyone argue that murder and wet willies are equally bad. Thus eliminating the binary and declaring all violence to fall on a continuum is a simpler theory, and by Occam’s razor should be favoured until contrary evidence comes along.

The other major objection is that while not every human society agrees on what constitutes violence, all of them agree that physical violence is violence. Sometimes this objection can be quite subtle:

Albeit rare, there are cases of violence occurring without rights being violated. This point has been made by Audi (1971, p. 59): ‘[while] in the most usual cases violence involves the violation of some moral right …there are also cases, like wrestling and boxing, in which even paradigmatic violence can occur without the violation of any moral right’.

*Bufacchi, Vittorio. “Two concepts of violence.” Political Studies Review 3.2 (2005): 193-204.*

That quote only works if you think wrestling is paradigmatic, something everyone agrees counts as violence. Wrestling fans would disagree, and either point to the hardcore training and co-operation involved or the efforts made to prevent injury, depending on which fandom you were querying. Societies definitely disagree on what physical acts count as violence, and even within a single country physical acts that are considered horrifically immoral to many today were perfectly acceptable to many a century ago. This pragmatic argument can also be turned on its head, by pointing out that if violence is binary then we wouldn’t expect a correlation between (for example) hostile views of women and violence towards women. If a violence continuum exists, however, such a correlation must exist.

Studies using Glick and Fiske’s (1996) Ambivalent Sexism Inventory, which contains different subscales for benevolent and hostile sexism, support this idea. Studies have found that greater endorsement of hostile sexism predicted more positive attitudes toward violence against a female partner (Forbes, Jobe, White, Bloesch, & Adams-Curtis, 2005; Sakalli, 2001). Other studies of IPV among college samples have found that men with more hostile sexist attitudes were more likely to have committed verbal aggression (Forbes et. al., 2004) and sexual coercion (Forbes & Adams-Curtis, 2001; Forbes et al., 2004).

Allen, Christopher T., Suzanne C. Swan, and Chitra Raghavan. “Gender symmetry, sexism, and intimate partner violence.” Journal of interpersonal violence 24.11 (2009): 1816-1834.

At this point in the post, though, I was supposed to pump the breaks a little. People have certain ideas in mind when you say “violence,” I’d say, and would likely equivocate between physical and non-physical violence. This would poison the well. Of course you can’t change language or create awareness by sitting on your hands, so EssenceOfThought were 100% in the right in arguing Rationality Rules was a violent transphobe, but at the same time I wasn’t willing to join in. I needed more time to think about it. After finishing that paragraph, I’d title this post “Rationality Rules is a ‘Violent’ Transphobe” and punch the Publish button.

But now that I’ve finished gathering my sources and writing this post, I *have* had time to think about it. I cannot find a good reason to reject the violence-as-intentional-rights-violation definition, in particular I cannot come up with a superior alternative. Rationality Rules argues that the rights of some transgender people should be restricted, via special pleading. As I point out at that link, Stephen Woodford is aware of the argument from human rights, so he cannot claim his restriction is being done out of ignorance. That gives us proof of intent.

So no quote marks are necessary: I too believe Rationality Rules is a violent transphobe, for the definitions and reasons above.

Consider this scenario.

ME: “Hey, thanks for coming over! If you’re thirsty, I’ve got your choice of Pepsi, A&W Root Beer, and Mountain Dew.”

YOU: “I’d prefer Mountain Dew, but I’ll take anything.”

ME: “Gotcha, I’ll be right back!”

ME: [leaves, then returns with a Pepsi]

YOU: “Thanks. Too bad you ran out of Mountain Dew.”

ME: “Oh no, I’ve got tonnes.”

YOU: “… but it was too tough to reach, right?”

ME: “No, the Dew was right next to the Pepsi.”

Have I done anything wrong here? No, at least technically. You said you were fine with any soft drink, and I gave you a soft drink. At the same time, though, you expressed a preference for Mountain Dew over the other choices. I could be forgiven for ignoring your preference if I wasn’t able to fulfill it, or doing so would have been inconvenient for me, but in this fictional scenario both of those were off the table. My choice to hand you a Pepsi instead of a Mountain Dew reveals something about me, most likely that I think other people should prefer Pepsi over other soft drinks. You could point to the ordering of my list as further evidence: I’d be more likely to list my preferred option first, as it would be more prominent in my mind than the other choices, and then rattle off others as they came to me. This isn’t strong evidence, but you’d be justified in suspecting my motives as you enjoyed your Pepsi.

robbe de Boer: Hey EOT, I had a question not related to this video, what pronouns do you prefer? I couldn’t find anything real quick.

EssenceOfThought: They. But I’m fine with both she and he as well. It’s all on the channel description/Facebook page ‘About’ section. 😛

By merely existing, EssenceOfThought has set up a very similar situation. They have a pronoun preference, and when listing their pronoun choices put “he” last, but also say they aren’t offended if you pick another reasonable one. The difficulty in typing two extra letters is practically zero, and “they” as a singular pronoun has been in the English language for six centuries, so there isn’t any obstacle to its use beyond your hang-ups. Hell, even the guy who became famous for refusing to use transgender people’s pronouns is perfectly capable of using singular “they.”

Jordan Peterson: I don’t recognize another person’s right to decide what words I’m going to use, especially when the wordstheywant me to use, first of all, are non-standard elements of the English language and they are constructs of a small coterie of ideologically motivated people.Theymight have a point but I’m not going to saytheirwords forthem.

So if we encounter someone calling EssenceOfThought “he,” we’re justified in raising an eyebrow. While they’re not technically in the wrong, the use of “he” is suggestive that they’d overrule someone’s pronoun preference if they thought they could get away with it.

Steve McRae: Essence of thought bullies Rachel Oates and demonstrateshimselfto probably one of the worse humans ever to be on Twitter or YouTube.[12:00]

Noel Plum: The truth of the matter is, is that effectively in saying whathe’ssaying, Essence of Thought has said, to the majority of people who have sided withhim, you are transphobes. Your position is, transphobic unless you adopt this position that anyone who identifies as a woman gets to compete in this category, then you are holding a transphobic position. Andhismasterstroke is that it seems to have done the trick, and nobody’s arguing – none ofhissupporters are arguing withhim.

Rachel Oates: In regards to Essence of Thought calling the police and claiming to be the ‘only one’ actually helping me. It didn’t.Hecalled the police. Who apparently turned up at my old flat, broke the door down and then they wasted hours and many resources trying to find me. Meanwhile, I was at home with my friends around me, having all my Youtube friends send me love and support and check in on me. My family phoned me. People were there for me, helping.

The analogy isn’t a perfect fit. In there, I asked for your choice and received it. EssenceOfThought will mention their preferred pronoun if asked, but does not put it on blast nor do they bother to correct people who don’t pick their preference. Ignorance is more of an option than the analogy presents.

Provided, of course, these people were ignorant. Some of them claim to care about transgender people, though. They should know not to screw up someone’s pronouns, and thus be willing to do a little extra legwork to get things right. Even if they only use “he” because their friends and peers do so, that means their social circle is overwhelmingly dominated by transphobic people or people with a high tolerance to transphobia.

Rachel Oates: Also, I’m really sorry if I got EoT’s pronouns wrong in this thread – I’ve heard different things about which pronouns they prefer & may have slipped up here.

A good way to rule out ignorance is to correct them on EoT’s preferred pronouns. If that person responds with something like this…

Noel Plum: If he drops the “he” then I will drop it too. As it stands he accepts he, she or they. I couldn’t give a toss which he “prefers” as i don’t like him.

… then you’re pretty justified in believing the “cloaked transphobia” hypothesis.

“He” is also a strange choice given how EssenceOfThought presents. You see someone with no facial hair and long flowing locks in front of a transgender flag, and you immediately jump to “he?” C’mon, even Rationality Rules splits the baby and uses “she.” These people didn’t settle on “he” by accident, their choice reveals something about their internal opinion of transgender people, their peer group, or their ignorance.

And it isn’t very refreshing.

Apparently I know the solar system very well?

I attended a lecture on Carl Sagan, hosted by the Atheist Society of Calgary, and part of the event was a trivia challenge. While I wasn’t the only person at my table offering answers, my answers seemed to be the ones most consistently endorsed by the group. Assisted by some technical issues, our team wound up with a massive lead over the second-place finisher. The organizer from ASC surprised us all by saying everyone at our table could pick up a free T-shirt. I wasn’t terribly keen on wearing their logo, but I wandered over to the merch table anyway.

Sitting among the other designs was one that stopped me cold.

I was browsing YouTube videos on PyMC3, as one naturally does, when I happened to stumble on this gem.

Tech has spent millions of dollars in efforts to diversify workplaces. Despite this, it seems after each spell of progress, a series of retrograde events ensue. Anti-diversity manifestos, backlash to assertive hiring, and sexual misconduct scandals crop up every few months, sucking the air from every board room. This will be a digest of research, recent events, and pointers on women in STEM.

Lorena A. Barba *really* knows her stuff; the entire talk is a rapid-fire accounting of claims and counterclaims, aimed to directly appeal to the male techbros who need to hear it. There was a lot of new material in there, for me at least. I thought the only well-described matriarchies came from the African continent, but it turns out the Algonquin also fit that bill. Some digging turns up a rich mix of gender roles within First Nations peoples, most notably the Iroquois and Hopi. I was also depressed to hear that the R data analysis community is better at dealing with sexual harassment than the skeptic/atheist community.

But what really grabbed my ears was the section on gender quotas. I’ve long been a fan of them on logical grounds: if we truly believe the sexes are equal, then if we see unequal representation we know discrimination is happening. By forcing equality, we greatly reduce network effects where one gender can team up against the other. Worried about an increase in mediocrity? At worst that’s a temporary thing that disappears once the disadvantaged sex gets more experience, and at best the overall quality will actually go up. The research on quotas has advanced quite a bit since that old Skepchick post. Emphasis mine.

In 1993, Sweden’s Social Democratic Party centrally adopted a gender quota and imposed it on all the local branches of that party (…). Although their primary aim was to improve the representation of women, proponents of the quota observed that

the reform had an impact on the competence of men. Inger Segelström (the chair of Social Democratic Women in Sweden (S-Kvinnor), 1995–2003) made this point succinctly in a personal communication:At the time, our party’s quota policy of mandatory alternation of male and female names on all party lists became informally known as

the crisis of the mediocre man…We study the selection of municipal politicians in Sweden with regard to their competence, both theoretically and empirically. Moreover, we exploit the Social Democratic quota as a shock to municipal politics and ask how it altered the competence of that party’s elected politicians, men as well as women, and leaders as well as followers.

Besley, Timothy. “Gender Quotas and the Crisis of the Mediocre Man: Theory and Evidence from Sweden.” THE AMERICAN ECONOMIC REVIEW 107, no. 8 (2017): 39.

We can explain this with the benefit of hindsight: if men can rely on the “old boy’s network” to keep them in power, they can afford to slack off. If other sexes cannot, they have to fight to earn their place. These are all social effects, though; if no sex holds a monopoly on operational competence in reality, the net result is a handful of brilliant women among a sea of iffy men. Gender quotas severely limit the social effects, effectively kicking out the mediocre men to make way for average women, and thus *increase* the average competence.

As tidy as that picture is, it’s wrong in one crucial detail. Emphasis again mine.

These estimates show that the overall effect mainly reflects an improvement in the selection of men. The coefficient in column 4 means that a 10-percentage-point larger quota bite (just below the cross-sectional average for all municipalities) raised the proportion of competent men by 4.4 percentage points. Given an average of 50 percent competent politicians in the average municipality (by definition, from the normalization),

this corresponds to a 9 percent increase in the share of competent men.For women, we obtain a negative coefficient in the regression specification without municipality trends, but a positive coefficient with trends. In neither case, however, is the estimate significantly different from zero, suggesting that

the quota neither raised nor cut the share of competent women. This is interesting in view of the meritocratic critique of gender quotas, namely that raising the share of women through a quota must necessarily come at the price of lower competence among women.

Increasing the number of women *does not also increase the number of incompetent women*. When you introduce a quota, apparently, *everyone* works harder to justify being there. The only people truly hurt by gender quotas are mediocre men who rely on the Peter Principle.

Alas, if that YouTube like ratio is any indication, there’s a lot of them out there.

I still can’t believe this post exists, given its humble beginnings.

The “women’s category” is, in my opinion, poorly named given our current climate, and so I’d elect a name more along the lines of the “Under 5 nmol/l category” (as in, under 5 nanomoles of testosterone per litre), but make no mistake about it, the “woman’s category” is not based on gender or identity, or even genitalia or chromosomes… it’s based on hormone levels and the absence of male puberty.

The above comment wasn’t in Rationality Rules’ latest transphobic video, it was just a casual aside by RR himself in the YouTube comment section. He’s obiquely doubled-down via Twitter (hat tip to Essence of Thought):

Of course, just as I support trans men competing in all “men’s categories” (poorly named), women who have not experienced male puberty competing in all women’s sport (also poorly named) and trans women who have experienced male puberty competing in long-distance running.

To further clarify, I think that we must rename our categories according to what they’re actually based on. It’s not right to have a “women’s category” and yet say to some trans women (who are women!) that they can’t compete within it; it should be renamed.

The proposal itched away at me, though, because I knew it was testable.

There is a need to clarify hormone profiles that may be expected to occur after competition when antidoping tests are usually made. In this study, we report on the hormonal profile of 693 elite athletes, sampled within 2 h of a national or international competitive event. These elite athletes are a subset of the cross-sectional study that was a component of the GH-2000 research project aimed at developing a test to detect abuse with growth hormone.

Healy, Marie-Louise, et al. “Endocrine profiles in 693 elite athletes in the postcompetition setting.” Clinical endocrinology 81.2 (2014): 294-305.

The GH-2000 project had already done the hard work of collecting and analyzing blood samples from athletes, so checking RR’s proposal was no tougher than running some numbers. There’s all sorts of ethical guidelines around sharing medical info, but fortunately there’s an easy shortcut: ask one of the scientists involved to run the numbers for me, and report back the results. Aggregate data is much more resistant to de-anonymization, so the ethical concerns are greatly reduced. The catch, of course, is that I’d have to find a friendly researcher with access to that dataset. About a month ago, I fired off some emails and hoped for the best.

I wound up much, much better than the best. *I got full access to the dataset!!* You don’t get handed an incredible gift like this and merely use it for a blog post. In my spare time, I’m flexing my Bayesian muscles to do a re-analysis of the above paper, while also looking for observations the original authors may have missed. Alas, that means my slow posting schedule is about to crawl.

But in the meantime, we have a question to answer.

(Click here to show the code)

```
import numpy as np
import pandas as pd
dataset = pd.read_csv('dataset.minimal.tsv',sep='\t')
mask_afab = dataset['Gender']==2
mask_amab = dataset['Gender']==1
print( "{:24} = {}".format("Total Assigned-female Athletes", np.sum(mask_afab)) )
print( "{:24} = {:.2f} cm".format(" Height, Mean", np.mean( dataset['height'][mask_afab] )) )
print( "{:24} = {:.2f} cm".format(" Height, Std.Dev", np.std( dataset['height'][mask_afab] )) )
print( "{:24} = {:.2f} kg".format(" Weight, Mean", np.mean( dataset['weight'][mask_afab] )) )
print( "{:24} = {:.2f} kg".format(" Weight, Std.Dev", np.std( dataset['weight'][mask_afab] )) )
print( "{:24} = {:.2f} kg".format(" Body Fat, Mean", np.mean( (dataset['weight']*dataset['body-fat']*.01)[mask_afab] )) )
print( "{:24} = {:.2f} kg".format(" Body Fat, Std.Dev", np.std( (dataset['weight']*dataset['body-fat']*.01)[mask_afab] )) )
print( "{:24} = {:.2f} nmol/L".format(" Testosterone, Mean", np.mean( dataset['Testo'][mask_afab] )) )
print( "{:24} = {:.2f} nmol/L".format(" Testosterone, Std.Dev", np.std( dataset['Testo'][mask_afab] )) )
print( "{:24} = {:.2f} nmol/L".format(" Testosterone, Max", np.max( dataset['Testo'][mask_afab] )) )
print( "{:24} = {:.2f} nmol/L".format(" Testosterone, Min", np.min( dataset['Testo'][mask_afab] )) )
print()
print( "{:24} = {}".format("Total Assigned-male Athletes", np.sum(mask_amab) ) )
print( "{:24} = {:.2f} cm".format(" Height, Mean", np.mean( dataset['height'][mask_amab] )) )
print( "{:24} = {:.2f} cm".format(" Height, Std.Dev", np.std( dataset['height'][mask_amab] )) )
print( "{:24} = {:.2f} kg".format(" Weight, Mean", np.mean( dataset['weight'][mask_amab] )) )
print( "{:24} = {:.2f} kg".format(" Weight, Std.Dev", np.std( dataset['weight'][mask_amab] )) )
print( "{:24} = {:.2f} kg".format(" Body Fat, Mean", np.mean( (dataset['weight']*dataset['body-fat']*.01)[mask_amab] )) )
print( "{:24} = {:.2f} kg".format(" Body Fat, Std.Dev", np.std( (dataset['weight']*dataset['body-fat']*.01)[mask_amab] )) )
print( "{:24} = {:.2f} nmol/L".format(" Testosterone, Mean", np.mean( dataset['Testo'][mask_amab] )) )
print( "{:24} = {:.2f} nmol/L".format(" Testosterone, Std.Dev", np.std( dataset['Testo'][mask_amab] )) )
print( "{:24} = {:.2f} nmol/L".format(" Testosterone, Max", np.max( dataset['Testo'][mask_amab] )) )
print( "{:24} = {:.2f} nmol/L".format(" Testosterone, Min", np.min( dataset['Testo'][mask_amab] )) )
```

Total Assigned-female Athletes = 239 Height, Mean = 171.61 cm Height, Std.Dev = 7.12 cm Weight, Mean = 64.27 kg Weight, Std.Dev = 9.12 kg Body Fat, Mean = 13.19 kg Body Fat, Std.Dev = 3.85 kg Testosterone, Mean = 2.68 nmol/L Testosterone, Std.Dev = 4.33 nmol/L Testosterone, Max = 31.90 nmol/L Testosterone, Min = 0.00 nmol/L Total Assigned-male Athletes = 454 Height, Mean = 182.72 cm Height, Std.Dev = 8.48 cm Weight, Mean = 80.65 kg Weight, Std.Dev = 12.62 kg Body Fat, Mean = 8.89 kg Body Fat, Std.Dev = 7.20 kg Testosterone, Mean = 14.59 nmol/L Testosterone, Std.Dev = 6.66 nmol/L Testosterone, Max = 41.00 nmol/L Testosterone, Min = 0.80 nmol/L

The first step is to get a basic grasp on what’s there, via some crude descriptive statistics. It’s also useful to compare these with the original paper, to make sure I’m interpreting the data correctly. Excusing some minor differences in rounding, the above numbers match the paper.

The only thing that stands out from the above, to me, is the serum levels of testosterone. At least one source says the mean of these assigned-female athletes is higher than the normal range for their non-athletic cohorts. Part of that may simply be because we don’t have a good idea of what the normal range is, so it’s not uncommon for each lab to have their own definition of “normal.” This is even worse for those assigned female, since their testosterone levels are poorly studied; note that my previous link collected the data of over a million “men,” but doesn’t mention “women” once. Factor in inaccurate test results and other complicating factors, and “normal” is quite poorly-defined.

Still, Rationality Rules is either convinced those complications are irrelevant, or ignorant of them. And, to be fair, that 5nmol/L line implicitly sweeps a lot of them under the rug. Let’s carry on, then, and look for invalid data. “Invalid” covers everything from missing data, to impossible data, and maybe even data we think might be made inaccurate due to measurement error. I consider a concentration of zero testosterone as invalid, even though it may technically be possible.

(Click here to show the code)

```
t_number = dataset['Testo'] >= 0
t_max = np.max(dataset['Testo'][t_number])
t_min = np.min(dataset['Testo'][t_number])
t_valid = (dataset['Testo'] > 0.5) & np.isfinite(dataset['Testo'])
print( "{:52} = {}".format("Total Assigned-male Athletes w/ T levels >= 0", np.sum(mask_amab & t_number) ) )
print( "{:52} = {}".format(" w/ T levels <= 0.5", np.sum(mask_amab & (dataset['Testo']<=0.5)) ) )
print( "{:52} = {}".format(" w/ T levels == 0", np.sum(mask_amab & (dataset['Testo']==0)) ) )
print( "{:52} = {}".format(" w/ missing T levels", np.sum(mask_amab & np.isnan(dataset['Testo'])) ) )
print( "{:52} = {}".format(" that I consider valid", np.sum(mask_amab & t_valid)) )
print()
print( "{:52} = {}".format("Total Assigned-female Athletes w/ T levels >= 0", np.sum(mask_afab & t_number)) )
print( "{:52} = {}".format(" w/ T levels <= 0.5", np.sum(mask_afab & (dataset['Testo']<=0.5)) ) )
print( "{:52} = {}".format(" w/ T levels == 0", np.sum(mask_afab & (dataset['Testo']==0)) ) )
print( "{:52} = {}".format(" w/ missing T levels", np.sum(mask_afab & np.isnan(dataset['Testo'])) ) )
print( "{:52} = {}".format(" that I consider valid", np.sum(mask_afab & t_valid)) )
```

Total Assigned-male Athletes w/ T levels >= 0 = 446 w/ T levels <= 0.5 = 0 w/ T levels == 0 = 0 w/ missing T levels = 8 that I consider valid = 446 Total Assigned-female Athletes w/ T levels >= 0 = 234 w/ T levels <= 0.5 = 5 w/ T levels == 0 = 1 w/ missing T levels = 5 that I consider valid = 229

Fortunately for us, the losses are pretty small. 229 datapoints is a healthy sample size, so we can afford to be liberal about what we toss out. Next up, it would be handy to see the data in chart form.

(Click here to show the code)

```
# %matplotlib notebook # makes the plots interactive, but only one can be active
%matplotlib inline
import matplotlib.pyplot as pp
pp.rcParams['figure.dpi'] = 96 # MUST SET THIS FIRST
pp.rcParams['figure.figsize'] = [9.5, 6]
bins = 9
pp.hist( np.log(dataset['Testo'][mask_afab & t_valid]), bins, density=1, facecolor='black', alpha=0.2)
pp.hist( np.log(dataset['Testo'][mask_amab & t_valid]), bins, density=1, facecolor='green', alpha=0.2)
pp.legend(['aFab','aMab'], loc=0)
pp.title('Testosterone, elite athletes')
pp.xlabel('nmol/L')
pp.xticks(np.linspace(-2,4,9), ["{:.1f}".format(np.exp(x)) for x in np.linspace(-2,4,9)])
pp.yticks([])
pp.axvline(np.log(0.5),0,1)
pp.axvline(np.log(5),0,1)
# Source: https://www.exeterlaboratory.com/test/testosterone/
pp.fill( np.log([29,8.6,8.6,29]), [1.2,1.2,0,0], facecolor='green', alpha=0.05 )
pp.fill( np.log([1.68,.3,.3,1.68]), [1.2,1.2,0,0], facecolor='black', alpha=0.05 )
pp.show()
```

I've put vertical lines at both the 0.5 and 5 nmol/L cutoffs. There's a big difference between categories, but we can see clouds on the horizon: a substantial number of assigned-female athletes have greater than 5 nmol/L of testosterone in their bloodstream, while a decent number of assigned-male athletes have less. How many?

(Click here to show the code)

```
mask_gt_5nmol = t_valid & (dataset['Testo'] > 5)
mask_lt_5nmol = t_valid & (dataset['Testo'] < 5)
mask_eq_5nmol = t_valid & (dataset['Testo'] == 5)
print("Segregating Athletes by Testosterone")
table = {"Concentration":["> 5nmol/L","< 5nmol/L","= 5nmol/L"],
"aFab":[sum(mask_gt_5nmol & mask_afab),sum(mask_lt_5nmol & mask_afab),sum(mask_eq_5nmol & mask_afab)],
"aMab":[sum(mask_gt_5nmol & mask_amab), sum(mask_lt_5nmol & mask_amab), sum(mask_eq_5nmol & mask_amab)]}
print(pd.DataFrame(table).to_string(index=False))
print()
print("{:.1f}% of assigned-female athletes have > 5nmol/L".format(100.*sum(mask_gt_5nmol & mask_afab)/sum(t_valid & mask_afab)))
print("{:.1f}% of assigned-male athletes have < 5nmol/L".format(100.*sum(mask_lt_5nmol & mask_amab)/sum(t_valid & mask_amab)))
print("{:.1f}% of athletes with > 5nmol/L are assigned-female".format(100.*sum(mask_gt_5nmol & mask_afab)/sum(mask_gt_5nmol)))
print("{:.1f}% of athletes with < 5nmol/L are assigned-male".format(100.*sum(mask_lt_5nmol & mask_amab)/sum(mask_lt_5nmol)))
```

Segregating Athletes by Testosterone Concentration aFab aMab > 5nmol/L 19 417 < 5nmol/L 210 26 = 5nmol/L 0 3 8.3% of assigned-female athletes have > 5nmol/L 5.8% of assigned-male athletes have < 5nmol/L 4.4% of athletes with > 5nmol/L are assigned-female 11.0% of athletes with < 5nmol/L are assigned-male

Looks like anywhere from 6-8% of athletes have testosterone levels that cross Rationality Rules' line. For comparison, maybe 1-2% of the general public has some level of gender dysphoria, though estimating exact figures is hard in the face of widespread discrimination and poor sex-ed in schools. Even that number is misleading, as the number of transgender athletes is substantially lower than 1-2% of the athletic population. The share of transgender athletes is irrelevant to this dataset anyway, as it was collected between 1996 and 1999, when no sporting agency had policies that allowed transgender athletes to openly compete.

That 6-8%, in other words, is entirely cisgender. This echoes one of Essence Of Thought's arguments: RR's 5nmol/L policy has far more impact on cis athletes than trans athletes, which could have catastrophic side-effects. *Could* is the operative word, though, because as of now we don't know anything about these athletes. Do >5nmol/L assigned-female athletes have bodies more like >5nmol/L assigned-male athletes than <5nmol/L assigned-female athletes? If so, then there's no problem. Equivalent body types are competing against each other, and outcomes are as fair as could be reasonably expected.

What, then, counts as an "equivalent" body type when it comes to sport?

One reasonable measure of equivalence is height. It's one of the stronger sex differences, and height is also correlated with longer limbs and greater leverage. Whether that's relevant to sports is debatable, but height and correlated attributes dominate Rationality Rules' list.

[19:07] In some events - such as long-distance running, in which hemoglobin and slow-twitch muscle fibers are vital - I think there's a strong argument to say no, [transgender women who transitioned after puberty] don't have an unfair advantage, as the primary attributes are sufficiently mitigated. But in most events, and especially those in which

height, width, hip size, limb length, muscle mass, and muscle fiber type are the primary attributes - such as weightlifting, sprinting, hammer throw, javelin, netball, boxing, karate, basketball, rugby, judo, rowing, hockey, and many more - my answer is yes, most do have an unfair advantage.

Fortunately for both of us, most athletes in the dataset have a "valid" height, which I define as being at least 30cm tall.

(Click here to show the code)

```
height_valid = dataset['height'] > 30
print("Out of {:3} athletes, {} have valid height data.".format(len(height_valid), sum(height_valid)) )
bins = 9
pp.hist( dataset['height'][mask_afab & height_valid], bins, density=1, facecolor='black', alpha=0.2)
pp.hist( dataset['height'][mask_amab & height_valid], bins, density=1, facecolor='green', alpha=0.2)
pp.legend(['aFab','aMab'], loc=0)
pp.title('Height, elite athletes')
pp.xlim([145,215])
pp.xlabel("cm")
pp.yticks([])
# source: https://ourworldindata.org/human-height, Germany, 1976
pp.axvline(166.3,0,1, color='k', alpha=0.2)
pp.axvline(np.mean(dataset['height'][mask_afab & height_valid]),0,1, color='k')
pp.axvline(179.3,0,1, color='g', alpha=0.2)
pp.axvline(np.mean(dataset['height'][mask_amab & height_valid]),0,1, color='g')
pp.show()
```

Out of 693 athletes, 678 have valid height data.

The faint vertical lines are for the mean adult height of Germans born in 1976, which should be a reasonable cohort to European athletes that were active between 1996 and 1999, while the darker lines are each category's mean. Athletes seem slightly taller than the reference average, but only by 2-5cm. The amount of overlap is also surprising, given that height is supposed to be a major sex difference. We actually saw *less* overlap with testosterone! Finally, the height distribution isn't quite Gaussian, there's a subtle bias towards the taller end of the spectrum.

Height is a pretty crude metric, though. You could pair any athlete with a non-athlete of the same height, and there's no way the latter would perform as well as the former. A better measure of sporting ability would be muscle mass. We shouldn't use the absolute mass, though: bigger bodies have more mass and need more force to accelerate as smaller bodies do, so height and muscle mass are correlated. We need some sort of dimensionless scaling factor which compensates.

And we have one! It's called the Body Mass Index, or BMI.

$$ BMI = \frac w {h^2}, $$

where \(w\) is a person's mass in kilograms, and \(h\) is a person's height in metres. Unfortunately, BMI is quite problematic. Partly that's because it is a crude measure of obesity. But part of that is because there are two types of tissue which can greatly vary, body fat and muscle, yet both contribute equally towards BMI.

That's all fixable. For one, some of the athletes in this dataset had their body fat measured. We can subtract that mass off, so their weight consists of tissues that are strongly correlated with height plus one that is fudgable: muscle mass. For two, we're not assessing these individual's health, we only want a dimensionless measure of muscle mass relative to height. For three, we're not comparing these individuals to the general public, so we're not restricted to using the general BMI formula. We can use something more accurate.

The oddity is the appearance of that exponent 2, though our world is three-dimensional. You might think that the exponent should simply be 3, but that doesn't match the data at all. It has been known for a long time that people don't scale in a perfectly linear fashion as they grow. I propose that a better approximation to the actual sizes and shapes of healthy bodies might be given by an exponent of 2.5. So here is the formula I think is worth considering as an alternative to the standard BMI:

$$ BMI' = 1.3 \frac w {h^{2.5}} $$

I can easily pop body fat into Nick Trefethen's formula, and get a better measure of relative muscle mass,

$$ \overline{BMI} = 1.3 \frac{ w - bf }{h^{2.5}}, $$

where \(bf\) is total body fat in kilograms. Individuals with excess muscle mass, relative to what we expect for their height, will have a high \(\overline{BMI}\), and vice-versa. And as we saw earlier, muscle mass is another of Rationality Rules' determinants of sporting performance.

Time for more number crunching.

(Click here to show the code)

```
BMI_adj = 1.3*(dataset['weight']*(100. - dataset['body-fat']))*0.01/((dataset['height']*.01)**(2.5))
BMI_adj_valid = BMI_adj > 1
print( "Out of {:3} athletes, {} have valid adjusted BMIs.".format(len(BMI_adj_valid), sum(BMI_adj_valid)) )
print( " {} have valid weights.".format(sum(dataset['weight'] > 10)) )
print( " {} have valid body fat percentages.".format(sum(dataset['body-fat'] >= 0)) )
print()
print( "{:24} = {}".format("Total Assigned-female Athletes", np.sum(mask_afab)) )
print( "{:24} = {}".format(" total with valid adjusted BMI", np.sum(mask_afab & BMI_adj_valid)) )
print( "{:24} = {:.2f}".format(" adjusted BMI, Mean", np.mean( BMI_adj[mask_afab & BMI_adj_valid] )) )
print( "{:24} = {:.2f}".format(" adjusted BMI, Std.Dev", np.std( BMI_adj[mask_afab & BMI_adj_valid] )) )
print( "{:24} = {:.2f}".format(" adjusted BMI, Median", np.median( BMI_adj[mask_afab & BMI_adj_valid] )) )
print()
print( "{:24} = {}".format("Total Assigned-male Athletes", np.sum(mask_amab)) )
print( "{:24} = {}".format(" total with valid adjusted BMI", np.sum(mask_amab & BMI_adj_valid)) )
print( "{:24} = {:.2f}".format(" adjusted BMI, Mean", np.mean( BMI_adj[mask_amab & BMI_adj_valid] )) )
print( "{:24} = {:.2f}".format(" adjusted BMI, Std.Dev", np.std( BMI_adj[mask_amab & BMI_adj_valid] )) )
print( "{:24} = {:.2f}".format(" adjusted BMI, Median", np.median( BMI_adj[mask_amab & BMI_adj_valid] )) )
```

Out of 693 athletes, 227 have valid adjusted BMIs. 663 have valid weights. 241 have valid body fat percentages. Total Assigned-female Athletes = 239 total with valid adjusted BMI = 86 adjusted BMI, Mean = 16.98 adjusted BMI, Std.Dev = 1.21 adjusted BMI, Median = 16.96 Total Assigned-male Athletes = 454 total with valid adjusted BMI = 141 adjusted BMI, Mean = 20.56 adjusted BMI, Std.Dev = 1.88 adjusted BMI, Median = 20.28

The bad news is that most of this dataset lacks any information on body fat, which really cuts into our sample size. The good news is that we've still got enough to carry on. It also looks like there's a strong sex difference, and the distribution is pretty clustered. Still, a chart would help clarify the latter point.

(Click here to show the code)

```
bins = 9
pp.hist( np.log(BMI_adj)[mask_afab & BMI_adj_valid], bins, density=1, facecolor='black', alpha=0.2)
pp.hist( np.log(BMI_adj)[mask_amab & BMI_adj_valid], bins, density=1, facecolor='green', alpha=0.2)
pp.legend(['aFab','aMab'], loc=0)
pp.title('Adjusted BMI, elite athletes')
pp.xticks(np.linspace(np.log(14),np.log(28),8),
["{:.1f}".format(np.exp(x)) for x in np.linspace(np.log(14),np.log(28),8)])
pp.yticks([])
pp.axvline(np.log(np.mean( BMI_adj[mask_afab & BMI_adj_valid] )),0,1, color='k')
pp.axvline(np.log(np.mean( BMI_adj[mask_amab & BMI_adj_valid] )),0,1, color='g')
pp.show()
```

Whoops! There's more overlap and skew than I thought. Even in logspace, the results don't look Gaussian. We'll have to remember that for the next step.

Just looking at charts isn't going to solve this question, we need to do some sort of hypothesis testing. Fortunately, all the pieces I need are here. We've got our hypothesis, for instance:

Athletes with exceptional testosterone levels are more like athletes of the same sex but with typical testosterone levels, than they are of other athletes with a different sex but similar testosterone levels.

If you know me, you know that I'm all about the Bayes, and that gives us our methodology.

- Fit a model to a specific metric for assigned-female athletes with less than 5nmol/L of serum testosterone.
- Fit a model to a specific metric for assigned-male athletes with more than 5nmol/L of serum testosterone.
- Apply the first model to the test group, calculating the overall likelihood.
- Apply the second model to the test group, calculating the overall likelihood.
- Sample the probability distribution of the Bayes Factor.

"Metric" is one of height or \(\overline{BMI}\), while "test group" is one of assigned-female athletes with >5nmol/L of serum testosterone or assigned-male athletes with <5nmol/L of serum testosterone. The Bayes Factor is simply

$$ \text{Bayes Factor} = \frac{ p(E \mid H_1) \cdot p(H_1) }{ p(E \mid H_2) \cdot p(H_2) } = \frac{ p(H_1 \mid E) }{ p(H_2 \mid E) }, $$

which means we need two hypotheses, not one. Fortunately, I've phrased the hypothesis to make it easy to negate: athletes with exceptional testosterone levels are *less* like athletes of the same sex but with typical testosterone levels, than they are of other athletes with a different sex but similar testosterone levels. We'll call this new hypothesis \(H_2\), and the original \(H_1\). Bayes factors greater than 1 mean \(H_1\) is more likely than \(H_2\), and vice-versa.

Calculating all that would be easy if I was using Stan or PyMC3, but I ran into problems translating the former's probability distributions into charts, and I don't have any experience with the latter. My next choice, emcee, forces me to manually convolve two posterior distributions. Annoying, but not difficult.

That just leaves one thing left: what models are we going to use? The obvious choice for height is the Gaussian distribution, as from previous research we know it's a great model.

(Click here to show the code)

```
import emcee
nwalkers, nsamples = 128, 6
models = dict()
import os
import scipy.stats as spst
ndim = 2
def lnLike_gaussian( theta, x ):
mu, sigma = theta
return np.sum( spst.norm( mu, sigma ).logpdf( x ) )
def lnPrior_gaussian( theta ):
mu, sigma = theta
if sigma <= 0: # standard deviation must be positive
return -np.inf
return -2*np.log(sigma) # favor lower standard deviations
def lnProb_gaussian( theta, x ):
temp = lnPrior_gaussian( theta )
if temp == -np.inf:
return temp
return temp + lnLike_gaussian( theta, x )
x = dataset['height'][mask_lt_5nmol & height_valid & mask_afab]
pos = [np.array([150,15]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)] # rough estimate
print("Fitting the height of lT aFab athletes to a Gaussian distribution ...")
lnprob = None
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnProb_gaussian, threads=os.cpu_count(), args=[ x ] )
model_mean = np.mean( pos, axis=0 )
print("{:6}: ({:5f}) mu={:7f}, sigma={:7f}".format(0, lnProb_gaussian( model_mean, x ), *model_mean))
for it in range(nsamples):
global pos, lnprob
pos, lnprob, _ = sampler.run_mcmc( pos, 64, storechain=False )
model_mean = np.mean( pos, axis=0 ) # fairer than going with the maximal likelihood
print("{:6}: ({:5f}) mu={:7f}, sigma={:7f}".format((it+1) * 64, lnProb_gaussian( model_mean, x ), *model_mean))
best = np.argmax( lnprob )
print("{:>6}: ({:5f}) mu={:7f}, sigma={:7f}".format("ML", lnProb_gaussian( pos[best], x ), *pos[best]))
model_median = np.median( pos, axis=0 )
print("{:>6}: ({:5f}) mu={:7f}, sigma={:7f}".format("median", lnProb_gaussian( model_median, x ), *model_median))
models["lT_aFab_height"] = pos # store the posterior for later use
```

Fitting the height of lT aFab athletes to a Gaussian distribution ... 0: (-980.322471) mu=150.000819, sigma=15.000177 64: (-710.417497) mu=169.639051, sigma=8.579088 128: (-700.539260) mu=171.107358, sigma=7.138832 192: (-700.535241) mu=171.154151, sigma=7.133279 256: (-700.540692) mu=171.152701, sigma=7.145515 320: (-700.552831) mu=171.139668, sigma=7.166857 384: (-700.530969) mu=171.086422, sigma=7.094077 ML: (-700.525284) mu=171.155240, sigma=7.085777 median: (-700.525487) mu=171.134614, sigma=7.070993

Alas, emcee also lacks a good way to assess model fitness. One crude metric is look at the progression of the mean fitness; if it grows and then stabilizes around a specific value, as it does here, we've converged on something. Another is to compare the mean, median, and maximal likelihood of the posterior; if they're about equally likely, we've got a fuzzy caterpillar. Again, that's also true here.

As we just saw, though, charts are a better judge of fitness than a handful of numbers.

(Click here to show the code)

```
def plotWrapper( func, theta, x ):
retVal = list()
for data in x:
retVal.append( np.exp( func( theta, data ) ) )
return retVal
minVal = min(dataset['height'][height_valid])
maxVal = max(dataset['height'][height_valid])
x = np.linspace(minVal,maxVal,255)
bins = 9
pp.hist( dataset['height'][mask_afab & height_valid], bins, density=1, facecolor='black', alpha=0.2)
pp.legend(['lT aFab'], loc=0)
pp.plot( x, plotWrapper(lnLike_gaussian, np.mean(models['lT_aFab_height'],axis=0), x ), 'k' )
pp.title('Height, elite athletes')
pp.xlim([145,215])
pp.xlabel("cm")
pp.yticks([])
# source: https://ourworldindata.org/human-height, Germany, 1976
pp.axvline(166.3,0,1, color='k', alpha=0.2)
pp.axvline(np.mean(dataset['height'][mask_afab & height_valid]),0,1, color='k')
pp.show()
```

If you were wondering why I didn't make much of a fuss out of the asymmetry in the height distribution, it's because I've already seen this graph. A good fit isn't necessarily *the best* though, and I might be able to get a closer match by incorporating the sport each athlete played.

(Click here to show the code)

```
sport_names = {
1: 'Power lifting',
2: 'Basketball',
3: 'Football',
4: 'Swimming',
5: 'Marathon',
6: 'Canoeing',
7: 'Rowing',
8: 'Cross-country skiing',
9: 'Alpine skiing',
10: 'Weight lifting',
11: 'Judo',
12: 'Bandy',
13: 'Ice Hockey',
14: 'Handball',
15: 'Track and field'}
print("{:^48}".format("Assigned-female Athletes"))
print("{:^24} {:^23}".format("sport","below/above 171cm"))
for sport in pd.Categorical( dataset['sport'] ).categories:
below = sum(dataset['sport'][mask_afab & height_valid & (dataset['height'] < 171)] == sport) above = sum(dataset['sport'][mask_afab & height_valid & (dataset['height'] >= 171)] == sport)
print("{:>24}: {:2} /{:2}".format(sport_names[sport], below, above))
```

Assigned-female Athletes sport below/above 171cm Power lifting: 1 / 0 Basketball: 2 /12 Football: 0 / 0 Swimming: 41 /49 Marathon: 0 / 1 Canoeing: 1 / 0 Rowing: 9 /13 Cross-country skiing: 8 / 1 Alpine skiing: 11 / 1 Weight lifting: 7 / 0 Judo: 0 / 0 Bandy: 0 / 0 Ice Hockey: 0 / 0 Handball: 12 /17 Track and field: 22 /27

Basketball attracts tall people, unsurprisingly, while skiing seems to attract shorter people. This could be the cause of that asymmetry. It's no guarantee that I'll actually get a better fit, though, as I'm also dramatically cutting the number of datapoints to fit to. The model's uncertainty must increase as a result, and that may be enough to dilute out any increase in fitness. I'll run those numbers for the paper, but for now the Gaussian model I have is plenty good.

(Click here to show the code)

```
x = dataset['height'][mask_gt_5nmol & height_valid & mask_amab]
pos = [np.array([150,15]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)] # rough estimate
print("Fitting the height of hT aMab athletes to a Gaussian distribution ...")
lnprob = None
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnProb_gaussian, threads=os.cpu_count(), args=[ x ] )
model_mean = np.mean( pos, axis=0 )
print("{:6}: ({:5f}) mu={:7f}, sigma={:7f}".format(0, lnProb_gaussian( model_mean, x ), *model_mean))
for it in range(nsamples):
global pos, lnprob
pos, lnprob, _ = sampler.run_mcmc( pos, 64, storechain=False )
model_mean = np.mean( pos, axis=0 ) # fairer than going with the maximal likelihood
print("{:6}: ({:5f}) mu={:7f}, sigma={:7f}".format((it+1) * 64, lnProb_gaussian( model_mean, x ), *model_mean))
best = np.argmax( lnprob )
print("{:>6}: ({:5f}) mu={:7f}, sigma={:7f}".format("ML", lnProb_gaussian( pos[best], x ), *pos[best]))
model_median = np.median( pos, axis=0 )
print("{:>6}: ({:5f}) mu={:7f}, sigma={:7f}".format("median", lnProb_gaussian( model_median, x ), *model_median))
models["hT_aMab_height"] = pos # store the posterior for later use
```

Fitting the height of hT aMab athletes to a Gaussian distribution ... 0: (-2503.079578) mu=150.000061, sigma=15.001179 64: (-1482.315571) mu=179.740851, sigma=10.506003 128: (-1451.789027) mu=182.615810, sigma=8.620333 192: (-1451.748336) mu=182.587979, sigma=8.550535 256: (-1451.759883) mu=182.676004, sigma=8.546410 320: (-1451.746697) mu=182.626918, sigma=8.538055 384: (-1451.747266) mu=182.580692, sigma=8.534070 ML: (-1451.746074) mu=182.591047, sigma=8.534584 median: (-1451.759295) mu=182.603231, sigma=8.481894

We get the same results when fitting the model to >5 nmol/L assigned-male athletes. The log likelihood, that number in brackets, is a lot lower for these athletes, but that number is roughly proportional to the number of samples. If we had the same degree of model fitness but doubled the number of samples, we'd expect the log likelihood to double. And, sure enough, this dataset has roughly twice as many assigned-male athletes as it does assigned-female athletes.

(Click here to show the code)

```
x = np.linspace(minVal,maxVal,255)
bins = 9
pp.hist( dataset['height'][mask_lt_5nmol & mask_afab & height_valid], bins, density=1, facecolor='black', alpha=0.2)
pp.hist( dataset['height'][mask_gt_5nmol & mask_amab & height_valid], bins, density=1, facecolor='green', alpha=0.2)
pp.legend(['lT aFab','hT aMab'], loc=0)
pp.plot( x, plotWrapper(lnLike_gaussian, np.mean(models['lT_aFab_height'],axis=0), x ), 'k' )
pp.plot( x, plotWrapper(lnLike_gaussian, np.mean(models['hT_aMab_height'],axis=0), x ), 'g' )
pp.title('Height, elite athletes')
pp.xlim([145,215])
pp.xlabel("cm")
pp.yticks([])
# source: https://ourworldindata.org/human-height, Germany, 1976
pp.axvline(166.3,0,1, color='k', alpha=0.2)
pp.axvline(np.mean(dataset['height'][mask_lt_5nmol & mask_afab & height_valid]),0,1, color='k')
pp.axvline(179.3,0,1, color='g', alpha=0.2)
pp.axvline(np.mean(dataset['height'][mask_gt_5nmol & mask_amab & height_valid]),0,1, color='g')
pp.show()
```

The updated charts are more of the same.

Unfortunately, adjusted BMI isn't nearly as tidy. I don't have any prior knowledge that would favour a particular model, so I wound up testing *five* candidates: the Gaussian, Log-Gaussian, Gamma, Weibull, and Rayleigh distributions. All but the first needed an offset parameter to get the best results, which has the same interpretation as last time.

(Click here to show the code)

```
ndim = 2
x = BMI_adj[mask_gt_5nmol & BMI_adj_valid & mask_amab]
pos = [np.array([15,5]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)] # rough estimate
print("Fitting the adjusted BMI of hT aMab athletes to a Gaussian distribution ...")
lnprob = None
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnProb_gaussian, threads=os.cpu_count(), args=[ x ] )
model_mean = np.mean( pos, axis=0 )
print("{:6}: ({:5f}) mu={:7f}, sigma={:7f}".format(0, lnProb_gaussian( model_mean, x ), *model_mean))
for it in range(nsamples):
global pos, lnprob
pos, lnprob, _ = sampler.run_mcmc( pos, 64, storechain=False )
model_mean = np.mean( pos, axis=0 ) # fairer than going with the maximal likelihood
print("{:6}: ({:5f}) mu={:7f}, sigma={:7f}".format((it+1) * 64, lnProb_gaussian( model_mean, x ), *model_mean))
best = np.argmax( lnprob )
print("{:>6}: ({:5f}) mu={:7f}, sigma={:7f}".format("ML", lnProb_gaussian( pos[best], x ), *pos[best]))
model_median = np.median( pos, axis=0 )
print("{:>6}: ({:5f}) mu={:7f}, sigma={:7f}".format("median", lnProb_gaussian( model_median, x ), *model_median))
models['hT_aMab_BMI_gaussian'] = pos
```

Fitting the adjusted BMI of hT aMab athletes to a Gaussian distribution ... 0: (-410.901047) mu=14.999563, sigma=5.000388 384: (-256.474147) mu=20.443497, sigma=1.783979 ML: (-256.461460) mu=20.452817, sigma=1.771653 median: (-256.477475) mu=20.427138, sigma=1.781139

(Click here to show the code)

```
ndim = 3
def lnLike_loggaussian( theta, x ):
mu, sigma, off = theta
x_adj = x - off
if np.any( x_adj < 0 ):
return -np.inf
return np.sum( -.5*( ((x_adj-mu)/sigma)**2 ) - np.log( x_adj*sigma ) )
def lnPrior_loggaussian( theta ):
mu, sigma, off = theta
if (mu < 0) or (sigma <= 0):
return -np.inf
if (off < 0) or (off > 25):
return -np.inf
return -2*np.log(sigma)
def lnProb_loggaussian( theta, x ):
temp = lnPrior_loggaussian( theta )
if temp == -np.inf:
return temp
return temp + lnLike_loggaussian(theta, x)
pos = [np.array([7,2,10]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)] # rough estimate
print("Fitting the adjusted BMI of hT aMab athletes to a Log-Gaussian distribution ...")
lnprob = None
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnProb_loggaussian, threads=os.cpu_count(), args=[ x ] )
model_mean = np.mean( pos, axis=0 )
print("{:6}: ({:5f}) mu={:7f}, sigma={:7f}, off={:7f}".format(0, \
lnProb_loggaussian( model_mean, x ), *model_mean))
for it in range(nsamples):
global pos, lnprob
pos, lnprob, _ = sampler.run_mcmc( pos, 64, storechain=False )
model_mean = np.mean( pos, axis=0 )
print("{:6}: ({:5f}) mu={:7f}, sigma={:7f}, off={:7f}".format((it+1) * 64, \
lnProb_loggaussian( model_mean, x ), *model_mean))
best = np.argmax( lnprob )
print("{:>5}: ({:5f}) mu={:7f}, sigma={:7f}, off={:7f}".format("ML", \
lnprob[best], *pos[best]) )
model_median = np.median( pos, axis=0 )
print("{:6}: ({:5f}) mu={:7f}, sigma={:7f}, off={:7f}".format("median", \
lnProb_loggaussian( model_median, x ), *model_median))
models['hT_aMab_BMI_loggaussian'] = pos
```

Fitting the adjusted BMI of hT aMab athletes to a Log-Gaussian distribution ... 0: (-629.141577) mu=6.999492, sigma=2.001107, off=10.000768 384: (-290.910651) mu=3.812746, sigma=1.789607, off=16.633741 ML: (-277.119315) mu=3.848383, sigma=1.818429, off=16.637382 median: (-288.278918) mu=3.795675, sigma=1.778238, off=16.637076

(Click here to show the code)

```
import scipy as sp
ndim = 3
def lnLike_gammaoffset(theta, x):
alpha, beta, off = theta
x_adj = x - off
if np.any( x_adj < 0 ):
return -np.inf
return np.sum( alpha*np.log(beta) - sp.special.loggamma(alpha) + (alpha-1)*np.log(x_adj) - beta*x_adj )
lnPrior_gammaoffset = lnPrior_loggaussian # the two are similar enough to reuse
def lnProb_gammaoffset( theta, x ):
temp = lnPrior_gammaoffset( theta )
if temp == -np.inf:
return temp
return temp + lnLike_gammaoffset(theta, x)
pos = [np.array([20,3.,10.]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
print("Fitting the adjusted BMI of hT aMab athletes to a Gamma distribution ...")
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnProb_gammaoffset, threads=os.cpu_count(), args=[ x ] )
gammaoffset_mean = np.mean( pos, axis=0 )
print("{:5}: ({:5f}) alpha={:7f}, beta={:7f}, off={:7f}".format(0,
lnProb_gammaoffset( gammaoffset_mean, x ), *gammaoffset_mean))
for it in range(nsamples):
global pos, lnprob
pos, lnprob, _ = sampler.run_mcmc( pos, 64, storechain=False )
model_mean = np.mean( pos, axis=0 )
print("{:6}: ({:5f}) alpha={:7f}, beta={:7f}, off={:7f}".format((it+1) * 64, \
lnProb_gammaoffset( model_mean, x ), *model_mean))
best = np.argmax( lnprob )
print("{:6}: ({:5f}) alpha={:7f}, beta={:7f}, off={:7f}".format("ML", \
lnprob[best], *pos[best]) )
model_median = np.median( pos, axis=0 )
print("{:6}: ({:5f}) alpha={:7f}, beta={:7f}, off={:7f}".format("median", \
lnProb_gammaoffset( model_median, x ), *model_median))
models['hT_aMab_BMI_gamma'] = pos
```

Fitting the adjusted BMI of hT aMab athletes to a Gamma distribution ... 0: (-564.227696) alpha=19.998389, beta=3.001330, off=9.999839 384: (-256.999252) alpha=15.951361, beta=2.194827, off=13.795466 ML : (-248.056301) alpha=8.610936, beta=1.673886, off=15.343436 median: (-249.115483) alpha=12.411010, beta=2.005287, off=14.410945

(Click here to show the code)

```
ndim = 3
def lnLike_weibulloffset(theta, x):
k, beta, off = theta
x_adj = x - off
if np.any( x_adj < 0 ):
return -np.inf
return np.sum( np.log(k*beta) + (k-1)*np.log(x_adj*beta) - (x*beta)**k )
lnPrior_weibulloffset = lnPrior_loggaussian
def lnProb_weibulloffset( theta, x ):
temp = lnPrior_weibulloffset( theta )
if temp == -np.inf:
return temp
return temp + lnLike_weibulloffset(theta, x)
pos = [np.array([8,.1,1.]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
print("Fitting the adjusted BMI of hT aMab athletes to a Weibull distribution ...")
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnProb_weibulloffset, threads=os.cpu_count(), args=[ x ] )
weibull_mean = np.mean( pos, axis=0 )
print("{:5}: ({:5f}) k={:7f}, beta={:7f}, off={:7f}".format(0,
lnProb_weibulloffset( weibull_mean, x ), *weibull_mean))
for it in range(nsamples):
global pos, lnprob
pos, lnprob, _ = sampler.run_mcmc( pos, 64, storechain=False )
model_mean = np.mean( pos, axis=0 )
print("{:>5}: ({:5f}) k={:7f}, beta={:7f}, off={:7f}".format((it+1) * 64, \
lnProb_weibulloffset( model_mean, x ), *model_mean))
best = np.argmax( lnprob )
print("{:>5}: ({:5f}) k={:7f}, beta={:7f}, off={:7f}".format("ML", \
lnprob[best], *pos[best]) )
model_median = np.median( pos, axis=0 )
print("{:>5}: ({:5f}) k={:7f}, beta={:7f}, off={:7f}".format("median", \
lnProb_weibulloffset( model_median, x ), *model_median))
models['hT_aMab_BMI_weibull'] = pos
```

Fitting the adjusted BMI of hT aMab athletes to a Weibull distribution ... 0: (-48865.772268) k=7.999859, beta=0.099877, off=0.999138 384: (-271.350390) k=9.937527, beta=0.046958, off=0.019000 ML: (-270.340284) k=9.914647, beta=0.046903, off=0.000871 median: (-270.974131) k=9.833793, beta=0.046947, off=0.011727

(Click here to show the code)

```
ndim = 2
def lnLike_rayleighoffset(theta, x):
tau, off = theta
x_adj = x - off
if np.any( x_adj < 0 ):
return -np.inf
return np.sum( np.log(x_adj*tau) - .5*x_adj*x_adj*tau )
def lnPrior_rayleighoffset( theta ):
tau, off = theta
if (tau <= 0) or (off < 0):
return -np.inf
return 0
def lnProb_rayleighoffset( theta, x ):
temp = lnPrior_rayleighoffset( theta )
if temp == -np.inf:
return temp
return temp + lnLike_rayleighoffset(theta, x)
pos = [np.array([.5,10]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
print("Fitting the adjusted BMI of hT aMab athletes to a Rayleigh distribution ...")
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnProb_rayleighoffset, threads=os.cpu_count(), args=[ x ] )
rayleigh_mean = np.mean( pos, axis=0 )
print("{:5}: ({:5f}) tau={:7f}, off={:7f}".format(0,
lnProb_rayleighoffset( rayleigh_mean, x ), *rayleigh_mean))
for it in range(nsamples):
global pos, lnprob
pos, lnprob, _ = sampler.run_mcmc( pos, 64, storechain=False )
rayleigh_mean = np.mean( pos, axis=0 )
print("{:5}: ({:5f}) tau={:7f}, off={:7f}".format((it+1)*64,
lnProb_rayleighoffset( rayleigh_mean, x ), *rayleigh_mean))
best = np.argmax( lnprob )
print("{:>5}: ({:5f}) tau={:7f}, off={:7f}".format("ML", lnprob[best], *pos[best]))
rayleigh_median = np.median( pos, axis=0 )
print("{:5}: ({:5f}) tau={:7f}, off={:7f}".format("median",
lnProb_rayleighoffset( rayleigh_median, x ), *rayleigh_median))
models['hT_aMab_BMI_rayleigh'] = pos
```

Fitting the adjusted BMI of hT aMab athletes to a Rayleigh distribution ... 0: (-3378.099000) tau=0.499136, off=9.999193 384: (-254.717778) tau=0.107962, off=16.378780 ML: (-253.012418) tau=0.110751, off=16.574934 median: (-253.092584) tau=0.108740, off=16.532576

(Click here to show the code)

```
minVal = min(BMI_adj[BMI_adj_valid])
maxVal = max(BMI_adj[BMI_adj_valid])
x = np.linspace(minVal,maxVal,255)
bins = 9
pp.hist( BMI_adj[mask_gt_5nmol & mask_amab & BMI_adj_valid], bins, density=1, facecolor='green', alpha=0.2)
pp.plot( x, plotWrapper(lnLike_gaussian, np.median(models['hT_aMab_BMI_gaussian'],axis=0), x ), 'k' )
pp.plot( x, plotWrapper(lnLike_loggaussian, np.median(models['hT_aMab_BMI_loggaussian'],axis=0), x ), 'r', alpha=.2 )
pp.plot( x, plotWrapper(lnLike_gammaoffset, np.median(models['hT_aMab_BMI_gamma'],axis=0), x ), 'g' )
pp.plot( x, plotWrapper(lnLike_weibulloffset, np.median(models['hT_aMab_BMI_weibull'],axis=0), x ), 'b', alpha=.2 )
pp.plot( x, plotWrapper(lnLike_rayleighoffset, np.median(models['hT_aMab_BMI_rayleigh'],axis=0), x ), 'y' )
pp.legend(['Gaussian','log-Gaussian + Offset','Gamma + Offset', 'Weibull + Offset','Rayleigh + Offset','high-T aMab'], loc=0)
pp.title('Adjusted BMI, elite athletes')
pp.xlim([minVal,maxVal])
pp.ylim([0,.3])
pp.yticks([])
pp.show()
```

Looks like the Gamma distribution is the best of the bunch, though only if you use the median or maximal likelihood of the posterior. There must be some outliers in there that are tugging the mean around. Visually, there isn't too much difference between the Gaussian and Gamma fits, but the Rayleigh seems artificially sharp on the low end. It's a bit of a shame, the Gamma distribution is usually related to rates and variance so we don't have a good reason for applying it here, other than "it fits the best." We might be able to do better with a per-sport Gaussian distribution fit, but for now I'm happy with the Gamma.

Time to fit the other pool of athletes, and chart it all.

(Click here to show the code)

```
x = BMI_adj[mask_lt_5nmol & BMI_adj_valid & mask_afab]
ndim = 3
pos = [np.array([20,3.,10.]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
print("Fitting the adjusted BMI of lT aFab athletes to a Gamma distribution ...")
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnProb_gammaoffset, threads=os.cpu_count(), args=[ x ] )
gammaoffset_mean = np.mean( pos, axis=0 )
print("{:5}: ({:5f}) alpha={:7f}, beta={:7f}, off={:7f}".format(0,
lnProb_gammaoffset( gammaoffset_mean, x ), *gammaoffset_mean))
for it in range(nsamples):
global pos, lnprob
pos, lnprob, _ = sampler.run_mcmc( pos, 64, storechain=False )
model_mean = np.mean( pos, axis=0 )
print("{:6}: ({:5f}) alpha={:7f}, beta={:7f}, off={:7f}".format((it+1) * 64, \
lnProb_gammaoffset( model_mean, x ), *model_mean))
best = np.argmax( lnprob )
print("{:6}: ({:5f}) alpha={:7f}, beta={:7f}, off={:7f}".format("ML", \
lnprob[best], *pos[best]) )
model_median = np.median( pos, axis=0 )
print("{:6}: ({:5f}) alpha={:7f}, beta={:7f}, off={:7f}".format("median", \
lnProb_gammaoffset( model_median, x ), *model_median))
models['lT_aFab_BMI_gamma'] = pos
```

Fitting the adjusted BMI of lT aFab athletes to a Gamma distribution ... 0: (-127.467934) alpha=20.000007, beta=3.000116, off=9.999921 384: (-128.564564) alpha=15.481265, beta=3.161022, off=12.654149 ML : (-117.582454) alpha=2.927721, beta=1.294851, off=14.713479 median: (-120.689425) alpha=11.961847, beta=2.836153, off=13.008723

(Click here to show the code)

```
x = np.linspace(minVal,maxVal,255)
bins = 9
pp.hist( BMI_adj[mask_lt_5nmol & mask_afab & BMI_adj_valid], bins, density=1, facecolor='black', alpha=0.2)
pp.hist( BMI_adj[mask_gt_5nmol & mask_amab & BMI_adj_valid], bins, density=1, facecolor='green', alpha=0.2)
pp.legend(['lT aFab','hT aMab'], loc=0)
pp.plot( x, plotWrapper(lnLike_gammaoffset, np.median(models['lT_aFab_BMI_gamma'],axis=0), x ), 'k', alpha=.5 )
pp.plot( x, plotWrapper(lnLike_gammaoffset, np.median(models['hT_aMab_BMI_gamma'],axis=0), x ), 'g', alpha=.5 )
pp.legend(['Gamma + Offset, lT aFab', 'Gamma + Offset, hT aMab','low-T aFab','high-T aMab'], loc=0)
pp.title('Adjusted BMI, elite athletes')
pp.xlim([minVal,maxVal])
pp.yticks([])
pp.show()
```

Those models look pretty reasonable, though the upper end of the assigned-female distribution could be improved on. It's a good enough fit to get some answers, at least.

It's easier to combine step 3, applying the model, with step 5, calculating the Bayes Factor, when writing the code. The resulting Bayes Factor has a probability distribution, as the uncertainty contained in the posterior contaminates it.

(Click here to show the code)

```
logBF_gt5_afab_height = list()
gt5_afab_height = dataset['height'][mask_gt_5nmol & height_valid & mask_afab]
numer_gt5_afab_height = [lnProb_gaussian( pos, gt5_afab_height ) for pos in models['lT_aFab_height']]
denom_gt5_afab_height = [lnProb_gaussian( pos, gt5_afab_height ) for pos in models['hT_aMab_height']]
for laf in numer_gt5_afab_height:
for ham in denom_gt5_afab_height:
logBF_gt5_afab_height.append( laf-ham )
print("Summary of the BF distribution, for the height of >5nmol/L aFab athletes")
percentiles = np.percentile( logBF_gt5_afab_height, [5,16,50,84,95] )
geo_mean = np.exp(np.mean( logBF_gt5_afab_height ))
temp = np.exp(logBF_gt5_afab_height)
mean = np.mean( temp )
favour = np.sum( temp > 1 ) / len(temp)
decisive = np.sum( temp > 19 ) / len(temp)
print("{:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10}".format(
"n","mean","geo.mean","5%","16%","50%","84%","95%"))
print("{:>10} {:10.2f} {:10.2f} {:10.2f} {:10.2f} {:10.2f} {:10.2f} {:10.2f}".format(
len(gt5_afab_height), mean, geo_mean, *np.exp(percentiles) ))
print()
print("Percentage of BF's that favoured the primary hypothesis: {:.2f}%".format(favour*100.))
print("Percentage of BF's that were 'decisive': {:.2f}%".format(decisive*100.))
bins = 21
pp.hist( logBF_gt5_afab_height, bins, density=1, facecolor='black', alpha=0.2)
pp.legend(['high-T aFab'], loc=0)
pp.axvline( x=percentiles[2], color='r' )
pp.axvline( x=percentiles[1], color='r', alpha=.5 )
pp.axvline( x=percentiles[3], color='r', alpha=.5 )
pp.axvline( x=percentiles[0], color='r', alpha=.1 )
pp.axvline( x=percentiles[4], color='r', alpha=.1 )
pp.xticks( np.linspace(-2,6,5), ["{:.1f}".format(x) for x in np.exp(np.linspace(-2,6,5))] )
pp.yticks([])
pp.title('Bayes factor, height, >5nmol/L aFab athletes')
pp.show()
```

Summary of the BF distribution, for the height of >5nmol/L aFab athletes n mean geo.mean 5% 16% 50% 84% 95% 19 10.64 5.44 0.75 1.76 5.66 17.33 35.42 Percentage of BF's that favoured the primary hypothesis: 92.42% Percentage of BF's that were 'decisive': 14.17%

That looks a lot like a log-Gaussian distribution. The arthithmetic mean fails us here, thanks to the huge range of values, so the geometric mean and median are better measures of central tendency.

The best way I can interpret this result is via an eight-sided die: our credence in the hypothesis that >5nmol/L aFab athletes are more like their >5nmol/L aMab peers than their <5nmol/L aFab ones is similar to the credence we'd place on rolling a one via that die, while our credence on the primary hypothesis is similar to rolling any other number *except* one. About 92% of the calculated Bayes Factors were favourable to the primary hypothesis, and about 16% of them crossed the 19:1 threshold, a close match for the asserted evidential bar in science.

That's strong evidence for a mere 19 athletes, though not quite conclusive. How about the Bayes Factor for the height of <5nmol/L aMab athletes?

(Click here to show the code)

```
logBF_lt5_amab_height = list()
lt5_amab_height = dataset['height'][mask_lt_5nmol & height_valid & mask_amab]
numer_lt5_amab_height = [lnProb_gaussian( pos, lt5_amab_height ) for pos in models['hT_aMab_height']]
denom_lt5_amab_height = [lnProb_gaussian( pos, lt5_amab_height ) for pos in models['lT_aFab_height']]
for laf in numer_lt5_amab_height:
for ham in denom_lt5_amab_height:
logBF_lt5_amab_height.append( laf-ham )
print("Summary of the BF distribution, for the height of <5nmol/L aMab athletes") percentiles = np.percentile( logBF_lt5_amab_height, [5,16,50,84,95] ) geo_mean = np.exp(np.mean( logBF_lt5_amab_height )) temp = np.exp(logBF_lt5_amab_height) mean = np.mean( temp ) favour = np.sum( temp > 1 ) / len(temp)
decisive = np.sum( temp > 19 ) / len(temp)
print("{:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10}".format(
"n","mean","geo.mean","5%","16%","50%","84%","95%"))
print("{:>10} {:10.2e} {:10.2e} {:10.2e} {:10.2e} {:10.2e} {:10.2e} {:10.2e}".format(
len(lt5_amab_height), mean, geo_mean, *np.exp(percentiles) ))
print()
print("Percentage of BF's that favoured the primary hypothesis: {:.2f}%".format(favour*100.))
print("Percentage of BF's that were 'decisive': {:.2f}%".format(decisive*100.))
bins = 21
pp.hist( logBF_lt5_amab_height, bins, density=1, facecolor='black', alpha=0.2)
pp.legend(['low-T aMab'], loc=0)
pp.axvline( x=percentiles[2], color='r' )
pp.axvline( x=percentiles[1], color='r', alpha=.5 )
pp.axvline( x=percentiles[3], color='r', alpha=.5 )
pp.axvline( x=percentiles[0], color='r', alpha=.1 )
pp.axvline( x=percentiles[4], color='r', alpha=.1 )
pp.xticks( np.linspace(30,55,7), ["{:.1e}".format(x) for x in np.exp(np.linspace(30,55,7))], rotation='vertical' )
pp.yticks([])
pp.title('Bayes factor, height, <5nmol/L aMab athletes')
pp.show()
```

Summary of the BF distribution, for the height of <5nmol/L aMab athletes n mean geo.mean 5% 16% 50% 84% 95% 26 4.67e+21 3.49e+18 5.67e+14 2.41e+16 5.35e+18 4.16e+20 4.61e+21 Percentage of BF's that favoured the primary hypothesis: 100.00% Percentage of BF's that were 'decisive': 100.00%

Wow! Even with 26 data points, our primary hypothesis was *extremely* well supported. Betting against that hypothesis is like betting a particular person in the US will be hit by lightning three times in a single year!

That seems a little *too* favourable to my view, though. Did something go wrong with the mathematics? The simplest check is to graph the models against the data they're evaluating.

(Click here to show the code)

```
minVal = min(dataset['height'][height_valid])
maxVal = max(dataset['height'][height_valid])
x = np.linspace(minVal,maxVal,255)
bins = 8
pp.hist( lt5_amab_height, bins, density=1, facecolor='green', alpha=0.2)
pp.plot( x, plotWrapper(lnLike_gaussian, np.mean(models['lT_aFab_height'],axis=0), x ), 'k' )
pp.plot( x, plotWrapper(lnLike_gaussian, np.mean(models['hT_aMab_height'],axis=0), x ), 'g' )
pp.legend(['low-T aFab', 'high-T aMab', 'low-T aMab'], loc=0)
pp.yticks([])
# pp.ylim([0,.45])
pp.title('Height, elite athletes, <5nmol/L aMab athletes')
pp.show()
```

Nope, the underlying data genuinely is a better fit for the high-testosterone aMab model. But *that* good of a fit? In linear space, we multiply each of the individual probabilities to arrive at the Bayes factor. That's equivalent to raising the geometric mean to the *n*th power, where *n* is the number of athletes. Since *n* = 26 here, even a geometric mean barely above one can generate a big Bayes factor.

(Click here to show the code)

```
temp = [lnProb_gaussian( np.median(models['hT_aMab_height'],axis=0), x ) - \
lnProb_gaussian( np.median(models['lT_aFab_height'],axis=0), x ) for x in lt5_amab_height]
print( "{}th root of the median Bayes factor of the high-T aMab model applied to low-T aMab athletes: {:.4f}".format(len(lt5_amab_height), \
np.exp(percentiles[2]/len(temp))) )
print( "{}th root of the Bayes factor for the median marginal: {:.4f}".format(len(lt5_amab_height), \
np.exp(np.mean(temp))) )
```

26th root of the median Bayes factor of the high-T aMab model applied to low-T aMab athletes: 5.2519 26th root of the Bayes factor for the median marginal: 3.6010

Note that the Bayes factor we generate by using the median of the marginal for each parameter isn't as strong as the median Bayes factor in the above convolution. That's simply because I'm using a small sample from the posterior distribution. Keeping more samples would have brought those two values closer together, but also greatly increased the amount of computation I needed to do to generate all those Bayes factors.

With that check out of the way, we can move on to \(\overline{BMI}\).

(Click here to show the code)

```
logBF_gt5_afab_BMI = list()
gt5_afab_BMI_invalid = [0,0]
gt5_afab_BMI = BMI_adj[mask_gt_5nmol & BMI_adj_valid & mask_afab]
numer_gt5_afab_BMI = [lnProb_gammaoffset( pos, gt5_afab_BMI ) for pos in models['lT_aFab_BMI_gamma']]
denom_gt5_afab_BMI = [lnProb_gammaoffset( pos, gt5_afab_BMI ) for pos in models['hT_aMab_BMI_gamma']]
for laf in numer_gt5_afab_BMI:
for ham in denom_gt5_afab_BMI:
if not np.isfinite(laf):
gt5_afab_BMI_invalid[0] += 1
continue
elif not np.isfinite(ham):
gt5_afab_BMI_invalid[1] += 1
continue
else:
logBF_gt5_afab_BMI.append( laf-ham )
print("Summary of the BF distribution, for the adjusted BMI of >5nmol/L aFab athletes")
percentiles = np.percentile( logBF_gt5_afab_BMI, [5,16,50,84,95] )
geo_mean = np.exp( np.mean(logBF_gt5_afab_BMI) )
temp = np.exp(logBF_gt5_afab_BMI)
mean = np.mean( temp )
favour = np.sum( temp > 1 ) / len(temp)
decisive = np.sum( temp > 19 ) / len(temp)
print("{:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10}".format(
"n","mean","geo.mean","5%","16%","50%","84%","95%"))
print("{:>10} {:10.2e} {:10.2e} {:10.2e} {:10.2e} {:10.2e} {:10.2e} {:10.2e}".format(
len(gt5_afab_BMI), mean, geo_mean, *np.exp(percentiles) ))
print()
print("Percentage of BF's that favoured the primary hypothesis: {:.2f}%".format(favour*100.))
print("Percentage of BF's that were 'decisive': {:.2f}%".format(decisive*100.))
print("Percentage of non-finite probabilities, when applying the low-T aFab model to high-T aFab athletes: {:.2f}%".format(\
gt5_afab_BMI_invalid[0]*100./(sum(gt5_afab_BMI_invalid)+len(logBF_gt5_afab_BMI)) ))
print("Percentage of non-finite probabilities, when applying the high-T aMab model to high-T aFab athletes: {:.2f}%".format(\
gt5_afab_BMI_invalid[1]*100./(sum(gt5_afab_BMI_invalid)+len(logBF_gt5_afab_BMI)) ))
bins = 20
pp.hist( logBF_gt5_afab_BMI, bins, density=1, facecolor='black', alpha=0.2)
pp.legend(['high-T aFab'], loc=0)
pp.axvline( x=percentiles[2], color='r' )
pp.axvline( x=percentiles[1], color='r', alpha=.5 )
pp.axvline( x=percentiles[3], color='r', alpha=.5 )
pp.axvline( x=percentiles[0], color='r', alpha=.1 )
pp.axvline( x=percentiles[4], color='r', alpha=.1 )
pp.xticks( np.linspace(0,50,7), ["{:.1e}".format(x) for x in np.exp(np.linspace(0,50,7))], rotation='vertical' )
pp.yticks([])
pp.title('Bayes factor, BMI, >5nmol/L aFab athletes')
pp.show()
```

Summary of the BF distribution, for the adjusted BMI of >5nmol/L aFab athletes n mean geo.mean 5% 16% 50% 84% 95% 4 1.70e+12 1.06e+05 2.31e+02 1.60e+03 4.40e+04 3.66e+06 3.99e+09 Percentage of BF's that favoured the primary hypothesis: 100.00% Percentage of BF's that were 'decisive': 99.53% Percentage of non-finite probabilities, when applying the low-T aFab model to high-T aFab athletes: 0.00% Percentage of non-finite probabilities, when applying the high-T aMab model to high-T aFab athletes: 10.94%

This distribution is much stranger, with a number of extremely high BF's that badly skew the mean. The offset contributes to this, with 7-12% of the model posteriors for high-T aMab athletes assigning a zero percent likelihood to an adjusted BMI. Those are excluded from the analysis, but they suggest the high-T aMab model poorly describes high-T aFab athletes.

Our credence in the primary hypothesis here is similar to our credence that an elite golfer *will not* land a hole-in-one on their next shot. That's surprisingly strong, given we're only dealing with *four* datapoints. More data may water that down, but it's unlikely to overcome that extreme level of credence.

(Click here to show the code)

```
logBF_lt5_amab_BMI = list()
lt5_amab_BMI_invalid = [0,0]
lt5_amab_BMI = BMI_adj[mask_lt_5nmol & BMI_adj_valid & mask_amab]
numer_lt5_amab_BMI = [lnProb_gammaoffset( pos, lt5_amab_BMI ) for pos in models['hT_aMab_BMI_gamma']]
denom_lt5_amab_BMI = [lnProb_gammaoffset( pos, lt5_amab_BMI ) for pos in models['lT_aFab_BMI_gamma']]
for laf in numer_lt5_amab_BMI:
for ham in denom_lt5_amab_BMI:
if not np.isfinite(laf):
lt5_amab_BMI_invalid[0] += 1
continue
elif not np.isfinite(ham):
lt5_amab_BMI_invalid[1] += 1
continue
else:
logBF_lt5_amab_BMI.append( laf-ham )
print("Summary of the BF distribution, for the adjusted BMI of <5nmol/L aMab athletes") percentiles = np.percentile( logBF_lt5_amab_BMI, [5,16,50,84,95] ) geo_mean = np.exp( np.mean(logBF_lt5_amab_BMI) ) temp = np.exp(logBF_lt5_amab_BMI) mean = np.mean( temp ) favour = np.sum( temp > 1 ) / len(temp)
decisive = np.sum( temp > 19 ) / len(temp)
print("{:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10}".format(
"n","mean","geo.mean","5%","16%","50%","84%","95%"))
print("{:>10} {:10.2e} {:10.2e} {:10.2e} {:10.2e} {:10.2e} {:10.2e} {:10.2e}".format(
len(lt5_amab_BMI), mean, geo_mean, *np.exp(percentiles) ))
print()
print("Percentage of BF's that favoured the primary hypothesis: {:.2f}%".format(favour*100.))
print("Percentage of BF's that were 'decisive': {:.2f}%".format(decisive*100.))
print("Percentage of non-finite probabilities, when applying the high-T aMab model to low-T aMab athletes: {:.2f}%".format(\
lt5_amab_BMI_invalid[0]*100./(sum(lt5_amab_BMI_invalid)+len(logBF_lt5_amab_BMI)) ))
print("Percentage of non-finite probabilities, when applying the low-T aFab model to low-T aMab athletes: {:.2f}%".format(\
lt5_amab_BMI_invalid[1]*100./(sum(lt5_amab_BMI_invalid)+len(logBF_lt5_amab_BMI)) ))
bins = 20
pp.hist( logBF_lt5_amab_BMI, bins, density=1, facecolor='black', alpha=0.2)
pp.legend(['low-T aMab'], loc=0)
pp.axvline( x=percentiles[2], color='r' )
pp.axvline( x=percentiles[1], color='r', alpha=.5 )
pp.axvline( x=percentiles[3], color='r', alpha=.5 )
pp.axvline( x=percentiles[0], color='r', alpha=.1 )
pp.axvline( x=percentiles[4], color='r', alpha=.1 )
pp.xticks( np.linspace(20,100,10), ["{:.1e}".format(x) for x in np.exp(np.linspace(20,100,10))], rotation='vertical' )
pp.yticks([])
pp.title('Bayes factor, BMI, <5nmol/L aMab athletes')
pp.show()
```

Summary of the BF distribution, for the adjusted BMI of <5nmol/L aMab athletes n mean geo.mean 5% 16% 50% 84% 95% 9 6.64e+35 2.07e+22 4.05e+12 4.55e+16 6.31e+21 7.72e+27 9.81e+32 Percentage of BF's that favoured the primary hypothesis: 100.00% Percentage of BF's that were 'decisive': 100.00% Percentage of non-finite probabilities, when applying the high-T aMab model to low-T aMab athletes: 0.00% Percentage of non-finite probabilities, when applying the low-T aFab model to low-T aMab athletes: 0.00%

The hypotheses' Bayes factor for the adjusted BMI of low-testosterone aMab athletes is much better behaved. Even here, the credence is above three-lightning-strikes territory, pretty decisively favouring the hypothesis.

Our final step would normally be to combine all these individual Bayes factors into a single one. That involves multiplying them all together, however, and a small number multiplied by a very large one is an even larger one. It isn't worth the effort, the conclusion is pretty obvious.

Our primary hypothesis is on quite solid ground: Athletes with exceptional testosterone levels are more like athletes of the same sex but with typical testosterone levels, than they are of other athletes with a different sex but similar testosterone levels. If we divide up sports by testosterone level, then, roughly 6-8% of assigned-male athletes will wind up in the <5 nmol/L group, and about the same share of assigned-female athletes will be in the >5 nmol/L group. Note, however, that it doesn't follow that 6-8% of those in the <5 nmol/L group will be assigned-male. About 41% of the athletes at the 2018 Olymics were assigned-female, for instance. If we fix the rate of exceptional testosterone levels at 7%, and assume PyeongChang's rate is typical, a quick application of Bayes' theorem reveals

$$ \begin{align} p( \text{aMab} \mid \text{<5nmol/L} ) &= \frac{ p( \text{<5nmol/L} \mid \text{aMab} ) p( \text{aMab} ) }{ p( \text{<5nmol/L} \mid \text{aMab} ) p( \text{aMab} ) + p( \text{<5nmol/L} \mid \text{aFab} ) p( \text{aFab} ) } \\ {} &= \frac{ 0.07 \cdot 0.59 }{ 0.07 \cdot 0.59 + 0.93 \cdot 0.41 } \\ {} &\approx 9.8\% \end{align} $$

If all those assumptions are accurate, about 10% of <5 nmol/L athletes will be assigned-male, more-or-less matching the number I calculated way back at the start. In sports where performance is heavily correlated with height or \(\overline{BMI}\), then, the 10% of assigned-male athletes in the <5 nmol group will heavily dominate the rankings. The odds of a woman earning recognition in this sport are negligible, leading many of them to drop out. This increases the proportion of men in that sport, leading to more domination of the rankings, more women dropping out, and a nasty feedback loop.

Conversely, about 5% of >5nmol/L athletes will be assigned-female. In a heavily-correlated sport, those women will be outclassed by the men and have little chance of earning recognition for their achievements. They have no incentive to compete, so they'll likely drop out or avoid these sports as well.

In events where physicality has less or no correlation with sporting performance, these effects will be less pronounced or non-existent, of course. But this still translates into fewer assigned-female athletes competing than in the current system.

But it gets worse! We'd also expect an uptick in the number of assigned-female athletes doping, primarily with testosterone inhibitors to bring themselves just below the 5nmol/L line. Alternatively, high-testosterone aFab athletes may inject large doses of testosterone to bulk up and remain competitive with their assigned-male competitors.

By dividing up testosterone levels into only two categories, sporting authorities are implicitly stating that everyone within those categories is identical. A number of athletes would likely go to court to argue that boosting or inhibiting testosterone should be legal, provided they do not cross the 5nmol/L line. If they're successful, then either the rules around testosterone usage would be relaxed, or sporting authorities would be forced to subdivide these groups further. This would lead to an uptick in testosterone doping among *all athletes*, not just those assigned female.

Notice that assigned-male athletes don't have the same incentives to drop out, and in fact the low-testosterone subgroup may even be encouraged to compete as they have an easier path to sporting fame and glory. Sports where performance is heavily correlated with height or \(\overline{BMI}\) will come to be dominated by men.

[1:15] In a nutshell, I find the arguments and logic that currently permit transgender women to compete against biological women to be remarkably flawed, and I’m convinced that unless quickly rectified, this will KILL women’s sports.

[14:00] I don’t want to see the day when women’s athletics is dominated by Y chromosomes, but without a change in policy, that is precisely what’s going to happen.

It's rather astounding. Transgender athletes are a not a problem, on several levels; as I've pointed out before, they've been allowed to compete in the category they identify for over a decade in some places, and yet no transgender athlete has come to dominate *any* sport. The Olympics has held the door open since 2004, and not a single transgender athlete has ever openly competed as a transgender athlete. Rationality Rules, like other transphobes, is forced to cherry-pick and commit lies of omission among a handful of examples, inflating them to seem more significant than they actually are.

In response to this non-existent problem, Rationality Rules' proposed solution would accomplish the very thing he wants to avoid! You don't get that turned around if you're a rational person with a firm grasp on the science.

No, this level of self-sabotage is only possible if you're a clueless bigot who's ignorant of the relevant science, and so frightened of transgender people that your critical thinking skills abandon you. The vast difference between what Rationality Rules claims the science says, and what his own citations say, must be because he knows that if he puts on a good enough act nobody will check his work. Everyone will walk away assuming he's rational, rather than a scared, dishonest loon.

It's hard to fit any other conclusion to the data.

Oh HEY, didn’t see you there! Sorry for the absence, but life caught up to me again. I was doggedly plugging away on my next post when I stumbled across a comic that must have been tailor-made for my current headspace, and it simply *had* to jump the queue.

With the benefit of hindsight, I can see another omission from Rationality Rules’ latest transphobic video. In his citations, he cites two sporting bodies: the International Association of Athletics Federations and the Australian Sports Anti-Doping Authority. He relies heavily on the former, which is strange. The World Medical Association has condemned the IAAF’s policies on intersex and transgender athletes as “contrary to international medical ethics and human rights standards.” The IAAF has defended itself, in part, by arguing this:

The IAAF is not a public authority, exercising state powers, but rather a private body exercising private (contractual) powers. Therefore, it is not subject to human rights instruments such as the Universal Declaration of Human Rights or the European Convention on Human Rights.

Which is A) not a good look, and B) false. If you won’t take my word on that last one, maybe you’ll take the UN’s? [Read more…]

Essence of Thought has published a timeline of the Rationality Rules affair. If you’re missed any of the last five months, it’ll bring you up to speed.

Cripes, has it been that long already?! I had a look through my archives, and all but two of my posts over the last two months have been focused on Rationality Rules, and even those two were about transphobia. I know, I know, the constant drumbeat is getting a bit repetitive and boring. But there’s a reason for it.

[11:31] Now, some of the walkouts had formed a support group, which I was later added to, and reading through their accounts is truly horrifying. Many discussed the abuse they suffered thanks to Woodford and his audience. There are numerous discussions on how their sleep was impacted, about how they’re having to see psychiatrists and other specialists. I’ve even seen [a post?] discussing suicide in relation to what had occurred. That’s the level of severity we are talking about with this issue: people discussing suicide. That’s the damage Woodford and his supporters have caused this one group, this one organization.

I don’t have any way to verify this part, but some of it tracks with comments I’ve read elsewhere, the claims have remained consistent over time, and it would explain why ACA members seem willing to talk to Essence of Thought despite the ocean between them.

One thing I do know: the odds of anyone holding Rationality Rules responsible are basically zero. Some big names in the atheo-skeptic sphere, such as Matt Dillahunty and AronRa, either agree with RR or don’t care enough to do their homework. The ACA tried to do the right thing, but it appears RR supporters elected themselves into a majority on the ACA’s board, possibly breaking the rules in the process, and promptly started kissing their abuser’s ass.

In order to remove any ambiguity in the following statement, I wish to make clear that the ACA earnestly and sincerely apologizes to Stephen Woodford (Rationality Rules) for vilifying his character and insinuating that he is opposed to the LGBTQIA+ community. The Board of Directors has officially retracted our original statement.

Rationality Rules was so confident nobody would take him to task, his “improved” video contains the same arguments as his “flawed” one. And honestly, he was right; I’ve seen this scenario play out often enough within this community to know that we try to bury our skeletons, that we treat our minorities like shit, that we “skeptics” are just as prone to being blind followers as the religious/woo crowds we critique. And just like all those other times, I cope by writing words until I get sick of the topic. Sometimes, that takes a while.

This is one of those “a while” times. If it helps, I’m actively trying to avoid covering topics other people already have, and elevating the voices of others to break up the monotony.

- «Previous Page
- 1
- 2
- 3
- 4
- 5
- …
- 11
- Next Page»