Ana içeriğe geç

R’de Bootstrap Eğitimi

Bu eğitimde, farklı türde bootstrap güven aralıkları elde etmek için boot paketini nasıl kullanacağınızı öğreneceksiniz.
· 12 dk. oku

Yapay Zeka ile Keşfedin

Bootstrap nedir?

Bootstrap, bir anakütle hakkında çıkarım yapmak için örneklem verilerini kullanan bir yöntemdir. Bradley Efron bunu ilk kez 1979’da bu makalede tanıttı. Bootstrap, örneklem verilerinden yerine koyarak örnekleme yapmaya dayanır. Bu teknik, herhangi bir istatistiğin standart hatasını tahmin etmek ve onun için bir güven aralığı (GA) elde etmek amacıyla kullanılabilir. GA’nın kapalı bir formu olmadığında veya çok karmaşık olduğunda bootstrap özellikle kullanışlıdır.

n elemandan oluşan bir örneklemimiz olduğunu varsayalım: X = {x1, x2, …, xn} ve bazı bir istatistik için GA ile ilgileniyoruz: T = t(X). Bootstrap çerçevesi basittir. Sadece aşağıdaki şemayı R kez tekrarlarız: i. tekrarda, mevcut örneklemden yerine koyarak n eleman örnekleriz (bazıları birden fazla kez seçilecektir). Bu yeni örnekleme i. bootstrap örneklemi, Xi diyelim ve istenen istatistiği hesaplayalım: Ti = t(Xi).

Sonuç olarak, istatistiğimizin R değerini elde ederiz: T1, T2, …, TR. Bunlara T’nin bootstrap gerçekleşmeleri ya da T’nin bootstrap dağılımı deriz. Buna dayanarak T için GA hesaplayabiliriz. Bunu yapmanın birkaç yolu vardır. Yüzdelikleri almak en kolay görünenidir.

Uygulamada Bootstrap

(Bir kez daha) iyi bilinen iris veri setini kullanalım. İlk birkaç satıra göz atın:

 head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

Diyelim ki medyan Sepal.Length, medyan Sepal.Width ve bu ikisi arasındaki Spearman sıra korelasyon katsayısı için GA’ları bulmak istiyoruz. R’nin boot paketini ve adı... boot olan bir fonksiyonu kullanacağız. Bu gücü kullanmak için, yeniden örneklenmiş verilerden istatistiğimizi(lerimizi) hesaplayan bir fonksiyon oluşturmamız gerekir. En az iki argümana sahip olmalıdır: bir dataset ve bir bootstrap örneklemi oluşturmak için veri setinden seçilen öğelerin indices vektörü.

Aynı anda birden fazla istatistik için GA hesaplamak istersek, fonksiyonumuz bunları tek bir vektör olarak döndürmelidir.

Bizim örneğimiz için şu şekilde görünebilir:

 library(boot)

foo <- function(data, indices){
  dt<-data[indices,]
  c(
    cor(dt[,1], dt[,2], method='s'),
    median(dt[,1]),
    median(dt[,2])
  )
}

foo, data içinden istenen öğeleri (indices içinde numaraları bulunan) seçer ve ilk iki sütunun korelasyon katsayısını hesaplar (method='s', Spearman sıra katsayısını seçer; method='p' Pearson katsayısı ile sonuçlanır) ve bunların medyanlarını döndürür.

Ekstra argümanlar da ekleyebiliriz; örneğin, kullanıcının korelasyon katsayısı türünü seçmesine izin verelim:

 foo <- function(data, indices, cor.type){
  dt<-data[indices,]
  c(
    cor(dt[,1], dt[,2], method=cor.type),
    median(dt[,1]),
    median(dt[,2])
  )
}

Bu eğitimi daha genel hale getirmek için, foo’nun ikinci sürümünü kullanacağım.

Artık boot fonksiyonunu kullanabiliriz. Ona veri setinin adını, az önce oluşturduğumuz fonksiyonu, tekrar sayısını (R) ve fonksiyonumuzun ek argümanlarını (örneğin cor.type) sağlamalıyız. Aşağıda, örneğin tekrarlanabilirliği için set.seed kullanıyorum.

 set.seed(12345)
myBootstrap <- boot(iris, foo, R=1000, cor.type='s')

boot fonksiyonu, adı... (evet, doğru bildiniz!) boot olan bir sınıfa ait bir nesne döndürür. İki ilginç öğesi vardır. $t, bootstrap prosedürüyle üretilen istatistiğimizin R değerlerini (T’nin bootstrap gerçekleşmeleri) içerir:

 head(myBootstrap$t)

