2
$\begingroup$

In R, how can I perform a two sample t-test in the case of equal variances?

My actual problem is to compare different algorithms for length measurements in machine vision.

I have studied Applied Statistics and Probability for Engineers by Montgomery & Runger and I have tried to implement the test in R and compare the results to the R command t.test.

This is my code:

# Douglas C. Montgomery, George C. Runger
# Applied Statistics and Probability for Engineers
# Third Edition
# 10-3.1 Hypotheses Tests for a Difference in Means, Variances Unknown, p.337
# Case 1: sigma_1^2 = sigma_2^2 = sigma^2
#
# See also: EXAMPLE 10-5, p.339
two_sample_t_test_equal_variance <- function(x1,x2,Delta_0,alpha)
{
   n1 <- length(x1)
   n2 <- length(x2)
   dof <- n1+n2-2
   S1_squared <- var(x1)
   S2_squared <- var(x2)
   #pooled estimator:
   Sp_squared <- ((n1-1)*S1_squared + (n2-1)*S2_squared)/(n1+n2-2) 
   x1_bar <- mean(x1)
   x2_bar <- mean(x2)
   Sp <- sqrt(Sp_squared)
   t0 <- (x1_bar - x2_bar - Delta_0)/(Sp*sqrt(1/n1+1/n2))
   t_half_alpha_dof <- -qt(alpha/2,dof)
   reject_H0 <- (t0 > t_half_alpha_dof || t0 < - t_half_alpha_dof)
   if ( reject_H0 ) {
      cat("Reject H0 (alpha =",alpha,").\nH0 is mu1 - mu2 =", Delta_0, "\n\n" )
   } else {
      cat("Cannot reject H0 (alpha =",alpha,").\nH0 is mu1 - mu2 =", Delta_0, "\n\n" )
   }

   test_result <- t.test(x1, x2, alternative="two.sided", mu=Delta_0, 
                          paired=FALSE, var.equal = TRUE )

   accept_H0_by_R <- (test_result$p.value > alpha)
   if ( reject_H0 != !accept_H0_by_R ) {
      cat("WARNING: R t.test gives a different answer for accepting H0.\n" )
   }

   rel_error_warn <- 0.01
   if ( abs(test_result$statistic - t0)/abs(test_result$statistic) > rel_error_warn ) {
      cat("WARNING: t0 relative error is >", rel_error_warn, "\n"  )
   }
}

# data from EXAMPLE 10-5, p.339
cat1<-c(91.50, 94.18, 92.18, 95.39, 91.79, 89.07, 94.72, 89.21)
cat2<-c(89.19, 90.95, 90.46, 93.21, 97.19, 97.04, 91.07, 92.75)
two_sample_t_test_equal_variance(cat1,cat2,0,0.05)

Is my invocation of t.test correct for the case at hand?

Is there any error in the code?

I tagged the question with "homework" even if I am not at school since more than 10 years and I am self-studying statistics for job-related reasons.

Thank you. Alessandro

$\endgroup$
  • $\begingroup$ Are you aware that there is already a function in R doing that ? Try help(t.test). $\endgroup$ – steffen Jun 29 '11 at 10:15
  • $\begingroup$ @steffen Yes, I am aware. I included in the code my invocation of t.test and I would like to know whether it is correct. $\endgroup$ – Alessandro Jacopson Jun 29 '11 at 10:22
  • $\begingroup$ @uvts_csv Sorry, did not read the question properly :( $\endgroup$ – steffen Jun 29 '11 at 10:29
2
$\begingroup$

If you were still in school then you might have a couple points taken off for saying that you accept the null hypothesis. The purists (well frequentist purists) will always say that we never accept the null, just fail to reject it.

As a style thing, generally test functions will return an object with the test statistic, p-value, etc. and not print anything. Then a print method will be used to print the results nicely. But for learning or simple use what you have done is fine (I would go the other route if you are planning on building on this, or doing more of your own tests).

In your last cat statement you use relative_error_warning, but I don't see it defined anywhere, did you mean rel_error_warn? Not that that line is ever likely to be run.

Everything looks correct, both your calculations and the running of t.test, I would expect the only differences you ever see to be rounding error (or due to handling of missing values).

$\endgroup$
  • $\begingroup$ +1, thank you very much for the "frequentist purists" clarification. I fixed the last cat statement. $\endgroup$ – Alessandro Jacopson Jun 30 '11 at 7:00
  • $\begingroup$ I also changed "Accept null" into "Cannot reject". $\endgroup$ – Alessandro Jacopson Jun 30 '11 at 19:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.