Replace t-test with modified z-test - anyone ever seen this? I just read about a procedure I'd never heard of, and was wondering if anyone has had any experience with it.  In Donald Berry's "Statistics" book, he presents an alternative to doing a t-test by modifying the sample standard deviation, with a factor that increases with lower-N:
$s\rightarrow s\left(1+\frac{20}{n^2}\right)$
and then doing a simple z-test.  Has anyone ever seen this, or seen a derivation of this?  It seems to empirically match the t-distribution test at the 99% level, and overestimates the standard deviation for lower percentiles.
I doubt that it would make a lot of practical difference, but if justified, it might help students in intro stats classes. 
thoughts?
 A: For any given significance level, a simpler formula does better--quite well, in fact.
What's really going on here is that we want to approximate quantiles of the t-distribution of $n-1$ degrees of freedom by adjusting the corresponding quantiles of the standard Normal distribution.  A multiplicative adjustment is a natural one to try (rather than an additive one, especially if you have ever spent time staring at tables of these critical values).  So what we would like to know is how the ratio of quantiles $t_{n,1-\alpha}/z_{1-\alpha}$ varies with $n$ for smallish values of $\alpha$ around $1$% or $5$%.
A good way to obtain the answer is to compute some of these ratios and plot them against $n$.  Using standard methods of exploratory data analysis, I find that 
$$f(n, \alpha) = 1 / \left(t_{n-1,1-\alpha}/z_{1-\alpha} - 1\right)$$
varies in a beautifully linear fashion with $n$ for a wide range of small $n$ (which is where any such adjustment actually matters).  Not only that, the variation is linear regardless of the value of $\alpha$, too.  The evidence is abundantly clear in this plot of $f(n,\alpha)$ versus $n$ for $\alpha$ ranging from $1/8 = 0.125$ down to $1/2^{10}\approx 0.001$ by factors of $1/2$:

In this figure, colors distinguish values of $\alpha$.  The slopes get small as $\alpha$ decreases.
Knowing this, you can pick your favorite value of $\alpha$, such as $\alpha = 5$%, find the formula of the corresponding line (visually if you like or using least squares for more precision), and fiddle a little with the results to obtain a pleasing formula.  For instance, with $\alpha=5$%, an intercept of $-2/3$ and slope of $1$ work fine.  Specifically, this means the standard deviation should be multiplied by $1 + 1/(n-2/3))$.  Here's the evidence:

A plot of the residuals (not shown) indicates this approximation is accurate to about one part in $300$ except for 1 and 2 degrees of freedom, where it errs by $4$% and $-1.5$%, respectively.  This is pretty good considering the Student t quantiles themselves range from $3.84$ down to $1.03$ times the Standard Normal quantile.

Given that the actual behavior of quantiles of the Student t distribution is proportional to $1/n$ rather than $1/n^2$, I see little pedagogical value in using the $1/n^2$ adjustment proposed in the question: at best it might be useful as an approximation, but it provides no valid insight into the relationship between the Student t and standard Normal distributions.

Here is R code used to produce these plots.
#
# Approximate qt by qnorm.
#
y <- function(x, q) {
  sapply(q, function(r) 1/(qt(r, x) / qnorm(r) - 1))
}
#
# Plot the approximation.
#
x <- 1:30                                           # Degrees of freedom
q <- 1 - exp(-seq(log(8), 7, log(2)))               # Quantiles
k <- length(q)
data <- y(x, q)                                     # Approximations
fits <- apply(data, 2, function(z) lm(z ~ x)$coeff) # Fitted lines
colors <- hsv(1 - (1:k - 1/2)/k, .8, .8)
plot(t(matrix(rep(x,k), ncol=k)), t(data), 
     pch=19, col=colors,
     xlab="df", ylab="1/(t/z - 1)")
tmp <- sapply(1:k, function(i) abline(fits[,i], col=colors[i]))
#
# Figure 2: Evaluate the approximation for a particular quantile
#
plot(qt(.95, x),qnorm(.95) *(1 + 1 / (x-2/3)), pch=19, col="Red",
     log="xy", main="Estimated versus actual Student t quantile")
abline(0,1)
#
# Plot the relative residuals
# (The approximation typically is better than 2 sig figs except for 
# 1 and 2 df.)
#
plot(qt(.95, x) / (qnorm(.95) *(1 + 1 / (x-2/3))), pch=19, col="Red",
     log="y", main="Approximation residuals")
abline(h=1)

