# Heteroskedasticity and standard deviation

I am looking at state-wide data (entire population) of a school's grade as a function of the school's poverty index. The data appears to me to be an unconditional heteroskedastic distribution. I am drawing a line of best fit (LBF) through the data using linear regression (I changed this from 2-degree poly). What I am attempting to do is look at how each school is doing compared to how it was ‘predicted’ to do. In past years, a different test was used and the result was that the data was clearly linear, and moderately heteroskedastic. I then used SD as a measure so that I could look across multiple years.

The question I have is how can I find SD across a LBF through heteroskedastic data? I just clustered the data based on Chernick's reply. You can see what it now looks like. Very interesting.

-
 Could you post a plot of the data? – jbowman Aug 15 '12 at 18:25

So y is school grade and x is poverty index. If you have repeated x values or several x values close together you could compute the average squared distance from the observed y to the fitted y is each group. These could serve as variance estimates for the residuals that you can plot vs x (picking a central x for each group) to see how the variance is changing with x.

-
 Correct on XY variables. Your suggestion is what I was thinking. I was wondering if there was a specified approach. As you might expect, the data points start off close together between X=15-40, with Y between 90-100 (low poverty have higher scores). There appears to be another section that spreads apart more between X=40-75, and an even larger spread between X=75-100. But it seems really smooth – Gregory MacDougall Aug 15 '12 at 17:39 I do not know of any formal approach to this. – Michael Chernick Aug 15 '12 at 17:48 Let's see if I can upload the graphic. I divided the data into 3 sections – Gregory MacDougall Aug 27 '12 at 16:20

The shape of the heteroscedasticity in the scatterplot is a strong indicator that a nonlinear re-expression of $y$, the school grade point total, will cure the problem.

A good choice at the outset would be the "folded logarithm" or, equivalently, a logit-like formula

$$z = \log(y / (100 + \varepsilon - y))$$

where $\varepsilon$ is a small positive value, perhaps around $1/2$ or $1$, chosen to avoid singular calculations.

To demonstrate this, I simulated data like the ones shown in the question. The x-coordinates appear concentrated at the right side, so (in R) I created 500 values by scaling a Beta$(4,1)$ distribution to the range $[0,100]$:

set.seed(17)
x <- 100 * rbeta(500, 4, 1)


Next I posited that the logits of the school scores, $z$, would have a linear relationship to the poverty indexes $x$ plus normally-distributed random error. I chose an intercept and slope to reproduce the apparent fit and a standard deviation of the error to reproduce the scatter:

z <- 4 - .025*x + rnorm(length(x), sd=1)
logistic <- function(y) 1 / (1 + exp(-y))
y <- 100 * logistic(z)


Let's look at the plot to see how successful these formulae are at emulating the data:

par(mfrow=c(1,1))
plot(x,y, main="Simulated Data")


That looks quite similar to the plot included in the question.

Now let's display the re-expressed scores:

logit <- function(x) log(x/(100-x))
plot(x, logit(y), ylab="Logit(y)", main="Re-expressed Data")
lines(lowess(logit(y) ~ x), col="Red", lwd=2)


To the uninitiated, this scatterplot may appear heteroscedastic, but that is an illusion brought about by the concentration of $x$ values at the right: as is clear from the code, these data are perfectly homoscedastic. To make this clear, construct a wandering schematic plot by slicing $x$ at key quantiles and displaying side-by-side boxplots of the residuals:

thresh <- c(0, 1/16, 1/8, 1/4, 1/2, 3/4, 7/8, 15/16, 1)
boxplot(residuals(lm(logit(y) ~ x)) ~  cut(x, quantile(x, probs=thresh)),
main="Residuals of logit(y)~x")


These plots show the residuals have essentially the same dispersion across the entire range of $x$ values, demonstrating the homoscedasticity of the re-expressed scores.

Finally--to drive these points home--let's compare least-squares fits of the original scores and of their re-expressed values. First, the original scores:

summary(lm(y ~ x))

Estimate Std. Error t value Pr(>|t|)
(Intercept) 101.8218     2.6770   38.04  < 2e-16 ***
x            -0.2154     0.0327   -6.58  1.2e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12 on 498 degrees of freedom
Multiple R-squared: 0.0801, Adjusted R-squared: 0.0782
F-statistic: 43.4 on 1 and 498 DF,  p-value: 1.16e-10


Standard diagnostic plots (not shown here) are terrible, showing strongly heteroscedastic residuals, non-normal residuals, etc.

Now for the re-expressed data:

summary(lm(logit(y) ~ x))

Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.75714    0.21751   17.27  < 2e-16 ***
x           -0.02157    0.00266   -8.12  3.8e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.975 on 498 degrees of freedom
Multiple R-squared: 0.117,  Adjusted R-squared: 0.115
F-statistic: 65.9 on 1 and 498 DF,  p-value: 3.76e-15


The original coefficients of $4$ and $-0.025$ are adequately estimated as $3.75$ and $-0.022$, respectively. The residual standard error of $0.975$ is close to the intended SD of $1$. All diagnostics are beautiful.

-
 These analyses (except for the boxplots) are straightforward to do in Excel, too: but any statistical package is easier, faster, more reliable, and offers many more analytical capabilities. – whuber♦ Aug 27 '12 at 22:10