Welch's t-test

Given two lists of data, calculate the p-Value^[1] used for null hypothesis testing.

This uses Welch's t-test, which assumes that the variances between the two sets are not equal. Welch's t-test statistic can be computed thus:

$t\quad =\quad {\;{\overline {X}}_{1}-{\overline {X}}_{2}\; \over {\sqrt {\;{s_{1}^{2} \over N_{1}}\;+\;{s_{2}^{2} \over N_{2}}\quad }}}$

and the degrees of freedom, $\nu$ can be approximated:

$\nu \quad \approx \quad {{\left(\;{s_{1}^{2} \over N_{1}}\;+\;{s_{2}^{2} \over N_{2}}\;\right)^{2}} \over {\quad {s_{1}^{4} \over N_{1}^{2}\nu _{1}}\;+\;{s_{2}^{4} \over N_{2}^{2}\nu _{2}}\quad }}$

The p-value, $p$ , can be computed as a cumulative distribution function

$p=1-{\frac {1}{2}}I_{\frac {\nu }{t^{2}+\nu }}\left({\frac {\nu }{2}},{\frac {1}{2}}\right)$

where I is the incomplete beta function. This is the same as:

$p=1-{\frac {1}{2}}{\frac {B({\frac {\nu }{t^{2}+\nu }};{\frac {\nu }{2}},{\frac {1}{2}})}{B({\frac {\nu }{2}},{\frac {1}{2}})}}$

Keeping in mind that

$B(x;a,b)=\int _{0}^{x}t^{a-1}\,(1-t)^{b-1}\,\mathrm {d} t.\!$

and

