Women’s shoe sizes

The research question that we wish to answer is the following: what is the proportion of women with shoe size greater or equal to UK size 9? We want to know how unusual this is. Is it more or less unusual than the proportion of all trans people in the general population?

Katy demonstrated that the number of UK cis women with shoe size greater or equal to UK size 9 is overwhelmingly larger than the number of all trans women in the UK in total.

Katy found a paper with some data on the range of women’s shoe size:

https://www.researchgate.net/publication/232851673_Foot_shape_of_older_people_Implications_for_shoe_design

The sample size is a bit low: 154 women were reported as being surveyed. However, data for four of the women appears to be missing, and Katy hypothesised that either the measurements were not recorded, or the measurements fall outside the range of histogrammed data. It also cannot be absolutely ruled-out that the authors screwed-up – we can’t tell, because they don’t provide the raw data or their analysis code.

The small sample size isn’t actually a problem, provided we recognise that the statistic of interest – the proportion of women with shoe size greater or equal to UK size 9 – is a point estimate which has an uncertainty. We have methods for quantifying that uncertainty that we can use to determine a range of plausible values that the statistic might take.

More of a problem is that the data might not be representative of the general population. The paper is about the foot shape of older people. The first line of the abstract states: “Foot problems and deformities, which are highly prevalent in older people, may affect foot anthropometrics and result in older people having difficulty finding shoes that fit.”

Moreover, the methods section reports the following: “The sample included 158 men (age: 71.1±6.3\(\,\)years; height: 1.73±0.1\(\,\)m; mass: 84.2±13.8\(\,\)kg; BMI: 28.2±3.6\(\,\)kg/m²) and 154 women (age: 71.7±7.0\(\,\)years; height: 1.59±0.1\(\,\)m; mass: 72.4±16.3\(\,\)kg; BMI: 28.6±6.1\(\,\)kg/m²).”

The age of those sampled, and the suggestion that they may have foot deformities that affect shoe size, suggests that the people who were sampled may not be representative of the general population. Might this affect Katy’s conclusion? To answer this question, we found another (much larger) data set, which can be found here:

http://stat.pugetsound.edu/hoard/projectDetails.aspx?id=53

The data is described as follows:

"Each semester in all of the Introductory Statistics sections at the University of Puget Sound, a survey is given to the students during the first week of class. The survey is voluntary and is used as an example data set throughout the class.

The data set given here is a compilation of much of the data collected in the period from the Fall 2002 semester through the Spring 2008 semester. Values that have been determined to be incorrect (such as 8-foot tall students) have been removed from this data set."

Let’s load and examine the data; there are several columns, but we’re only interested in student gender and their shoe size (and, later, height), which we assume to be self-reported:

# install.packages("tidyverse") # Uncomment to install
suppressPackageStartupMessages( # STFU!
  suppressWarnings( # STFU!
    require(tidyverse, quietly = TRUE)
  )
)
# Note: the printed link on the webpage is dead, but the hyperlink still works!
suppressMessages({ # STFU!
  raw <- read_csv("http://stat.pugetsound.edu/hoard/datasets/introStats2002-2008.csv")
})

raw %>% select(gender, shoeSize, height) -> dat # Drop some columns
dat %>% select(gender, shoeSize) %>% glimpse() # Show some stats

## Observations: 2,433
## Variables: 2
## Error in get(genname, envir = envir) : object 'vec_ptype2' not found
## $ gender   <chr> "M", "M", "M", "M", "M", "M", "F", "M", "F", "F", "F"...
## $ shoeSize <dbl> 12.0, 11.0, 13.0, 12.0, 11.0, 10.5, 5.5, 12.0, 8.5, 7...

We have no information on whether there are any trans people in the data, but we assume that if there are any, they’re represented in the same proportion as occurs in the general population, so they won’t skew the results. Interestingly, there are ten students with no gender recorded! (What could this mean?)

table(dat$gender, useNA = "always")

## 
##    F    M <NA> 
## 1394 1029   10

