R Dersleri - Ali Onur Gitmez

Ana Sayfa Ders 1 Ders 2 Ders 3 Ders 4 Ders 5 Ders 6 Ders 7 Ders 8 Ders 9 Ders 10 Ders 11 Ders 12 Ders 13 Ders 14 Ders 15

R Dersi 15

Computational Methods

Bu son kisimdaysa gecen hafta causal inference’a giris yaptigimiz gibi sosyal bilimlerde sikca kullanilan computational metodlara bir giris yapip orneklerine bakacagiz. Bu yontemler

• Web Scraping
• Text Analysis
• Sentiment Analysis
• Network Analysis
• Spatial Analysis
• Machine Learning

olacak. Bunlarin hepsi ozellikle son donemlerde yaygin bir sekilde sosyal bilimler alaninda kullaniliyor ve calisilan konuya gore de onemleri yuksek. Bundan dolayı uygulamalı olarak temel bir giriş yapmak önemli olacaktır.

Diğer paketler

library(tidyverse)
library(hms)
library(stringr)
library(writexl)
library(readxl)
library(lubridate)

Scraping paketleri

library(rvest)
library(xml2)
library(XML)
library(httr)

Text paketleri

library(tidytext)
library(stopwords)
library(ggwordcloud)
library(seededlda)
library(glmnet)
library(caret)

Web Scraping

Bir siteden veri cekmek icin once robots.txt kismina bakariz. Orada hangi sayfaların çekmeye uygun olup olmadığını görebiliriz. Bir örnek bile bunu deneyelim.

Sabah Gazetesi’nin sayfasından verileri çekeceğimizi varsayalım. Bunun için öncelikle robots.txt kısmına bakarız. Bunun sonrasında verileri çekebileceğimiz sayfayı bulmamız lazım. Bunun için ya sayfa içinde arama yapmak bir çözüm olabilir. Verileri çekmek için farklı yöntemler kullanabiliriz. Bunlardan ilki HTML ile çekmek olacak. Neredeyse websitesi HTML formatı üstünde bir yapıya kurulu olur ve bunu kullanarak site içinden resim, yazı veya tablo çekebilirsiniz.

h1 ve p gibi kısımlara node deriz. Bunlar veriyi çekmemiz için en önemli kısımlar çünkü bütün veriler bunların altında saklanıyor. Sayfa yapısını inspect ederek veya belli uzantılar aracılığıyla bakarak neyin hangi yapıda olduğunu anlayabiiriz. Şimdi bir örnek ile görelim.

url <- "http://books.toscrape.com/index.html"
webpage <- read_html(url)
titles <- webpage %>%
  html_nodes("h3 a") %>%
  html_attr("title")
titles

##  [1] "A Light in the Attic"                                                                          
##  [2] "Tipping the Velvet"                                                                            
##  [3] "Soumission"                                                                                    
##  [4] "Sharp Objects"                                                                                 
##  [5] "Sapiens: A Brief History of Humankind"                                                         
##  [6] "The Requiem Red"                                                                               
##  [7] "The Dirty Little Secrets of Getting Your Dream Job"                                            
##  [8] "The Coming Woman: A Novel Based on the Life of the Infamous Feminist, Victoria Woodhull"       
##  [9] "The Boys in the Boat: Nine Americans and Their Epic Quest for Gold at the 1936 Berlin Olympics"
## [10] "The Black Maria"                                                                               
## [11] "Starving Hearts (Triangular Trade Trilogy, #1)"                                                
## [12] "Shakespeare's Sonnets"                                                                         
## [13] "Set Me Free"                                                                                   
## [14] "Scott Pilgrim's Precious Little Life (Scott Pilgrim #1)"                                       
## [15] "Rip it Up and Start Again"                                                                     
## [16] "Our Band Could Be Your Life: Scenes from the American Indie Underground, 1981-1991"            
## [17] "Olio"                                                                                          
## [18] "Mesaerion: The Best Science Fiction Stories 1800-1849"                                         
## [19] "Libertarianism for Beginners"                                                                  
## [20] "It's Only the Himalayas"

prices <- webpage %>%
  html_nodes(".price_color") %>%
  html_text() %>%
  sub("£", "", .) %>% 
  as.numeric()   
prices

##  [1] 51.77 53.74 50.10 47.82 54.23 22.65 33.34 17.93 22.60 52.15 13.99 20.66
## [13] 17.46 52.29 35.02 57.25 23.88 37.59 51.33 45.17

ratings <- webpage %>%
  html_nodes(".star-rating") %>%
  html_attr("class") %>%
  sub("star-rating ", "", .)
ratings

##  [1] "Three" "One"   "One"   "Four"  "Five"  "One"   "Four"  "Three" "Four" 
## [10] "One"   "Two"   "Four"  "Five"  "Five"  "Five"  "Three" "One"   "One"  
## [19] "Two"   "Two"

books_df <- data.frame(
  Title = titles,
  Price = prices,
  Rating = ratings,
  stringsAsFactors = FALSE
)

Benzer şekilde websitelerindeki tabloları da çekebiliriz. Burada tabloyu sitede gördüğümüz şekliyle çektiğini görüyoruz.

url2 <- "https://en.wikipedia.org/wiki/List_of_largest_companies_in_the_United_States_by_revenue"
webpage2 <- read_html(url2)
tables <- webpage2 %>% html_nodes("table.wikitable")
companies_table <- tables[[1]] %>% html_table()

Xpath ile de veri çekebiliriz. xpath, html ve xml dökümanları ile çalışmak için güçlü bir tool. Bunu HTML lokasyonlarını öğrenmek için de her objenin kendi yerini öğrenmek için de kullanabiliriz.

url4 <- "http://quotes.toscrape.com/"
webpage4 <- read_html(url4)
quotes <- webpage4 %>%
  html_nodes(xpath = '/html/body/div[1]/div[2]/div[1]/div[1]/span[1]') %>%
  html_text() %>%
  gsub("^\\s+|\\s+$", "", .)

Bu şekilde manüel seçimler yapmak yeriine toplu bir şekilde de veri çekebiliriz. Sabah Gazetesi öğreneğine dönecek olursak:

sitemap_urls <- c("https://www.sabah.com.tr/sitemaparchives/post/2023-1.xml",
                  "https://www.sabah.com.tr/sitemaparchives/post/2023-2.xml",
                  "https://www.sabah.com.tr/sitemaparchives/post/2023-3.xml",
                  "https://www.sabah.com.tr/sitemaparchives/post/2022-12.xml")

2023 yılının ilk 3 ayında Sabah gazetesinde yayınlanmış haberleri çekelim.

df_sabah <- map_dfr(sitemap_urls, function(url) {
  sabah_data <- read_xml(url)
  sabah_urls <- xml_text(xml_find_all(sabah_data, "//d1:loc"))
  sabah_dates <- xml_text(xml_find_all(sabah_data, "//d1:lastmod"))
  sabah_links <- tibble(urls = sabah_urls, dates = sabah_dates)
  sabah_links_2 <- sabah_links %>%
    filter(str_detect(dates, "^2023-(01|02|03)")) %>%
    mutate(dates = substr(dates, 1, 10))
  return(sabah_links_2)
})
num_stories <- nrow(df_sabah)
sprintf("In the first 3 months of 2023, Sabah published %d stories.", num_stories)

## [1] "In the first 3 months of 2023, Sabah published 47228 stories."

Text Analysis

Sonra bu haberlerin kategorizasyonunu yapalım

df_sabah <- df_sabah %>%
  mutate(category = str_extract(urls, "(?<=sabah.com.tr\\/)[^\\/]+")) %>%
  mutate(category = case_when(
    category %in% c("dunya") ~ "World",
    category %in% c("egitim") ~ "Education",
    category %in% c("ekonomi", "finans") ~ "Economy",
    category %in% c("gundem") ~ "Politics",
    category %in% c("kultur-sanat", "magazin", "medya") ~ "Entertainment",
    category %in% c("spor") ~ "Sports",
    category %in% c("yazarlar") ~ "Opinion",
    TRUE ~ "Others"
  ))

df_sabah %>%
  count(category) %>%
  mutate(category = reorder(category, -n)) %>%
  ggplot(aes(x = category, y = n, fill = category)) +
  geom_bar(stat = "identity") +
  ggtitle("Number of News by Category in Sabah Newspaper 2023-01/03") +
  xlab("Category") +
  ylab("Count") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +  
  theme(legend.position = "none")

Haber başlıklarında kaç kere Kılıçdaroğlu kaç kere Erdoğan geçmiş diye bakalım

count_df <- df_sabah %>%
  mutate(names = str_extract_all(urls, "erdogan|kilicdaroglu")) %>%
  unnest(names) %>%
  count(names)
sabah_stories <- df_sabah %>%
  filter(str_detect(urls, "kilicdaroglu|erdogan")) %>%
  distinct(urls, .keep_all = TRUE)
invisible(any(duplicated(sabah_stories$urls)))
paste0("Erdogan appears ", count_df$n[1], 
       " times, and Kilicdaroglu appears ", count_df$n[2], " times.")

## [1] "Erdogan appears 1020 times, and Kilicdaroglu appears 314 times."

Haber içerikleri analizini yapalım

sabah_stories <- sabah_stories %>%
  filter(category == "Politics")
sabah_stories_combined <- data.frame(title = rep(NA, nrow(sabah_stories)),
                                     publish_date = rep(NA, nrow(sabah_stories)),
                                     main_text = rep(NA, nrow(sabah_stories)),
                                     erdogan = rep(NA, nrow(sabah_stories)),
                                     url = rep(NA, nrow(sabah_stories)),
                                     stringsAsFactors = FALSE)
sabah_stories_combined <- sabah_stories_combined %>%
  mutate(publish_date = sabah_stories$dates) %>%
  mutate(urls = sabah_stories$urls) %>%
  mutate(erdogan = ifelse(grepl('erdogan', urls), 1, 0)) %>%
  select(title, publish_date, main_text, erdogan, urls)

Hem haber başlığını hem de textini çıkaralım

for (i in seq_along(sabah_stories$urls)) {
  # get HTML content
  html <- read_html(sabah_stories$urls[i])
  
  # Extract title
  title <- html %>% html_nodes("h1.pageTitle") %>% html_text()
  
  # Extract main text
  main_text_p <- html %>%
    html_nodes("p") %>%
    html_text()
  
  if (length(main_text_p) > 0) {  
    main_text <- main_text_p
  } else { 
    main_text_newsBox <- html %>%
      html_nodes(".newsBox.selectionShareable") %>%
      html_text()
    main_text <- main_text_newsBox
  }
  
  
  sabah_stories_combined$title[i] <- title
  sabah_stories_combined$main_text[i] <- paste(main_text, collapse = "\n")
  
  
  Sys.sleep(7)
}

Kaydedelim ve tekrar çağıralım

write_csv(sabah_stories_combined, "sabahstories.csv")

sabah_stories_combined <- read_csv("sabahstories.csv")
sabah_stories_combined$publish_date <- as.Date(sabah_stories_combined$publish_date)

Günlük olarak çıkan haberlere bakalım

counts <- sabah_stories_combined %>%
  mutate(date = as.Date(publish_date, format = "%Y-%m-%d")) %>%
  group_by(date) %>%
  summarise(erdogan_count = sum(grepl("erdogan", urls, ignore.case = TRUE)),
            kilicdaroglu_count = sum(grepl("kilicdaroglu", urls, ignore.case = TRUE))) %>%
  ungroup()

Başlıklardan noktalamaları çıkaralım, küçük harfe çevirelim ve içinde kılıcdaroğlu geçenleri ayrı bir değişkene alalım.

sabah_stories_combined$title2 <- 
  str_to_lower(str_replace_all(sabah_stories_combined$title, 
                               "[[:punct:]]", ""))
sabah_stories_combined <- sabah_stories_combined %>% 
  mutate(kilicdaroglu = ifelse(grepl('kilicdaroglu', urls), 1, 0))

Stopwords ekleyelim. Bunu R’ın içinde Türkçe stopwords olmadığı için ayrı paketle yapıyoruz.

tr_stopwords <- stopwords::stopwords("tr", source = "stopwords-iso")
tr_stopwords <- tr_stopwords[tr_stopwords != "iyi"]

Sentiment Analysis

Erdoğan geçen başlıklardaki kelimeleri sayalım

erdogan_words <- sabah_stories_combined %>%
  filter(erdogan == 1) %>%
  select(title2) %>%
  unnest_tokens(word, title2) %>%
  filter(!word %in% c("son", "dakika", 
                      tr_stopwords, "recep", "tayyip", "daki̇ka")
         & !grepl("erdoğan", word)) %>%
  count(word, sort = TRUE)

Aynısını Kılıçdaroğlu için yapalım

kilicdaroglu_words <- sabah_stories_combined %>%
  filter(kilicdaroglu == 1) %>%
  select(title2) %>%
  unnest_tokens(word, title2) %>%
  filter(!word %in% c("son", "dakika", 
                      tr_stopwords, "kemal", "daki̇ka") &
           !grepl("kılıçdaroğlu", word)) %>%
  count(word, sort = TRUE)

İkisi için de en çok kullanılan 15 kelimeyi çıkaralım.

top_erdogan_words <- erdogan_words %>%
  slice_head(n = 15)
top_kilicdaroglu_words <- kilicdaroglu_words %>%
  slice_head(n = 15)

erdogan_plot <- ggplot(top_erdogan_words, 
                       aes(x = reorder(word, n), y = n)) +
  geom_bar(stat = "identity", fill = "orange") +
  labs(title = "Most Common Words in Erdogan Titles", 
       y = "Word", x = "Count") +
  theme_minimal() +
  theme(axis.text.y = element_text(angle = 0, hjust = 1)) +
  coord_flip()
ggsave("erdogan.png", dpi = 300)

kilicdaroglu_plot <- ggplot(top_kilicdaroglu_words, 
                            aes(x = reorder(word, n), y = n)) +
  geom_bar(stat = "identity", fill = "red") +
  labs(title = "Most Common Words in Kilicdaroglu Titles", 
       y = "Word", x = "Count") +
  theme_minimal() +
  theme(axis.text.y = element_text(angle = 0, hjust = 1)) +
  coord_flip()
ggsave("kilicdaroglu.png", dpi = 300)

Bunu world cloud olarak da yapabiliriz.

kilicdaroglu_wordcloud <- ggplot(top_kilicdaroglu_words, 
                                 aes(label = word, size = n)) +
  geom_text_wordcloud() +
  scale_size_area(max_size = 15) +
  theme_void() +
  labs(title = "Most Common Words in Kilicdaroglu Titles")
kilicdaroglu_wordcloud

erdogan_wordcloud <- ggplot(top_erdogan_words, 
                                 aes(label = word, size = n)) +
  geom_text_wordcloud() +
  scale_size_area(max_size = 30) +
  theme_void() +
  labs(title = "Most Common Words in Erdogan Titles")
erdogan_wordcloud

Text analizi kısmında textin kendisine ve içinde kelime olarak geçtiğine odaklanmıştık ancak içerik olarak nelerden bahsettiklerine hiç bakmadık.

sentiment <- read_xlsx("sentiment.xlsx")

## New names:
## • `` -> `...1`

Önce verimizi alıp analiz için hazırlayalım

dfm <- dfm(corpus(sentiment$fulltext) %>% 
             tokens(remove_punct = TRUE) %>%
             tokens_tolower() %>%
             tokens_remove(tr_stopwords) %>%
             tokens_wordstem(language = "tr")) 

Haberleri konulara bölelim ve beraber bulundukları haberlere göre ortlama olarak nasıl ayrıldıklarına bakalım.

set.seed(94)
lda_model <- textmodel_lda(dfm, k = 10)
top_terms <- terms(lda_model, n= 10)
top_terms_df <- data.frame(top_terms)
top_terms_df_1 <- top_terms_df[ , 1:5]
top_terms_df_2 <- top_terms_df[ , 6:10]
gt_table_1 <- gt::gt(top_terms_df_1) %>% 
  gt::tab_header(title = "Top Terms per Topic (1-5)")
gt_table_2 <- gt::gt(top_terms_df_2) %>% 
  gt::tab_header(title = "Top Terms per Topic (6-10)")

gt_table_1

