Vizualizacija u R
Amra
September 13, 2016
Ovo je R Markdown dokument. Markdown je jednostavno formatirana sintaksa za pisanje HTML, PDF i MS Word dokumenta.
Za više detalja posjetite http://rmarkdown.rstudio.com.
Ovo je cetvrta nedelja i radimo malo napredniju vizualizaciju u odnosu na trecu nedelju. Za to nam treba paket ggplot2. Ako nemamo instaliran ggplot onda je dobro vrijeme da to sada uradimo:
install.packages("ggplot2")
Vizualizacija u R u koristeći paket ggplot2
Za više o paketu ggplot2 u consoli ukucajte: ??ggplot2.
Prvo učitavamo bazu u Rmarkdown, na isti način kao i u Rstudio consoli samo sto to radimo u chunks (gornji desni ugao ili istovremeno pritisneite ctrl+alt+i).
Prije svake analize, potrebno je da ukucamo summary() ili str (), da vidimo s čime raspolažemo.
str(WHO)## 'data.frame': 194 obs. of 13 variables:
## $ Country : chr "Afghanistan" "Albania" "Algeria" "Andorra" ...
## $ Region : chr "Eastern Mediterranean" "Europe" "Africa" "Europe" ...
## $ Population : int 29825 3162 38482 78 20821 89 41087 2969 23050 8464 ...
## $ Under15 : num 47.4 21.3 27.4 15.2 47.6 ...
## $ Over60 : num 3.82 14.93 7.17 22.86 3.84 ...
## $ FertilityRate : num 5.4 1.75 2.83 NA 6.1 2.12 2.2 1.74 1.89 1.44 ...
## $ LifeExpectancy : int 60 74 73 82 51 75 76 71 82 81 ...
## $ ChildMortality : num 98.5 16.7 20 3.2 163.5 ...
## $ CellularSubscribers : num 54.3 96.4 99 75.5 48.4 ...
## $ LiteracyRate : num NA NA NA NA 70.1 99 97.8 99.6 NA NA ...
## $ GNI : num 1140 8820 8310 NA 5230 ...
## $ PrimarySchoolEnrollmentMale : num NA NA 98.2 78.4 93.1 91.1 NA NA 96.9 NA ...
## $ PrimarySchoolEnrollmentFemale: num NA NA 96.4 79.4 78.2 84.5 NA NA 97.5 NA ...Raspolazemo sa 194 obzervacije i 13 varijabli. Imamo numericke varijable, kategoricke, i integer (broj stanovnika = $Population). Baza se odnosi na demografske, ekonomske karakteristike i socijlane karakteristike po zemljama (194 zemlje). Izvor baze je Svjetska zdravstvena organizacija.
Prije nego što krenemo da radimo grafički prikaz moramo da znamo šta želimo da prikažemo. Tako sada želimo da vidmo odnos između dvije varijable a to je FertilityRate tj. stopa fertiliteta i GNI index.
Nakon instalacije paketa otvaramo paket sa funkcijom library ().
Inace ggplot2 je paket koji pravi grafove u Ru u slojevima i zbog toga dodajemo pluseve.
library (ggplot2) #warning=FALSE smo stavili tako da nemamo warninge u rezultatima rmarkdownaScaterplot
Scatterplots: je vrsta grafova koja stavlja u odnos dvije numericke varijable.
plot = ggplot(WHO, aes(x= GNI, y = FertilityRate)) #aes je skracenica za aestethics.
scatterplot = plot + geom_point() #geom_point znaci da su tacke odnosno skaterplot.
scatterplot 
I na kraju mozete da mijenjate boje i uljepšavate grafove.
plot + geom_point(color= "blue", size = 3, shape = 17) + ggtitle ("Fertility Rate vs. Gross National Income") #size i shape odnose se na vrstu i velicinu *tackica*.
Spasavamo scaterplot u varijablu fertilityGNIplot. Nakon toga mozemo i da ga spasavamo, i prikazemo kad zelimo ali i da mu dodajemo neke druge karakteristike.
fertilityGNIplot = plot+geom_point(color="darkred", size=3, shape=8)+ggtitle ("FertilityRate vs. GNI")Ako želimo da obojimo tačkice u odnosu na regiju (note Regija je kategorička farijabla) onda je kod:
ggplot(WHO, aes (x=GNI, y=FertilityRate, color=Region))+geom_point()
A kada je numerička varijabla, grafikon izgleda…
ggplot(WHO, aes (x=GNI, y=FertilityRate, color=LifeExpectancy))+geom_point()
Dodavanje linija scatterplotu: tj. regresionih linija
Fertility rate i Under 15 su u nelinearonoj vezi i iz tog razloga ga logujemo, da bismo dobili linearnu vezu
library (gridExtra) # da odredimo raspored grafova ako zelimo vise odjednom
NL <- ggplot(WHO,aes(x=FertilityRate, y=Under15))+ geom_point()
# a mozemo i da logujemo pa dobijemo linearnu vezu:
L<- ggplot(WHO,aes(x=log(FertilityRate), y=Under15))+ geom_point()
grid.arrange(NL, L, ncol=2)
Posto je linearna veza mozemo da dodamo linearni model
model = lm(Under15~log(FertilityRate), data = WHO)
summary (model) # provjerava signifikantnost modela##
## Call:
## lm(formula = Under15 ~ log(FertilityRate), data = WHO)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.3131 -1.7742 0.0446 1.7440 7.7174
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.6540 0.4478 17.09 <2e-16 ***
## log(FertilityRate) 22.0547 0.4175 52.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.65 on 181 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.9391, Adjusted R-squared: 0.9387
## F-statistic: 2790 on 1 and 181 DF, p-value: < 2.2e-16Dodavanje regresione linije:
ggplot(WHO,aes(x=log(FertilityRate), y=Under15))+ geom_point() + stat_smooth(method = "lm", level = 0.99, color="orange") #level, mi mozemo da odredimo koji nivo signifikantnosti zelimo
Vise linija na jednom grafu ili kombinacija linija tacki…
ggplot(WHO, aes(x= FertilityRate, y= values, color=variable))+geom_line(aes(y = Under15, col = "Under15",group=1)) + geom_line(aes(y = Over60, col = "Over60", group=2))+geom_point(aes(x=FertilityRate,y=ChildMortality, col="ChildMortality" ))
Bazirajmo se samo na Evropu
WHOEurope<- WHO[WHO$Region=="Europe",]Sada uradimo opet scatterplot i pokusajmo da uradimo labels, ali i ne za sve nego pojedine, tj. samo tamo gdje je vise starih od mladih.
ggplot(WHOEurope, aes(x= FertilityRate, y= ChildMortality, label=Country))+ geom_point() +geom_text(aes(label=ifelse(Over60>Under15,as.character(Country),'')),hjust=0,vjust=0) #hjust i vjust su justification tj. poravnanje teksta
Taj graf nam daje mnogo outlayera, tj. imamo mnogo zemalja koje imaju dosta visoku stopu djecijeg mortaliteta ali i visoku stopu Fertiliteta, sto nam otezava analizu. Prvo provjerimo da li najveci outlayer ima smisla i ko je to>
which.max(WHOEurope$ChildMortality) #dobijemo redni broj## [1] 47WHOEurope[47,1] #47 je broj reda a 1 kolone i dobijemo drzavu## [1] "Tajikistan"Sada uradimo novu tabelu sa WHO evropom i ChildMortality<15
WHOEurope15<- WHOEurope[WHOEurope$ChildMortality<15,] # note: ne moramo da stavljamo navodnike za brojevekada Ukucamo slijedeci kod imamo nekoliko varijabli prikazanih. Fertilitya rate vs. chiled mortaltity, + oni koji imaju stariju populaciju vise od mladje. Tako vidimo da BiH je u dosta nepovoljnoj situaciji, ima najvisi fertitlitet, child mortalitiy dosta visok s obyirom na niyak fertilitet. U grupi je zemalja sa maltom, srbijom ali oni imaju veci fertilitet+stanovnistvo stari.
library (gridExtra) #potrebno je za grid arrange
plotWHO1 = ggplot(WHOEurope15, aes(x= FertilityRate, y= ChildMortality, label=Country))+ geom_point() + geom_text(aes(label=ifelse(Over60>Under15,as.character(Country),'')),hjust=0,vjust=0)
plotWHO2 = ggplot(WHOEurope15, aes(x= FertilityRate, y= ChildMortality, label=Country))+ geom_point() + geom_text(aes(label=ifelse(Over60<=Under15,as.character(Country),'')),hjust=0,vjust=0) #hjust i vjust se odnose na tekst na grafiku kao labels
grid.arrange (plotWHO1, plotWHO2,ncol=2)
Posto smo spasili plot pod plotWHO dodacemo i zadnje elemente, naziv…
plotWHOTitle = plotWHO1 + ggtitle("Stopa fertiliteta i Mortalitet kod djece 2014.g.") + xlab("Stopa Fertiliteta") + ylab ("Mortalitet kod djece")
plotWHOTitle
dodajemo jos boje..
plotWHOTitle + geom_point(colour = "grey70", size = 2)
Barplotovi
library (ggplot2)
ggplot (WHOEurope15, aes(x= Country, y=FertilityRate, fill=Country)) + geom_bar(stat = "identity")# identity koristimo kada nije counts u Y osi, to jest da se varijabla ponavlja više puta npr BiH onda bi bili counts i bilo bi bez stat = "identity"
Ako ne zelimo da nam prikaze sve zemlje u bazi WHOEurope15, nego npr samo prvih deset zemalja onda kucamo:
ggplot (WHOEurope15[1:10,], aes(x= Country, y=FertilityRate, fill=Country)) + geom_bar(stat = "identity")+ theme(legend.position="none")
Ako zelimo da poredamo po velicini onda prvo sortiramo dataframe (Sortiranje smo radili na pocetku prva ili druga semdica): sortiramo data frame
sortWHOEurope15 <- WHOEurope15[order(WHOEurope15$FertilityRate),]Nakon toga napravimo Order zemalja prema fertility Rate jer ggplot2 automatski reda zemlje po abecednom redu.
orderWHOEurope15 <- transform(sortWHOEurope15, Country=reorder(Country, -FertilityRate))ggplot (orderWHOEurope15[1:10,], aes(x= Country, y=FertilityRate, fill=Country)) + geom_bar(stat = "identity")
Trenutno imamo dva problema sa ovim grafom. Nazivi zemalja se preklapaju i imamo legendu a ona nam nije potrebna.
ggplot (orderWHOEurope15[1:15,], aes(x= Country, y=FertilityRate, fill=Country)) + geom_bar(stat = "identity") + theme(axis.title.x = element_blank(), axis.text.x = element_text(angle = 45,hjust=1),legend.position="none") #Legend position izbacuje legendu a ostavlja boje. 
Barplot sa standard Error - napredni dio
Barplot sa standard error su česti kada imamo uzorak. da bismo dobili takav graf prvo moramao da razmjemo sta je standard erorr. Uzimamo prosjek varijable određenog sampla, nakon toga izračunamo njegovu standardnu devijaciju sa sd(), broj observacija lenght(), i onda dobijemo standard error. Radicemo na podacima jednog uzorka pacijenata, istrazivanje provedeno u bolnici. Bazu cemo nazvati kardio.
kardio <- read.csv2("~/alma master/amra0705.csv")
head(kardio[,1:6]) #provjeravamo varijable i brojeve## TG ID pol starost bolest kormobiditeti
## 1 1 Pacijent1 0 39 A.psoriatica 2
## 2 2 Pacijent2 0 67 A.psoriatica 2
## 3 3 Pacijent3 1 66 A.psoriatica 2
## 4 4 Pacijent4 0 33 A.psoriatica 2
## 5 5 Pacijent5 0 25 A.psoriatica 2
## 6 6 Pacijent6 1 28 Arthr.reactiva 1tail (kardio[,1:6]) #provjeramamo kraj tabele i vidimo da nije potpuna## TG ID pol starost bolest kormobiditeti
## 113 NA NA NA
## 114 NA kormobiditeti NA NA
## 115 NA Ne 1 NA
## 116 NA metabolièka 2 NA
## 117 NA pulmonarna 3 NA
## 118 NA nervna 4 NAkardio <- kardio[1:110,] #uzimamo samo sto je baza
tapply (kardio$SE, kardio$pol, mean, na.rm=TRUE)## 0 1
## 45.87753 43.04762Sada izracunamo varijable koje nam trebaju: mean, length, sd, se i napravimo bazu sa ovim podacima. mi cemo racunati za sedimentaciju po spolu:
NSE = length(kardio$SE)
meanSE = aggregate(SE~pol, data = kardio, mean) #class (meanSE) je vektor
sdSE = aggregate(SE~pol, data = kardio, sd)
seSE <- sdSE/sqrt(NSE)Sada izvucemo vektore za graf tj. intersuje nas mean i s.error sedimentacije po spolu.
seSE<-unlist(seSE$SE)
final_data<-cbind(meanSE,seSE)
final_data$pol <- as.character(final_data$pol) #class() je pokazivao da je $pol integer, a on je u stvari character, jer graf je bio drugaciji, legenda je bila "valer plave boje"
ggplot(final_data, aes (x=pol, y=SE, fill=pol)) + geom_bar(position=position_dodge(), stat="identity")+geom_errorbar(aes(ymin=SE-seSE, ymax=SE+seSE), width=.1) 
Ostali grafovi
Sada npr. mozete dodati linije a kasnije i boje grafikonu. Ipak nije primjereno da stavljamo linije na ovim varijablama. Linije su primjerene kada posmatramo godine (trendove), ali radi primjera i mogucnosti sada stavljamo:
Linijski
plot1 = ggplot(WHO, aes(x= Region, y=FertilityRate))
Linijski= plot1 + geom_line() + ggtitle("Linijski")Barplot jednostavni
Bargraf = plot1 + geom_bar(stat="identity") + ggtitle ("Bargraf")Distribucija
plot2 = ggplot(WHO, aes(x=FertilityRate))
distribucija= plot2 + geom_density() + ggtitle("Distribucija")Histogram
plot2 = ggplot(WHO, aes(x=FertilityRate))
histogram = plot2 + geom_histogram() + ggtitle ("Histogram, banwidth=30")grid.arrange (Linijski, Bargraf, distribucija, histogram,ncol=2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Poredjenje distribucija
plot3 = ggplot(WHO, aes(x=FertilityRate, fill=Region))
kompar.poredjen = plot3 + geom_density() + ggtitle("Distribucija komparativno")
kompar.poredj_viz = plot3 + geom_density(alpha=0.25) + ggtitle("Preglednije poredjenje uz `geom_density(alpha=0.25)`")grid.arrange (kompar.poredjen,kompar.poredj_viz, ncol=2)
legend.postion i legend.justification pozicioniraju legendu.
Tako npr. vidimo da je R automatski postavio legendu van grafa sto u ovom uporednom slucaju dovodi do suzavanja grafova i njihove nepreglednosti. Radi toga dodadajemo theme () sa odredjivanjem pozicije legende.
kompar.poredjen = plot3 + geom_density() + ggtitle("Distribucija komparativno")+ theme(legend.position=c(1,1), legend.justification =c(1,1))
kompar.poredj_viz = plot3 + geom_density(alpha=0.25) + ggtitle("Preglednije poredjenje uz alpha=0.25")+ theme(legend.position=c(1,1), legend.justification =c(1,1))grid.arrange (kompar.poredjen,kompar.poredj_viz, ncol=2)
Word Clouds
MAPS
library (maps)
library (ggmap)
sarajevo = get_map(location = "sarajevo", zoom=15)## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=sarajevo&zoom=15&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=sarajevo&sensor=falseggmap(sarajevo)
Mapa svijeta
Ako zelimo da neki indikator predstavimo na mapi svijeta, ili mapi pojedinog kontitenta onda nam za to trebal world_map.
Prvo formirajmo world_map
world_map = map_data("world")
head (world_map,3)## long lat group order region subregion
## 1 -69.89912 12.45200 1 1 Aruba <NA>
## 2 -69.89571 12.42300 1 2 Aruba <NA>
## 3 -69.94219 12.43853 1 3 Aruba <NA>head (WHO,3)## Country Region Population Under15 Over60
## 1 Afghanistan Eastern Mediterranean 29825 47.42 3.82
## 2 Albania Europe 3162 21.33 14.93
## 3 Algeria Africa 38482 27.42 7.17
## FertilityRate LifeExpectancy ChildMortality CellularSubscribers
## 1 5.40 60 98.5 54.26
## 2 1.75 74 16.7 96.39
## 3 2.83 73 20.0 98.99
## LiteracyRate GNI PrimarySchoolEnrollmentMale
## 1 NA 1140 NA
## 2 NA 8820 NA
## 3 NA 8310 98.2
## PrimarySchoolEnrollmentFemale
## 1 NA
## 2 NA
## 3 96.4Sada cemo da mergamo ove dvije baze. Vidimo da su $region iz world_map i $Country iz WHO iste varijable. Novu bazu cemo nazvati world_map1
world_map1 = merge(world_map, WHO, by.x = "region", by.y = "Country")
head(world_map1,3)## region long lat group order subregion
## 1 Afghanistan 68.21397 31.80737 2 122 <NA>
## 2 Afghanistan 66.23125 29.86572 2 156 <NA>
## 3 Afghanistan 69.18691 31.83809 2 111 <NA>
## Region Population Under15 Over60 FertilityRate
## 1 Eastern Mediterranean 29825 47.42 3.82 5.4
## 2 Eastern Mediterranean 29825 47.42 3.82 5.4
## 3 Eastern Mediterranean 29825 47.42 3.82 5.4
## LifeExpectancy ChildMortality CellularSubscribers LiteracyRate GNI
## 1 60 98.5 54.26 NA 1140
## 2 60 98.5 54.26 NA 1140
## 3 60 98.5 54.26 NA 1140
## PrimarySchoolEnrollmentMale PrimarySchoolEnrollmentFemale
## 1 NA NA
## 2 NA NA
## 3 NA NASada nacrtajmo mapu Under15. :
ggplot(world_map, aes(x=long, y=lat, group=group, fill=Under15))+geom_polygon(fill="white", color="black") + coord_map("mercator")
Ova mapa nam cudno izgleda. Najverovatnije da nam nije pravilo poslozen. Napravimo order.
world_map1=world_map1[order(world_map1$group,world_map1$order),]
ggplot(world_map1, aes(x=long, y=lat, group=group, fill=Under15))+geom_polygon() + coord_map("mercator")
Ni ova mapa nije zavrsena, tj. vidimo da je poprilicno cudna. To je zbog toga sto nam imena nekih drzava nisu jednaka u map_data i u WHO database. Tako npr. u mapa_data Iran je nazvan Iran a u WHO bazi podataka nazvan je "Iran (Islamic Republic of)". Posto to sada znamo zamjenimo u WHO nazovimo Iran i opet merdzajmo..
Ponovimo merdjanje, order i plot.
Ako imamo situaciju da imamo dane u grafikonu. R nam reda dane po abecedi. Da bismo to ispraviliprvo nam dani u sedmici moraju biti formatirani kao datum. nakon toga kucamo:
baza$dani = factor (baza$dani, order=TRUE, level=c("Sunday", "Monday", ..."Saturday")).
Graficki prikaz intervala povjerenja ANOVE - NAPREDNI DIO
Stavicemo odnos izmedju bolesti i trajanja bolesti. Prvo provjerimo nivoe facotra bolest.
levels(kardio$bolest)## [1] "" "A.psoriatica" "Arthr.reactiva" "M.Bechterew"
## [5] "MBVT" "RA-" "Ra+" "RA+"
## [9] "SSP"Vidimo da imamo pogresno napisano Ra+ umjesto RA+ i da imamo jedan prazan level.
To izmjenimo:
kardio$bolest[kardio$bolest=="Ra+"] <- "RA+"
table (kardio$bolest)##
## A.psoriatica Arthr.reactiva M.Bechterew MBVT
## 0 5 7 3 16
## RA- Ra+ RA+ SSP
## 16 0 56 7Ispustimo prazne nivoe:
kardio$bolest_relevel <- droplevels(kardio$bolest, droplevels.factor = c("Ra+", " "))Uradimo anovu:
anova=aov(kardio$trajanje~kardio$bolest_relevel)Predstavimo u tabelarnom pregledu:
library(xtable)
library(knitr)
kable (xtable (anova))| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| kardio$bolest_relevel | 6 | 431.6107 | 71.93512 | 1.293746 | 0.2671077 |
| Residuals | 99 | 5504.6157 | 55.60218 | NA | NA |
Iako je p-value veci od 0.05 predstavicemo anovu i na grafikonu:
par (mar=c(4.1,8.1,2.1,2.1)) #ovo su marigine: c(dno, lijeva, vrh, desno)
plot(TukeyHSD(anova),las=2,cex.axis=0.7)
Interpretacija. Vidimo da je vecina CI dosta vlika posebno kod M.Bechterew i A.psoriatica. To je radi malog uzorka. Osnovno sto nam ovaj graf pokazuje je da svaki interval povjerenja sijece razliku u mean kad je mean 0. Dakle radi intervala povjerenja mozemo da kazema da niti izmedju bilo koje dvije bolesti ne postoji razlika u trajanju bolesti.
Graficki prikaz Heatmap
Npr. ako zelimo da vidimo koja nam se marka automobila u kojoj opstini najvise puta registrovala preglednije nam je to uraditi sa heatmap nego tabelarno jer ima dosta podataka.
Uzmimo nas primjer registracije automobila.
registracija <- read.csv("~/private/Academy/Week3/prvi_put_registrovana_vozila_12_2015.csv", sep=";")
heat <- as.data.frame(table(registracija$MJESTO.REGISTRACIJE.dec.2015, registracija$MARKA.VOZILA))
heat <- heat[139:24797,]heat_10 <- heat[heat$Freq>=10,] #uzmimo samo one opstine u kojima je pojedina marka registrovana bar 10 puta
ggplot (heat_10, aes(x=Var2, y=Var1)) + geom_tile (aes(fill=Freq)) + theme (axis.text.x = element_text(angle = 45,hjust=1,vjust=1, size = 5)) + theme (axis.text.y = element_text(size = 5)) + scale_fill_gradient (low="white", high = "red")