2. HAFTANIN TEKRARI

2.Hafta derste üç ana başlık altında çalıştık:

  1. veri setinin R ortamına aktarılması ve betimsel incelenmesi,

  2. eksik veri analizi ile eksik verilerle başa çıkma yöntemleri,

  3. uç değer tespiti, normallik varsayımı kontrolü ve çoklu bağlantı analizi

Veri İnceleme

İlk aşamada haven paketi ile SPSS formatındaki SCREEN.sav dosyası R’a aktardık. Veri setini tanımak için head(), summary() ve psych::describe() fonksiyonlarıyla temel yapıyı incelendik. Betimsel istatistikleri daha düzenli sunmak adına gtsummary, vtable ve skimr gibi farklı paketler de denedik. Sonra da DataExplorer::create_report() ile veri setinin otomatik profilleme raporu oluşturduk. Rapor ile verinin genel yapısını, dağılımlarını ve olası sorunlu noktaları hızlıca görmemizi sağladı.

library(haven)
library(dplyr)
library(psych)

screen <- read_sav("SCREEN.sav")
head(screen)
summary(screen)
describe(screen[,-1])

Eksik Veri Analizi

İkinci aşamada naniar paketi ile eksik veri durumu detaylıca inceledik. any_na(), n_miss(), prop_miss(), miss_var_summary() ve miss_var_table() fonksiyonları ile eksik verilerin hangi değişkenlerde yoğunlaştığı belirledik. SCREEN verisinde özellikle INCOME değişkeninde 26 (%5.59) ve ATTHOUSE değişkeninde 1 (%0.22) eksik gözlem bulunduğunu tespit ettik. Çünkü yüzdeyi öğrendikten sonra ona göre kayıp veri ile baş edeceğiz (Örn. MCAR ve %5 altında ise sileceğiz gibi.).

Eksik verilerin mekanizmasını anlamak için mcar_test() ile Little’ın MCAR testi uyguladık Sonuç anlamlı çıktı (p = .044), yani eksiklikler tamamen rastlantısal değil (MCAR değilmiş MAR ve MNAR tespit edemiyorduk zaten tam olarak).

Eksik veriyle başa çıkma yöntemlerinden üçünü uyguladık ve denedik:

  • Listwise deletion (na.omit()): Eksik verili satırların tamamen silinmesi; örneklem kaybına yol açar.
  • Ortalama ile doldurma: Eksik değerlerin değişken ortalamasıyla değiştirilmesi; standart sapmayı düşürür.
  • Çoklu atama (Multiple Imputation): mice paketi ile 5 tamamlanmış veri seti (m=5) oluşturularak PMM (Predictive Mean Matching) yöntemi uygulandı. with() ve pool() fonksiyonları ile sonuçlar birleştirildi. (Şu an en kabul göreni buydu)
library(naniar)
miss_var_summary(screen)
miss_var_table(screen)
mcar_test(data = screen[, c(2,3,4,5,7,8)])

screen$INCOME[is.na(screen$INCOME)] <- mean(screen$INCOME, na.rm = TRUE)

library(mice)
imputed_data <- mice(screen, m = 5, maxit = 50, method = 'pmm', seed = 500)
fit <- with(imputed_data, lm(TIMEDRS ~ ATTHOUSE + INCOME))

pooled_results <- pool(fit)
summary(pooled_results)

Uç Değerler ve Normallik

Üçüncü aşamada veri temizliği ve varsayım kontrolleri yapmışız derste tam takip edememiştim. Tek değişkenli uç değerler için outliers::scores() fonksiyonu ile z-puanları hesaplanıyor. ±3 sınırını aşan değerler inceledik.. ggplot2 ile histogram ve boxplot grafikleri çizdik (Hatırladığım en pratik kod kümesi); plotly ile interaktif boxplot’lar üzerinden uç değerlerin ID numaraları belirlemişiz. Çok değişkenli uç değerler için mahalanobis() uzaklığı hesapladık α = .001 düzeyinde Ki-kare kritik değeri belirlenerek bu değeri aşan 9 gözlemi veriden çıkarıldı. Çoklu bağlantı için VIF değerleri hesaplandı. Tüm değişkenlerde VIF < 5 olduğu görüldü, dolayısıyla çoklu bağlantı sorunu yoktu.

