Logarithms and heteroskedasticity¶

Here’s a question from the Reddit statistics forum.

Is it correct to use logarithmic transformation in order to mitigate heteroskedasticity?

For my studies I gathered data on certain preferences across a group of people. I am trying to figure out if I can pinpoint preferences to factors such as gender in this case.

I used mixed ANOVA analysis with good success however one of my hypothesis came up with heteroskedasticity when doing Levene’s test. [I’ve been] breaking my head all day on how to solve this. I’ve now used logarithmic transformation to all 3 test results and run another Levene’s. When using the media value the test now results [in] homoskedasticity, however interaction is no longer significant?

Is this the correct way to deal with this problem or is there something I am missing? Thanks in advance to everyone taking their time to help.

Although the question is about ANOVA, I’m going to reframe it in terms of regression, for two reasons:

• Discussion of heteroskedasticity is clearer in the context of regression.

• For many problems, a regression model is better than ANOVA anyway.

What is heteroskedasticity?¶

Linear regression is based on a model of the data-generating process where the independent variable y is the sum of

• A linear function of x with an unknown slope and intercept, and

• Random values drawn from a Gaussian distribution with mean 0 and unknown standard deviation, sigma, that does not depend on x.

If, contrary to the second assumption, sigma depends on x, the data-generating process is heteroskedastic. Some amount of heteroskedasticity is common in real data, but most of the time it’s not a problem, because:

• Heteroskedasticity doesn’t bias the estimated parameters of the model, only the standard errors,

• Even very strong heteroskedasticity doesn’t affect the standard errors by much, and

• For practical purposes we don’t need standard errors to be particularly precise.

To demonstrate, let’s generate some data. First, we’ll draw xs from a normal distribution.

To generate heteroskedastic data, I’ll use interpolate to construct a function where sigma depends on x.

To generate strong heteroskedasticity, I’ll vary sigma over a wide range.

Here’s what sigma looks like as a function of x.

Now we can generate ys with variable values of sigma.

If we make a scatter plot of the data, we see a cone shape that indicates heteroskedasticity.

Now let’s fit a model to the data.

Here’s what the fitted line looks like.

If we plot the absolute values of the residuals, we can see the heteroskedasticity more clearly.

Testing for heteroskedasticity¶

OP mentions using the Levene test for heteroskedasticity, which is used to test whether sigma is different between groups. For continuous values of x and y, we can use the Breusch-Pagan Lagrange Multiplier test:

Or White’s Lagrange Multiplier test:

Both tests produce small p-values, which means that if we generate a dataset by a homoskedastic process, there is almost no chance it would have as much heteroskedasticity as the dataset we generated.

If you have never heard of either of these tests, don’t panic — neither had I under I looked them up for this example. And don’t worry about remembering them, because you should never use them again. Like testing for normality, testing for heteroskedasticity is never useful.

Why? Because in almost any real dataset, you will find some heteroskedasticity. So if you test for it, there are only two possible results:

• If the heteroskedasticity is small and you don’t have much data, you will fail to reject the null hypothesis.

• If the heteroskedasticity is large or you have a lot of data, you will reject the null hypothesis.

Either way, you learn nothing — and in particular, you don’t learn the answer to the question you actually care about, which is whether the heteroskedasticity is so large that the effect on the standard errors is large enough that you should care.

And the answer to that question is almost always no.

Should we care?¶

The dataset we generated has very large heteroskedasticity. Let’s see how much effect that has on the results. Here are the standard errors from simple linear regression:

Now, there are several ways to generate standard errors that are robust in the presence of heteroskedasticity. One is the Huber-White estimator, which we can compute like this:

Another is to use Huber regression.

Another is to use quantile regression.

And one more option is a wild bootstrap, which resamples the residuals by multiplying them by a random sequence of 1 and -1. This way of resampling preserves heteroskedasticity, because it only changes the sign of the residuals, not the magnitude, and it maintains the relationship between those magnitudes and x.

We can use wild_bootstrap to generate a sample from the sampling distributions of the intercept and slope.

The standard deviation of the sampling distributions is the standard error.

