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.

10. En küçük kareler#

Bu bölüm ve sonraki bölüm verilere model uydurma düşüncesini tanıtır. Bu bağlamda model, değişkenler arasındaki ilişkinin —örneğin düz bir çizgi— matematiksel açıklaması ile rastgele değişimin —örneğin normal dağılım— açıklamasından oluşur.

Bir modelin verilere uyduğunu söylerken genellikle model ile veri arasındaki uzaklıkları, yani hataları, en aza indirdiğini kastederiz. Model uydurmanın en yaygın yollarından biriyle başlayacağız: karesel hatalar toplamını en aza indiren en küçük kareler uyumu.

Önce aynı anda yalnızca iki değişkenle çalışan modelleri ele alacağız. Sonraki bölüm ikiden fazla değişkeni işleyebilen modelleri tanıtacak.

Bu notebook’u Colab’da çalıştırmak için buraya tıklayın.

Kod hücresi içeriğini gizle

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/nb/thinkstats.py")

Kod hücresi içeriğini gizle

try:
    import empiricaldist
except ImportError:
    %pip install empiricaldist

Kod hücresi içeriğini gizle

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from thinkstats import decorate

10.1. En küçük kareler uyumu#

İlk örnek olarak Bölüm 8’deki senaryoya dönelim. Antarktika’da yerel penguen popülasyonlarını inceleyen bir araştırmacı olduğunuzu varsayın. Veri toplama kapsamında bir penguen örneklemi yakalıyor, ölçüp tartıyor ve zarar vermeden serbest bırakıyorsunuz.

Penguenleri doğru ölçüm yapmaya yetecek kadar tartıda tutmanın zor olduğunu kısa sürede öğrenirdiniz. Bazı penguenler için yüzgeç ve gaga ölçümlerimiz olduğunu, ancak ağırlıklarının olmadığını varsayalım. Eksik verileri diğer ölçümlerle doldurup dolduramayacağımıza bakalım; bu sürece eksik değer atama denir.

Antarktika’daki Palmer İstasyonu araştırmacılarının 2007–2010 arasında topladığı verilerle ağırlık ve ölçümler arasındaki ilişkiyi inceleyerek başlayacağız. Veriler ücretsizdir; indirme yönergeleri bu bölümün notebook’unda yer alır.

Aşağıdaki hücre verileri Allison Horst tarafından oluşturulan bir depodan indirir.

Horst AM, Hill AP, Gorman KB (2020). palmerpenguins: Palmer Archipelago (Antarctica) penguin data. R paketi, sürüm 0.1.0. https://allisonhorst.github.io/palmerpenguins/. doi: 10.5281/zenodo.3960218.

Veriler şu makaleye temel olan araştırma kapsamında toplanmıştır: Gorman KB, Williams TD, Fraser WR (2014). Ecological sexual dimorphism and environmental variability within a community of Antarctic penguins (genus Pygoscelis). PLoS ONE 9(3):e90081. https://doi.org/10.1371/journal.pone.0090081

download(
    "https://raw.githubusercontent.com/allisonhorst/palmerpenguins/c19a904462482430170bfe2c718775ddb7dbb885/inst/extdata/penguins_raw.csv"
)

Verileri read_csv ile okuyabiliriz.

penguins = pd.read_csv("penguins_raw.csv").dropna(subset=["Body Mass (g)"])
penguins.shape
(342, 17)

Veri kümesi 151 Adélie pengueninin ölçümlerini içerir. Bu verileri içeren satırları query ile seçebiliriz.

adelie = penguins.query('Species.str.startswith("Adelie")')
len(adelie)
151

Bir Adélie pengueninin yüzgeç uzunluğunu bildiğimizi varsayalım; ağırlığını ne kadar iyi tahmin edebileceğimize bakalım. Önce ilgili sütunları DataFrame’den seçeceğiz.

xvar = "Flipper Length (mm)"
yvar = "Body Mass (g)"

flipper_length = adelie[xvar]
body_mass = adelie[yvar]

Bu nicelikler arasındaki ilişkiyi gösteren saçılım grafiği şöyledir.

plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
decorate(xlabel=xvar, ylabel=yvar)
_images/3f28135b7ff2230ba387ca5d467155d6d1b46e4f0bdafe88dffd3d14e8c11f42.png

