4
$\begingroup$

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")

enter image description here

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")

enter image description here

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:

enter image description here

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") 

enter image description here

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!

$\endgroup$
4
  • 2
    $\begingroup$ Hi! Which are these "given weight ranges" are they fixed? Furthermore, what do you mean by "mitigate the risk of p-hacking" ? Here you are doing EDA, I don't see any issues with p-hacking. $\endgroup$
    – utobi
    Jan 25, 2023 at 21:51
  • $\begingroup$ @ utobi: Hi! Thank you for your reply! By "given weight ranges" - consider "percentiles" (e.g. how many men in the bottom 10% of all weights have asthma compared to how many men in the top 10% of all weights have asthma". $\endgroup$
    – stats_noob
    Jan 25, 2023 at 21:55
  • $\begingroup$ I have heard that making arbitrary ranges within the data can lead to p-hacking ...I have been trying to see if there are ways to avoid this (e.g. stats.stackexchange.com/questions/602844/…) $\endgroup$
    – stats_noob
    Jan 25, 2023 at 21:55
  • $\begingroup$ Since there are no fixed specific values of the weight variable conditionally on which you want to compute asthma proportions, I think logistic regression is the most reasonable choice. $\endgroup$
    – utobi
    Jan 26, 2023 at 11:14

3 Answers 3

2
+50
$\begingroup$

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: Plot with GAM smoothed curve

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.

$\endgroup$
11
  • $\begingroup$ @ Björn: thank you so much for your answer! I really like the approach you have taken! Its very unique and creative, and I think it is answering the exact question I am interested in! $\endgroup$
    – stats_noob
    Jan 27, 2023 at 3:06
  • $\begingroup$ Two questions I had about this: $\endgroup$
    – stats_noob
    Jan 27, 2023 at 3:06
  • $\begingroup$ Question 1: Why have you chosen a GAM/Spline model? Could a regular logistic regression have worked as well? For example: model <- glm(asthma ~ Weight, data = my_data, family = binomial()) $\endgroup$
    – stats_noob
    Jan 27, 2023 at 3:08
  • $\begingroup$ Question 2: What was the significance of using "y=asthma+case_when(asthma==1 & gender=="female" ~ -0.05, asthma==0 & gender!="female" ~ 0.05, TRUE ~ 0)"? Why did you not just use "y = asthma"? $\endgroup$
    – stats_noob
    Jan 27, 2023 at 3:10
  • $\begingroup$ Thank you so much for your everything! $\endgroup$
    – stats_noob
    Jan 27, 2023 at 3:11
6
$\begingroup$

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)

enter image description here

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)

enter image description here

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)

enter image description here

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)

enter image description here

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.

$\endgroup$
3
  • 1
    $\begingroup$ @ utobi: thank you for your answer! $\endgroup$
    – stats_noob
    Jan 25, 2023 at 21:58
  • 1
    $\begingroup$ @ utobi: thank you so much for this update! $\endgroup$
    – stats_noob
    Jan 28, 2023 at 5:37
  • $\begingroup$ @utobi: very good answer. However, your violin plot command , yields strange results i.e. geom_density() in place of geom_boxplot(). Any insight as to why this might be? I just a side-by-side plot with straight up and down lines in the plots $\endgroup$
    – Andy
    Jan 30, 2023 at 8:25
1
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ @ Ben: Thank you for your reply! Originally, I was interested in visualizing the "asthma rates"- but now I realize that this might not be "objectively" possible given the data I have? $\endgroup$
    – stats_noob
    Jan 24, 2023 at 4:35
  • $\begingroup$ If you want to visualise asthma rates of males and females then it would be a simple bar chart with two bars. $\endgroup$
    – Ben
    Jan 24, 2023 at 4:56
  • $\begingroup$ @ Ben : thank you for your reply! I am interested in learning about ways to visualize the proportion of people who have asthma in given weight ranges ... and mitigate the risk of engaging in p-hacking while creating these weight ranges $\endgroup$
    – stats_noob
    Jan 25, 2023 at 21:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.