Friday, March 1, 2013

Overlapping Histogram in R

While preparing a class exercise involving the use of overlaying of histogram, I searched Google on possible article or discussion on the said topic. Luckily, I found a blog where the author demonstrated an R function to create an overlapping histogram. However, a comment from a guy also showed the same output using transparency. Below were the sample codes that can be used to generate overlapping histogram in R as based on the blog and the viewers comment.

                                                                                                                                     



                                                                                                                        
Here are the codes:

#Random numbers
h2<-rnorm(1000,4)
h1<-rnorm(1000,6)

# Histogram Grey Color
hist(h1, col=rgb(0.1,0.1,0.1,0.5),xlim=c(0,10), ylim=c(0,200), main="Overlapping Histogram")
hist(h2, col=rgb(0.8,0.8,0.8,0.5), add=T)
box()

# Histogram Colored (blue and red)
hist(h1, col=rgb(1,0,0,0.5),xlim=c(0,10), ylim=c(0,200), main="Overlapping Histogram", xlab="Variable")
hist(h2, col=rgb(0,0,1,0.5), add=T)
box()

Friday, October 26, 2012

NSCB Sexy Stats Version 2

This was a revised version of my previous post about the NSCB article. With the suggestion from Tal Galili, below were the new pie charts and the R codes to produce these plots by directly scrapping the data from the webpage using XML and RColorBrewer pagkage.

Unemployment by Age Group

Unemployment by Gender

Unemployment by Civil Status

Unemployment by Educational Level




##################################################################################

library (XML)
library(RColorBrewer)
 

url<-"http://www.nscb.gov.ph/sexystats/2012/Filipinoversion/SS20121022_joblessness_filver.asp"
unemployment<-readHTMLTable(url, header=T, which=2,stringsAsFactors=F)
agegroup<-unemployment[3:8,c(1,3,5)]
gender<-unemployment[12:13,c(1,3,5)]
civil<-unemployment[17:20,c(1,3,5)]
education<-unemployment[25:31,c(1,3,5)]


#Copy to clipboard
Education    Y2006    Y2009
Elementary     42.5    37.9
High School    47.7    52.2
College    9.7    10


educ<-read.table("clipboard", header=T, sep="\t")

colnames(agegroup)<-c("Age.Group","Y2006","Y2009")
colnames(gender)<-c("Gender","Y2006","Y2009")
colnames(civil)<-c("Civil.Status","Y2006","Y2009")
colnames(education)<-c("Education","Y2006","Y2009")

agegroup$Age.Group[6]<-"65 & Up"
agegroup$Y2006<-as.numeric(agegroup$Y2006)
agegroup$Y2009<-as.numeric(agegroup$Y2009)

gender$Gender[2]<-"Female"
gender$Gender[1]<-"Male"
gender$Y2006<-as.numeric(gender$Y2006)
gender$Y2009<-as.numeric(gender$Y2009)

civil$Y2006<-as.numeric(civil$Y2006)
civil$Y2009<-as.numeric(civil$Y2009)
cs<-c("Single","Married","Widowed", "Divorced")

win.graph(w=14.3,h=7)
par(mfrow=c(1,2), oma=c(1,0,1,1) , mar=c(1,1,0,1))

#Chart 1
pie(agegroup$Y2006,label=agegroup$Age.Group, col=brewer.pal(6,"Set1"), border="white")
par(new=TRUE) 
pie(c(1), labels=NA, border='white', radius=0.4)
text(0,0,labels="Percent\nUnemployment\nby Age Group\nYear 2006", cex=1.5, font=2)

pie(agegroup$Y2009,label=agegroup$Age.Group, col=brewer.pal(6,"Set1"), border="white")
par(new=TRUE) 
pie(c(1), labels=NA, border='white', radius=0.4)
text(0,0,labels="Percent\nUnemployment\nby Age Group\nYear 2009", cex=1.5, font=2)
text(0.5,-1, "Data Source: NSCB\nCreated by: ARSalvacion", adj=c(0,0), cex=0.7)

#Chart 2
pie(gender$Y2006,label=gender$Gender, col=brewer.pal(2,"Set1"), border="white", cex=1.5)
par(new=TRUE) 
pie(c(1), labels=NA, border='white', radius=0.4)
text(0,0,labels="Percent\nUnemployment\nby Gender\nYear 2006", cex=1.5, font=2)