İlişkili görünüyorlar; korelasyon katsayısını hesaplayarak ilişkinin gücünü nicelleştirebiliriz.

np.corrcoef(flipper_length, body_mass)[0, 1]
np.float64(0.4682016942179394)

Korelasyon yaklaşık 0,47’dir; yani daha uzun yüzgeçli penguenler daha ağır olma eğilimindedir. Yüzgeç uzunluğunu bildiğimizde ağırlığı daha doğru tahmin edebilmemiz yararlıdır; ancak yalnızca korelasyon bu tahmini nasıl yapacağımızı söylemez. Bunun için bir en iyi uyum doğrusuna gereksinimimiz vardır.

“En iyi” doğruyu tanımlamanın birçok yolu vardır. Bu tür veriler için yaygın bir seçim, ortalama karesel hatayı (MSE) en aza indiren doğru olan doğrusal en küçük kareler uyumudur.

SciPy, en küçük kareler uyumunu hesaplayan linregress işlevini sağlar. Adı bu tür bir modelin diğer adı olan doğrusal regresyonun kısaltmasıdır. Argümanları sırasıyla x ve y değerleridir.

from scipy.stats import linregress

result = linregress(flipper_length, body_mass)
result
LinregressResult(slope=np.float64(32.831689751150094), intercept=np.float64(-2535.8368022002524), rvalue=np.float64(0.4682016942179394), pvalue=np.float64(1.3432645947790076e-09), stderr=np.float64(5.076138407990822), intercept_stderr=np.float64(964.7984274994061))

Sonuç, uydurulan doğrunun eğimini ve kesişimini, ayrıca birazdan inceleyeceğimiz başka bilgileri içeren bir LinregressResult nesnesidir. Eğim yaklaşık 32,8’dir; yüzgeç uzunluğundaki her ek milimetre, vücut ağırlığında 32,8 gramlık artışla ilişkilidir.

Kesişim -2535 gramdır; ölçülen ağırlık negatif olamayacağından anlamsız görünebilir. Eğimi ve kesişimi kullanıp uydurulan doğruyu ortalama yüzgeç uzunlığında değerlendirirsek daha anlamlı olabilir.

x = flipper_length.mean()
y = result.intercept + result.slope * x
x, y
(np.float64(189.95364238410596), np.float64(3700.662251655629))

Ortalama yüzgeç uzunluğuna, yaklaşık 190 mm’ye sahip bir penguen için beklenen vücut ağırlığı yaklaşık 3700 gramdır.

Aşağıdaki işlev linregress sonucunu ve bir xs dizisini alıp her x değeri için uydurulan doğru üzerindeki noktayı bulur.

def predict(result, xs):
    ys = result.intercept + result.slope * xs
    return ys

predict adı burada garip görünebilir. Gündelik dilde öngörü çoğu zaman gelecekte olacak bir şeyle ilgilidir; ancak regresyon bağlamında uydurulan doğru üzerindeki noktalara da öngörü denir.

Bir yüzgeç uzunluğu aralığında doğru üzerindeki noktaları hesaplamak için predict kullanabiliriz.

fit_xs = np.linspace(np.min(flipper_length), np.max(flipper_length))
fit_ys = predict(result, fit_xs)

Verilerin saçılım grafiğiyle birlikte uydurulan doğru şöyledir.

plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)
plt.plot(fit_xs, fit_ys, color="C1")
decorate(xlabel=xvar, ylabel=yvar)
_images/da49d28ee033a69169fc186a65345eeab123a31e230aa40feb28a75e2076e491.png

Beklendiği gibi uydurulan doğru verinin merkezinden geçer ve eğilimi izler. Bazı öngörüler doğrudur; ancak birçok veri noktası doğrudan uzaktadır. Öngörülerin ne kadar iyi ya da kötü olduğunu anlamak için her noktanın doğruya dikey uzaklığı olan öngörü hatasını hesaplayabiliriz. Aşağıdaki işlev artıklar da denen bu hataları hesaplar.

def compute_residuals(result, xs, ys):
    fit_ys = predict(result, xs)
    return ys - fit_ys

