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.
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, we need to make some assumption about variation in the attractiveness measurement. 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.
But is 0.5 a reasonable value? To answer that question, it helps to have data.
The data at hand come from wave 3 of the Add Health study. I already downloaded those data awhile ago in preparing my earlier analysis, but I've had difficulty working with the numbers. I can't remember why, but there's no single dataset and no single codebook. It's kind of a mess. I spent a few minutes going through to see if there were any repeated measurements of attractiveness there, but I wasn't sure where to look.
So I googled *reliability of add health attractiveness measure*. I was hoping to find the names of the relevant variables in the survey, and that didn't come up in the search. What did appear was a bunch of research papers on the topic. I looked at a few of them and then came across this one: "Ugly criminals," by Naci Mocan and Erdal Tekin. I'm not so interested in the topic of the paper---I'll leave it to others to decide if its conclusions are reasonably supported by its data---; what's relevant here is that they discuss the Add Health attractiveness measure. Which is indeed how it came up in my search.
They looked at data from respondents whose attractiveness had been recorded in two waves. That's exactly what we need for a reliability analysis:
There's only one mild annoyance here: they pooled categories 1 and 2, and they pooled categories 3 and 4. So we don't have the full 5 x 5 table, only a 3 x 3 collapsing of it.
But we can deal with that. We'll just simulate and see what happens.
First we write a couple of simple functions:
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)))
}
We set up the parameters for the latent data:
mu_obs <- c(3.87, 3.69) s_latent <- 0.53 s_obs <- c(0.82, 0.52)
Then we simulate two independent measurements on each of 14654 people---that's the number of cases in the above table:
y_latent <- rnorm(N, 0, s_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
y_obs[,j] <- round_and_bound(rnorm(N, mu_obs[j] + y_latent, s_obs[j]))
}
And then we look at the two margins and the table showing the joint distribution of the two measurements:
print(table(clumped(observed_1))) print(table(clumped(observed_2))) print(table(clumped(observed_1), clumped(observed_2)))
I also looked at the unclumped data but I won't bother showing you that here.
Here are the results:
> print(table(clumped(observed_1)))
1 2 3
668 11186 2800
> print(table(clumped(observed_2)))
1 2 3
723 11221 2710
> print(table(clumped(observed_1), clumped(observed_2)))
1 2 3
1 227 437 4
2 490 9292 1404
3 6 1492 1302
These don't quite match the table above. The simulations have not enough "unattractive" people and too many "very attractive." Also they two waves in the simulations are essentially the same (by construction, as we simulated them from the same distribution), but the real data show a few more unattractive people and a lot more very attractive people in wave 1.
There would be a way to set up my simulation to alter its parameters and get the closest possible match to the observed two-way table---I guess it wouldn't be hard to program this in Stan, actually! But here I took the lazy way out and just tinkered with the numbers to get things to work. It took about five tries of playing with the inputs, and eventually I converged on this:
sd_measurement <- 0.5 mu_latent <- 3.8 sd_latent <- 0.67 latent_beauty <- rnorm(n_sim, mu_latent, sd_latent) observed_1 <- round_and_bound(rnorm(n_sim, latent_beauty, sd_measurement)) observed_2 <- round_and_bound(rnorm(n_sim, latent_beauty - 0.1, sd_measurement*0.6))
As you can see, I hacked the simulations to give the wave 2 measurement a lower mean and lower variance.
Here's the result:
> print(table(clumped(observed_1)))
1 2 3
845 10839 2970
> print(table(clumped(observed_2)))
1 2 3
731 11852 2071
> print(table(clumped(observed_1), clumped(observed_2)))
1 2 3
1 330 514 1
2 400 9727 712
3 1 1611 1358
OK, now I have the margins pretty good, but the joint table doesn't have enough people whose measurements shifted between the lowest and highest categories.
I played around some more with the parameters of the model but I couldn't get the counts to line up with the data.
So I bit the bullet and did what I should've done from the start, which was to write a Stan program to find a best fit. That was its own struggle---as is often the case, I learn a lot from these little examples---so I'll devote a separate post to it.
Here's what I ended up using:
mu_obs <- c(3.87, 3.69)
s_latent <- 0.53
s_obs <- c(0.82, 0.52)
y_latent <- rnorm(N, 0, s_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
y_obs[,j] <- round_and_bound(rnorm(N, mu_obs[j] + y_latent, s_obs[j]))
}
And here's the result, which is still not perfect (as will be discussed in my followup post) but is enough for my goals here:
> for (j in 1:2) print(table(clumped(y_obs[,j])))
1 2 3
1225 9684 3745
1 2 3
815 11857 1982
> print(table(clumped(y_obs[,1]), clumped(y_obs[,2])))
1 2 3
1 208 990 27
2 563 8141 980
3 44 2726 975
Which is far from a perfect match to the actual table of responses from the two waves:
> print(N_table)
[,1] [,2] [,3]
[1,] 208 640 40
[2,] 501 9879 464
[3,] 51 1439 1432
But it's close enough, as we're just trying to make a general point about attenuation of effects under measurement error.