Top Terms per Topic (1-5)
topic1	topic2	topic3	topic4	topic5
kılıçdaroğlu	erdoğa	em	kemal	erdoğa
genel	millet	erdoğa	bay	yıl
parti	seç	çocuk	millet	ülke
kemal	genç	dedi	işçi	millet
başka	başka	ev	iş	büyük
chp	ak	hesap	başka	son
kılıçdaroğlu'n	ülke	medya	ifade	başka
aday	söz	atık	erdoğa	ala
hdp	an	eş	ada	yer
i̇yi̇	parti	erdoğan'	diyecek	dünya

gt_table_2

Top Terms per Topic (6-10)
topic6	topic7	topic8	topic9	topic10
başka	ülke	depre	başka	su
erdoğa	türki	erdoğa	cumhurbaşkan	açılış
recep	türk	vatandaş	terör	lira
tayyip	dünya	başka	medya	yol
cumhurbaşka	uluslararas	çalışma	erdoğan'	başka
genel	konu	bölge	i̇sveç'	yüz
ziyaret	el	konut	karar	proje
erdoğan'	barış	felaket	örgüt	tünel
görüşme	halk	afet	eyle	i̇stanbul
kabul	ilişki	yer	i̇sveç	yatır

Bu sefer 5 konu altınada toplarlayıp farkına bakalım. Ne değişti?

set.seed(94)
lda_model2 <- textmodel_lda(dfm, k = 5)
top_terms2 <- terms(lda_model2, n = 10)
top_terms2_df <- data.frame(top_terms2)
gt_table_3 <- gt::gt(top_terms2_df) %>% 
  gt::tab_header(title = "Top Terms per Topic")
gt_table_3

Top Terms per Topic
topic1	topic2	topic3	topic4	topic5
başka	erdoğa	erdoğa	başka	millet
parti	depre	yıl	erdoğa	erdoğa
genel	başka	başka	cumhurbaşka	ülke
kılıçdaroğlu	vatandaş	ülke	tayyip	genç
kemal	bölge	dünya	recep	söz
chp	çalışma	su	türk	siz
aday	konut	açılış	erdoğan'	an
kılıçdaroğlu'n	felaket	erdoğan'	görüşme	zama
cumhurbaşka	devlet	proje	türki	türki
6	afet	tarih	kabul	dedi

topics <- topics(lda_model2)
sentiment$Topic <- topics
topic_labels <-  c("Erdogan","Turkey", "Kilicdaroglu", "Earthquake", "Election" )
sentiment$Topic <- factor(sentiment$Topic, labels = topic_labels)
topic_counts <- table(sentiment$Topic)
topic_counts <- as.data.frame(topic_counts)

Top Terms per Topic
topic1	topic2	topic3	topic4	topic5
an	hdp	konut	emine	hizmet
terör	masa	felaket	türk	i̇stanbul
i̇sveç	akşener	depremzede	görüşme	dünya
14	i̇yi̇	kahramanmaraş	dünya	alan
bay	adaylık	ev	kabul	proje
genç	lı	afet	halk	su
mayıs	partisi	hatay	uluslararası	kadın
siyasi	koalisyon	il	mesaj	sahip
masa	isim	yardım	ilet	hayat
örgüt	görüşme	bura	atık	yatırım

Konu başına haç haber düştüğünü bu şekilde de görebiliriz

Cluster Frequencies

Cluster	Frequency
1	7
2	35
3	6
4	64
5	480

Kelimeleri pozitif ne negatif olarak bölelim. Bunu şimdi öğrenmeyeceğimiz lasso regresyonu diye bir yöntemle yapıyoruz ve böylece hangi kelimelerin neyi tahmin ettiğini görebiliyoruz. Bu kelimeler negatif veya pozitif olarak grupladığımız haberlerde en çok geçen kelimeler.

predictors <- convert(dfm, to = 'matrix')
response <- sentiment$erdogan
set.seed(94)
lasso <- cv.glmnet(predictors, response, family = "binomial", alpha = 1)
tmp_coeffs <- coef(lasso, s = "lambda.min")
coeffs_df <- data.frame(name = tmp_coeffs@Dimnames[[1]][tmp_coeffs@i + 1], coefficient = tmp_coeffs@x)
coeffs_df <- coeffs_df[coeffs_df$name != "(Intercept)",]
coeffs_df <- coeffs_df[order(-abs(coeffs_df$coefficient)),]
top_positive_words <- head(coeffs_df[coeffs_df$coefficient > 0,], 10)
top_negative_words <- head(coeffs_df[coeffs_df$coefficient < 0,], 10)

Negatif haberlerde en çok geçen kelimeler

Top Positive Words
Word	Coefficient
ereğli'	7.2696803
beledi̇yeni̇n	2.5043161
başvura	2.4157332
tazmina	2.3183030
toplanmış	1.7814972
daires	1.0560055
yama	0.9433218
tayyip	0.6628557
etkileşim	0.6190165
onardi	0.4788912

Pozitif haberlerde en çok geçen kelimeler

top_negative_words <- gt::gt(top_negative_words) %>%
  gt::tab_header(title = "Top Negative Words") %>%
  gt::cols_label(name = "Word", coefficient = "Coefficient")
top_negative_words

Top Negative Words
Word	Coefficient
kılıçdaroğlu’n	-4.749593
sevigen'	-4.154531
ağırlayarak	-3.702199
flört	-3.103029
pervasız	-2.786683
kılıçdaroğlu’na	-2.731575
tüm	-2.563847
mensubiyet	-2.469966
çürüdü	-2.436991
açiklama	-2.348500

İsimleri çıkarıp bir daha bakalım

coeffs_df <- coeffs_df[-grep('erdoğa|kılıçdaroğlu|tayyip|recep|kemal', coeffs_df$name), ]
coeffs_df <- coeffs_df[order(-abs(coeffs_df$coefficient)),]
top_words <- head(coeffs_df, 10)

Bunu yaparken dikkat edilmesi gereken noktayı hatırlatmak gerekirse: robots.txt dışına çıkmak etik olarak görülmez ve zaten bir çok sitede buna yönelik engeller vardır. Sürekli veri çekmeye çalışırsanız da sizi IP üstünden geçici veya kalıcı olarak engelleyebilir. Bu da ileride herhangi bir çalışma yapmamızın önüne geçer.

Bazı sitelerse doğrudan scrape edilmeyi istemezler ancak application programming interface (API) aracılığıyla buna izin verirler. Bunlar ücretli veya ücretsiz olabilir ancak düzenli bir veri çekimine olanak tanırlar.

Daha advanced sentiment analiz icin R’in paketleri var ancak bunlar Ingilizce ile calisiyor. Turkce yapmak icinse Python kullanmamiz lazim. Neyse ki R kullanirken Python cagirip isimizi yapitrabiliyoruz.

Once Python islemlerini yapalim

from transformers import AutoModelForSequenceClassification, AutoTokenizer, pipeline
import pandas as pd
import warnings
from openpyxl import Workbook
warnings.filterwarnings('ignore', category=FutureWarning)
warnings.filterwarnings('ignore', category=UserWarning)
model = AutoModelForSequenceClassification.from_pretrained("savasy/bert-base-turkish-sentiment-cased")
tokenizer = AutoTokenizer.from_pretrained("savasy/bert-base-turkish-sentiment-cased")
sa = pipeline("sentiment-analysis", tokenizer=tokenizer, model=model, device=-1)
input_file = "/Users/alionurgitmez/Desktop/R Dersi/Hafta 15/sabah_PS2.xlsx"
output_file = "after_analysis.xlsx"
df = pd.read_excel(input_file)
print(df.head())

##    Unnamed: 0  ...                                lemmatized_fulltext
## 0           0  ...  Başkan Erdoğan dan Ali İhsan Destici taziye me...
## 1           1  ...  Emine Erdoğan vatandaş böyle seslen Duyarlı he...
## 2           2  ...  Başkan Erdoğan Yeni Azerbaycan Partisi Gençler...
## 3           3  ...  Son dakika Denizli emekçi kadın buluş Başkan E...
## 4           4  ...  O fotoğraf kahraman konuş Huriye teyze Erdoğan...
## 
## [5 rows x 8 columns]

Burada sentiment paketini calistiyoruz

sentiment_results = df['title'].apply(lambda x: sa(x)[0])
df['sentiment_label'] = [result['label'] for result in sentiment_results]
df['sentiment_score'] = [result['score'] for result in sentiment_results]
df.to_excel(output_file, index=False)