Yüzgeç uzunluğunun işlevi olarak vücut kütlesi artıklarının dağılımı şöyledir.

residuals = compute_residuals(result, flipper_length, body_mass)

Örnek olarak veri kümesindeki ilk penguenin sonuçlarına bakabiliriz.

x = flipper_length[0]
y = predict(result, x)
x, y
(np.float64(181.0), np.float64(3406.699042757915))

Seçilen penguenin yüzgeç uzunluğu 181 mm ve öngörülen vücut kütlesi 3407 gramdır. Şimdi gerçek kütleye bakalım.

body_mass[0], residuals[0]
(np.float64(3750.0), np.float64(343.30095724208513))

Bu penguenin gerçek kütlesi 3750 gram, öngörü çıkarıldıktan sonra kalan artık ise 343 gramdır.

Karesel artıkların ortalaması öngörülerin ortalama karesel hatasıdır (MSE).

mse = np.mean(residuals**2)
mse
np.float64(163098.85902884742)

Bu sayı tek başına pek anlamlı değildir. Belirleme katsayısını hesaplayarak daha anlamlı hâle getirebiliriz.

10.2. Belirleme katsayısı#

Bir penguenin ağırlığını tahmin etmek istediğinizi varsayın. Yüzgeç uzunluğunu biliyorsanız tahmininizi en küçük kareler uyumuna dayandırabilir; MSE ile tahminlerinizin ortalama doğruluğunu nicelleştirebilirsiniz.

Peki yüzgeç uzunluğunu bilmiyorsanız ne tahmin ederdiniz? MSE’yi en aza indirmesi bakımından en iyi stratejinin ortalamayı tahmin etmek olduğu ortaya çıkar. Her zaman ortalamayı tahmin edersek öngörü hataları ortalamadan sapmalardır.

deviations = body_mass - np.mean(body_mass)

MSE ise karesel sapmaların ortalamasıdır.

np.mean(deviations**2)
np.float64(208890.28989956583)

Ortalama karesel sapmanın varyans olduğunu anımsayabilirsiniz.

np.var(body_mass)
np.float64(208890.28989956583)

Dolayısıyla kütlelerin varyansını her zaman ortalamayı tahmin ettiğimizde oluşan MSE, artıkların varyansını ise regresyon doğrusunu kullandığımızda oluşan MSE olarak düşünebiliriz. Bu varyansların oranını hesaplayıp 1’den çıkarırsak sonuç, tahminlerimizde yüzgeç uzunluğunu kullandığımızda MSE’nin ne kadar azaldığını gösterir.

Aşağıdaki işlev teknik adı belirleme katsayısı olan bu değeri hesaplar. Değer \(R^2\) ile gösterildiğinden çoğu kişi buna “R kare” der.

def coefficient_of_determination(ys, residuals):
    return 1 - np.var(residuals) / np.var(ys)

Örnekte \(R^2\) yaklaşık 0,22’dir; yani uydurulan doğru MSE’yi %22 azaltır.

R2 = coefficient_of_determination(body_mass, residuals)
R2
np.float64(0.21921282646854878)

Belirleme katsayısı \(R^2\) ile korelasyon katsayısı \(r\) arasında ilişki vardır. Gösterimden tahmin edebileceğiniz gibi \(r^2 = R^2\)’dir.

\(R^2\)’nin karekökünü hesaplayarak bunun doğru olduğunu gösterebiliriz.

r = np.sqrt(R2)
r
np.float64(0.46820169421793933)

Ardından daha önce hesapladığımız korelasyonla karşılaştırabiliriz.

corr = np.corrcoef(flipper_length, body_mass)[0, 1]
corr
np.float64(0.4682016942179394)

Kayan nokta yaklaşımından kaynaklanan küçük bir fark dışında aynıdırlar.

linregress de bu değeri hesaplar ve RegressionResult nesnesinde bir öznitelik olarak döndürür.

result.rvalue
np.float64(0.4682016942179394)

Belirleme ve korelasyon katsayıları çoğunlukla aynı bilgiyi taşır, ancak farklı yorumlanır:

  • Korelasyon, ilişkinin gücünü -1 ile 1 arasında nicelleştirir.

  • \(R^2\), uydurulan doğrunun MSE’yi azaltma yeteneğini nicelleştirir.