pie(gender$Y2009,label=gender$Gender, col=brewer.pal(2,"Set1"), border="white", cex=1.5)
par(new=TRUE) 
pie(c(1), labels=NA, border='white', radius=0.4)
text(0,0,labels="Percent\nUnemployment\nby Gender\nYear 2009", cex=1.5, font=2)
text(0.5,-1, "Data Source: NSCB\nCreated by: ARSalvacion", adj=c(0,0), cex=0.7)

#Chart 3
pie(civil$Y2006,label=cs, col=brewer.pal(4,"Dark2"), border="white")
par(new=TRUE) 
pie(c(1), labels=NA, border='white', radius=0.4)
text(0,0,labels="Percent\nUnemployment\nby Civil  Status\nYear 2006", cex=1.5, font=2)

pie(civil$Y2009,label=cs, col=brewer.pal(4,"Dark2"), border="white")
par(new=TRUE) 
pie(c(1), labels=NA, border='white', radius=0.4)
text(0,0,labels="Percent\nUnemployment\nby Civil  Status\nYear 2009", cex=1.5, font=2)
text(0.5,-1, "Data Source: NSCB\nCreated by: ARSalvacion", adj=c(0,0), cex=0.7)

#Chart 4
pie(educ$Y2006,label=educ$Education, col=brewer.pal(3,"Dark2"), border="white", cex=1.5)
par(new=TRUE) 
pie(c(1), labels=NA, border='white', radius=0.4)
text(0,0,labels="Percent\nUnemployment by\nEducational Level\nYear 2006", cex=1.5, font=2)

pie(educ$Y2009,label=educ$Education, col=brewer.pal(3,"Dark2"), border="white", cex=1.5)
par(new=TRUE) 
pie(c(1), labels=NA, border='white', radius=0.4)
text(0,0,labels="Percent\nUnemployment by\nEducational Level\nYear 2009", cex=1.5, font=2)
text(0.5,-1, "Data Source: NSCB\nCreated by: ARSalvacion", adj=c(0,0), cex=0.7)


##################################################################################


Thursday, October 25, 2012

NSCB Sexy Statistics (Unemployment)


Recently, my friend posted on her Facebook account about the article published by the National Statistical Coordination Board (NSCB) about poverty and unemployment in the country.  Looking at the report I saw a lot of tables, so I thought why not use R to make some graphs out of these tables (with some inspiration from the mazmascience).


Percent Unemployed Per Age Group




Percent Unemployed by Sex Group



 Percent Unemployed by Civil Status





Wednesday, August 8, 2012

Rainfall Amount Flooding Quezon City Philippines

The rainfall received by Quezon City, Philippines was almost double of what the city normally receive for the entire month of August, causing flooding and land slide to various villages in the area.




August 6-7 Rainfall on Metro Manila

Majority of Metro Manila is affected by floods. Looking at the hourly data from PAGASA weather  station located at Bicutan, Taguig, the graph below will gave the viewers of this blog on the rainfall situation in Manila from August 6 (12:00am)-7(11:15 pm), 2012.

Tuesday, August 7, 2012

Provincial Monthly Rainfall of the Philippines from WORLDCLIM

Preparing for a future conference on climate change, I downloaded and extracted average monthly rainfall in the Philippines from worldclim.org. Using maptools, raster, and animation package in R, I produced an animation of average monthly rainfall of the country.


Monday, August 6, 2012

Provincial Map using GADM

This blog demonstrates how to produce political/provincial boundary map (below) using R maptools and raster packages.



## Load required packages
library(maptools)
library(raster)

 ## Download data from gadm.org
adm <- getData('GADM', country='PHL', level=2)
mar<-(adm[adm$NAME_1=="Marinduque",])plot(mar, bg="dodgerblue", axes=T)

##Plot downloaded data
plot(mar, lwd=10, border="skyblue", add=T)
plot(mar,col="green4", add=T)
grid()
box()
invisible(text(getSpPPolygonsLabptSlots(mar), labels=as.character(mar$NAME_2), cex=1.1, col="white", font=2))
mtext(side=3, line=1, "Provincial Map of Marinduque", cex=2)
mtext(side=1, "Longitude", line=2.5, cex=1.1)
mtext(side=2, "Latitude", line=2.5, cex=1.1)
text(122.08,13.22, "Projection: Geographic\nCoordinate System: WGS 1984\nData Source: GADM.org\nCreated by: ARSsalvacion", adj=c(0,0), cex=0.7, col="grey20")