A nonlinear effect plus measurement error in x becomes close to linear: A simulation study

It’s time for another one of Bob’s favorite sort of post, which is when I use simulation to demonstrate a statistical principle.

1. Flogging a dead horse

The issue in question came up the other day in comments to my post, This is how we do modern frequentist statistics: Using fake-data simulation to understand what can happen in a study. I discussed a famously flawed paper from 2007 that had claimed to find a strong relationship between parents’ physical attractiveness and the sex of their children. It had been apparent at the time (ok, apparent to quantitative social scientists, not to the hapless Freakonomics team) that for statistical reasons the study was hopeless—it was based on a sample of about 3000 cases, and to detect any reasonable effect it would be necessary to have orders of magnitude more people in the study: maybe with 3,000,000 cases you’d be able to see some signal.

That n = 3000 is not enough can be a surprise: after all, 1500 people are enough for a national poll, so why is 3000 not enough to study some sociobiological pattern? The answers to why 3000 is not nearly enough are:
1. Effect size;
2. The standard error scales like 1/sqrt(n);
3. Noisy measurement.

Regarding the first point: a survey of 1500 people is enough (if it’s a representative sample) to estimate national opinion to within +/-3 percentage points. Which is great. But, based on the literature on sex ratio variation, we can be pretty sure that any differences based on attractiveness will be less than half a percentage point, probably less than a tenth of a percentage point.

Regarding the second point: to go from an uncertainty of 3 percentage points to an uncertainty of 0.03 percentage points will require the sample size to increase by a factor of 100^2. So, instead of 3000, you’d need 30,000,000. But, hey, maybe the underlying difference is larger than you think, which is why I’m saying you might have a shot with n = 3 million.

Regarding the third point: “attractiveness” is not a clearly defined construct, and the measurement (based on a survey interviewer’s one-time interaction with the respondent) is noisy. This will attenuate any differences.

It’s that last point I’ll demonstrate in this post. Using simulation.

2. What happens to a nonlinear effect when you add random noise to the predictor?

My immediate motivation here came from a discussion in comments. In my simulation, I estimated the association in the data between attractiveness and sex ratio by fitting a linear regression. But commenter Sandro wrote:

Not sure linear regression makes the most sense for every pattern because it imposes linearity. what if, for example, there is threshold and the pattern is flat to the left of the threshold and flat at a higher level to right of the threshold?

I replied:

This attractiveness measure is noisy. So even if there’s an underlying threshold effect if measured based on some latent or idealized “true attractiveness,” the relationship with measured attractiveness would be pulled toward linearity by the measurement (and also any coefficient would be attenuated, which is one reason I’m sure that any true underlying difference would be tiny).

Sandro responded:

I don’t dispute the idea that linearity is a useful default, but it is not necessarily the most powerful test. I don’t understand the “pulled towards linearity” argument. Here is where I am coming from: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=7149560

And I replied:

I guess the easiest way to show this would be to construct a simulation which starts with some relationship between latent or “true” attractiveness, and then adds random error to latent attractiveness to create measured attractiveness. I think that, if the underlying relationship is monotonic and then random error is large, that linear regression will be a robust and efficient approach. Again, this is moot for this particular example, but more generally it’s an interesting question.

Moot it may be, but I decided to do the simulation anyway, because why not. It’s a better morning’s work than submitting 258 papers to SSRN.

3. Setting up the model

This simulation is a little different than what came before. For my earlier post I simulated hypothetical future data from the beauty-and-sex-ratio study. For my new post I’m simulating data from a hypothetical strong relationship between two variables. Here I’m entirely focused on the question of what happens to a nonlinear (but monotonic) pattern when independent measurement error is added to the predictor. To make the patterns clear, I will assume a large effect; indeed, I’m only concerned here with the effect size, not its variation.

The one thing I’ll preserve from the previous simulation is the distribution of attractiveness measurements in the sample: 2%, 5%, 45%, 37%, and 11% in the five attractiveness categories.

I want to set up a model where there’s an assumed true nonlinear relationship between some outcome of interest and a continuous underlying latent attractiveness variable.