Now let’s put all of the result in a table.

The standard errors we get from different methods are notably different, but the differences probably don’t matter.

First, I am skeptical of the results from Huber regression. With this kind of heteroskedasticity, the standard errors should be larger than what we get from OLS. I’m not sure what’s the problem is, and I haven’t bothered to find out, because I don’t think Huber regression is necessary in the first place.

The results from bootstrapping and the Huber-White estimator are the most reliable — which suggests that the standard errors from quantile regression are too big.

In my opinion, we don’t need esoteric methods to deal with heteroskedasticity. If heteroskedasticity is extreme, consider using wild bootstrap. Otherwise, just use ordinary least squares.

Now let’s address OP’s headline question, “Is it correct to use logarithmic transformation in order to mitigate heteroskedasticity?”

Log transform help?¶

In some cases, a log transform can reduce or eliminate heteroskedasticity. However, there are several reasons this is not a good idea in general:

• As we’ve seen, heteroskedasticity is not a big problem, so it usually doesn’t require any mitigation.

• Taking a log transform of one or more variables in a regression model changes the meaning of the model — it hypothesizes a relationship between the variables that might not make sense in context.

• Anyway, taking a log transform doesn’t always help.

To demonstrate the last point, let’s see what happens if we apply a log transform to the dependent variable:

Here’s what the scatter plot looks like after the transform.

Here’s what we get if we fit a model to the data.

And here’s the fitted line.

If we plot the absolute values of the residuals, we can see that the log transform did not entirely remove the heteroskedasticity.

Which we can confirm by running the tests again (which we should never do).

The p-values are bigger, which suggests that the log transform mitigated the heteroskedasticity a little. But if the goal was to eliminate heteroskedasticity, the log transform didn’t do it.

Discussion¶

To summarize:

• Heteroskedasticity is common in real datasets — if you test for it, you will often find it, provided you have enough data.

• Either way, testing does not answer the question you really care about, which is whether the heteroskedasticity is extreme enough to be a problem.

• Plain old linear regression is robust to heteroskedasticity, so unless it is really extreme, it is probably not a problem.

• Even in the worst case, heteroskedasticity does not bias the estimated parameters — it only affects the standard errors — and we don’t need standard errors to be particularly precise anyway.

• Although a log transform can sometimes mitigate heteroskedasticity, it doesn’t always succeed, and even if it does, it’s usually not necessary.

• A log transform changes the meaning of the regression model in ways that might not make sense in context.

So, use a log transform if it makes sense in context, not to mitigate a problem that’s not much of a problem in the first place.

Data Q&A: Answering the real questions with Python

Copyright 2024 Allen B. Downey

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

Combining Risks¶

Here’s a question from the Reddit statistics forum.

Bit of a weird one but I’m hoping you’re the community to help. I work in children’s residential care and I’m trying to find a way of better matching potential young people together.

The way we calculate individual risk for a child is

risk = likelihood + impact (R=L+I), so

L4 + I5 = R9

That works well for individuals but I need to work out a good way of calculating a combined risk to place children [in] the home together. I’m currently using the [arithmetic] average but I don’t feel that it works properly as the average is always lower then the highest risk.

I’ll use a fairly light risk as an example, running away from the home. (We call this MFC missing from care) It’s fairly common that one of the kids will run away from the home at some point or another either out of boredom or frustration. If young person A has a risk of 9 and young person B has a risk of 12 the the average risk of MFC in the home would be 10.5

HOWEVER more often then not having two young people that go MFC will often result in more episodes as they will run off together, so having a lower risk rating doesn’t really make sense. Adding the two together to 21 doesn’t really work either though as the likelihood is the thing that increases not necessarily the impact.

I’m a lot better at chasing after run away kids then I am mathematics so please help 😂.

Here’s one way to think about this question: based on background knowledge and experience, OP has qualitative ideas about what happens when we put children at different risks together, and he is looking for a statistical summary that is consistent with these ideas.

