Skip to content
Snippets Groups Projects
Commit 3ab2279c authored by rhaase's avatar rhaase
Browse files

added another potential solution (which is not working)

parent 91fd7ece
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
11_Hypothesis_testing/exercise_otsu/Otsu_Method_Comparison.ipynb
+ 95
− 11
View file @ 3ab2279c
%% Cell type:code id: tags:
``` python
import pandas as pd
# Load data
table_otsu_python = pd.read_csv("otsu_python.csv", delimiter=',')
table_otsu_imagej = pd.read_csv("otsu_imagej.csv", delimiter=',')
# read out columns
measurement_python = table_otsu_python["# Area"];
measurement_imagej = table_otsu_imagej["Area"];
print(measurement_python)
print(measurement_imagej)
```
%% Output
0 199.127197
1 469.055176
2 782.623291
3 990.295410
4 1100.158691
5 1150.970459
6 1173.553467
7 1184.692383
8 1190.338135
9 1185.302734
10 1189.270020
11 1192.169189
12 1206.970215
13 1242.980957
14 1277.008057
15 1309.814453
16 1334.838867
17 1371.765137
18 1375.122070
19 1315.155029
20 1172.790527
21 917.663574
22 587.005615
23 224.304199
24 128.479004
Name: # Area, dtype: float64
0 197.601
1 468.292
2 780.640
3 988.770
4 1099.396
5 1148.682
6 1172.485
7 1183.167
8 1188.507
9 1183.472
10 1187.134
11 1189.117
12 1205.750
13 1238.708
14 1275.330
15 1307.678
16 1331.329
17 1369.019
18 1373.444
19 1313.019
20 1171.265
21 915.527
22 584.717
23 222.931
24 127.411
Name: Area, dtype: float64
%% Cell type:markdown id: tags:
First we determine the mean of both measurementsand compare them
%% Cell type:code id: tags:
``` python
import numpy as np
mean_python = np.mean(measurement_python)
mean_imagej = np.mean(measurement_imagej)
print("Python: " + str(mean_python))
print("ImageJ: " + str(mean_imagej))
print("Difference: " + str(mean_python - mean_imagej))
```
%% Output
Python: 1010.858154296875
ImageJ: 1008.93564
Difference: 1.9225142968749651
%% Cell type:markdown id: tags:
Then we have a look at the histograms for the measurement
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
import numpy as np
def draw_histogram(data):
counts, bins = np.histogram(data, bins=10, range=(0,1400))
plt.hist(bins[:-1], bins, weights=counts)
#plt.axis([0, 10, 0, 4])
plt.show()
draw_histogram(measurement_python)
draw_histogram(measurement_imagej)
```
%% Output
%% Cell type:markdown id: tags:
We should also plot them against each other in a scatter plot to get a first impression on they relate to each other
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
# plot our data
plt.plot(measurement_python, measurement_imagej, "*")
# plot another line which corresponds to identidy
plt.plot([0, 1400], [0, 1400])
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
Now let's have a look at the Bland-Altman plot
%% Cell type:code id: tags:
``` python
# A function for drawing Bland-Altman plots
# source https://stackoverflow.com/questions/16399279/bland-altman-plot-in-python
import matplotlib.pyplot as plt
import numpy as np
def bland_altman_plot(data1, data2, *args, **kwargs):
data1 = np.asarray(data1)
data2 = np.asarray(data2)
mean = np.mean([data1, data2], axis=0)
diff = data1 - data2 # Difference between data1 and data2
md = np.mean(diff) # Mean of the difference
sd = np.std(diff, axis=0) # Standard deviation of the difference
plt.scatter(mean, diff, *args, **kwargs)
plt.axhline(md, color='gray', linestyle='--')
plt.axhline(md + 1.96*sd, color='gray', linestyle='--')
plt.axhline(md - 1.96*sd, color='gray', linestyle='--')
```
%% Cell type:code id: tags:
``` python
# draw a Bland-Altman plot
bland_altman_plot(measurement_python, measurement_imagej)
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
The cofindence interval
%% Cell type:code id: tags:
``` python
data1 = measurement_python
data2 = measurement_imagej
mean = np.mean([data1, data2], axis=0)
diff = data1 - data2 # Difference between data1 and data2
md = np.mean(diff) # Mean of the difference
sd = np.std(diff, axis=0) # Standard deviation of the difference
CI = [md - 2 * sd, md+2*sd]
print(CI)
```
%% Output
[0.3143669053595717, 3.530661688390393]
%% Cell type:code id: tags:
``` python
from statsmodels.stats.weightstats import ttost_ind
pval = ttost_ind(measurement_python, measurement_imagej, low=-10, upp=10)
pval
```
%% Output
(0.47051036223078757,
(0.1097795345763236, 0.4565209578658961, 48.0),
(-0.0743754713943518, 0.47051036223078757, 48.0))
%% Cell type:code id: tags:
``` python
# source: https://gist.github.com/josef-pkt/3900314
'''Test of Equivalence and Non-Inferiority
currently only TOST for paired sample
Application for example bioequivalence
Author: Josef Perktold
License: BSD-3
'''
import numpy as np
from scipy import stats
def tost_paired(y, x, low, upp, transform=None):
'''test of (non-)equivalence for paired sample
TOST: two one-sided t tests
null hypothesis x - y < low or x - y > upp
alternative hypothesis: low < x - y < upp
If the pvalue is smaller than a threshold, then we reject the hypothesis
that there is difference between the two samples larger than the one
given by low and upp.
Parameters
----------
y, x : array_like
two paired samples
low, upp : float
equivalence interval low < x - y < upp
transform : None or function
If None (default), then the data is not transformed. Given a function
sample data and thresholds are transformed. If transform is log the
the equivalence interval is in ratio: low < x / y < upp
Returns
-------
pvalue : float
pvalue of the non-equivalence test
t1, pv1 : tuple of floats
test statistic and pvalue for lower threshold test
t2, pv2 : tuple of floats
test statistic and pvalue for upper threshold test
Notes
-----
tested on only one example
uses stats.ttest_1samp which doesn't have a real one-sided option
'''
if transform:
y = transform(y)
x = transform(x)
low = transform(low)
upp = transform(upp)
t1, pv1 = stats.ttest_1samp(x - y, low)
t2, pv2 = stats.ttest_1samp(x - y, upp)
return max(pv1, pv2)/2., (t1, pv1 / 2.), (t2, pv2 / 2.)
```
%% Cell type:code id: tags:
``` python
tolerance = 0.1
tost_paired(measurement_python, measurement_imagej, -tolerance, tolerance)
```
%% Output
(3.0607159809231846e-11,
(-11.104032129950907, 3.0607159809231846e-11),
(-12.32257205021279, 3.601584280854005e-12))
%% Cell type:code id: tags:
``` python
from scipy import stats
def equivalence_ttest_rel(x1, x2, limit):
alpha = 0.05
# print(x1)
# print(x2)
# print(x1 - x2)
# less than test:
t1, p_value1 = stats.ttest_1samp(x1 - x2, limit)
if p_value1 / 2 < alpha and t1 < 0:
print("Reject null-hypothesis: Difference is above limit")
# greater than test:
t2, p_value2 = stats.ttest_1samp(x1 - x2, -limit)
if p_value2 / 2 < alpha and t2 > 0:
print("Reject null-hypothesis: Difference is below -limit")
# return the maximum p_value from both test
# representing the worst case scenario
return [
[t1, p_value1],
[t2, p_value2]
]
tolerance = 10
equivalence_ttest_rel(measurement_python, measurement_imagej, tolerance)
```
%% Output
Reject null-hypothesis: Difference is above limit
Reject null-hypothesis: Difference is below -limit
[[-49.213693923012265, 1.289351642345751e-25],
[72.64029810317597, 1.19906617345889e-29]]
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment