Visualizing Categorical Response Data I have a dataset that contains the "weight", the "gender" and if someone has "asthma" or not.
I am interested in learning about ways to visualize the proportion of people who have asthma in different weight ranges (e.g. what percent of "heavy" people have asthma vs what percent of "non-heavy" people what asthma).
Below, I used some R Code to generate a sample dataset (with some patterns incorporated):
library(dplyr)
library(purrr)
library(ggplot2)
set.seed(123)

my_data1 = data.frame(Weight =  rnorm(500,100,100), asthma = sample(c(0,1), prob = c(0.7,0.3), replace=TRUE, size= 500))
my_data2 = data.frame(Weight = rnorm(500, 200, 50),  asthma = sample(c(0,1), prob = c(0.3,0.7), replace=TRUE, size= 500))
my_data_a = rbind(my_data1, my_data2)
my_data_a$gender =  "male"


my_data1 = data.frame(Weight =  rnorm(500,100,100), asthma = sample(c(0,1), prob = c(0.7,0.3), replace=TRUE, size= 500))
my_data2 = data.frame(Weight = rnorm(500, 200, 50),  asthma = sample(c(0,1), prob = c(0.3,0.7), replace=TRUE, size= 500))
my_data_b = rbind(my_data1, my_data2)
my_data_b$gender =  "female"


my_data = rbind(my_data_a, my_data_b)

I am interested in learning how about different ways to visualize and analyze this data. For example, I can make a plot for this dataset:
ggplot(my_data, aes(x=Weight, y=asthma, color=gender)) +
  geom_point() + 
    ggtitle("Relationship Between Weight and Asthma") + scale_color_manual(values = c("black", "red"), name = "Gender", labels = c("Male", "Female")) +
    theme(legend.position = "bottom")


Unfortunately, this plot does not seem to be very informative.
I thought of two different Approaches ("Approach 1" and "Approach 2") that might be useful in visualizing this data.
Approach 1: Even though in Statistics we are advised never to "bin" a continuous variable as it can lead to p-hacking, add bias and lose information in the process (i.e. inherently arbitrary) - as an example, I binned this data into increasing groups of 10-percentiles and took the average asthma rate for each percentile. The resulting graph looks something like this:
final = my_data %>%
    arrange(Weight) %>%
    mutate(ntile = ntile(Weight, 10)) %>%
    group_by(ntile, gender) %>%
    summarise(mean = mean(asthma))  

ggplot(final, aes(x = ntile, y = mean, color = gender)) + 
  geom_line() +
  xlab("Weight Bin Number") + 
  ylab("Average Asthma Rate") + 
  ggtitle("Relationship Between Weight Bins and Average Asthma Rate") + scale_color_manual(values = c("black", "red"), name = "Gender", labels = c("Male", "Female")) +
  theme(legend.position = "bottom")


Now, it seems like we can see some sort of pattern, regardless of gender, heavier people tend to have asthma at higher rates than less heavy people. However, different binning schemes will inherently produce different patterns. For example:

Here, we can see that too few bins can really flatten out the trend and that too many bins make it too difficult to make out any patterns at all!
Approach 2: I thought of another approach in which I can compare the distributions of Weight by Asthma in different Genders:
a = ggplot(my_data[my_data$gender == "male",], aes(x = Weight, colour = asthma)) +
    geom_density() + 
    ggtitle("Males : Distribution of Weight by Asthma") 

b = ggplot(my_data[my_data$gender == "female",], aes(x = Weight, colour = asthma)) +
    geom_density() + 
    ggtitle("Females : Distribution of Weight by Asthma") 


Initially, I had thought that "Approach 2" might be less likely to distort and add bias to the visualizations compared to "Approach 1" since I am not required to make a decision on how to bin variables - but then I realized that these density visualizations in "Approach 2" are likely based on some kernel smoothing method (e.g. kernel options https://rdrr.io/r/stats/density.html) which will induce some bias (though I am not sure if as strongly as "Approach 1"). For what it's worth, visually, we can see that people with asthma are more concentrated in regions corresponding to heavier weights compared to people without asthma.
In the end, it seems like fitting a Logistic Regression model to this dataset is likely a better choice as it does not directly require the user to make some decision on binning the data - but in the end, are there other statistical methods that are better suited for visualizing these data?
Just to clarify - I am interested in learning about ways to visualize the proportion of people who have asthma in different weight ranges ... and mitigate the risk of engaging in p-hacking.
Thanks!

Another option to consider is to use a smoothed relationship on some suitable scale. I.e. based on some model that suitably transforms proportions to ensure that they are between 0 and 1 - e.g. via a logit-link function - and ensures that the model is flexible enough to display the trends in the data rather than some enforced functional forms such as a linear relationship. One option could be a GAM for a binomial outcome that uses some splines (one can then of course play around with different splines, different smoothing parameters, different models, etc.) or of course you could do something even more non-parametric (to really stay flexible/reflect the data rather than an overly restrictive model). Of course, some things would be easier to plot than others, e.g. GAMs integrate nicely with ggplot2::geom_smooth in R. One could fit some much more sophisticated model, but might then have to get the models predictions and explicitly provide them as data to plot them.
Here's an example implementation:
library(tidyverse)
set.seed(123)

my_data1 = data.frame(Weight =  rnorm(500,100,100), asthma = sample(c(0,1), prob = c(0.7,0.3), replace=TRUE, size= 500))
my_data2 = data.frame(Weight = rnorm(500, 200, 50),  asthma = sample(c(0,1), prob = c(0.3,0.7), replace=TRUE, size= 500))
my_data_a = rbind(my_data1, my_data2)
my_data_a$gender =  "male"


my_data1 = data.frame(Weight =  rnorm(500,100,100), asthma = sample(c(0,1), prob = c(0.7,0.3), replace=TRUE, size= 500))
my_data2 = data.frame(Weight = rnorm(500, 200, 50),  asthma = sample(c(0,1), prob = c(0.3,0.7), replace=TRUE, size= 500))
my_data_b = rbind(my_data1, my_data2)
my_data_b$gender =  "female"


my_data = rbind(my_data_a, my_data_b) %>% tibble()

my_data %>%
  ggplot(aes(x=Weight, y=asthma+case_when(asthma==1 & gender=="female" ~ -0.05, asthma==0 & gender!="female" ~ 0.05, TRUE ~ 0), col=gender, fill=gender, label=gender)) +
  theme_bw(base_size=18) +
  geom_smooth(method = "gam", 
              formula = y ~ s(x, bs = "bs"),
              method.args=list(family="binomial", link="logit"),
              se=TRUE, level=0.95) +
  # theme(legend.position="none") +
  #directlabels::geom_dl(method="smart.grid", list(directlabels::dl.trans(x=x-0.1))) + # didn't work so well, but would be ideal, if it did
  geom_jitter(width = 0, height=0.02, alpha=0.5) +
  coord_cartesian(ylim=c(0,1)) +
  scale_y_continuous(label=scales::percent_format(1), breaks=seq(0, 1, 0.25)) +
  ggthemes::scale_color_colorblind() +
  ggthemes::scale_fill_colorblind() +
  xlab("Weight (kg)") +
  ylab("Individuals with asthma (%)")

This produces the following plot:

I explicitly wrote down the se=TRUE, level=0.95 options to make it clear what is being shown as shaded areas around the curves (i.e. two-sided 95% confidence intervals). That these are 95% ones is of course a completely arbitrary choice, but it is useful to visualize the uncertainty about the fitted curve, because one might e.g. otherwise overinterpret the ends of the curve for females.