The arithmetic mean probably makes sense as a starting point, but it clashes with the observation that if you put two children together who are high risk, they interact in ways that increase the risk. For example, if we put together children with risks 9 and 12, the average is 10.5, and OP’s domain knowledge says that’s too low — it should be more than 12.

At the other extreme, I’ll guess that putting together two low risk children might be beneficial to both — so the combined risk might be lower than either.

And that implies that there is a neutral point somewhere in the middle, where the combined risk is equal to the individual risks.

To construct a summary statistic like that, I suggest a weighted sum of the arithmetic and geometric means. That might sound strange, but I’ll show that it has the properties we want. And it might not be as strange as it sounds — there’s a reason it might make sense.

Weighted sum of means¶

The following function computes the arithmetic mean of a sequence of values, which is the sum divided by n.

The following function computes the geometric mean of a sequence, which is the product raised to the power 1/n.

And the following function computes the weighted sum of the arithmetic and geometric means. The constant k determines how much weight we give the geometric mean.

The value 4 determines the neutral point. So if we put together two people with risk 4, the combined average is 4.

Above the neutral point, there is a penalty if we put together two children with higher risks.

In that case, the combined risk is higher than the individual risks. Below the neutral point, there is a bonus if we put together two children with low risks.

In that case, the combined risk is less than the individual risks.

If we combine low and high risks, the discrepancy brings the average down a little.

In the example OP presented, where we put together two people with high risk, the penalty is substantial.

If that penalty seems too high, we can adjust the weight, k, and the neutral point accordingly.

This behavior extends to more than two people. If everyone is neutral, the result is neutral.

If you add one person with higher risk, there’s a moderate penalty, compared to the arithmetic mean.

With two higher risk people, the penalty is bigger.

And with three it is even bigger.

Does this make any sense?¶

The idea behind this suggestion is logistic regression with an interaction term. Let me explain where that comes from. OP explained:

The way we calculate individual risk for a child is

risk = likelihood + impact (R=L+I), so

L4 + I5 = R9

At first I thought it was strange to add a likelihood and an impact score, Thinking about expected value, I thought it should be the product of a probability and a cost. But if both are on a log scale, adding these logarithms is like multiplying probability by impact on a linear scale, so that makes more sense.

And if the scores are consistent with something like a log-odds scale, we can see a connection with logistic regression. If r1 and r2 are risk scores, we can imagine a regression equation that looks like this, where p is the probability of an outcome like “missing from care”:

logit(p) = a r1 + b r2 + c r1 r2

In this equation, logit(p) is the combined risk score, a, b, and c are unknown parameters, and the product r1 r2 is an interaction term that captures the tendency of high risks to amplify each other.

With enough data, we could estimate the unknown parameters. Without data, the best we can do is chose values that make the results consistent with expectations.

Since r1 and r2 are interchangeable, they have to have the same parameter. And since the whole risk scale has an unspecified zero point, we can set it a and b to 1/2. Which means there is only one parameter left, the weight of the interaction term.

logit(p) = (r_1 + r2) / 2 + k r1 r2

Now we can see that the first term is the arithmetic mean and the second term is close to the geometric mean, but without the square root.

So the function I suggested — the weighted sum of arithmetic and weighted means — is not identical to the logistic model, but it is motivated by it.

With this rationale in mind, we might consider a revision: rather than add the likelihood and impact scores, and then compute the weighted sum of means, it might make more sense to separate likelihood and impact, compute the weighted sum of the means of the likelihoods, and then add back the impact.

Computing by hand¶

In case the Python code makes it hard to see what’s going on, let’s work an example by hand. Suppose r1 is 9 and r2 is 12.

Here’s the arithmetic mean.

Here’s the geometric mean.

And here’s how we combine them.

Discussion¶

This question got my attention because OP is working on a challenging and important problem — and they provided useful context. It’s an intriguing idea to define something that is intuitively like an average, but is not always bounded between the minimum and maximum of the data.

If we think strictly about generalized means, that’s not possible. But if we think in terms of logarithms, regression, and interaction terms, we find a way.

Data Q&A: Answering the real questions with Python

Copyright 2024 Allen B. Downey

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

