Principal Component Analysis (PCA) R ile nasıl yapılır?


Egzersiz

PCA ile ilgili olarak unuttuklarımızı gözden gözgen geçirmek için basit bir egzersiz tasarladım. Bu örnek hipotetik bir veri seti üzerinde tasarlanmıştır. Durumu tam kavradıktan sonra bunu gerçek veri seti üzerinde uygulayacağız.

Memleket ve Dil Örneği

Diyelim ki bir grup öğrencinin farazi bir okulda Türkçe, Farsça, Arapça ve Süryanice dillerinden aldıkları dil puanları ve memleket bilgileri yer alsın. Bu veri setinde özellikle Van’lı olanları Türkçe ve İngilizcelerinin iyi Farsçalarının da azıcık iyi olmasına özen gösterdik. Mardinlilerin ise Arapçalarının iyi, Türkçe e İngilizcelerinin ise Van’lılara göre kötü olacağını düşündük. Daha sonra Siirtlileri ise Türkçe ve İngilizcelerinin kötü ancak Arapça’da en az Mardinliler kadar iyi olmalarını sağladık. Son olarak Süryanice’de en usta olanların Siirt’liler, Mardinliler ise azıcık Süryaniceye hakim olsun.

Şimdi PCA yaparak bu 5 dil yerine tek bir gösterge ortaya çıkararak tüm veri setinde memleket değişkenini en yüksek ölçüde açıklamayaçalışacağız.

Veri setimiz şuydu:

data = read.csv("ornek2.csv")
rownames(data) <- data$student
data[2:7]
##       memleket turkce farsca arapca ingilizce suryanice
## suat       van     90     60     10        90         0
## garb       van     88     70     11        88         0
## vahi       van     89     71     13        93         0
## mahi       mrd     70     10     90        66        11
## meli       mrd     65     11     98        67         7
## fahro      srt     40     40     89        60        90
## abdo       srt     40     40     89        60        90
## hako       srt     40     40     89        60        90

PCA Uygulaması

data.pca <- prcomp(data[3:7])
data.pca
## Standard deviations (1, .., p=5):
## [1] 62.697451 32.779140  3.550695  1.399083  0.912195
## 
## Rotation (n x k) = (5 x 5):
##                  PC1         PC2        PC3        PC4        PC5
## turkce    -0.3574760  0.11677660 -0.2578622  0.1165373 -0.8823267
## farsca    -0.2328702 -0.56494423  0.7378478 -0.1834692 -0.2202935
## arapca     0.6056701  0.49571857  0.5167242  0.1554323 -0.3102639
## ingilizce -0.2295579 -0.07426912  0.1540874  0.9442049  0.1628538
## suryanice  0.6312284 -0.64494065 -0.3135936  0.1925512 -0.2240209

PC1 adlı komponentimizin %62 oranında tüm veri setini açıkladığını görüyoruz, PC2 adlı komponent ise %32 oranında açıklıyor. Sadece bu iki komponenti kullarak tüm veri setini 62+32 = %94 oranında açıklayabiliriz. Bunun için 5 dille ilgili parametrelere gerek kalmaz.

#library("devtools")
#install_github("vqv/ggbiplot")

Görselleştirme

data.pca adlı değişkenimize tanımladığımız komponentlerden en güçlü iki komponentimizi grafiğe yerleştiriyoruz. Görülen okun yönü ise her bir parametrenin (dil) ekseninde yer aldığı komponente hizmet ettiğini gösteriyor, ok ne kadar uzunsa o parametre o komponente o kadar hizmet ediyor (açıklıyor). Buna göre yatay eksende PCA1 komponentini en çok arapça ve süryanice açıklıyor. Aynı komponenti Türkçe, İngilize ve Farça bilgisi de açıklıyor ancak ters yönde.

Öte yandan Süryanic bilgisi ile Farsça bilgisi ise PCA2’yi açıklıyor bu bilgilerin aksi istikamete ise Arapça yer alıyor. Tüm bu ilişkiler yani hangi komponentin hangi parametre ile ne kadar ilgili olduğu zaten PCA tablosunda görülebiliyor.

