Think Stats kitabının üçüncü baskısını Bookshop.org ve Amazon üzerinden edinebilirsiniz (bunlar satış ortaklığı bağlantılarıdır). Ücretsiz çevrimiçi sürümü beğeniyorsanız bana bir kahve ısmarlamayı düşünebilirsiniz.
13. Sağkalım analizi#
Sağkalım analizi, bir şeyin ne kadar sürdüğünü inceleme yöntemidir. Sıklıkla insan ömrünü incelemek için kullanılır; ancak mekanik ve elektronik parçaların dayanma sürelerine, herhangi bir olay gerçekleşene kadar geçen zamana, hatta uzaydaki aralıklara da uygulanabilir.
Ampullerin ömrüyle ilgili basit bir örnekle başlayıp ilk evlenme yaşı ve bunun ABD’de son 50 yıldaki değişimini ele alacağız.
Bu notebook’u Colab’da çalıştırmak için buraya tıklayın.
13.1. Sağkalım fonksiyonları#
Sağkalım analizinin temel kavramı, incelenen evrende belirli bir süreden uzun dayananların oranını gösteren sağkalım fonksiyonudur. İlk örnekte ampul ömürlerinin sağkalım fonksiyonunu hesaplayacağız.
2007’deki bir deneyin verilerini kullanacağız. Araştırmacılar 50 yeni ampul takıp sürekli açık bıraktı. Her 12 saatte ampulleri denetleyip sönenlerin ömrünü kaydettiler ve 50 ampul de sönene kadar deneyi sürdürdüler. İndirme yönergeleri bu bölümün notebook’unda yer alır.
Aşağıdaki hücre, burada belgelenen verileri indirir.
Veri kaynağı:
V.J. Menon ve D.C. Agrawal, “Renewal Rate of Filament Lamps: Theory and Experiment”. Journal of Failure Analysis and Prevention, Aralık 2007, s. 421, Tablo 2. DOI: 10.1007/s11668-007-9074-9
download(
"https://gist.github.com/epogrebnyak/7933e16c0ad215742c4c104be4fbdeb1/raw/c932bc5b6aa6317770c4cbf43eb591511fec08f9/lamps.csv"
)
Verileri şöyle okuyabiliriz.
df = pd.read_csv("lamps.csv", index_col=0)
df.tail()
| h | f | K | |
|---|---|---|---|
| i | |||
| 28 | 1812 | 1 | 4 |
| 29 | 1836 | 1 | 3 |
| 30 | 1860 | 1 | 2 |
| 31 | 1980 | 1 | 1 |
| 32 | 2568 | 1 | 0 |
h sütunu saat cinsinden ömürleri, f sütunu her h değerinde sönen ampul sayısını içerir. Ömür dağılımını göstermek için değerleri bir Pmf nesnesine koyup normalize edeceğiz.
from empiricaldist import Pmf
pmf_bulblife = Pmf(df["f"].values, index=df["h"])
pmf_bulblife.normalize()
np.int64(50)
Her h değerinde veya öncesinde sönenlerin oranını gösteren CDF’yi make_cdf ile hesaplayabiliriz. Örneğin ampullerin %78’i 1656 saat veya öncesinde sönmüştür.
cdf_bulblife = pmf_bulblife.make_cdf()
cdf_bulblife[1656]
np.float64(0.7800000000000002)
Sağkalım fonksiyonu her h değerinden sonra sönenlerin oranıdır ve CDF’nin tümleyenidir. Şöyle hesaplayabiliriz.
complementary_cdf = 1 - cdf_bulblife
complementary_cdf[1656]
np.float64(0.21999999999999975)
Ampullerin %22’si 1656 saatten sonra sönmüştür.
empiricaldist, sağkalım fonksiyonunu temsil eden Surv nesnesini ve bunu oluşturan make_surv yöntemini sağlar.
surv_bulblife = cdf_bulblife.make_surv()
surv_bulblife[1656]
np.float64(0.21999999999999997)
CDF ile sağkalım fonksiyonunu çizersek tüm h değerlerinde toplamlarının 1 olduğunu, yani birbirlerinin tümleyeni olduklarını görürüz.
cdf_bulblife.plot(ls="--", label="CDF")
surv_bulblife.plot(label="Sağkalım")
decorate(xlabel="Ampul ömrü (saat)", ylabel="Olasılık")
Bu anlamda CDF ile sağkalım fonksiyonu eşdeğerdir; biri verilirse diğeri hesaplanabilir. Ancak sağkalım analizinde sağkalım eğrileri daha yaygındır. Sağkalım eğrisini hesaplamak bir sonraki önemli kavrama, tehlike fonksiyonuna doğru bir adımdır.
13.2. Tehlike fonksiyonu#
Ampul veri kümesindeki her h, h saatinde biten 12 saatlik aralığı, yani “h aralığını” temsil eder. Ampulün h aralığına kadar dayandığını bildiğimizi ve bu aralıkta sönme olasılığını istediğimizi varsayın. Sağkalım fonksiyonu h sonrasına kalanların, PMF ise h içinde sönenlerin oranını gösterir. Toplamları, h aralığında sönebilecek, yani “risk altındaki” ampullerin oranıdır. Örneğin 1656 aralığında ampullerin %26’sı risk altındadır.
at_risk = pmf_bulblife + surv_bulblife
at_risk[1656]
np.float64(0.25999999999999995)
Tüm ampullerin %4’ü 1656 aralığında sönmüştür.
pmf_bulblife[1656]
np.float64(0.04)
Tehlike, pmf_bulblife ile at_risk oranıdır.
hazard = pmf_bulblife / at_risk
hazard[1656]
np.float64(0.15384615384615388)
1656 aralığına kadar dayanan ampullerin yaklaşık %15’i bu aralıkta sönmüştür.
Tehlike fonksiyonunu kendimiz hesaplamak yerine empiricaldist içindeki Hazard nesnesini ve make_hazard yöntemini kullanabiliriz.
hazard_bulblife = surv_bulblife.make_hazard()
hazard_bulblife[1656]
np.float64(0.15384615384615397)
Ampullerin tehlike fonksiyonu şöyledir.
hazard_bulblife.plot()
decorate(xlabel="Ampul ömrü (saat)", ylabel="Tehlike")
Tehlikenin bazı yerlerde daha yüksek olduğunu görebiliriz; ancak özellikle az verili aralıklarda bu görselleştirme yanıltıcı olabilir. Daha iyi bir seçenek, tehlikelerin kümülatif toplamı olan kümülatif tehlike fonksiyonunu çizmektir.
cumulative_hazard = hazard_bulblife.cumsum()
cumulative_hazard.plot()
decorate(xlabel="Ampul ömrü (saat)", ylabel="Kümülatif tehlike")
Sönme olasılığının yüksek olduğu yerde kümülatif tehlike eğimi dik, düşük olduğu yerde yatıktır. Tehlikenin 1500–2000 saat arasında en yüksek olduğunu görürüz. Sonrasında azalır; ancak bu sonuç alışılmadık ölçüde uzun ömürlü tek bir ampule dayandığından başka verilerde farklı olabilir.
Sağkalım ve tehlike fonksiyonlarının genel fikrini daha kapsamlı bir veri kümesine uygulayalım.
13.3. Evlilik verileri#
Birçok ülkede insanlar geçmişe göre daha geç evleniyor ve daha fazla kişi evlenmeden kalıyor. ABD’deki eğilimleri incelemek için sağkalım araçlarını ve Ulusal Aile Büyümesi Araştırması (NSFG) verilerini kullanacağız.
Önceki bölümlerdeki NSFG gebelik dosyasında bildirilen her gebelik için bir satır vardır. Bu bölümde katılımcıların kendileri hakkındaki bilgileri içeren katılımcı dosyasıyla çalışacağız.
1982–2019 arasında yapılan dokuz araştırmanın yanıtlarını birleştirip evlilikle ilgili verileri seçtim. Alt kümeyi indirme yönergeleri notebook’tadır.
Aşağıdaki hücre, 1982–2019 arasındaki çeşitli NSFG döngülerini birleştirerek oluşturduğum CSV dosyasını indirir.
Hazırlama ayrıntıları bu notebook’tadır.
filename = "marriage_nsfg_female.csv.gz"
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/" + filename)
Verileri şöyle okuyabiliriz.
resp = pd.read_csv("marriage_nsfg_female.csv.gz")
resp.shape
(70183, 34)
Alt kümede 70.000’den fazla katılımcı için birer satır ve yaş/evlilikle ilgili şu değişkenler vardır:
cmbirth: Herkes için bilinen doğum tarihi.cmintvw: Herkes için bilinen görüşme tarihi.cmmarrhx: Uygunsa ve biliniyorsa ilk evlilik tarihi.evrmarry: Görüşmeden önce evlendiyse 1, aksi durumda 0.
İlk üçü Aralık 1899’dan beri geçen tam ay sayısı olan “yüzyıl ayları”yla kodlanmıştır; 1, Ocak 1900’dür.
Kuşak değişimini incelemek için katılımcıları doğdukları on yıla göre gruplayacağız. Aşağıdaki işlev cmbirth değerinden doğum on yılını hesaplar; // ile 10’a bölüp aşağı yuvarlar.
month0 = pd.to_datetime("1899-12-31")
def decade_of_birth(cmbirth):
date = month0 + pd.DateOffset(months=cmbirth)
return date.year // 10 * 10
İşlevi apply ile kullanıp her katılımcının doğum on yılını hesaplayarak cohort sütununa atayabiliriz. Bu bağlamda kohort, doğduğu on yıl gibi ortak özelliği bulunan ve analizde grup olarak ele alınan kişilerdir.
value_counts her kohorttaki kişi sayısını gösterir.
from thinkstats import value_counts
resp["cohort"] = resp["cmbirth"].apply(decade_of_birth)
value_counts(resp["cohort"])
cohort
1930 325
1940 3608
1950 10631
1960 14953
1970 16438
1980 14271
1990 8552
2000 1405
Name: count, dtype: int64
Veri kümesinde 1950’lerden 1980’lere her on yılda doğan 10.000’den fazla, daha erken ve geç on yıllarda daha az katılımcı vardır.
Şimdi uygunsa ilk evlilik yaşını ve görüşme yaşını hesaplayacağız.
resp["agemarr"] = (resp["cmmarrhx"] - resp["cmbirth"]) / 12
resp["age"] = (resp["cmintvw"] - resp["cmbirth"]) / 12
Başlangıç için bir DataFrame ile kohort listesini alıp her kohortu sağkalım fonksiyonunu temsil eden Surv nesnesine eşleyen sözlük döndüren işlevi kullanacağız. Her kohort için ilk evlilik yaşlarını seçip Surv.from_seq çağırır. dropna=False, evlenmemiş kişileri de sağkalım fonksiyonuna katar.
from empiricaldist import Surv
def make_survival_map(resp, cohorts):
surv_map = {}
grouped = resp.groupby("cohort")
for cohort in cohorts:
group = grouped.get_group(cohort)
surv_map[cohort] = Surv.from_seq(group["agemarr"], dropna=False)
return surv_map
İşlevi şöyle kullanırız.
cohorts = [1980, 1960, 1940]
surv_map = make_survival_map(resp, cohorts)
1940’lar, 1960’lar ve 1980’lerde doğanların sonuçları şöyledir.
import matplotlib.pyplot as plt
import cycler
# Varsayılan renk döngüsünü çıkar
default_colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
# İstenen çizgi biçemlerini tanımla
linestyles = ["--", "-", "-.", "--", "-", "-.", "--", "-", "-.", "--"]
# Biçem ve renklerde doğru sırayla dolaşılmasını sağla
custom_cycler = cycler.cycler(color=default_colors) + cycler.cycler(
linestyle=linestyles
)
# Yeni döngüyü uygula
plt.rc("axes", prop_cycle=custom_cycler)
for cohort, surv in surv_map.items():
surv.plot(label=f"{cohort}s")
ylim = [-0.05, 1.05]
decorate(xlabel="Yaş", ylabel="Hiç evlenmeme olasılığı", ylim=ylim)
Sonuçları olduğu gibi alırsak eski kuşakların daha genç evlendiğini ve daha çoğunun sonunda evlendiğini gösterir. Ancak henüz yorumlamamalıyız; doğru değiller. İki sorun vardır:
Bölüm 1’de açıklandığı gibi NSFG bazı grupları bilinçli olarak fazla örnekleyen tabakalı örnekleme kullanır.
Sağkalım fonksiyonunu bu biçimde hesaplamak henüz evlenmemiş kişileri doğru hesaba katmaz.
İlk sorun için ağırlıklı bootstrap, ikincisi için Kaplan–Meier tahmini kullanacağız. Yeniden örneklemeyle başlayalım.
13.4. Ağırlıklı bootstrap#
NSFG’deki finalwgt sütunu, her katılımcının evrende temsil ettiği kişi sayısını belirten örnekleme ağırlığını içerir. Tabakalı örneklemeyi düzeltmek için yeniden örnekleme sırasında bu ağırlıkları kullanabiliriz. Aşağıdaki işlev bir DataFrame ile ağırlık sütununun adını alır; ağırlıkları hesaba katarak satırları yeniden örnekleyip yeni DataFrame döndürür.
def resample_rows_weighted(df, column="finalwgt"):
n = len(df)
weights = df[column]
return df.sample(n, weights=weights, replace=True)
Veri kümesi araştırmanın “döngü” denen çeşitli yinelemelerinden katılımcılar içerir. Yeniden örneklemek için kişileri döngüye göre gruplamalı, her grubu yeniden örnekleyip grupları birleştirmeliyiz. Aşağıdaki işlev bunu yapar.
def resample_cycles(resp):
grouped = resp.groupby("cycle")
samples = [resample_rows_weighted(group) for _, group in grouped]
return pd.concat(samples)
Başlangıç olarak verileri bir kez yeniden örnekleyeceğiz.
# Her çalıştırmada aynı sonuçları almak için rastgele sayı üretecini tohumla
np.random.seed(1)
sample = resample_cycles(resp)
Rastgele örneklemeden kaynaklanan değişimi görmek için daha sonra veriyi birkaç kez yeniden örnekleyeceğiz.
Aşağıdaki şekil yeniden örneklemeli sonuçları, noktalı çizgilerle gösterilen yeniden örneklemesiz sonuçlarla karşılaştırır.
for label, surv in surv_map.items():
surv.plot(ls=":", color="gray", alpha=0.6)
survs_resampled = make_survival_map(sample, cohorts)
for label, surv in survs_resampled.items():
surv.plot(label=label)
decorate(xlabel="Yaş", ylabel="Hiç evlenmeme olasılığı", ylim=ylim)
plt.rc("axes", prop_cycle=plt.rcParamsDefault["axes.prop_cycle"])
Yeniden örneklemeli ve örneklemesiz sonuçlar arasındaki fark büyüktür; doğru sonuç için tabakalı örneklemeyi düzeltmemiz gerektiğini gösterir.
Şimdi ikinci soruna, eksik veriye geçelim.
13.5. Tehlike fonksiyonlarının tahmini#
Ampul örneğinde 50 ampulün de ömrünü bildiğimiz için sağkalım fonksiyonunu doğrudan, tehlike fonksiyonunu da ondan hesaplayabiliriz.
Evlilik örneğinde görüşmeden önce evlenmiş kişilerin ilk evlilik yaşını biliyoruz. Hiç evlenmemiş olanların gelecekte hangi yaşta evleneceğini ya da evlenip evlenmeyeceğini bilmiyoruz.
Bu tür eksik verilere sansürlü denir. Terim garip görünebilir; bilgi kasıtlı değil, geleceği bilmediğimiz için gizlidir.
Yine de kısmi bilgi vardır: araştırma sırasında evli olmayan birinin evlenme yaşı, evlenecekse, şimdiki yaşından büyük olmalıdır.
Bu bilgiyle tehlike fonksiyonunu, ondan da sağkalım fonksiyonunu tahmin edebiliriz. Buna Kaplan–Meier tahmini denir.
Göstermek için yeniden örneklenmiş veriden bir kohort seçeceğim.
resp60 = sample.query("cohort == 1960")
Araştırma sırasında evli olanların ilk evlilik yaşlarını seçeceğiz. Bunlardan 9921 tane vardır; “tamamlanmış” olgular diyeceğiz.
complete = resp60.query("evrmarry == 1")["agemarr"]
complete.count()
np.int64(9921)
Evlenmemiş katılımcıların görüşme yaşını seçeceğiz. Bunlardan 5468 tane vardır; “süren” olgular diyeceğiz.
ongoing = resp60.query("evrmarry == 0")["age"]
ongoing.count()
np.int64(5468)
Tehlike fonksiyonunu tahmin etmek için her yaşa kadar evlenmemiş herkesi içeren “risk altındaki” toplam olgu sayısını hesaplayacağız. Her yaştaki tamamlanmış ve süren olguları sayan FreqTab nesneleri kullanışlıdır.
from empiricaldist import FreqTab
ft_complete = FreqTab.from_seq(complete)
ft_ongoing = FreqTab.from_seq(ongoing)
Örneğin 58 katılımcı ilk kez 25 yaşında evlendiğini bildirmiştir.
ft_complete[25]
np.int64(58)
Başka 5 katılımcı 25 yaşında görüşülmüş ve hiç evlenmediğini bildirmiştir.
ft_ongoing[25]
np.int64(5)
Bu FreqTab nesnelerinden her yaşı aşan tamamlanmış ve süren olguların sayısını içeren normalize edilmemiş Surv nesneleri hesaplayabiliriz.
surv_complete = ft_complete.make_surv()
surv_ongoing = ft_ongoing.make_surv()
Örneğin 2848 kişi 25 yaşından sonra evlendiğini bildirmiştir.
surv_complete[25]
np.int64(2848)
25 yaşından sonra görüşülen ve hiç evlenmemiş 2273 kişi vardır.
surv_ongoing[25]
np.int64(2273)
Hesapladığımız dört sayının toplamı 25 yaşında evlenebilecek, yani risk altındaki katılımcıların sayısıdır. “Risk altında” terimi, hastalık ya da ölüm riskini anlattığı tıptaki sağkalım analizinden kalmadır; genellikle olumlu bir dönüm noktası kabul edilen evlilik bağlamında uyumsuz görünebilir. Şöyle hesaplarız.
at_risk = ft_complete[25] + ft_ongoing[25] + surv_complete[25] + surv_ongoing[25]
at_risk
np.int64(5184)
Bunlardan 25 yaşında gerçekten evlenenlerin sayısı ft_complete[25]’tir. Tehlikeyi şöyle hesaplayabiliriz.
hazard = ft_complete[25] / at_risk
hazard
np.float64(0.011188271604938271)
Tek bir yaşta tehlike fonksiyonunu böyle hesaplarız. Tüm yaşlar için Index.union yöntemiyle ft_complete ile ft_ongoing içindeki bütün yaşları sıralı bir Pandas Index içinde birleştireceğiz.
ts = pd.Index.union(ft_complete.index, ft_ongoing.index)
Şimdi ts içindeki yaşları dört nesnede arayarak her yaştaki risk altındaki kişi sayısını hesaplayabiliriz.
at_risk = ft_complete(ts) + ft_ongoing(ts) + surv_complete(ts) + surv_ongoing(ts)
Son olarak her yaşın tehlikesini hesaplayıp sonuçları Hazard nesnesine koyabiliriz.
from empiricaldist import Hazard
hs = ft_complete(ts) / at_risk
hazard = Hazard(hs, ts)
Kümülatif tehlike fonksiyonu şöyledir.
hazard.cumsum().plot()
decorate(xlabel="Yaş", ylabel="Kümülatif tehlike")
20–30 yaş arasında en diktir; evli olmayan kişinin evlenme “riskinin” bu yaşlarda en yüksek olduğunu gösterir. Sonrasında yatıklaşır, yani tehlike giderek azalır.
13.6. Sağkalım fonksiyonlarının tahmini#
Sağkalım fonksiyonundan tehlikeyi hesaplamayı biliyoruz; şimdi ters yöne gidelim.
Tehlike fonksiyonu, daha önce evlenmemiş olma koşuluyla her yaşta evlenme olasılığını gösterir. Tümleyeni her yaşta evlenmeden kalma olasılığıdır.
Belirli t yaşını “aşmak” için t dâhil her yaşta evlenmeden kalmak gerekir. Bunun olasılığı tümleyen tehlike fonksiyonunun çarpımıdır; şöyle hesaplarız.
ps = (1 - hazard).cumprod()
Hazard nesnesinin bu hesabı yapan make_surv yöntemi vardır.
surv = hazard.make_surv()
Sonuç, tabakalı örneklemeyi düzeltip sansürlü veriyi işlemeyen önceki sonuçla (noktalı çizgi) karşılaştırıldığında şöyledir.
survs_resampled[1960].plot(ls=":", color="gray", label="yeniden örneklenmiş")
surv.plot(label="Kaplan-Meier")
decorate(xlabel="Yaş", ylabel="Hiç evlenmeme olasılığı", ylim=ylim)
Sansürlü veriyi doğru işlemenin ne kadar önemli olduğunu görebiliriz.
1986 tarihli ünlü bir dergi yazısı buna benzer bir sağkalım fonksiyonuna dayanıyordu. Newsweek, 40 yaşındaki evlenmemiş bir kadının evlenmektense “bir terörist tarafından öldürülme olasılığının daha yüksek” olduğunu bildirdi. İddia yaygınlaşıp popüler kültüre girdi; fakat hem hatalı analize dayandığı için o zaman yanlıştı hem de sürmekte olan kültürel değişimler nedeniyle daha da yanlış çıktı. Newsweek 2006’da hatasını kabul eden başka bir yazı yayımladı.
Bu yazı, dayandığı istatistikler ve tepkiler hakkında okumayı öneririm. İstatistiksel analizi özenle yapma, sonuçları uygun kuşkuculukla yorumlama ve kamuya doğru, dürüst biçimde sunma etik yükümlülüğünü anımsatmalıdır.
Snopes’un konuyla ilgili iyi bir yazısı vardır.
Aşağıdaki işlev Kaplan–Meier tahmininin adımlarını kapsar. Tamamlanmış ve süren olguların sağkalım sürelerini alıp Hazard nesnesi döndürür.
def estimate_hazard(complete, ongoing):
"""Kaplan–Meier tahmini."""
ft_complete = FreqTab.from_seq(complete)
ft_ongoing = FreqTab.from_seq(ongoing)
surv_complete = ft_complete.make_surv()
surv_ongoing = ft_ongoing.make_surv()
ts = pd.Index.union(ft_complete.index, ft_ongoing.index)
at_risk = ft_complete(ts) + ft_ongoing(ts) + surv_complete(ts) + surv_ongoing(ts)
hs = ft_complete(ts) / at_risk
return Hazard(hs, ts)
Aşağıdaki işlev bir katılımcı grubundan sağkalım sürelerini çıkarır, tehlikeyi estimate_hazard ile tahmin eder ve karşılık gelen sağkalım fonksiyonunu hesaplar.
def estimate_survival(group):
"""Sağkalım fonksiyonunu tahmin eder."""
complete = group.query("evrmarry == 1")["agemarr"]
ongoing = group.query("evrmarry == 0")["age"]
hf = estimate_hazard(complete, ongoing)
sf = hf.make_surv()
return sf
Bu işlevlerle yakında sağkalım fonksiyonlarının güven aralıklarını hesaplayacağız. Önce Kaplan–Meier tahmini yapmanın başka bir yoluna bakalım.
13.7. Lifelines#
lifelines adlı Python paketi Kaplan–Meier tahminlerini hesaplayan işlevlerin de bulunduğu sağkalım analizi araçları sağlar.
Aşağıdaki hücre gerekirse lifelines paketini kurar.
try:
import lifelines
except ImportError:
%pip install lifelines
Önceki sonucun doğruluğunu bununla onaylayabiliriz. Önce sağkalım fonksiyonunu estimate_survival ile hesaplayacağız.
surv = estimate_survival(resp60)
Ardından lifelines ile hesaplayacağız. Önce verileri paketin istediği biçime getirelim.
complete = complete.dropna()
durations = np.concatenate([complete, ongoing])
event_observed = np.concatenate([np.ones(len(complete)), np.zeros(len(ongoing))])
Şimdi KaplanMeierFitter nesnesi oluşturup verileri uydurabiliriz.
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()
kmf.fit(durations=durations, event_observed=event_observed)
<lifelines.KaplanMeierFitter:"KM_estimate", fitted with 15389 total observations, 5468 right-censored observations>
Uyumdan sonra plot ile tahmini sağkalım fonksiyonunu ve güven aralığını gösterebiliriz. Ancak güven aralığı tabakalı örneklemeyi düzeltmediği için burada doğru değildir.
kmf.plot()
decorate(xlabel="Yaş", ylabel="Hiç evlenmeme olasılığı", ylim=ylim)
Bizim sağkalım fonksiyonumuzdan farklı olarak lifelines sonucu 0’dan başlar; geri kalanı kayan nokta hatası sınırında aynıdır.
ps = kmf.survival_function_["KM_estimate"].drop(0)
np.allclose(ps, surv)
True
Sonraki bölümde tabakalı örneklemeyi hesaba katan güven aralıklarını ağırlıklı yeniden örneklemeyle hesaplayacağız.
13.8. Güven aralıkları#
Kaplan–Meier tahminimiz veri kümesinin tek bir yeniden örneklemesine dayanır. Rastgele örneklemeden kaynaklanan değişimi görmek için analizi birkaç yeniden örneklemeyle çalıştırıp sonuçları çizeceğiz.
Aşağıdaki işlev bir DataFrame ve kohort listesi alır, her kohortun sağkalım fonksiyonunu tahmin eder ve tam sayı kohortları Surv nesnelerine eşleyen sözlük döndürür.
make_survival_map ile aynıdır; ancak sansürlü veri yoksa çalışan Surv.from_seq yerine Kaplan–Meier kullanan estimate_survival çağırır.
def estimate_survival_map(resp, cohorts):
"""Kohortları Surv nesnelerine eşleyen sözlük oluşturur."""
surv_map = {}
grouped = resp.groupby("cohort")
for cohort in cohorts:
group = grouped.get_group(cohort)
surv_map[cohort] = estimate_survival(group)
return surv_map
Aşağıdaki döngü veri kümesinin 101 rastgele yeniden örneklemesini üretir ve tahmini sağkalım fonksiyonlarını içeren 101 sözlükten oluşan liste oluşturur.
cohorts = [1940, 1950, 1960, 1970, 1980, 1990]
surv_maps = [estimate_survival_map(resample_cycles(resp), cohorts) for i in range(101)]
Sonuçları çizmek için sözlük listesini, bir tam sayı kohortu ve renk dizesini alan aşağıdaki işlevi kullanacağız. Sözlüklerde dolaşıp kohortun sağkalım fonksiyonunu neredeyse saydam çizgiyle çizer; yeniden örneklemeler arasındaki değişimi görselleştirmenin bir yoludur.
def plot_cohort(surv_maps, cohort, color):
"""Tek bir kohortun sonuçlarını çizer."""
survs = [surv_map[cohort] for surv_map in surv_maps]
for surv in survs:
surv.plot(color=color, alpha=0.05)
x, y = surv.index[-1], surv.iloc[-1]
plt.text(x + 1, y, f"{cohort}s", ha="left", va="center")
1940’lardan 1990’lara doğum kohortlarının sonuçları şöyledir.
colors = [f"C{i}" for i in range(len(cohorts))]
for cohort, color in zip(cohorts, colors):
plot_cohort(surv_maps, cohort, color)
xlim = [8, 55]
decorate(xlabel="Yaş", ylabel="Hiç evlenmeme olasılığı", xlim=xlim, ylim=ylim)
Bu görselleştirme keşif için yeterlidir; ancak çizgiler bulanık, bazı etiketler üst üstedir. Yayına hazır şekil için daha fazla çalışma gerekir; şimdilik basit tutacağız.
Şu örüntüler görülüyor:
1940’larda doğan kadınlar en erken evlendi. 1950’ler ve 1960’lar daha geç evlendi, ancak yaklaşık aynı oran evlenmeden kaldı.
1970’lerde doğanlar hem daha geç evlendi hem önceki kohortlardan daha yüksek oranda evlenmeden kaldı.
1980’ler ve 1990’lar daha da geç evleniyor ve daha yüksek oranlarda evlenmeden kalacak gibi görünüyor; gelecekte örüntü değişebilir.
Daha fazlası için NSFG’nin sonraki veri yayımını beklemeliyiz.
13.9. Beklenen kalan ömür#
Bir dağılım verildiğinde geçen sürenin işlevi olarak beklenen kalan ömrü hesaplayabiliriz. Örneğin gebelik süreleri dağılımından doğuma kadar beklenen süreyi bulabiliriz. Göstermek için NSFG gebelik verilerini kullanacağız.
Aşağıdaki hücreler veri dosyalarını indirir ve okumak için gereken statadict paketini kurar.
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/nb/nsfg.py")
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/2002FemPreg.dct")
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/2002FemPreg.dat.gz")
try:
import statadict
except ImportError:
%pip install statadict
Verileri okumak ve ilk bebeklerle diğerlerini ayırmak için get_nsfg_groups kullanacağız.
from nsfg import get_nsfg_groups
live, firsts, others = get_nsfg_groups()
Verinin tek bir yeniden örneklemesiyle başlayacağız.
sample = resample_rows_weighted(live, "finalwgt")
Gebelik sürelerinin PMF’si şöyledir.
pmf_durations = Pmf.from_seq(sample["prglngth"])
Gebeliğin 36. haftasının başında olduğumuzu varsayın. En yaygın sürenin 39 hafta olduğunu anımsayarak kalan sürenin 3–4 hafta olmasını bekleriz. Tahmini kesinleştirmek için dağılımda 36 haftaya eşit ya da büyük değerleri seçebiliriz.
t = 36
is_remaining = pmf_durations.qs >= t
Şimdi yalnızca bu değerleri içeren ve şimdiki zamanı 0’a getirecek biçimde sola kaydırılmış yeni bir Pmf oluşturacağız.
ps = pmf_durations.ps[is_remaining]
qs = pmf_durations.qs[is_remaining] - t
pmf_remaining = Pmf(ps, qs)
Pmf değerlerinin alt kümesini seçtiğimiz için olasılıkların toplamı artık 1 değildir; normalize ederek düzeltebiliriz.
pmf_remaining.normalize()
np.float64(0.9155006558810669)
haftanın başındaki kalan süre dağılımı şöyledir.
pmf_remaining.bar(label="36. hafta")
decorate(xlabel="Kalan süre (hafta)", ylabel="PMF")
Bu dağılımın ortalaması beklenen kalan süredir.
pmf_remaining.mean()
np.float64(3.2145671641791043)
Aşağıdaki işlev bu adımları kapsayıp belirli bir Pmf ve t zamanı için kalan süre dağılımını hesaplar.
def compute_pmf_remaining(pmf, t):
"""Kalan sürenin dağılımı."""
is_remaining = pmf.qs >= t
ps = pmf.ps[is_remaining]
qs = pmf.qs[is_remaining] - t
pmf_remaining = Pmf(ps, qs)
pmf_remaining.normalize()
return pmf_remaining
Aşağıdaki işlev gebelik sürelerinin Pmf nesnesini alıp 36.–43. haftaların başındaki beklenen kalan süreyi hesaplar.
def expected_remaining(pmf):
index = range(36, 44)
expected = pd.Series(index=index)
for t in index:
pmf_remaining = compute_pmf_remaining(pmf, t)
expected[t] = pmf_remaining.mean()
return expected
Verinin tek bir yeniden örneklemesinin sonuçları şöyledir.
expected = expected_remaining(pmf_durations)
expected
36 3.214567
37 2.337714
38 1.479095
39 0.610133
40 0.912517
41 0.784211
42 0.582301
43 0.589372
dtype: float64
Rastgele örneklemeden kaynaklanan değişimi görmek için analizi birkaç yeniden örneklemeyle çalıştırıp sonuçları çizebiliriz.
for i in range(21):
sample = resample_rows_weighted(live, "finalwgt")
pmf_durations = Pmf.from_seq(sample["prglngth"])
expected = expected_remaining(pmf_durations)
expected.plot(color="C0", alpha=0.1)
decorate(
xlabel="Gebelik haftası", ylabel="Beklenen kalan süre (hafta)", ylim=[0, 3.4]
)
ile 39. haftalar arasında beklenen kalan süre azalır ve 39. haftanın başında yaklaşık 0,6 haftadır. Sonra eğri yatıklaşır. 40. haftanın başında hâlâ yaklaşık 0,6 haftadır —hatta biraz daha yüksek— ve 41., 42. ve 43. haftaların başında neredeyse aynıdır. Bebeğin doğumunu heyecanla bekleyenler için bu davranış oldukça acımasız görünüyor.
13.10. Sözlük#
sağkalım analizi (survival analysis): Genellikle ömür veya sürelere odaklanarak ilgilenilen olaya kadar geçen zamanı açıklama ve öngörme yöntemleri.
sağkalım fonksiyonu (survival function): \(t\) zamanını \(t\) sonrasına sağkalma olasılığına eşleyen fonksiyon.
tehlike fonksiyonu (hazard function): \(t\)’ye kadar sağ kalan olgular arasında olayı \(t\)’de yaşayanların oranını gösteren fonksiyon.
kümülatif tehlike fonksiyonu: Tehlike fonksiyonunun görselleştirmede yararlı kümülatif toplamı.
ağırlıklı bootstrap: Temsil edici örneklemi simüle ederek tabakalı örneklemeyi düzeltmek için örnekleme ağırlıkları kullanan yeniden örnekleme.
sansürlü veri (censored data): İlgilenilen olay henüz gerçekleşmediği ya da gözlenmediği için yalnızca kısmen bilinen veri.
Kaplan–Meier tahmini: Sansürlü gözlemler içeren verilerde sağkalım ve tehlike fonksiyonlarını tahmin etme yöntemi.
kohort (cohort): Doğum on yılı gibi ortak özellikleri olan ve grup olarak analiz edilen denekler.
13.11. Alıştırmalar#
13.11.1. Alıştırma 13.1#
Bu bölümün yöntemleriyle evlilik süresinin tehlike ve sağkalım fonksiyonlarını tahmin edebiliriz. Basitlik için yalnızca ilk evlilikleri ve ayrılık ya da ölüm yerine bitiş olayı olarak boşanmayı ele alacağız.
NSFG’deki cmdivorcx, uygunsa ilk evliliğin boşanma tarihini yüzyıl ayıyla içerir. Boşanmayla biten evliliklerin ve hâlen süren evliliklerin bugüne kadarki süresini hesaplayın.
Tamamlanmış olgular için
cmdivorcxilecmmarrhxarasındaki süreyi hesaplayın. İki değer de geçerliyse ilk evlilik boşanmayla bitmiştir.Süren olgular için yalnızca bir kez evlenmiş ve hâlâ evli kişileri seçin. Evlilik sayısını veren
fmarnoile gayriresmî medeni durumu kodlayanrmaritaldeğişkenlerini kullanabilirsiniz; 1 evli demektir.
Bazı değerler yaklaşık olduğundan az sayıda negatif fark bulabilirsiniz; bir yıldan fazla olmamalıdır.
Evlilik süresinin tehlike ve sağkalım fonksiyonlarını tahmin edin. Kümülatif tehlikeyi çizip boşanma tehlikesinin ne zaman en yüksek olduğunu bulun. Sağkalım fonksiyonunu çizip evliliklerin ne kadarının boşanmayla bittiğini belirleyin.
# Verinin tek bir yeniden örneklemesini kullanmanızı öneririm
sample = resample_cycles(resp)
13.11.2. Alıştırma 13.2#
2012’de Güney Kaliforniya Üniversitesi’nden demograflar 1800’lerin ve 1900’lerin başında İsveç’te doğanların yaşam beklentisini tahmin etti. 0–91 yaşlarında, o yaşa kadar sağ kalanlar arasındaki ölüm oranı olan yaşa özgü ölüm hızını —tehlike fonksiyonunu— tahmin ettiler.
Makalelerindeki şeklin verilerini çevrimiçi grafik sayısallaştırıcıyla çıkarıp CSV dosyasına kaydettim. İndirme yönergeleri notebook’tadır.
Veri kaynağı: Beltrán-Sánchez, H., Crimmins, E. M. ve Finch, C. E. (2012). Early cohort mortality predicts the rate of aging in the cohort: a historical analysis. Journal of Developmental Origins of Health and Disease, 3(5), 380–386.
Aşağıdaki hücre verileri indirir.
download(
"https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/mortality_rates_beltran2012.csv"
)
Verileri şöyle yükleyebiliriz.
mortality = pd.read_csv("mortality_rates_beltran2012.csv", header=[0, 1]).dropna()
Aşağıdaki işlev 0–99 arasındaki her yaş için yaklaşık ölüm hızlarını içeren tehlike fonksiyonu oluşturmak üzere verilere ara değer biçer.
from scipy.interpolate import interp1d
from empiricaldist import Hazard
def make_hazard(ages, rates):
"""Bir Series üzerinde ara değer biçerek Hazard fonksiyonu oluşturur.
series: Series
döndürür: Hazard
"""
interp = interp1d(ages, rates, fill_value="extrapolate")
xs = np.arange(0, 100)
ys = np.exp(interp(xs))
return Hazard(ys, xs)
Şimdi şöyle bir Hazard nesnesi oluşturabiliriz.
ages = mortality["1800", "X"].values
rates = mortality["1800", "Y"].values
hazard = make_hazard(ages, rates)
Ölüm hızları şöyle görünür.
hazard.plot()
decorate(xlabel="Yaş (yıl)", ylabel="Tehlike")
Bu hızlardan sağkalım fonksiyonu oluşturmak için make_surv, karşılık gelen CDF için make_cdf kullanıp sonuçları çizin.
Ardından ömür dağılımını temsil eden Pmf için make_pmf kullanıp grafiğini çizin. Son olarak 0–99 arasındaki her yaşta ortalama kalan ömrü compute_pmf_remaining ile hesaplayıp çizin.
Kalan ömür eğrisinde sezgiye aykırı bir örüntü görmelisiniz: yaşamın ilk yıllarında kalan ömür artar. 1800’lerin başında bebek ölümleri çok yüksek olduğundan daha büyük çocuğun beklenen ömrü daha küçük çocuktan uzundu. Yaklaşık 5 yaştan sonra yaşam beklentisi beklediğimiz örüntüye döner; gençlerin yaşlılardan uzun yaşaması beklenir.
Konuyla ilgileniyorsanız Probably Overthinking It kitabımın çeşitli istatistik alanlarından benzer sezgiye aykırı sonuçlar sunan 5. bölümünü sevebilirsiniz.
https://probablyoverthinking.it adresinden edinilebilir.
Think Stats: Python ile Keşifsel Veri Analizi, 3. Baskı
Telif hakkı 2024 Allen B. Downey
Kod lisansı: MIT Lisansı
Metin lisansı: Creative Commons Atıf-Gayriticari-AynıLisanslaPaylaş 4.0 Uluslararası