Estimation with Small Samples¶

Here’s a question from the Reddit statistics forum.

Hey, so imagine I only have 6 samples from a value that has a normal distribution. Can I estimate the range of likely distributions from those 6?

Let’s be more specific. I’m considering the accuracy of a blood testing device. I took 6 samples of my blood at the same time from the same vein and gave them to the machine. The results are not all the same (as expected), indicating the device’s inherent level of imprecision.

So, I’m wondering if there’s a way to estimate the range of possibilities of what I would see if I could give 100 or 1000 samples?

I’m comfortable assuming a normal distribution around the “true” value. Is there any stats method to guesstimate the range of likely values for sigma? Or would I just need to drain my blood dry to get 1000 samples to figure that out?

Fyi, not a statistician.

Because the sample size is so small, this question cries out for a Bayesian approach. Why?

• Bayesian methods do a good job of taking advantage of background information and extracting as much information as possible from the data,

• With a small sample size, the uncertainty of the result is large, so it is important to quantify it, and

• The motivating question is explicitly about making probabilistic predictions, which is what Bayesian methods do and classical methods don’t.

As an example, OP suggested looking at blood potassium (K+) levels, with the following data:

4.0, 3.9, 4.0, 4.0, 4.7, 4.2 (mmol/L)

OP also explained that the blood samples were collected “continuously from a single puncture … within seconds, not minutes,” so we don’t expect the true level to change much between samples. In that case, the variation in the measurements would largely reflect the precision of the testing device.

With that assumption, let’s see what we can do.

Bayesian estimation¶

Here’s the data.

Now we need prior distributions for mu and sigma. For the prior distribution of mu we can use the distribution in the general population. The “normal range” for potassium — which is usually the 5th and 95th percentiles of the population — is 3.5 to 5.4 mmol/L. So we can construct a normal distribution with these percentiles.

Now I’ll make a Pmf object that represents this distribution.

For the prior distribution of sigma, we could use some background information about the device. For example, based on similar devices, what values of sigma would be expected? I’ll use a gamma distribution to construct a prior PMF, but it could be anything.

This distribution represents the background knowledge that we expect the standard deviation to be less than 2, and we’d be surprised by anything much higher than that.

I’ll use this function to make a Pandas DataFrame to represent the joint prior.

I’ll use the following function to make a contour plot of the prior.

The update¶

To use the data to update the prior, we have to compute the likelihood of the data for each possible pair of mu and sigma. Can can do that by creating a 3-D mesh with the possible values of mu and sigma, and the observed values of the data.

Now we can evaluate the normal distribution for each data point and each pair of mu and sigma.

The likelihood of each pair is the product of the densities for the data points.

The unnormalized posterior is the product of the prior and the likelihood.

We can normalize it like this.

And here’s what the result looks like.

Here’s the posterior distribution of sigma compared to the prior.

The posterior mean of sigma is about 0.4, somewhat higher than the standard deviation of the data.

Predictions¶

The last part of OP’s question is “So, I’m wondering if there’s a way to estimate the range of possibilities of what I would see if I could give 100 or 1000 samples?”

We can simulate a future experiment with larger sample size by drawing random pairs from the joint posterior distribution and generating simulated measurements. That’s easier to do if we stack the posterior PMF.

Now we can use NumPy’s choice function to choose a random pair of mu and sigma. For each random pair, we can generate 1000 simulated measurements from a normal distribution, and plot the distribution of the results.

Each line shows the distribution of a possible sample of 1000 measurements. You can see that they vary in both location and spread, due to the uncertainty represented by the posterior distribution of mu and sigma.

Discussion¶

The normal distribution is probably a good enough model for this data, but for some pairs of mu and sigma in the prior, the normal distribution extends into negative values with non-negligible probability — and for measurements like these, negative values are nonsensical.

An alternative would be to run the whole calculation with the logarithms of the measurements. In effect, this would model the distribution of the data as lognormal, which is generally a good choice for measurements like these.

In this example, the sample size is small, so the choice of the prior is important. I suspect there is more background information we could use to make a better choice for the prior distribution of sigma.