##             [,1] [,2] [,3]
## [1,] -0.26405188 5.70    3
## [2,] -0.12973299 5.80    3
## [3,] -0.07972066 5.75    3
## [4,] -0.16122705 6.00    3
## [5,] -0.20664808 6.00    3
## [6,] -0.12221170 5.80    3

$t0, istatistiğimizin orijinal, tam veri setindeki değerlerini içerir:

 myBootstrap$t0

## [1] -0.1667777  5.8000000  3.0000000

boot nesnesini konsola yazdırmak daha fazla bilgi verir:

 myBootstrap

##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = iris, statistic = foo, R = 1000, cor.type = "s")
##
##
## Bootstrap Statistics :
##       original       bias    std. error
## t1* -0.1667777  0.002546391  0.07573983
## t2*  5.8000000 -0.013350000  0.10295571
## t3*  3.0000000  0.007900000  0.02726414

original, $t0 ile aynıdır. bias, bootstrap gerçekleşmelerinin ortalaması (yani $t’den gelen), T’nin bootstrap tahmini ile orijinal veri setindeki değer ($t0’daki) arasındaki farktır.

 colMeans(myBootstrap$t)-myBootstrap$t0

## [1]  0.002546391 -0.013350000  0.007900000

std. error, bootstrap tahmininin standart hatasıdır ve bootstrap gerçekleşmelerinin standart sapmasına eşittir.

 apply(myBootstrap$t,2,sd)

## [1] 0.07573983 0.10295571 0.02726414

Farklı Bootstrap GA Türleri

GA ile başlamadan önce, bootstrap gerçekleşmelerinin dağılımına göz atmak her zaman değerlidir. foo içinde hesaplanan hangi istatistiğe bakmak istediğimizi belirten index ile plot fonksiyonunu kullanabiliriz. Burada index=1 sepal uzunluğu ile genişliği arasındaki Spearman korelasyon katsayısı, index=2 sepal uzunluğunun medyanı ve index=3 sepal genişliğinin medyanıdır.

plot(myBootstrap, index=1)
bootstrap gerçekleşmelerinin dağılımı

Bootstrap korelasyon katsayılarının dağılımı oldukça normale benziyor. Bunun için GA bulalım. boot.ci kullanabiliriz. Varsayılan olarak yüzde 95 GA verir, ancak conf parametresiyle değiştirilebilir.

 boot.ci(myBootstrap, index=1)

## Warning in boot.ci(myBootstrap, index = 1): bootstrap variances needed for
## studentized intervals

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, index = 1)
##
## Intervals :
## Level      Normal              Basic         
## 95%   (-0.3178, -0.0209 )   (-0.3212, -0.0329 )  
##
## Level     Percentile            BCa          
## 95%   (-0.3007, -0.0124 )   (-0.3005, -0.0075 )  
## Calculations and Intervals on Original Scale

boot.ci 5 tür bootstrap GA sağlar. Bunlardan biri olan studentize aralık özeldir. Bootstrap varyansı tahminine ihtiyaç duyar. Bunu sağlamadık, bu yüzden R bir uyarı yazdırdı: bootstrap variances needed for studentized intervals. Varyans tahminleri ikinci düzey bootstrap ile veya (daha kolay) jackknife tekniğiyle elde edilebilir. Bu konu bu eğitimin kapsamının biraz dışında, bu yüzden kalan dört tür bootstrap GA’ya odaklanalım.

Hepsini görmek istemezsek, type argümanında ilgili olanları seçebiliriz. Olası değerler norm, basic, stud, perc, bca veya bunların bir vektörüdür.

 boot.ci(myBootstrap, index=1, type=c('basic','perc'))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, type = c("basic", "perc"), index = 1)
##
## Intervals :
## Level      Basic              Percentile     
## 95%   (-0.3212, -0.0329 )   (-0.3007, -0.0124 )  
## Calculations and Intervals on Original Scale

boot.ci fonksiyonu, adı... (evet, tahmin ettiniz!) bootci olan bir sınıfın nesnesini oluşturur. Öğeleri, type argümanında kullanılan GA türleriyle aynı adları taşır. $norm, güven düzeyini ve GA sınırlarını içeren 3 elemanlı bir vektördür.

 boot.ci(myBootstrap, index=1, type='norm')$norm

##      conf                       
## [1,] 0.95 -0.3177714 -0.02087672

$basic, $stud, $perc ve $bca ise GA’yı hesaplamak için kullanılan yüzdelikleri de içeren 5 elemanlı vektörlerdir (buna birazdan döneceğiz):

 boot.ci(myBootstrap, index=1, type='basic')$basic

##      conf                                    
## [1,] 0.95 975.98 25.03 -0.3211981 -0.03285178

Biraz gösterim (üzgünüm!)

Farklı GA türlerinin ne olduğunu hissetmek için biraz gösterim tanıtmamız gerekiyor. Şöyle ki:

  • t⋆ bootstrap tahmini (bootstrap gerçekleşmelerinin ortalaması) olsun,
  • t0 orijinal veri setindeki istatistik değerimiz olsun,
  • se⋆ bootstrap tahmininin standart hatası olsun,
  • b bootstrap tahmininin yanlılığı olsun, b = t⋆ − t0
  • α güven düzeyi olsun, tipik olarak α = 0.95,
  • zα standart normal dağılımın $1-\frac \alpha 2$-kantili olsun,
  • θα bootstrap gerçekleşmelerinin dağılımının α-yüzdeliği olsun.

Yüzdelik (Percentile) GA

Yukarıdaki gösterimle, yüzdelik GA şöyledir:

(θ(1 − α)/2, θ1 − (1 − α)/2)


yani sadece ilgili yüzdelikleri alır. Başka bir şey yok.

Normal GA

Tipik bir Wald GA şöyle olurdu:

t0 ± zα ⋅ se⋆

ancak bootstrap durumunda, bunu yanlılık için düzeltmeliyiz. Böylece şu hale gelir:

$$ t_0 - b \pm z_\alpha \cdot se^\star \\ 2t_0 - t^\star \pm z_\alpha \cdot se^\star$$

Temel (Basic) GA

Yüzdelik GA genel olarak önerilmez çünkü tuhaf kuyruklu dağılımlarda kötü performans gösterebilir. Temel GA (aynı zamanda pivotal veya empirical GA olarak da adlandırılır) buna çok daha dayanıklıdır. Bunun mantığı, her bir bootstrap yinelemesi ile t0 arasındaki farkları hesaplamak ve bunların dağılımının yüzdeliklerini kullanmaktır. Tüm ayrıntılar örneğin L. Wasserman’ın All of Statistics kitabında bulunabilir.

Temel GA için nihai formül şöyledir:

(2t0 − θ1 − (1 − α)/2, 2t0 − θ(1 − α)/2)

BCα GA

BCα, bias-corrected, accelerated ifadesinden gelir. Formülü çok karmaşık olmasa da biraz sezgisizdir, bu yüzden atlayacağım. Ayrıntılar ilginizi çekiyorsa Thomas J. DiCiccio ve Bradley Efron’un makalesine bakın.

Yöntemin adında geçen hızlandırma, bootstrap gerçekleşmelerinin belirli yüzdeliklerinin kullanılmasını gerektirir. Bazen bunlar uç yüzdelikler, muhtemel aykırı değerler olabilir. BCα bu gibi durumlarda kararsız olabilir.

Petal genişliği medyanı için BCα GA’ya bakalım. Orijinal veri setinde bu medyan tam olarak 3’tür.

 boot.ci(myBootstrap, index=3)

## Warning in boot.ci(myBootstrap, index = 3): bootstrap variances needed for
## studentized intervals

## Warning in norm.inter(t, adj.alpha): extreme order statistics used as
## endpoints

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, index = 3)
##
## Intervals :
## Level      Normal              Basic         
## 95%   ( 2.939,  3.046 )   ( 2.900,  3.000 )  
##
## Level     Percentile            BCa          
## 95%   ( 3.0,  3.1 )   ( 2.9,  2.9 )  
## Calculations and Intervals on Original Scale
## Warning : BCa Intervals used Extreme Quantiles
## Some BCa intervals may be unstable

BCα GA’sı (2.9, 2.9) elde ediyoruz. Bu tuhaf, ama neyse ki R bize extreme order statistics used as endpoints diye uyardı. Burada tam olarak ne olduğuna bakalım:

 plot(myBootstrap, index=3)
t'nin histogramı

Bootstrap gerçekleşmelerinin dağılımı alışılmadık. Büyük çoğunluğu (yüzde 90’dan fazlası) 3’tür.

 table(myBootstrap$t[,3])

##
##  2.9 2.95    3 3.05  3.1 3.15  3.2
##    1    1  908   24   63    1    2

Bu gibi bir durumda, BCα yöntemindeki “hızlandırma” kolayca uç değerlere “sıçrayabilir”. Burada, BCα GA’sının temelini oluşturan yüzdelik GA (3, 3.1)’dir. Şimdi, bunu “sola” doğru genişlettiğimizde 0,002-yüzdelik kullanmamız gerekir. “Sağa” genişleme de uçlara gider: 0,997-yüzdelik.

