3
$\begingroup$

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.

Improve this question
$\endgroup$
  • $\begingroup$ I have created a working solution for the problem myself. I'll probably add an answer later. (Although, I would still love to see what other people would do) $\endgroup$ – Eff Mar 14 '19 at 8:11
  • $\begingroup$ You have 12 equations in only $2+6$ unknowns so unless you know the exact true population proportions your system of equations will not have a solution. The problem would be more interesting (and perhaps more realistic) if you instead only had multinomial data on the proportions in which case you could estimate the parameters (6 thresholds, $\Delta$ and $\rho^2$) by maximum likelihood. $\endgroup$ – Jarle Tufto Mar 14 '19 at 10:41
  • $\begingroup$ @JarleTufto I'm not sure I understand your point exactly. I ended up finding a "solution" (at least, I think it works?). It's given as an answer below. It seems to work decently--for this data set at least. If I have more equations than unknowns, that should simply be more information, and I can get a solution which is the "best fit" (won't fit exactly with each point, necessarily). $\endgroup$ – Eff Mar 14 '19 at 12:02
  • $\begingroup$ I'm not trying to find values that will fit exactly, necessarily. Because I'm looking for a method that makes estimates based on data which can be noisy. $\endgroup$ – Eff Mar 14 '19 at 12:05
  • 2
    $\begingroup$ Exactly, the data is noisy. How to best take this into account depends on the exact details of the sampling process. If the proportions are based on counts of number of individuals within each of the 6 categories in a sample taken from the population, then what I'm saying is that these counts follow a multinomial distribution. If this is the case, I can post an answer. $\endgroup$ – Jarle Tufto Mar 14 '19 at 12:34
3
$\begingroup$

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
Improve this answer
$\endgroup$
  • $\begingroup$ Thanks a lot for the answer (+1)! I will try and implement it later, and test it out. $\endgroup$ – Eff Mar 14 '19 at 13:49
  • $\begingroup$ Hmm.. Perhaps I'm missing something, but the convergence doesn't appear to be very stable. I try writing exactly your code, and get the same result as you. But if I try and put in my data, i.e. I set z <- rnorm(n, c(161,175)[sex], c(7.1,8.5)[sex]) then, with $n = 500$, I get $\Delta = 1.8978$ and $\rho = 1.9822$. With $n = 2000$, I get $\Delta = 1.8800$ and $\rho = 1.3274$. With $n = 5000$, I get $\Delta = 1.7998$ and $\rho = 1.0331$. With $n = 10000$, I get $\Delta = 1.9752$ and $\rho = 1.2202$. With $n = 20000$, I get $\Delta = 0.4666$ and $\rho = 0.4780$. Am I misunderstanding something? $\endgroup$ – Eff Mar 14 '19 at 21:01
  • $\begingroup$ These seems like wildly different estimates. I am still using set.seed(1) for easy replicability. But I am correct in setting z <- rnorm(n, c(161,175)[sex], c(7.1,8.5)[sex]) to set $\mu_F = 161$ and $\mu_M = 175$ and $\sigma_F = 7.1$ and $\sigma_M = 8.5$, right? $\endgroup$ – Eff Mar 14 '19 at 21:03
  • $\begingroup$ Yes, the code may not be very robust. There are awkward constraints on the parameters so you may want to do the optimisation with respect to a different parameterisation of the model. Also, reasonable starting values for the thresholds can be computed directly from the counts in males. $\endgroup$ – Jarle Tufto Mar 14 '19 at 21:43
  • $\begingroup$ Btw, this is almost like ordinal regression (see e.g. function MASS::polr) except that the variance of the underlying latent normally distributed variable is not constant. $\endgroup$ – Jarle Tufto Mar 14 '19 at 21:45
2
$\begingroup$

(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):

Plot of best fit

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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