An 19 minute recording of a discussion around plotting in R.
You can watch the recording (18 minutes 43 seconds)
or download (right click and choose the save/download link option) the recording to your computer via the link:
march21.m4v (.m4v format, 88 MB)
and read the parts of the script below (code available when viewing the full article)
if(!(require(dplyr))){
install.packages("dplyr")
require(dplyr)
}
gftext <- "
angle,actual_n,expected_n
0,3735,4171
10,3723,4265
20,4501,5267
30,3512,3947
40,2529,2981
50,2091,2336
60,1599,1763
70,393,469
80,0,0
90,0,0
100,376,469
110,1595,1762
120,1993,2337
130,2542,2980
140,3481,3947
150,4479,5266
160,3724,4265
170,3760,4170
180,3580,3758
190,4565,4215
200,6150,5172
210,4574,3955
220,3449,3049
230,2760,2431
240,2014,1862
250,504,500
260,0,0
270,0,0
280,497,500
290,1993,1861
300,2749,2430
310,3535,3050
320,4630,3954
330,6329,5173
340,4947,4215
350,3970,3758"
gf1 <- read.csv(text=gftext)
#######################
# make a bar graph
barplot(gf1$actual_n,names.arg = gf1$angle, cex.names=0.7,
las=2, cex.axis=0.7,
xlab="10 degree angle of sun (lower bound)",
ylab="Count of Earthquakes")
#####################
# decide to calculate confidence intervals
if(!(require(binom))){
install.packages("binom")
require(binom)
}
ci_brackets <- gf1 %>% mutate(grand_total=sum(actual_n)) %>%
rowwise() %>% mutate(ci_lower = binom.confint(actual_n, grand_total, method=c("wilson"), conf.level = .999999999997440)[1,5] * grand_total,
ci_upper = binom.confint(actual_n, grand_total, method=c("wilson"), conf.level = .999999999997440)[1,6] * grand_total)
# because in base plotting we are drawing points, lines, and shapes in a corodinate system,
# I express the confidence intervals as polygons
gf_ymax <- max(ci_brackets$ci_upper)
gf_ymin <- min(ci_brackets$ci_lower)
# create a blank graph
plot(c(0,360),c(gf_ymin,gf_ymax), type="n", frame.plot=FALSE,
xlab="10 degree angle of sun (lower bound)",
ylab="Count of Earthquakes")
#draw a bunch of polygons
suppress_return <- apply(ci_brackets[,c(1,5,6)], 1, function(x){polygon(x=c(x[1], x[1]+10, x[1]+10, x[1]),
y=c(x[2], x[2], x[3], x[3]), col="#0000FF44",
border=NA)})
# draw a bunch of lines
suppress_return <- apply(ci_brackets[,c(1,2)], 1, function(x){lines(x=c(x[1], x[1]+10),
y=c(x[2], x[2]))})
# add two legends
legend(x=50, y=7000, legend=c("total"), lwd=c(1), col=c("#000000"), bty="n", cex=0.8)
legend(x=50, y=6500, legend=c("7 Sigma Conf. Int."),border=NA , fill=c("#0000FF44"), bty="n", cex=0.8)
#####################
# decide to add expected values
gf_ymax <- max(ci_brackets$ci_upper)
gf_ymin <- min(ci_brackets$ci_lower)
plot(c(0,360),c(gf_ymin,gf_ymax), type="n", frame.plot=FALSE,
xlab="10 degree angle of sun (lower bound)",
ylab="Count of Earthquakes")
suppress_return <- apply(ci_brackets[,c(1,5,6)], 1, function(x){polygon(x=c(x[1], x[1]+10, x[1]+10, x[1]),
y=c(x[2], x[2], x[3], x[3]), col="#0000FF44",
border=NA)})
suppress_return <- apply(ci_brackets[,c(1,2)], 1, function(x){lines(x=c(x[1], x[1]+10),
y=c(x[2], x[2]))})
# add line for expected data
suppress_return <- apply(ci_brackets[,c(1,3)], 1, function(x){lines(x=c(x[1], x[1]+10),
y=c(x[2], x[2]), lty=3)})
legend(x=50, y=7400, legend=c("total","expected"), lty=c(1,3), col=c("#000000"), bty="n", cex=0.8)
legend(x=50, y=6400, legend=c("7 Sigma Conf. Int."),border=NA , fill=c("#0000FF44"), bty="n", cex=0.8)
#####################
# then decide to normalise on the expected values by subtraction in the graph
gf_ymax <- max(ci_brackets$ci_upper - ci_brackets$expected_n)
gf_ymin <- min(ci_brackets$ci_lower - ci_brackets$expected_n)
plot(c(0,360),c(gf_ymin,gf_ymax), type="n", frame.plot=FALSE,
xlab="10 degree angle of sun (lower bound)",
ylab="Difference in count of Earthquakes\nto expected count")
suppress_return <- apply(ci_brackets[,c(1,5,6,3)], 1, function(x){polygon(x=c(x[1], x[1]+10, x[1]+10, x[1]),
y=c(x[2], x[2], x[3], x[3]) - x[4], col="#0000FF44",
border=NA)})
suppress_return <- apply(ci_brackets[,c(1,2,3)], 1, function(x){lines(x=c(x[1], x[1]+10),
y=c(x[2], x[2]) - x[3])})
suppress_return <- apply(ci_brackets[,c(1,3,3)], 1, function(x){lines(x=c(x[1], x[1]+10),
y=c(x[2], x[2]) - x[3], lty=3)})
legend(x=50, y=1400, legend=c("total","expected"), lty=c(1,3), col=c("#000000"), bty="n", cex=0.8)
legend(x=50, y=1000, legend=c("7 Sigma Conf. Int."),border=NA , fill=c("#0000FF44"), bty="n", cex=0.8)
#####################
# But as angles, I can swap to a round corordinate system
# this could painfully be done by hand,
# but I am going to use the plotrix package
if(!(require(plotrix))){
install.packages("plotrix")
require(plotrix)
}
# put the differences in the data
ci_brackets <- ci_brackets %>% mutate(diff_act= actual_n - expected_n,
diff_CI_U = ci_upper - expected_n,
diff_CI_L = ci_lower - expected_n)
ci_brackets$border = 2
# because I need to draw polygons, I need a copy of the data for each edge
ci_brackets2 <- ci_brackets
ci_brackets2$angle <- ci_brackets2$angle + 10
ci_brackets2$border = 1
graphdata <- bind_rows(ci_brackets,ci_brackets2) %>% arrange(angle,border)
# inside and outside limits on the graph
limits=2 * max(abs(c(graphdata$diff_CI_L, graphdata$diff_CI_U)))
#plot the graph
polar.plot(graphdata$diff_CI_U, polar.pos=graphdata$angle, radial.lim=c(-1*limits,limits),
labels = "", main=NULL,lwd=0.5, rp.type="p",
show.grid.labels=FALSE, show.grid=FALSE, mar=c(0,0,0,0),
grid.col="#0000FF44", line.col="#0000FF44", poly.col="#0000FF44")
polar.plot(graphdata$diff_CI_L, polar.pos=graphdata$angle, add=TRUE, radial.lim=c(-1*limits,limits),
line.col="white", lwd=0.5, rp.type="p", poly.col="white")
polar.plot(rep(0,nrow(graphdata)), polar.pos=graphdata$angle, add=TRUE,radial.lim=c(-1*limits,limits),
rp.type="p", lty=4)
polar.plot(rep(-500,nrow(graphdata)), polar.pos=graphdata$angle, add=TRUE,radial.lim=c(-1*limits,limits),
rp.type="p", lty=1, line.col="#00000044")
polar.plot(rep(500,nrow(graphdata)), polar.pos=graphdata$angle, add=TRUE,radial.lim=c(-1*limits,limits),
rp.type="p", lty=3, line.col="#00000044")
#######
lines(c(limits* -1.4,limits* -1.2), c(0,0))
lines(c(limits*1.2,limits*1.4), c(0,0))
text(limits* -1.4,0, label="sunset\n180", cex=0.7)
text(limits* 1.4,0, label="sunrise\n0", cex=0.7)
### Now to have a legend, lets divide up the plot area
layout(matrix(c(1,1,2), nrow=1, ncol=3, byrow = TRUE))
#plot graph 1
polar.plot(graphdata$diff_CI_U, polar.pos=graphdata$angle, radial.lim=c(-1*limits,limits),
labels = "", main=NULL,lwd=0.5, rp.type="p",
show.grid.labels=FALSE, show.grid=FALSE, mar=c(0,0,0,0),
grid.col="#0000FF44", line.col="#0000FF44", poly.col="#0000FF44")
polar.plot(graphdata$diff_CI_L, polar.pos=graphdata$angle, add=TRUE, radial.lim=c(-1*limits,limits),
line.col="white", lwd=0.5, rp.type="p", poly.col="white")
polar.plot(rep(0,nrow(graphdata)), polar.pos=graphdata$angle, add=TRUE,radial.lim=c(-1*limits,limits),
rp.type="p", lty=4)
polar.plot(rep(-500,nrow(graphdata)), polar.pos=graphdata$angle, add=TRUE,radial.lim=c(-1*limits,limits),
rp.type="p", lty=1, line.col="#00000044")
polar.plot(rep(500,nrow(graphdata)), polar.pos=graphdata$angle, add=TRUE,radial.lim=c(-1*limits,limits),
rp.type="p", lty=3, line.col="#00000044")
lines(c(limits* -1.4,limits* -1.2), c(0,0))
lines(c(limits*1.2,limits*1.4), c(0,0))
text(limits* -1.4,0, label="sunset\n180")
text(limits* 1.4,0, label="sunrise\n0")
#plot graph 2
plot(c(0,1),c(0,1), type="n", frame.plot=FALSE, axes=FALSE, xlab="", ylab="")
legend("center", legend=c("Expected Occurance", "500 above expected", "500 below expected", "7 Sigma Conf. Int."),
bty = "n", lty=c(2,3,1,1), lwd=c(1,1,1,8), col=c("#000000","#00000088","#00000088","#0000FF44"))
par(mfrow=c(1,1))
#########################
## picture within picture
# if you are tampering with par, save the current settings
# otherwise you need to restart at the end
oldpar <- par()
gftext <- "
lat,long,loc,col,val
-45.88,170.50,DN,2,26
-43.53,172.63,CH,3,20
-41.29,174.78,WL,4,18
"
gfdata <- read.csv(text=gftext)
if(!(require("maps"))){
install.packages("maps")
require(maps)
}
if(!(require("mapdata"))){
install.packages("mapdata")
require(mapdata)
}
par(fig=c(0,1,0,1), new=TRUE)
map("nzHires")
points(x=gfdata$long, y=gfdata$lat,col=gfdata$col, pch=19, cex=3)
par(oma= c(0,0,0,0))
par(mar= c(0,0,1,0))
par(fig=c(0.4,0.6,0.6,0.9), new=TRUE)
barplot(gfdata$val, col=gfdata$col, border="#888888")
par(fig=oldpar$fig)
par(mar=oldpar$mar)
par(oma=oldpar$oma)
### ggplot
# Core idea is that you are building a figure in which the data is doing a job
if(!(require(ggplot2))){
install.packages("ggplot2")
require(ggplot2)
}
data(iris)
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, colour=Species, shape=Species)) +
geom_point()
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, colour=Species, shape=Species)) +
geom_point() + stat_smooth()
####
#igraph used it network analysis
if(!(require(igraph))){
install.packages("igraph")
require(igraph)
}
ig <- graph( edges=c(1, 2, 2, 3, 3, 1, 1, 4), n=4, directed=F )
plot(ig)
ig <- graph( edges=c(1, 2, 2, 3, 3, 1, 1, 4), n=6, directed=F )
plot(ig)
## also diagrammeR