The nonlinear model isn’t hard to come up with. The hard thing is creating a distribution for the latent attractiveness, under the constraint that it be roughly consistent with the data.

I’ll assume a symmetric model, a normal distribution on the measurement scale. This might seem like a strong assumption but it seems reasonable enough given the data, which have an average that’s between 3 and 4. The tails are asymmetric but that’s just because the scale has 5 as an upper bound.

So I have a two-stage model:

x_latent[i] ~ normal(mu_latent, s_latent), for i=1,…,n
x_obs[i] ~ normal(x_latent[i], s_obs), for i=1,…,n
x_discrete[i] = round_and_bound(x_obs[i]), for i=1,…,n, where “round_and_bound” refers to rounding the continuous x_obs to the nearest integer between 1 and 5.

To simulate from this model, we need the hyperparameters mu_latent, s_latent, and s_obs—but with only one measurement per person, we can’t separately identify s_latent and s_obs; all we can estimate is sqrt(s_latent^2 + s_obs^2).

To separately identify s_latent and s_obs, we need to make some assumption about variation in the attractiveness measurement.

My first step is to compute the mean and sd of the observed data:

n <- round(2792*c(.02,.05,.45,.37,.11))
mean_ratings <- sum(n*(1:5))/sum(n)
sd_ratings <- sqrt(sum(n*((1:5) - mean_ratings)^2)/sum(n))

These come to 3.5 and 0.83. Assuming random zero-mean measurement error, the latent attractiveness distribution should have the same mean but a lower variance. We should subtract 1/12 to account for the rounding (as this is approximately equivalent to adding an independent uniform(-0.5,0.5) random variable, which has a variance of 1/12) and then more to account for measurement error.

But how much to subtract? For this, Attractiveness seems pretty subjective to me (also it can vary a lot over time, as we know from those high-school yearbook photos of people who later became movie stars). I gotta assume something here, so let's assume that the standard deviation of the measurements is 0.5. I'm picking this because 0.5 is a large value, but it's not as large as the sd of 0.83 in the data. Indeed, 0.5^2 + 1/12 = 0.33 and 0.83^2 = 0.69, so I'm assuming that about half of the variance in the data is coming from measurement error. If s_obs = 0.5 and the sd of the rounded data is 0.83, then s_latent should be approximately sqrt(0.83^2 - 0.5^2 - 1/12) = 0.6.

So here's a model, which I'll express in R:

round_and_bound <- function(a, lower=1, upper=5){
  pmin(upper, pmax(lower, round(a)))
}
N <- sum(n)
mu_latent <- 3.5
s_latent <- 0.6
s_obs <- 0.5
x_latent <- rnorm(N, mu_latent, s_latent)
x_obs <- rnorm(N, x_latent, s_obs)
x_discrete <- round_and_bound(x_obs)

Again, x_latent is the vector of respondents' latent attractiveness (as implicitly defined by the model), x_obs is the vector of survey interviewers' continuous perceptions of the respondents' attractiveness, and x_discrete is the recorded attractiveness on a 1-5 scale.

Let's check the results:

> print(round(table(x_discrete)/sum(n), 2))
x_discrete
   1    2    3    4    5 
0.01 0.09 0.40 0.40 0.10

OK, not perfect, but close enough. We'll get back to this model-fitting problem in a future post. For now, we'll go with the simulations we have here.

4. The prediction model

Now that we have a model of latent and recorded attractiveness, we can simulate some data.

Here's what we're gonna do.

First we'll suppose a linear relationship between latent attractiveness and some hypothetical outcome of interest. Then we'll look at what the corresponding relationship is between recorded attractiveness and that outcome. It should still be close to linear.

Next we'll do the same thing, but with a latent nonlinear model, and it turns out that this will be closer to linear. That is, adding noise to the measurement will make the E(y) vs. x relationship closer to linear.

We'll simulate the hypothetical data as follows:

x_latent and x_discrete are the 2972 attractiveness values simulated above. The latent values are all between 0 and 6.

There will be corresponding values y generated, first by specifying E(y|x_latent), then adding noise:

