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.
11. Çoklu regresyon#
Önceki bölümdeki doğrusal en küçük kareler uyumu, regresyon örneğidir. Regresyon, yanıt değişkenleri ya da bağımlı değişkenler denen bir değişken kümesiyle açıklayıcı değişkenler ya da bağımsız değişkenler denen başka bir küme arasındaki ilişkiyi modelleme problemidir.
Önceki örneklerde bir yanıt ve bir açıklayıcı değişken vardı; buna basit regresyon denir. Bu bölümde yine tek yanıt değişkeniyle, fakat birden çok açıklayıcı değişken kullanan çoklu regresyona geçiyoruz. Birden çok yanıt değişkeni varsa çok değişkenli regresyon söz konusudur; bu kitapta onu ele almayacağız.
Bu notebook’u Colab’da çalıştırmak için buraya tıklayın.
11.1. StatsModels#
Önceki bölümde en küçük kareler uyumlarını hesaplamak için SciPy’nin linregress işlevini kullandık. Bu işlev basit regresyon yapar, çoklu regresyon yapmaz. Bunun için çeşitli regresyon ve analiz türleri sağlayan StatsModels paketini kullanacağız.
İlk örnek olarak penguen verilerini incelemeyi sürdüreceğiz.
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 yüklerken sütunlara boşluk içermeyen adlar vermek için aşağıdaki sözlüğü kullanacağız; böylece StatsModels ile kullanımları kolaylaşır.
columns = {
"Body Mass (g)": "mass",
"Flipper Length (mm)": "flipper_length",
"Culmen Length (mm)": "culmen_length",
"Culmen Depth (mm)": "culmen_depth",
}
Şimdi verileri yükleyebilir, kütlesi eksik satırları çıkarabilir ve sütunları yeniden adlandırabiliriz.
penguins = (
pd.read_csv("penguins_raw.csv")
.dropna(subset=["Body Mass (g)"])
.rename(columns=columns)
)
penguins.shape
(342, 17)
Veri kümesi üç penguen türünü içerir. Yalnızca Adélie penguenleriyle çalışacağız.
adelie = penguins.query('Species.str.startswith("Adelie")').copy()
len(adelie)
151
Önceki bölümde penguenlerin yüzgeç uzunlukları ile ağırlıkları arasında en küçük kareler uyumu hesapladık.
flipper_length = adelie["flipper_length"]
body_mass = adelie["mass"]
Anımsatmak için bunu linregress ile şöyle yapmıştık.
from scipy.stats import linregress
result_linregress = linregress(flipper_length, body_mass)
result_linregress.intercept, result_linregress.slope
(np.float64(-2535.8368022002524), np.float64(32.831689751150094))
StatsModels iki arayüz (API) sağlar. Yanıt ve açıklayıcı değişkenleri belirtmek için Patsy formül dilini kullanan “formül” API’sini kullanacağız. Aşağıdaki formül, yanıt değişkeni mass değerinin açıklayıcı flipper_length değişkeninin doğrusal işlevi olduğunu belirtir.
formula = "mass ~ flipper_length"
Bu formülü verilerle birlikte StatsModels içindeki ols işlevine verebiliriz.
import statsmodels.formula.api as smf
model = smf.ols(formula, data=adelie)
type(model)
statsmodels.regression.linear_model.OLS
ols adı “ordinary least squares”, yani olağan en küçük karelerin kısaltmasıdır; işlevin en yaygın varsayımlar altında en küçük kareler uyumu hesapladığını belirtir.
Sonuç, modeli temsil eden bir OLS nesnesidir. Model, değişkenler arasındaki ilişkinin basitleştirilmiş bir açıklamasıdır. Bu örnekte model doğrusaldır; yanıt değişkeninin açıklayıcı değişkenlerin doğrusal bir işlevi olduğunu varsayar.
fit yöntemi modeli verilere uydurur ve sonucu içeren bir RegressionResults nesnesi döndürür.
result_ols = model.fit()
# Teknik olarak bir RegressionResultsWrapper nesnesidir
type(result_ols)
statsmodels.regression.linear_model.RegressionResultsWrapper
RegressionResults çok fazla bilgi içerdiğinden thinkstats, şimdilik gereken bilgileri gösteren bir işlev sağlar.
from thinkstats import display_summary
display_summary(result_ols)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -2535.8368 | 964.798 | -2.628 | 0.009 | -4442.291 | -629.382 |
| flipper_length | 32.8317 | 5.076 | 6.468 | 0.000 | 22.801 | 42.862 |
| R-kare: | 0.2192 |
İlk sütunda modelin katsayıları olan kesişim ve eğim bulunur. Bunların linregress ile elde ettiğimiz katsayılarla aynı olduğunu doğrulayabiliriz.
result_linregress.intercept, result_linregress.slope
(np.float64(-2535.8368022002524), np.float64(32.831689751150094))
İkinci sütunda katsayıların standart hataları vardır; bunlar da linregress sonuçlarıyla aynıdır.
result_linregress.intercept_stderr, result_linregress.stderr,
(np.float64(964.7984274994061), np.float64(5.076138407990822))
Sonraki sütun p-değerlerini hesaplamakta kullanılan \(t\) istatistiklerini bildirir; p-değerleri P>|t| etiketli sonraki sütunda olduğundan bunları görmezden gelebiliriz. flipper_length p-değeri 0’a yuvarlanmıştır; şöyle gösterebiliriz.
result_ols.pvalues["flipper_length"]
np.float64(1.3432645947789658e-09)
linregress işlevinin aynı sonucu hesapladığını da doğrulayabiliriz.
result_linregress.pvalue
np.float64(1.3432645947790076e-09)
p-değeri çok küçüktür. Ağırlık ile yüzgeç uzunluğu arasında gerçekte ilişki olmasaydı tahmin edilen kadar büyük bir eğimi şans eseri görmemiz çok düşük olasılıklı olurdu.
[0.025 ve 0.975] etiketli son iki sütun, kesişim ile eğimin %95 güven aralıklarını bildirir. Eğimin %95 GA’sı [22,8; 42,9]’dur.
Son satır modelin yaklaşık 0,22 olan \(R^2\) değerini bildirir. Yalnızca ortalama ağırlık yerine yüzgeç uzunluğunu kullanarak ağırlığı öngörürsek MSE’yi yaklaşık %22 azaltabiliriz.
Basit korelasyondan elde edilen \(R^2\), korelasyon katsayısı \(r\)’nin karesidir. Bu nedenle ols ile hesaplanan rsquared değerini linregress ile hesaplanan rvalue karesiyle karşılaştırabiliriz.
result_ols.rsquared, result_linregress.rvalue**2
(np.float64(0.219212826468549), np.float64(0.2192128264685488))
Kayan nokta yaklaşımından doğan küçük fark dışında aynıdırlar.
Çoklu regresyona geçmeden önce açıklayıcı değişken olarak culmen_length kullanan bir basit regresyon daha hesaplayalım (culmen, gaganın üst sırtıdır).
formula = "mass ~ culmen_length"
result = smf.ols(formula, data=adelie).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 34.8830 | 458.439 | 0.076 | 0.939 | -870.998 | 940.764 |
| culmen_length | 94.4998 | 11.790 | 8.015 | 0.000 | 71.202 | 117.798 |
| R-kare: | 0.3013 |
Eğimin p-değeri yine çok küçüktür; kütle ile gaga sırtı uzunluğu arasında gerçekte ilişki olmasaydı bu kadar büyük bir eğimi şans eseri görmemiz düşük olasılıklı olurdu. Kesişimle ilişkili p-değerinin büyük olması sorun değildir; kesişimin sıfır olup olmadığıyla ilgilenmiyoruz. Bu modelde kesişimin sıfıra yakın olması yalnızca rastlantıdır.
Bu modelin \(R^2\) değeri yaklaşık 0,30’dur. Açıklayıcı değişken olarak gaga sırtı uzunlığını kullanırsak MSE azalması yüzgeç uzunluğuna göre biraz daha yüksektir (yüzgeç uzunluğunda 0,22). Şimdi ikisini birleştirelim.
11.2. Çoklu regresyona geçiş#
Kütlenin hem yüzgeç hem gaga sırtı uzunluğunun doğrusal işlevi olduğu çoklu regresyon modelinin Patsy formülü şöyledir.
formula = "mass ~ flipper_length + culmen_length"
Modeli verilere uydurmanın sonucu aşağıdadır.
result = smf.ols(formula, data=adelie).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -3573.0817 | 866.739 | -4.122 | 0.000 | -5285.864 | -1860.299 |
| flipper_length | 22.7024 | 4.742 | 4.787 | 0.000 | 13.331 | 32.074 |
| culmen_length | 76.3402 | 11.644 | 6.556 | 0.000 | 53.331 | 99.350 |
| R-kare: | 0.3949 |
Bu modelin üç katsayısı vardır: kesişim ve iki eğim. Yüzgeç uzunluğuyla ilişkili eğim 22,7’dir; gaga sırtı uzunluğu aynı kaldığında yüzgeci bir milimetre daha uzun olan penguenin 22,7 gram daha ağır olmasını bekleriz. Benzer biçimde yüzgeç uzunluğu sabitken gagası bir milimetre daha uzun olan penguenin 76,3 gram daha ağır olmasını bekleriz.
İki eğimin p-değeri de küçüktür; iki açıklayıcı değişkenin katkısının şans eseri oluşması düşük olasılıklıdır.
\(R^2\) 0,39’dur; yalnızca gaga sırtı uzunluğu (0,30) ya da yalnızca yüzgeç uzunluğu (0,22) kullanan modellerden yüksektir. İki değişkene dayanan öngörüler tek başına birine dayananlardan daha iyidir.
Ancak umduğumuz kadar iyi değildir. Yüzgeç uzunluğu MSE’yi %22, gaga sırtı uzunluğu %30 azaltıyorsa ikisi neden toplam %52 azaltmaz? Çünkü açıklayıcı değişkenler birbiriyle korelasyonludur.
from thinkstats import corrcoef
corrcoef(adelie, "flipper_length", "culmen_length")
np.float64(0.3257847151651593)
Yüzgeci daha uzun olan penguenin gagası da ortalama olarak daha uzundur. Açıklayıcı değişkenler birbirleri hakkında, dolayısıyla yanıt değişkeni hakkında aynı bilgilerin bir bölümünü taşır. Modele açıklayıcı değişken eklediğimizde \(R^2\) iyileşmesi yalnızca yeni değişkenin sağladığı yeni bilgiyi yansıtır.
Üçüncü değişken olarak gaga sırtı derinliğini eklediğimizde de aynı örüntüyü görürüz.
formula = "mass ~ flipper_length + culmen_length + culmen_depth"
result = smf.ols(formula, data=adelie).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -4341.3019 | 795.117 | -5.460 | 0.000 | -5912.639 | -2769.964 |
| flipper_length | 17.4215 | 4.385 | 3.973 | 0.000 | 8.756 | 26.087 |
| culmen_length | 55.3676 | 11.133 | 4.973 | 0.000 | 33.366 | 77.369 |
| culmen_depth | 140.8946 | 24.216 | 5.818 | 0.000 | 93.037 | 188.752 |
| R-kare: | 0.5082 |
Bu modelin dört katsayısı vardır. Tüm p-değerleri küçüktür; her açıklayıcı değişkenin katkısının şans eseri oluşması düşük olasılıklıdır. \(R^2\) yaklaşık 0,51’dir; iki değişkenli önceki modelden (0,39) ve tek değişkenli modellerden (0,22 ve 0,30) daha iyidir.
Ancak gaga sırtı derinliği diğer iki ölçümle korelasyonlu olduğundan ek iyileşme yine umulandan küçüktür.
[
corrcoef(adelie, "culmen_depth", "flipper_length"),
corrcoef(adelie, "culmen_depth", "culmen_length"),
]
[np.float64(0.3076201793966855), np.float64(0.3914916918358764)]
Bu örnek çoklu regresyonun yaygın kullanımını, daha iyi öngörüler için birden çok açıklayıcı değişkeni birleştirmeyi gösterir. Başka bir yaygın kullanım, bir değişken kümesinin katkısını diğer kümenin katkısını kontrol ederken nicelleştirmektir.
11.3. Kontrol değişkenleri#
Bölüm 4’te ilk bebeklerin ortalama olarak diğerlerinden daha hafif olduğunu gördük. Bölüm 9’da doğum ağırlığının annenin yaşıyla korelasyonlu olduğunu; daha yaşlı annelerin bebeklerinin ortalama olarak daha ağır olduğunu gördük.
Bu sonuçlar ilişkili olabilir. İlk kez anne olanlar diğer annelerden daha gençse —ki olasıdır— bebeklerinin neden daha hafif olduğunu bu açıklayabilir. Anne yaşını kontrol ederken ilk bebeklerle diğerleri arasındaki doğum ağırlığı farkını tahmin ederek bu varsayımı çoklu regresyonla sınayabiliriz.
NSFG verilerini indirme yönergeleri bu bölümün notebook’unda yer alır.
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, canlı doğumları seçmek ve ilk bebekler ile diğerlerini gruplamak için get_nsfg_groups kullanabiliriz.
from nsfg import get_nsfg_groups
live, firsts, others = get_nsfg_groups()
Geçerli doğum ağırlığı, doğum sırası ve anne yaşı bulunan satırları dropna ile seçeceğiz.
valid = live.dropna(subset=["agepreg", "birthord", "totalwgt_lb"]).copy()
Doğum ağırlığının yaşla korelasyonlu olduğunu doğrulamak ve doğrusal ilişki varsayımıyla eğimi tahmin etmek için StatsModels kullanabiliriz.
formula = "totalwgt_lb ~ agepreg"
result_age = smf.ols(formula, data=valid).fit()
display_summary(result_age)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 6.8304 | 0.068 | 100.470 | 0.000 | 6.697 | 6.964 |
| agepreg | 0.0175 | 0.003 | 6.559 | 0.000 | 0.012 | 0.023 |
| R-kare: | 0.004738 |
Eğim küçüktür: yılda yalnızca 0,0175 pound. Yaşları on yıl farklı iki annenin bebeklerinin ağırlıklarının 0,175 pound farklı olmasını bekleriz. Ancak p-değeri küçüktür; gerçekte ilişki olmasa bu eğim düşük olasılıklı olurdu.
\(R^2\) de küçüktür; anne yaşı öngörü değişkeni olarak pek yararlı değildir. Annenin yaşını bilmek bebeğin ağırlığını öngörme yeteneğimizi neredeyse hiç iyileştirmez.
Küçük p-değeri ile küçük \(R^2\) birlikteliği sık karıştırılır; ilişki istatistiksel olarak anlamlıysa öngörücü olması gerektiği düşünülür. Bu örnek çelişki olmadığını gösterir: ilişki anlamlı, ancak öngörü için yararsız olabilir. Görselleştirdiğimizde nedenini göreceğiz. Önce ilgili sütunları seçelim.
totalwgt = valid["totalwgt_lb"]
agepreg = valid["agepreg"]
Uyum doğrusunu hesaplamak için result_age içinden kesişim ve eğimi çıkarabilirdik; ancak RegressionResults nesnesinin predict yöntemini kullanabiliriz. Önce bir agepreg değer aralığı hesaplayacağız.
agepreg_range = np.linspace(agepreg.min(), agepreg.max())
predict kullanmak için açıklayıcı değişkenlerin değerlerini bir DataFrame içine koymalıyız.
df = pd.DataFrame({"agepreg": agepreg_range})
DataFrame sütunlarının adları açıklayıcı değişkenlerle aynı olmalıdır. Ardından bunu predict işlevine verebiliriz.
fit_ys = result_age.predict(df)
Sonuç öngörülen değerleri içeren bir Series’dir. Verilerin saçılım grafiğiyle birlikte şöyle görünür.
plt.scatter(agepreg, totalwgt, marker=".", alpha=0.1, s=5)
plt.plot(agepreg_range, fit_ys, color="C1", label="doğrusal model")
decorate(xlabel="Anne yaşı", ylabel="Doğum ağırlığı (libre)")
Uyum doğrusunun eğimi küçük olduğundan en genç ve en yaşlı annelerin beklenen doğum ağırlıkları arasındaki farkı güçlükle görürüz. Her anne yaşında doğum ağırlığındaki değişim çok daha büyüktür.
Şimdi ilk bebeklerin diğerlerinden daha hafif olduğunu StatsModels ile doğrulayacağız. İlk bebeklerde True, diğerlerinde False olan bir Boole Series oluşturup is_first adlı sütun olarak ekleyeceğiz.
valid["is_first"] = valid["birthord"] == 1
from thinkstats import value_counts
# sonuçları denetle
value_counts(valid["is_first"])
is_first
False 4675
True 4363
Name: count, dtype: int64
Yanıt değişkeni doğum ağırlığı, açıklayıcı değişkeni is_first olan modelin formülü şöyledir. Patsy dilinde C ile değişken adının çevresindeki parantez, değişkenin kategorik olduğunu; doğum ağırlığı gibi ölçüm değil “ilk bebek” gibi kategori temsil ettiğini belirtir.
formula = "totalwgt_lb ~ C(is_first)"
Şimdi her zamanki gibi modeli uydurup sonuçları gösterebiliriz.
result_first = smf.ols(formula, data=valid).fit()
display_summary(result_first)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 7.3259 | 0.021 | 356.007 | 0.000 | 7.286 | 7.366 |
| C(is_first)[T.True] | -0.1248 | 0.030 | -4.212 | 0.000 | -0.183 | -0.067 |
| R-kare: | 0.00196 |
Sonuçlardaki C(is_first)[T.True], is_first değişkeninin kategorik olduğunu ve katsayının True değeriyle ilişkili olduğunu belirtir. True önündeki T, kontrollü deney dilindeki “treatment”ı gösterir; ilk bebekler müdahale, diğerleri referans grubudur. Bu tanımlar keyfîdir; tersini de seçebilirdik. Sonuçları yorumlamak için hangisinin hangisi olduğunu bilmeliyiz.
Kesişim yaklaşık 7,3’tür; referans grubunun ortalama ağırlığı 7,3 pounddur. is_first katsayısı -0,12’dir; müdahale grubunun, yani ilk bebeklerin ortalaması 0,12 pound daha düşüktür. İki sonucu doğrudan hesaplayarak denetleyebiliriz.
others["totalwgt_lb"].mean()
np.float64(7.325855614973262)
diff_weight = firsts["totalwgt_lb"].mean() - others["totalwgt_lb"].mean()
diff_weight
np.float64(-0.12476118453549034)
StatsModels bu katsayıların yanı sıra p-değerlerini, güven aralıklarını ve \(R^2\)’yi de hesaplar. İlk bebeklerle ilişkili p-değeri küçüktür; gruplar arasındaki fark istatistiksel olarak anlamlıdır. \(R^2\) küçüktür; bebeğin ağırlığını tahmin ederken ilk bebek olup olmadığını bilmek fazla yardımcı olmaz.
Şimdi farkın anne yaşındaki farktan kaynaklanmasının makul olup olmadığına bakalım. İlk kez anne olanlar ortalama 3,6 yıl daha gençtir.
diff_age = firsts["agepreg"].mean() - others["agepreg"].mean()
diff_age
np.float64(-3.5864347661500275)
Doğum ağırlığının yaşın işlevi olarak eğimi yılda 0,0175 pounddur.
slope = result_age.params["agepreg"]
slope
np.float64(0.0174538514718027)
Eğimi yaş farkıyla çarparsak anne yaşından kaynaklanması beklenen doğum ağırlığı farkını elde ederiz.
slope * diff_age
np.float64(-0.06259709972169203)
Sonuç 0,063 pounddur; gözlenen farkın yaklaşık yarısıdır. Dolayısıyla gözlenen doğum ağırlığı farkı anne yaşı farkıyla kısmen açıklanabilir.
Çoklu regresyonla anne yaşı ve ilk bebek katsayılarını aynı anda tahmin edebiliriz.
formula = "totalwgt_lb ~ agepreg + C(is_first)"
result = smf.ols(formula, data=valid).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 6.9142 | 0.078 | 89.073 | 0.000 | 6.762 | 7.066 |
| C(is_first)[T.True] | -0.0698 | 0.031 | -2.236 | 0.025 | -0.131 | -0.009 |
| agepreg | 0.0154 | 0.003 | 5.499 | 0.000 | 0.010 | 0.021 |
| R-kare: | 0.005289 |
is_first katsayısı -0,0698’dir; anne yaşı farkı hesaba katıldığında ilk bebekler ortalama 0,0698 pound daha hafiftir. Bu, anne yaşını hesaba katmadan elde edilen farkın yaklaşık yarısıdır.
p-değeri 0,025’tir; yine anlamlı kabul edilir, ancak bu büyüklükteki farkın şans eseri oluşmasını eleyemediğimiz sınır bölgesindedir.
Model anne yaşından kaynaklanan ağırlık farkını hesaba kattığı için anne yaşını kontrol ettiğini söyleyebiliriz. Ancak ağırlık ile anne yaşı ilişkisinin doğrusal olduğunu varsayar. Bunun doğru olup olmadığına bakalım.
11.4. Doğrusal olmayan ilişkiler#
agepreg katkısının doğrusal olmayıp olmadığını denetlemek için veri kümesine agepreg değerlerinin karelerini içeren yeni bir sütun ekleyebiliriz.
valid["agepreg2"] = valid["agepreg"] ** 2
Şimdi hem doğrusal hem karesel ilişki içeren bir model tanımlayabiliriz.
formula = "totalwgt_lb ~ agepreg + agepreg2"
Modeli her zamanki gibi uydurabiliriz.
result_age2 = smf.ols(formula, data=valid).fit()
display_summary(result_age2)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 5.5720 | 0.275 | 20.226 | 0.000 | 5.032 | 6.112 |
| agepreg | 0.1186 | 0.022 | 5.485 | 0.000 | 0.076 | 0.161 |
| agepreg2 | -0.0019 | 0.000 | -4.714 | 0.000 | -0.003 | -0.001 |
| R-kare: | 0.00718 |
Karesel agepreg2 teriminin p-değeri çok küçüktür; doğum ağırlığı hakkında şanstan beklenenden fazla bilgi sağladığını düşündürür. Modelin \(R^2\) değeri 0,0072’dir; doğrusal modelin 0,0047 değerinden yüksektir.
agepreg ve agepreg2 katsayılarını tahmin ederek verilere aslında bir parabol uyduruyoruz. Bunu görmek için RegressionResults ile bir anne yaşı aralığı için öngörüler üretebiliriz.
Önce agepreg_range temelinde agepreg ve agepreg2 sütunları içeren geçici bir DataFrame oluşturacağız.
df = pd.DataFrame({"agepreg": agepreg_range})
df["agepreg2"] = df["agepreg"] ** 2
Şimdi DataFrame’i predict yöntemine verip bir öngörü Series’i elde edebiliriz.
fit_ys = result_age2.predict(df)
Verilerin saçılım grafiğiyle birlikte uydurulan parabol şöyle görünür.
plt.scatter(agepreg, totalwgt, marker=".", alpha=0.1, s=5)
plt.plot(agepreg_range, fit_ys, color="C1", label="karesel model")
decorate(xlabel="Anne yaşı", ylabel="Doğum ağırlığı (libre)")
Eğrilik hafiftir; ancak en genç ve en yaşlı annelerin bebeklerinin daha hafif, ortadakilerin daha ağır olduğunu düşündürür.
Karesel model ilişkiyi doğrusal modelden daha iyi yakalar; anne yaşından kaynaklanan doğum ağırlığı farkını daha etkili biçimde açıklayabilir. is_first değişkenini karesel modele eklediğimizde ne olduğuna bakalım.
formula = "totalwgt_lb ~ agepreg + agepreg2 + C(is_first)"
result = smf.ols(formula, data=valid).fit()
display_summary(result)
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 5.6923 | 0.286 | 19.937 | 0.000 | 5.133 | 6.252 |
| C(is_first)[T.True] | -0.0504 | 0.031 | -1.602 | 0.109 | -0.112 | 0.011 |
| agepreg | 0.1124 | 0.022 | 5.113 | 0.000 | 0.069 | 0.155 |
| agepreg2 | -0.0018 | 0.000 | -4.447 | 0.000 | -0.003 | -0.001 |
| R-kare: | 0.007462 |
Anne yaşını daha etkili kontrol ettiğimizde ilk bebeklerle diğerleri arasındaki tahmini fark 0,0504 pounddur; yalnızca doğrusal modeldeki 0,0698’den küçüktür. is_first p-değeri 0,109’dur; gruplar arasında kalan farkın şanstan kaynaklanması makuldür.
Doğum ağırlığı farkının en azından kısmen, belki tamamen anne yaşı farkıyla açıklandığı sonucuna varabiliriz.
11.5. Lojistik regresyon#
Doğrusal regresyonda yanıt değişkeninin beklenen değeri açıklayıcı değişkenlerin ağırlıklı toplamı ile kesişimden oluşur. Bu model doğum ağırlığı gibi sürekli yanıtlar için uygundur; sayım ya da kategori gibi kesikli yanıtlar için uygun değildir.
Bu yanıtlar için genelleştirilmiş doğrusal modelleri (GLM) kullanabiliriz:
Yanıt bir sayımsa Poisson regresyonu.
İki kategoriliyse lojistik regresyon.
İkiden fazla kategoriliyse çok terimli lojistik regresyon.
Kategoriler sıralanabiliyorsa sıralı lojistik regresyon.
Kitapta yalnızca en yaygın olan lojistik regresyonu ele alacağız. Penguen verileriyle ağırlık ve diğer ölçümlerden penguenin erkek mi dişi mi olduğunu belirlemeye çalışacağız.
StatsModels, lojistik regresyon yapan logit işlevini sağlar; adı tanımdaki matematiksel işlevden gelir. Kullanmadan önce yanıt değerlerini 0 ve 1 olacak biçimde dönüştürmeliyiz.
adelie["y"] = (adelie["Sex"] == "MALE").astype(int)
adelie["y"].value_counts()
y
0 78
1 73
Name: count, dtype: int64
Yanıt değişkeni y, açıklayıcı değişkeni mass olan basit bir modelle başlayacağız. Modeli şöyle oluşturup uydururuz; disp=False argümanı uyum süreci iletilerini gizler.
model = smf.logit("y ~ mass", data=adelie)
result = model.fit(disp=False)
Sonuçlar şöyledir.
display_summary(result)
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -25.9871 | 4.221 | -6.156 | 0.000 | -34.261 | -17.713 |
| mass | 0.0070 | 0.001 | 6.138 | 0.000 | 0.005 | 0.009 |
| Sözde R-kare: | 0.5264 |
Belirleme katsayısı \(R^2\) lojistik regresyona uygulanmaz; “sözde \(R^2\)” olarak kullanılan birkaç seçenek vardır. Bu modelinki yaklaşık 0,526’dır; tek başına fazla anlam taşımaz, ancak modelleri karşılaştırmakta kullanacağız.
mass katsayısı pozitiftir; daha ağır penguenlerin erkek olma olasılığı daha yüksektir. Katsayıların geri kalanını yorumlamak güçtür; öngörüleri çizerek modeli daha iyi anlayabiliriz. Bir mass değer aralığı içeren DataFrame oluşturup predict ile öngörü Series’i hesaplayacağız.
mass = adelie["mass"]
mass_range = np.linspace(mass.min(), mass.max())
df = pd.DataFrame({"mass": mass_range})
fit_ys = result.predict(df)
Her öngörülen değer, ağırlığının işlevi olarak penguenin erkek olma olasılığıdır. Şöyle görünür.
plt.plot(mass_range, fit_ys)
decorate(xlabel="Kütle (g)", ylabel="Olasılık(erkek)")
En hafif penguenler neredeyse kesinlikle dişi, en ağırları büyük olasılıkla erkektir. Ortada, 3750 gramlık bir penguenin erkek ya da dişi olma olasılığı yaklaşık eşittir.
Diğer ölçümleri açıklayıcı değişken olarak eklersek ne olduğuna bakalım.
formula = "y ~ mass + flipper_length + culmen_length + culmen_depth"
model = smf.logit(formula, data=adelie)
result = model.fit(disp=False)
display_summary(result)
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -60.6075 | 13.793 | -4.394 | 0.000 | -87.642 | -33.573 |
| mass | 0.0059 | 0.001 | 4.153 | 0.000 | 0.003 | 0.009 |
| flipper_length | -0.0209 | 0.052 | -0.403 | 0.687 | -0.123 | 0.081 |
| culmen_length | 0.6208 | 0.176 | 3.536 | 0.000 | 0.277 | 0.965 |
| culmen_depth | 1.0111 | 0.349 | 2.896 | 0.004 | 0.327 | 1.695 |
| Sözde R-kare: | 0.6622 |
Bu modelin sözde \(R^2\) değeri 0,662’dir; önceki 0,526’dan yüksektir. Ek ölçümler cinsiyetleri ayıran ek bilgi içerir.
Gaga sırtı uzunluğu ile derinliğinin p-değerleri küçüktür; şanstan beklenenden fazla bilgi sağlarlar. Yüzgeç uzunluğunun p-değeri büyüktür; ağırlık ile gaga boyutları biliniyorsa yüzgeç uzunluğunun ek bilgi sağlamadığını düşündürür.
Öngörüleri incelemek için bir kütle dizisi ile belirli bir culmen_length değeri alan aşağıdaki işlevi kullanacağız. Diğer ölçümleri ortalamalarına ayarlar, kütlenin işlevi olarak olasılıkları hesaplar ve çizer.
def plot_predictions(mass_range, culmen_length, **options):
"""Kütlenin işlevi olarak öngörülen olasılıkları çizer."""
df = pd.DataFrame({"mass": mass_range})
df["flipper_length"] = adelie["flipper_length"].mean()
df["culmen_length"] = culmen_length
df["culmen_depth"] = adelie["culmen_depth"].mean()
fit_ys = result.predict(df)
plt.plot(mass_range, fit_ys, **options)
Ortalamanın bir standart sapma üstü, ortalama ve bir standart sapma altı olmak üzere üç culmen_length değerinin sonuçları şöyledir.
culmen_length = adelie["culmen_length"]
m, s = culmen_length.mean(), culmen_length.std()
plot_predictions(mass_range, m + s, ls="--", label="Ortalamanın üzerinde gaga sırtı uzunluğu")
plot_predictions(mass_range, m, alpha=0.5, label="Ortalama gaga sırtı uzunluğu")
plot_predictions(mass_range, m - s, ls=":", label="Ortalamanın altında gaga sırtı uzunluğu")
decorate(xlabel="Kütle (g)", ylabel="Olasılık(erkek)")
Basit modelde gördüğümüz gibi ağır penguenlerin erkek olma olasılığı daha yüksektir. Ayrıca herhangi bir ağırlıkta gagası daha uzun penguenin erkek olma olasılığı daha yüksektir.
Bu model göründüğünden daha yararlıdır; verilerin toplandığı özgün araştırma makalesindeki modellere benzer. Araştırmanın ana konusu, erkek ve dişi bedenlerin farklılık derecesi olan eşeysel dimorfizmdir. Dimorfizmi nicelleştirmenin bir yolu ölçümlerle erkek ve dişileri sınıflandırmaktır. Dimorfizmi yüksek türlerde sınıflandırmaların daha doğru olmasını bekleriz.
Özgün makale: 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
Yöntemi sınamak için aynı modeli başka bir türde deneyelim. Adélie penguenlerine ek olarak veri kümesi 123 Gentoo pengueninin ölçümlerini içerir. Bunları seçmek için aşağıdaki işlevi kullanacağız.
def get_species(penguins, species):
df = penguins.query(f'Species.str.startswith("{species}")').copy()
df["y"] = (df["Sex"] == "MALE").astype(int)
return df
gentoo = get_species(penguins, "Gentoo")
len(gentoo)
123
Lojistik regresyon modelinin sonuçları şöyledir.
formula = "y ~ mass + flipper_length + culmen_length + culmen_depth"
model = smf.logit(formula, data=gentoo)
result = model.fit(disp=False)
display_summary(result)
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -173.9123 | 62.326 | -2.790 | 0.005 | -296.069 | -51.756 |
| mass | 0.0105 | 0.004 | 2.948 | 0.003 | 0.004 | 0.017 |
| flipper_length | 0.2839 | 0.183 | 1.549 | 0.121 | -0.075 | 0.643 |
| culmen_length | 0.2734 | 0.285 | 0.958 | 0.338 | -0.286 | 0.833 |
| culmen_depth | 3.0843 | 1.291 | 2.389 | 0.017 | 0.554 | 5.614 |
| Sözde R-kare: | 0.848 |
Sözde \(R^2\) değeri 0,848’dir; Adélie penguenlerinde elde ettiğimiz 0,662’den yüksektir. Gentoo penguenleri fiziksel ölçümlerle daha doğru sınıflandırılabilir; bu da daha dimorfik olduklarını düşündürür.
11.6. Sözlük#
regresyon: Verilere bir model uyduran katsayıları tahmin etme yöntemi.
yanıt değişkenleri (response variables): Regresyon modelinin öngörmeye çalıştığı, bağımlı değişkenler de denen değişkenler.
açıklayıcı değişkenler (explanatory variables): Modelin yanıtları öngörmek için kullandığı, bağımsız değişkenler de denen değişkenler.
basit regresyon: Bir yanıt ve bir açıklayıcı değişkenli regresyon.
çoklu regresyon: Tek yanıt, birden çok açıklayıcı değişkenli regresyon.
katsayılar: Regresyon modelindeki kesişim ve açıklayıcı değişkenlerin tahmini eğimleri.
kategorik değişken: Genellikle sayısal olmayan, kesikli değer kümesinden birini alabilen değişken.
kontrol değişkeni: Açıklayıcı değişkenin doğrudan etkisini dolaylı etkiden ayırmak için regresyona eklenen değişken.
genelleştirilmiş doğrusal modeller: Açıklayıcı ve yanıt değişkenleri arasındaki farklı matematiksel ilişkilere dayanan regresyon modelleri.
lojistik regresyon: Yanıt değişkeninin yalnızca iki olası değeri olduğunda kullanılan genelleştirilmiş doğrusal model.
11.7. Alıştırmalar#
11.7.1. Alıştırma 11.1#
Erkek bebekler kız bebeklerden daha mı ağırdır? NSFG verilerini yeniden kullanacağız.
Yanıt değişkeni totalwgt_lb, kategorik açıklayıcı değişkeni babysex olan doğrusal regresyon modeli uydurun; bu değişken erkeklerde 1, kızlarda 2’dir. Tahmini ağırlık farkı nedir? İstatistiksel olarak anlamlı mı? Annenin yaşını kontrol ederseniz anne yaşı gözlenen farkın bir bölümünü ya da tamamını açıklar mı?
11.7.2. Alıştırma 11.2#
Trivers–Willard hipotezi, birçok memelide cinsiyet oranının annenin yaşı, bedeni, sağlığı ve sosyal statüsü gibi “anne koşullarına” bağlı olduğunu öne sürer. Bazı insan çalışmaları bu etkiyi göstermiştir, ancak sonuçlar karışıktır.
Bkz. https://en.wikipedia.org/wiki/Trivers-Willard_hypothesis.
Annenin yaşı ile erkek çocuk doğurma olasılığı arasında ilişki olup olmadığına bakalım. Bebeğin cinsiyetini yanıt, annenin yaşını açıklayıcı değişken alan lojistik regresyon modeli uydurun. Yaşlı annelerin erkek çocuk doğurma olasılığı daha mı yüksek, daha mı düşük? Anne yaşının karesel modelini kullanırsanız ne olur?
babysex değişkenini lojistik regresyonda yanıt olarak kullanmak için erkeklerde 1, kızlarda 0 olacak biçimde yeniden kodlayacağız.
valid["y"] = (valid["babysex"] == 1).astype(int)
11.7.3. Alıştırma 11.3#
Adélie penguenleri için ağırlığı flipper_length, culmen_depth ve kategorik Sex değişkenlerinin işlevi olarak öngören doğrusal regresyon modeli uydurun. Yüzgeç uzunluğu ile gaga sırtı derinliğini kontrol ettiğimizde erkekler ne kadar daha ağırdır? culmen_depth ortalamasına sabitlenmişken erkek ve dişi penguenler için bir yüzgeç uzunluğu aralığında öngörüler üretip çizin.
11.7.4. Alıştırma 11.4#
Sakallı penguenlerin, modelin sözde \(R^2\) değeriyle nicelleştirildiğinde veri kümesindeki diğer türlerden daha mı az, daha mı çok dimorfik olduğuna bakalım. Sakallı penguenleri get_species ile seçin; cinsiyeti yanıt, dört ölçümü açıklayıcı değişken alan lojistik regresyon modeli uydurun. Sözde \(R^2\) diğer modellerle nasıl karşılaştırılıyor?
chinstrap = get_species(penguins, "Chinstrap")
len(chinstrap)
68
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ı