avaialable from (https://rpubs.com/staszkiewicz/EX_7_NA_PL)

Wprowadzenie

Spróbujemy sprawdzić gdzie fizycznie są umieszczone Banki w USa.

Do tego potrzebujemy adresy banków (głównie siedzib) oraz inforamcje o ich umiejscowieniu geograficznym w postaci długości i szerokości geograficznej (geocodowanie).

Przestrzenne obrazowanie Polega na przypisaniu obiektowi miejsca w przestrzeni. Sposód wielu pakietów dostępnych do geocodowania, my skorzystamy z tmaptools

#rm(list = ls()) # Gdyby była potrzeba wyczyszczenia środowiska
#install.packages("tmaptools") # jeśli po raz pierwszy instalujemy bibliotekę

library(tmaptools) # ładowanie pakietu

## Warning: pakiet 'tmaptools' został zbudowany w wersji R 4.1.3

Jak zwykle, wykorzystamy dane pierwotne z artykułu Audit fee and communication sentiment. Badania Ekonomiczne-Ekonomska Istraživanja. https://doi.org/10.1080/1331677X.2021.1985567 Staszkiewicz and Karkowska (2021).

# z bazowych funkcji systemu wczytamy klasyczny plik z csv ale w taki sposób
#   że wybierzemy z okienka umiejscowienie pliku na włanym komuterze 
# dlaego zagnieżdzamy polecenie "file.choose()"

bank <-  read.csv(file.choose())

Podstawy

przypomnijmy sobie strukturę naszego plikut str(bank) Mamy dane adresowe rozbite do kilku kolumn zmiennych a mianowicie: Companyx, Bus_Strreet_1x, Bus_Strreet_2x, Cityx, Countryx, State_Codex,State_Namex,Region, Zipx - te wszyskie zmienne są ciągami alafanumerycznymi (znakowymi). Nasz plik ma 5356 obeserwacji za okres table(bank$Year_Ended) od 2000 do 2016. Ponieważ siedziby banków raczej zbyt dynamicznie się nie zmieniają to dla naszych potrzeb zanalizujemy rok 2016. Wpierw potrzebujemy ograniczyć nasz zbiór tylko do roku 2016

bank2016<- bank[bank$Year_Ended == 2016,]

bank2016<- bank[bank$Year_Ended == 2016,]

Proszę zauważyć, iż potrzebujemy przecinka, by wskazć, że wybieramy wszystkie kolumy.

Porównajmy oba obiekty bank i bank2016 otóż, mamy po 197 zmiennych w obu ramkach danych lecz zredukowaliśmy liczbę obserwacji (w tym przypaku banków) od 359.

Spróbujmy dla ustalenie uwagi wykonać geokodowanie jednego adresu - pierwszego banku w naszej nowej ramce danych tj. bank2016

View(bank2016[1,])

Co możemy zrobić to połączyć poszczególne elementy adresu w jeden obiekt i przypisać takiem adresowi długość i szerokość geograficzną.

połączenie adresu

ASB<-bank2016[1,] # definujemy dane o banku

2 Na prostym wektorze dane o adresie mamy w kolumnach od 7 do 14, lecz dla naszych potrzeb waży są nastepujące elementy, nazwa stanu, region i kod pocztowy tj. pozycjie od 12 do 14.

więc wykorzystamy funkcję paste() z argumentem collapse =’’ by skleić odrębne komórki do jednej i przypiszemy je do obiektu adresASB w taki sposób:

adresASB<-paste(ASB[c(12:14)],collapse=' ')

mając już adre potrzebujemy do niego przypisać długość i szerokość geograficzną ponieważ uruchomiliśmy biblotekę library(tmaptools) to możemy bezpośrednio zgeokodować adresASB

library(tmaptools)
gk<-geocode_OSM(adresASB)
gk

## $query
## [1] "WISCONSIN US Midwest 54301"
## 
## $coords
##         x         y 
## -88.49582  44.24636 
## 
## $bbox
##      xmin      ymin      xmax      ymax 
## -88.49677  44.24636 -88.49532  44.24637

Proszę zauważyć, iż geokodowanie zwraca nam obiekt w postacji listy

class(gk)

## [1] "list"

i jak popatrzym na strukturę tej listy to mamy następujące elementy

str(gk)

## List of 3
##  $ query : chr "WISCONSIN US Midwest 54301"
##  $ coords: Named num [1:2] -88.5 44.2
##   ..- attr(*, "names")= chr [1:2] "x" "y"
##  $ bbox  : 'bbox' Named num [1:4] -88.5 44.2 -88.5 44.2
##   ..- attr(*, "names")= chr [1:4] "xmin" "ymin" "xmax" "ymax"
##   ..- attr(*, "crs")=List of 2
##   .. ..$ input: chr "EPSG:4326"
##   .. ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223"| __truncated__
##   .. ..- attr(*, "class")= chr "crs"

zapytanie: $query
długość i szerokość $coords
oraz powierzchnia $bbox

W tym przykładzie nie będziemy korzystać z powierzchni, potrzebujemy jedynie długość (x) i szerokość (y) geograficzną.

Więc dodajmy x i y jako zmienne do naszej ramki. Tutaj proszę poeksperymentować w jaki spób adresujemy w liście elemnty: np gk, gk[1], gk$query,gk$query[1]

i połączymy do ASB nasze x i y. ASB<-c(ASB,gk$coords[1:2])

operacje na całej ramce danych

Ten sam tok rozumowania możem powtórzyć na całej ramce danych i przypiszmy nasze wynki do danych

bank2016$adr <- paste(bank2016$Countyx,bank2016$Region, bank2016$Zipx, sep=” “)

bank2016$adr <- paste(bank2016$Region, bank2016$Zipx, sep=" ")

Teraz zgeokodujemy nasz adres, ten proces będzie trochę trwał i niesetety zwróci nam kilka błędów. Zapiszmy wyniki do obiektu XY

## No results found for "US Southwest 75201".

## No results found for "US Mid Atlantic 14203".

## No results found for "US Mid Atlantic 15212".

## No results found for "US Mid Atlantic 18964".

## No results found for "US Mid Atlantic 26003".

## No results found for "US Mid Atlantic 41502-2947".

## No results found for "US Mid Atlantic 12302".

## No results found for "US Mid Atlantic 17604".

## No results found for "US Southwest 75701".

## No results found for "US Mid Atlantic 15907".

## No results found for "US Mid Atlantic 15701".

## No results found for "US Mid Atlantic 7632".

## No results found for "US Mid Atlantic 19335".

## No results found for "US Mid Atlantic 15222-2401".

## No results found for "US Mid Atlantic 7470".

## No results found for "US Mid Atlantic 17059-0066".

## No results found for "US Mid Atlantic 17325".

## No results found for "US Mid Atlantic 17740".

## No results found for "US Mid Atlantic 12801".

## No results found for "US Mid Atlantic 15701".

## No results found for "US Mid Atlantic 22657".

## No results found for "US Mid Atlantic 13214".

## No results found for "US Mid Atlantic 17201-0819".

## No results found for "US Mid Atlantic 25313".

## No results found for "US Mid Atlantic 25301".

## No results found for "US Mid Atlantic 16830".

## No results found for "US Mid Atlantic 18603".

## No results found for "US Mid Atlantic 16933".

## No results found for "US Mid Atlantic 11545".

## No results found for "US Mid Atlantic 22853".

## No results found for "US Mid Atlantic 23663".

## No results found for "US Mid Atlantic 24541".

## No results found for "US Mid Atlantic 18951-9005".

## No results found for "US Mid Atlantic 14902".

## No results found for "US Mid Atlantic 21550".

## No results found for "US Mid Atlantic 13815".

## No results found for "US Midwest 43502".

## No results found for "US Mid Atlantic 24060".

## No results found for "US Mid Atlantic 19010".

## No results found for "US Mid Atlantic 17403".

## No results found for "US Mid Atlantic 16901".

## No results found for "US Mid Atlantic 26836".

## No results found for "US Mid Atlantic 20832".

## No results found for "US Mid Atlantic 17257".

## No results found for "US Mid Atlantic 19801".

## No results found for "US Mid Atlantic 10013".

## No results found for "US Mid Atlantic 19102".

## No results found for "US Mid Atlantic 40206".

## No results found for "US Mid Atlantic 14006".

## No results found for "US Mid Atlantic 11932".

## No results found for "US Mid Atlantic 7438".

## No results found for "US Mid Atlantic 20601".

## No results found for "US Mid Atlantic 16373".

## No results found for "US Mid Atlantic 24605".

## No results found for "US Mid Atlantic 14569".

## No results found for "US Mid Atlantic 21404".

## No results found for "US Mid Atlantic 17061".

## No results found for "US Mid Atlantic 22611".

## No results found for "US Mid Atlantic 23219".

## No results found for "US Mid Atlantic 25702".

## No results found for "US Mid Atlantic 21227".

## No results found for "US Midwest 45631".

## No results found for "US Mid Atlantic 11590".

## No results found for "US Mid Atlantic 15237".

## No results found for "US Mid Atlantic 23181".

## No results found for "US Mid Atlantic 8809".

## No results found for "US Mid Atlantic 40202".

## No results found for "US Mid Atlantic 11556".

## No results found for "US Mid Atlantic 40362-0157".

## No results found for "US Mid Atlantic 08753-8396".

## No results found for "US Mid Atlantic 11201".

## No results found for "US Mid Atlantic 14850".

## No results found for "US Mid Atlantic 24210".

## No results found for "US Mid Atlantic 18431".

## No results found for "US Mid Atlantic 10027-4512".

## No results found for "US Mid Atlantic 7432".

## No results found for "US Mid Atlantic 7416".

## No results found for "US Mid Atlantic 22482".

## No results found for "US Mid Atlantic 21601-3013".

## No results found for "US Mid Atlantic 18512".

## No results found for "US Mid Atlantic 42440".

## No results found for "US Mid Atlantic 7921".

## No results found for "US Mid Atlantic 18503".

## No results found for "US Mid Atlantic 10901".

## No results found for "US Mid Atlantic 12414".

## No results found for "US Mid Atlantic 42103".

## No results found for "US Southwest 75201".

## No results found for "US Mid Atlantic 20186".

## No results found for "US Midwest 45830".

## No results found for "US Midwest 50010".

## No results found for "US Mid Atlantic 8512".

## No results found for "US Mid Atlantic 24260".

## No results found for "US Mid Atlantic 20191".

## No results found for "US Mid Atlantic 7306".

## No results found for "US Mid Atlantic 7002".

## No results found for "US Mid Atlantic 20716".

## No results found for "US Mid Atlantic 26554-2777".

## No results found for "US Mid Atlantic 23113".

## No results found for "US Mid Atlantic 41702".

## No results found for "US Mid Atlantic 8080".

## No results found for "US Mid Atlantic 23233".

## No results found for "US Mid Atlantic 22911".

## No results found for "US Mid Atlantic 8901".

## No results found for "US Mid Atlantic 14048".

## No results found for "US Mid Atlantic 7724".

## No results found for "US Mid Atlantic 40223".

## No results found for "US Midwest 54701".

## No results found for "US Mid Atlantic 15219".

## No results found for "US Mid Atlantic 18360".

## No results found for "US Mid Atlantic 21043".

## No results found for "US Mid Atlantic 7024".

## No results found for "US Mid Atlantic 18966".

## No results found for "US Mid Atlantic 17522-0457".

## No results found for "US Mid Atlantic 18017".

## No results found for "US Mid Atlantic 24011".

## No results found for "US Mid Atlantic 16365".

## No results found for "US Mid Atlantic 10016".

## No results found for "US Mid Atlantic 19610".

## No results found for "US Mid Atlantic 7095".

## No results found for "US Southwest 75225".

## No results found for "US Mid Atlantic 41101".

## No results found for "US Mid Atlantic 11753".

## No results found for "US Mid Atlantic 19301".

## No results found for "US Mid Atlantic 21286".

## No results found for "US Midwest 46037".

## No results found for "US Southwest 75069".

## No results found for "US Mid Atlantic 22911".

## No results found for "US Mid Atlantic 19145".

## No results found for "US Mid Atlantic 17110".

## No results found for "US Mid Atlantic 7078".

## No results found for "US Mid Atlantic 15320".

## No results found for "US Mid Atlantic 13126".

## No results found for "US Mid Atlantic 19103".

## No results found for "US Mid Atlantic 7004".

## No results found for "US Mid Atlantic 21050".

## No results found for "US Mid Atlantic 24091".

## No results found for "US Mid Atlantic 10462".

## No results found for "US Mid Atlantic 15237".

Geokodowanie jest dość żmudnym procesem i czasami wymaga iteracyjnego podejścia dla naszych celów, po prostu pominiemy te obiekty, dla których nie ustaliśmy X i Y.

nowa<-merge(bank2016,XY, by.x = "adr",by.y="query")
# by połączyć ramki, trzeba wskazć kolumnę (mny) w pierwszej i drugiej, wg których będziemy łączyć ramki danch

Proszę zauważyć że wielkośći obiektów po połączeniu będie inna Bank20166 miał 259 obserawcji, XY 221 zaś po mapowaniu nowa 243 obserwacje. Dlaczego?

Wizualizacja danych na mapie

Wykorzystamy pakiet maps i ggplot by zwizualizować skąd pochodzą banki (gdzie mają siedziby), które notowane są na rynku publiczny w USA w roku 2016

# właczenie biblioteki
library(maps)

## Warning: pakiet 'maps' został zbudowany w wersji R 4.1.3

library(ggmap)

## Ładowanie wymaganego pakietu: ggplot2

## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.1.3

## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.

## Please cite ggmap if you use it! See citation("ggmap") for details.

# kontur świata
mapWorld <- borders("world", colour="gray50", fill="white")

# kontur stanów - nie potrzebny
#states<- borders("state",colour="gray50") 

#wygenrowanie kontury
mpw<-ggplot()+mapWorld
# dodanie punków z naszej bazy

mpw <- mpw+ geom_point(data=nowa, aes(x=lon, y=lat) ,color="blue",alpha=0.5, size=2)

#wyświetlenie zbudowanej grafiki

mpw

Zadanie

Pokaż banki mające siedzib tylko w USA, na czerwono jeśli był zintegrowany raport , na niebiesko jeśli takiego raportu nie było, na czerwono zaznacz banki audytowane przez KMPG zas różnymi kształtami, różne typy emitentów (Accelerated, not-Accelerated itd). Wynik powienien być podoby jak ten poniżej

# mpw<-ggplot()+mapWorld
# mpw <- mpw+ geom_point(data=nowa, aes(x=lon, y=lat), ,color=ifelse(nowa$Is_Integrated_Audit=="Yes","Red","Blue"),alpha=0.5, size=2, shape =factor(nowa$Filer_Status), show.legend = "Yes" )
# mpw
# 
# # to jest wersja, gdzie mapowanie jest na poziomie osi
# 
# mpw<-ggplot(data=nowa, aes(x=lon, y=lat, ,color=ifelse(nowa$Is_Integrated_Audit=="Yes","Red","Blue"),alpha=0.5, size=2, shape =factor(nowa$Filer_Status), show.legend = "Yes" ))+mapWorld
# mpw <- mpw+ geom_point()
# mpw
# 
# #do tego dodamy ograniczenie dla USA
# mpw<-ggplot(data=nowa, aes(x=lon, y=lat, ,color=ifelse(nowa$Is_Integrated_Audit=="Yes","Red","Blue"),alpha=0.5, size=2, shape =factor(nowa$Filer_Status), show.legend = "Yes" ))+mapWorld
# mpw <- mpw+ geom_point()
# mpw<- mpw + xlim(-130,-60) + ylim(20, 50)
# mpw
# 
# # i jeszcze zróżnicujemy wielść `size` w zależności od tego czy bank był badany przez KPMG czy nie
# mpw<-ggplot(data=nowa, aes(x=lon, y=lat, ,color=ifelse(Is_Integrated_Audit=="Yes","Red","Blue"),alpha=0.5, size=ifelse(Auditorx=="KPMG LLP","6","2"), shape =factor(Filer_Status), show.legend = "Yes" ))+mapWorld
# mpw <- mpw+ geom_point()
# mpw<- mpw + xlim(-130,-60) + ylim(20, 50)
# mpw
# 
# 
# dodatkowo chcemy wiedzieć które z banków jakim były typem emitenta `factor(nowa$Filer_Status)`
# 
# 
# 
# 
# 
# # teraz skupimy się tylko na USA podając szergokość i długość geograficzną świata którą chcemy przybliżyć
# mpw + xlim(-130,-60) + ylim(20, 50)

Zadanie 2

Spróbujmy terza połączyć siedzibę audytorów z siedzibą banku. Do tego potrzebujemy geokodowanie audytorów, ponieważ nie mamy adresów, spróbujemy tego dokonać na podstawie nazwy i wpiszemy to do obiekut audGC, tj. audGC<-geocode_OSM(nowa$Auditorx) a następnie połączym do ramki danych tytle tylko że nie możemy skorzystać z poprzedniej metody np. nowaAU<-merge(nowa,audGC, by.x = "Auditorx",by.y="query") bo byśmy zmulitplikowali obseracje z uwagi na to że w obiekcie audGC mamy takie same nazwy zmiennych jakie już są w obiekcie nowa, stąd wygenerujemy nowy obiekt, będący ramką danych któRy anzwiem audLL i zmienimy nazwy zmeinnych lat i log na Alat i Alot a mianowicie audLL<-data.frame(aud =audGC$query,Alat= audGC$lat,Alon =audGC$lon) a dopiero teraz połączymy - no i tak nie wyszło dlaczego? Bo w augGC są powtarzające się rekordy

nrow(audGC)==length(unique(audGC$query)) stąd potrzebujem zredukować audGC jedynie do rekordów unikatowych liczba rzędów w obiekcienrow(audGC) jest inne niż liczab unikatowych wartości length(unique(audGC$query)) Dlateko połączymy tylko unikatowe wartości nowaAU<-merge(nowa,unique(audLL), by.x = "Auditorx",by.y="aud") W tym przypadku zgeokodowaliśmy zaledwie 88 obseracji dla 15 firm arudytowych levels(factor(nowaAU$Auditorx)) Mając już przygotowane dane chcemy połączyć graficzne banki z ich firmami audytowymi.

#
#Te rozwiązanie jest zrobione na podstawie galerii i kodu z tej strony: 
# https://www.data-to-viz.com/story/MapConnection.html

# przypiszę dane do obiektu data

#data <-nowaAU
# 
# # Download NASA night lights image
# download.file("https://www.nasa.gov/specials/blackmarble/2016/globalmaps/BlackMarble_2016_01deg.jpg", 
# destfile = "IMG/BlackMarble_2016_01deg.jpg", mode = "wb")
# 
# # Load picture and render
# library(jpeg) # readJPGEG
# library(grid)# rasterGrob
# earth <- readJPEG("IMG/BlackMarble_2016_01deg.jpg", native = TRUE)
# earth <- rasterGrob(earth, interpolate = TRUE)
# 
# # Count how many times we have each unique connexion + order by importance
# library(dplyr) # trzeba wywołać bibliotekę, bo z jakiegoś powodu nie chce działać
# summary=data %>% 
#   dplyr::count(lat,lon,Auditorx) %>%
#   arrange(n)
# 
# # A function that makes a dateframe per connection (we will use these connections to plot each lines)
# 
# library(geosphere) # for gcIntermediate
# 
# data_for_connection=function( lon,lat, Alon, Alat, group){
#   inter <- gcIntermediate(c(lon, lat), c(Alon, Alat), n=50, addStartEnd=TRUE, breakAtDateLine=F)             
#   inter=data.frame(inter)
#   inter$group=NA
#   diff_of_lon=abs(lon) + abs(Alon)
#   if(diff_of_lon > 180){
#     inter$group[ which(inter$lon>=0)]=paste(group, "A",sep="")
#     inter$group[ which(inter$lon<0)]=paste(group, "B",sep="")
#   }else{
#     inter$group=group
#   }
#   return(inter)
# }
# 
# # Utworzenie kompletnej ramki danych z punktami wszystkich linii, które mają być wykonane.
#  
# data_ready_plot=data.frame()
# for(i in c(1:nrow(summary))){
#   tmp=data_for_connection(summary$lon[i], summary$lat[i], summary$Alon[i], summary$Alat[i] , i)
#   tmp$homecontinent=summary$Auditorx[i]
#   tmp$n=summary$n[i]
#   data_ready_plot=rbind(data_ready_plot, tmp)
# }
# data_ready_plot$homecontinent=factor(data_ready_plot$Auditorx, levels= levels(as.factor(data_ready_plot$Auditorx)))
# 
# # Plot
# library(ggplot2)
# p <- ggplot() +
#   annotation_custom(earth, xmin = -180, xmax = 180, ymin = -90, ymax = 90) +
#   geom_line(data=nowaAU, aes(x=lon, y=lat,  colour="yellow", alpha=1), size=0.6) +
#   scale_color_brewer(palette="Set3") +
#   theme_void() +
#   theme(
#         legend.position="none",
#         panel.background = element_rect(fill = "black", colour = "black"), 
#         panel.spacing=unit(c(0,0,0,0), "null"),
#         plot.margin=grid::unit(c(0,0,0,0), "cm"),
#   ) +
#   ggplot2::annotate("text", x = -150, y = -45, hjust = 0, size = 11, label = paste("Where surfers travel."), color = "white") +
#   ggplot2::annotate("text", x = -150, y = -51, hjust = 0, size = 8, label = paste("data-to-viz.com | NASA.gov | 10,000 #surf tweets recovered"), color = "white", alpha = 0.5) +
#   #ggplot2::annotate("text", x = 160, y = -51, hjust = 1, size = 7, label = paste("Cacedédi Air-Guimzu 2018"), color = "white", alpha = 0.5) +
#   xlim(-180,180) +
#   ylim(-60,80) +
#   scale_x_continuous(expand = c(0.006, 0.006)) +
#   coord_equal() 
# 
# # Save at PNG
# ggsave("IMG/Surfer_travel.png", width = 36, height = 15.22, units = "in", dpi = 90)

Staszkiewicz, Piotr, and Renata Karkowska. 2021. “Audit Fee and Banks’ Communication Sentiment.” Economic Research-Ekonomska Istraživanja, October, 1–21. https://doi.org/10.1080/1331677x.2021.1985567.

EX. 7 Analiza przestrzenna - firmy i zarządzający

Piotr Staszkiewicz

23 01 2022

Wprowadzenie

Podstawy

operacje na całej ramce danych

Wizualizacja danych na mapie

Zadanie