📔1. Öğrenme Günlüğüm

Regresyon ve korelasyonla ilgili netleştirdiğim ilk şeyi possum veri seti üzerinden örneklendirmem gerekirse; regresyon yönlüdür korelasyon simetriktir. Yani regresyonda “taill ile totlngth’i açıklayabilir miyim?” diye soruyoruz; korelasyonda ise “taill ile totlngth birlikte ne kadar değişiyor?” diye bakıyoruz. Regresyon yönlü olduğu için x ile y’nin yeri değişince anlam değişir ama korelasyonda x ile y’nin yer değiştirmesi anlamı değiştirmez. Bu nedenle regresyon grafiğinde x ve y’ye atadığımız değişkenlere de ekstra dikkat etmemiz gerekiyor. Regresyon doğrusunu saçılım grafiğinin üzerine eklemek için geom_smooth(method = "lm", se =....) yazıyorduk. TRUE yazarsak regresyon çizginin etrafında güven aralığı çıkar FALSE yazarsak sadece çigiyi görürüz. Bunu önceki derslerde de öğrenmiştim. Bunun üzerine eklediğim şeyse, saçılım grafiğinin üzerine doğrudan denklemi ekleyebiliyor oluşumdu. Denklemi grafiğe eklemek için ggpubr paketindeki stat_regline_equation() fonksiyonunu kullanmamız gerektiğini öğrendim. Çizgi, güven aralığı ve Pearson korelasyon katsayısını (r) görmek için de ggscatter(..., add="reg.line", conf.int=TRUE, cor.coef=TRUE) kullanmak gerekiyor.

Sürekli bir değişkeni kesitlere göre karşılaştırmam gerekiyorsa önce onu aralıklara bölmem gerektiğini, bunu da grup <- cut(xxx, breaks = k) şeklinde yapmam gerektiğini öğrendim. Ardından her kesitte bağımlı değişkenin dağılımını geom_boxplot() ile göstermem, ham noktalar kaybolmasın diye kutuların üstüne geom_jitter(width = 0.x) eklemem ve aykırıların iki kez çizilmemesi için geom_boxplot(outlier.alpha = 0) ayarını yapmam gerektiğini öğrendim. Böylece saçılım–kutu grafiği bağlantısı kurulabiliyor. Kutu genel özeti, jitter ise tek tek gözlemleri görünür kılıyor.

İki değişken arasında kıvrımlı bir ilişki varsa düz lm çizgisi ilişkiyi iyi temsil etmez onun yerine geom_smooth(method = "loess", se = FALSE) ile eğriyi daha doğru gösterebiliriz. Görselde belirgin bir desen yoksa çizgi eklemek yerine “ilişki yok” demek daha doğru olur, bu durumda cor(x, y) genelde 0’a yakındır. Desen parabolikse doğrusal çizgiyle ilikişiyi gösteremeyiz. Bu durumda dönüşüm ya da eğrisel bir yaklaşım (ör. geom_smooth(method="loess")) kullanmalıyız. Grafikte belirgin bir aykırı nokta varsa annotate() ile işaretlemek gerekebilir çünkü çizgiyi ve yorumu etkileyebilir. İki değişken arasındaki ilişkinin her zaman doğrusal olmadığını saçılım grafiğinde tuhaf/kıvrımlı desenler ya da bazen hiç ilişki görülmeyebileceğini öğrendim. Böyle durumlarda dönüşüm yapmak (özellikle log) ilişkiyi daha doğrusal hale getirmeye ve x arttıkça büyüyen yayılımı toparlamaya yardımcı oluyor. Bunu R’da iki yolla uygulayabiliyoruz coord_trans(x = "log10", y = "log10") grafiğin koordinatlarını dönüştürüyor; scale_x_log10() + scale_y_log10() ise eksenleri log ölçeğe alıp grafiği düzgün ayarladığı için raporlamada daha temiz duruyor. Bu yelpaze şeklindeki görselleşmenin de sabit olmayan v aryans olduğunu öğrendim.