Ayrıca \(R^2\) her zaman pozitiftir; korelasyonun pozitif mi negatif mi olduğunu göstermez.

10.3. MSE’nin en aza indirilmesi#

Daha önce en küçük kareler uyumunun MSE’yi en aza indiren doğru olduğunu söyledim. Bunu kanıtlamayacağız; ancak kesişime ve eğime küçük rastgele değerler ekleyip MSE’nin kötüleşip kötüleşmediğine bakarak sınayabiliriz.

intercept = result.intercept + np.random.normal(0, 1)
slope = result.slope + np.random.normal(0, 1)

Sınamayı çalıştırmak için intercept ve slope özniteliklerine sahip bir nesne gerekir; types modülünün sağladığı SimpleNamespace nesnesini kullanacağız.

from types import SimpleNamespace

fake_result = SimpleNamespace(intercept=intercept, slope=slope)
fake_result
namespace(intercept=np.float64(-2536.2119265008996),
          slope=np.float64(33.19402369403873))

Bu nesneyi compute_residuals işlevine verip MSE’yi artıklardan hesaplayabiliriz.

fake_residuals = compute_residuals(fake_result, flipper_length, body_mass)
fake_mse = np.mean(fake_residuals**2)

Sonucu en küçük kareler doğrusunun MSE’siyle karşılaştırırsak her zaman daha kötü olduğunu görürüz.

mse, fake_mse, fake_mse > mse
(np.float64(163098.85902884742), np.float64(167790.04789653933), np.True_)

MSE’yi en aza indirmek yararlıdır, ancak “en iyi”nin tek tanımı değildir. Hataların mutlak değerlerini ya da her noktadan doğruya en kısa uzaklık olan “toplam hatayı” en aza indirmek başka seçeneklerdir. Bazı bağlamlarda fazla yüksek tahmin etmek fazla düşük tahmin etmekten daha iyi ya da daha kötü olabilir. Bu durumda her artık için maliyet fonksiyonu hesaplayıp toplam maliyeti en aza indirmek isteyebilirsiniz.

Yine de en küçük kareler uyumu, esas olarak hesaplama verimliliği nedeniyle bu seçeneklerden çok daha yaygındır. Aşağıdaki işlev bunun nasıl yapıldığını gösterir.

def least_squares(xs, ys):
    xbar = np.mean(xs)
    ybar = np.mean(ys)

    xdev = xs - xbar
    ydev = ys - ybar

    slope = np.sum(xdev * ydev) / np.sum(xdev**2)
    intercept = ybar - slope * xbar

    return intercept, slope

Bu işlevi sınamak için yine yüzgeç uzunluğu ile vücut kütlesini kullanacağız.

intercept, slope = least_squares(flipper_length, body_mass)
intercept, slope
(np.float64(-2535.8368022002524), np.float64(32.831689751150094))

linregress ile aynı sonuçları elde ettiğimizi doğrulayabiliriz.

np.allclose([intercept, slope], [result.intercept, result.slope])
True

Hesaplama verimliliğinin probleme en uygun yöntemi seçmekten daha önemli olduğu dönemde MSE’yi en aza indirmek anlamlıydı. Artık durum böyle olmadığından karesel artıkların gerçekten en aza indirilmesi gereken şey olup olmadığını düşünmeye değer.

10.4. Tahmin#

slope ve intercept parametreleri bir örnekleme dayalı tahminlerdir. Diğer tahminler gibi temsil edici olmayan örneklemeye, ölçüm hatasına ve rastgele örnekleme değişkenliğine açıktırlar. İlk ikisinin etkisini nicelleştirmek zordur; rastgele örneklemenin etkisini nicelleştirmek daha kolaydır.

Bunu yapmanın bir yolu bootstrap yeniden örneklemesidir: örneklemi evrenin tamamı gibi kabul edip gözlenen veriden yerine koyarak özgünle aynı büyüklükte yeni örneklemler çekeriz. Aşağıdaki işlev bir DataFrame alır, satırları sample yöntemiyle yeniden örnekler ve yeni bir DataFrame döndürür.

def resample(df):
    n = len(df)
    return df.sample(n, replace=True)

Aşağıdaki işlev bir DataFrame alır, en küçük kareler uyumunu bulur ve uydurulan doğrunun eğimini döndürür.

def estimate_slope(df):
    xs, ys = df["Flipper Length (mm)"], df["Body Mass (g)"]
    result = linregress(xs, ys)
    return result.slope

Bu işlevlerle çok sayıda simüle edilmiş veri kümesi üretip her birinin eğimini hesaplayabiliriz.

# Her çalıştırmada aynı sonuçları almak için rastgele sayı üretecini tohumla

np.random.seed(1)
resampled_slopes = [estimate_slope(resample(adelie)) for i in range(1001)]

Sonuç, eğimin örnekleme dağılımından bir örneklemdir. Şöyle görünür.

from thinkstats import plot_kde

plot_kde(resampled_slopes)
decorate(xlabel="Uydurulan doğrunun eğimi (g/mm)", ylabel="Yoğunluk")
_images/b938c79244b4eaa74453f2505fa8d2ced90f9a220ca754c9e1eb019fbe1c0b68.png
ci90 = np.percentile(resampled_slopes, [5, 95])
print(result.slope, ci90)
32.831689751150094 [25.39604591 40.21054526]

Dolayısıyla tahmini eğimi 33 gram/mm ve %90 GA’yı [25, 40] gram/mm olarak bildirebiliriz.

Tahminin standart hatası örnekleme dağılımının standart sapmasıdır.

stderr = np.std(resampled_slopes)
stderr
np.float64(4.570238986584832)

linregress ile elde ettiğimiz RegressionResult, dağılımın biçimine ilişkin bazı varsayımlara dayanarak standart hatanın bir yaklaşımını sağlar.

result.stderr
np.float64(5.076138407990822)

Yeniden örneklemeyle hesapladığımız standart hata biraz küçüktür; ancak fark uygulamada muhtemelen önemli değildir.

10.5. Belirsizliğin görselleştirilmesi#

Veri kümesini her yeniden örneklediğimizde farklı bir uyum doğrusu elde ederiz. Doğrulardaki değişimi görmek için bunların tümünü bir döngüyle çizebiliriz. Aşağıdaki işlev yeniden örneklenmiş bir DataFrame alır, en küçük kareler uyumunu hesaplar ve bir xs dizisi için öngörülen değerleri üretir.

def fit_line(df, fit_xs):
    xs, ys = df["Flipper Length (mm)"], df["Body Mass (g)"]
    result = linregress(xs, ys)
    fit_ys = predict(result, fit_xs)
    return fit_ys

Kullanacağımız xs dizisi şöyledir.

xs = adelie["Flipper Length (mm)"]
fit_xs = np.linspace(np.min(xs), np.max(xs))

Verilerin saçılım grafiğiyle birlikte uydurulan doğrular şöyle görünür.

plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)

for i in range(101):
    fit_ys = fit_line(resample(adelie), fit_xs)
    plt.plot(fit_xs, fit_ys, color="C1", alpha=0.05)

decorate(xlabel=xvar, ylabel=yvar)
_images/1bf348c41e42617c9f7a736611bfff684a0fb7b06a118d8e94d560176fbe427a.png

Ortaya yakın yerde doğrular birbirine yakındır; uçlarda daha fazla ayrılırlar.

Doğruların değişkenliğini göstermenin başka bir yolu, öngörülen her değer için %90 güven aralığı çizmektir. Bunu, uydurulan doğruları dizilerden oluşan bir listede toplayarak yapabiliriz.

fitted_ys = [fit_line(resample(adelie), fit_xs) for i in range(1001)]

Bu dizi listesini her uyum doğrusu için bir satırı ve xs değerlerinin her biri için bir sütunu olan iki boyutlu bir dizi gibi düşünebiliriz.

Her xs değeri için karşılık gelen ys değerlerinin 5., 50. ve 95. yüzdeliklerini bulmak üzere axis=0 argümanıyla percentile kullanabiliriz.

low, median, high = np.percentile(fitted_ys, [5, 50, 95], axis=0)

Şimdi fill_between ile 5. ve 95. yüzdelikler arasındaki %90 GA bölgesini; ayrıca her sütundaki medyanı ve verilerin saçılım grafiğini çizeceğiz.

plt.scatter(flipper_length, body_mass, marker=".", alpha=0.5)

plt.fill_between(fit_xs, low, high, color="C1", lw=0, alpha=0.2)
plt.plot(fit_xs, median, color="C1")

decorate(xlabel=xvar, ylabel=yvar)
_images/d6b7424d47a32ebfb7408d84b9d67140f8e38d0560bc9d37872f6e1df01ebd4d.png

Rastgele örneklemeden kaynaklanan uyum doğrusu değişkenliğini göstermenin en sevdiğim yolu budur.

10.6. Dönüşüm#

Verilere doğru uydurmadan önce değerlerin karelerini, kareköklerini ya da logaritmalarını hesaplayarak değişkenlerden birini veya ikisini dönüştürmek bazen yararlıdır. Göstermek için Bölüm 5’te açıklanan BRFSS boy ve ağırlıklarını kullanacağız.

Aşağıdaki hücre BRFSS verilerini indirir.

download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/CDBRFS08.ASC.gz")

BRFSS verilerini şöyle yükleyebiliriz.

from thinkstats import read_brfss

brfss = read_brfss()

Ardından geçerli verisi olan satırları bulup boy ve ağırlık sütunlarını seçeceğiz.

valid = brfss.dropna(subset=["htm3", "wtkg2"])
heights, weights = valid["htm3"], valid["wtkg2"]

En küçük kareler uyumunun eğimini ve kesişimini linregress ile hesaplayabiliriz.

result_brfss = linregress(heights, weights)
result_brfss.intercept, result_brfss.slope
(np.float64(-82.65926054409942), np.float64(0.9570745850332297))

Eğim yaklaşık 0,96’dır; yani boydaki 1 santimetrelik artış ortalama olarak ağırlıktaki yaklaşık 1 kilogramlık artışa karşılık gelir. Bir xs aralığının öngörülen değerlerini üretmek için yine predict kullanabiliriz.

fit_xs = np.linspace(heights.min(), heights.max())
fit_ys = predict(result_brfss, fit_xs)

Saçılım grafiğini çizmeden önce boylara ve ağırlıklara titreşim eklemek yararlıdır.

from thinkstats import jitter

jittered_heights = jitter(heights, 2)
jittered_weights = jitter(weights, 1.5)

\(x\) ekseninin sınırlarını seçmek için boyların ortalaması ile standart sapmasını kullanacağız.

m, s = heights.mean(), heights.std()
xlim = m - 4 * s, m + 4 * s
ylim = 0, 200

Titreşim eklenmiş verilerin saçılım grafiği ve uydurulan doğru şöyledir.

plt.scatter(jittered_heights, jittered_weights, alpha=0.01, s=0.1)
plt.plot(fit_xs, fit_ys, color="C1")
decorate(xlabel="Boy (cm)", ylabel="Ağırlık (kg)", xlim=xlim, ylim=ylim)
_images/66c13f007c17659bd9738aa73852ecd9de9b4c725f72555ca8b6c10bbe09cd77.png

Uydurulan doğru saçılım grafiğinin en yoğun bölümünden geçmez. Bunun nedeni ağırlıkların normal dağılım izlememesidir. Bölüm 5’te gördüğümüz gibi yetişkin ağırlıkları büyük değerlere doğru çarpık olan lognormal dağılım izleme eğilimindedir; bu değerler doğruyu yukarı çeker.

Başka bir kaygı nedeni artıkların aşağıdaki dağılımıdır.

residuals = compute_residuals(result_brfss, heights, weights)
from thinkstats import make_pmf

pmf_kde = make_pmf(residuals, -60, 120)
pmf_kde.plot()

decorate(xlabel="Artık (kg)", ylabel="Yoğunluk")
_images/3333b2cf232999c5e7f179ffbf4ddf741ebfec259329556103ec86c9804a122f.png

Artıkların dağılımı sağa çarpıktır. Bu tek başına mutlaka sorun değildir; ancak en küçük kareler uyumunun değişkenler arasındaki ilişkiyi uygun biçimde açıklamadığını düşündürür.

Ağırlıklar lognormal dağılım izliyorsa logaritmaları normal dağılım izler. Ağırlık logaritmasını boyun işlevi olarak uydurursak ne olduğuna bakalım.

