Dünya Bankası veri tabanından enerji verileri çekmek ve otomatik model belirlenmesi¶
Barış Sanlı, www.barissanli.com barissanli2@gmail.com
Niye Excel değil de R? Gerçekte R veya Python veya başka bir dil tercihim özel olarak yok. Sadece bir dili kullananların sayısı ne kadar çok ise o kadar fazla soruya daha hızlı cevap bulabiliyoruz.
Excel değil de R sorusu "elma" ile "armut" kıyaslaması gibi gelebilir. Fakat bu derste Excel ile veri temizleme hazırlama, model seçimi gibi bir çok sorunun R ile ne kadar hızlı halledilebileceğini göreceğiz.
Bunun için Dünya Bankası veri bankasını kullanacağız. Sırası ile:
- Dünya Bankası veri setine erişmek için "wbstats" kütüphanesini kuracağız
- Veri kütüphanesinden, kişi başı gelir, nüfus ve enerji göstergelerini çekeceğiz
- Bu değişkenlerin birbiri ile ilintisine bakacağız
- Grafikleyeceğiz, hatta "The Economist" temasını kullanacağız
- Son olarak da otomatik doğrusal regresyon modeli belirleyerek, modeli kuracağız
"wbstats" kütüphanesi¶
wbstats, Dünya Bankası veri tabanına erişmeyi kolaylaştıran bir R kütüphanesi. Daha önce NOAA ve EIA veritabanlarında görüldüğü gibi API anahtarı almaya gerek yok.
Doğrudan veri talep ediyoruz.
Eğer kütüphane daha önce kurulmadıysa:
install.packages("wbstats")
ile bu kütüphaneyi kurabiliriz. Açıklamalar konusunda wbstats rehberine bakabilirsiniz. https://cran.r-project.org/web/packages/wbstats/wbstats.pdf
Eğer kütüphaneyi kurarken bir sorun yaşamadıysanız, her zaman olduğu gibi library komutu ile bu kütüphaneyi devreye alabiliriz.
Verileri çekerken ise sadece göstergelerin kodlarını vektör olarak indicator parametresine girmemiz yeterli. Bu çalışmada benim araştırma sorum:
"Küresel olarak elektrikte kömür kullanımı hangi parametrelerle nasıl ilişkili, kömürü az kullanmak bir kalkınmışlık göstergesi mi?"
sorusu. Dünya Bankasındaki tüm ülkelerin yıllara sair verileri üzerinden kömür, nüfus, ekonomik büyüme ve diğer kaynaklar ile ilişkisini bulmaya çalışacağız.
Önce Kütüphanemizi yükleyelim:
library(wbstats)
Aradığımız Veri kodlarını bulmak¶
Herhangi bir veriyi sistemden çekebilmek için veri kodunu bilmemiz gerekiyor. Mesela elektrik üretimindeki kömürün oranı "EG.ELC.COAL.ZS" gibi biraz garip bir kod. Dolayısıyla önce veritabanında doğru kodları bulmamız gerekiyor.
Bunun için wbsearch komutunu kullanacağız. Komuta aradığımız anahtar kelimeleri yazınca bizlere o konu ile ilgili tüm kodları getirecek. Burada indicatorID sütunu altında yer alan kod ise sistemden talep edeceğimz verinin kodu.
Şimdi elektrik ile ilgili veri kodlarını bulmaya çalışalım
wbsearch("electricity production")
Dünya Bankası Web sayfasından verileri bulmak¶
Dünya Bankası kütüphaneleri ile ilgili İngilizce sitelere de bakmak faydalı olabilir. Mesela : https://cengel.github.io/gearup2016/worldbank.html gibi.
Ana olarak WorldBank data catalog: https://data.worldbank.org/ adresinden erişilebilir. Burada zaten bir kaç kelime yazınca verinin adı otomatik geliyor. Veriye tıklayınca da açılan sayfanın adres kısmında verinin ismi gözüküyor.
Örneğin: Nüfus için https://data.worldbank.org/indicator/SP.POP.TOTL adresinin sonundaki "SP.POP.TOTL" veri kodudur. Benim hangi verileri kullanacağımı bulmak haftalar sürdü, sonra kodlarını bulmak dakikalar. Çünkü Dünya Bankası veri seti çok geniş. Bunun sonucunda:
- Kişi başı satın alma pariteli GSYİH "NY.GDP.PCAP.PP.KD"
- Nüfus "SP.POP.TOTL"
- Elektrik üretiminde doğalgaz oranı "EG.ELC.NGAS.ZS"
- Elektrik üretiminde kömür oranı "EG.ELC.COAL.ZS"
- Elektrik üretiminde yenilenebilir oranı "EG.ELC.RNEW.ZS"
- Toplam elektrik üretimi "4.1.1_TOTAL.ELECTRICITY.OUTPUT"
- Elektriğe şehirde erişim oranı "EG.ELC.ACCS.UR.ZS"
- Elektriğe kırsalda erişim oranı "EG.ELC.ACCS.RU.ZS"
- Toplam nihai enerji tüketimi "1.1_TOTAL.FINAL.ENERGY.CONSUM"
- Yenilenebilirin toplam enerji tüketimindeki oranı "2.1_SHARE.TOTAL.RE.IN.TFEC"
Tabii ki istediğim başka veriler de oldu, ama bu ders için çekemedim. Bunun veri çekme sınırından mı, verinin kendisinden mi kaynaklandığı başka bir problem.
Yukarıdaki verileri alarak wb_dat değişkenine alıyoruz. Verilerin gelmesi zaman alabilir, biraz bekleyin, veriler gelince de gelen verinin başlıklarına names komutu ile bir bakın
wb_dat <- wb(indicator = c("NY.GDP.PCAP.PP.KD", "EG.ELC.NGAS.ZS","EG.ELC.COAL.ZS",
"EG.ELC.ACCS.RU.ZS","EG.ELC.ACCS.UR.ZS", "EG.ELC.RNEW.ZS",
"2.1_SHARE.TOTAL.RE.IN.TFEC" ,
"1.1_TOTAL.FINAL.ENERGY.CONSUM", "SP.POP.TOTL",
"4.1.1_TOTAL.ELECTRICITY.OUTPUT" , "12.1_TD.LOSSES ",
"4.1.2_REN.ELECTRICITY.OUTPUT ", "11.1_THERMAL.EFFICIENCY "
))
names(wb_dat)
Veri setimizin üst kısmına bir bakalım. Neler varmış? Bir de ben internetten çektiğim veriyi daima başka bir değişkene de (wb2) aktarıyorum ki mevcut veriyi değiştirirken yanlışlık yaparken orjinal veriyi tekrar internete bağlanmadan yükleyebileyim.
head(wb_dat)
wb2<-wb_dat
Veri ile işlem yapmak ve sütunlara bölmek¶
Bir de Dünya Bankası veri tabanında yer alan ülkelerin isim ve bilgilerini alalım. Çünkü gelen veri seti, çok kompakt, yani kaç veri istersek isteyelim değerler hep "value" sütununda, veri ismi de "indicator" sütununda. Burada değişkenleri sütunlara yaymamız gerekecek.
Önce sistemdeki tüm ülkelerin isimlerini wbcountries() komutu ile çekerek başlıklarına names komutu ile bakalım.
wb_countries <- wbcountries()
names(wb_countries)
Verilerdeki bölge değişkenlerini iso2c altında birleştirelim
wb_dat <- merge(wb_dat, y = wb_countries[c("iso2c", "region")], by = "iso2c", all.x = TRUE)
wb_dat <- subset(wb_dat, region != "Aggregates") # this also removes NAs
Şimdi bu karışık kodlara ait verileri yeni isimlendirdiğimiz sütunlara çekelim. Öncelikle "NY.GDP.PCAP.PP.KD" kodunu GDP olarak değiştirelim. Bu yöntem yerine ayırıp sonra sütun isimlendirme de yapılabilir. Bunu web'den aldığım bir koddan değiştirerek uyguladım, daha güvenli bir yöntem olabilir.
Önce kodları isimlendirelim
wb_dat$indicatorID[wb_dat$indicatorID == "NY.GDP.PCAP.PP.KD"] <- "GDP"
wb_dat$indicatorID[wb_dat$indicatorID == "EG.ELC.NGAS.ZS"] <- "gas"
wb_dat$indicatorID[wb_dat$indicatorID == "EG.ELC.COAL.ZS"] <- "coal"
wb_dat$indicatorID[wb_dat$indicatorID == "EG.ELC.RNEW.ZS"] <- "renewables"
wb_dat$indicatorID[wb_dat$indicatorID == "EG.ELC.ACCS.UR.ZS"] <- "urbanaccess"
wb_dat$indicatorID[wb_dat$indicatorID == "EG.ELC.ACCS.RU.ZS"] <- "ruralaccess"
wb_dat$indicatorID[wb_dat$indicatorID == "SP.POP.TOTL"] <- "population"
wb_dat$indicatorID[wb_dat$indicatorID == "4.1.1_TOTAL.ELECTRICITY.OUTPUT"] <- "electricityGWh"
wb_dat$indicatorID[wb_dat$indicatorID == "2.1_SHARE.TOTAL.RE.IN.TFEC"] <- "reinTFEC"
wb_dat$indicatorID[wb_dat$indicatorID == "1.1_TOTAL.FINAL.ENERGY.CONSUM"] <- "TFEC"
Eğer yüklü değil ise "reshape2" kütüphanesini kurmamız gerekir bunun için de
install.packages("reshape2")
Ben de zaten kurulu olduğundan reshape2 kütüphanesinden dcast komutu ile indicatorID sütunundaki daha önce isimlendirdiğim tüm verileri karşılık gelen değerleri ile ayrı sütunlara çeviriyorum, yeni veri matrisini tekrar wb_dat olarak adlandırıyorum.
library(reshape2)
wb_dat <- dcast(wb_dat, iso2c + country + date ~ indicatorID, value.var = 'value')
Şimdi veriye baktığımız zaman, head(wb_dat) ile, ayrı sütunlarda her bir ülke için 1960 dan buyana tüm verilerin gruplandığını göreceğiz.
head(wb_dat)
Şimdi de veri setimizin bir özetine bakalım, görüldüğü üzere 12377 adet veri satırından oluşmakta, ama bunların bir kısmı boş ya da eksik veri. Ayrıca burada verilerimizin hangi aralıkta yer aldığı, ortamalarının ne olduğu da rahatlıkla görülebilir.
summary(wb_dat)
Veri setinden bir kısım seçmek¶
Hazır veri setimiz tam iken burada bir komut ile değişkenlerimizi diske saklayalım. Bunun için saveRDS komutu yeterli olacaktır. Komutun kullanımı çok basit, doğrudan komutu yazıp, değişkeni ve dosya ismini yazıp saklıyoruz.
saveRDS(wb_dat,"ulkeverileri.rds")
Tüm veri setimiz görüldüğü üzere oldukça geniş. Ben analiz için daha çok 2000 yılı sonrası verilerini kullanmak istiyorum, bu sebeple wb_dat verisinden 2000 yılı sonrası verilerini seçerek, wb_dat verisine geri yüklüyorum. Artık wb_dat da sadece 2000 yılından sonraki veriler var.
# yıl veya değer seçelim
wb_dat<-subset(wb_dat,date>2000)
Şimdi de yeni bir veri çerçevesi oluşturarak, sütunları tekrar isimlendirelim. Bu işleme gerek olmayabilir, ama olaki siz de benim gibi belki bir çok veri tabanından veri çekerseniz bu şekilde bir birleştirme gerekebilir. colnames komutu ile de sütunlara isimlerini verelim.
data_set isimli verisetimize baktığımız zaman ise bir sürü olmayan veri (NA)lar göreceğim. Hatırlarsanız yukarıda Andorra'nın 1960'dan başlayan veri setleri yoktu. En üstte hep bu ülkenin 1960lı yıllar verileri olduğundan bu veriler olmayabilir.
data_set<-data.frame(wb_dat$coal,wb_dat$gas, wb_dat$GDP,
wb_dat$ruralaccess, wb_dat$urbanaccess, wb_dat$renewables, wb_dat$TFEC, wb_dat$reinTFEC,
wb_dat$electricityGWh, wb_dat$population)
colnames(data_set)<-c("coal","gas","gdp","raccess","uaccess", "ren","TFEC","renTFEC","elec","pop")
head(data_set)
NA olmayan verileri bulmak¶
Normalde NA, yani olmayan verileri, birçok yöntemle elimine edebiliriz. Burada ben herhangi bir satırda NA var ise o satırın çıkmasını isteyeceğim. Bunun için de complete.cases komutunu kullanıyoruz. Tüm sütunlarında veri olan satırların numarasınız complete.cases(data_set) ile nasiz_veriler değişkenine atarak başlığına bakıyoruz.
nasiz_veriler<-data_set[complete.cases(data_set),]
head(nasiz_veriler)
Yeni sütunlar eklemek veya hesaplamak¶
Bu bilgisayar işlerinin sağı solu belli olmadığından, ben "nasiz_veriler" değişkeni üzerinde değil, değişkeni kopyaladığım nasiz_veriler2 değişkeni üzerinden işleme devam edeceğim. Bu çalışmalar esnasında genelde bir R defteri haftalarca değiştirildiğinden ben her aşamada veri yedekliyorum.
Ardından GDP ve TFEC ve daha sonra elec verilerinin logarithmalarını hesaplayacağım için bu sütunlarda 1'den büyük değerleri alıyorum. 0'dan büyük değerleri almam matematiksel olarak yeterli ama 1'den büyük alarak güvenilir tarafta kalıyorum. 3.sütun GDP, 7.sütun da TFEC(total final energy consumption, nihai enerji tüketimi). Bu iki sütunun da 1'den büyük olduğu satırları alıyorum.
Kodlama stili olarak ben genelde KISS, kullanıyorum, yani zekice birşey yapmamaya çalışıyorum. Sadece adım adım neyi değiştirince ne olduğu, nasıl geliştiğini görmeye çalışıyorum. Burada c(3,7) ve & ile mantıksal operatör de kullanılabilir belki. Ama ben bildiğim ve hatayı izleyebileceğim yolu tercih ediyorum
nasiz_veriler2<-nasiz_veriler
nasiz_veriler2<-nasiz_veriler2[nasiz_veriler2[,7]>1.0,]
nasiz_veriler2<-nasiz_veriler2[nasiz_veriler2[,3]>1.0,]
Şimdi de sütunlar hesaplayalım. R'da veri çerçevesinde yeni bir sütun eklemek kolay. Sadece yeni sütun ismini yazarak, içinde ne olması gerektiğini belirtiyorsunuz. İlk satırda görüleceği üzere, nasiz_veriler2 'deki TFEC sütununun logarithmasını alarak logTFEC diye bir sütun oluşturuyorum. 3.satırda 1 milyon ile çarpma sebebi, değişkenin GWh olması ben ise birey başına kWh elektrik üretimi vs hesaplayacağım. Bu sebeple GWh'i kWh'e çeviriyorum.
Tüm hesaplamalardan sonra değişkene tekrar bakalım
nasiz_veriler2$logTFEC<-log(nasiz_veriler2$TFEC)
nasiz_veriler2$logGDP<-log(nasiz_veriler2$gdp)
nasiz_veriler2$elecPerPop<-(nasiz_veriler2$elec/nasiz_veriler2$pop)*1000000
nasiz_veriler2$TFECPerPop<-nasiz_veriler2$TFEC/nasiz_veriler2$pop
nasiz_veriler2$logelec<-log(nasiz_veriler2$elec)
head(nasiz_veriler2)
Verilerin birbiri ile ilişkisine bakmak¶
İlinti ya da korelasyonlara bakmak için veri çerçevesini cor komutuna vermemiz yeterli. Her sütunun diğer sütunlar ile ilişkisine bakan cor komutu ile hangi değişken hangi diğer değişken ile ne kadar ilintili bulmak mümkün
cor(nasiz_veriler2)
Korelasyon tablosunu incelemek¶
Benim çalışmadaki amacım, herkesin iddia ettiği gibi petrolün, kömürün sonunun gelip gelmediği ve hangi değişkenlerden etkilendiği. "Mesela kömür gelişmede bir adım mıdır?" sorusuna cevap da, ilinti oranı çok düşük 0.02'tür. Yani geliştikçe daha az kömür tüketilir diye birşey söylenemez. Aksine hem ren yani elektrikteki yenilenebilir oranı ve renTFEC nihai enerji talebindeki yenilenebilir oranının katsayıları (negatif)'dir. Yani gelişmiş ülkeler yenilenebiliri sanki daha az kullanan ülkelerdir. Yalnız ilintide sebep sonuç mümkün değildir. Gelişmişlik mesela en çok kişi başı elektrik, enerji ve kırsal da elektriğe erişim ile ilintilidir.
Peki kömür kullanımına dair bir doğrusal regresyon yapsak hangi değişkenlerden kurmamız gerekir.
Bunun için otomatik bir model seçimi yapacağız.
Otomatik regresyon modeli seçme¶
Dersin bu kısmı biraz İLERİ SEVİYE. Eğer zor geliyorsa, bir sonraki kısımda grafikleme kısmına geçebilirsiniz.
Tüm değişken ve ilinti tablolarını oluşturduktan sonra, ülkelerin kömür kullanım oranları hangi değişkenlerle ilgili bunu bulmak için "leaps" paketini kullanabiliriz. Bunun için önce paketi yüklemek gereki
- install.packages("leaps")
Paketi yükledikten sonra asıl komutumuz regsubsets. Ben aşağıda bu kütüphane için kullandığım web adreslerini de verdim.
- regsubsets(coal ~ . , data=nasiz_veriler2 şeklinde devam eden aşağıdaki komutta kömür değişkeninin tüm veri setindeki diğer değişkenlerden hangilerine bağlı olduğunu parametrelerle hesaplatıyoruz.
#--------------------------------------
#best model selection for coal
#-------------------------------------
# web: https://rstudio-pubs-static.s3.amazonaws.com/2897_9220b21cfc0c43a396ff9abf122bb351.html
# install.packages("leaps")
require(leaps)
results<-regsubsets(coal ~ ., data = nasiz_veriler2, nvmax = 1000, nbest=1)
Şimdi ise sonuçları görelim. Öncelikle:
- summary(results) komutu ile sonuçların özetini alalım
- which.max(results.summary$adjr2) komutu ile Adjusted R2 (düzeltilmiş r-kare)'ye göre em iyi sonucu bul diyoruz
- results.by.adjr2 komutu ile hangi denklemin en iyi sonucu verdiğini ekrana yazıyoruz
- plot komutu ile de hangi değişkenler ile en yüksek R2 elde edildiğini çizdiriyoruz
results.summary<-summary(results)
results.by.adjr2 <- which.max(results.summary$adjr2)
results.by.adjr2
results.summary$adjr2[169]
results.summary$which[results.by.adjr2,]
Yukarıda kodlama açısından tüm değişkenleri ekleyerek sadece TFEC değişkenini kabul etmediğini gördük. TFEC değişkeninşi kabul etmeden elde edilen en yüksek yakınsama değeri 0.445'tir. Bunu da ekrana yazdırdıktan sonra, sonuçları ekrana çizelim.
Birşeye dikkat etmekte fayda var, otomatik seçimde modelin doğruluğu değil belirtilen yakınsamaya bakılıyor. Yani R2 veya BIC değerlerini minimum veya maksimum yapan değerler seçiliyor.
plot(results , scale = "adjr2", main = "Adjusted R^2")
Görüldüğü üzere en yüksek R2, TFEC değişkeni dışarıda bırakıldığı durumda elde ediliyor.
Verileri Grafiklendirmek için ggplot2 ve ggthemes¶
ggplot2 kütüphanesi grafikleme açısından en popüler kütüphanelerden biridir. Ben buna ek olarak "ggthemes" kütüphanesini de kullanıyorum. Arada raporlar daha iyi gözüksün diye bu kütüphanedeki "The Economist" veya "Wall Street Journal" temalarını seçiyorum.
Aşağıdaki grafikleme de "The Economist"'in beyaz temasını kullanacağız.
Yoksa, önce ggplot2 ve ggthemes'i kuralım
install.packages("ggplot2")
install.packages("ggthemes")
Şimdi de kurduğumuz kütüphaneleri kullanalım
library(ggplot2)
library(ggthemes)
Şimdi öncelikle veri setimizde gaz ve kömürün elektrikte payı 0'dan, kişi başı elektrik kullanımı da 1000 kWh'dan büyük olan veri satırları için işlem yapacağız. Bunun için subset komutunu kullanıyoruz. Filtreleme için ise "gas>1", mantıksal operatör olarak da "&" kullanılıyor. Yani 3 değişkende de belirli bir sınırın üzerini istiyoruz
2.satırda bu veri kümesinde,
- x ekseninde toplam elektrik tüketiminin logaritması
- y ekseninde kömürün elektrik üretimindeki payı
- renk farklılaştırması için de GSYİH(gayri safi yurt içi hasıla)'yı kullanıyoruz. Daha zengin daha açık renk
3.satırdaki theme_economist_white() komutu ile de meşhur derginin beyaz temasını seçiyoruz.
ggplot(subset(nasiz_veriler2, gas>1 & coal>1 & elecPerPop>1000),
aes(x = log(elec), y = coal, colour= logGDP) ) +
theme_economist_white() + geom_point()
Şimdi bir tane de elektrik tüketimine bakmadan gelişmişlik ile Toplam Nihai Enerji Tüketimindeki yenilenebilir oranına bakalım. Renk skalası da logaritmik toplam enerji tüketimi olsun
#ggplot(subset(wb_dat, gas>1 & coal>1 ), aes(x = log(GDP, y = log(TFEC), colour=ruralaccess ) ) + geom_point()
ggplot(subset(nasiz_veriler2, gas>1 & coal>1 ),
aes(x = logGDP, y = renTFEC, colour= logTFEC) ) +
theme_economist_white() + geom_point()
Sonuç¶
Bu derste daha karışık daha profesyonel bir analiz denedik. Herşeyi anlamak kolay olmayabilir. Fakat kısaca neler yaptığımıza bakalım
- Dünya Bankasından veri almayı gördük
- Bu verileri genişleterek düzgün sütunlu veri çerçeveleri oluşturduk
- Verilerin bibirbir ile ilişkilerinin incelendiği ilinti matrisi hesaplandı
- "leaps" paketi ile otomatik model değişkenleri seçtik
- Grafikleme ve profesyonel tema kullanımını inceledik
Sanıyorum bu kısımda bir çok soru oluşmuştur. Daha kısa kodlarla bu egzersiz yapılabilirmiydi? Sanırım evet Önerilerinizi bekliyorum barissanli2@gmail.com
EK¶
Otomatik model seçimi konusunda meraklı olanlar için aşağıda bir de meifly kütüphanesi ile otomatik model seçimi denemesi var:
#install.packages("meifly")
library(meifly)
# coal~gas+raccess+uaccess+ren+pop+logelec,data=nasiz_veriler2
# https://rstudio-pubs-static.s3.amazonaws.com/2897_9220b21cfc0c43a396ff9abf122bb351.html
fitall.out <- fitall(y = nasiz_veriler2$coal,
x = nasiz_veriler2[,c(2,3,4,5,6,7,8,9,10,11,12,13)])
fitall.out.aic <- t(sapply(fitall.out, extractAIC))
## Create an order list of increasing AIC
final.out.order <- order(fitall.out.aic[,2])
## Show the result for the best model
fitall.out[final.out.order][1]
summary(fitall.out[final.out.order[1]])