Unbalanced Four Core Genotype Analyses

This came out of a discussion related to a rat model designed to explore the same question at the Four Core Genotype mouse model: what happens if you don’t have all four arms? Are you hopelessly confounded? The short answer is no under some assumptions. We’ll explore that with some simple simulations.

To start with we’ll create a function to simulate the four core genotype model and then drop one of the arms.

library(tidyverse)
library(ggplot2)
library(broom)

generateFCGData <- function(
    XXM = 0,
    XXF = 0,
    XYM = 0,
    XYF = 0,
    std = 1,
    ngroup=10) {
  rbind(
    data.frame(group="XXM", gonad="T", chromosomes="XX", 
               volume = rnorm(ngroup, XXM, sd=std)),
    data.frame(group="XXF", gonad="O", chromosomes="XX", 
               volume = rnorm(ngroup, XXF, sd=std)),
    data.frame(group="XYM", gonad="T", chromosomes="XY", 
               volume = rnorm(ngroup, XYM, sd=std)),
    data.frame(group="XYF", gonad="O", chromosomes="XY", 
               volume = rnorm(ngroup, XYF, sd=std))
  )
}

Let’s assume that having male gonads increases the output variable by 2 standard deviations. That would look as follows:

ddd <- generateFCGData(XYM=2, XXM=2, ngroup = 500)
  
ggplot(ddd) + 
  aes(x=group, y=volume, fill=chromosomes) + 
  geom_boxplot() + 
  facet_wrap(~gonad, scales="free_x") +
  scale_fill_brewer(palette = "Set1") + 
  theme_classic()

Let’s run the stats three ways, first separately fitting gonads and chromosomes and then using an additive model incorporating both chromosomes and gonads:

summary( ( lfg <- lm(volume ~ gonad, ddd)))

Call:
lm(formula = volume ~ gonad, data = ddd)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3858 -0.6770 -0.0152  0.6524  3.9618 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.03784    0.03136  -1.206    0.228    
gonadT       2.09769    0.04435  47.293   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9918 on 1998 degrees of freedom
Multiple R-squared:  0.5282,    Adjusted R-squared:  0.5279 
F-statistic:  2237 on 1 and 1998 DF,  p-value: < 2.2e-16

summary( ( lfc <- lm(volume ~ chromosomes, ddd)))

Call:
lm(formula = volume ~ chromosomes, data = ddd)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8468 -1.0806  0.0336  1.0483  5.0089 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.009223   0.045660  22.103   <2e-16 ***
chromosomesXY 0.003567   0.064573   0.055    0.956    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.444 on 1998 degrees of freedom
Multiple R-squared:  1.527e-06, Adjusted R-squared:  -0.000499 
F-statistic: 0.003052 on 1 and 1998 DF,  p-value: 0.956

summary( ( lfgc <- lm(volume ~ gonad + chromosomes, ddd)))

Call:
lm(formula = volume ~ gonad + chromosomes, data = ddd)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3840 -0.6771 -0.0163  0.6512  3.9600 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -0.039623   0.038422  -1.031    0.303    
gonadT         2.097692   0.044366  47.282   <2e-16 ***
chromosomesXY  0.003567   0.044366   0.080    0.936    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9921 on 1997 degrees of freedom
Multiple R-squared:  0.5282,    Adjusted R-squared:  0.5277 
F-statistic:  1118 on 2 and 1997 DF,  p-value: < 2.2e-16

The output here is correct - no matter how analyzed, the statistics tell us that having testes increases the output by 2 (everything is in units of standard deviation here), and there is no effect of chromosomes.

Let’s now drop one of the arms, and assume that rather than having all four groups we do not have the XY female group. The exact same data would look like this:

ddd <- ddd %>% filter(group!="XYF")
ggplot(ddd) + 
  aes(x=group, y=volume, fill=chromosomes) + 
  geom_boxplot() + 
  facet_wrap(~gonad, scales="free_x") +
  scale_fill_brewer(palette = "Set1") + 
  theme_classic()

And we run the same statistics:

summary( ( lrg <- lm(volume ~ gonad, ddd)))

Call:
lm(formula = volume ~ gonad, data = ddd)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3858 -0.6819 -0.0288  0.6690  3.9618 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.03693    0.04433  -0.833    0.405    
gonadT       2.09678    0.05430  38.618   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9913 on 1498 degrees of freedom
Multiple R-squared:  0.4989,    Adjusted R-squared:  0.4985 
F-statistic:  1491 on 1 and 1498 DF,  p-value: < 2.2e-16