R’da Pearson korelasyonunu cor(x, y) ile hesaplayabiliyoruz. Fonksiyon NA görünce varsayılan olarak NA döndürüyor, bu yüzden use argümanını "pairwise.complete.obs" olarak ayarlamak gerekiyor. r pozitifse eğilim yukarıya doğru, negatifse eğilim aşağı doğru, sıfıra yakın ise de görüntü dağınık görünüyor, buradan da ilişki yok anlamını çıkarabiliyoruz.

Görselleştirmede ise elimde sayısal sütunlar olduğunda önce as.matrix() ile veri setini saysal matrise çevirip Hmisc paketindeki rcorr() fonksiyonunu kullanmam gerektiğini öğrendim. Bu fonksiyonla r ve p-değeri matrisleri elde ediliyor. Çıkan r ve p’yi doğrudan görselleştirmek için corrplot(r, ...) yazınca p > 0.05 olan hücreler boş bırakılıyor ve değişkenler benzerliğe göre gruplanıyor. methodkısmına circle veya square yazarak görseli değiştirebiliyoruz. Aynı işi ggplot ekosisteminde yapmak istersek de ggcorrplot(r, hc.order = TRUE, type = "lower", lab = TRUE) kullanabiliriz. Sayı etiketleri kutuların içine yazılır ve tema ayarları (başlık, yazı tipi, ölçek) daha kolay kontrol edilir. Tüm bu öğrendiklerime ek olarak psych paketindeki corPlot(..., stars = TRUE)fonksiyonu hoşuma gitti. stars = TRUE ile anlamlılık düzeylerini hücrelere koymak okumayı hızlandırıyor ve yorumlamayı kolaylaştırıyor. Beğendiğim bir diğer şey ise PerformanceAnalytics içindeki chart.Correlation(...) fonksiyonu, tam all in one ekranı gibi tek bakışta kısa rapor sunuyor.

R’da analiz yapmak ve görselleştirmek çok keyifli. Özellikle çıkan görseller SPSS görsellerinden daha estetik ve canlı duruyor. SPSS ile karşılaştırdığımda verilerle tek tek uğraşmak yerine küçük kodlarla istediğimi elde elebilmek daha kolay ve eğlenceli geliyor bana. SPSS gözümü korkutan bir uygulama ve orada ezbere bir şeyler yapıp sonuçları yorumluyorum ama R kullanırken yazdığım her kodun nedenini bilerek yazıyorum. Bu da ilerleyen zamanlarda analiz yaptığımda analizimden daha emin olmamı sağlayacak ve SPSS kullanırkenki karanlıkta ve gözlerim bağlı yürüyormuşum hissinden kurtulacağım bence :’)

📑 2. Ödev

2.1 Veri Analizi

2.1.1. Veri Setini Seçme

Bu ödevde, R’ın içinde hazır gelen airquality veri setini seçtim. Veri, 1973 yılı Mayıs–Eylül döneminde New York’ta günlük ölçülen hava kalitesi ve hava durumu değişkenlerinden oluşuyor. Toplam 153 gözlem ve 6 değişken var: Ozone (ozon düzeyi, ppb), Solar.R (güneş radyasyonu, Langley), Wind (rüzgâr hızı, mph), Temp (günlük maksimum sıcaklık, °F), ayrıca tarih bilgisi için Month (5–9) ve Day var. Veri günlük olduğu için mevsimsel-sıcaklık etkileri, rüzgârın dağıtıcı etkisi gibi meteorolojik durumları görselleştirebileceğimi düşündüm. Bu seti tercih etmemin temel nedeni, ödevin istediği gibi birden fazla sayısal değişken içermesi ve aralarındaki ilişkiyi hem saçılım grafikleri hem de korelasyon/korelogram ile incelemeye uygun olması.

2.1.2. Kullanılacak Paketleri Yükleme