$\mathrm {B} (x,y)={\dfrac {\Gamma (x)\,\Gamma (y)}{\Gamma (x+y)}}=\exp(\ln {\dfrac {\Gamma (x)\,\Gamma (y)}{\Gamma (x+y)}})=\exp((\ln(\Gamma (x))+\ln(\Gamma (y))-\ln(\Gamma (x+y)))\!$

$p$ can be calculated in terms of gamma functions and integrals more simply:

$p=1-{\frac {1}{2}}\times {\frac {\int _{0}^{\frac {\nu }{t^{2}+\nu }}r^{{\frac {\nu }{2}}-1}\,(1-r)^{-0.5}\,\mathrm {d} r}{\exp((\ln(\Gamma ({\frac {\nu }{2}}))+\ln(\Gamma (0.5))-\ln(\Gamma ({\frac {\nu }{2}}+0.5)))}}$

which simplifies to

$p=1-{\frac {1}{2}}\times {\frac {\int _{0}^{\frac {\nu }{t^{2}+\nu }}{\frac {r^{{\frac {\nu }{2}}-1}}{\sqrt {1-r}}}\,\mathrm {d} r}{\exp((\ln(\Gamma ({\frac {\nu }{2}}))+\ln(\Gamma (0.5))-\ln(\Gamma ({\frac {\nu }{2}}+0.5)))}}$

C

Works with: C99

This shows how pvalue can be calculated from any two arrays, using Welch's 2-sided t-test, which doesn't assume equal variance. The high "n" value ensures roughly 6 digits convergence with Simpson's integral approximation. <lang C>#include <stdio.h>

include <math.h>

double calculate_Pvalue (const double *array1, const size_t array1_size, const double *array2, const size_t array2_size) { if (array1_size <= 1) { return 1.0; } if (array2_size <= 1) { return 1.0; } double mean1 = 0.0, mean2 = 0.0; for (size_t x = 0; x < array1_size; x++) { mean1 += array1[x]; } for (size_t x = 0; x < array2_size; x++) { mean2 += array2[x]; } if (mean1 == mean2) { return 1.0; } mean1 /= array1_size; mean2 /= array2_size; double variance1 = 0.0, variance2 = 0.0; for (size_t x = 0; x < array1_size; x++) { variance1 += (array1[x]-mean1)*(array1[x]-mean1); } for (size_t x = 0; x < array2_size; x++) { variance2 += (array2[x]-mean2)*(array2[x]-mean2); } if ((variance1 == 0.0) && (variance2 == 0.0)) { return 1.0; } variance1 = variance1/(array1_size-1); variance2 = variance2/(array2_size-1); const double WELCH_T_STATISTIC = (mean1-mean2)/sqrt(variance1/array1_size+variance2/array2_size); const double DEGREES_OF_FREEDOM = pow((variance1/array1_size+variance2/array2_size),2.0)//numerator / ( (variance1*variance1)/(array1_size*array1_size*(array1_size-1))+ (variance2*variance2)/(array2_size*array2_size*(array2_size-1)) ); const double a = DEGREES_OF_FREEDOM/2, x = DEGREES_OF_FREEDOM/(WELCH_T_STATISTIC*WELCH_T_STATISTIC+DEGREES_OF_FREEDOM); const unsigned short int N = 65535; const double h = x/N; double sum1 = 0.0, sum2 = 0.0; for(unsigned short int i = 0;i < N; i++) {

     sum1 += (pow(h * i + h / 2.0,a-1))/(sqrt(1-(h * i + h / 2.0)));
     sum2 += (pow(h * i,a-1))/(sqrt(1-h * i));

} double return_value = ((h / 6.0) * ((pow(x,a-1))/(sqrt(1-x)) + 4.0 * sum1 + 2.0 * sum2))/(expl(lgammal(a)+0.57236494292470009-lgammal(a+0.5)));//0.5723... is lgammal(0.5) if ((isinf(return_value) != 0) || (return_value > 1.0)) { return 1.0; } else { return return_value; } } //------------------- int main(void) { const double d1[] = {27.5,21.0,19.0,23.6,17.0,17.9,16.9,20.1,21.9,22.6,23.1,19.6,19.0,21.7,21.4}; const double d2[] = {27.1,22.0,20.8,23.4,23.4,23.5,25.8,22.0,24.8,20.2,21.9,22.1,22.9,20.5,24.4}; const double d3[] = {17.2,20.9,22.6,18.1,21.7,21.4,23.5,24.2,14.7,21.8}; const double d4[] = {21.5,22.8,21.0,23.0,21.6,23.6,22.5,20.7,23.4,21.8,20.7,21.7,21.5,22.5,23.6,21.5,22.5,23.5,21.5,21.8}; const double d5[] = {19.8,20.4,19.6,17.8,18.5,18.9,18.3,18.9,19.5,22.0}; const double d6[] = {28.2,26.6,20.1,23.3,25.2,22.1,17.7,27.6,20.6,13.7,23.2,17.5,20.6,18.0,23.9,21.6,24.3,20.4,24.0,13.2}; const double d7[] = {30.02,29.99,30.11,29.97,30.01,29.99}; const double d8[] = {29.89,29.93,29.72,29.98,30.02,29.98}; const double x[] = {3.0,4.0,1.0,2.1}; const double y[] = {490.2,340.0,433.9};

printf("Test sets 1 p-value = %lf\n",calculate_Pvalue(d1,sizeof(d1)/sizeof(*d1),d2,sizeof(d2)/sizeof(*d2))); printf("Test sets 2 p-value = %lf\n",calculate_Pvalue(d3,sizeof(d3)/sizeof(*d3),d4,sizeof(d4)/sizeof(*d4))); printf("Test sets 3 p-value = %lf\n",calculate_Pvalue(d5,sizeof(d5)/sizeof(*d5),d6,sizeof(d6)/sizeof(*d6))); printf("Test sets 4 p-value = %lf\n",calculate_Pvalue(d7,sizeof(d7)/sizeof(*d7),d8,sizeof(d8)/sizeof(*d8))); printf("Test sets 5 p-value = %lf\n",calculate_Pvalue(x,sizeof(x)/sizeof(*x),y,sizeof(y)/sizeof(*y))); return 0; } </lang>

Output:

Test sets 1 p-value = 0.021378
Test sets 2 p-value = 0.148842
Test sets 3 p-value = 0.035972
Test sets 4 p-value = 0.090773
Test sets 5 p-value = 0.010751

Racket

Translation of: C

, producing the same output.

<lang racket>#lang racket (require math/statistics math/special-functions)

(define (p-value S1 S2 #:n (n 11000))

 (define σ²1 (variance S1 #:bias #t))
 (define σ²2 (variance S2 #:bias #t))
 (define N1 (sequence-length S1))
 (define N2 (sequence-length S2))
 (define σ²/sz1 (/ σ²1 N1))
 (define σ²/sz2 (/ σ²2 N2))
 
 (define degrees-of-freedom
   (/ (sqr (+ σ²/sz1 σ²/sz2))
      (+ (/ (sqr σ²1) (* (sqr N1) (sub1 N1)))
         (/ (sqr σ²2) (* (sqr N2) (sub1 N2))))))
 
 (define a (/ degrees-of-freedom 2))
 (define a-1 (sub1 a))
 (define x (let ((welch-t-statistic (/ (- (mean S1) (mean S2)) (sqrt (+ σ²/sz1 σ²/sz2)))))
             (/ degrees-of-freedom (+ (sqr welch-t-statistic) degrees-of-freedom))))
 (define h (/ x n))
 
 (/ (* (/ h 6)
       (+ (* (expt x a-1)
             (expt (- 1 x) -1/2))
          (* 4 (for/sum ((i (in-range 0 n)))
                 (* (expt (+ (* h i) (/ h 2)) a-1)
                    (expt (- 1 (+ (* h i) (/ h 2))) -1/2))))
          (* 2  (for/sum ((i (in-range 0 n)))
                  (* (expt (* h i) a-1) (expt (- 1 (* h i)) -1/2))))))
    (* (gamma a) 1.77245385090551610 (/ (gamma (+ a 1/2))))))

(module+ test

 (list
  (p-value (list 27.5 21.0 19.0 23.6 17.0 17.9 16.9 20.1 21.9 22.6 23.1 19.6 19.0 21.7 21.4)
           (list 27.1 22.0 20.8 23.4 23.4 23.5 25.8 22.0 24.8 20.2 21.9 22.1 22.9 20.5 24.4))
  
  (p-value (list 17.2 20.9 22.6 18.1 21.7 21.4 23.5 24.2 14.7 21.8)
           (list 21.5 22.8 21.0 23.0 21.6 23.6 22.5 20.7 23.4 21.8
                 20.7 21.7 21.5 22.5 23.6 21.5 22.5 23.5 21.5 21.8))
  
  (p-value (list 19.8 20.4 19.6 17.8 18.5 18.9 18.3 18.9 19.5 22.0)
           (list 28.2 26.6 20.1 23.3 25.2 22.1 17.7 27.6 20.6 13.7
                 23.2 17.5 20.6 18.0 23.9 21.6 24.3 20.4 24.0 13.2))
  
  (p-value (list 30.02 29.99 30.11 29.97 30.01 29.99)
           (list 29.89 29.93 29.72 29.98 30.02 29.98))
  
  (p-value (list 3.0 4.0 1.0 2.1)
           (list 490.2 340.0 433.9))))</lang>

Output:

(0.021378001462867013 0.14884169660532798 0.035972271029796624 0.09077332428567102 0.01075139991904718)

Tcl

Translation of: Racket

Works with: Tcl version 8.6

Library: Tcllib (Package: math::statistics)

Library: Tcllib (Package: math::special)

This is not particularly idiomatic Tcl, but perhaps illustrates some of the language's relationship with the Lisp family.

<lang Tcl>#!/usr/bin/tclsh

package require math::statistics package require math::special namespace path {::math::statistics ::math::special ::tcl::mathfunc ::tcl::mathop}

proc incf {_var {inc 1.0}} {

   upvar 1 $_var var
   if {![info exists var]} {
       set var 0.0
   }
   set var [expr {$inc + $var}]

}

proc sumfor {_var A B body} {

   upvar 1 $_var var
   set var $A
   set res 0
   while {$var < $B} {
       incf res [uplevel 1 $body]
       incr var
   }
   return $res

}

proc sqr {x} {expr {$x*$x}}

proc pValue {S1 S2 {n 11000}} {

   set σ²1 [var $S1]
   set σ²2 [var $S2]
   set N1  [llength $S1]
   set N2  [llength $S2]
   set σ²/sz1 [/ ${σ²1} $N1]
   set σ²/sz2 [/ ${σ²2} $N2]

   set d1 [/ [sqr ${σ²1}] [* [sqr $N1] [- $N1 1]]]
   set d2 [/ [sqr ${σ²2}] [* [sqr $N2] [- $N2 1]]]
   set DoF [/ [sqr [+ ${σ²/sz1} ${σ²/sz2}]] [+ $d1 $d2]]

   set a [/ $DoF 2.0]

   set welchTstat [/ [- [mean $S1] [mean $S2]] [sqrt [+ ${σ²/sz1} ${σ²/sz2}]]]
   set x [/ $DoF [+ [sqr $welchTstat] $DoF]]
   set h [/ $x $n]

   / [* [/ $h 6] \
        [+ [* [** $x [- $a 1]] \
              [** [- 1 $x] -0.5]] \
           [* 4 [sumfor i 0 $n {
                   * [** [+ [* $h $i] [/ $h 2]] [- $a 1]] \
                     [** [- 1 [* $h $i] [/ $h 2]] -0.5]}]] \
           [* 2 [sumfor i 0 $n {
                   * [** [* $h $i] [- $a 1]] [** [- 1 [* $h $i]] -0.5]}]]]] \
     [* [Gamma $a] 1.77245385090551610 [/ 1.0 [Gamma [+ $a 0.5]]]]

}

foreach {left right} {

   { 27.5 21.0 19.0 23.6 17.0 17.9 16.9 20.1 21.9 22.6 23.1 19.6 19.0 21.7 21.4 }
   { 27.1 22.0 20.8 23.4 23.4 23.5 25.8 22.0 24.8 20.2 21.9 22.1 22.9 20.5 24.4 }

   { 17.2 20.9 22.6 18.1 21.7 21.4 23.5 24.2 14.7 21.8 }
   { 21.5 22.8 21.0 23.0 21.6 23.6 22.5 20.7 23.4 21.8 20.7 21.7 21.5 22.5 23.6 21.5 22.5 23.5 21.5 21.8 }

   { 19.8 20.4 19.6 17.8 18.5 18.9 18.3 18.9 19.5 22.0 }
   { 28.2 26.6 20.1 23.3 25.2 22.1 17.7 27.6 20.6 13.7 23.2 17.5 20.6 18.0 23.9 21.6 24.3 20.4 24.0 13.2 }

   { 30.02 29.99 30.11 29.97 30.01 29.99 }
   { 29.89 29.93 29.72 29.98 30.02 29.98 }

   { 3.0 4.0 1.0 2.1 }
   { 490.2 340.0 433.9 }

} {

   puts [pValue $left $right]

} </lang>

Output:

0.021378001462853034
0.148841696604164
0.035972271029770915
0.09077332428458083
0.010751399918798182

↑ Template:Cite doi

[nature506-1] Template:Cite doi

[1]