summary( ( lrc <- lm(volume ~ chromosomes, ddd)))

Call:
lm(formula = volume ~ chromosomes, data = ddd)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8468 -0.9131  0.0039  0.8902  3.9573 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    1.00922    0.04139   24.38   <2e-16 ***
chromosomesXY  1.05511    0.07169   14.72   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.309 on 1498 degrees of freedom
Multiple R-squared:  0.1263,    Adjusted R-squared:  0.1257 
F-statistic: 216.6 on 1 and 1498 DF,  p-value: < 2.2e-16

summary( ( lrgc <- lm(volume ~ gonad + chromosomes, ddd)))

Call:
lm(formula = volume ~ gonad + chromosomes, data = ddd)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3814 -0.6814 -0.0285  0.6713  3.9573 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -0.036928   0.044347  -0.833    0.405    
gonadT         2.092301   0.062716  33.361   <2e-16 ***
chromosomesXY  0.008958   0.062716   0.143    0.886    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9916 on 1497 degrees of freedom
Multiple R-squared:  0.4989,    Adjusted R-squared:  0.4982 
F-statistic: 745.2 on 2 and 1497 DF,  p-value: < 2.2e-16

Losing one arm of the study means that the design is now unbalanced. Looking at the three different stats models shows that modelling only gonads gives the correct answer of a 2 SD increase, but modelling only chromosomes gives an incorrect answer of an increase of 1SD when having XY chromosomes. What saves us is covarying - modelling gonads + chromosomes again returns the correct answer, since the model simultaneously estimates the effects of gonads and chromosomes while controlling for the other.

Let’s show that graphically:

allStats <- rbind(
  tidy(lfg, conf.int = T) %>% mutate(model = "FCG: G"),
  tidy(lfc, conf.int = T) %>% mutate(model = "FCG: C"),
  tidy(lfgc, conf.int = T) %>% mutate(model = "FCG: G + C"),
  tidy(lrg, conf.int = T) %>% mutate(model = "uFCG: G"),
  tidy(lrc, conf.int = T) %>% mutate(model = "uFCG: C"),
  tidy(lrgc, conf.int = T) %>% mutate(model = "uFCG: G + C")
) %>% filter(term != "(Intercept)")

cols <- RColorBrewer::brewer.pal(2, "Set2")

Warning in RColorBrewer::brewer.pal(2, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

allStats %>% 
  ggplot() + 
  aes(x=term, y=estimate, colour=term) + 
  geom_point() + 
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high)) + 
  facet_wrap(.~model, scale="free_x", nrow=1) + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle=45, hjust=1)) + 
  scale_color_manual(values = cols[1:2]) + 
  geom_hline(yintercept = 2, color=cols[2]) +
  geom_hline(yintercept = 0, color=cols[1]) +
  labs(title="Model effects", caption="FCG = balanced four core genotypes; uFCG = unbalanced four core genotypes (no XYF); C = chromosomes; G = gonads") 

Now, can’t we just do a simple t-test comparing two groups with the desired contrasts which avoids the issue of more complex modelling? Yes, that works:

summary(lm(volume ~ group, data = ddd %>% filter(group %in% c("XXF", "XXM"))))

Call:
lm(formula = volume ~ group, data = ddd %>% filter(group %in% 
    c("XXF", "XXM")))

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3814 -0.6694 -0.0223  0.6893  2.9014 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.03693    0.04407  -0.838    0.402    
groupXXM     2.09230    0.06232  33.575   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9853 on 998 degrees of freedom
Multiple R-squared:  0.5304,    Adjusted R-squared:  0.5299 
F-statistic:  1127 on 1 and 998 DF,  p-value: < 2.2e-16

But do we pay a price in statistical power? Let’s compare the simple two group comparison to the linear model with all three groups (in the unbalanced FCG case) estimating both a gonad and a chromosome effect.