library("ggbiplot")
ggbiplot::ggbiplot(data.pca,labels= rownames(data),groups = data$memleket)

Bu grafiğe isimleri de eklediğimizde memleketlerine göre de renklendirdiğimizde sağ alt köşede üç mardinli arkadaşın da aybı yerde yer aldığını ve bu arkadaşların da bir küme olma nedeninin Süryanice bilgisi olduğunu başka bir şeye bakmaksızın anlıyoruz. Yani memlekete bilgisi olmasaydı dahi PCA bize bu arkadaşların ortak bir özelliğinin var olduğunu söylemiş oluyor. Aynı durum grafiğin üstünde siirtiler kaar olmasa da yine ayrı bir grubun var olduğunu gösteriyor. Solda is Van’lılar yer alıyor onlar ise Türkçe, İngilizce ve Farsça bilgisi ile göz dolduruyor 🙂

Sonuçlar

  • PCA ile çok fazla sayıda parametre daha az parametreye indirgenebilir. (PCA1, PCA2 gibi)
  • Bu indiregeme işlemi görselleştirildiğinde birbiri ile ilgili (benzer ya da tam tersi) parametreler (örneğin Türkçe, İngilizce ve Farsça’nın grubun bir kısmı için birbiri ile ilgili parametreler olması yani öğrencilerin memleketindeki şu veya bu nedenden ötürü Farsçası iyi olan kişinin İngilizcesinin de iyi olması gibi) ortaya çıkarılabilir.
  • GGBiplot aracı yardımı ile gözlemlerin (öğrencilerimiz) de grafiğe eklenmesi ile veri seti içinde gizli gruplar da ortaya çıkarılabilir. Örneğimizde memleket kolonu yer almasa bile biz bu öğrencilerinin bazılarının adını koymadığımız bir nedenden ötürü (memleket) grup olduklarını anlamış oluyoruz.
Reklamlar

Malthus Popülasyon Yasası Bir Örnek ve R Dili ile Modelleme (Adi Diferansiyel Denklemle)


Sistem dinamiklerinde ‘stocks’ olarak anılan her bir düğüm artış ve azalışa hazir olan bir olgudan başka bir şey değildir. Örneğin stock nüfus olduğunda nüfus için ‘inflow’ doğumlar, ‘outflow’ ölümlerdir. Amaç diyelim ki bugün nüfusunu bildiğimiz bir yerde 10 sene sonraki nüfusu hesaplamaktır. Bu hesaplama yapılırken elbette elimizde nüfusun nasıl artacağına dair bir fonksiyon olmalıdır. Bu fonksiyon literatürden veya uzmanlıktan gelen bir fonksiyon olabilir.

Girdilerimiz 1970 yılındaki nüfus: 3.7 milyar, artış oranımız: 0.01801. Başka bir şey yok. aux parametremize artış oranını, 1970 yılındaki nüfüsu sCustomers değişkenimize atadık. START<-1970; FINISH<-2008; STEP<-0.25 zaten kendi kendini açıklıyor. Daha sonra her yıl çeyreklik bazda (step) nüfus artışını ve grafikleri ortaya çıkarıyoruz.

Sistem dinamikleri elbette bundan ibaret değil. Malthus modeli gibi modellerin çıktılarını başka bir modele örneğin artan nüfusa göre gıda ihtiyacı fonksiyonuna bağlayarak fonkisyonlardan oluşan bir network yaratarak simülasyonlar yapıyoruz.

Malthus fonksiyonunun matematiksel gösterimi ve çözümü için şuraya bakabilirsin: https://math.stackexchange.com/questions/345242/differential-equation-word-problem-malthuss-law

Modeli Tasarlama

library(deSolve)
library(ggplot2)
require(gridExtra)
library(scales)

# Setup simulation times and time step
START<-1970; FINISH<-2008; STEP<-0.25

# Create time vector
simtime <- seq(START, FINISH, by=STEP)

Kademeli olarak çeyreklik bazda (0.25) yılları arttırdık

