Bu yazımızda belirlediğimiz regresyon modellerini veri setine dinamik olarak uygulayarak, R ile istatistiksel programlama bakış açısıyla regresyon analizi yapacağız. Oluşturulan modellerin özet sonuçlarını ve bilgi kriterlerini tek bir tabloda göreceğiz. Analizin kolay ve anlaşılabilir olabilmesi için yalnızca lineer regresyon modeli ile çalışılmıştır. Ancak bu kod bloğu ve mantığı, farklı modeller için de rahatlıkla uyarlanabilir.
Gerekli Kütüphaneler ve Tanımlamalar
- tidyverse: Birbiriyle uyumlu çalışan dplyr, tidyr, tibble gibi paketleri tek bir çatı altında toplar.
- magrittr: R’ın semantik yapısını geliştiren bir grup operatör sunar.
- broom: İstatiksel veri analizi nesnelerini düzenli veri seti yapısına dönüştürür.
Sistemimizde eksik olan kütüphanelerin kurulumunu yapalım.
install.packages(c("tidyverse","magrittr", "broom"))
Her bir kütüphaneyi R oturumumuzda aktif hale getirelim.
library(tidyverse)
library(magrittr)
library(broom)
Analize Hazırlık
R’ın örnek veri setlerinden olan mtcars veri setinden “mpg”, “disp”, “hp”, “drat”, “wt”, “qsec” sütunlarını seçerek analiz veri setini oluşturalım.
df <- mtcars %>%
select("mpg","disp","hp","drat","wt","qsec")
df %>% head()
## mpg disp hp drat wt qsec
## Mazda RX4 21.0 160 110 3.90 2.620 16.46
## Mazda RX4 Wag 21.0 160 110 3.90 2.875 17.02
## Datsun 710 22.8 108 93 3.85 2.320 18.61
## Hornet 4 Drive 21.4 258 110 3.08 3.215 19.44
## Hornet Sportabout 18.7 360 175 3.15 3.440 17.02
## Valiant 18.1 225 105 2.76 3.460 20.22
Veri setinin mpg değişkenini bağımlı değişken olarak belirleyelim. Geri kalan bağımsız değişkenlerle modeller oluşturarak bir liste yaratalım.
y <- "mpg"
x <- list( model1=c("disp"), model2=c("hp"), model3=c("drat"),
model4=c("wt"), mdoel5=c("qsec"),model6=c("disp","hp"),
mdoel7=c("disp","drat"), model8=c("disp","hp","drat"),
model9=c("disp","wt","qsec") )
x listesinde belirlenen modelleri sprintf fonksiyonunu kullanarak regresyon denklemi oluşturabiliriz. lapply fonksiyonu ile liste elemanı olan her bir modele sprintf fonksiyonunu uygulayalım ve modelFamily listesini oluşturalım. Listenin her bir elemanı regresyon denklemini içermektedir.
modelFamily <- lapply(x, function(z) sprintf("%s ~ %s", y,
ifelse(length(z)>1,
paste(z, collapse=" + "),
z)))
Analize geçmeden önce sonuç şablonlarımızı oluşturalım. Regresyon sonuçlarını özetleyeceğimiz tablo yapısını tribble fonksiyonu ile oluşturabiliriz. Ayrıca lm fonksiyonundan çıkan fit sonuçlarını tutacağımız boş bir liste yaratalım.
resTable <- tribble(~id, ~model, ~formula, ~r2, ~adjr2, ~Fstat, ~Fstat_p, ~AIC, ~BIC,
~logLik, ~tTest)
resModel <- list()
Analiz
Bu bölümde modelFamily listesinde belirtilen her bir modelin regresyon sonuçlarını for yapısı ile dinamik olarak elde edeceğiz.
try fonksiyonuyla modelleri fitlemeyi deneyerek error kontrolünü sağlayabiliriz. Böylelikle fitlemede yaşanan bir hata durumunda bir sonraki model için analize devam edebileceğiz.
bloom paketi içindeki glace fonksiyonunu kullanarak modellerin bir satırlık özetlerini veri seti yapısında alalım. Bu veri setinin bazı sütun değerleri ile analiz tablomuzu dolduracağız.
Aynı zamanda oluşturduğumuz resTable veri setinin tTest sütununda, Intercepts dışındaki tüm değişkenler için t-testi uygulayalım. Her bir değişkenin %95 güven seviyesinde istatistiksel olarak anlamlı olması durumunda TRUE, aksi durumda FALSE sonucu elde edilecektir.
for(i in 1:length(modelFamily)){
temp <- try(lm(as.formula(modelFamily[[i]]), data = df))
if(inherits(temp, "try-error")){
resTable %<>% add_row(id = i,
model = names(modelFamily[[i]]),
formula = modelFamily[[i]],
r2 = NA_real_,
adjr2 = NA_real_,
Fstat = NA_real_,
Fstat_p =NA_real_,
AIC = NA_real_,
BIC = NA_real_,
logLik = NA_real_,
tTest= NA)
} else {
tempGlance <- glance(temp)
resTable %<>% add_row(id = i,
model = names(modelFamily[i]),
formula = modelFamily[[i]],
r2 = tempGlance$r.squared,
adjr2 = tempGlance$adj.r.squared,
Fstat = tempGlance$statistic,
Fstat_p = tempGlance$p.value,
AIC = tempGlance$AIC,
BIC = tempGlance$BIC,
logLik = tempGlance$logLik,
tTest = all(summary(temp)$coefficients[-1, "Pr(>|t|)"]< 0.05))
resModel[[i]] <- temp
rm(tempGlance)
}
rm(temp)
}
Oluşturduğumuz 9 modelin özet sonuçlarını resTable tablosunda görelim.
resTable
## # A tibble: 9 x 11
## id model formula r2 adjr2 Fstat Fstat_p AIC BIC logLik tTest
## <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 1 model1 mpg ~ ~ 0.718 0.709 76.5 9.38e-10 170. 175. -82.1 TRUE
## 2 2 model2 mpg ~ ~ 0.602 0.589 45.5 1.79e- 7 181. 186. -87.6 TRUE
## 3 3 model3 mpg ~ ~ 0.464 0.446 26.0 1.78e- 5 191. 195. -92.4 TRUE
## 4 4 model4 mpg ~ ~ 0.753 0.745 91.4 1.29e-10 166. 170. -80.0 TRUE
## 5 5 mdoel5 mpg ~ ~ 0.175 0.148 6.38 1.71e- 2 205. 209. -99.3 TRUE
## 6 6 model6 mpg ~ ~ 0.748 0.731 43.1 2.06e- 9 169. 174. -80.3 FALSE
## 7 7 mdoel7 mpg ~ ~ 0.731 0.712 39.4 5.39e- 9 171. 177. -81.4 FALSE
## 8 8 model8 mpg ~ ~ 0.775 0.751 32.2 3.28e- 9 167. 174. -78.5 FALSE
## 9 9 model9 mpg ~ ~ 0.826 0.808 44.4 8.95e-11 159. 166. -74.4 FALSE
Oluşan resTable veri setine aşağıdaki filtreleme ve sıralamaları uygulayacağız.
1. Fstat_p değeri 0.05 ’den küçük (%95 güven seviyesinde istatistiksel olarak anlamlı)
2. tTest değeri TRUE (%95 güven seviyesinde tüm değişkenler istatistiksel olarak anlamlı)
3. AIC değerine göre küçükten büyüğe, ardından da adjusted R2 değerine göre büyükten küçüğe sıralı
Bu düzenlemeler sonrası ilk sıradaki model en iyi model olarak seçilecektir.
selectedTable <- resTable %>% filter(Fstat_p < 0.05, tTest) %>% arrange(AIC,desc(adjr2))
selectedTable
## # A tibble: 5 x 11
## id model formula r2 adjr2 Fstat Fstat_p AIC BIC logLik tTest
## <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 4 model4 mpg ~ ~ 0.753 0.745 91.4 1.29e-10 166. 170. -80.0 TRUE
## 2 1 model1 mpg ~ ~ 0.718 0.709 76.5 9.38e-10 170. 175. -82.1 TRUE
## 3 2 model2 mpg ~ ~ 0.602 0.589 45.5 1.79e- 7 181. 186. -87.6 TRUE
## 4 3 model3 mpg ~ ~ 0.464 0.446 26.0 1.78e- 5 191. 195. -92.4 TRUE
## 5 5 mdoel5 mpg ~ ~ 0.175 0.148 6.38 1.71e- 2 205. 209. -99.3 TRUE
Oluşturduğumuz selectedTable veri setinde öne çıkan modelin id numarasını, id sütununun ilk satır değeri ile alabiliriz. Bu model seçim kriterlerimize göre en iyi modeldir. Bu id numarası ile resModel listesinden alacağımız modelin özet bilgisini summary fonksiyonuyla görelim.
summary(resModel[[selectedTable$id[1]]] )
##
## Call:
## lm(formula = as.formula(modelFamily[[i]]), data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10
Not: lmtest paketi içerisinde yer alan testler kullanılarak model ile ilgili diagnostic check’ler yapılabilmektedir. Bu konu,farklı bir yazımızda detaylı olarak ele alınacaktır.
