Статистички софтвер 1: Визуализација података са ggplot2 пакетом

Стефан Малбашић, Математички факултет, Универзитет у Београду

Графичко представљање података

Графичко представљање података је метод који приказује вредности, извештаје, упоређује резултате или даје предвиђање у току одређених ситуација. То је визуелни приказ података кроз графиконе, а у наставку ћемо користити пакет ggplot2.

install.packages("ggplot2")
library(ggplot2)

Пакет ggplot2 је један од најкоришћенијих и најзаступљенијих пакета у програмском језику R. Има доста различиту синтаксу од синтаксе класичних графика, али највећа предност овог пакета је што се са мало кода праве естетски богати графици.

Пакет се заснива на тзв. “граматици графике”. Графици се праве “сабирањем”, односно додавањем компоненти једну по једну на график. Главни делови графика су:

  • подаци (data): скуп података који се приказује,

  • геометријски објекат (geom_*): тип геометријског објекта којим се приказују подаци,

  • естетски параметри (aes()): естетски атрибути који се могу доделити геометријским објектима.

Графике правимо пратећи следећи шаблон:

ggplot(data = <Podaci>) + 
  <geom_funkcija>(mapping = aes(<estetski parametri>))

Битно је поменути да не постоји подразумевани график у ggplot, већ се мора додати неки геометријски објекат (позивом неке geom_* функције). Позивом само ggplot(data=...) добија се празан график.

ggplot(data = cars)

Сада ћемо приказати неколико врста графика који се могу цртати.

Scatterplotovi

Ови графици приказују податке као тачке у равни. Геометријски објекат, тј. функција која се користи у ове сврхе је geom_point. Користићемо базу података mpg из пакета ggplot2. Прикажимо потрошњу аутомобила (hwy) у односу на запремину мотора (displ).

ggplot(data = mpg, aes(x = displ, y = hwy)) +
  geom_point()

Приметимо да постоји опадајући тренд, тј. што је већи мотор, потрошња горива расте, али имамо и неколико аутлајера од тренда на десном крају графика. Због јаснијег уочавања правилности на графику, додајмо естетске елементе и то наћимо кластере на основу класе возила (комби, SUV, мали ауто, спортски итд), која се налази у променљивој class у бази.

ggplot(data = mpg, aes(x = displ, y = hwy)) +
  geom_point(aes(color = class))

Сада закључујемо да су уочени аутлајери заправо спортски аутомобили који имају много велике моторе, а малу масу, па им је потрошња мања него код осталих.

Даље, можемо направити график тако да и величину тачке ставимо у зависности од тога колико цилиндара има мотор. Опет само додајемо естетски параметар size = cyl.

ggplot(data = mpg, aes(x = displ, y = hwy)) +
  geom_point(aes(color = class, size = cyl))

Јасно уочавамо да велику потрошњу имају мотори са више цилиндара.

Такође, можемо мењати облик тачкица прослеђивањем естетсског параметра shape. Нпр. у зависности од погона (предњи, задњи, сва 4 точка) имамо различит облик. Додатно, можемо повећати тачке додавањем аргумента size у самој функцији.

ggplot(data = mpg, aes(x = displ, y = hwy)) +
  geom_point(mapping = aes(color = class, shape = drv), size = 3)

Још једна од предности ggplot пакета је што се могу користити и неке статистичке функције. Нпр. можемо доцртати линију која описује тренд у подацима коришћењем функције geom_smooth, која прави одређену апроксимацију података (линеарни модел, полиномијални и сл).

ggplot(data = mpg, aes(x = displ, y = hwy)) +
  geom_point(mapping = aes(color = class, shape = drv), size = 3) + 
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Подразумевано је да фунцкија geom_smooth користи методу “loess” као апроксимацију која нам даје нелинеарну оцену тренда. Најједноставнија оцена тренда је да повучемо праву кроз податке коришћењем линеарне регресије. То радимо додавањем аргумента method = "lm".

