Vincent Schubert Malbas ®: May 2013

tutorials in python v R

UPDATED:
Here's the BioPython Tutorial and Cookbook --> BioPython
It's made by Jeff Chang and colleagues, and it seems to be a very useful reference for bioinformatics.

Actually, I have now found the Documentation website for Python 2.6 --> Python Docs
It probably has almost all basic functions defined in there. Of course, there's always StackOverflow...

------- ------- ------- ------- ------- ------- ------- ------- ------- ------- -------

For future reference, I'm keeping this Python tutorial URL --> Tutorials Point

Obviously I'm new to this. I'm keeping tabs on useful websites for Python learning.
I am not aware yet of a function that easily calls on the documentation from the Python GUI.
This makes my learning very arduous, a "trial-and-error" kind of thing.
#imstillhavingfunatleast

In contrast, R has "?" and "??" functions, which will search the appropriate documentation for that queried function. For example, I can type ?spplot, which will look up the term 'spplot' from the local library, i.e., files in the computer, or ??spplot, which will look up the term from the internet, i.e., from the CRAN website. #obviouslyimafan

If you happen to know how to search for Python functions more easily, do give me a heads up!
Currently, looking for the right function to write is a P.I.A. I'm not complaining though; I thank God that for some reason I'm having fun with this particular challenge.

senate election results

click on the image to enlarge, or click this link to view this infographic on visual.ly

These data are sourced from the 16th Official Canvass Report from the Comelec Rappler Mirror Server. Retrieved May 23, 2013 4:30am. http://bit.ly/12AmEyO

you may also download the high-res dataviz from dropbox --> http://bit.ly/1a9KBPU

The following source code for making the heatmaps is working and functional in R v. 3.0.0

#this code was tested in
# R version 3.0.0 (2013-04-03) -- "Masked Marvel"
# Copyright (C) 2013 The R Foundation for Statistical Computing
# Platform: i386-w64-mingw32/i386 (32-bit)

getwd()
# this is your working directory

# setwd("")
# you can set your personal working directory here

getwd()

#  ------------------------------------------------------------------------------------------------
# # LOAD VOTES DATASETS from COMELEC's 16th canvass report (OFFICIAL TALLY)
# Accessed from http://www.rappler.com/nation/politics/elections-2013/features/rich-media/29126-official-tally-votes-2013-senatorial-race
# Retrieved May 23, 2013 4:30am


data1 <- "http://dl.dropboxusercontent.com/u/7911075/election%202013%20csv/namebyregion.csv"
download.file(data1, destfile="./datasets/namebyregion.csv")

data2 <- "http://dl.dropboxusercontent.com/u/7911075/election%202013%20csv/namebyprovince.csv"
download.file(data2, destfile="./datasets/namebyprovince.csv")

data3 <- "http://dl.dropboxusercontent.com/u/7911075/election%202013%20csv/provincetotalbyname.csv"
download.file(data3, destfile="./datasets/provincetotalbyname.csv")

data4 <- "http://dl.dropboxusercontent.com/u/7911075/election%202013%20csv/regiontotalbyname.csv"
download.file(data4, destfile="./datasets/regiontotalbyname.csv")

list.files(".")
dateDownloaded <-date()
dateDownloaded

byregion <- read.csv("./datasets/namebyregion.csv", header=TRUE)
str(byregion)

byprovince <- read.csv("./datasets/namebyprovince.csv", header=TRUE)
str(byprovince)

tprovince <- read.csv("./datasets/provincetotalbyname.csv", header=TRUE)
str(tprovince)

tregion <- read.csv("./datasets/regiontotalbyname.csv", header=TRUE)
str(tregion)

#  ------------------------------------------------------------------------------------------------
# # LOAD COLOR PALETTES

library(RColorBrewer)

cols4v1 <- brewer.pal(5, "YlGnBu")
pal4v1 <- colorRampPalette(cols4v1)  # <- dark blue cyan white gradient

cols4v2 <- brewer.pal(5, "YlOrBr")
pal4v2 <- colorRampPalette(cols4v2)  # <- brown orange white gradient

cols4v3 <- brewer.pal(5, "PuRd")
pal4v3 <- colorRampPalette(cols4v3)  # <- purple pink white gradient