simNtimes <- function(nsims=100, effectSize=2, ngroup=20) {
  bFCG <- map(1:nsims, ~ generateFCGData(XYM=effectSize, 
                                         XXM=effectSize, 
                                         ngroup=ngroup))
  uFCG <- map(bFCG, ~ .x %>% filter(group!="XYF"))
  tgroup <- map(bFCG, ~ .x %>% filter(group %in% c("XXF", "XXM")))
  
  lbFCG <- map(bFCG, ~ lm(volume ~ gonad + chromosomes, data=.x))
  luFCG <- map(uFCG, ~ lm(volume ~ gonad + chromosomes, data=.x))
  ltgroup <- map(tgroup, ~ lm(volume ~ group, data=.x))
  
  tlbFCG <- map(lbFCG, tidy)
  tluFCG <- map(luFCG, tidy)
  tltgroup <- map(ltgroup, tidy)
  
  data.frame(
    bFCG = map_dbl(tlbFCG, ~.x %>% filter(term == "gonadT") %>% pull(p.value)),
    uFCG = map_dbl(tluFCG, ~.x %>% filter(term == "gonadT") %>% pull(p.value)),
    twogroup = map_dbl(tltgroup, ~.x %>% filter(term == "groupXXM") %>% pull(p.value))
  )
  
}

allSims <- map(seq(0, 2, by=0.1), ~ simNtimes(nsims=200, effectSize=.x, ngroup = 10) %>% mutate(d=.x))
powerDF <- map_dfr(allSims, ~ .x) %>% group_by(d) %>% summarise_at(vars(bFCG:twogroup), function(x) mean(x<0.05))

powerDF %>% 
  pivot_longer(bFCG:twogroup, values_to = "power") %>% 
  ggplot() + aes(x=d, y=power, colour=name) + 
  geom_smooth(method="gam", formula=y~s(x, k=8), se=F) + 
  scale_color_brewer("Model", palette = "Dark2", 
                     label=c("balanced FCG", 
                             "unbalanced FCG", 
                             "two group comparison")) + 
  theme_classic() + 
  xlab("Simulated effect size") + 
  labs(title="Power of gonad effect", caption ="Assuming n=10 per group")

Not really - this shows the power we lose when moving from a balanced FCG dataset to the unbalanced one missing one group. But the power is the same between the two ways of modelling the unbalanced dataset, indicating that in the case of the unbalanced design (i.e. only having three rather than all four groups of the FCG model) the gonad and chromosome terms of the statistical model essentially reduce down to a two-group comparison.

Let’s expand this to the FCG-like rat model (male arm only for now). There are four genotypes:

Genotype	Gonads	Sex Chromosomes	Sry
XX	Ovaries	XX	None
XX-Sry	Testes	XX	Transgene
XY	Testes	XY	WT
XY-Sry	Testes	XY	WT + Transgene

Once again we need a function to simulate from this dataset

generateFCGRatData <- function(
    XX = 0,
    XXSry = 0,
    XY = 0,
    XYSry = 0,
    std = 1,
    ngroup=10) {
  rbind(
    data.frame(group="XX", gonad="O", chromosomes="XX", sry=0,
               volume = rnorm(ngroup, XX, sd=std)),
    data.frame(group="XXSry", gonad="T", chromosomes="XX", sry=1,
               volume = rnorm(ngroup, XXSry, sd=std)),
    data.frame(group="XY", gonad="T", chromosomes="XY", sry=1,
               volume = rnorm(ngroup, XY, sd=std)),
    data.frame(group="XYSry", gonad="T", chromosomes="XY", sry=2,
               volume = rnorm(ngroup, XYSry, sd=std))
  )
}

Let’s see what a few different scenarios look like:

gonads <- generateFCGRatData(XX=0, XXSry=2, XY=2, XYSry=2, ngroup = 500) %>%
  mutate(contrast = "Gonad")
chrom <- generateFCGRatData(XX=0, XXSry=0, XY=2, XYSry=2, ngroup = 500) %>%
  mutate(contrast = "Chromosomes")
sry <- generateFCGRatData(XX=0, XXSry=2, XY=2, XYSry=4, ngroup = 500) %>%
  mutate(contrast = "Sry")
  
rbind(gonads, chrom, sry) %>%
ggplot() + 
  aes(x=group, y=volume, fill=chromosomes) + 
  geom_boxplot() + 
  facet_wrap(~gonad, scales="free_x") +
  scale_fill_brewer(palette = "Set1") + 
  facet_grid(contrast ~ .) + 
  theme_classic()

Can we appropriately detect the right effects in those models, focusing first on analyzing gonads and chromosomes?

summary(lm(volume ~ gonad + chromosomes, data=gonads))

Call:
lm(formula = volume ~ gonad + chromosomes, data = gonads)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4345 -0.7236  0.0055  0.6643  3.7615 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -0.002916   0.045251  -0.064    0.949    
gonadT         1.967447   0.063994  30.744   <2e-16 ***
chromosomesXY  0.048733   0.055421   0.879    0.379    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.012 on 1997 degrees of freedom
Multiple R-squared:  0.4233,    Adjusted R-squared:  0.4227 
F-statistic: 732.9 on 2 and 1997 DF,  p-value: < 2.2e-16