simtime
##   [1] 1970.00 1970.25 1970.50 1970.75 1971.00 1971.25 1971.50 1971.75
##   [9] 1972.00 1972.25 1972.50 1972.75 1973.00 1973.25 1973.50 1973.75
##  [17] 1974.00 1974.25 1974.50 1974.75 1975.00 1975.25 1975.50 1975.75
##  [25] 1976.00 1976.25 1976.50 1976.75 1977.00 1977.25 1977.50 1977.75
##  [33] 1978.00 1978.25 1978.50 1978.75 1979.00 1979.25 1979.50 1979.75
##  [41] 1980.00 1980.25 1980.50 1980.75 1981.00 1981.25 1981.50 1981.75
##  [49] 1982.00 1982.25 1982.50 1982.75 1983.00 1983.25 1983.50 1983.75
##  [57] 1984.00 1984.25 1984.50 1984.75 1985.00 1985.25 1985.50 1985.75
##  [65] 1986.00 1986.25 1986.50 1986.75 1987.00 1987.25 1987.50 1987.75
##  [73] 1988.00 1988.25 1988.50 1988.75 1989.00 1989.25 1989.50 1989.75
##  [81] 1990.00 1990.25 1990.50 1990.75 1991.00 1991.25 1991.50 1991.75
##  [89] 1992.00 1992.25 1992.50 1992.75 1993.00 1993.25 1993.50 1993.75
##  [97] 1994.00 1994.25 1994.50 1994.75 1995.00 1995.25 1995.50 1995.75
## [105] 1996.00 1996.25 1996.50 1996.75 1997.00 1997.25 1997.50 1997.75
## [113] 1998.00 1998.25 1998.50 1998.75 1999.00 1999.25 1999.50 1999.75
## [121] 2000.00 2000.25 2000.50 2000.75 2001.00 2001.25 2001.50 2001.75
## [129] 2002.00 2002.25 2002.50 2002.75 2003.00 2003.25 2003.50 2003.75
## [137] 2004.00 2004.25 2004.50 2004.75 2005.00 2005.25 2005.50 2005.75
## [145] 2006.00 2006.25 2006.50 2006.75 2007.00 2007.25 2007.50 2007.75
## [153] 2008.00
#Stock ve Aux değerlerini atadık. Müşteri sayısı 10.000
# Create stock and auxs
stocks  <- c(sCustomers=3712000000)
auxs    <- c(aGrowthFraction=0.018021)

Model basit bir fonksiyondur.

# Model function
model <- function(time, stocks, auxs){
  with(as.list(c(stocks, auxs)),{ 
    #Denklem: fRecruits işe değerimiz müşteri sayısı ile artış katsayısı
    #çarpımı kadardır
    fRecruits<-sCustomers*aGrowthFraction
    # Inflow
    dC_dt <- fRecruits
    return (list(c(dC_dt),
                 Recruits=fRecruits,
            GF=aGrowthFraction))   
  })
}
# Modeli ve değerleri Düzenli Difereansiyel Denkleme Sokalım
# Run simulation
o<-data.frame(ode(y=stocks, times=simtime, func = model, 
                  parms=auxs, method="euler"))

Model Sonuçları

Her yıl olacak işe alım ve kayıplar hesaplandı:

time<dbl>sCustomers<dbl>Recruits<dbl>GF<dbl>
1970.003712000000668939520.018021
1970.253728723488671953260.018021
1970.503745522319674980580.018021
1970.753762396834678021530.018021
1971.003779347372681076190.018021
1971.253796374277684144610.018021
1971.503813477892687226850.018021
1971.753830658563690322980.018021
1972.003847916638693433060.018021
1972.253865252464696557150.018021

Next12345616Previous1-10 of 153 rows

Modeli Görselleştirme

p1<-ggplot()+
  geom_line(data=o,aes(time,o$sCustomers,color="Customers"))+
  scale_y_continuous(labels = comma)+
  ylab("Stock")+
  xlab("Year") +
  labs(color="")+
  theme(legend.position="none")
p1

Görsel 1: Yıllara göre nüfusun artışı