Creating \(t\)-tests for use in R when only summary statistics are available

I. Creating a function to evaluate \(t\)-tests from summary statistics only

Neither native R nor its Mosaic package have versions of the \(t\)-tests for use when only the summary statistics are known. However, we can create them ourselves based on our statistics formulas.

One-Sample \(t\)-test

The one-sample \(t\)-statistic is given by \[t = \frac{\bar x - \mu_0}{\frac{s}{\sqrt{n}}}\] with p-values from the \(t\)-distribution with degrees of freedom \[df = n-1\]

Knowing the formulas, you should be able see how things are calculated in the function below once you realize the the function pt finds the exact p-value based a \(t\)-statistic and degrees of freedom. The function my.t.test returns the two-tailed p-value only. The parameter mu defaults to zero but may be set to match \(\mu_0\).

my.t.test = function(m, s, n, mu = 0)
{
  df = n-1
  se = s/sqrt(n)
  t = (m-mu)/(se)
  # Calculate left-tail p-value, then double it for two-tailed p-value
  p = 2*pt(-abs(t),df)
  dat = c(se, df, t, p)
  names(dat) = c("Std Error", "Deg Freedom", "t", "p-value")
return(dat) 
}

Example 1

Suppose we are testing the reliability of a machine that fills 12 ounce aluminum cans with a soft drink. After sampling, we find the following data:

\[\begin{array}{ccc}\bar x && s && n \\ \hline 11.9 && 0.1 && 50 \end{array}\]

We can test the hypothesis that the machine is underfilling the cans using the 0.05 level.

\[H_0 : \mu = 12 \\ H_0 : \mu < 12 \]

We put the summary statistics in, and make sure to set the value \(\mu_0 = 12\) as our last argument.

my.t.test(11.9, 0.2, 50, mu = 12)

   Std Error  Deg Freedom            t      p-value 
 0.028284271 49.000000000 -3.535533906  0.000900003 

If you run the same test on your TI graphing calculator, you will find the same results if you use a two-tailed test. Here, however, our alternative hypothesis specified a one-tailed test, but we can simply divide our p-value in half to get the correct value. Again, you can check this by comparing to you graphing calculator.

II. Creating a function to evaluate 2-sample \(t\)-tests from summary statistics only

For two-sample testing, the correct statistic is given by

\[t = \frac{\bar x_1 - \bar x_2}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}}\]

with conservative p-values from the \(t\)-distribution with degrees of freedom:

\[df = \min(n_1 - 1 , n_2 - 1)\]

However, most modern statistics apps use the Welch \(t\)-test which should be used whenever we assume unequal variances – which we should in basically all cases. We are testing for unequal means, so assuming equal variances is more than just counter-intuitive. The Welch version of Student’s (William Sealy Gosset’s) \(t\)-test adjusts the degrees of freedom for better accuracy when the standard deviations are, indeed, different. The degrees of freedom calculation is a bit messy, but don’t worry, we’ll just let R do the calculations for us.

\[df = \frac{\left(\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}\right)^2}{\frac{\left(\frac{s_1^2}{n_1}\right)^2}{n_1-1}+\frac{\left(\frac{s_2^2}{n_2}\right)^2}{n_2-1}}\]

Here’s the code. Notice the p-value is calculated by finding the area corresponding to a one-tailed test, and then doubled.

myWelch.t.test = function(m1,m2,s1,s2,n1,n2)
{
  se = sqrt( (s1^2/n1) + (s2^2/n2) )
  # welch-satterthwaite degrees of freedom
  df = ( (s1^2/n1 + s2^2/n2)^2 )/( (s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1) )
  t = (m1-m2)/se
  # Calculate left-tail p-value, then double it for two-tailed p-value
  p = 2*pt(-abs(t),df)
  dat = c(m1-m2, se, df, t, p)
  names(dat) = c("Mean Diff", "Std Error", "Deg Freedom", "t", "p-value")
return(dat) 
}

Example 2

Who’s more optimistic? A recent survey measured optimism levels and asked participants their biological sex. Test for a significant difference in Optimism based upon Sex at the 0.05 level given the data table below.

\[\begin{array}{lcc} &\textbf{Females} & \textbf{Males}\\ \hline \textbf{AVG} & 28.5 & 30.0 \\ \textbf{SD} & 4.9 & 5.7 \\ \textbf{N} & 130 & 74 \end{array}\]

We simply plug in the values to our Welch \(t\)-test function.

myWelch.t.test(28.5, 30, 4.9, 5.9, 130, 74)

   Mean Diff    Std Error  Deg Freedom            t      p-value 
 -1.50000000   0.80938107 130.21690377  -1.85326796   0.06610732 

Since the alternative hypothesis specified a two-tailed, the p-value given in the table is correct, and we fail to reject the null.

III. How to use these functions

If you would like to use these functions, just open this notebook in RStudio and run the code blocks associated with the two \(t\)-test functions. Executing their code blocks will load them into R’s working memory for the current session. Return to your notebook, and implement them as shown in the examples.

You do not need to copy the code blocks into your working notebook. Just execute them from this notebook during any session where you need them.