summary(lm(volume ~ gonad + chromosomes, data=chrom))

Call:
lm(formula = volume ~ gonad + chromosomes, data = chrom)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0717 -0.6588 -0.0195  0.6691  3.3288 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.01701    0.04404   0.386    0.699    
gonadT        -0.05503    0.06228  -0.884    0.377    
chromosomesXY  1.96188    0.05394  36.372   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9848 on 1997 degrees of freedom
Multiple R-squared:  0.4915,    Adjusted R-squared:  0.4909 
F-statistic: 964.9 on 2 and 1997 DF,  p-value: < 2.2e-16

Yes. What about treating sry as an ordered factor (i.e. assuming that a dose of 2 is more than a dose of 1, but not assuming that it is purely linear).

summary(lm(volume ~ sry + chromosomes, data=gonads %>% mutate(sry=as.ordered(sry))))

Call:
lm(formula = volume ~ sry + chromosomes, data = gonads %>% mutate(sry = as.ordered(sry)))

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4345 -0.7225  0.0037  0.6645  3.7615 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    1.32328    0.03991   33.15   <2e-16 ***
sry.L          1.42209    0.06400   22.22   <2e-16 ***
sry.Q         -0.78537    0.03695  -21.25   <2e-16 ***
chromosomesXY  0.02689    0.06400    0.42    0.674    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.012 on 1996 degrees of freedom
Multiple R-squared:  0.4234,    Adjusted R-squared:  0.4226 
F-statistic: 488.6 on 3 and 1996 DF,  p-value: < 2.2e-16

summary(lm(volume ~ sry + chromosomes, data=chrom %>% mutate(sry=as.ordered(sry))))

Call:
lm(formula = volume ~ sry + chromosomes, data = chrom %>% mutate(sry = as.ordered(sry)))

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0717 -0.6636 -0.0133  0.6644  3.3082 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -0.033370   0.038846  -0.859    0.390    
sry.L         -0.067958   0.062293  -1.091    0.275    
sry.Q          0.005698   0.035965   0.158    0.874    
chromosomesXY  1.982422   0.062293  31.824   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9849 on 1996 degrees of freedom
Multiple R-squared:  0.4916,    Adjusted R-squared:  0.4908 
F-statistic: 643.3 on 3 and 1996 DF,  p-value: < 2.2e-16

summary(lm(volume ~ sry + chromosomes, data=sry %>% mutate(sry=as.ordered(sry))))

Call:
lm(formula = volume ~ sry + chromosomes, data = sry %>% mutate(sry = as.ordered(sry)))

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0507 -0.6783 -0.0034  0.6878  3.2917 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.926884   0.040140  48.004   <2e-16 ***
sry.L         2.753428   0.064368  42.777   <2e-16 ***
sry.Q         0.008242   0.037163   0.222    0.825    
chromosomesXY 0.077437   0.064368   1.203    0.229    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.018 on 1996 degrees of freedom
Multiple R-squared:  0.6561,    Adjusted R-squared:  0.6556 
F-statistic:  1269 on 3 and 1996 DF,  p-value: < 2.2e-16

Still works but is a bit more unstable. Some interpretation aids to these models - when treating sry as an ordered factor, there are two possible trends - a linear trends, indicating that with every increase in sry there is an increase in the outcome, and a quadratic trend, indicating that going from 1 to 2 is not the same as going from 0 to 1. So in our third model where we artificially simulated an increase with every dose of sry there is a linear increase (sry.L), but no quadratic increase (sry.Q is non-significant). But in the gonad model there is both a linear and a quadratic effect, since having 2 sry copies (WT + transgene) has the same effect as having one sry copy.

What about power? First just looking at the gonads + chromosomes model:

simRatNtimes <- function(nsims=100, effectSize=2, ngroup=20) {
  ratG <- map(1:nsims, ~ generateFCGRatData(XX=0,
                                            XXSry=effectSize,
                                            XY=effectSize,
                                            XYSry=effectSize,
                                            ngroup=ngroup))
  ratC <- map(1:nsims, ~ generateFCGRatData(XX=0,
                                            XXSry=0,
                                            XY=effectSize,
                                            XYSry=effectSize,
                                            ngroup=ngroup))
  
  lratG <- map(ratG, ~ lm(volume ~ gonad + chromosomes, data=.x))
  lratC <- map(ratC, ~ lm(volume ~ gonad + chromosomes, data=.x))
  
  tlratG <- map(lratG, tidy)
  tlratC <- map(lratC, tidy)
  
  data.frame(
    gonads = map_dbl(tlratG, ~.x %>% filter(term == "gonadT") %>% pull(p.value)),
    chromosomes = map_dbl(tlratC, ~.x %>% filter(term == "chromosomesXY") %>% pull(p.value))
  )
  
}