log_weights = np.log10(weights)
result_brfss2 = linregress(heights, log_weights)
result_brfss2.intercept, result_brfss2.slope
(np.float64(0.9930804163932864), np.float64(0.005281454169417781))

Değişkenlerden birini dönüştürdüğümüz için eğim ile kesişimi yorumlamak daha zordur. Yine de uydurulan doğruyu predict ile hesaplayabiliriz.

fit_xs = np.linspace(heights.min(), heights.max())
fit_ys = predict(result_brfss2, fit_xs)

Ardından dönüştürülmüş verilerin saçılım grafiğiyle birlikte çizebiliriz.

jittered_log_weights = jitter(log_weights, 1.5)
plt.scatter(jittered_heights, jittered_log_weights, alpha=0.01, s=0.1)
plt.plot(fit_xs, fit_ys, color="C1")
decorate(xlabel="Boy (cm)", ylabel="Ağırlık (log10 kg)", xlim=xlim)
_images/2c1dac461e99d43a2f83d1ee8f02d9609605ab4ae575aa862614f8b812cc8463.png

Uydurulan doğru grafiğin en yoğun bölümünden geçer ve gerçek değerler doğrunun üstünde ve altında yaklaşık eşit uzaklıklara uzanır; dolayısıyla artıkların dağılımı yaklaşık simetriktir.

residuals = compute_residuals(result_brfss2, heights, log_weights)
pmf_kde = make_pmf(residuals, -0.6, 0.6)
pmf_kde.plot()

decorate(xlabel="Artık (kg)", ylabel="Yoğunluk")
_images/ee8bbfeb6fde4566153235c7e6a95f535b560bf0fabedd89eff43476ea94e589.png

Saçılım grafiğinin görünümü ve artık dağılımı, boy ile logaritmik dönüşümlü ağırlık arasındaki ilişkinin uydurulan doğruyla iyi açıklandığını düşündürür. İki regresyonun \(r\) değerlerini karşılaştırdığımızda boy ile logaritmik dönüşümlü ağırlıkların korelasyonunun biraz daha yüksek olduğunu görürüz.

result_brfss.rvalue, result_brfss2.rvalue
(np.float64(0.5087364789734595), np.float64(0.5317282605983427))

Bu, \(R^2\) değerinin de biraz daha yüksek olduğu anlamına gelir.

result_brfss.rvalue**2, result_brfss2.rvalue**2
(np.float64(0.2588128050383132), np.float64(0.282734943118939))

Boylardan ağırlıkları tahmin ediyorsak logaritmik dönüşümlü ağırlıklarla çalıştığımızda tahminler biraz daha iyidir.

Ancak veriyi dönüştürmek model parametrelerini yorumlamayı zorlaştırır; sonuçları sunmadan önce dönüşümü tersine çevirmek yararlı olabilir. Örneğin 10 tabanındaki logaritmanın tersi, 10 tabanında üs almadır. Ters dönüşümden sonra uydurulan doğrunun dönüştürülmemiş verilerle görünümü şöyledir.

plt.scatter(jittered_heights, jittered_weights, alpha=0.01, s=0.1)
plt.plot(fit_xs, 10**fit_ys, color="C1")
decorate(xlabel="Boy (cm)", ylabel="Ağırlık (kg)", xlim=xlim, ylim=ylim)
_images/d60cc7311c5f81f786a87194b6ffc01ec99df07276bc4c654d208fb4e2b930a0.png

Dönüştürülmüş verilere göre düz olan uyum doğrusu, özgün verilere göre eğridir.

10.7. Sözlük#

  • model: Regresyon bağlamında değişkenler arasındaki ilişkinin matematiksel açıklaması ile rastgele değişimin açıklaması.

  • eksik değer atama (imputation): Veri kümesindeki eksik değerleri tahmin edip doldurma süreci.

  • en iyi uyum doğrusu (line of best fit): Bir ilişkiyi belirli bir “en iyi” tanımına göre en iyi açıklayan doğru ya da eğri.

  • doğrusal regresyon (linear regression): En iyi uyum doğrusunu bulma yöntemi.

  • öngörü (prediction): En iyi uyum doğrusu üzerindeki nokta; regresyonda mutlaka geleceğe ilişkin değildir.

  • artık (residual): Gözlenen değer ile uyum doğrusunun öngördüğü değer arasındaki fark.

  • doğrusal en küçük kareler uyumu: Karesel artıklar toplamını en aza indiren doğru.

  • belirleme katsayısı (coefficient of determination): \(R^2\) ile gösterilen ve modelin veriye uyumunu nicelleştiren istatistik.

  • bootstrap yeniden örneklemesi: Örneklemi evren kabul edip yerine koyarak özgünle aynı büyüklükte yeni örneklemler çekme yöntemi.

