Welch's t-test: Difference between revisions

Your task is to discern whether or not the difference in means between the two sets is statistically significant and worth further investigation. P-values are significance tests to gauge the probability that the difference in means between two data sets is significant, or due to chance. A threshold level, alpha, is usually chosen, 0.01 or 0.05, where p-values below alpha are worth further investigation and p-values above alpha are considered not significant. The p-value is not considered a final test of significance, [http://www.nature.com/news/scientific-method-statistical-errors-1.14700 only whether the given variable should be given further consideration].

 

 

<math>t \quad = \quad {\; \overline{X}_1 - \overline{X}_2 \; \over \sqrt{ \; {s_1^2 \over N_1} \; + \; {s_2^2 \over N_2} \quad }} </math>

 

The <math>\ln(\Gamma(x))</math>, or <code>lgammal(x)</code> function is necessary for the program to work with large <code>a</code> values, as [http://rosettacode.org/wiki/Gamma_function Gamma functions] can often return values larger than can be handled by <code>double</code> or <code>long double</code> data types. The <code>lgammal(x)</code> function is standard in <code>math.h</code> with C99 and C11 standards.

 

=={{header|C}}==

This shows how pvalue can be calculated from any two arrays, using Welch's 2-sided t-test, which doesn't assume equal variance.

This is the equivalent of R's<code>t.test(vector1,vector2, alternative="two.sided", var.equal=FALSE)</code> and as such, it is compared against R's pvalues with the same vectors/arrays to show that the differences are very small (here 10^-14).

#include <math.h>

#include <stdlib.h>

return 0;

}

 

{{out}}

'''If''' your computer does not have <code>lgammal</code>, add the following function before <code>main</code> and replace <code>lgammal</code> with <code>lngammal</code> in the <code>calculate_Pvalue</code> function:

 

#include <math.h>

 

}

 

 

=={{header|Fortran}}==

Alternatively, the program shows the p-value computed using the IMSL '''BETAI''' function.

 

use tdf_int

implicit none

print *, t, df, p

print *, betai(df / (t**2 + df), 0.5d0 * df, 0.5d0)

 

'''Output'''

=== Using SLATEC ===

 

 

implicit none

integer :: n1, n2

call welch_ttest(4, x, 3, y, t, df, p)

print *, t, df, p

 

'''Output'''

 

<pre> -9.55949772193266 2.00085234885628 1.075156114978449E-002</pre>

 

=={{header|Go}}==

 

import (

func(r float64) float64 { return math.Pow(r, ν/2-1) / math.Sqrt(1-r) }) /

math.Exp(g1+g2-g3)

{{out}}

<pre>

Implementation:

 