library(outliers)
z.scores <- screen %>% select(2:5) %>% scores(type = "z") %>% round(2)

veri <- screen[, 1:5]
md <- mahalanobis(veri, center = colMeans(veri), cov = cov(veri))
cutoff <- qchisq(p = 1 - 0.001, df = ncol(veri))
ucdegerler <- which(md > cutoff)

model <- lm(SUBNO ~ TIMEDRS + ATTDRUG + ATTHOUSE + INCOME + RACE + MSTATUS, data = screen)

library(olsrr)

ols_vif_tol(model)

FARKLI VERİ SETİ İLE UYGULAMA: airquality

Derste öğrendiğimiz tüm aşamaları, R’ın kendi airquality veri seti üzerinde uygulayacağım (Bu veri setini kullanmayı chat önerdi çünkü eksik veri de içeren ama kontrol edebileceğim veri seti gelmedi aklıma. İlk dönemki derste proje ödevindeki kullandığım veri seti çok teferruatlıydı onu kullansaydım vei setindeki değişkenleri kavrama çalışana kadar kodları deneyemeyecektim.)

CHAT: “Bu veri seti New York’ta 1973 yılında Mayıs-Eylül ayları arasında günlük hava kalitesi ölçümlerini içermektedir. Ozon değişkeninin sağa çarpık olması ve doğal eksik veriler içermesi, derste işlediğiniz konuları denemek için uygun veri setlerinden birisidir.”

Paketlerin Yüklenmesi

library(dplyr)
library(psych)
library(naniar)
library(ggplot2)
library(plotly)
library(moments)
library(mice)
library(outliers)
library(corrplot)
library(olsrr)
library(knitr)

Verinin Yüklenmesi ve İncelenmesi

veri <- airquality
head(veri)
##   Ozone Solar.R Wind Temp Month Day
## 1    41     190  7.4   67     5   1
## 2    36     118  8.0   72     5   2
## 3    12     149 12.6   74     5   3
## 4    18     313 11.5   62     5   4
## 5    NA      NA 14.3   56     5   5
## 6    28      NA 14.9   66     5   6
str(veri)
## 'data.frame':    153 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...
describe(veri)
##         vars   n   mean    sd median trimmed   mad  min   max range  skew
## Ozone      1 116  42.13 32.99   31.5   37.80 25.95  1.0 168.0   167  1.21
## Solar.R    2 146 185.93 90.06  205.0  190.34 98.59  7.0 334.0   327 -0.42
## Wind       3 153   9.96  3.52    9.7    9.87  3.41  1.7  20.7    19  0.34
## Temp       4 153  77.88  9.47   79.0   78.28  8.90 56.0  97.0    41 -0.37
## Month      5 153   6.99  1.42    7.0    6.99  1.48  5.0   9.0     4  0.00
## Day        6 153  15.80  8.86   16.0   15.80 11.86  1.0  31.0    30  0.00
##         kurtosis   se
## Ozone       1.11 3.06
## Solar.R    -1.00 7.45
## Wind        0.03 0.28
## Temp       -0.46 0.77
## Month      -1.32 0.11
## Day        -1.22 0.72

CHAT: “Veri seti 153 gözlem ve 6 değişkenden oluşuyor. Değişkenler: Ozone (ozon konsantrasyonu, ppb), Solar.R (güneş radyasyonu, lang), Wind (rüzgar hızı, mph), Temp (sıcaklık, °F), Month ve Day.”

Eksik Veri Analizi

Eksik Veri Durumunun İncelenmesi