- In one simulation, we'll assume a linear relationship, E(y|x_latent) = x_latent/6. In the other simulation, we'll assume an S-curve, E(y|x_latent) = invlogit(2*(x_latent - 4)). We choose this curve because it is monotonic but strongly nonlinear in the range.

- My first thought was to add linear noise to the data. But I wanted to keep the range of the data between 0 and 1. Then the natural choice is binary data. But for this purpose it will be helpful to visualize with continuous data. So I'll take advantage of the constraint that E(y) is always between 0 and 1 and draw y from a beta distribution with 10 degrees of freedom, that is, y ~ beta(10*E(y|x_latent), 10*(1 - E(y|x_latent))).

Here's what we get:

The left column shows what happens when we assume a linear underlying relationship. As expected, the linear relationship is approximately preserved with the observed data.

The right column shows what happens when we assume a nonlinear underlying relationship. With the observed data, the relationship becomes much more linear.

In summary, adding measurement error attenuates the function (that is, the relation between y and x_recorded is weaker than the relation between y and x); also it brings a strongly nonlinear relation closer to linear.

This was not a surprise to me, but it's good to see it in a simulation.

# Add Health H3IR1 S35Q1 PHYSICAL ATTRACTIVENESS OF R-W3
# How physically attractive is the respondent?
# Taken from: National Longitudinal Study of Adolescent to Adult Health (Add Health), 1994-2018 [Public Use].
# in DS8

# Raw totals:  1.7% very unattractive, 4.9% unattractive, 45.3% about average, 36.7% attractive, 11.4% very attractive, N = 4877
# Section 24:  Respondent identification number (AID), Was the baby a boy or a girl (H3LB3), 1=boy (687 responses), 2 = girl (644)
# Kanazawa study:  N = 2972 ("Wave III respondents who have had at least one biological child")

library("arm")
library("cmdstanr")
linear_0 <- cmdstan_model("linear_0.stan", pedantic=TRUE)
linear_prior <- cmdstan_model("linear_prior.stan", pedantic=TRUE)
set.seed(123)

x <- seq(-2,2,1)
y <- c(50, 44, 50, 47, 56)
sexratio_data <- list(N=length(x), x=x, y=y)

display(lm(y~x))

fit_0 <- linear_0$sample(data=sexratio_data, refresh=0)
print(fit_0)
sims_0 <- fit_0$draws(format="df")

optimize_0 <- linear_0$optimize(data=sexratio_data)
print(optimize_0)
mle_0 <- optimize_0$mle()

fit_1 <-  linear_prior$sample(data=c(sexratio_data, mu_a=45.8, sigma_a=0.5, mu_b=0, sigma_b=0.2), refresh=0)
print(fit_1)
sims_1 <- fit_1$draws(format="df")

pdf("sexratio_simple.pdf", height=3, width=5)
par(mar=c(3,3,3,2), mgp=c(1.7,.5,0), tck=-.01)
plot(x+3, y, ylim=c(43, 57), xlab="Attractiveness of parent", ylab="Percentage of girl babies", bty="l", yaxt="n", main="From a survey of 3000 people",  pch=19, cex=1)
axis(2, c(45,50,55), paste(c(45,50,55), "%", sep=""))
dev.off()

pdf("sexratio_bayes.pdf", height=4, width=10)
par(mfrow=c(1,2), mar=c(3,3,3,2), mgp=c(1.7,.5,0), tck=-.01)
plot(x, y, ylim=c(43, 57), xlab="Attractiveness of parent", ylab="Percentage of girl babies", bty="l", yaxt="n", main="Least-squares estimate",  pch=19, cex=1)
axis(2, c(45,50,55), paste(c(45,50,55), "%", sep=""))
curve(mle_0["a"] + mle_0["b"]*x, col="blue", lwd=2, add=TRUE)
text(1, 49.2, paste("y = ", fround(mle_0["a"], 2), " + ", fround(mle_0["b"], 2), " x", sep=""), col="blue")

