Değişkenlik hipotezi#
Değişkenlik hipotezi, erkeklerin fiziksel ve psikolojik özelliklerde genel olarak kadınlardan daha değişken olduğu kuramıdır. Hafif söylemek gerekirse tartışmalı bir iddiadır.
Konuyu hızlıca incelemek ve bootstrap yeniden örneklemesini göstermek için erkeklerin boylarının kadınlardan daha değişken olup olmadığına bakalım. 400.000’den fazla katılımcının kendi bildirdiği boy ve ağırlıkları içeren 2022 BRFSS verisini kullanacağız.
try:
import empiricaldist
except ImportError:
!pip install empiricaldist
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
np.random.seed(42)
plt.rcParams["figure.dpi"] = 75
plt.rcParams["figure.figsize"] = [6, 3.5]
Veriler#
Aşağıdaki işlev, 1 Ağustos 2024’te BRFSS sitesinden indirilmiş SAS dışa aktarma biçimindeki veriyi okur. Kullanacağımız sütunları seçip çok daha küçük alt kümeyi HDF dosyasına yazar.
def write_hdf():
"""SAS dışa aktarma dosyasını okuyup bir HDF dosyası yazar."""
brfss = pd.read_sas("LLCP2022.XPT.gz")
columns = ["_SEX", "HTM4", "WTKG3", "_LLCPWT"]
brfss[columns].to_hdf("LLCP2022.hdf5", key="brfss", complevel=6)
Tüm veri kümesi yerine yalnızca alt kümeyi indirebiliriz.
from os.path import basename, exists
def download(url):
filename = basename(url)
if not exists(filename):
from urllib.request import urlretrieve
local, _ = urlretrieve(url, filename)
print("İndirildi: " + local)
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/LLCP2022.hdf5")
İlk birkaç satır şöyledir.
brfss = pd.read_hdf("LLCP2022.hdf5", key="brfss")
brfss.describe()
| _SEX | HTM4 | WTKG3 | _LLCPWT | |
|---|---|---|---|---|
| count | 445132.000000 | 416480.000000 | 403054.000000 | 445132.000000 |
| mean | 1.529942 | 170.269057 | 8307.447039 | 594.856344 |
| std | 0.499103 | 10.717750 | 2144.817270 | 1134.837415 |
| min | 1.000000 | 91.000000 | 2268.000000 | 0.020464 |
| 25% | 1.000000 | 163.000000 | 6804.000000 | 115.885991 |
| 50% | 2.000000 | 170.000000 | 8074.000000 | 274.632388 |
| 75% | 2.000000 | 178.000000 | 9525.000000 | 627.913694 |
| max | 2.000000 | 241.000000 | 29257.000000 | 54390.520926 |
Veri kümesi yaklaşık 200.000 erkek ve 220.000 kadının bildirdiği boyları içerir. Kaydedilen en kısa boylara sahip erkek ve kadın sayıları şöyledir.
xtab = pd.crosstab(brfss["HTM4"], brfss["_SEX"])
xtab.columns = ["Erkek", "Kadın"]
xtab.head()
| Erkek | Kadın | |
|---|---|---|
| HTM4 | ||
| 91.0 | 7 | 17 |
| 92.0 | 0 | 1 |
| 95.0 | 1 | 0 |
| 97.0 | 1 | 3 |
| 99.0 | 1 | 0 |
En uzun kaydedilen boylar da şöyledir.
xtab.tail()
| Erkek | Kadın | |
|---|---|---|
| HTM4 | ||
| 226.0 | 9 | 2 |
| 229.0 | 4 | 1 |
| 234.0 | 4 | 0 |
| 236.0 | 1 | 0 |
| 241.0 | 4 | 1 |
Kod kitabına göre 91 cm altındaki boylar 91, 241 cm üzerindekiler 241 cm yapılır. Yine de kalan bazı değerler muhtemelen hatalıdır. Dünyanın yaşayan en uzun kadını 234 cm olduğundan veri kümesindeki 241 cm büyük olasılıkla yanlış bildirilmiş ya da kaydedilmiştir.
Bu DataFrame’den bootstrap örneklemi çekmek için aşağıdaki işlevi kullanacağız.
def resample(df):
"""Bir bootstrap örneklemi çeker.
df: DataFrame
döndürür: DataFrame
"""
n = len(df)
return df.sample(n, replace=True, weights="_LLCPWT")
Sınamada kullanacağımız tek bir örneklem şöyledir.
sample = resample(brfss)
Ortalamalar arasındaki fark#
Erkeklerin daha değişken olup olmadığını sınamadan önce ortalama olarak daha uzun olduklarını doğrulayalım. Ortalama farkının örnekleme dağılımını bootstrap ile tahmin edeceğiz. Aşağıdaki yardımcı işlev bir dizinin ilk iki öğesi arasındaki farkı hesaplar.
def diff(seq):
"""Bir dizinin ilk iki öğesi arasındaki fark.
seq: dizi
döndürür: float
"""
return np.diff(seq)[0]
Bunu bir örneklemdeki ortalama farkını, yani “test istatistiğini” hesaplamak için kullanacağız.
def diff_means(sample):
"""Ortalama boy farkı (erkek eksi kadın).
sample: _SEX ve HTM4 sütunlarını içeren DataFrame
döndürür: float
"""
grouped_heights = sample.groupby("_SEX")["HTM4"]
return -diff(grouped_heights.mean())
Fark, pozitif değerin erkeklerin daha uzun olduğunu göstereceği biçimde tanımlıdır. Bootstrap örneklemiyle örnek şöyledir.
diff_means(sample)
np.float64(14.455796580284954)
Erkekler ortalama yaklaşık 14 cm daha uzundur.
Aşağıdaki işlev bir DataFrame ve test istatistiği hesaplayan işlev alır. Çok kez yeniden örnekler, her örneklemin istatistiğini hesaplar ve liste döndürür.
def sampling_dist(df, test_stat, iters=201):
"""Bootstrap örneklemleri üretip test istatistiğini hesaplar.
df: DataFrame
test_stat: DataFrame alan işlev
iters: örneklem sayısı
döndürür: float listesi
"""
return [test_stat(resample(df)) for i in range(iters)]
Şöyle çağırırız. Sonuç ortalama farklarının örnekleme dağılımından bir örneklemdir.
t1 = sampling_dist(brfss, diff_means)
np.mean(t1)
np.float64(14.433755828000592)
Aşağıdaki işlev KDE ile örnekleme dağılımının yoğunluğunu tahmin eder, güven aralığını hesaplar ve dağılımı, ortalamayı ve %95 GA’yı gösteren grafik üretir.
from scipy.stats import gaussian_kde
from empiricaldist import Pmf
def plot_sampling_distribution(data, format_str=".2f", **options):
"""Örnekleme dağılımını, ortalamayı ve güven aralığını çizer.
data: örnekleme dağılımından örneklem
format_str: değerleri biçimlendirmek için kullanılır
options: plt.plot'a iletilen seçenekler
"""
# ortalamayı ve güven aralığını hesapla
m, s = np.mean(data), np.std(data)
qs = np.linspace(m - 4 * s, m + 4 * s, 501)
ps = gaussian_kde(data)(qs)
pmf = Pmf(ps, qs)
pmf.normalize()
low, high = pmf.make_cdf().inverse([0.025, 0.975])
# yoğunluğu çiz
pmf.plot(**options)
inside = pmf[(pmf.qs >= low) & (pmf.qs <= high)]
plt.fill_between(inside.index, inside, color="gray", alpha=0.2)
y = pmf.max()
ax = plt.gca()
# metinleri ve çizgileri ekle
plt.text(m, 0.4 * y, "ortalama", fontsize=14, ha="center")
plt.text(m, 0.05 * y, f"%95 GA", fontsize=14, ha="center")
plt.plot([m, m], [0, y], ls=":", color="gray")
plt.plot([low, high], [0, 0], ls="-", color="gray")
# adjust xticks
ticks = [low, m, high]
labels = [format(tick, format_str) for tick in ticks]
ax.tick_params(axis="x", labelsize=14)
plt.xticks(ticks, labels)
plt.yticks([])
# remove spines
for spine in ax.spines.values():
spine.set_visible(False)
Ortalama farkının örnekleme dağılımı şöyledir.
plot_sampling_distribution(t1)
plt.title("Ortalama farkı (cm)", fontsize=14)
plt.tight_layout()
plt.savefig("variability1.png", dpi=300)
Ortalama farkı yaklaşık 14,4 cm’dir. Veri toplama sürecini çok kez simüle edince fark biraz değişir, ancak sıfıra yaklaşmaz.
Aşağıdaki işlev bu dağılımdan bir değerin 0’ı aşma olasılığını, yani erkeklerin ortalama daha uzun olduğu hipotezinin p-değerini tahmin eder.
from scipy.stats import norm
def p_value(t):
m, s = np.mean(t), np.std(t)
if m > 0:
return norm.cdf(0, m, s)
else:
return norm.sf(0, m, s)
p-değeri kayan nokta aritmetiğiyle hesaplanamayacak kadar sıfıra yakındır; gerçekten ihmal edilebilir.
p_value(t1)
np.float64(0.0)
Aynı evrenden bu büyüklükte başka bir örneklem alırsak kadınların ortalama olarak daha uzun çıkması uygulamada olanaksızdır.
Örneklemde gördüğümüz farkın şans eseri oluşmasının çok düşük olasılıklı olduğu sonucuna da makul biçimde varabiliriz. Hipotez testlerini çok biçimsel yorumlayan bazı kişiler bunu sevmez, ancak görüşümün arkasındayım.
Standart sapmalar arasındaki fark#
Şimdi asıl soruya gelelim. Erkekler daha değişkense boylarının standart sapmasının daha yüksek olmasını bekleriz. Bir örneklem alıp standart sapma farkını hesaplayan işlev yeterlidir.
def diff_stds(sample):
"""Standart sapma farkı (erkek eksi kadın).
sample: `_SEX` ve `HTM4` sütunlarını içeren DataFrame
döndürür: float
"""
grouped_heights = sample.groupby("_SEX")["HTM4"]
return -diff(grouped_heights.std())
Tek yeniden örneklemeye göre erkek boylarının standart sapması daha yüksek görünüyor.
diff_stds(sample)
np.float64(0.7755433526023632)
sampling_dist işlevini farklı test istatistiğiyle yeniden kullanabiliriz.
t2 = sampling_dist(brfss, diff_stds)
np.mean(t2)
np.float64(0.8125943377103155)
Ortalaması ve %95 GA’sıyla örnekleme dağılımı şöyledir.
plot_sampling_distribution(t2)
plt.title("Standart sapma farkı (cm)", fontsize=14)
plt.tight_layout()
plt.savefig("variability2.png", dpi=300)
p_value(t2)
np.float64(3.3148844524687797e-169)
p-değeri ancak hesaplanabilecek kadar büyüktür; yine de gözlenen farkın şanstan kaynaklanmasının çok düşük olasılıklı olduğu sonucuna varacak kadar küçüktür.
Ancak standart sapmaları karşılaştırmak erkeklerin daha değişken olup olmadığını değerlendirmek için anlamlı olmayabilir. Erkekler ortalama daha uzundur; mutlak olarak daha fazla değişmeleri göreli olarak da daha fazla değiştikleri anlamına gelmez. Ortalamaya göre değişim daha ilgiliyse değişim katsayısını ele alabiliriz.
Değişim katsayıları arasındaki fark#
Değişim katsayısı standart sapmanın ortalamaya oranıdır; değişkenliği göreli olarak nicelleştirir. Aşağıdaki işlev erkeklerle kadınların değişim katsayısı farkını hesaplar.
def diff_cvs(sample):
"""DK farkı (erkek eksi kadın).
sample: `_SEX` ve `HTM4` sütunlarını içeren DataFrame
döndürür: float
"""
grouped_heights = sample.groupby("_SEX")["HTM4"]
means = grouped_heights.mean()
stds = grouped_heights.std()
return -diff(stds / means)
Farkların örnekleme dağılımından bir örneklem şöyledir.
t3 = sampling_dist(brfss, diff_cvs)
np.mean(t3)
np.float64(0.000673958474897113)
Örnekleme dağılımı ve %95 GA şöyledir.
plot_sampling_distribution(t3, format_str="0.4f")
plt.title("Değişim katsayısı farkı", fontsize=14)
plt.tight_layout()
plt.savefig("variability3.png", dpi=300)
GA sıfırı içermez; fark istatistiksel olarak anlamlıdır. p-değeri de küçüktür.
p_value(t3)
np.float64(4.91233303730651e-05)
Ancak uygulamada fark çok küçüktür. Değişim katsayıları yaklaşık 0,048’dir.
grouped_heights = sample.groupby("_SEX")["HTM4"]
means = grouped_heights.mean()
stds = grouped_heights.std()
CVs = stds / means
CVs
_SEX
1.0 0.048558
2.0 0.048107
Name: HTM4, dtype: float64
DK farkını katsayıların yüzdesi olarak yazarsak yalnızca yaklaşık %1,4’tür.
np.mean(t3) / CVs * 100
_SEX
1.0 1.387945
2.0 1.400971
Name: HTM4, dtype: float64
Bu ölçüye göre erkek boyları biraz daha değişkendir; farkın uygulamada sonucu olması beklenmez.
Ayrıca standart sapma ile değişim katsayısı aykırı değerlerden etkilenir ve veri kümesindeki uçların muhtemelen hatalı olduğunu gördük. Görünürdeki fark veri hatalarının sonucu olabilir.
Daha dirençli istatistikle yeniden deneyelim.
Dirençli DK farkı#
Değişim katsayısına alternatif olarak daha dirençli değişkenlik ölçüsü kullanabiliriz. Seçeneklerden biri medyan mutlak sapmanın (MAD) medyana oranıdır. Aşağıdaki işlev bunu hesaplar.
def robust_cv(series):
"""Medyan mutlak sapma.
series: sayı serisi
döndürür: float
"""
m = series.median()
deviations = series - m
return deviations.abs().median() / m
Aşağıdaki işlev dirençli DK farkını hesaplar. İlk satır verilere titreşim ekler; aksi durumda her yeniden örneklemede aynı sonuçları alırdık.
def diff_robust_cv(sample):
"""Dirençli DK farkı (erkek eksi kadın).
sample: `_SEX` ve `HTM4` sütunlarını içeren DataFrame
döndürür: float
"""
sample["HTM4"] += np.random.normal(0, 1, size=len(sample))
grouped_heights = sample.groupby("_SEX")["HTM4"]
rcvs = grouped_heights.apply(robust_cv)
return -diff(rcvs)
Örneklemin dirençli DK farkı negatiftir; kadın boylarının daha değişken olduğunu düşündürür.
diff_robust_cv(sample)
np.float64(-0.0025971928221166715)
Örnekleme dağılımı şöyledir.
t4 = sampling_dist(brfss, diff_robust_cv)
np.mean(t4)
np.float64(-0.002502441445045549)
plot_sampling_distribution(t4, format_str="0.4f")
plt.title("Dirençli DK farkı", fontsize=14)
plt.tight_layout()
plt.savefig("variability4.png", dpi=300)
p-değeri de şöyledir.
p_value(t4)
np.float64(6.93074075784399e-122)
Dirençli DK’ye göre kadınlar daha değişkendir; fakat fark yine uygulamada önemsizdir.
Tartışma#
Bu örnek bootstrap yeniden örneklemesinin gücünü gösterir: birden çok test istatistiğini incelemek ve hipoteze en uygun olanı değerlendirmek kolaydır.
Boyların standart sapmasını karşılaştırırsak erkekler daha değişkendir; ama daha uzun olduklarından bu şaşırtıcı değildir. Büyük şeyler üreten süreçler çoğu zaman daha değişken şeyler de üretir.
Farklı büyüklükleri karşılaştırmak için değişim katsayısı daha iyi olabilir. Buna göre erkekler yine biraz daha değişkendir; ancak fark uygulamada önemsizdir.
Sonuç veri hatalarının ürünü olabilir. DK’nin dirençli alternatifi MAD/medyan kullanılırsa kadınlar daha değişkendir; fakat yine önemli bir fark değildir.
Sonuç olarak insan boyu dağılımı değişkenlik hipotezine güçlü destek sağlamaz.
ÇDA / medyan farkı#
Bir dirençli değişkenlik ölçüsü daha.
def diff_relative_iqr(sample):
"""ÇAA / medyan farkı (erkek eksi kadın).
sample: `_SEX` ve `HTM4` sütunlarını içeren DataFrame
döndürür: float
"""
sample["HTM4"] += np.random.normal(0, 1, size=len(sample))
grouped_heights = sample.groupby("_SEX")["HTM4"]
medians = grouped_heights.median()
# sonuçlar kullandığımız kantillere bağlıdır
iqrs = grouped_heights.quantile(0.75) - grouped_heights.quantile(0.25)
#iqrs = grouped_heights.quantile(0.95) - grouped_heights.quantile(0.05)
return -diff(iqrs / medians)
diff_relative_iqr(sample)
np.float64(-0.0031180661772347157)
t5 = sampling_dist(brfss, diff_relative_iqr)
np.mean(t5)
np.float64(-0.004621824998497872)
plot_sampling_distribution(t5, format_str="0.4f")
plt.title("ÇAA / medyan farkı", fontsize=14)
plt.tight_layout()
plt.savefig("variability5.png", dpi=300)
p_value(t5)
np.float64(6.773833573551112e-98)