Visualizing Categorical Response Data

Question

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!

Answer 1

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.

Answer 2

Note: Option 1 was supposed to address the OP's original question. Indeed, this is similar to Ben's answer (to the original question). Later the OP updated his question by asking something slightly different. Therefore, Option 2 is my answer to this updated question.

Option 1. If you want to make comparisons in terms of weight, you can plot side-by-side distributions split by a given condition. Here, I am summarising the distributions using boxplots via the ggplot2 package of R; an alternative could be to use violins, as Ben suggests in his answer.

If you want to compare the presence vs absence of asthma, for each type of gender, you could do

my_data$asthma <- factor(my_data$asthma)

## males and females in separate panels
ggplot(my_data, aes(col=asthma, x=asthma, y=Weight))+
  geom_boxplot()+
  facet_wrap(~gender)

If, on the other hand, you wish to compare gender, for each condition of asthma, then you can do

# asthma 0 and asthma 1 in separate panels
ggplot(my_data, aes(col = gender, x=gender, y=Weight))+
  geom_boxplot()+
  facet_wrap(~asthma)

To obtain violins instead of boxplots in the above code just replace geom_boxplot() with geom_density().

Option 2 (updated question) The idea is to compute a $2 \times 2$ frequency table for asthma vs gender for individuals having weight less than or equal to a given value. By varying this weight value, one obtains a series of $2\times 2$ frequency tables, which can then be plotted by, e.g. a barplot. The R code below shows how to perform these calculations and how to plot the results. Again, depending on the emphasis one can split the windows by gender or by asthma condition. I'm choosing 20 equally-distant weight values but you can change it at your will.

my_data$asthma <- factor(my_data$asthma, labels = c("no", "yes"))

# the number of weight values, e.g. the number of bars (see below)
n <- 20
# cut weights in n intervals
weight_cuts <- seq(from= range(my_data$Weight)[1],
                   to= range(my_data$Weight)[2], 
                   length.out=n)

hist(my_data$Weight)
points(weight_cuts, y=rep(0,n), pch=20)

# computes the distribution of asthma conditionally for unit
# with weight weight_thersh or lower.  
wcond_table <- function(weight_thersh, dataset) {
  dd <- dataset[dataset$Weight<= weight_thersh,]
  print(table(dd[,-1]))
  matrix(table(dd[,-1]), ncol=2, nrow=2)
}

wcond_table(weight_cuts[3], my_data)

# build the dataset of conditional counts
dd <- t(sapply(weight_cuts, wcond_table, 
               dataset = my_data))
dd2 <- data.frame( noasthma_female= dd[,1], 
                   noasthma_male = dd[,2], 
                   asthma_female = dd[,3],
                   asthma_male = dd[,4],
                   weight = weight_cuts)
# reshape it to be suitable for ggplot2
dd_long <- reshape(dd2, direction = "long", 
                   varying = list(1:4),
                   timevar = "asthma_gender",
                   times = c("noasthma_female", "noasthma_male",
                             "asthma_female", "asthma_male"),
                   v.names = "counts")
dd_long$gender = factor(c(rep("female", 20), rep("male", 20),
                          rep("female", 20), rep("male", 20)))
dd_long$asthma = factor(c(rep("noasthma", 40), 
                            rep("asthma", 40)))

Now, for the plots, if you want a window for every kind of gender, that is if you want to compare asthma vs non-asthma within a given gender you may consider this

ggplot(dd_long, aes(y = counts, 
                    x = weight, 
                    fill = asthma))+
  geom_bar(stat = "identity")+
  facet_wrap(~gender)

On the other hand, if you want to compare gender for a given condition of asthma, then you may consider this.

ggplot(dd_long, aes(y = counts, 
                    x = weight, 
                    fill = gender))+
  geom_bar(stat = "identity")+
  facet_wrap(~asthma)

As per my last comment, if you want to leverage yourself from binning the weight, you'll have to use a logistic regression or a GAM as suggested in another answer. But, since this is an exploratory data analysis, I do not see any issue with what you call p-hacking and in my opinion, "hand-made" calculations give precious insight into the data.

Answer 3

You have four groups of interest: males with asthma; males without asthma; females with asthma; and females without asthma. For each group you want to visualise the distribution of their weights. Given this objective, I would recommend you consider using a violin plot showing the densities for the four groups next to each other, with appropriate labelling/colouring of the densities to make them easily discernible. Please read the linked material to review what a violin plot does and how to construct it. You can find instructions on making a violin plot in ggplot2 here.