Burada kaydettigimiz veriyi R’a aliyoruz. Boylece R’a geri donmus olduk. Ikisi de ayni anda calisabiliyor.

## [1] "Average positive sentiment score of Kilicdaroglu articles is: 0.85"
## [1] "Average positive sentiment score of Erdogan articles is: 0.91"

Gorsellestirmesini de yapalim

sabah_stories_positive$date <- as.Date(sabah_stories_positive$date)

average_sentiment2 <- sabah_stories_positive %>%
  mutate(week = as.integer(format(date, "%U"))) %>%
  group_by(erdogan, week) %>%
  summarize(average_sentiment = mean(sentiment_score))

ggplot(average_sentiment2, aes(x = week, y = average_sentiment, color = factor(erdogan))) +
  geom_line() +
  labs(x = "Week", y = "Average Sentiment", color = "Leader") +
  scale_color_manual(values = c("red", "orange"), labels = c("Kilicdaroglu", "Erdogan")) +
  scale_x_continuous(breaks = unique(average_sentiment2$week)) +
  theme_minimal()

Spatial Analysis

Coğrafi analiz özellikle IR alanında ve conflict çalışmalarında gelişmekte olan bir analiz yöntemi. Siyasi ve fiziksel harita veya uydu bilgilerini kullanarak çıkarımlar yapmamıza olanak sağlıyor. Text analizi gibi çok geniş bir alan olduğu için kısa bir giriş yapacağım. Öncelikle bu alanda zorunlu olan paketleri yükleyelim.

library(sf)
library(raster)
library(spData)
library(viridis)
library(scales)

Örnek bir veriyi hazırlayıp bakalım. Bu veri North Carolina sudden infant death syndrome (SIDS) oranlarına bakıyor. Aynı zamanda bu veriyi görselleştirelim. st olarak belirttiğimiz için buna göre görselleştirme yaptı.

nc <- st_read(system.file("shapes/sids.shp", package = "spData"))

## Reading layer `sids' from data source 
##   `/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/spData/shapes/sids.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 22 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## CRS:           NA

plot(nc["SID74"]) 

Her bölge çevresine 10 km bir güvenli alan oluşturalım.

buffered_nc <- st_buffer(nc, dist = 10000)
plot(buffered_nc, max.plot = 22)

Kesişen alanları bulalım. Bu ilk verimiz ve 10 km koyduğumuz veri arasında kesişen alanları gösteriyor.

intersection_nc <- st_intersection(nc, buffered_nc)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

plot(intersection_nc)

## Warning: plotting the first 10 out of 44 attributes; use max.plot = 44 to plot
## all

Farklı bir görselleştirme alternatifi olarak da bunu kullanabiliriz.

library(tmap)
tm_shape(nc) +
  tm_polygons(col = "SID74", palette = "Blues", title = "SIDS Rate 1974") +
  tm_layout(legend.outside = TRUE)

Haritalar ile çalışmanın farklı bir örneğini de Türkiye seçimleri üstünden görebiliriz. Bunu ilçe bazlı ve il bazlı seçim sonuçları üstünden yapabiliriz.

library(haven)
library(tidyverse)
library(TRmaps)
ilce_election <- read_dta("Elections_wide.dta")
ilce_election <- ilce_election %>% dplyr::select(adm1_tr, adm2_tr, akp_2018P, chp_2018P, mhp_2018P, 
                                          iyi_2018P, hdp_2018P, sp_2018P)

ilce_election$adm2_tr <- gsub(".*MERKEZ*", "MERKEZ", ilce_election$adm2_tr)

convert_to_lower <- function(str) {
  converted_str <- str_to_sentence(str, locale = "tr")  
  return(converted_str)
}
ilce_election$adm1_tr <- convert_to_lower(ilce_election$adm1_tr)
ilce_election$adm2_tr <- convert_to_lower(ilce_election$adm2_tr)

ilce_election$il_ilce <- paste(ilce_election$adm1_tr, ilce_election$adm2_tr, sep = "-")

ilce_election <- ilce_election %>% dplyr::select(-adm1_tr, -adm2_tr)

tr_ilce <- as.data.frame(st_as_sf(tr_ilce))

tr_ilce <- tr_ilce %>% dplyr::select(il_ilce, geometry, tuik_no, Shape_Area, Shape_Leng)

tr_ilce <- tr_ilce[tr_ilce$il_ilce != "Hakkari-Derecik", ]

ilce_election <- ilce_election %>%
  mutate(il_ilce = ifelse(il_ilce == "Ankara-Kazan", "Ankara-Kahramankazan", il_ilce))
ilce_election <- ilce_election %>%
  mutate(il_ilce = ifelse(il_ilce == "Samsun-19 mayıs", "Samsun-19 Mayıs", il_ilce))
ilce_election <- ilce_election %>%
  mutate(il_ilce = ifelse(il_ilce == "Kırıkkale-Bahşili", "Kırıkkale-Bahşılı", il_ilce))

final_ilce_election <- merge(ilce_election, tr_ilce, by.x = "il_ilce", by.y = "il_ilce")

final_ilce_election$cumhur <- final_ilce_election$akp_2018P + final_ilce_election$mhp_2018P
final_ilce_election$millet <- final_ilce_election$chp_2018P + final_ilce_election$iyi_2018P + final_ilce_election$sp_2018P
final_ilce_election$hdp <- final_ilce_election$hdp_2018P 


final_ilce_election$kazanan <- ifelse(final_ilce_election$hdp > final_ilce_election$millet & final_ilce_election$hdp > final_ilce_election$cumhur, "HDP",
                                      ifelse(final_ilce_election$millet > final_ilce_election$hdp & final_ilce_election$millet > final_ilce_election$cumhur, "MILLET",
                                             ifelse(final_ilce_election$cumhur > final_ilce_election$hdp & final_ilce_election$cumhur > final_ilce_election$millet, "CUMHUR", "ESIT")))

final_ilce_election_sf <- st_as_sf(final_ilce_election)
my_colors <- c("HDP" = "#8f1ca5", "MILLET" = "red", "CUMHUR" = "orange", "CHP" = "firebrick1", 
               "MHP" = "darkred", "IYI" = "cyan", "SAADET" = "maroon", "ESIT" = "darkgrey")

ggplot(final_ilce_election_sf) + 
  geom_sf(aes(fill = kazanan)) +
  scale_fill_manual(values = my_colors) +
  labs(fill = "Kazanan") + 
  theme_void() +
  ggtitle("Tum Ittifaklar Beraber") 

Spatial data ile bir analiz öğreneğini inceleyelim. Öncelikle kullanılacak verileri alıp dönüşüm yapıyoruz.

url <- "https://www2.census.gov/geo/tiger/TIGER2019/COUNTY/tl_2019_us_county.zip"
download.file(url, destfile = "us_counties.zip")
unzip("us_counties.zip")

counties <- st_read('/Users/alionurgitmez/Desktop/R Dersi/Hafta 15/tl_2019_us_county.shp')

## Reading layer `tl_2019_us_county' from data source 
##   `/Users/alionurgitmez/Desktop/R Dersi/Hafta 15/tl_2019_us_county.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3233 features and 17 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.43979
## Geodetic CRS:  NAD83

counties_transformed <- st_transform(counties, 2163)

covid_url <- "https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-counties.csv"
covid_data <- read.csv(covid_url)
write_csv(covid_data, "covid_data.csv")

covid_data <- read_csv("covid_data.csv")

## Rows: 2502832 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (2): county, state
## dbl  (3): fips, cases, deaths
## date (1): date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

covid_latest <- covid_data %>%
  group_by(fips) %>%
  filter(date == max(date)) %>%
  ungroup() %>%
  dplyr::select(fips, cases, deaths)
covid_latest$fips <- as.character(covid_latest$fips)

counties_covid <- counties_transformed %>%
  left_join(covid_latest, by = c("GEOID" = "fips"))
counties_covid <- counties_covid %>%
  mutate(cases_per_100k = (cases / as.numeric(ALAND)) * 100000)

Vakaların çok yüksek olduğu yerlere bakalım

hotspots <- counties_covid %>%
  arrange(desc(cases_per_100k)) %>%
  slice_head(n = 10)
print(hotspots[, c("NAME", "STATEFP", "cases", "cases_per_100k")])

## Simple feature collection with 10 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1930000 ymin: -428834.5 xmax: 2352825 ymax: 140953.9
## Projected CRS: NAD27 / US National Atlas Equal Area
##                    NAME STATEFP  cases cases_per_100k
## 1               Suffolk      25 226017      149.81600
## 2                Hudson      34 179135      149.72731
## 3          Philadelphia      42 317808       91.38143
## 4  District of Columbia      11 143943       90.90732
## 5            Alexandria      51  32647       84.40414
## 6             Arlington      51  46198       68.61124
## 7                 Essex      34 222558       68.15279
## 8         Manassas Park      51   3828       58.54358
## 9                Nassau      36 425386       57.70325
## 10                Union      34 151057       56.75346
##                          geometry
## 1  MULTIPOLYGON (((2322932 130...
## 2  MULTIPOLYGON (((2144411 -12...
## 3  MULTIPOLYGON (((2074643 -24...
## 4  MULTIPOLYGON (((1955476 -40...
## 5  MULTIPOLYGON (((1956985 -41...
## 6  MULTIPOLYGON (((1952303 -40...
## 7  MULTIPOLYGON (((2122970 -12...
## 8  MULTIPOLYGON (((1930000 -42...
## 9  MULTIPOLYGON (((2171046 -10...
## 10 MULTIPOLYGON (((2119662 -14...

Yakındaki yerler üstünden vakaların yüksek olduğu yerleri belirleyelim

high_case_counties <- counties_covid %>%
  filter(cases_per_100k > quantile(cases_per_100k, 0.9, na.rm = TRUE))

neighbors <- st_join(counties_covid, high_case_counties, join = st_touches, left = FALSE)

clusters <- neighbors %>%
  group_by(GEOID.x) %>%
  summarise(cluster_size = n())

Şimdiyse elimizdeki verileri görselleştirelim:

library(viridis)
library(scales)
ggplot(data = counties_covid) +
  geom_sf(aes(fill = cases_per_100k), color = NA) +
  scale_fill_viridis(trans = "log", labels = comma, 
                     name = "Cases per 100,000\n(log scale)") +
  theme_minimal() +
  theme(legend.position = "right") +
  labs(title = "COVID-19 Cases per 100,000 Population by US County",
       subtitle = "Data as of latest available date",
       caption = "Source: New York Times COVID-19 Data")

Bu tip spatial veriye vektör verisi diyorduk. Başka bir tip spatial verisi ise raster dediğimiz ve fiziki haritalarda da görebileceğimiz harita tipi.

library(raster)
library(terra)
elevation <- raster("raster.tif")
print(elevation)

## class      : RasterLayer 
## dimensions : 10812, 10812, 116899344  (nrow, ncol, ncell)
## resolution : 9.259259e-05, 9.259259e-05  (x, y)
## extent     : -106.0006, -104.9994, 38.99944, 40.00056  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=NAD83 +no_defs 
## source     : raster.tif 
## names      : raster

summary(elevation)

##           raster
## Min.    1555.488
## 1st Qu. 2363.127
## Median  2744.982
## 3rd Qu. 3063.692
## Max.    4317.352
## NA's       0.000

plot(elevation, main = "Elevation in Colorado", xlab = "Longitude", ylab = "Latitude")

Buradan fiziki bir analiz de yapmamız mümkün. Mesela birazdan bahseceğim verideki bigi yükseklik yerine eğim kullanarak analiz yapabiliriz.

slope <- terrain(elevation, opt = "slope", unit = "degrees")
aspect <- terrain(elevation, opt = "aspect", unit = "degrees")
par(mfrow = c(1, 2))
plot(slope, main = "Slope in Colorado", xlab = "Longitude", ylab = "Latitude")
plot(aspect, main = "Aspect in Colorado", xlab = "Longitude", ylab = "Latitude")

Terror example

library(rnaturalearth)
library(maps)
library(rnaturalearthdata)
library(TRmaps)
library(Hmisc)
library(gt)
library(tidyverse)
load("pkkattacks.rdata")

pkkattacks <- dataNewC %>% dplyr::select(attack_p22, nkill_n_p22, nwound_n_p22, margin,
                                   unemployment, literacy, infant_mort_perthousand, attack_p1,
                                   subsidy, curfew, border, district, X_coordinate, Y_coordinate, population, district)
pkkattacks$curfew[pkkattacks$curfew == "TRUE"] <- 1
pkkattacks$border[pkkattacks$border == "TRUE"] <- 1

Yükseklik verisini çekelim. Uzun süreceği için ben öncede çekmiştim. Onu kulalanacağım.

elevation_url <- "https://maps.googleapis.com/maps/api/elevation/json"
geocoding_url <- "https://maps.googleapis.com/maps/api/geocode/json"
api_key <- "AIzaSyDxC-P1tX7vAZ6JVnevvflwZ0YmH6sQn4"
elevation_data <- c()
location_names <- c()
for (i in 1:nrow(pkkattacks)) {
  
  
  lat <- dataNewC[i, "Y_coordinate"]
  lng <- dataNewC[i, "X_coordinate"]
  
  
  elevation_params <- list(
    locations = paste(lat, lng, sep = ","),
    key = api_key
  )
  
  
  elevation_response <- GET(elevation_url, query = elevation_params)
  
  
  elevation <- content(elevation_response)$results[[1]]$elevation
  
  
  geocoding_params <- list(
    latlng = paste(lat, lng, sep = ","),
    key = api_key
  )
  
  
  geocoding_response <- GET(geocoding_url, query = geocoding_params)
  
  
  location_name <- content(geocoding_response)$results[[1]]$formatted_address
  
  
  elevation_data <- c(elevation_data, elevation)
  location_names <- c(location_names, location_name)
}

library(readxl)
pkkattacks_yeni <- read_xlsx("pkkattacks.xlsx")

## New names:
## • `district_id` -> `district_id...2`
## • `province.x` -> `province.x...3`
## • `province.x` -> `province.x...11`
## • `province.y` -> `province.y...44`
## • `province.y` -> `province.y...48`
## • `district_id` -> `district_id...56`
## • `province.x2` -> `province.x2...57`
## • `province.x2` -> `province.x2...66`
## • `province.y2` -> `province.y2...99`
## • `province.y2` -> `province.y2...103`

pkkattacks_yeni <- pkkattacks_yeni %>% dplyr::select(attack_p22, nkill_n_p22, nwound_n_p22, margin, unemployment, literacy, infant_mort_perthousand, attack_p1, subsidy, curfew, border, district, X_coordinate, Y_coordinate, population, district, elevation, location_name)

Sınıra olan uzaklıkları hesaplayalım

countries <- ne_countries(returnclass = "sf")
syria_border <- subset(countries, name == "Syria")
iraq_border <- subset(countries, name == "Iraq")
iran_border <- subset(countries, name == "Iran")
cities_sf <- st_as_sf(pkkattacks_yeni, coords = c("X_coordinate", "Y_coordinate"), crs = 4326)
pkkattacks_yeni$syria <- st_distance(cities_sf, syria_border)
pkkattacks_yeni$iran <- st_distance(cities_sf, iran_border)
pkkattacks_yeni$iraq <- st_distance(cities_sf, iraq_border)
pkkattacks_yeni <- pkkattacks_yeni %>%
  mutate(syria = as.numeric(str_remove(syria, "\\[m\\]")) / 1000)
pkkattacks_yeni <- pkkattacks_yeni %>%
  mutate(iran = as.numeric(str_remove(iran, "\\[m\\]")) / 1000)
pkkattacks_yeni <- pkkattacks_yeni %>%
  mutate(iraq = as.numeric(str_remove(iraq, "\\[m\\]")) / 1000)
pkkattacks_yeni <- pkkattacks_yeni %>%
  mutate(border_distance = pmin(syria, iraq, iran, na.rm = TRUE))

Harita yapalım. İlk olarak yükseklik haritası olsun.

map <- pkkattacks_yeni %>% dplyr::select(district, X_coordinate, Y_coordinate, elevation, population, attack_p22)

pkkattacks_sf <- st_as_sf(map, coords = c("X_coordinate", "Y_coordinate"), crs = 4326)

turkey_map <- ne_states(country = "Turkey", returnclass = "sf")

color_scale <- terrain.colors(n = 7)
color_scale[7] <- rgb(102, 51, 0, maxColorValue = 255)
districtmap <- ggplot() +
  geom_sf(data = turkey_map, fill = "transparent", color = "black") +
  geom_sf(data = pkkattacks_sf, aes(color = elevation, size = population)) +
  scale_color_gradientn(colors = color_scale, name = "Elevation") +
  labs(title = "Map of Districts", color = "Elevation", size = "Population") +
  guides(size = "none") +
  theme_bw()
districtmap

Saldırıların haritası

attackmap <- ggplot() +
  geom_sf(data = turkey_map, fill = "transparent", color = "black") +
  geom_sf(data = pkkattacks_sf[pkkattacks_sf$attack_p22 != 0,], aes(color = elevation, size = attack_p22)) +
  scale_color_gradientn(colors = color_scale, name = "Elevation") +
  labs(title = "Map of Attacks", color = "Elevation", size = "Attacks") +
  guides(size = FALSE) +
  theme_bw()

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

attackmap

Daha detaylı bir harita yapalım

tr_ilce_sf <- st_as_sf(tr_ilce)
attackmap2 <- ggplot() + 
  geom_sf(data = tr_ilce_sf) +
  geom_sf(data = turkey_map, fill = "transparent", color = "black") +
  geom_sf(data = pkkattacks_sf[pkkattacks_sf$attack_p22 != 0,], aes(color = elevation, size = attack_p22)) +
  scale_color_gradientn(colors = color_scale, name = "Elevation") +
  labs(title = "Map of Attacks", color = "Elevation", size = "Attacks") +
  guides(size = "none") +
  theme_bw()
attackmap2

Network Analysis

Burada bahsettiğim diğer alanlar gibi Network de çok hızlıca gelişen bir alan. Ancak çok fazla bilimsel derinliği bir alan olduğu için çok detaylıca bahsetmem mümkün değil. Ondan dolayı temel analiz kısmına odaklanıp, görselleştirme ve basit kavramları anlatacağım. Zaten başlıbaşına programı olan bir alan.

library(igraph)
library(igraphdata)
data(karate)

Network alanının çok temel kavramları var:

• Nodes(Veritces): Network’ün unitleridir.
• Edges(Links): Nodelar arasındaki ilişkiyi gösterir.

Farklı network tipleri: - Directed vs undirected: İki node arasındaki edge’in yani ilişkinin bir yönü olup olmadığını gösterir. - Weighted vs unweighted: İki node arasındaki ilişkinin gücünü gösterir.

Basit network

g <- graph(edges = c(1, 2, 1, 3, 2, 4, 3, 4, 4, 5), directed = FALSE)
plot(g, vertex.label = V(g)$name, edge.arrow.size = 0.5)

Direted network

g_directed <- graph(edges = c(1, 2, 1, 3, 2, 4, 3, 4, 4, 5), directed = TRUE)
plot(g_directed, vertex.label = V(g_directed)$name, edge.arrow.size = 0.5)

Weighted network

g_weighted <- graph(edges = c(1, 2, 1, 3, 2, 4, 3, 4, 4, 5), directed = FALSE)
E(g_weighted)$weight <- c(1, 2, 1.5, 2.5, 3)
plot(g_weighted, vertex.label = V(g_weighted)$name, edge.width = E(g_weighted)$weight)

Sized network

g_sized <- graph(edges = c(1, 2, 1, 3, 2, 4, 3, 4, 4, 5), directed = FALSE)
vertex_degree <- degree(g_sized)
plot(g_sized, vertex.label = V(g_sized)$name, vertex.size = vertex_degree * 10, edge.arrow.size = 0.5)

Adjacency Matrix ile nodeların diğer nodlara olan ilişkisine bakabiliriz.

as_adjacency_matrix(g)

## 5 x 5 sparse Matrix of class "dgCMatrix"
##               
## [1,] . 1 1 . .
## [2,] 1 . . 1 .
## [3,] 1 . . 1 .
## [4,] . 1 1 . 1
## [5,] . . . 1 .

Dataframe ile de network grafik oluşturabiliriz

edges <- data.frame(from = c("A", "A", "B", "C", "D"),
                    to = c("B", "C", "D", "E", "E"))
g_df <- graph_from_data_frame(edges, directed = TRUE)
plot(g_df)

plot(g, vertex.label = V(karate)$name, edge.arrow.size = 0.5)

## This graph was created by an old(er) igraph version.
##   Call upgrade_graph() on it to use with the current igraph version
##   For now we convert it on the fly...

Önemli bazı metrikleri çıkartma. Centrality her node’un diğer kaç node’a bağlı olduğunu gösterir. Böylece bir node’un ne kadar önemli olup olmadığını görebiliriz.

degree(g)

## [1] 2 2 2 3 1

Betweenness ise bir node’un diğer iki nod arasında kaç defa en kısa yolu sağladığını gösterir.

betweenness(g)

## [1] 0.5 1.0 1.0 3.5 0.0

Closeness ise bir node’un diğer nodelar ile ne kadar yakın olduğuna bakar. Böylece diğer nodelar ile ne kadar interact ettiğini görebiliriz. Ortalama bir değer olarak alınır. Yüksek değer daha merkezi anlamına gelir.

closeness(g)

## [1] 0.1428571 0.1666667 0.1666667 0.2000000 0.1250000

Diameter ise iki netwoek arasındaki maksiumum uzaklığı hesaplar. Kaç edge geçerek bu değere ulaşacağımızı gösterir.

diameter(g)

## [1] 3

Ortalama mesafe ise nodelar arasındaki ortalama edge sayısına bakar. Yine merkezi bir node olduğuna bakmamız açısından önemlidir.

mean_distance(g, directed = FALSE)

## [1] 1.6

Daha gelişmiş bir grafik oluşturması için:

library(tidygraph)
library(ggraph)
V(g)$name <- c("A", "B", "C", "D", "E")
g_tbl <- as_tbl_graph(g)
ggraph(g_tbl, layout = 'fr') +
  geom_edge_link(aes(edge_alpha = 0.5)) +
  geom_node_point(aes(size = degree(g_tbl)), color = 'skyblue') +
  geom_node_text(aes(label = name), repel = TRUE) +
  theme_void()

Cliqueler nodeların bir alt grubudur. Nodeların diğer nodelara doğrudan bağlı olduğu gruplara denir.

cliques_found <- cliques(g)
print(cliques_found)

## [[1]]
## + 1/5 vertex, named, from f88f7e7:
## [1] D
## 
## [[2]]
## + 1/5 vertex, named, from f88f7e7:
## [1] A
## 
## [[3]]
## + 1/5 vertex, named, from f88f7e7:
## [1] E
## 
## [[4]]
## + 2/5 vertices, named, from f88f7e7:
## [1] D E
## 
## [[5]]
## + 1/5 vertex, named, from f88f7e7:
## [1] C
## 
## [[6]]
## + 2/5 vertices, named, from f88f7e7:
## [1] A C
## 
## [[7]]
## + 2/5 vertices, named, from f88f7e7:
## [1] C D
## 
## [[8]]
## + 1/5 vertex, named, from f88f7e7:
## [1] B
## 
## [[9]]
## + 2/5 vertices, named, from f88f7e7:
## [1] A B
## 
## [[10]]
## + 2/5 vertices, named, from f88f7e7:
## [1] B D

Louvalin algoritması kullanarak community yani birbiri ile yakından ilişkili olan grupları görebiliriz. Bunu özellikle ilişkiler içinde grup bulmak istiyorsak kullanırız. Cliquelerin akine kendi içinde güçlü bağlantıya ama aralarında zayıf bağlantıya sahip olamalarıdır.

community <- cluster_louvain(g)
plot(community, g)

Edge betweenness ile de birbirine bağlı olan nodeları gruplandırabiliriz. Bu özellikle büyük verisetlerinde işe yarar.

eb_communities <- cluster_edge_betweenness(g)
plot(eb_communities, g)

Yeni node ve edge ekleme

g <- add_vertices(g, 1, name = "New Node")
g <- add_edges(g, c("New Node", "A"))
plot(g, vertex.label = V(g)$name, edge.arrow.size = 0.5)

Edge çıkartma

g <- delete_edges(g, E(g, P = c("New Node", "A")))
plot(g, vertex.label = V(g)$name, edge.arrow.size = 0.5)

Subset etme

sub_g <- induced_subgraph(g, vids = V(g)[degree(g) > 1])
plot(sub_g)

Communityleri daha büyük bir veristeinde de görüntüleyelim

community <- cluster_edge_betweenness(karate)

## Warning in cluster_edge_betweenness(karate): At
## vendor/cigraph/src/community/edge_betweenness.c:498 : Membership vector will be
## selected based on the highest modularity score.

membership(community)

##    Mr Hi  Actor 2  Actor 3  Actor 4  Actor 5  Actor 6  Actor 7  Actor 8 
##        1        1        2        1        3        3        3        1 
##  Actor 9 Actor 10 Actor 11 Actor 12 Actor 13 Actor 14 Actor 15 Actor 16 
##        4        2        3        1        1        2        4        4 
## Actor 17 Actor 18 Actor 19 Actor 20 Actor 21 Actor 22 Actor 23 Actor 24 
##        3        1        4        1        4        1        4        5 
## Actor 25 Actor 26 Actor 27 Actor 28 Actor 29 Actor 30 Actor 31 Actor 32 
##        5        5        6        5        2        6        4        4 
## Actor 33   John A 
##        4        4

plot(community, karate, vertex.label = V(karate)$name, main = "Communities in Karate Network")

Modularity böyle daha büyük verilerde verinin communitylere bölünme derecesini inceler.

modularity(community)

## [1] 0.345299

Clustering coefficent network içindeki nodeların ne derece cluster olduğuna bakar

transitivity(karate, type = "global")

## [1] 0.2556818

Bunu her nod için de yapabiliriz

transitivity(karate, type = "local")

##     Mr Hi   Actor 2   Actor 3   Actor 4   Actor 5   Actor 6   Actor 7   Actor 8 
## 0.1500000 0.3333333 0.2444444 0.6666667 0.6666667 0.5000000 0.5000000 1.0000000 
##   Actor 9  Actor 10  Actor 11  Actor 12  Actor 13  Actor 14  Actor 15  Actor 16 
## 0.5000000 0.0000000 0.6666667       NaN 1.0000000 0.6000000 1.0000000 1.0000000 
##  Actor 17  Actor 18  Actor 19  Actor 20  Actor 21  Actor 22  Actor 23  Actor 24 
## 1.0000000 1.0000000 1.0000000 0.3333333 1.0000000 1.0000000 1.0000000 0.4000000 
##  Actor 25  Actor 26  Actor 27  Actor 28  Actor 29  Actor 30  Actor 31  Actor 32 
## 0.3333333 0.3333333 1.0000000 0.1666667 0.3333333 0.6666667 0.5000000 0.2000000 
##  Actor 33    John A 
## 0.1969697 0.1102941

Networkün yoğunluğuna da bakabiliriz. Bu ise potansiyel bağlantıların ne kadarının gerçekten bağlantı oluşturduğuna bakar. Böylece networkün potansiyel ilişkisine ne kadar ulaştığına bakabiliriz.

edge_density(karate)

## [1] 0.1390374

Eigenvector ise bütün nodelara bir değer atar. Bu atamayı da daha merkezde olan veya daha çok bağlantısı bulnan nodelara daha yüksek değer atayarak yapar.

eigenvector_centrality <- evcent(g)$vector

## Warning: `evcent()` was deprecated in igraph 2.0.0.
## ℹ Please use `eigen_centrality()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

print(eigenvector_centrality)

##         A         B         C         D         E  New Node 
## 0.7807764 0.8337830 0.8337830 1.0000000 0.4682132 0.0000000

Alternatif olarak başka paketler kullanarak da logo oluşturabiliriz.

library(tidygraph)
library(ggraph)
network_data <- tibble::tribble(
  ~from, ~to,
  "A", "B",
  "B", "C",
  "C", "D",
  "D", "A"
)

g_tbl <- as_tbl_graph(network_data)

ggraph(g_tbl, layout = "fr") + 
  geom_edge_link() + 
  geom_node_point() + 
  geom_node_text(aes(label = name), vjust = 1.5) +
  theme_void()

Bunu kendi kullandığım bir veri üstünden anlatayım. Bu veri Twitter ile çalıştığım bir veri. Ülkelerin arasındaki troll saldırılarını, övme veya eleştirmeyi gösteriyor ve kaç defa bu derecede bir saldırı olduğunu belirtiyor. Böylece hem edge width, hem direction hem onun tipi hem de node boyutu ile ilgili analiz yapabiliriz. Şimdi veriyi alıp analiz edelim.

library(readxl)
library(igraph)
library(ggraph)
library(tidyverse)
edges <- read_xlsx("edges.xlsx")
nodes <- read_xlsx("nodes.xlsx")
edges <- edges %>%
  mutate(attack_type = case_when(
    attack_type == "P" ~ "green",
    attack_type == "N" ~ "red"
  ))
used_nodes <- unique(c(edges$source, edges$target))
filtered_nodes <- nodes %>% filter(country %in% used_nodes)
graph <- graph_from_data_frame(d = edges, vertices = filtered_nodes, directed = TRUE)
V(graph)$country <- filtered_nodes$country

ggraph(graph, layout = "fr") + 
  geom_edge_link(aes(color = attack_type, width = attack_density), 
                 arrow = arrow(length = unit(3, "mm"), type = "open"), 
                 end_cap = circle(3, 'mm'), 
                 edge_alpha = 0.7, 
                 lineend = "round") + 
  geom_edge_loop(aes(color = attack_type, width = attack_density), 
                 arrow = arrow(length = unit(3, "mm"), type = "open"), 
                 end_cap = circle(3, 'mm'), 
                 edge_alpha = 0.7, 
                 lineend = "round") + 
  geom_node_point(aes(size = node_size), color = "blue", fill = "lightblue", shape = 21, stroke = 1) + 
  geom_node_text(aes(label = V(graph)$country), repel = TRUE, size = 5, color = "black") + 
  scale_size_continuous(range = c(4, 12)) + 
  scale_edge_width_continuous(range = c(0.5, 2)) + 
  scale_edge_color_manual(values = c("green" = "green", "red" = "red")) + 
  theme_void() + 
  theme(legend.position = "none")

Machine Learning

Bilgisayarlara veriden öğrenmeyi ve böylece gelecekte benzer bir durumla karşılaşılmasıdurumunda neler olabileceğini öğretmemize ve bunları tahmin etmemize yarayan bir yöntemdir. Temel olarak 2 tip vardır.

• Supervised: Label ettiğimiz bir veriyi kullanarak, parametreleri bilgisayara vererek öğrettiğimiz modeldir.

• Unsupervised: Label olamadan verileri gruplama üstüne kuruludur. Bir değişken belirtmeden ortaklık aranarak verileri gruplar ve bu şekilde ayrımı sağlar.

Bu alanın çok klasik bir verisi titanic verisidir. Titanic yolcularının belli parametrelere göre hayatta kalıp kalmadığını tahmin etmeki için kullanılır. Bu örnek supervised learninge bir örnek olacak. Önemli olan bazı değişkenler:

• Survived: Yolcunun hayatta kalıp kalmadığına dair veri. 0 ve 1 değeri alır.
• Pclass: Yolcunun hangi sınıfta olduğunu gösterir. 1, 2, 3 değerlerini alır.
• Sex: Yolcunun cinsiyeti
• Age: Yolcunun yaşı
• Fare: Yolcunun bilete ödedeiği ücret

library(titanic)
library(caret)
library(randomForest)
library(e1071)
library(nnet)
library(xgboost)
data("titanic_train")

Önce missing değerlere bakalım.

colSums(is.na(titanic_train))

## PassengerId    Survived      Pclass        Name         Sex         Age 
##           0           0           0           0           0         177 
##       SibSp       Parch      Ticket        Fare       Cabin    Embarked 
##           0           0           0           0           0           0

median(titanic_train$Age, na.rm = TRUE)

## [1] 28

Age columında 177 missing değer var. 891 veri içinde çok yüksek. Onun için impute dediğimiz yöntemi kullanalım ve median ile dolduralım. Böylece median değişmeden veriyi elde etmiş olduk.

titanic_train$Age[is.na(titanic_train$Age)] <- median(titanic_train$Age, na.rm = TRUE)
colSums(is.na(titanic_train))

## PassengerId    Survived      Pclass        Name         Sex         Age 
##           0           0           0           0           0           0 
##       SibSp       Parch      Ticket        Fare       Cabin    Embarked 
##           0           0           0           0           0           0

