Simpson Paradoksu: Palmer Penguins Örneği

Bu çalışmanın amacı Simpson Paradoksunu, palmerpenguins paketindeki penguins veri seti üzerinde incelemektir.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(palmerpenguins) #verinin olduğu kütüphaneyi çağırma
## 
## Attaching package: 'palmerpenguins'
## The following objects are masked from 'package:datasets':
## 
##     penguins, penguins_raw

1. Veri Setini Düzenleme

1.1. Veri setindeki değişken isimlerini değiştirme

İlk adımda veri setindeki değişkenlerin isimlerini Türkçe yapacağız.

penguenler0 <- penguins

penguenler <- penguenler0 %>% #r'da pipe operatörü fonksiyonları birbirine bağlar.
  rename( tur = species, #rename kullanarak sütunların adları değiştirilebilir
          ada = island,
          gaga_uzunlugu = bill_length_mm,
          gaga_derinligi = bill_depth_mm,
          yuzgec_uzunlugu = flipper_length_mm,
          vucut_kutlesi = body_mass_g,
          cinsiyet = sex,
          yil = year
          )
glimpse(penguenler) #veri setine bakmak için glimpse fonksiyonu kullanırız
## Rows: 344
## Columns: 8
## $ tur             <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adelie…
## $ ada             <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgersen,…
## $ gaga_uzunlugu   <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, 42…
## $ gaga_derinligi  <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, 20…
## $ yuzgec_uzunlugu <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186, …
## $ vucut_kutlesi   <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, 42…
## $ cinsiyet        <fct> male, female, female, NA, female, male, female, male, …
## $ yil             <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, …

1.2 Eksik verileri bulma ve silme

1.2.1 Eksik verileri bulmadan önceki veri setinin boyutlarını belirleme:

onceki_boyut <- c(satir = nrow(penguenler), sutun = ncol(penguenler)) #c fonksiyonu verilen değeri tek. bir vektörde birleştirir. burada satır ve sütun sayılarınıı tek bir vektörde vermiş olduk.
print(onceki_boyut) #print çıktıları basmamızı, yazdırmamızı sağlar.
## satir sutun 
##   344     8

Eksik verileri çıkarmadan önce veri setimizde 344 satır ve 8 sütun bulunmakta.

1.2.2. Eksik verileri silme:

penguenler_temiz <- na.omit(penguenler) #na.omit veri setindeki eksik değerlerin olduğu tüm satırları siler.
print(penguenler)
## # A tibble: 344 × 8
##    tur    ada       gaga_uzunlugu gaga_derinligi yuzgec_uzunlugu vucut_kutlesi
##    <fct>  <fct>             <dbl>          <dbl>           <int>         <int>
##  1 Adelie Torgersen          39.1           18.7             181          3750
##  2 Adelie Torgersen          39.5           17.4             186          3800
##  3 Adelie Torgersen          40.3           18               195          3250
##  4 Adelie Torgersen          NA             NA                NA            NA
##  5 Adelie Torgersen          36.7           19.3             193          3450
##  6 Adelie Torgersen          39.3           20.6             190          3650
##  7 Adelie Torgersen          38.9           17.8             181          3625
##  8 Adelie Torgersen          39.2           19.6             195          4675
##  9 Adelie Torgersen          34.1           18.1             193          3475
## 10 Adelie Torgersen          42             20.2             190          4250
## # ℹ 334 more rows
## # ℹ 2 more variables: cinsiyet <fct>, yil <int>

1.2.3 Eksik verileri sildikten sonraki veri setinin boyutları:

sonraki_boyut <- c(satir = nrow(penguenler_temiz), sutun = ncol(penguenler_temiz)) #nrow ve ncol satır ve sütunların toplam sayısını verir.
print(sonraki_boyut)
## satir sutun 
##   333     8

Eksik verileri temizledikten sonra 333 satır ve 8 sütun kalmıştır. 11 satır silinmiştir.

1.3. BMI değişkenini oluşturma ve betimsel istatistikler

