confidence interval for 2-sample t test with scipy

Question

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?

EDIT: Scipy is not mandatory. Python is mandatory.

Answer 1

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

Answer 2

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.

Answer 3

I found my own question from a couple years ago, and let me add now a very simple answer in Python. There is no need to stick to Scipy - I didn't know that back then. Instead, install the Pingouin library.

https://pingouin-stats.org/

Then:

import pingouin as pg
res = pg.ttest(ts1, ts2, paired=False)
print(res)

Output:

               T  dof alternative     p-val          CI95%   cohen-d   BF10  \
T-test -1.832542   18   two-sided  0.083467  [-2.15, 0.15]  0.819538  1.225   

           power  
T-test  0.411029