This example demonstrates the use of grid methods for Bayesian statistics. For many common statistical questions, grid methods like this are simple and fast to compute, especially with libraries like NumPy and SciPy. And they make it possible to answer many useful probabilistic questions.

Data Q&A: Answering the real questions with Python

Copyright 2024 Allen B. Downey

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

Sample Size Selection¶

Here’s a question from the Reddit statistics forum.

Hi Redditors, I am a civil engineer trying to solve a statistical problem for a current project I have. I have a pavement parking lot 125,000 sf in size. I performed nondestructive testing to render an opinion about the areas experiencing internal delimitation not observable from the surface. Based on preliminary testing, it was determined that 9% of the area is bad, and 11% of the total area I am unsure about (nonconclusive results if bad or good), and 80% of the area is good. I need to verify all areas using destructive testing, I will take out slabs 2 sf in size. my question is how many samples do I need to take from each area to confirm the results with 95% confidence interval?

There are elements of this question that are not clear, and OP did not respond to follow-up questions. But the question is generally about sample size selection, so let’s talk about that.

If the parking lot is 125,000 sf and each sample is 2 sf, we can imagine dividing the total area into 62,500 test patches. Of those, some unknown proportion are good and the rest are bad.

In reality, there is probably some spatial correlation — if a patch is bad, the nearby patches are more likely to be bad. But if we choose a sample of patches entirely at random, we can assume that they are independent. In that case, we can estimate the proportion of patches that are good and quantify the precision of that estimate by computing a confidence interval.

Then we can choose a sample size that meets some requirement. For example, we might want the 95% confidence interval needs to be smaller than a given threshold, or we might want to bound the probability that the proportion falls below some threshold.

But let’s start by estimating proportions and computing confidence intervals.

The beta-binomial model¶

Based on preliminary testing, it is likely that the proportion of good patches is between 80% and 90%. We can take advantage of that information by using a beta distribution as a prior and updating it with the data. Here’s a prior distribution that seems like a reasonable choice, given the background information.

The prior mean is at the low end of the likely range, so the results will be a little conservative. Here’s what the prior looks like.

This prior leaves open the possibility of values below 80% and greater than 90%, but it assigns them lower probabilities.

Now let’s generate a hypothetical dataset to see what the update looks like. Suppose the actual percentage of good patches is 90%, and we sample n=10 of them.

And suppose that, in line with expectations, 9 out of 10 tests are good.

Under the beta-binomial model, computing the posterior is easy.

Here’s how we run the update.

The posterior mean is 85%, which is half way between the prior mean and the proportion observed in the data.

Here’s what the posterior distribution looks like, compared to the prior.

Given the posterior distribution, we can use ppf, which computes the inverse CDF, to compute a confidence interval.

Here’s the result for this example.

With a sample size of only 10, the confidence interval is still quite wide — that is, the estimate of the proportion is not precise.

More data?¶

Now let’s run the same analysis with a sample size of n=100.

With a larger sample size, the posterior mean is closer to the proportion observed in the data. And the posterior distribution is narrower, which indicates greater precision.

The confidence interval is much smaller.

If we need more precision than that, we can increase the sample size more. If we don’t need that much precision, we can decrease it.

With some math, we could compute the sample size algorithmically, but a simple alternative is to run this analysis with different sample sizes until we get the results we need.

But what about that prior?¶

Some people don’t like using Bayesian methods because they think it is more objective to ignore perfectly good background information, even in cases like this where they come from preliminary testing that is clearly applicable.

To satisfy them, we can run the analysis again with a uniform prior, which is not actually more objective, but it seems to make people happy.

The mean of the uniform prior is 50%, so it is more pessimistic. Here’s the update with n=10.

Now let’s compare the posterior distributions with the informative prior and the uniform prior.

With the informative prior, the posterior distribution is a little narrower — an estimate that uses background information is more precise.

Let’s do the same thing with n=100.

With a larger sample size, the choice of the prior has less effect — the posterior distributions are almost the same.

Discussion¶

