Welch's t-test: Difference between revisions

 

'''Task Description'''<br>

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 program, for example, pvalue.c, can be compiled by

 

 

or

 

 

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

#include <math.h>

 

return 1.0;

return 1.0;

}

}

}

}

}

}

/

(

);

}

return 1.0;

}

 

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 x[] = {3.0,4.0,1.0,2.1};

const double y[] = {490.2,340.0,433.9};

 

return 0;

}

 

{{out}}

<pre>Test sets 1 p-value = 0.021378

Test sets 2 p-value = 0.148842

 

'''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|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|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>

 

=={{header|Python}}==

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]))

=={{header|R}}==

d1 <- c(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)

d2 <- c(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)

x <- c(3.0,4.0,1.0,2.1)

y <- c(490.2,340.0,433.9)

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

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

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

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

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

 

{{out}}

</pre>

 

=={{header|Racket}}==

 

(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|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>