median(titanic_train$Age, na.rm = TRUE)

## [1] 28

Cinsiyet değişkenini dummy olarak kodlayalım.

titanic_train$Sex <- ifelse(titanic_train$Sex == "male", 1, 0)

Verinin basit olması için diğer columnları drop edelim.

titanic_train <- titanic_train %>% select(Survived, Pclass, Sex, Age)

Şu an verimiz analize hazır durumda. Şimdi machine learning kısmına başlayabiliriz. Öncelikle verimizi ikiye bölebiliriz. Titanic verisi çoktan ikiye bölünmüş olarak geldiği için bu aşamayı pas geçeceğiz ama verimiz tek halde gelseydi bu bölmeyi yapmak gerekecekti. Yine de nasıl yapılacağını görelim.

set.seed(123)
train_index <- createDataPartition(titanic_train$Survived, p = 0.7, list = FALSE)
train_data <- titanic_train[train_index, ]
test_data <- titanic_train[-train_index, ]

Modelimizi çalıştıralım. Gördüğümüz üzere model aslında lojistik regresyon. Sonuçlara bakalım.

model_log <- train(Survived ~ ., data = train_data, method = "glm", family = "binomial")

## Warning in train.default(x, y, weights = w, ...): You are trying to do
## regression and your outcome only has two possible values Are you trying to do
## classification? If so, use a 2 level factor as your outcome column.

summary(model_log)

## 
## Call:
## NULL
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.989959   0.547376   9.116  < 2e-16 ***
## Pclass      -1.168618   0.142408  -8.206 2.28e-16 ***
## Sex         -2.657501   0.225041 -11.809  < 2e-16 ***
## Age         -0.040044   0.009014  -4.442 8.90e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 833.37  on 623  degrees of freedom
## Residual deviance: 564.30  on 620  degrees of freedom
## AIC: 572.3
## 
## Number of Fisher Scoring iterations: 5

Tahminlerimizi inceleyelim.

pred_log_prob <- predict(model_log, test_data)
pred_log <- ifelse(pred_log_prob > 0.5, 1, 0)
confusionMatrix(as.factor(pred_log), as.factor(test_data$Survived))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 133  27
##          1  34  73
##                                           
##                Accuracy : 0.7715          
##                  95% CI : (0.7164, 0.8205)
##     No Information Rate : 0.6255          
##     P-Value [Acc > NIR] : 2.334e-07       
##                                           
##                   Kappa : 0.5191          
##                                           
##  Mcnemar's Test P-Value : 0.4424          
##                                           
##             Sensitivity : 0.7964          
##             Specificity : 0.7300          
##          Pos Pred Value : 0.8313          
##          Neg Pred Value : 0.6822          
##              Prevalence : 0.6255          
##          Detection Rate : 0.4981          
##    Detection Prevalence : 0.5993          
##       Balanced Accuracy : 0.7632          
##                                           
##        'Positive' Class : 0               
## 

Accuracy için alternatif bir yöntemle bakalım

accuracy_log <- mean(pred_log == test_data$Survived)
print(accuracy_log)

## [1] 0.7715356

Özellikle verinin balansının tam olmadığı durumlarda precision, recall ve f1 gibi değerleri kullanabiliriz. - Precision predicted pozitiflerin ne kadar doğru olduğunu görmemize yardımcı olur. - Recall veya sensitivity actual pozitiflerin ne kadar doğru tahmin edildiğini gösterir. - F-1 score ise bu ikisi arasında bir denge kurar.

precision_log <- posPredValue(as.factor(pred_log), as.factor(test_data$Survived), positive = "1")
recall_log <- sensitivity(as.factor(pred_log), as.factor(test_data$Survived), positive = "1")
f1_log <- 2 * ((precision_log * recall_log) / (precision_log + recall_log))
cat("Precision: ", precision_log, "\n")

## Precision:  0.682243

cat("Recall: ", recall_log, "\n")

## Recall:  0.73

cat("F1-Score: ", f1_log, "\n")

## F1-Score:  0.705314

Son olarak ROC Curve ve AUC değerlerine bakarak da gücünü görebiliriz. Curve sol üst köşeye ne kadar yakınsa o kadar iyi bir modele sahibiz diyebiliriz.

library(pROC)
pred_log_prob <- predict(model_log, test_data, type = "raw")
roc_curve <- roc(test_data$Survived, pred_log_prob)
plot(roc_curve, col = "blue", main = "ROC Curve for Logistic Regression")

auc_log <- auc(roc_curve)
cat("AUC: ", auc_log, "\n")

## AUC:  0.8449401

AUC yorumlamasını şu şekilde yaparız.

• AUC = 0.5 ise bu rastgele seçim yapmaktan farksızdır
• AUC > 0.7 ise iyi bir model olduğunu gösterir
• AUC > 0.9 ise mükemmele yakın bir model olduğunu gösterir.

Bu veri Bostondaki evlerin değerini tahmin etmek için kullanılacak.

library(MASS)
data("Boston")
colSums(is.na(Boston))

##    crim      zn   indus    chas     nox      rm     age     dis     rad     tax 
##       0       0       0       0       0       0       0       0       0       0 
## ptratio   black   lstat    medv 
##       0       0       0       0

Burada hiç missing değer olmadığını görüyoruz.

Verimizi yine ikiye bölelim

set.seed(123)
library(caTools)
split <- sample.split(Boston$medv, SplitRatio = 0.7)
train_data <- subset(Boston, split == TRUE)
test_data <- subset(Boston, split == FALSE)

Analizi yapalım

model_lm <- lm(medv ~ ., data = train_data)
summary(model_lm)

## 
## Call:
## lm(formula = medv ~ ., data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.5797  -2.6006  -0.5752   2.1161  18.1208 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  33.765376   5.622702   6.005 4.75e-09 ***
## crim         -0.117755   0.033906  -3.473 0.000579 ***
## zn            0.044972   0.014919   3.014 0.002761 ** 
## indus         0.040452   0.066552   0.608 0.543693    
## chas          2.457784   0.963830   2.550 0.011194 *  
## nox         -19.711311   4.168890  -4.728 3.28e-06 ***
## rm            4.618992   0.440164  10.494  < 2e-16 ***
## age          -0.008863   0.013928  -0.636 0.524931    
## dis          -1.559286   0.219516  -7.103 6.76e-12 ***
## rad           0.226005   0.072663   3.110 0.002021 ** 
## tax          -0.010045   0.004202  -2.390 0.017351 *  
## ptratio      -0.933469   0.142952  -6.530 2.29e-10 ***
## black         0.005653   0.002966   1.906 0.057479 .  
## lstat        -0.500263   0.055431  -9.025  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.418 on 353 degrees of freedom
## Multiple R-squared:  0.7905, Adjusted R-squared:  0.7828 
## F-statistic: 102.5 on 13 and 353 DF,  p-value: < 2.2e-16

Tahminlerimizi çıkaralım

predictions <- predict(model_lm, test_data)
mse <- mean((test_data$medv - predictions)^2)
rmse <- sqrt(mse)
rss <- sum((test_data$medv - predictions)^2)
tss <- sum((test_data$medv - mean(test_data$medv))^2)
r_squared <- 1 - (rss/tss)
cat("MSE: ", mse, "\n")

## MSE:  31.88877

cat("RMSE: ", rmse, "\n")

## RMSE:  5.647014

cat("R-squared: ", r_squared, "\n")

## R-squared:  0.5366081

R-squared’in ne olduğunu öğrenmiştik. Predictorlarımız dependent variable içindeki değişimin % kaçını ölçüyor.

MSE actual değer ve fitted değerler arasındaki farka bakar. Bunun düşük olması daha iyidir.

RMSE ise hatanın oranını dependent değişken ile aynı ünitlerde olduğu bir birimdir. MSE’ye göre interpret etmesi de bundan dolayı daha kolaydır. Mesela bu önreğiizde RMSE 5-6 arası. Bu da bizim tahminlerimizin gerçek değerlere göre 5-6 dolar arası saptığını gösteriyor.

İyi bir model düşük RMSE, MSE ve yüksek R-squared değerleine sahip olmalıdır. Bu değerler özellikle iki modeli karşılaştırırken hangisinin daha iyi olduğunu anlamamıza yardımcı olur.

SON