cols4v4 <- brewer.pal(5, "RdPu")
pal4v4 <- colorRampPalette(cols4v4)  # <- violet pink white gradient

cols8 <- brewer.pal(8, "Set3")
pal8 <- colorRampPalette(cols8)

cols12 <- brewer.pal(12, "Set3")
pal12 <- colorRampPalette(cols12)

#  ------------------------------------------------------------------------------------------------
# # DRAW HEATMAPS FOR REGION

row.names(byregion) <- byregion$X
byregion <- byregion[,2:19]
byregion_matrix <- data.matrix(byregion)

regions_heatmap1 <- heatmap(byregion_matrix, Rowv=NA, Colv=NA, col = cm.colors(256), scale="column", margins=c(5,10)) # <-- pink, blue
regions_heatmap2 <- heatmap(byregion_matrix, Rowv=NA, Colv=NA, col = heat.colors(256), scale="column", margins=c(5,10)) # <-- too much red
regions_heatmap3 <- heatmap(byregion_matrix, Rowv=NA, Colv=NA, col = pal4v2(5), scale="column", margins=c(5,10)) # <-- yellow brown gradient
regions_heatmap4 <- heatmap(byregion_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
regions_heatmap5 <- heatmap(byregion_matrix, Rowv=NA, Colv=NA, col = pal12(5), scale="column", margins=c(5,10)) # <- messy looking

regions_heatmap6 <- heatmap(byregion_matrix, Rowv=NA, Colv=NA, col = pal4v4(5), scale="column", margins=c(5,10)) # <-- purle pink gradient
# dev.copy2pdf(file="./output/votesperregion.pdf", height = 8, width = 11)

#  ------------------------------------------------------------------------------------------------
# # DRAW HEATMAPS PER PROVINCE

str(byprovince)
names(byprovince)

#  REGION 1 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region1
provincesregion1 <- byprovince[,33:37]
provincesregion1$X <- byprovince$X
str(provincesregion1)
row.names(provincesregion1) <- provincesregion1$X
provincesregion1 <- provincesregion1[,1:5]
provincesregion1_matrix <- data.matrix(provincesregion1)

#draw heatmap for provinces in region1
region1_heatmap <- heatmap(provincesregion1_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
dev.copy2pdf(file="./output/votesregion1.pdf", height = 8, width = 11)

#  REGION 2 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region2
provincesregion2 <- byprovince[,38:43]
provincesregion2$X <- byprovince$X
str(provincesregion2)
row.names(provincesregion2) <- provincesregion2$X
provincesregion2 <- provincesregion2[,1:6]
provincesregion2_matrix <- data.matrix(provincesregion2)

#draw heatmap for provinces in region2
region2_heatmap <- heatmap(provincesregion2_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion2.pdf", height = 8, width = 11)

#  REGION 3 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region3
provincesregion3 <- byprovince[,44:51]
provincesregion3$X <- byprovince$X
str(provincesregion3)
row.names(provincesregion3) <- provincesregion3$X
provincesregion3 <- provincesregion3[,1:8]
provincesregion3_matrix <- data.matrix(provincesregion3)

#draw heatmap for provinces in region3
region3_heatmap <- heatmap(provincesregion3_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion3.pdf", height = 8, width = 11)

#  REGION 4A ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region4a
provincesregion4a <- byprovince[,52:57]
provincesregion4a$X <- byprovince$X
str(provincesregion4a)
row.names(provincesregion4a) <- provincesregion4a$X

provincesregion4a <- provincesregion4a[,1:6]
provincesregion4a_matrix <- data.matrix(provincesregion4a)

#draw heatmap for provinces in region4a
region4a_heatmap <- heatmap(provincesregion4a_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion4a.pdf", height = 8, width = 11)

#  REGION4B ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region4b
provincesregion4b <- byprovince[,58:63]
provincesregion4b$X <- byprovince$X
str(provincesregion4b)
row.names(provincesregion4b) <- provincesregion4b$X

provincesregion4b <- provincesregion4b[,1:6]
provincesregion4b_matrix <- data.matrix(provincesregion4b)

#draw heatmap for provinces in region4b
region4b_heatmap <- heatmap(provincesregion4b_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion4b.pdf", height = 8, width = 11)

#  REGION5 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region5
provincesregion5 <- byprovince[,64:70]
provincesregion5$X <- byprovince$X
str(provincesregion5)
row.names(provincesregion5) <- provincesregion5$X

provincesregion5 <- provincesregion5[,1:7]
provincesregion5_matrix <- data.matrix(provincesregion5)

#draw heatmap for provinces in region5
region5_heatmap <- heatmap(provincesregion5_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion5.pdf", height = 8, width = 11)

#  REGION6 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region6
provincesregion6 <- byprovince[,71:79]
provincesregion6$X <- byprovince$X
str(provincesregion6)
row.names(provincesregion6) <- provincesregion6$X

provincesregion6 <- provincesregion6[,1:9]
provincesregion6_matrix <- data.matrix(provincesregion6)

#draw heatmap for provinces in region6
region6_heatmap <- heatmap(provincesregion6_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion6.pdf", height = 8, width = 11)

#  REGION7 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region7
provincesregion7 <- byprovince[,80:86]
provincesregion7$X <- byprovince$X
str(provincesregion7)
row.names(provincesregion7) <- provincesregion7$X

provincesregion7 <- provincesregion7[,1:7]
provincesregion7_matrix <- data.matrix(provincesregion7)

#draw heatmap for provinces in region7
region7_heatmap <- heatmap(provincesregion7_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion7.pdf", height = 8, width = 11)

#  REGION8 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region8
provincesregion8 <- byprovince[,87:93]
provincesregion8$X <- byprovince$X
str(provincesregion8)
row.names(provincesregion8) <- provincesregion8$X

provincesregion8 <- provincesregion8[,1:7]
provincesregion8_matrix <- data.matrix(provincesregion8)

#draw heatmap for provinces in region8
region8_heatmap <- heatmap(provincesregion8_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion8.pdf", height = 8, width = 11)

#  REGION9 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region9
provincesregion9 <- byprovince[,94:98]
provincesregion9$X <- byprovince$X
str(provincesregion9)
row.names(provincesregion9) <- provincesregion9$X

provincesregion9 <- provincesregion9[,1:5]
provincesregion9_matrix <- data.matrix(provincesregion9)

#draw heatmap for provinces in region9
region9_heatmap <- heatmap(provincesregion9_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion9.pdf", height = 8, width = 11)

#  REGION10 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region10
provincesregion10 <- byprovince[,99:106]
provincesregion10$X <- byprovince$X
str(provincesregion10)
row.names(provincesregion10) <- provincesregion10$X

provincesregion10 <- provincesregion10[,1:7]
provincesregion10_matrix <- data.matrix(provincesregion10)

#draw heatmap for provinces in region10
region10_heatmap <- heatmap(provincesregion10_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion10.pdf", height = 8, width = 11)

#  REGION11 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region11
provincesregion11 <- byprovince[,107:112]
provincesregion11$X <- byprovince$X
str(provincesregion11)
row.names(provincesregion11) <- provincesregion11$X

provincesregion11 <- provincesregion11[,1:6]
provincesregion11_matrix <- data.matrix(provincesregion11)

#draw heatmap for provinces in region11
region11_heatmap <- heatmap(provincesregion11_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion11.pdf", height = 8, width = 11)

#  REGION12 ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region12
provincesregion12 <- byprovince[,113:117]
provincesregion12$X <- byprovince$X
str(provincesregion12)
row.names(provincesregion12) <- provincesregion12$X

provincesregion12 <- provincesregion12[,1:5]
provincesregion12_matrix <- data.matrix(provincesregion12)

#draw heatmap for provinces in region12
region12_heatmap <- heatmap(provincesregion12_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion12.pdf", height = 8, width = 11)

#  REGION13 CARAGA ------------------------------------------------------------------------------------------------
# for example, we take the provinces in region13
provincesregion13 <- byprovince[,10:15]
provincesregion13$X <- byprovince$X
str(provincesregion13)
row.names(provincesregion13) <- provincesregion13$X

provincesregion13 <- provincesregion13[,1:6]
provincesregion13_matrix <- data.matrix(provincesregion13)

#draw heatmap for provinces in region13
region13_heatmap <- heatmap(provincesregion13_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesregion13.pdf", height = 8, width = 11)

#  NCR  ------------------------------------------------------------------------------------------------
# for example, we take the provinces in ncr
provincesncr <- byprovince[,16:32]
provincesncr$X <- byprovince$X
str(provincesncr)
row.names(provincesncr) <- provincesncr$X

provincesncr <- provincesncr[,1:17]
provincesncr_matrix <- data.matrix(provincesncr)

#draw heatmap for provinces in ncr
ncr_heatmap <- heatmap(provincesncr_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesncr.pdf", height = 8, width = 11)

#  CAR  ------------------------------------------------------------------------------------------------
# for example, we take the provinces in car
provincescar <- byprovince[,2:9]
provincescar$X <- byprovince$X
str(provincescar)
row.names(provincescar) <- provincescar$X

provincescar <- provincescar[,1:8]
provincescar_matrix <- data.matrix(provincescar)

#draw heatmap for provinces in car
car_heatmap <- heatmap(provincescar_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votescar.pdf", height = 8, width = 11)

#  ARMM  ------------------------------------------------------------------------------------------------
# for example, we take the provinces in armm
provincesarmm <- byprovince[,118:123]
provincesarmm$X <- byprovince$X
str(provincesarmm)

row.names(provincesarmm) <- provincesarmm$X
provincesarmm <- provincesarmm[,1:6]
provincesarmm_matrix <- data.matrix(provincesarmm)

#draw heatmap for provinces in armm
armm_heatmap <- heatmap(provincesarmm_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesarmm.pdf", height = 8, width = 11)

#  OAV OVERSEAS ABSENTEE VOTING COUNTRIES  ------------------------------------------------------------------------------------------------
# for example, we take the provinces in oav
provincesoav <- byprovince[,124:146]
provincesoav$X <- byprovince$X
str(provincesoav)

row.names(provincesoav) <- provincesoav$X
provincesoav <- provincesoav[,1:23]
provincesoav_matrix <- data.matrix(provincesoav)

#draw heatmap for provinces in oav
oav_heatmap <- heatmap(provincesoav_matrix, Rowv=NA, Colv=NA, col = pal4v1(5), scale="column", margins=c(5,10)) # <-- yellow blue gradient
# dev.copy2pdf(file="./output/votesoav.pdf", height = 8, width = 11)


# # ---------------------------
# END OF CODE

dataviz on antenatal care

in a previous post, i showed you how to make choropleth maps, particularly both world and Philippine maps, using R (see previous post). In this case we used the data we collected from the websites of the World Health Organization and the Philippine Department of Health.

improving antenatal care infographic

infographic newborn screening

what is newborn screening infographic

ironman3 v startrek2

this news report was brought to you by xkcd

i hate doing reviews; one can argue that films don't need them.
reviews exist because of movies, not vice versa.

don't get me wrong; star trek is pretty solid work from the reliable jj abrams.
it had a tight script, crazy visuals, character buildup, and some bromance on the side.

that being said, i agonized for 2 hours inside the cinema.
i was hoping to find salvation for startrek2 after much anticipation.
i wanted that one break-out moment that can topple iron man 3.
there was none.
to be fair, my pre-expectations were quite high for star trek, and a bit low for iron man.
there was a very high chance that it will fail my expectations.

as a trekkie, it pains me to say that ironman3 has edged out startrek2 in the emotions department.
where did my good ol' jj abrams go?
ironman3 felt, surprisingly, more human.

you can do the review on your own.
think about 3 themes:
infallibility: which felt more vulnerable?
sense of threat: what was the main goal of the antagonist? was it killing or was it just escaping?
consequences: do you really think that the protagonist will die?

having said that, startrek is best seen on 3D and with fellow trekkies.
live long and prosper, b*tches.

choropleth maps in R

below is the source code for making the choropleths for the infographic about newborns. you can run the following code in R. By the way, I'm using version 3.0.0; some library packages may have a different syntax and output in older versions of R.

getwd()
# this is your working directory on your computer

# DOWNLOAD THE FOLLOWING DATA FROM THE INTERNET
provinces <- "http://dl.dropboxusercontent.com/u/7911075/csv%20newborn/provinces.csv"
download.file(provinces, destfile="./provinces.csv")

phattended <- "http://dl.dropboxusercontent.com/u/7911075/csv%20newborn/dataPH2003-2008attended.csv"
download.file(phattended, destfile="./dataPH2003-2008attended.csv")

phbirthshome <- "http://dl.dropboxusercontent.com/u/7911075/csv%20newborn/dataPH2003-2008facility.csv"
download.file(phbirthshome, destfile="./dataPH2003-2008facility.csv")

countries <- "http://dl.dropboxusercontent.com/u/7911075/csv%20newborn/datacountries2005-2011attended.csv"
download.file(countries, destfile="./datacountries2005-2011attended.csv")

# this is a large file from the gadm website, might take a while to finish downloading
phil1 <- "http://gadm.org/data/rda/PHL_adm1.RData"
download.file(phil1, destfile="./phil1.RData")

list.files(".")
dateDownloaded <-date()
dateDownloaded

#  ------------------------------------------------------------------------------------------------
# # LOAD GEO MAPS

library(sp)

con <- "./phil1.RData"
print(load(con))
close(con)
str(gadm, max.level=2)
# phil1.RData has 82 observations of 16 variables

names(gadm)
#  [1] "ID_0"       "ISO"        "NAME_0"     "ID_1"       "NAME_1"    
# [6] "VARNAME_1"  "NL_NAME_1"  "HASC_1"     "CC_1"       "TYPE_1"    
# [11] "ENGTYPE_1"  "VALIDFR_1"  "VALIDTO_1"  "REMARKS_1"  "Shape_Leng"
# [16] "Shape_Area"

#  ------------------------------------------------------------------------------------------------
# # LOAD PROVINCES DATASET
provinces <- read.csv("./provinces.csv")
str(provinces)

#  ------------------------------------------------------------------------------------------------
# # LOAD DATA BIRTHS ATTENDED BY SKILLED PERSONNEL
phattended <- read.csv("./dataPH2003-2008attended.csv", header=TRUE)
str(phattended)

#calculate percent births with skilled attendants, i.e., doctors, nurses, and midwives only
phattended$skilled2008percent <- (phattended$X2008doctorspercent + phattended$X2008nursespercent + phattended$X2008midwivespercent)
phattended$skilled2007percent <- (phattended$X2007doctorspercent + phattended$X2007nursespercent + phattended$X2007midwivespercent)
phattended$skilled2006percent <- (phattended$X2006doctorspercent + phattended$X2006nursespercent + phattended$X2006midwivespercent)
phattended$skilled2005percent <- (phattended$X2005doctorspercent + phattended$X2005nursespercent + phattended$X2005midwivespercent)
phattended$skilled2004percent <- (phattended$X2004doctorspercent + phattended$X2004nursespercent + phattended$X2004midwivespercent)
phattended$skilled2003percent <- (phattended$X2003doctorspercent + phattended$X2003nursespercent + phattended$X2003midwivespercent)
#example: barplot(phattended$X2008midwivespercent, col=unique(phattended$REGION), xlab="", ylab="", cex.axis=0.8)

#multiply by 100
phattended$skilled2008percent <- phattended$skilled2008percent*100
phattended$skilled2007percent <- phattended$skilled2007percent*100
phattended$skilled2006percent <- phattended$skilled2006percent*100
phattended$skilled2005percent <- phattended$skilled2005percent*100
phattended$skilled2004percent <- phattended$skilled2004percent*100
phattended$skilled2003percent <- phattended$skilled2003percent*100

#  ------------------------------------------------------------------------------------------------
# # LOAD DATA PHL HOMEBIRTHS, NORMAL SPONTANEOUS DELIVERY
phbirthshome <- read.csv("./dataPH2003-2008facility.csv", header=TRUE)
str(phbirthshome)
names(phbirthshome)

#  ------------------------------------------------------------------------------------------------
# # FIRST MERGE: INNER JOIN FOR DATASETS OF % HOME BIRTHS & % ATTENDED BY SKILLED PERSONNEL
phbirthdata <- merge(phattended, phbirthshome, by="REGION")
str(phbirthdata)
 
#  ------------------------------------------------------------------------------------------------
# # SECOND MERGE: INNER JOIN TO PROVINCE DATASET
provattended1 <- merge(provinces, phbirthdata, by="REGION")
str(provattended1)
# write.csv(provattended1, file="provattended1.csv")

# SORT MERGED DATASET TO ARRANGE ACCDG TO GADM MAP's ORDER
provattended1.sorted <- provattended1[order(provattended1$PROVINCEORDER) , ]
# write.csv(provattended1.sorted, file="provattended1.sorted.csv")

# ------------------------------------------------------------------------------------------------------------
# # NOW WE CAN MERGE THE COMBINED DATASET TO THE GADM MAP

gadm$region <- as.factor(provattended1.sorted$REGION)
gadm$capital <- as.factor(provattended1.sorted$CAPITAL)
gadm$popln2010 <- as.numeric(provattended1.sorted$POPLN2010)
gadm$landarea <- as.numeric(provattended1.sorted$LANDAREA.km2)
gadm$popdense <- as.numeric(provattended1.sorted$POPDENS2010)
gadm$livebirths2008 <- as.numeric(provattended1.sorted$X2008total)
gadm$skilled2008 <- as.numeric(provattended1.sorted$skilled2008percent)
gadm$homebirths2008 <- as.numeric(provattended1.sorted$X2008HOMEpercent)
gadm$hospbirths2008 <- as.numeric(provattended1.sorted$X2008HOSPpercent)

library(Hmisc)
library(RColorBrewer)

cols4v1 <- brewer.pal(5, "YlGnBu")
pal4v1 <- colorRampPalette(cols4v1)  # <- dark blue cyan white gradient

cols4v2 <- brewer.pal(5, "YlOrBr")
pal4v2 <- colorRampPalette(cols4v2)  # <- brown orange white gradient

cols4v3 <- brewer.pal(5, "PuRd")
pal4v3 <- colorRampPalette(cols4v3)  # <- purple pink white gradient

cols8 <- brewer.pal(8, "Set3")
pal8 <- colorRampPalette(cols8)

cols12 <- brewer.pal(12, "Set3")
pal12 <- colorRampPalette(cols12)

str(gadm$region)
str(gadm$capital)
str(gadm$popln2010)
str(gadm$landarea)
str(gadm$popdense)
str(gadm$livebirths2008)
str(gadm$skilled2008)
str(gadm$homebirths2008)

# --------------------------------------------------------------------
#  PARSE REGIONS INTO LOW-AVERAGE-HIGH

gadm$g4area <- cut2(gadm$landarea, g = 4)
table(gadm$g4area)

gadm$g4popln <- cut2(gadm$popln2010, g = 4)
table(gadm$g4popln)

gadm$g4popdense <- cut2(gadm$popdense, g = 4)
table(gadm$g4popdense)

gadm$g4livebirths2008 <-cut2(gadm$livebirths2008, g = 4)
table(gadm$g4livebirths2008)

gadm$g4skilled2008 <-cut2(gadm$skilled2008, g=4)
table(gadm$g4skilled2008)

gadm$g4homebirths2008 <- cut2(gadm$homebirths2008, g=4)
table(gadm$g4homebirths2008)

gadm$g4hospbirths2008 <- cut2(gadm$hospbirths2008, g=4)
table(gadm$g4hospbirths2008)

gadm$g8popln <- cut2(gadm$popln2010, g = 8)
table(gadm$g8popln)

gadm$g8popdense <- cut2(gadm$popdense, g = 8)
table(gadm$g8popdense)

colg4skilled4v1 = pal4v1(length(levels(gadm$g4skilled2008)))
spplot(gadm, "g4skilled2008", col.regions=colg4skilled4v1, main="% Births with Skilled Attendants, 2008 by Region")
# dev.copy2pdf(file="PH percent births attended 2008.pdf", height =11, width = 8)

colg4home4v2 = pal4v2(length(levels(gadm$g4homebirths2008)))
spplot(gadm, "g4homebirths2008", col.regions=colg4home4v2, main="% Births at Home, Normal Spontaneous Delivery, 2008 by Region")
# dev.copy2pdf(file="PH percent births at home 2008.pdf", height =11, width = 8)

colg4hosp4v3 = pal4v3(length(levels(gadm$g4hospbirths2008)))
spplot(gadm, "g4hospbirths2008", col.regions=colg4hosp4v3, main="% Births at Hospital, Normal Spontaneous Delivery, 2008 by Region")
# dev.copy2pdf(file="PH percent births at hospital 2008.pdf", height =11, width = 8)

colpopln4v1 = pal4v1(length(levels(gadm$g4popln))) # <- dark blue cyan white gradient
colpopln4v2 = pal4v2(length(levels(gadm$g4popln))) # <- brown orange white gradient 
colpopln4v3 = pal4v3(length(levels(gadm$g4popln))) # <- purple pink white gradient
# spplot(gadm, "g4popln", col.regions=colpopln4v1, main="Population by Province")

colregions = pal12(length(levels(gadm$region)))
# spplot(gadm, "region", col.regions=colregions, main="Regions in the Philippines")
# dev.copy2pdf(file="PH regions.pdf", height =11, width = 8)

colpopln8 = pal8(length(levels(gadm$g8popln)))
# spplot(gadm, "g8popln", col.regions=colpopln8, main="Population by Province")
# dev.copy2pdf(file="PH population per province.pdf", height =11, width = 8)

colpopdense8 = pal8(length(levels(gadm$g8popdense)))
# spplot(gadm, "g8popdense", col.regions=colpopdense8, main="Population Density by Province")
# dev.copy2pdf(file="PH population density per province.pdf", height =11, width = 8)

collivebirth4v2 = pal4v2(length(levels(gadm$g4livebirths2008)))
# spplot(gadm, "g4livebirths2008", col.regions=collivebirth4v2, main="Live Births 2008")
# dev.copy2pdf(file="live births 2008.pdf", height =11, width = 8)

# ------------------------------------------------------------------------------------------
# # MAKE WORLD MAP

countries <- read.csv("./datacountries2005-2011attended.csv", header=TRUE)
str(countries)
countries$skilled100 <- countries$birthsattendedpercent*100
countries$g8skilled100 <- cut2(countries$skilled100, g = 8)
table(countries$g8skilled100)


library(Hmisc)
library(RColorBrewer)

cols4v1 <- brewer.pal(5, "YlGnBu")
pal4v1 <- colorRampPalette(cols4v1)  # <- dark blue cyan white gradient

cols4v2 <- brewer.pal(5, "YlOrBr")
pal4v2 <- colorRampPalette(cols4v2)  # <- brown orange white gradient

cols4v3 <- brewer.pal(5, "PuRd")
pal4v3 <- colorRampPalette(cols4v3)  # <- purple pink white gradient

cols8 <- brewer.pal(8, "Set3")
pal8 <- colorRampPalette(cols8)

cols8v2 <- brewer.pal(8, "YlOrRd")
pal8v2 <- colorRampPalette(cols8v2)

cols8v3 <- brewer.pal(8, "YlGnBu")
pal8v3 <- colorRampPalette(cols8v3)

cols8v4 <- brewer.pal(8, "RdPu")
pal8v4 <- colorRampPalette(cols8v4)

cols12 <- brewer.pal(12, "Set3")
pal12 <- colorRampPalette(cols12)

library(rworldmap)

mapattended <- joinCountryData2Map(countries, joinCode = "ISO3", nameJoinColumn = "ISO3V10")
str(mapattended, max.level=2)

par(mai=c(0,0,0.2,0),xaxs="i",yaxs="i")
mapParams <- mapCountryData(mapattended, nameColumnToPlot="skilled100", catMethod = "pretty", numCats = 8, colourPalette = cols8v3, mapTitle="", addLegend=FALSE)
do.call( addMapLegend, c(mapParams, legendWidth=0.5, legendMar = 2, legendLabels="all", legendIntervals="data"))
# dev.copy2pdf(file="countries births skilled attendants 2005-2011 v2.pdf", height =8, width = 11)

Navigation