In the data, US shoe sizes are reported, which are different for men and women. To put these shoe sizes on the same scale, we convert them to UK shoe sizes, which are nominally the same for men and women. To do this, we’ll need to get a conversion table from the internet:

# We need a shoe size conversion table; let's scrape one from the internet:
suppressMessages( # STFU!
  require(rvest, quietly = TRUE)
)
webpage <- read_html("https://www.blitzresults.com/en/shoe-size-conversion/")
tbls <- html_nodes(webpage, "table") # Scape tables from the HTML
head(tbls) 

## {xml_nodeset (6)}
## [1] <table id="tablepress-334" class="tablepress tablepress-id-334 table ...
## [2] <table id="tablepress-332" class="tablepress tablepress-id-332 table ...
## [3] <table id="tablepress-333" class="tablepress tablepress-id-333 table ...
## [4] <table id="tablepress-351" class="tablepress tablepress-id-351 table ...
## [5] <table id="tablepress-351-no-2" class="tablepress tablepress-id-351  ...
## [6] <table id="tablepress-352" class="tablepress tablepress-id-352 table ...

tbls[4] %>% # It's the 4th table we need
  html_table(fill = TRUE) %>% # Returns a list of data frames
  first() -> conversion # Take the first element in the list
conversion # Ta-da!

From this we see that the conversion rule appears to be:

• US Women’s to UK: subtract 2
• US Men’s to UK: subtract 0.5

However, other tables like this and this indicate that we should subtract 1 from the US Men’s size to get to the equivalent UK size. Another website indicates a half size difference between UK sizes for male and female shoes. Yet another suggests subtracting 2.5 from the US Women’s size to the equivalent UK size. Different manufacturers use different lasts to construct their shoes, and sizing may vary accordingly. Nevertheless, despite this inconsistency, we note that the difference between conversion rules is at most only half a size in UK measure.

Here are the counts for each shoe size, in US measure, grouped by gender (“NA” represents missing data):

with(dat, table(gender, shoeSize, useNA = "always"))

##       shoeSize
## gender   3 3.5   4   5 5.5 5.75   6 6.25 6.5   7 7.5   8 8.25 8.5 8.75   9
##   F      1   1   2  18  19    1  67    1  68 140 188 229    1 209    1 187
##   M      0   0   0   0   0    0   1    0   2   1   3  18    0  12    0  54
##   <NA>   0   0   0   0   0    0   0    0   0   1   1   2    0   2    0   0
##       shoeSize
## gender 9.25 9.5  10 10.5  11 11.5  12 12.5  13 13.5  14  15  16  17 <NA>
##   F       0  77 114   22  32    1  10    0   2    0   0   0   0   0    3
##   M       1  70 111  145 172   75 166   27 106    5  40   8   3   1    8
##   <NA>    0   1   3    0   0    0   0    0   0    0   0   0   0   0    0

Not everyone reported their shoe size, so we’ll remove the individuals who didn’t from the data:

dat %>% filter(!is.na(shoeSize)) -> dat

We’re not sure how to convert a person’s US shoe size to a UK shoe size in the case that there is no gender marker recorded for them, so for the sake of simplicity we’ll also remove these individuals from the data:

dat %>% filter(!is.na(gender)) -> dat

Now we convert everyone’s sizes to UK shoe size:

# Now we need to use our conversion rule to convert everyone's shoe size to UK,
# which nominally isn't different for men and women:
dat %>% 
  mutate(UKshoeSize = shoeSize - if_else(gender == "F", 2, 0.5)) %>% 
  {.} -> dat

We’ll now only refer to shoe size in UK measure, even though the individuals in the sample reported their shoe size in US measure.

Here’s the frequency table for our data with shoe size in UK measure:

with(dat, table(gender, UKshoeSize))

##       UKshoeSize
## gender   1 1.5   2   3 3.5 3.75   4 4.25 4.5   5 5.5   6 6.25 6.5 6.75   7
##      F   1   1   2  18  19    1  67    1  68 140 188 229    1 209    1 187
##      M   0   0   0   0   0    0   0    0   0   0   1   2    0   1    0   3
##       UKshoeSize
## gender 7.5   8 8.5 8.75   9 9.5  10 10.5  11 11.5  12 12.5  13 13.5 14.5
##      F  77 114  22    0  32   1  10    0   2    0   0    0   0    0    0
##      M  18  12  54    1  70 111 145  172  75  166  27  106   5   40    8
##       UKshoeSize
## gender 15.5 16.5
##      F    0    0
##      M    3    1