any_na(veri)
## [1] TRUE
n_miss(veri)
## [1] 44
prop_miss(veri)
## [1] 0.0479
miss_var_summary(veri)
## # A tibble: 6 × 3
##   variable n_miss pct_miss
##   <chr>     <int>    <num>
## 1 Ozone        37    24.2 
## 2 Solar.R       7     4.58
## 3 Wind          0     0   
## 4 Temp          0     0   
## 5 Month         0     0   
## 6 Day           0     0
miss_var_table(veri)
## # A tibble: 3 × 3
##   n_miss_in_var n_vars pct_vars
##           <int>  <int>    <dbl>
## 1             0      4     66.7
## 2             7      1     16.7
## 3            37      1     16.7
miss_case_table(veri)
## # A tibble: 3 × 3
##   n_miss_in_case n_cases pct_cases
##            <int>   <int>     <dbl>
## 1              0     111     72.5 
## 2              1      40     26.1 
## 3              2       2      1.31

Ozone değişkeninde ve Solar.R değişkeninde eksik gözlem bulunuyor. Wind, temp, month ve day değişkenlerinde ise hiç eksiklik yokmuş.

Eksik Verinin Görselleştirilmesi

vis_miss(veri) + theme(axis.text.x = element_text(angle = 80))

vis_miss(veri) + theme(axis.text.x = element_text(angle = 100))

gg_miss_upset(veri)

UpSet grafiğinde eksik verilerin büyük kısmının sadece Ozone değişkeninde tek başına gerçekleştiğini, bir kısmının ise hem Ozone hem Solar.R’da birlikte eksik olduğunu görüyoruz. Bu önemli çünkü belki problem durumuna göre ikili doldurmalar yapılabilir ama tabii o zaman farklı veri setleri ile alt problemelere bakıldığı için geçerli olmayacaktır, bakacağız.

MCAR Testi

mcar_test(veri)
## # A tibble: 1 × 4
##   statistic    df p.value missing.patterns
##       <dbl> <dbl>   <dbl>            <int>
## 1      35.1    14 0.00142                4

DERS NOT: Little’ın MCAR testinin sonucu p değerine göre yorumlanır. p > .05 ise eksiklikler tamamen rastlantısaldır (MCAR); p < .05 ise veriler MCAR değildir ve eksikliklerin sistematik bir yapıda olabileceği düşünülmelidir.

Eksik Veri ile Başa Çıkma

Yöntem 1: Listwise Deletion

veri_listwise <- na.omit(veri)
cat("Orijinal veri:", nrow(veri), "gözlem\n")
## Orijinal veri: 153 gözlem
cat("Listwise sonrası:", nrow(veri_listwise), "gözlem\n")
## Listwise sonrası: 111 gözlem
cat("Kaybedilen gözlem:", nrow(veri) - nrow(veri_listwise), "\n")
## Kaybedilen gözlem: 42

Yöntem 2: Ortalama ile Doldurma

veri_ort <- veri
cat("Orijinal Ozone SD:", round(sd(veri$Ozone, na.rm = TRUE), 2), "\n")
## Orijinal Ozone SD: 33
veri_ort$Ozone[is.na(veri_ort$Ozone)] <- mean(veri_ort$Ozone, na.rm = TRUE)
veri_ort$Solar.R[is.na(veri_ort$Solar.R)] <- mean(veri_ort$Solar.R, na.rm = TRUE)

cat("Ort doldurma sonrası Ozone SD:", round(sd(veri_ort$Ozone), 2), "\n")
## Ort doldurma sonrası Ozone SD: 28.7
any_na(veri_ort)
## [1] FALSE

Beklendiği gibi, ortalama ile doldurma standart sapmayı düşürdü. Bu yöntem varyansı küçülttüğü için ideal değildi ama medyanda daha az etki etmesi lazımdı onu deneyeceğim hemen altta: (Olmadı aynı çıktı ?)

veri_medyan <- veri

veri_medyan$Ozone[is.na(veri_medyan$Ozone)] <- mean(veri_medyan$Ozone, na.rm = TRUE)
veri_medyan$Solar.R[is.na(veri_medyan$Solar.R)] <- mean(veri_medyan$Solar.R, na.rm = TRUE)