'a b steps'=. 3{.y,128

size=. (b - a)%steps

hi=. v%(t^2)+v

(F f. simpson integrate lo,hi) % 0.5 B v%2

 

<code>integrate</code> and <code>simpson</code> are from the [[Numerical_integration#J|Numerical integration]] task.

 

Data for task examples:

d2=: 27.1 22 20.8 23.4 23.4 23.5 25.8 22 24.8 20.2 21.9 22.1 22.9 20.5 24.4

d3=: 17.2 20.9 22.6 18.1 21.7 21.4 23.5 24.2 14.7 21.8

d8=: 29.89 29.93 29.72 29.98 30.02 29.98

d9=: 3 4 1 2.1

 

Task examples:

0.021378

d3 p2_tail d4

0.0907733

d9 p2_tail da

 

=={{header|Java}}==

Using the '''[http://commons.apache.org/proper/commons-math/ Apache Commons Mathematics Library]'''.

 

public class WelchTTest {

System.out.println("p = " + res[2]);

}

 

'''Result'''

df = 2.0008523488562844

p = 0.010751561149784485</pre>

 

=={{header|Julia}}==

{{works with|Julia|0.6}}

 

 

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]

ttest = UnequalVarianceTTest(y1, y2)

println("\nData:\n y1 = $y1\n y2 = $y2\nP-value for unequal variance TTest: ", round(pvalue(ttest), 4))

 

{{out}}

=={{header|Kotlin}}==

This program brings in code from other tasks for gamma functions and integration by Simpson's rule as Kotlin doesn't have these built-in:

 

typealias Func = (Double) -> Double

println(f.format(p2Tail(da7, da8)))

println(f.format(p2Tail(x, y)))

 

{{out}}

0.010751

</pre>

 

=={{header|Octave}}==

{{trans|Stata}}

y = [490.2,340.0,433.9];

n1 = length(x);

ans =

 

 

=={{header|PARI/GP}}==

B3(x,a,b)=a--;b--;intnum(r=0,x,r^a*(1-r)^b)

Welch2(u,v)=my(m1=vecsum(u)/#u, m2=vecsum(v)/#v, v1=var(u,m1), v2=var(v,m2), s=v1/#u+v2/#v, t=(m1-m2)/sqrt(s), nu=s^2/(v1^2/#u^2/(#u-1)+v2^2/#v^2/(#v-1))); B3(nu/(t^2+nu),nu/2,1/2)/B2(nu/2,1/2);

{{out}}

<pre>%1 = 0.010751561149784496723954539777213062928</pre>

 

=={{header|Perl}}==

{{trans|Sidef}}

use List::Util qw(sum);

use Math::AnyNum qw(gamma pi);

 

my ($A, $B) = @_;

 

my ($left, $right) = splice(@tests, 0, 2);

print p_value($left, $right), "\n";

{{out}}

<pre>

</pre>

 

{{trans|C}}

 

}

 

 

 

 

 

 

 

 

 

 

 

 

{{out}}

</pre>

 

 

=== Using NumPy & SciPy ===

import scipy as sp

import scipy.stats

 

welch_ttest(np.array([3.0, 4.0, 1.0, 2.1]), np.array([490.2, 340.0, 433.9]))

 

else:

temp = qq - ai

if ns == 0:

rx = xx

temp = qq - ai

if ns == 0:

else:

temp = psq

 

 

 

 

 

 

 

=={{header|R}}==

 

printf <- function(...) cat(sprintf(...))

results <- t.test(z1,z2, alternative="two.sided", var.equal=FALSE)

printf("%.15g\n", results$p.value);

 

The output here is used to compare against C's output above.

{{trans|C}}

 

(require math/statistics math/special-functions)

 

(p-value (list 3.0 4.0 1.0 2.1)

 

{{out}}

<pre>(0.021378001462867013 0.14884169660532798 0.035972271029796624 0.09077332428567102 0.01075139991904718)</pre>

 

=={{header|SAS}}==

{{trans|Stata}}

input value group @@;

cards;

class group;

var value;

 

'''Output'''

Implementation in IML:

 

use tbl;

read all var {value} into x where(group=1);

p = 2*probt(-abs(t), df);

print t df p;

 

'''Output'''

 

=={{header|Scala}}==

 

object WelchTTest extends App {

println(s"\nSuccessfully completed without errors. [total ${scala.compat.Platform.currentTime - executionStart} ms]")

 

{{trans|Stata}}

 

Scilab will print a warning because the number of degrees of freedom is not an integer. However, the underlying implementation makes use of the [http://www.netlib.org/random/ dcdflib] Fortran library, which happily accepts a noninteger df.

 

y = [490.2,340.0,433.9];

n1 = length(x);

df = (v1/n1+v2/n2)^2/(v1^2/(n1^2*(n1-1))+v2^2/(n2^2*(n2-1)));

[p, q] = cdft("PQ", -abs(t), df);

 

'''Output'''

 

=={{header|Sidef}}==

[A.len, B.len].all { _ > 1 } || return 1

 

tests.each_slice(2, {|left, right|

say p_value(left, right)

{{out}}

<pre>

Notice that here we use the option '''unequal''' of the '''ttest''' command, and not '''welch''', so that Stata uses the Welch-Satterthwaite approximation.

 

clear

svmat double a

 

di r(p)

 

The computation can easily be implemented in Mata. Here is how to compute the t statistic (t), the approximate degrees of freedom (df) and the p-value (p).

 

x = select(a[., 1], a[., 2] :== 1)

y = select(a[., 1], a[., 2] :== 2)

+----------------------------------------------+

1 | -9.559497722 2.000852349 .0107515611 |

 

=={{header|Tcl}}==

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

 

 

package require math::statistics

puts [pValue $left $right]

}

 

{{out}}

0.09077332428458083

0.010751399918798182

</pre>

 

=={{header|zkl}}==

{{trans|C}}

if (array1.len()<=1 or array2.len()<=1) return(1.0);

 

foreach x in (cof){ ser+=(x/(y+=1)); }

return((2.5066282746310005 * ser / x).log() - tmp);

T(T(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),

T(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)),

 

foreach x,y in (testSets)

{{out}}

<pre>