\[\mathrm{BMI} = \frac{\text{Vücut kütlesi}}{\text{Yüzgeç uzunluğu}} \] olacak şekilde yeni bir değişken tanımlayacağız

penguenler_temiz <- penguenler_temiz %>% mutate(bmi = vucut_kutlesi/ yuzgec_uzunlugu) #mutate yeni sütun ekler veya veri setindeki bir sütunu değiştirmemizi sağlar. Burada bmi adında yeni bir sütun ekledik.

glimpse(penguenler_temiz)
## Rows: 333
## Columns: 9
## $ tur             <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adelie…
## $ ada             <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgersen,…
## $ gaga_uzunlugu   <dbl> 39.1, 39.5, 40.3, 36.7, 39.3, 38.9, 39.2, 41.1, 38.6, …
## $ gaga_derinligi  <dbl> 18.7, 17.4, 18.0, 19.3, 20.6, 17.8, 19.6, 17.6, 21.2, …
## $ yuzgec_uzunlugu <int> 181, 186, 195, 193, 190, 181, 195, 182, 191, 198, 185,…
## $ vucut_kutlesi   <int> 3750, 3800, 3250, 3450, 3650, 3625, 4675, 3200, 3800, …
## $ cinsiyet        <fct> male, female, female, female, male, female, male, fema…
## $ yil             <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, …
## $ bmi             <dbl> 20.71823, 20.43011, 16.66667, 17.87565, 19.21053, 20.0…

Veri setimize bmi adında yeni bir sütun eklenmiştir. Veri setimizde toplamda 333 satır ve 9 sütun bulunmaktadır.

BMI sütununun ortalama, standart sapma, minimum ve maksimum değerlerini bulacağız.

bmi_betimsel <- penguenler_temiz %>% 
  group_by(tur) %>% #group_by veri setini gruplara ayırmamızı sağlar. Biz burada türe göre grupladık ve bundan sonraki tüm işlemleri türe göre yaptırmış olduk.
  summarise(n = n(), ortalama = mean(bmi), standart_sapma= sd(bmi), minimum= min(bmi), maksimum= max(bmi)) #summarise veri çerçevesini özetlememizi sağlar. Özet istatistik çıkarır.

print(bmi_betimsel)
## # A tibble: 3 × 6
##   tur           n ortalama standart_sapma minimum maksimum
##   <fct>     <int>    <dbl>          <dbl>   <dbl>    <dbl>
## 1 Adelie      146     19.5           2.18    15.2     25.3
## 2 Chinstrap    68     19.0           1.60    14.1     22.9
## 3 Gentoo      119     23.4           1.88    19.0     28.5

Kutu grafiği:

ggplot(penguenler_temiz, aes(x = tur, y = bmi, fill = tur)) + #ggplot grafik oluşturmamızı sağlar. ggplotta yazdığımız fonksiyonları pipe operatörüyle değil, + ile bağlamamız gerekir.
  geom_boxplot() + #geom_ hangi grafik türünü seçeceğimizi belirtir.
  theme_minimal() + #theme_ grafiğin temasını ayarlamamızı sağlar.
  labs(title = "Türlere Göre BMI", x = "Tür", y = "BMI") #labs grafiğin başlığını, x ve y eksenlerinin adlarını yazmamızı sağlar.

Betimsel istatistiklere baktığımızda türlere göre ortalama BMI değerlerinde farklılıklar görülmektedir. En yüksek BMI Gentoo’da, sonra Adelie’de, en düşük Chinstrap’dadır. Adelie ile Chinstrap’ın ortalamaları birbirine oldukça yakındır. Grafik de betimsel istatistik sonuçlarını destekler niteliktedir.

2. Görev 1: Toplam Düzeyde Analiz

Eksik gözlemleri çıkarılmış veri seti üzerinde gaga uzunluğu ve gaga derinliği arasındaki ilişkiyi inceleyeceğiz ve grafiğe bir regresyon doğrusu ekleyeceğiz. Sonrasında basit regresyon denklemi kuracağız.