cat("Medyan doldurma sonrası Ozone SD:", round(sd(veri_medyan$Ozone), 2), "\n")
## Medyan doldurma sonrası Ozone SD: 28.7
any_na(veri_medyan)
## [1] FALSE
cat("Ort doldurma sonrası Ozone SD:", round(sd(veri_ort$Ozone), 10), "\n")
## Ort doldurma sonrası Ozone SD: 28.7
cat("Medyan doldurma sonrası Ozone SD:", round(sd(veri_medyan$Ozone), 10), "\n")
## Medyan doldurma sonrası Ozone SD: 28.7

Yöntem 3: Çoklu Atama (Multiple Imputation-mice ile)

imputed <- mice(veri, m = 5, maxit = 50, method = 'pmm', seed = 123, print = FALSE)
md.pattern(veri, rotate.names = TRUE)

##     Wind Temp Month Day Solar.R Ozone   
## 111    1    1     1   1       1     1  0
## 35     1    1     1   1       1     0  1
## 5      1    1     1   1       0     1  1
## 2      1    1     1   1       0     0  2
##        0    0     0   0       7    37 44
fit_imp <- with(imputed, lm(Ozone ~ Solar.R + Wind + Temp))
pooled <- pool(fit_imp)
summary(pooled)
##          term estimate std.error statistic   df  p.value
## 1 (Intercept) -63.7747   21.5203     -2.96 50.3 4.64e-03
## 2     Solar.R   0.0619    0.0218      2.84 63.0 6.01e-03
## 3        Wind  -2.9529    0.6373     -4.63 38.7 4.00e-05
## 4        Temp   1.5716    0.2411      6.52 44.9 5.37e-08

Çoklu atama ile oluşturulan 5 veri seti üzerinden regresyon modeli kurdum ve pool() fonksiyonu ile sonuçlar birleştirdim. PMM yöntemi, sürekli değişkenler için en sık önerilen atama yöntemi olarak geçiyor

Temiz Veriyi Hazırlama

Bundan sonraki analizler için eksik verileri ortalama ile doldurup devam edeceğim

veri_temiz <- veri
veri_temiz$Ozone[is.na(veri_temiz$Ozone)] <- mean(veri_temiz$Ozone, na.rm = TRUE)
veri_temiz$Solar.R[is.na(veri_temiz$Solar.R)] <- mean(veri_temiz$Solar.R, na.rm = TRUE)

Uç Değer Analizi

Tek Değişkenli Uç Değerler (Z-Puanları)

sayisal_degiskenler <- veri_temiz %>% select(Ozone, Solar.R, Wind, Temp)

z_scores <- scores(sayisal_degiskenler, type = "z") %>% round(2)
head(z_scores)
##   Ozone Solar.R  Wind  Temp
## 1 -0.04    0.05 -0.73 -1.15
## 2 -0.21   -0.77 -0.56 -0.62
## 3 -1.05   -0.42  0.75 -0.41
## 4 -0.84    1.44  0.44 -1.68
## 5  0.00    0.00  1.23 -2.31
## 6 -0.49    0.00  1.40 -1.26
data.frame(
  Degisken = names(sayisal_degiskenler),
  Min_Z = round(apply(z_scores, 2, min), 2),
  Max_Z = round(apply(z_scores, 2, max), 2)
) %>% kable()
Degisken Min_Z Max_Z
Ozone Ozone -1.43 4.39
Solar.R Solar.R -2.03 1.68
Wind Wind -2.34 3.05
Temp Temp -2.31 2.02

Histogram ve Boxplot ile Görselleştirme

