Male-female height difference estimation from thresholds

Question

Introduction. Assume that two populations $A$ and $B$ are distributed with normal distributions $N(\mu_A, \sigma_A^2)$ and $N(\mu_B, \sigma_B^2)$.

This is a general problem, but as an example, I will consider height distributions of men and women. Assume that the heights of men and women are normally distributed with the following means and standard deviations:

$$\begin{align} & \textbf{Mean} & \textbf{St.d.}\\ \textbf{Female}\quad & 161\,\text{cm} & 7.1\,\text{cm}\\ \textbf{Male}\quad& 175\, \text{cm} & 8.5\,\text{cm} \end{align}$$

Importantly, I don't know these parameter values, but I wish to infer/estimate them (Read on). Instead, what I have is a table of the proportion of men and women exceeding some (unknown) threshold. For example, I have manually created some thresholds. The following table contains the proportion of women and men, respectively, who exceed a certain threshold of height:

$$\begin{align} \textbf{Threshold} &\quad \textbf{Female} & \textbf{Male}\\ \geq158\,\text{cm} &\quad0.6637 & 0.9772\\ \geq 162\,\text{cm} & \quad 0.4440 & 0.9369\\ \geq 166\,\text{cm} & \quad 0.2406 & 0.8552 \\ \geq 170\,\text{cm} & \quad 0.1025 & 0.7218\\ \geq 174\,\text{cm} & \quad 0.0336 & 0.5468\\ \geq 178\,\text{cm} & \quad 0.0083 & 0.3621\end{align}$$

I only have the actual proportion values (0.6636, 0.9772, ... etc) but I don't have the threshold values (158 cm, 162 cm, etc) as they are unknown.

Question. Given such a threshold table, how do I find the group mean difference (in SD units) and variance ratio? (Assuming that they are normally distributed)

Results. Given I have manually constructed the data, I can say what the results should approximately be. If we standardize the female distribution, then the mean difference is

$$\Delta := \frac{\mu_M-\mu_F}{\sigma_F} = \frac{175-161}{7.1} \approx 1.972 \quad\text{female standard deviation units} $$

The variance ratio is $$\rho^2 := \frac{\sigma_M^2}{\sigma_F^2} = \frac{8.5^2}{7.1^2} \approx 1.433$$

Which are the two values that I'm looking for (mean difference and variance ratio).

---

Thoughts. My thinking is to first reduce the problem from four unknown variables ($\mu_A$, $\mu_B$, $\sigma_A^2$, $\sigma_B^2$) to a problem of two unknowns. This can be done by standardizing one group to be distributed with mean 0 and variance 1. Then I am simply looking for the mean difference (say, $\Delta$) and the variance ratio.

Then I probably have to calculate some density function or cumulative density function for the difference, and make a least squares optimization problem with the table (There can be noise in the data, hence I am thinking a least squares minimization). But I'm unsure exactly how.

Answer 1

Here is some R code simulating data, fitting a model by standard maximum likelihood methods assuming that the counts follow a multinomial distribution, and comparing observed and fitted expected observed counts.

# Simulate the data generating process
set.seed(1)
n <- 500
sex <- rep(c(1,2),each=n/2)
z <- rnorm(n, c(170,180)[sex], c(10,15)[sex])
threshold <- seq(150,190,len=5)
group <- 0
for (j in 1:length(threshold))
  group <- group + ifelse(z > threshold[j], 1, 0)
x <- table(sex, group)

# Function computing 2x6 matrix of multinomial probabilities from model parameters
pfn <- function(par) {
  delta <- par[1]
  rho <- par[2]
  threshold <- c(-Inf,par[-(1:2)],Inf)
  rbind(diff(pnorm(threshold)),           
        diff(pnorm((threshold-delta)/rho)))
}

# Negative log likelihood assuming that counts have a multinomial distributions
nll <- function(par, x) {
  p <- pfn(par)
  ll <- 0
  for (i in 1:2)
    ll <- ll + dmultinom(x[i,], prob=p[i,], log=TRUE)
  -ll
}

# Initial guess of parameter values
start <- c(delta=0,rho=1,threshold=seq(-1,1,len=5))
# Maximise the log likelihood
fit <- optim(start, nll, x=x, hessian=TRUE)

As you can see, the parameter estimates are reasonably close to the true values but this is not so surprising given the large sample of $n=500$. You may also want to test goodness-of-fit perhaps using the Pearson's chi-square test.

> # Estimates and asymptotic standard errors
> est <- cbind(fit$par,sqrt(diag(solve(fit$hessian))))
> colnames(est) <- c("Estimate", "SE")
> est
               Estimate         SE