Sample size analysis is a good thing to do when you are designing experiments, because it requires you to

• Make a model of the data-generating process,

• Generate hypothetical data, and

• Specify ahead of time what analysis you plan to do.

It also gives you a preview of what the results might look like, so you can think about the requirements. If you do these things before running an experiment, you are likely to clarify your thinking, communicate better, and improve the data collection process and analysis.

Sample size analysis can also help you choose a sample size, but most of the time that’s determined by practical considerations, anyway. I mean, how many holes do you think they’ll let you put in that parking lot?

Data Q&A: Answering the real questions with Python

Copyright 2024 Allen B. Downey

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

Likert scale analysis¶

Here’s a question from the Reddit statistics forum.

I have collected data regarding how individuals feel about a particular program. They reported their feelings on a scale of 1-5, with 1 being Strongly Disagree, 2 being Disagree, 3 being Neutral, 4 being Agree, and 5 being Strongly Agree.

I am looking to analyze the data for averages responses, but I see that a basic mean will not do the trick. I am looking for very simple statistical analysis on the data. Could someone help out regarding what I would do?

It sounds like OP has heard the advice that you should not compute the mean of values on a Likert scale. The Likert scale is ordinal, which means that the values are ordered, but it is not an interval scale, because the distances between successive points are not necessarily equal.

For example, if we imagine that “Neutral” maps to 0 on a number line and “Agree” maps to 1, it’s not clear where we should place “Strongly agree”. And an appropriate mapping might not be symmetric — for example, maybe the people who choose “Strongly agree” are content, but the people who choose “Strongly disagree” are angry. In an arithmetic mean, they would cancel each other out — but that might erase an important distinction.

Nevertheless, I think an absolute prohibition on computing means is too strong. I’ll show some examples where I think it’s a reasonable thing to do — but I’ll also suggest alternatives that might be better.

Political views¶

As an example, I’ll use data from the General Social Survey (GSS), which I have resampled to correct for stratified sampling.

The variable we’ll start with is polviews, which contains responses to this question:

We hear a lot of talk these days about liberals and conservatives. I’m going to show you a seven-point scale on which the political views that people might hold are arranged from extremely liberal–point 1–to extremely conservative–point 7. Where would you place yourself on this scale?

This is not a Likert scale, specifically, but it is ordinal. Here is the distribution of responses:

Here’s the mapping from numerical values to the options that were shown to respondents.

As always, it’s a good idea to visualize the distribution before we compute any summary statistics. I’ll use Pmf from empiricaldist to compute the PMF of the values.

Here’s what it looks like.

The modal value is “Moderate” and the distribution is roughly symmetric, so I think it’s reasonable to compute a mean. For example, suppose we want to know how self-reported political alignment has changed over time. We can compute the mean in each year like this:

And plot it like this.

I used LOWESS to plot a local regression line, which makes long-term trends easier to see. It looks like the center of mass trended toward conservative from the 1970s into the 1990s and trended toward liberal since then.

When you compute any summary statistic, you lose information. To see what we might be missing, we can use a normalized cross-tabulation to compute the distribution of responses in each year.

And we can use a heat map to visualize the results.

Based on the heat map, it seems like the general shape of the distribution has not changed much — so the mean is probably a good way to make comparisons over time.

However, it is hard to interpret the mean in absolute terms. For example, in the most recent data, the mean is about 4.1.

Since 4.0 maps to “Moderate”, we can say that the center of mass is slightly on the conservative side of moderate, but it’s hard to say what a difference of 0.1 means on this scale.

As an alternative, we could add up the percentage who identify as conservative or liberal, with or without an adverb.

And plot those percentages over time.

Or we could plot the difference in percentage points.

This figure shows the same trends we saw by plotting the mean, but the y-axis is more interpretable — for example, we could report that, at the peak of the Reagan era, conservatives outnumbered liberals by 10-15 percentage points.

Standard deviation¶

Suppose we are interested in polarization, so we want to see if the spread of the distribution has changed over time. Would it be OK to compute the standard deviation of the responses? As with the mean, my answer is yes and no.

First, let’s see what it looks like.