allRatSims <- map(seq(0, 2, by=0.1), ~ simRatNtimes(nsims=200, effectSize=.x, ngroup = 10) %>% mutate(d=.x))
powerDF <- map_dfr(allRatSims, ~ .x) %>% group_by(d) %>% summarise_at(vars(gonads:chromosomes), function(x) mean(x<0.05))

powerDF %>% 
  pivot_longer(gonads:chromosomes, values_to = "power") %>% 
  ggplot() + aes(x=d, y=power, colour=name) + 
  geom_smooth(method="gam", formula=y~s(x, k=8), se=F) + 
  scale_color_brewer("Model", palette = "Dark2") + 
  theme_classic() + 
  xlab("Simulated effect size") + 
  labs(title="Power of FCG-like rat study", 
       subtitle="Modelling gonads + chromosomes",
       caption ="Assuming n=10 per group")

Given the unbalanced design for gonads there is some power differential, in that it is easier to detect a chromosome effect than a gonad effect.

And if we model sry as an ordered factor?

simRatNtimesOF <- function(nsims=100, effectSize=2, ngroup=20) {
  ratG <- map(1:nsims, ~ generateFCGRatData(XX=0,
                                            XXSry=effectSize,
                                            XY=effectSize,
                                            XYSry=effectSize,
                                            ngroup=ngroup) %>%
    mutate(sry = as.ordered(sry))) 
  ratC <- map(1:nsims, ~ generateFCGRatData(XX=0,
                                            XXSry=0,
                                            XY=effectSize,
                                            XYSry=effectSize,
                                            ngroup=ngroup) %>%
    mutate(sry = as.ordered(sry))) 
  ratS <- map(1:nsims, ~ generateFCGRatData(XX=0,
                                            XXSry=effectSize,
                                            XY=effectSize,
                                            XYSry=effectSize*2,
                                            ngroup=ngroup) %>%
    mutate(sry = as.ordered(sry))) 
  
  lratG <- map(ratG, ~ lm(volume ~ sry + chromosomes, data=.x))
  lratC <- map(ratC, ~ lm(volume ~ sry + chromosomes, data=.x))
  lratS <- map(ratS, ~ lm(volume ~ sry + chromosomes, data=.x))
  
  tlratG <- map(lratG, tidy)
  tlratC <- map(lratC, tidy)
  tlratS <- map(lratS, tidy)
  
  data.frame(
    gonads = map_dbl(tlratG, ~.x %>% filter(term == "sry.L") %>% pull(p.value)),
    chromosomes = map_dbl(tlratC, ~.x %>% filter(term == "chromosomesXY") %>% pull(p.value)),
    sryL = map_dbl(tlratS, ~.x %>% filter(term == "sry.L") %>% pull(p.value)),
    sryQ = map_dbl(tlratS, ~.x %>% filter(term == "sry.Q") %>% pull(p.value))
  )
  
}

allRatSims <- map(seq(0, 2, by=0.1), ~ simRatNtimesOF(nsims=200, effectSize=.x, ngroup = 10) %>% mutate(d=.x))
powerDF <- map_dfr(allRatSims, ~ .x) %>% group_by(d) %>% summarise_at(vars(gonads:sryQ), function(x) mean(x<0.05))

powerDF %>% 
  pivot_longer(gonads:sryQ, values_to = "power") %>% 
  ggplot() + aes(x=d, y=power, colour=name) + 
  geom_smooth(method="gam", formula=y~s(x, k=8), se=F) + 
  scale_color_brewer("Model", palette = "Dark2") + 
  theme_classic() + 
  xlab("Simulated effect size") + 
  labs(title="Power of FCG-like rat study", 
       subtitle="Modelling sry + chromosomes",
       caption ="Assuming n=10 per group")

Here the gonads model looks at the linear part of the sry term when there is a change in the presence of sry but no dosage effect. sryL looks at the linear term when there is a dosage effect, which is now the most powered term (since there are three possible steps), and sryQ looks at the quadratic effect when there is a simulated dosage effect which is rightly 0.