delta       1.023136895 0.13996277
rho         1.579971333 0.12695298
threshold1 -2.016226101 0.16087921
threshold2 -1.052349248 0.09170107
threshold3  0.002726207 0.07482051
threshold4  1.001152575 0.08985068
threshold5  1.940070686 0.14024179
> 
> # Expected counts under fitted model
> pfn(fit$par)*apply(x,1,sum)
     threshold1 threshold2 threshold3 threshold4 threshold5          
[1,]   5.472042   31.10788   88.69198   85.13397   33.04774  6.546387
[2,]   6.799301   16.82235   41.17612   58.81452   56.17775 70.209964
> # Observed counts
> x
   group
sex  0  1  2  3  4  5
  1  4 32 97 77 33  7
  2  9 16 36 67 55 67

Answer 2

^{(I created a solution to the problem myself, so I give it as an answer here for anyone who is interested.)}

Method

First, we have the data that we wish to estimate group differences from. These indicate the proportions of women and men, respectively, that exceed a certain threshold. For example, when $66.36\%$ of women are taller than some threshold, then $97.72\%$ of men are taller than that same threshold. Solution code is given in R. Here is the data from the original question:

P_F = c(0.6636, 0.4440, 0.2406, 0.1025, 0.0336, 0.0083)
P_M = c(0.9772, 0.9369, 0.8552, 0.7218, 0.5468, 0.3621)

First, we standardize the female height distribution. They have mean 0 and standard deviation 1. From this, we can calculate the underlying threshold values $t_1, t_2, ..., t_n$ (where $n = 6$ is the number of thresholds, i.e. the length of one of the vectors above). We do this by solving for $t_i$ such that

$$\int_{t_1}^\infty f(t;\mu = 0,\sigma = 1)\,\mathrm{d}t = 0.6636$$ where $f(t; \mu, \sigma)$ is the normal density function with mean $\mu$ and standard deviation $\sigma$. We do this for all values in P_F (this was just the first one). Here is the code:

f0 = function(x){ # Standard normal density function
     1/sqrt(2*pi)*exp(-x^2/2)
}
# Find threshold values
thresholds = rep(0, length(P_F))
for(i in 1:length(P_F)){
     F0_temp = function(t){
     -P_F[i] + integrate(f0,t,Inf)$value
     }
thresholds[i] = uniroot(F0_temp, lower = -10, upper = 10)$root
}

Once we have the threshold values $t_1, t_2, ..., t_n$ we want to find values $\Delta$ and $\rho$ such that $f(t,\Delta,\rho)$ is the best fit to the data P_M. First we create the function that we wish to fit

# Create function to fit
F_fit = function(t,Delta, rho){
     f0_full = function(x, mu = Delta, sigma = rho){
          1/sqrt(2*pi*sigma^2)*exp(-(x-mu)^2/(2*sigma^2))
     }
     outp = rep(0, length(t))
     for(i in 1:length(outp)){
          outp[i] = integrate(f0_full,t[i],Inf)$value
     }
     return(outp)
}

Next we find the best fit

# Fit using nonlinear least squares
fitLS = nls(P_M ~ F_fit(thresholds, Delta, rho), start = list(Delta = 1.6, rho = 1.2))

The found parameters are then given by

# Store parameters
Delta = unname(coef(fitLS)[1]) # Mean difference (in SD units)
rho = unname(coef(fitLS)[2]) # Standard deviation ratio

Results

We can check if the results come close to the expected values (correct values), with $\Delta$ being the mean difference and $\rho^2$ being the variance ratio.

> Delta
[1] 1.972079
> rho^2
[1] 1.434806

As we can see, they are very close to the expected values $1.972$ and $1.433$ (see original question). The method seems to perform decently (at least in this case).

Plot

I now create a plot with the following code

# Create plot
linspac = seq(0,1, length.out = 100) # Line of equality
plot(linspac, linspac, type = "l", main = "Best fit", xlab = "Proportion of females exceeding threshold", ylab = "Proportion of males exceeding threshold")
points(P_F, P_M, col = "red") # Data points

tvals = seq(-5,5,length.out = 100)
xvals = rep(0,length(tvals))
yvals = rep(0,length(tvals))

for(i in 1:length(tvals)){ # Create curve
     xvals[i] = F_fit(tvals[i], 0, 1)
     yvals[i] = F_fit(tvals[i], Delta, rho)
}

lines(xvals, yvals, col = "blue") # Plot curve

Here is the plot from the code (annotated in paint):