The standard deviation is easy to compute, and it makes it easy to see the long-term trend. If we interpret the spread of the distribution as a measure of polarization, it looks like it has increased in the last 30 years.

But it is not easy to interpret this result in context. If it increased from about 1.35 to 1.5, is that a lot? It’s hard to say.

As an alternative, let’s compute the mean absolute deviation (MAD), which we can think of like this: if we choose two people at random, how much will they differ on this scale, on average?

A quick way to estimate MAD is to draw two samples from the responses and compute the mean pairwise distance.

The result is about 1.5 points, which is bigger than the distance from moderate to slightly conservative, and smaller than the distance from moderate to conservative.

Rather than sampling, we can compute MAD deterministically by forming the joint distribution of response pairs and computing the expected value of the distances. For this computation, it is convenient to use NumPy functions for outer product and outer difference.

Again, the result is about 1.5 points.

Now we can see how this value has changed over time.

Here’s the result, along with the standard deviation.

The two figures tell the same story — polarization is increasing. But the MAD is easier to interpret. In the 1970s, if you chose two people at random, they would differ by less than 1.5 points on average. Now the difference would be almost 1.7 points. Considering that the difference between a moderate and a conservative is 2 points, it seems like we should still be able to get along.

I think MAD is more interpretable than standard deviation, but it is based on the same assumption that the points on the scale are equally spaced. In most cases, that’s not an assumption we can easily check, but in this example, maybe we can.

For another project, I selected 15 questions in the GSS where conservatives and liberals are most likely to disagree, and used them to estimate the number of conservative responses from each respondent. The following figure shows the average number of conservative responses to the 15 questions for each point on the self-reported scale.

The result is close to a straight line, which suggests that the assumption of equal spacing is not bad in this case.

When is the mean bad?¶

In the examples so far, computing the mean and standard deviation of a scale variable is not necessarily the best choice, but it could be a reasonable choice. Now we’ll see an example where it is probably a bad choice.

The variable homosex contains responses to this question:

What about sexual relations between two adults of the same sex–do you think it is always wrong, almost always wrong, wrong only sometimes, or not wrong at all?

If the wording of the question seems loaded, remember that many of the core questions in the GSS were written in the 1970s. Here is the encoding of the responses.

And here are the value counts.

Before we do anything else, let’s look at the distribution of responses.

There are several reasons it’s a bad idea to summarize this distribution by computing the mean. First, one of the responses is not ordered. If we include “Other” in the mean, the result is meaningless.

If we exclude “Other”, the remaining responses are ordered, but arguably not evenly spaced on a spectrum of opinion.

If we compute a mean and report that the average response is somewhere between “Sometimes wrong” and “Almost always wrong”, that is not an effective summary of the distribution.

The distribution of results is strongly bimodal — most people are either accepting of homosexuality or not. And that suggests a better way to summarize the distribution: we can simply report the fraction of respondents who choose one extreme or the other.

I’ll start by creating a binary variable that is 1 for respondents who chose “Not wrong at all”, 0 for the other responses, and NaN for people who were not asked the question, did not respond, or chose “Other”.

The mean of this variable is the fraction of respondents who chose “Not wrong at all”.

Now we can see how this fraction has changed over time.

The percentage of people who accept homosexuality was almost unchanged in the 1970s and 1980s, and began to increase quickly around 1990. For a discussion of this trend, and similar trends related to racism and sexism, you might be interested in Chapter 11 of Probably Overthinking It.

Discussion¶

Computing the mean of an ordinal variable can be a quick way to make comparisons between groups or show trends over time. The computation implicitly assumes that the points on the scale are equally spaced, which is not true in general, but in many cases it is close enough.

However, even when the mean (or standard deviation) is a reasonable choice, there is often an alternative that is easier to interpret in context.

It’s always good to look at the distribution before choosing a summary statistic. If it’s bimodal, the mean is probably not the best choice.

Data Q&A: Answering the real questions with Python

Copyright 2024 Allen B. Downey

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

Appendix: Quantifying discord¶