ggplot(data = mpg, aes(x = displ, y = hwy)) +
  geom_point(mapping = aes(color = class, shape = drv), size = 3) + 
  geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

Линијски графици

За цртање временски уређених података, попут временских серија, стандардни графици су линијски. То је линија која приказује кретање неке променљиве кроз време.

data("economics", package = "ggplot2")

ggplot(economics, aes(x = date, y = psavert)) +
  geom_line(color = "darkgreen") +
  geom_smooth() + 
  labs(x = "Mesec", y = "Stopa stednje") 
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Хистограм

Груба оцена расподеле узорка се може визуализовати хистограмима, где бројимо колико елемената узорка упада у одговарајуће категорије. Користимо геометријски објекат geom_histogram.

Користићемо податке о температури из базе weather из пакета nycflights13, при чему ћемо узети само податке за 25. децембар 2013. Такође, претворићемо температуру из фаренхајта у целзијусе.

library(nycflights13)

# dodajemo kolonu u kojoj je temperatura po celzijusovoj skali
weather$celsius <- (weather$temp - 32) * 5 / 9

weather_today <- weather[weather$month == 12 & weather$day == 25, ]
ggplot(weather_today, aes(x = celsius)) + # dovoljna nam je samo x osa za histogram
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Подразумевано се формира 30 категорија, па је хистограм доста неинформативан у конкретном случају. Можемо променити број категорија или ширину хистограма.

ggplot(weather_today, aes(x = celsius)) + 
  geom_histogram(bins = 10)

Сада немамо проблем са бројем категорија, али хистограм није јасно видљив. Зато ћемо да урадимо бољу визуализацију: поставићемо да границе измећу категориај буду беле боје, а колоне плаве.

ggplot(weather_today, aes(x = celsius)) + 
  geom_histogram(bins = 10, color = "white", fill = "steelblue")

Упоредимо овај хистограм са оним који добијамо директно из R-а.

hist(weather_today$celsius)

У R-у је подразумевано да се број категорија рачуна по формули \(\lceil \log_2n + 1\rceil\), што се добија позивом nclass.Sturges функције.

bin_count <- nclass.Sturges(weather_today$celsius) 

ggplot(weather_today, aes(x = celsius)) + 
  geom_histogram(bins = bin_count, color = "white", fill = "steelblue")

Приказ више графика на истом прозору

Желимо да прикажемо хистограме за различите месеце из претходонг примера са температуром. За то користиму функцију facet_wrap. На графику испод можемо приметити да се разликују расподеле температуре по месецима, што је и очекивано. Расподеле у зимским месецима су померене улево у односу на летње.

ggplot(weather, aes(x = celsius)) +
  geom_histogram(color = "white", fill = "steelblue") +
  facet_wrap(~ month)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

Кутијасти дијаграми (Boxplot)

У неким случајевима бољи начин поређења расподела могу нам омогућити кутијасти дијаграми. Boxplot представља графички приказ 5 сумарних статистика - минимум, 1. квартил, медијана, 3. квартил и максимум. Те вредности добијамо покретањем функције summary на (нумеричком) узорку.

summary(weather_today$celsius)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -7.200  -5.600  -3.300  -3.693  -2.200   0.000

Boxplot се генерише функцијом geom_boxplot. Када цртамо график за један узорак, довољно је дати само y-осу.

ggplot(weather_today, aes(y = celsius)) +
  geom_boxplot()

Веома корисно својство функције geom_boxplot је то што може да користи и x-осу, која мора бити категоричка, тј. фактор. Погледајмо резултат ако додамо месец, претворен у категоричку променљиву.

ggplot(weather, aes(y = celsius, x = factor(month))) +
  geom_boxplot()
## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).

Bar plot

Bar plot или тракасти дијаграм служи за графички приказ категоричких променљивих. На њима се приказује колико постоји података у бази из одговарајуће категорије. Разликују се од хистограма по томе што не сврставају податке по категоријама у зависности од тога ком интервалу припадају њихове нумеричке вредности, већ просто броје чланове категорија. Тракасте дијаграме цртамо позивом функције geom_bar. Такође је довољно дати само x-осу.