Sonuçların Yeniden Üretilmesi

Bazen bootstrap yinelemelerini yeniden oluşturmak gerekir. Bunu yapmak için R’yi kullanabiliyorsak sorun yoktur. set.seed fonksiyonu bu sorunu çözer. Bootstrap yinelemelerini başka yazılımlarda yeniden oluşturmak çok daha karmaşık olacaktır. Dahası, büyük R için, R’de yeniden hesaplama da bir seçenek olmayabilir (örneğin zaman kısıtı nedeniyle).

Bu sorunun üstesinden, her bir bootstrap örneklemini oluşturan orijinal veri seti öğelerinin indekslerini kaydederek gelebiliriz. boot.array fonksiyonu ( indices=T argümanıyla) bunu yapar.

 tableOfIndices<-boot.array(myBootstrap, indices=T)

Her satır bir bootstrap örneklemidir. Örneğin, ilk örneklemimiz şu öğeleri içerir:

 tableOfIndices[1,]

##   [1] 109  12 143  65  41  28 105  62  23 102  37  34  55  53  38   3  11
##  [18] 146  31 122  85 141 118 116 117  78 100  13  98  49  59  88  24  56
##  [35] 136   4  59  67 126 118 104 101  70  13  12  19  21 149  73 133  67
##  [52]  86   6  88 131 105 121 145 121  83  70  68  71 111  76  73 122 116
##  [69]  85 144 114 139  89 124  64  27  81  78  58  39  17  50  37 117  76
##  [86]  97 114  64  17  93 141  65 124  80 137  54  37  57  70 128  10  55
## [103]  13  53  59  45 116  14  77 118 108 138  50  78  49 104  49   1  85
## [120]  28  43  82  74  64  33  55  32  59  62 112   8  11  96  30  14  30
## [137]  38  85  66  85  97 107   4  18  76  35  31 133  27  69

indices=F (varsayılan) ayarlandığında, “Orijinal veri setinin her bir öğesi her bir bootstrap örnekleminde kaç kez göründü?” sorusunun cevabını alırız. Örneğin, ilk örneklemde: veri setinin 1. öğesi bir kez, 2. öğe hiç, 3. öğe bir kez, 4. öğe iki kez vb. görünmüştür.

 tableOfAppearances<-boot.array(myBootstrap)
tableOfAppearances[1,]

##   [1] 1 0 1 2 0 1 0 1 0 1 2 2 3 2 0 0 2 1 1 0 1 0 1 1 0 0 2 2 0 2 2 1 1 1 1
##  [36] 0 3 2 1 0 1 0 1 0 1 0 0 0 3 2 0 0 2 1 3 1 1 1 4 0 0 2 0 3 2 1 2 1 1 3
##  [71] 1 0 2 1 0 3 1 3 0 1 1 1 1 0 5 1 0 2 1 0 0 0 1 0 0 1 2 1 0 1 1 1 0 2 2
## [106] 0 1 1 1 0 1 1 0 2 0 3 2 3 0 0 2 2 0 2 0 1 0 1 0 0 1 0 2 0 0 1 1 1 1 0
## [141] 2 0 1 1 1 1 0 0 1 0

Bu tür bir tablo, bootstrap gerçekleşmesini R dışında yeniden oluşturmamıza olanak tanır. Ya da set.seed kullanmak istemediğimiz ve tüm hesaplamayı tekrar yapmak istemediğimizde R’de de.

 onceAgain<-apply(tableOfIndices, 1, foo, data=iris, cor.type='s')

Sonuçların aynı olup olmadığını kontrol edelim:

 head(t(onceAgain))

##             [,1] [,2] [,3]
## [1,] -0.26405188 5.70    3
## [2,] -0.12973299 5.80    3
## [3,] -0.07972066 5.75    3
## [4,] -0.16122705 6.00    3
## [5,] -0.20664808 6.00    3
## [6,] -0.12221170 5.80    3

head(myBootstrap$t)

##             [,1] [,2] [,3]
## [1,] -0.26405188 5.70    3
## [2,] -0.12973299 5.80    3
## [3,] -0.07972066 5.75    3
## [4,] -0.16122705 6.00    3
## [5,] -0.20664808 6.00    3
## [6,] -0.12221170 5.80    3

all(t(onceAgain)==myBootstrap$t)

## [1] TRUE

Evet, aynılar!

Eğer R’de Makine Öğrenimi hakkında daha fazla bilgi edinmek isterseniz, DataCamp’in Machine Learning Toolbox kursunu alın ve R ile Makine Öğrenimi başlangıç eğiticisine göz atın.

Konular

İlgili