# Mean substitution - skewness and kurtosis

I am doing the review questions from Thompson: foundations of behavioral statistics, chapter 4, question 2. I cannot seem to conceptualize the correct answer. Would appreciate any help as I need to understand this review for the final exam at the end of term :)

"For a given dataset, how will using mean substitution for missing values impact the coefficient of skewness for data that were initially skewed? How will using mean substitution for missing values impact the coefficient of kurtosis? How does mean substitution differentially impact post-substitution means, SD, skewness, and kurtosis?"

Thank you!

• Have you considered actually doing this with some small datasets? There's nothing like working it out for yourself to really get the point. Sketching suitable graphs of the data and the mean-substituted data (use histograms, perhaps) can also be instructive. If you are mathematically inclined, you can work out fairly simple formulas for the post-substitution statistics. This is especially easy once you recognize that shifting the mean does not alter SD, skewness, or kurtosis, so you can work with datasets that have means of zero.
– whuber
Jan 21 at 20:52
• I did run a mock dataset on SPSS but I am finding such small differences. Not sure if I am doing it right :( Jan 22 at 0:27
• If the dataset is large, the differences will be small.
– whuber
Jan 22 at 3:40

Let's consider what these statistics are, and then get into what mean substitution does.

### Definitions

Ignoring the missing values for a moment, you have a dataset of numbers $$x_1, x_2, \ldots, x_n,$$ not all of them equal. The SD $$\sigma,$$ skewness $$\beta,$$ and kurtosis $$\kappa$$ of this dataset all depend on the residuals relative to the mean

$$m = (x_1 + x_2 + \cdots + x_n)/n.$$

These residuals are

$$y_i = x_i - m.$$

They sum to zero (because the $$x_i$$ and the $$n$$ subtracted copies of $$m$$ exactly balance). In these terms there are various formulas for things called "standard deviation" (SD), "skewness," and "kurtosis." Because the details don't matter, let's express these formulas as

$$\mu_k = \frac{1}{C_{n,k}} \left(y_1^k + y_2^k + \cdots + y_n^k\right)$$

where $$k=2$$ gives the square of the SD (that is, the variance), $$k=3$$ gives the third central moment, and $$k=4$$ gives the fourth central moment. The denominators $$C_{n,k}$$ increase with $$n,$$ are close to $$n,$$ and (on a relative basis) approach $$n$$ as $$n$$ grows large. See https://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm for a rundown of some of the formulas to see the details I'm avoiding here.

Finally, the statistics in question are

$$SD = \sigma= \sqrt{\mu_2};\quad \beta = \frac{\mu_3}{\sigma^3} = \frac{\mu_3}{\mu_2^{3/2}};\quad \kappa = \frac{\mu_4}{\sigma^4}= \frac{\mu_4}{\mu_2^{2}}.$$

### Analysis

Consider those missing values and suppose there are $$a$$ of them. "Mean substitution" introduces additional "data" $$x_{n+1}=x_{n+2}=\cdots=x_{n+1}=m$$ all equal to the original mean. Obviously the new mean is the same as the old, which implies all the residuals are the same. The new residuals are $$y_{n+i} = x_{n+i}-m = m-m= 0.$$ Consequently, mean substitution does not change any of the sums of powers in the definitions of the central moments $$\mu_k.$$

However, those pesky denominators do change: where originally $$C_{n,k} \approx n$$ is used, with the mean-substituted data $$C_{n+a,k}\approx n+a$$ will be used in its place. As a result,

Mean-substitution of $$a$$ missing values multiplies the original central moment $$\mu_k$$ by $$C_{n,k}/C_{n+a,k}.$$

Notice this implies a decrease in the size of every central moment. This is the crux of the matter.

### Results

Plugging the preceding result into the formulas for the statistics and using primes for the new values, we see

\begin{aligned} \sigma^\prime &= \sqrt{\frac{C_{n,2}}{C_{n+a,2}}}\,\sigma &\approx \sqrt{\frac{n}{n+a}}\,\sigma &\lt \sigma;\\ \beta^\prime &= \frac{C_{n,3}}{C_{n+a,3}}\left(\frac{C_{n+a,2}}{C_{n,2}}\right)^{3/2} \beta &\approx \sqrt{\frac{n+a}{n}}\beta;\\ \kappa^\prime &= \frac{C_{n,4}}{C_{n+a,4}}\left(\frac{C_{n+a,2}}{C_{n,2}}\right)^2 \kappa &\approx \frac{n+a}{n}\kappa & \gt \kappa. \end{aligned}

That is,

• The SD $$\sigma$$ decreases because a bunch of zero residuals have been thrown into the mix;

• So the skewness (if not originally zero) increases because the decrease in $$\sigma^3$$ in the denominator overwhelms the decrease in $$\mu_3$$ from the inclusions of (cubes) of those zero residuals;

• And the kurtosis increases even more from the decrease in $$\sigma^4.$$

By "increase" in skewness I mean an increase in its size: its sign remains the same.

Notice that these are universal patterns on a relative scale: that is, we can predict exactly how $$\sigma^\prime/\sigma,$$ etc. will behave as a function of $$n$$ and $$a$$ once we have selected our denominators $$C_{n,k}$$ for computing these statistics. With $$C_{n,k}=n$$ in all cases (these are the empirical formulas) and $$n=4$$ non-missing values, here is what these plots look like (with positive skewness):

For those who like to experiment, here is the R code to produce the illustration. It is written to parallel the mathematical exposition here rather than for efficiency.

mu <- function(x, k) mean((x - mean(x))^k)             # Central moment of x of order k
moment.std <- function(x, k) mu(x, k) / mu(x, 2)^(k/2) # Standardized central moment of x
sigma <- function(x) sqrt(mu(x, 2))
beta <- function(x) moment.std(x, 3)
kappa <- function(x) moment.std(x, 4)

x <- c(0,1,2,5) # Any set of values of nonzero skewness will serve as an illustration
a <- 0:10       # Number of missing values to examine
stats <- sapply(a, function(a) {
z <- c(x, rep(mean(x), a))                       # Mean substitution of a values
c(Sigma=sigma(z), Beta=beta(z), Kappa=kappa(z))  # The three statistics
})
stats <- stats/stats[, 1] # Relative values of the statistics

par(mfrow=c(1,3))
plot(a, stats["Sigma",], ylab="", main=expression(sigma * "'" /sigma))
plot(a, stats["Beta",], ylab="", main=expression(beta * "'" /beta))
plot(a, stats["Kappa",], ylab="", main=expression(kappa * "'" /kappa))
par(mfrow=c(1,1))