10.8. Alıştırmalar#

10.8.1. Alıştırma 10.1#

Bu bölümde penguen ağırlığını yüzgeç uzunluğunun işlevi olarak en küçük karelerle uydurduk. Veri kümesindeki iki ölçümü daha ele alabiliriz: gaga sırtı uzunluğu ve derinliği.

Ağırlığı gaga sırtı uzunluğunun işlevi olarak en küçük karelerle uydurun. Bu değişkenlerin saçılım grafiğini ve uydurulan doğruyu çizin.

RegressionResult nesnesinin rvalue özniteliğine göre değişkenlerin korelasyonu nedir? Belirleme katsayısı nedir? Ağırlığı hangisi daha iyi tahmin eder: gaga sırtı uzunluğu mu, yüzgeç uzunluğu mu?

xvar = "Culmen Length (mm)"
yvar = "Body Mass (g)"

culmen_length = adelie[xvar]
body_mass = adelie[yvar]

10.8.2. Alıştırma 10.2#

Bu bölümde bir uyum doğrusunun eğiminin örnekleme dağılımına yaklaşmak için yeniden örnekleme kullandık. Kesişimin örnekleme dağılımına da aynı biçimde yaklaşabiliriz:

  1. Yeniden örneklenmiş bir DataFrame alan, penguen ağırlığını yüzgeç uzunluğunun işlevi olarak en küçük karelerle uyduran ve kesişimi döndüren estimate_intercept işlevini yazın.

  2. İşlevi adelie veri kümesinin çok sayıda yeniden örneklenmiş sürümüyle çağırıp kesişimleri toplayın.

  3. Kesişimin örnekleme dağılımını plot_kde ile çizin.

  4. Standart hatayı ve %90 güven aralığını hesaplayın.

  5. Yeniden örneklemeden elde ettiğiniz standart hatanın RegressionResult içindeki intercept_stderr özniteliğiyle tutarlı olduğunu denetleyin; biraz daha küçük olabilir.

10.8.3. Alıştırma 10.3#

Bir kişinin Vücut Kütle İndeksi (VKİ), kilogram cinsinden ağırlığının metre cinsinden boyunun karesine bölümüdür. BRFSS veri kümesinde boyları santimetreden metreye çevirdikten sonra VKİ’yi şöyle hesaplayabiliriz.

heights_m = heights / 100
bmis = weights / heights_m**2

Bu tanımda boyların başka bir kuvveti değil karesi alınır; çünkü istatistiğin erken dönemlerinde ortalama ağırlığın yaklaşık olarak boyun karesiyle orantılı arttığı gözlenmiştir.

Bunun doğru olup olmadığını BRFSS verileri, en küçük kareler uyumu ve biraz matematikle sınayabiliriz. Ağırlığın, bilinmeyen \(a\) kuvvetine yükseltilmiş boyla orantılı olduğunu varsayın:

\[w = b h^a\]

Burada \(w\) ağırlık, \(h\) boy ve \(b\) bilinmeyen orantı sabitidir. İki tarafın logaritmasını alırsak:

\[\log w = \log b + a \log h\]

Dolayısıyla logaritmik dönüşümlü ağırlıkları logaritmik dönüşümlü boyların işlevi olarak en küçük karelerle uydurursak doğrunun eğimi bilinmeyen \(a\) üssünü tahmin eder.

Boy ve ağırlıkların logaritmalarını hesaplayın. İki dönüşümde de aynı olmak koşuluyla herhangi bir taban kullanabilirsiniz. En küçük kareler uyumunu hesaplayın. Eğim 2’ye yakın mı?

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ı