ggplot(veri_temiz, aes(x = Ozone)) +
  geom_histogram(aes(y = ..density..), bins = 25, fill = "pink", alpha = 0.7) +
  geom_density(fill = "pink", alpha = 0.3) +
  geom_vline(xintercept = mean(veri_temiz$Ozone), color = "red", linetype = "dashed") +
  annotate("text", x = mean(veri_temiz$Ozone) + 10, y = 0.015,
           label = paste("Ort =", round(mean(veri_temiz$Ozone), 1)), color = "red") +
  theme_minimal() +
  labs(title = "Ozone Değişkeninin Dağılımı", x = "Ozone (ppb)", y = "Yoğunluk")

Ozone değişkeni derste incelediğimiz TIMEDRS’e benzer şekilde sağa çarpık bir dağılım sergiliyor.

ggplot(veri_temiz, aes(y = Ozone)) +
  geom_boxplot(fill = "pink", alpha = 0.5) +
  theme_minimal() +
  labs(title = "Ozone boxplot", y = "Ozone")

plot_ly(y = veri_temiz$Ozone, type = 'box', name = "Ozone") %>%
  layout(title = 'Ozone interaktiff boxplot')
ggplot(veri_temiz, aes(x = factor(Month), y = Ozone, fill = factor(Month))) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Aylara göre Ozone", x = "Ay", y = "Ozone", fill = "Ay")

Çok Değişkenli Uç Değerler (Mahalanobis uızaklığı)

md <- mahalanobis(sayisal_degiskenler,
                  center = colMeans(sayisal_degiskenler),
                  cov = cov(sayisal_degiskenler))

alpha <- 0.001
kritik_deger <- qchisq(p = 1 - alpha, df = ncol(sayisal_degiskenler))
cat("Ki-kare kritik değeri:", round(kritik_deger, 2), "\n")
## Ki-kare kritik değeri: 18.5
uc_indexler <- which(md > kritik_deger)
cat("Çok değişkenli uç değer sayısı:", length(uc_indexler), "\n")
## Çok değişkenli uç değer sayısı: 1
if (length(uc_indexler) > 0) {
  veri_temiz[uc_indexler, c("Ozone", "Solar.R", "Wind", "Temp")] %>% kable()
}
Ozone Solar.R Wind Temp
117 168 238 3.4 81 
if (length(uc_indexler) > 0) {
  veri_son <- veri_temiz[-uc_indexler, ]
  cat("Uç değerler çıkarıldıktan sonra:", nrow(veri_son), "gözlem\n")
} else {
  veri_son <- veri_temiz
  cat("Çıkarılacak uç değer bulunamadı.\n")
}
## Uç değerler çıkarıldıktan sonra: 152 gözlem

Normallik Varsayımı

Çarpıklık ve Basıklık

carpiklik_df <- data.frame(
  Degisken = c("Ozone", "Solar.R", "Wind", "Temp"),
  Carpiklik = c(skewness(veri_son$Ozone), skewness(veri_son$Solar.R),
                skewness(veri_son$Wind), skewness(veri_son$Temp)),
  Basiklik = c(kurtosis(veri_son$Ozone), kurtosis(veri_son$Solar.R),
               kurtosis(veri_son$Wind), kurtosis(veri_son$Temp))
)
carpiklik_df %>% mutate(across(where(is.numeric), ~round(., 2))) %>% kable()
Degisken Carpiklik Basiklik
Ozone 1.13 4.14
Solar.R -0.42 2.11
Wind 0.36 3.09
Temp -0.37 2.55

QQ Plot

ggplot(veri_son, aes(sample = Ozone)) +
  geom_qq() +
  geom_qq_line(color = "pink") +
  theme_minimal() +
  labs(title = "Ozone QQ Plot")

ggplot(veri_son, aes(sample = Wind)) +
  geom_qq() +
  geom_qq_line(color = "pink") +
  theme_minimal() +
  labs(title = "Wind QQ Plot")

Logaritmik Dönüşüm

Ozone değişkeni sağa çarpık olduğu için log dönüşüm uygulayalım (derste TIMEDRS için yapmıştık):

veri_son$log_Ozone <- log(veri_son$Ozone + 1)