Let’s plot the data:

ggplot(dat) +
  geom_histogram(
    aes(x = UKshoeSize, fill = gender), 
    alpha = 0.6,
    position = "identity",
    binwidth = 0.5, 
    center = 0
  ) +
  scale_x_continuous(
    breaks = seq(0, 17.5, 0.5)
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  xlab("UK shoe size") +
  ylab("Count") +
  ggtitle("Distribution of UK shoe size")

We see the data is distributed very similarly to the data used in Katy’s analysis.

In this sample, what proportion of women have shoe size greater or equal to UK size 9? To answer this, we filter the data to select only the women and calculate the proportion. We’ll make a function to do this that we can use again later:

dat %>% filter(gender == "F") -> women # Data for just the women
# A function that returns the proportion of women with size >= 9
test.size <- 9 # The UK shoe size we are testing
prop <- function(data, indices) {
  d <- data[indices,] # Allows boot to select sample
  sum(d$UKshoeSize >= test.size, na.rm = TRUE) / length(d$gender) # Proportion
}
# Apply this function to all rows of the data sample (just women):
result <- prop(women, 1:nrow(women)) # Result for this data sample
scales::percent(result)

## [1] "3.24%"

Therefore we see that about 3% of women in this sample have a shoe size greater or equal to UK size 9. However, we’d like to know what the uncertainty on this point estimate is. Could we have just been lucky with our sample data? If we randomly selected the same number of women from the same population, we would expect to get a similar result, but not the exact same result. Our results have a degree of uncertainty to them which is due to random variation.

Could it be that the true proportion of women with shoe size greater or equal to UK size 9 is much smaller in the general population? To answer this, it might seem that we would need data for the entire population. However, there are techniques that allow us to do meaningful statistical inference based on the data that we already have. One of these methods is the classical nonparametric bootstrap, and it is this method that we shall use now. This is a statistical method that was developed in the 1970’s. The paper presenting the method has over 42 thousand citations. The bootstrap is a very common tool for statistical inference.

We assume that the sample that we have is representative of the general population, and we have no reason to suspect that this is not the case. We have no reason to suspect that the distribution of shoe sizes of people who go to college is significantly different than those who do not go to college. What we do during bootstrapping is sample repeatedly (with replacement) from the data that we have and calculate our statistic of interest on each sample. We then have a distribution of the statistic of interest, which we can then use to estimate a range of plausible values for it. Let’s do that now, using the function we made earlier to simplify the process, and plot the resultant distribution:

# install.packages("boot") # Uncomment to install
require(boot, quietly = TRUE)
set.seed(42) # Set random number seed for reproducibility
n_resamples <- 1000000 # Number of resamples to make
# Calculate proportions on many bootstrap resamplings:
women.boot <- boot(data = women, statistic = prop, R = n_resamples)
women.boot

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = women, statistic = prop, R = n_resamples)
## 
## 
## Bootstrap Statistics :
##       original       bias    std. error
## t1* 0.03235083 8.621136e-06 0.004743474

props <- data.frame("t" = women.boot$t) # Proportions for each resampling

ggplot(props) +
  geom_histogram(
    aes(x = t, y = ..density..), 
    fill = "skyblue",
    bins = 25
  ) +
  scale_x_continuous(labels = scales::percent) +
  theme_minimal() +
  xlab("Percentage of women with shoe size greater or equal to UK size 9") +
  ylab("Count") +
  ggtitle("Sampling distribution of percentage of women with large shoe sizes")

We can then use this distribution to estimate a 99% confidence interval for the proportion of women with shoe sizes greater or equal to UK size 9:

q <- quantile(props$t, probs = c(0.005, 0.995)) # Quantiles for a 99% CI
scales::percent(q) # Express 99% CI as rounded percentages