plot(x, y, ylim=c(43, 57), xlab="Attractiveness of parent", ylab="Percentage of girl babies", bty="l", yaxt="n", main="Bayes estimate with informative prior",  pch=19, cex=1)
axis(2, c(45,50,55), paste(c(45,50,55), "%", sep=""))
curve(median(sims_1$a) + median(sims_1$b)*x, col="blue", lwd=2, add=TRUE)
text(1, 45, paste("y = ", fround(median(sims_1$a), 2), " + ", fround(median(sims_1$b), 2), " x", sep=""), col="blue")
dev.off()

# Simulation

n <- round(2972*c(.02,.05,.45,.37,.11))
p <- 0.488
N <- 1000
x_sim <- array(NA, c(N, length(n)))
for (k in 1:N){
  x_sim[k,] <- rbinom(length(n), n, p)
}

pdf("sexratio_rep_1.pdf", height=5, width=8)
par(mfrow=c(4,5))
par(mar=c(2.5,3,.5,.5), mgp=c(1.5,.3,0), tck=-.01)
for (k in 1:20){
  plot((-2):2, 100*x_sim[k,]/n, ylim=c(43, 57), pch=20, xaxt="n", xlab=if (k>15) "Attractiveness of parent" else "", ylab=if (k%%5==1) "% girls" else "", yaxt="n", bty="l")
  if (k%%5==1)
    axis(2, c(45, 50, 55), c("45%", "50%", "55%"))
  else
    axis(2, c(45, 50, 55), rep("", 3))
  if (k > 15)
    axis(1, (-2):2)
  else
    axis(1, (-2):2, rep("", 5))
  abline(lm(I(100*x_sim[k,]/n) ~ I((-2):2)), col="blue")
}
dev.off()

# Attractiveness and measurement error

round_and_bound <- function(a, lower=1, upper=5){
  pmin(upper, pmax(lower, round(a)))
}

clumped <- function(a) {
  ifelse(a==1 | a==2, 1, ifelse(a==3 | a==4, 2, ifelse(a==5, 3, NA)))
}

mean_ratings <- sum(n*(1:5))/sum(n)
sd_ratings <- sqrt(sum(n*((1:5) - mean_ratings)^2)/sum(n))

N <- sum(n)
mu_latent <- 3.5
s_latent <- 0.6
s_obs <- 0.5
x_latent <- rnorm(N, mu_latent, s_latent)
x_obs <- rnorm(N, x_latent, s_obs)
x_discrete <- round_and_bound(x_obs)

print(round(table(x_discrete)/sum(n), 2))

# Assume a relationship

linear_curve <- function(x) {
  x/6
}
s_curve <- function(x) {
  invlogit(2*(x-4))
}
response_curve <- list(linear=linear_curve, nonlinear=s_curve)

pdf("latent_error.pdf", width=6, height=4.5)
par(mfcol=c(2,2), mar=c(3,3,2,2), mgp=c(1.2,.2,0), tck=-.01)

for (k in 1:2) {
  p <- response_curve[[k]](x_latent)
  y <- rbeta(N, 10*p, 10*(1-p))
  plot(x_latent, y, ylim=c(0, 1), xlab="Latent beauty", ylab="Outcome", yaxt="n", pch=20, cex=.1, main=paste("Latent", names(response_curve)[k], "model"), cex.axis=.9, cex.lab=.9, cex.main=1)
  axis(2, c(0, 0.5, 1))
  curve(response_curve[[k]](x), col="white", lwd=3, add=TRUE)
  curve(response_curve[[k]](x), col="blue", lwd=2, add=TRUE)
  plot(x_discrete + runif(N, -.1, .1), y, ylim=c(0, 1), xlab="Recorded beauty (jittered)", ylab="Outcome", yaxt="n", pch=20, cex=.1, main=paste("Expected values given recorded beauty"), cex.axis=.9, cex.lab=.9, cex.main=.8)
  axis(2, c(0, 0.5, 1))
  Ey <- rep(NA, 5)
  sy <- rep(NA, 5)
  for (i in 1:5) {
    in_bin <- x_discrete==i
    Ey[i] <- mean(response_curve[[k]](x_latent[in_bin]))
    sy[i] <- sd(response_curve[[k]](x_latent[in_bin]))
  }
  lines(1:5, Ey, col="white", lwd=3)
  lines(1:5, Ey, col="red", lwd=2)
}
dev.off()