cat("Dönüşüm öncesi:\n")
## Dönüşüm öncesi:
describe(veri_son$Ozone)
##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 152 41.3 26.9   42.1      38 28.4   1 135   134 1.12     1.08 2.18
cat("\nDönüşüm sonrası:\n")
## 
## Dönüşüm sonrası:
describe(veri_son$log_Ozone)
##    vars   n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
## X1    1 152 3.53 0.71   3.76    3.57 0.64 0.69 4.91  4.22 -0.69     0.86 0.06
par(mfrow = c(1, 2))

hist(veri_son$Ozone, col = "pink", main = "Orijinal Ozone", xlab = "Ozone")

hist(veri_son$log_Ozone, col = "pink", main = "Log(Ozone + 1)", xlab = "Log Ozone")

Log dönüşüm sonrası çarpıklık belirgin şekilde azalmış, dağılım normale daha çok yaklaştı.

Korelasyon ve Çoklu Bağlantı

Korelasyon Matrisi

kor_matris <-
  cor(veri_son[, c("Ozone", "Solar.R", "Wind", "Temp")])


round(kor_matris, 2) %>% kable()
Ozone Solar.R Wind Temp
Ozone 1.00 0.31 -0.52 0.64
Solar.R 0.31 1.00 -0.05 0.26
Wind -0.52 -0.05 1.00 -0.46
Temp 0.64 0.26 -0.46 1.00 
corrplot(kor_matris, method = "number", type = "upper",
         tl.col = "darkred", tl.srt = 45,
         title = "Korelasyon matrisi", mar = c(0,0,1,0))

pairs(veri_son[, c("Ozone", "Solar.R", "Wind", "Temp")],
      main = "Değişkenler saçılım diyagramı")

VIF

model <- lm(Ozone ~ Solar.R + Wind + Temp, data = veri_son)
ols_vif_tol(model) %>% kable(digits = 2)
Variables Tolerance VIF
Solar.R 0.93 1.08
Wind 0.78 1.28
Temp 0.73 1.37

VIF değerlerinin 5’in altında olması çoklu bağlantı sorunu olmadığını gösterirdi.

YANSITMA

2.Haftaki derste, çok değişkenli istatistik analizlerinin ön koşulu olan veri hazırlama sürecini öğrendik. Öğrendiğim en önemli şey, analiz yapmadan önce verinin “analiz edilebilir” duruma getirilmesinin ne kadar sistematik bir süreç olduğuydu. Özellikle eksik veri konusunda dikkatimi çeken nokta, listwise deletion’ın ne kadar tehlikeli olabileceğiydi. airquality verisinde bile 42 gözlem kaybediliyordu bu toplam verinin yaklaşık %27’si demek. mice paketi ile yapılan çoklu atama çok daha güvenilir bir yöntem olarak öne çıkıyor; bireyin diğer yanıtlarına en çok benzeyen kişilere bakılarak eksik veriler tahmin ediliyor. Uç değer tespitinde tek değişkenli (z-puanları) ve çok değişkenli yaklaşımları birlikte kullanmanın gerekliliğini kavradım. Bir gözlem tek değişkende normal sınırlarda olsa bile, birden fazla değişken birlikte ele alındığında uç değer olarak ortaya çıkabiliyor. Normallik testleri konusunda, Kolmogorov-Smirnov gibi örneklem büyüklüğünden çok etkilenen testler yerine Jarque-Bera testinin ve çarpıklık/basıklık değerlerinin daha güvenilir olduğunu öğrendim ama Jarqur-Bera testini deneyemedim haftaya onu da deneyeceğim. Yorumlamasını becerememem gibi geldi ve basıklık çarpıklık alışkın olduğum ve kolay olduğu için onlara baktım. Sağa çarpık dağılımlarda log dönüşümün işe yaradığını uygulamalı olarak gördüm. naniar paketinin eksik veri analizindeki kullanım kolaylığı gerçekten etkileyici. Özellikle vis_miss() ve gg_miss_upset() ile eksik verilerin dağılımını görselleştirmek, sayısal tablolardan çok daha anlaşılır.