ggplot(penguenler_temiz, aes(x = gaga_uzunlugu, y = gaga_derinligi, color = tur)) +
geom_point() + 
geom_smooth(method = "lm", se = TRUE, color = "black") + #geom_smooth grafikteki noktaların uyum çizgisini çizmemizi sağlar. basit regresyonu grafikte geom_smooth kullanarak çizeriz.
theme_linedraw()+
labs(title = "Gaga Uzunluğu ve Gaga Derinliği Arasındaki İlişki",
x = "Gaga uzunluğu", y = "Gaga derinliği")
## `geom_smooth()` using formula = 'y ~ x'

Basit regresyon denklemi:

reg_penguenler <- lm(gaga_derinligi ~ gaga_uzunlugu, data = penguenler_temiz) #lm linear model demektir. doğrusal regresyon kurmamızı sağlar.
summary(reg_penguenler)
## 
## Call:
## lm(formula = gaga_derinligi ~ gaga_uzunlugu, data = penguenler_temiz)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1548 -1.4291  0.0122  1.3994  4.5004 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   20.78665    0.85417  24.335  < 2e-16 ***
## gaga_uzunlugu -0.08233    0.01927  -4.273 2.53e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.92 on 331 degrees of freedom
## Multiple R-squared:  0.05227,    Adjusted R-squared:  0.04941 
## F-statistic: 18.26 on 1 and 331 DF,  p-value: 2.528e-05

Grafiğe göre gaga uzunluğu arttıkça gaga derinliği azalmaktadır. Regresyon eğrisi aşağı yönlü olduğu için gaga uzunluğu ve gaga derinliği arasındaki ilişki negatif yönlüdür. Regresyon denkleminin sonuçlarında çıkan gaga_uzunlugu satırındaki estimate değeri eğimi vermektedir. Değer -0.08233 çıkmış ve p<.05 olduğundan bu değer istatistiksel olarak anlamlıdır. Toplam düzeyde gaga uzunluğu ve gaga derinliği arasındaki ilişki negatiftir.

3. Görev 2: Tür Bazında Analiz

3.1. Verilerin türlere göre ayrılması

adelie    <- filter(penguenler_temiz, tur == "Adelie") #filter veri çerçevesindeki satırları seçmemizi sağlar. Burada tür sütunundaki sadece adelie yazan satırları seçtik.
chinstrap <- filter(penguenler_temiz, tur == "Chinstrap")
gentoo    <- filter(penguenler_temiz, tur == "Gentoo")

3.2. Her tür için aynı scatterplot oluşturma