ggplot(mpg, aes(x = class)) +
  geom_bar(fill = "forestgreen", color = "black")

Некада нећемо да пребројимо елементе неке категорије, већ имамо у подацима за одређену категоричку променљиву додељену неку вредност. Нпр. у подацима gapminder имамо државе и њихове популације, па можемо да нацртамо популацију по државама. То радимо функцијом geom_col којој поред x-осе додамо и y-осу са одговарајућим вредностима.

library(gapminder)

# napravimo kratak spisak drzava koj cemo gledati
countries <- c("Bosnia and Herzegovina", "Croatia", "Montenegro", "Serbia", "Slovenia")

gapminder2007 <- gapminder[gapminder$year == 2007, ]
gapminder_small <- gapminder2007[gapminder2007$country %in% countries, ]

ggplot(gapminder_small, aes(x = country, y = pop)) +
  geom_col(fill = "forestgreen", color = "black")

Тракасти дијаграми могу бити наслагани, да приказују више вредности за исту категорију. Нацртаћемо тракасти дијаграм где се види популација за 1997, 2002. и 2007. годину.

gapminder567 <- gapminder[(gapminder$year %in% c(1997, 2002, 2007)) &
                            (gapminder$country %in% countries), ]

ggplot(gapminder567, aes(x = country, y = pop)) +
  geom_col(color = "black", mapping = aes(fill = factor(year)))

Ови графици нису лаки за читање јер је тешко упоредити величине правоугаоника разних боја. Кориснија алтернатива се добије када се дода аргумент position = "dodge".

ggplot(gapminder567, aes(x = country, y = pop)) +
  geom_col(color = "black", mapping = aes(fill = factor(year)), position = "dodge")

Све наведено за geom_col може се одрадити за geom_bar ако се ради пребројавање одговарајућих категорија. Нпр, можемо приказати број летова из пакета nycflights13 по компанијама у зависности од полазишта.

ggplot(flights, aes(x = carrier, fill = origin)) +
  geom_bar(position = "dodge")

График густине

График густине представља расподелу узорка у континуалном интервалу или временском периоду. То је заправо верзија хистограма која глатко исцртава вредности како би се избегли груби прелази.

#Podešavamo alfa transparentnost i bojimo gustine u odnosu na broj cilindara 
 
ggplot(mpg, aes(cty)) + 
  geom_density(aes(fill = factor(cyl)), alpha = 0.8) + 
  labs(title= "Density plot", 
       subtitle = "City Mileage Grouped by Number of cylinders",
       caption = "Source: mpg",
       x = "City Mileage",
       fill = "# Cylinders")

Кружни дијаграм

Кружни дијаграм или питица је кружни статистички график који је подељен на парчиће. Површина сваког парчића је пропорционална количини коју представља.

#Sada ćemo piticom da predstavimo odnos broja vozila svih ponudjenih klasa
 
pie <- ggplot(mpg, aes(x = "", fill = factor(class))) + 
   geom_bar(width = 1) +
   theme(axis.line = element_blank(), #izbacujemo x-osu i y-osu 
         plot.title = element_text(hjust = 0.5)) + #podešavamo i centriramo naslov
   labs(fill="class", 
        x = NULL, 
        y = NULL, 
        title = "Pie Chart of class", 
        caption = "Source: mpg")
 
pie + coord_polar(theta = "y")

Хијерархијски дентограм

Помоћу њих групишемо објекте, који су слични у неком смислу, у групе које називамо кластерима и правимо хијерархијску везу између њих.

#U narednom primeru ćemo hijerarhijski predstaviti broj uhapšenih kriminalaca na 100000 ljudi u svim državama SAD-a.

library(ggdendro)
theme_set(theme_bw())
 
hc <- hclust(dist(USArrests), "ave") #hijerarhijski klastering, metoda "average"
 
ggdendrogram(hc, rotate = T, size = 2)