showing significance of gender differences in control vs treatment groups

Question

I have 4 groups with independent samples : Female Control Female Treatment Male Control Male Treatment

When I compare Control to Treatment within the same sex, I find that Female Control to Female Treatment is significantly different whereas Male Control to Male Treatment is indifferent.

I'm confused about two things:

1. Do I need to also perform multiple comparison tests (e.g., one way anova with post hoc correction)?

2. If outcome of the multiple comparison test is negative, does it negate my initial findings?

3. Do I really need to use post hoc corrections? I don't really get the rationale for the post hoc corrections. For example, Bonferroni correction divides the p-value threshold by the number of comparisons, which is 6 in my case. However, from the standpoint of my hypotheses, some pairwise comparisons doesn't make any sense to me .

4. If I claim that treatment acts differently in each sex, how can I prove that?

Answer 1

To summarize: You collected data from a (randomized?) experiment to compare treatment vs control. And you have information about a single covariate, gender.

What you have done to analyze the data:

You split the participants by gender and compared treatment and control separately for each group.

What you could do instead:

Model all your data points with a single model that includes main terms for treatment and gender plus a treatment-gender interaction.

Regression is a good place to start. The details of what kind of regression is most appropriate depend on the outcome about which you provide no information.

A regression model can handle gracefully the complexities you are worried about:

• It will estimate all parameters of interest at the same time, so you don't have to make corrections for multiple comparisons afterward.
• It will quantify the evidence that the treatment effect is different in males and females.

Once you fit the model (and verify that it fits the data well), you can perform various comparisons of interests. These are also called contrasts.

The interpretation of the model coefficients depends on how the categorical variables are encoded in the design matrix. This is a big topic of its own; you can learn more about it here.

Instead I show how to estimate contrasts; they are the same no matter the encoding. If that weren't true, the conclusions from an experiment would depend on an arbitrary choice such as the reference level for a categorical variable.

I use R but this analysis can be performed with any statistical software.

library("contrast")
library("tidyverse")

set.seed(1234)

# Simulate data
n <- 100
y <- rnorm(n)
# The data is not balanced in either categorical predictor
gender <- ifelse(runif(n) < .4, "m", "f")
treatment <- ifelse(runif(n) < .33, "t", "c")

# Fit model
fit <- lm(y ~ gender * treatment)


# I jump into making comparisons right away because the coefficients depend on the encoding.


# Compare the two genders in the control group.
contrast(
  fit,
  list(gender = "f", treatment = "c"),
  list(gender = "m", treatment = "c")
)
#> lm model parameter contrast
#> 
#>    Contrast      S.E.      Lower     Upper     t df Pr(>|t|)
#>  -0.3970983 0.2506247 -0.8945844 0.1003879 -1.58 96   0.1164

# The contrast estimate is a difference of sample averages, as we expect.
# But now we also know its standard error and p-value.
mean(y[gender == "f" & treatment == "c"]) - mean(y[gender == "m" & treatment == "c"])
#> [1] -0.3970983


# And here is how to compare the treatment and control (across both genders)
contrast(
  fit,
  list(gender = c("f", "m"), treatment = "t"),
  list(gender = c("f", "m"), treatment = "c"),
  type = "average"
)
#> lm model parameter contrast
#> 
#>     Contrast      S.E.      Lower     Upper     t df Pr(>|t|)
#> 1 -0.1556399 0.2252886 -0.6028342 0.2915544 -0.69 96   0.4913

# This contrast estimate is *not* a difference of sample averages
mean(y[treatment == "t"]) - mean(y[treatment == "c"])
#> [1] -0.1136949
# but a difference in average effects
(
  (mean(y[gender == "f" & treatment == "t"]) + mean(y[gender == "m" & treatment == "t"])) / 2
    -
  (mean(y[gender == "f" & treatment == "c"]) + mean(y[gender == "m" & treatment == "c"])) / 2
)
#> [1] -0.1556399

^{Created on 2022-03-26 by the reprex package (v2.0.1)}