ggplot(penguenler_temiz, aes(gaga_uzunlugu, gaga_derinligi, color = tur)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  theme_minimal() +
  labs(title = "Tür Bazında Gaga Uzunluğu ve Gaga Derinliği",
       x = "Gaga uzunluğu", y = "Gaga derinliği", color = "Tür")
## `geom_smooth()` using formula = 'y ~ x'

3.3. Tür bazında basit doğrusal regresyon

reg_adelie <- lm(gaga_derinligi ~ gaga_uzunlugu, data = penguenler_temiz %>% 
                   filter(tur == "Adelie"))
summary(reg_adelie)
## 
## Call:
## lm(formula = gaga_derinligi ~ gaga_uzunlugu, data = penguenler_temiz %>% 
##     filter(tur == "Adelie"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1487 -0.7926 -0.0842  0.5550  3.4990 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   11.48771    1.37010   8.385 4.23e-14 ***
## gaga_uzunlugu  0.17668    0.03521   5.018 1.51e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.129 on 144 degrees of freedom
## Multiple R-squared:  0.1489, Adjusted R-squared:  0.1429 
## F-statistic: 25.18 on 1 and 144 DF,  p-value: 1.515e-06
reg_chinstrap <- lm(gaga_derinligi ~ gaga_uzunlugu, data = penguenler_temiz %>% 
                   filter(tur == "Chinstrap"))
summary(reg_chinstrap)
## 
## Call:
## lm(formula = gaga_derinligi ~ gaga_uzunlugu, data = penguenler_temiz %>% 
##     filter(tur == "Chinstrap"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.65742 -0.46033 -0.01862  0.61473  1.69801 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    7.56914    1.55053   4.882 6.99e-06 ***
## gaga_uzunlugu  0.22221    0.03168   7.015 1.53e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8659 on 66 degrees of freedom
## Multiple R-squared:  0.4271, Adjusted R-squared:  0.4184 
## F-statistic: 49.21 on 1 and 66 DF,  p-value: 1.526e-09
reg_gentoo <- lm(gaga_derinligi ~ gaga_uzunlugu, data = penguenler_temiz %>% 
                   filter(tur == "Gentoo"))
summary(reg_gentoo)
## 
## Call:
## lm(formula = gaga_derinligi ~ gaga_uzunlugu, data = penguenler_temiz %>% 
##     filter(tur == "Gentoo"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.57143 -0.52974 -0.04479  0.45417  2.96109 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     5.1210     1.0583   4.839 4.02e-06 ***
## gaga_uzunlugu   0.2076     0.0222   9.352 7.34e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7491 on 117 degrees of freedom
## Multiple R-squared:  0.4277, Adjusted R-squared:  0.4229 
## F-statistic: 87.45 on 1 and 117 DF,  p-value: 7.337e-16

3.4. Sonuçların karşılaştırılması

Grafiklere bakıldığında üçünde de eğim yukarı yönlü yani pozitiftir. Regresyon sonuçlarına baktığımızda:

Adelie: eğim = +0.177, p< .05, pozitif ve anlamlı,

Chinstrap: eğim = +0.222, p< .05, pozitif ve anlamlı,

Gentoo: eğim = +0.208, p< .05, pozitif ve anlamlı.

Türlere teker teker baktığımızda gaga uzunluğu ve gaga derinliği arasındaki ilişki her birinde pozitif ve anlamlıdır. Tür düzeyinde elde edilen ilişki toplam düzeyde elde edilen ilişkiyle çelişmektedir.

4. Görev 3: Yorum ve Tartışma

Toplam veriye baktığımızda gaga uzunluğu ile gaga derinliği arasındaki regresyon çizgisi aşağı doğru gitmekte ve basit doğrusal modelle de bu doğrulanmaktadır. Eğim -.082 ve p<.05’tir. Yani toplam düzeyde gaga uzunluğu ve gaga derinliği arasında negatif ve anlamlı bir ilişki vardır. Veriyi türlere ayırdığımızda ise üç türde de gaga uzunluğu ve gaga derinliği arasında ilişki pozitif ve anlamlı bir ilişki görülmektedir. Toplam analizde görülen negatif eğilimin nedeni, türler arasındaki farklardır. Bazı türlerde gaga genel olarak daha uzun ama derin değil, bazılarında ise daha kısa ama daha derindir. Bu üç grubu tek bir modele koyduğumuzda noktalar türlere göre kümelenmektedir ve bu kümeler arası farkı tek bir doğruyla özetlemeye çalıştığımız için çizgi negatif yönde görünmektedir. Türleri ayrı incelediğimizde ise, aynı tür içindeki bireyler birbirine daha benzer olduğundan gerçek eğilimin gaga uzadıkça derinlik de artar şeklinde olduğu görülmektedir.

Bu durum Simpson Paradoksuna uygun bir örnek çünkü alt gruplarda pozitif yönlü olan ilişki, gruplar birleştirilince negatif oluyor. Biz bu analizi yaparken sadece toplam düzeye baksaydık gaga uzadıkça derinlik azalıyor diyerek yanlış yorum yapmış olacaktık. Analize tür bilgisini dahil edince ilişkinin aslında pozitif olduğunu görmüş olduk. Alt gruplara göre analiz yapmak daha doğru ve güvenilir sonuçlar elde etmemizi sağlar.