Gerekli analiz ve görselleştirmeler için dplyr, ggplot2, tidyverse, ggpubr, patchwork, corrplot, GGally, psych, PerformanceAnalytics ve Hmiscpaketlerini yükledim ve kütüphanelerini çağırdım. Kütüphaneleri çağırma sırasındaggpubrpaketinde bir hata aldım.carData paketi olmadığı için ggpubr kütüphanesini çağıramıyordum. Bu yüzden carData paketini de yükledim ve öyle çalıştırabildim.

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(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.1     ✔ stringr   1.5.2
## ✔ lubridate 1.9.4     ✔ tibble    3.3.0
## ✔ purrr     1.1.0     ✔ tidyr     1.3.1
## ✔ readr     2.1.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggpubr)
library(patchwork)
library(corrplot)
## corrplot 0.95 loaded
library(GGally)
library(psych)
## 
## Attaching package: 'psych'
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(PerformanceAnalytics)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## 
## Attaching package: 'PerformanceAnalytics'
## 
## The following object is masked from 'package:graphics':
## 
##     legend
library(Hmisc) 
## 
## Attaching package: 'Hmisc'
## 
## The following object is masked from 'package:psych':
## 
##     describe
## 
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, units

2.1.3. Veri Setini Yükleme ve İnceleme

Kütüphaneleri çağırdıktan sonra R içerisindeki àirquality` veri setini kendim kullanacağım şekilde hava_kalitesi olarak globale ekledim.

data(airquality)
hava_kalitesi0 <- airquality

Veri setini inceleme için dim() ve glimpse()fonksiyonlarını kullandım. Veri setinde 153 gözlem ve 6 değişken olduğu ve bunlardan 5 tanesinin integer yani tam sayı, 1 tanesinin ise double yani ondalık sayılardan oluştuğu görülmüştür.

dim(hava_kalitesi0) 
## [1] 153   6
glimpse(hava_kalitesi0)
## Rows: 153
## Columns: 6
## $ Ozone   <int> 41, 36, 12, 18, NA, 28, 23, 19, 8, NA, 7, 16, 11, 14, 18, 14, …
## $ Solar.R <int> 190, 118, 149, 313, NA, NA, 299, 99, 19, 194, NA, 256, 290, 27…
## $ Wind    <dbl> 7.4, 8.0, 12.6, 11.5, 14.3, 14.9, 8.6, 13.8, 20.1, 8.6, 6.9, 9…
## $ Temp    <int> 67, 72, 74, 62, 56, 66, 65, 59, 61, 69, 74, 69, 66, 68, 58, 64…
## $ Month   <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ Day     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,…

2.1.4. Veri Setini Türkçeleştirme

rename() fonksiyonu kullanarak veri setindeki sütunların isimlerini türkçeleştirdikten sonra colnames() ile kontrol ettim.

hava_kalitesi <- hava_kalitesi0 %>%
  rename(
    ozon = Ozone,
    gunes_radyasyonu = Solar.R,
    ruzgar_hizi = Wind,
    sicaklik = Temp,
    ay = Month,
    gun = Day)

colnames(hava_kalitesi)
## [1] "ozon"             "gunes_radyasyonu" "ruzgar_hizi"      "sicaklik"        
## [5] "ay"               "gun"

2.1.5. Eksik Veri Kontrolü ve Temizleme

Eksik verileri kontrol ettiğimizde ozon sütununda 37, gunes_radyasyonu sütununda ise 7 eksik gözlem olduğu tespit edilmiştir. Bu gözlemleri drop_na() fonksiyonu kullanarak temizleyip yeni veri setini oluşturdum. Geriye 111 gözlem kaldı.

colSums(is.na(hava_kalitesi))
##             ozon gunes_radyasyonu      ruzgar_hizi         sicaklik 
##               37                7                0                0 
##               ay              gun 
##                0                0
hava_kalitesi_temiz <- drop_na(hava_kalitesi, ozon, gunes_radyasyonu)
glimpse(hava_kalitesi_temiz)
## Rows: 111
## Columns: 6
## $ ozon             <int> 41, 36, 12, 18, 23, 19, 8, 16, 11, 14, 18, 14, 34, 6,…
## $ gunes_radyasyonu <int> 190, 118, 149, 313, 299, 99, 19, 256, 290, 274, 65, 3…
## $ ruzgar_hizi      <dbl> 7.4, 8.0, 12.6, 11.5, 8.6, 13.8, 20.1, 9.7, 9.2, 10.9…
## $ sicaklik         <int> 67, 72, 74, 62, 65, 59, 61, 69, 66, 68, 58, 64, 66, 5…
## $ ay               <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ gun              <int> 1, 2, 3, 4, 7, 8, 9, 12, 13, 14, 15, 16, 17, 18, 19, …

2.2 Değişkenler Arası İlişkilerin İncelenmesi

Korelasyon ve regresyon yapabilmek için sayısal değerlerin kullanılması gerektiğinden analizlerde ozon, gunes_radyasyonu, ruzgar_hizi ve sicaklik değişkenleriyle çalıştım.

2.2.1. Saçılım Grafiği

Sıcaklık ve ozon ilişkisini incelemk için saçılım grafiği yaptım.

ggplot(hava_kalitesi_temiz, aes(x = sicaklik, y = ozon)) +
  geom_point(color = "steelblue", alpha = 0.8, shape = 16) +
  labs(title = "\nSıcaklık ile Ozon Arasındaki İlişki\n",
       x = "\nSıcaklık\n",
       y = "\nOzon\n") +
  theme_minimal()

Grafik bize sıcaklık arttıkça ozon derişiminin de arttığını gösteriyor. Yani ilişkinin yönü pozitif görünüyor. Noktalar tamamen çizgi gibi dizilmiyor ama yukarı doğru belirgin bir eğilim var.

2.2.2. Regresyon Grafiği

İlişkiyi daha iyi görebilmek için regresyon eğrisi ekledim.

ggplot(hava_kalitesi_temiz, aes(x = sicaklik, y = ozon)) +
  geom_point(color = "steelblue", alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", color = "black", se = TRUE) +
  labs(title = "\nSıcaklık–Ozon İlişkisi (LM çizgisiyle)\n",
       x = "\nSıcaklık\n",
       y = "\nOzon\n") +
  theme_light()
## `geom_smooth()` using formula = 'y ~ x'

Regresyon çizgisi de aradaki ilişkinin pozitif olduğunu gösteriyor.

2.2.3. Kutu Grafiğiyle Bağlantı

Sıcaklığı beş aralığa bölüp her dilimde kaç gün var diye baktım.

hava_kalitesi_kutulu <- hava_kalitesi_temiz %>%
  mutate(sicaklik_cut = cut(sicaklik, breaks = 5, include.lowest = TRUE))

table(hava_kalitesi_kutulu$sicaklik_cut)
## 
## [57,65] (65,73] (73,81] (81,89] (89,97] 
##      14      22      34      28      13

Saçılım grafiği çizdim:

ggplot(hava_kalitesi_kutulu, aes(x = sicaklik_cut, y = ozon)) +
  geom_point() +
  labs(title = "Ozon - Sıcaklık Dilimleri Saçılım Grafiği ",
       x = "Sıcaklık Dilimleri", y = "Ozon") +
  theme_classic()

x artık faktör olduğu için noktalar dikey şeritler halinde görünmektedir. Daha sıcak dilimlerde ozon değerleri daha yüksek yoğunlaşma eğiliminde diyebiliriz.

Kutu grafiğine çevirdim.

ggplot(hava_kalitesi_kutulu, aes(x = sicaklik_cut, y = ozon)) +
  geom_boxplot() +
  labs(title = "Ozon - Sıcaklık Dilimleri Kutu Grafiği",
       x = "Sıcaklık Dilimleri", y = "Ozon") +
  theme_classic()

Kutu grafiklerine bakınca, sıcaklık arttıkça ozonun ortancası da yukarı kayıyor. Üst tarafın açılması da sıcak günlerde daha yüksek ozon görülebildiğini gösteriyor.

Kutu + Jitter yaptım.

ggplot(hava_kalitesi_kutulu, aes(x = sicaklik_cut, y = ozon)) +
  geom_boxplot(outlier.alpha = 0) + 
  geom_jitter(color = "skyblue", width = 0.2) +
  labs(title = "Ozon - Sıcaklık Dilimleri Kutu + Jitter",
       x = "Sıcaklık Dilimleri", y = "Ozon") +
  theme_classic()

sicaklik sütununu dilimleyip ozon’u y eksenine koydum. Tek başına nokta grafiği dağınıktı. Kutu grafiği oluşturup medyanları gördüm üstüne jitter ekleyince ham veri de görünür oldu. Sıcaklık arttıkça kutuların medyanı yükseliyor, noktalar da yukarı toplanıyor. Yani sıcaklık arttıkça ozon da artıyor, ilişki pozitif yönlü.

2.2.4. Korelasyon

2.2.4.1 Hmisc ve Corrplot ile grafik oluşturma

library(Hmisc)
library(corrplot)

korelasyon_hava_kalitesi <-select( hava_kalitesi_temiz, ozon, gunes_radyasyonu, ruzgar_hizi, sicaklik)

corr_res <- rcorr(as.matrix(korelasyon_hava_kalitesi))
r <- corr_res$r

corrplot(r, method = "square", type = "lower", order = "hclust",
         addgrid.col = "darkgray", tl.srt = 45)

Korelasyon değerlerini ekleme:

library(ggcorrplot)

ggcorrplot(r,
           hc.order = TRUE,
           type = "lower",
           lab = TRUE)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggcorrplot package.
##   Please report the issue at <https://github.com/kassambara/ggcorrplot/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Anlamlılık değerlerini ekleme:

library(psych)

corPlot(korelasyon_hava_kalitesi, stars = TRUE)

corrplot ısı haritasını kırmızı negatif, mavi pozitif olarak okudum ve kare büyüdükçe ilişkinin gücünün arttığını gördüm. Bu grafikte sıcaklık ile ozon en güçlü pozitif, güneş radyasyonu ile ozon orta düzeyde pozitif, rüzgar hızı ile ozon ise belirgin biçimde negatif görünüyor. Aynı tabloyu ggcorrplot ile hücrelere korelasyon katsayılarını yazarak da kontrol ettim, okumak daha rahattı ve sonuçlar değişmedi; sıcaklık ile ozon yüksek pozitif, rüzgar ile ozon orta ile güçlü arası negatif, güneş ile ozon orta pozitif kaldı. psych paketindeki corPlot grafiğinde ilişkilerin çoğu istatistiksel olarak anlamlı çıktı.

PerformanceAnalytics ve GGally ile grafik oluşturma

library(PerformanceAnalytics)

hava_kalitesi_temiz1 <- hava_kalitesi_temiz %>%
  select(ozon, gunes_radyasyonu, ruzgar_hizi, sicaklik)

chart.Correlation(hava_kalitesi_temiz1, histogram = TRUE, pch = 19)

library(GGally)

ggpairs(
  hava_kalitesi_temiz1,
  title = "Hava Kalitesi: Sayısal Değişkenler Arasındaki İlişkiler",
  upper = list(continuous = wrap("cor", size = 4)),          
  lower = list(continuous = wrap("smooth", alpha = 0.3)),    
  diag  = list(continuous = wrap("densityDiag", alpha = 0.5)))

İki grafiğe birlikte bakıldığında ozonla sıcaklık arasındaki ilişki en belirgin olanı; noktalar sıcaklık arttıkça yukarı doğru toparlanıyor, eğim net pozitif ve güç olarak güçlü düzeyde. Ozonla rüzgâr hızı ters gidiyor; rüzgâr yükseldikçe ozon düşüyor, bu da orta-güçlü bir negatif ilişki gibi duruyor. Ozonla güneş radyasyonu arasında ise pozitif ama daha zayıf bir bağ var; saçılım daha dağınık. Diyagonaldeki dağılımlarda ozon sağa çarpık, birkaç çok yüksek gün var; bu uç değerler ilişkilerin görünüşünü yer yer etkiliyor. Sıcaklığın dağılımı tek tepeye yakın, rüzgâr daha yayvan, güneş radyasyonunda da geniş bir yayılım var. Rüzgâr ile sıcaklık arasında orta düzeyde negatif bir görüntü, güneş radyonu ile rüzgâr arasında ise belirgin bir desen yok gibi. Kısaca biçim çoğunlukla doğrusal, yön sıcaklık ve güneşle pozitif, rüzgârla negatif; güç olarak en net çift ozon-sıcaklık, onu ozon-rüzgâr izliyor.