4
from scipy import stats 
import numpy as np
ts1 = np.array([11,9,10,11,10,12,9,11,12,9])
ts2 = np.array([11,13,10,13,12,9,11,12,12,11])

r = stats.ttest_ind(ts1, ts2, equal_var=False)
print(r.statistic, r.pvalue)

The null hypothesis is that the averages are equal. This code will give me the t statistic and the P-value.

But what would be a simple way to calculate the 95% confidence interval for the difference in average? Is there a way to do that with scipy?

Florin Andrei
  • 163
  • 1
  • 6

2 Answers2

6

It's a very good detailed answer provided by @BruceET. So if you want to do it python, you have to calculate the pooled standard error. I moved the code from this link and you can see it gives you something similar to Bruce's answer:

import numpy as np
from scipy.stats import ttest_ind
from scipy.stats import t
import pandas as pd

def welch_ttest(x1, x2,alternative):
    
    n1 = x1.size
    n2 = x2.size
    m1 = np.mean(x1)
    m2 = np.mean(x2)
    
    v1 = np.var(x1, ddof=1)
    v2 = np.var(x2, ddof=1)
    
    pooled_se = np.sqrt(v1 / n1 + v2 / n2)
    delta = m1-m2
    
    tstat = delta /  pooled_se
    df = (v1 / n1 + v2 / n2)**2 / (v1**2 / (n1**2 * (n1 - 1)) + v2**2 / (n2**2 * (n2 - 1)))
    
    # two side t-test
    p = 2 * t.cdf(-abs(tstat), df)
    
    # upper and lower bounds
    lb = delta - t.ppf(0.975,df)*pooled_se 
    ub = delta + t.ppf(0.975,df)*pooled_se
  
    return pd.DataFrame(np.array([tstat,df,p,delta,lb,ub]).reshape(1,-1),
                         columns=['T statistic','df','pvalue 2 sided','Difference in mean','lb','ub'])

We run this function, i named the lower and upper bounds of the 95% CI as lb and ub.. You can simply modify them in the function:

welch_ttest(ts1,ts2,"equal")


    T statistic     df  pvalue 2 sided  Difference in mean  lb  ub
0   -1.832542   17.90031    0.08356     -1.0    -2.146912   0.146912
StupidWolf
  • 4,494
  • 3
  • 10
  • 26
  • Isn't pooled standard error only applicable when the variances are assumed to be the same in the two groups? The question implies that variances are different (`equal_var=False`) – Aleksandar Jovanovic Nov 11 '20 at 15:43
  • Thanks for this answer, should the df var be plugged into the lower and upper bounds in the t ppf? – Francisco Feb 07 '21 at 10:00
  • oh yes you are right, i did not put in df in the t.ppf part. it should be used in both upper and lower, like ```t.ppf(0.975,df)*pooled_se``` I have edited my answer. thanks for that. – StupidWolf Feb 07 '21 at 10:13
4

Not sure about Scripy. Maybe there's a Scripy help site that will show the code. [Perhaps this.]

In R, a 95% CI is part of t.test output, where the Welch version of the 2-sample t test is the default (and argument var.eq=T gets you the pooled test).

ts1 = c(11,9,10,11,10,12,9,11,12,9)
ts2 = c(11,13,10,13,12,9,11,12,12,11)
t.test(ts1, ts2)

        Welch Two Sample t-test

data:  ts1 and ts2
t = -1.8325, df = 17.9, p-value = 0.08356
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.1469104  0.1469104
sample estimates:
mean of x mean of y 
     10.4      11.4

Because the 95% CI includes $0$ the 2-sided test does not reject $H_0: \mu_1=\mu_2$ at the 5% level.

The 95% margin of error is $t^*\sqrt{\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2}},$ where $t^*$ cuts probability $0.025=2.5\%$ from the upper tail of Student's t distribution with degrees of freedom $\nu^\prime$ as found from the Welch formula involving sample variances and sample sizes. [Here, $\nu^\prime = 17.9,$ in some software rounded down to an integer. One always has $\min(n_1-1,n_2-1) \le \nu^\prime \le n_1+n_2-2.]$

me = qt(.975, 17.9)*sqrt(var(ts1)/10+var(ts2)/10); me
[1] 1.146912
pm=c(-1,1)
-1 + pm*me
[1] -2.1469118  0.1469118

It's always a good idea to keep the actual formulas in mind, even if one hopes to use them only rarely.

BruceET
  • 47,896
  • 2
  • 28
  • 76