We’ve seen that mean absolute deviation (MAD) can quantify the average level of disagreement between two people in a population. As a generalization, it can also quantify the discord between two groups.

I might come back and expand this example later.

Testing percentiles¶

Here’s a question from the Reddit statistics forum.

I have two different samples (about 100 observations per sample) drawn from the same population (or that’s what I hypothesize; the populations may in fact be different). The samples and population are approximately normal in distribution.

I want to estimate the 85th percentile value for both samples, and then see if there is a statistically significant difference between these two values. I cannot use a normal z- or t-test for this, can I? It’s my current understanding that those tests would only work if I were comparing the means of the samples.

As an extension of this, say I wanted to compare one of these 85th percentile values to a fixed value; again, if I was looking at the mean, I would just construct a confidence interval and see if the fixed value fell within it…but the percentile stuff is throwing me for a loop.

This is […] related to a research project I’m working on (in my job).

There are two questions here. The first is about testing a difference in percentiles between two groups. The second is about the difference between a percentile from an observed sample and an expected value.

We’ll answer the first question with a permutation test, and we’ll answer the second in two ways: bootstrap resampling and a Gaussian model.

Data¶

Since OP didn’t provide a dataset, we have to generate one. I’ll draw two samples from Gaussian distributions with the same standard deviation and different means.

Here’s what the distributions of the groups look like.

If we compute the 85th percentile in both groups, we see a difference, as expected.

Now let’s see if a difference of that size would be likely if the two samples were actually drawn from the same distribution.

Testing the Difference¶

When we test a difference between two groups, the usual model of the null hypothesis is that the groups are actually identical. If that’s true, the two samples came from the same distribution, so we can combine them into a single sample.

Now we can simulate the null hypothesis by permutation — that is, by shuffling the pooled data and splitting it into two groups with the same sizes as the originals. The following function generates two samples under this assumption, and returns the difference in their 85th percentiles.

If we call it many times, the result is a sample from the distribution of differences under the null hypothesis.

Here’s what it looks like, with a vertical line at the observed difference. I’m plotting it with cut=0 so the estimated density doesn’t extend past the minimum and maximum of the data.

The distribution is multimodal, which is the result of selecting a moderately high percentile from a moderately small dataset — the diversity of the results is limited. However, in this example we are interested in the tails of the distribution, so multimodality is not a problem.

To estimate a one-sided p-value, we can compute the fraction of the sample that exceeds the actual difference.

Or, for a two-sided p-value, we can compute the fraction of the sample that exceeds the actual difference in absolute value.

In this example, the result of the one-sided test would be considered significant at the 5% significance level, but the two-sided test would not. So which is it?

I think it’s not worth worrying about. My interpretation of the results is the same either way: they are inconclusive. Under the null hypothesis, a difference as big as the one we saw would be unlikely, but we can’t rule out the possibility that the groups are identical — or nearly so — and the apparent difference is due to random variation.

Testing a fixed value¶

Now let’s turn to the second question. Suppose we have reason to think that the actual value of the 85th percentile is 12.3, and we would like to know whether the data contradict this hypothesis.

We’ll test group1 first. Here’s the 85th percentile of group1 and its difference from the expected value.

Let’s see if a difference of this magnitude is likely to happen under the null hypothesis. One way to model the null hypothesis is to create a dataset that is similar to the observed data, but where the 85th percentile is exactly as expected. We can do that by shifting the observed data by the observed difference.

The 85th percentile of the shifted data is the expected value, exactly.

Now, to generate samples under the null hypothesis, we can use the following function, which takes a sample, shifts it to have the expected value of the 85th percentile, generates a bootstrap resample of the shifted values, and returns the difference between the 85th percentile of the sample and the expected value.

If we call this function many times, we get a sample of the differences we expect under the null hypothesis.

The following function shows the distribution of the sample with a vertical line at the observed value.

Without computing a p-value, we can see that a difference as big as actual_diff1 is entirely plausible under the null hypothesis. We can confirm that by computing a one-sided p-value.

So the observed difference in the first group is not statistically significant. Now let’s do the same thing for the second group.