## [1] "2.08%" "4.53%"

We can therefore be 99% confident that the proportion of women in the general population with shoe size greater or equal to UK size 9 is between approximately 2% to 4.5%. A survey¹ in 2016 estimated that 0.6% of US adults identify as transgender. The data in this analysis and the survey data corresponds to US adults, so we know that comparisons between the two are appropriate. Even if the proportion of women in the population with shoe size \(\geq\) UK size 9 is only 2%, this is still much larger than the percentage of all trans people in the US. Moveover, there is no reason to believe that this is not true more generally.

Clearly the percentage of all trans people is significantly smaller than the estimated proportion of women in the general population with shoe size greater or equal to UK size 9.

Consequently, Katy’s original conclusion is shown to be correct. We just accounted for the statistical uncertainties to give that conclusion a much more robust foundation.

At this point it might be tempting to estimate what women’s shoe size is so rare that the percentage of women with shoe size equal or greater to it is similar to the percentage of all trans people in the general population. However, this isn’t a question we can answer with this data because there are not enough data for women with large feet in this sample. Nevertheless, it is interesting to note that there are two women with UK size 11 in this data.

dat %>% 
  mutate(heightcm = 2.54 * height) %>% # Convert inches to cm
  {.} -> dat

dat %>% 
  drop_na(heightcm) %>% # Drop individuals with missing height
  drop_na(gender) %>% # Drop individuals with no reported gender
  ggplot() +
  geom_histogram(
    aes(x = heightcm, fill = gender), 
    alpha = 0.6,
    position = "identity",
    binwidth = 2.5, 
    center = 0
  ) +
  scale_x_continuous(
    breaks = seq(110, 220, 10)
  ) +
  theme_minimal() +
  xlab("Height (cm)") +
  ylab("Count") +
  ggtitle("Distribution of height")

dat %>% 
  filter(gender == "F") %>% # Select women
  drop_na(height) %>% # Drop women with missing heights
  {.} -> women # Data for just the women

# A function that returns the proportion of women with height >= 6 ft
test.height <- 6*12 # Height in inches (6 ft)
prop <- function(data, indices) {
  d <- data[indices,] # Allows boot to select sample
  sum(d$height >= test.height, na.rm = TRUE) / length(d$gender) # Proportion
}
# Apply this function to all rows of the data sample (just women):
result <- prop(women, 1:nrow(women)) # Result for this data sample
scales::percent(result)

## [1] "1.84%"

# Calculate proportions on many bootstrap resamplings:
women.boot <- boot(data = women, statistic = prop, R = n_resamples)
women.boot

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = women, statistic = prop, R = n_resamples)
## 
## 
## Bootstrap Statistics :
##       original       bias    std. error
## t1* 0.01839588 1.936718e-06 0.003644574

props <- data.frame("t" = women.boot$t) # Proportions for each resampling

ggplot(props) +
  geom_histogram(
    aes(x = t, y = ..density..), 
    fill = "aquamarine",
    bins = 25
  ) +
  scale_x_continuous(labels = scales::percent) +
  theme_minimal() +
  xlab("Percentage of women as tall as 6 ft or taller") +
  ylab("Count") +
  ggtitle("Sampling distribution of percentage of tall women")

q <- quantile(props$t, probs = c(0.005, 0.995)) # Quantiles for a 99% CI
scales::percent(q) # Express 99% CI as rounded percentages

## [1] "0.96%" "2.87%"

dat %>% 
  drop_na(heightcm) %>% # Drop individuals with missing height measurement
  ggplot() +
  geom_point(
    aes(
      x = heightcm,
      y = UKshoeSize,
      fill = gender,
      colour = gender
    ),
    alpha = 0.6
  ) +
  scale_x_continuous(
    breaks = seq(110, 220, 10)
  ) +
  scale_y_continuous(
    breaks = seq(0, 17.5, 1)
  ) +
  theme_minimal() +
  xlab("Height (cm)") +
  ylab("UK shoe size") +
  ggtitle("UK shoe size versus height")