# R code for chapter 57 Publication quality graphics
# Book details:
# Health Science Statistics using R and R Commander
# by Robin Beaumont
# Paperback - ISBN 9781907904318
# http://www.amazon.co.uk/Health-Science-Statistics-using-Commander/dp/190790431X
# or direct from publisher at: 
# http://www.scionpublishing.com/
# book website: http://www.robin-beaumont.co.uk/rbook
###########################################################
#
# install.packages("ggplot2",dependencies=TRUE)
library(ggplot2)
# 57.1.3	Boxplots
# boxplots
# weightgain
data(anorexia, package="MASS")
anorexia$weightgain <- with(anorexia, Postwt - Prewt)
anorexia

qplot(Treat, weightgain,  geom="boxplot", data=anorexia)
# same as:
ggplot(data=anorexia, aes(x=Treat, y=weightgain)) + 
geom_boxplot()

# add notches
qplot(Treat, weightgain,  geom="boxplot", data=anorexia, notch=TRUE,varwidth=TRUE)

######## boxplots for multiple variables
# get data locally - ignore:
# mydataframe <-  read.table("D:/web_sites_mine/HIcourseweb new/book2data/wilcoxon_paired_ex1.dat",
#  header=TRUE, sep="\t", na.strings="NA", dec=".", strip.white=TRUE)
# get data on server:
   mydataframe <- 
   read.delim("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/wilcoxon_paired_ex1.dat") 
 mydataframe  
mydataframe$difference <- mydataframe$post - mydataframe$pre
mydataframe$difference
library(reshape2)
# melt all three columns into one, where the melt function 
# automatically renames the two new columns, variable and value
# place it in a new dataframe called longdf
longdf <- melt(mydataframe)
longdf
library(ggplot2)

qplot(variable, value,  geom="boxplot", data=longdf, 
main="Example from comparing pre-post medians chapter")
# or
ggplot(data=longdf) + 
geom_boxplot(aes(x=variable,y=value))

#######
# add jitter
qplot(variable, value, geom=c("boxplot", "jitter"), data=longdf)
# add mean
ggplot(data=longdf, aes(x=variable, y=value)) + 
geom_boxplot() +
stat_summary(fun.y="mean", geom="point", shape=23, size=10, fill="black")

############ online example using facet - for factor levels
data <- data.frame(v1=rnorm(100),v2=rnorm(100),v3=rnorm(100), v4=rnorm(100))
library(reshape)
meltData <- melt(data)
boxplot(data=meltData, value~variable)
# or even then use ggplot2 package to make things nicer
library(ggplot2)
p <- ggplot(meltData, aes(factor(variable), value)) 
p + geom_boxplot() + facet_wrap(~variable, scale="free")
################
# 57.1.4	 Scatterplots – simple
# Scatterplot simple
# get data from url:
mydataframe <- 
read.delim("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/correlation_vo2max.dat")
########## get data locally - ignore
# vo2max data correlation_vo2max.dat in book2data folder
# mydataframe <- read.delim("D:\\web_sites_mine\\HIcourseweb new\\book2data\\correlation_vo2max.dat", header=TRUE)
names(mydataframe)
# [1] "id"       "age"      "sex"      "season"   "height"   "weight"  
 #[7] "time"     "tricep1"  "tricep2"  "tricep3"  "bicep1"   "bicep2"  
#[13] "bicep3"   "subscap1" "subscap2" "subscap3" "supra1"   "supra2"  
#[19] "supra3"   "vjwith"   "vjwout"   "flexboth" "flexleft" "flexrigh"
#[25] "griprigh" "gripleft" "abdcurln" "abdcurlt" "mftlev"   "mftshut" 
#[31] "mftvomax" "mst1"     "mst2"     "mst3"     "mst4"     "mst5"    
#[37] "mst6"     "mst7"     "mst8" 
   
library(ggplot2)
   qplot(abdcurln, height, data = mydataframe, main="height v. abdominal curl count")
qplot(abdcurln, height, data = mydataframe, 
main="height v. abdominal curl count")
# or
ggplot(data=mydataframe, aes(x=abdcurln, y=height)) + 
geom_point() +
labs(title = "height v. abdominal curl count")

# subgrouping
qplot(abdcurln, height, data = mydataframe,
main="height v. abdominal curl count", 
shape = factor(sex))

# or
ggplot(data=mydataframe, aes(x=abdcurln, y=height)) + 
geom_point() +
labs(title = "height v. abdominal curl count") +
aes(shape = factor(sex))
#
# enlarging in qplot
qplot(abdcurln, height, data = mydataframe,
main="height v. abdominal curl count", 
shape = factor(sex),
size=I(5))
or 
ggplot(data=mydataframe, aes(x=abdcurln, y=height)) + 
geom_point(size=5) +
labs(title = "height v. abdominal curl count") +
aes(shape = factor(sex)) +
scale_shape(solid = FALSE)

# 57.1.5	Linear Regression lines
######### regression lines
qplot(abdcurln, height, data = mydataframe, method="lm",
shape = factor(sex),
geom = c("point", "smooth") )
# or
ggplot(data=mydataframe, aes(x=abdcurln, y=height)) + 
aes(shape = factor(sex)) +
geom_point(size=5) +
stat_smooth( method="lm") 
#################
# 57.1.6	Scatterplot matrix
#scatterplot matrix using the car package
library(car)
 scatterplotMatrix(~abdcurln+height+mftvomax+weight, reg.line=lm, 
  smooth=FALSE, spread=FALSE, span=0.5, id.n=0, diagonal = 'histogram', 
  data=mydataframe)
##################
###### if determined to use GGally:
# GGally has a turn if there are any ns in data
str(mydataframe) 

reduced <- mydataframe[, c("abdcurln","height","mftvomax","weight")]
reduced
str(mydataframe_complete) 
############ ignore next 2 lines
## line of code below defunct plotmatrix last used in version 0.9.2:
# plotmatrix(mydataframe[, c("abdcurln","height","mftvomax","weight")])
############
install.packages("GGally")
library(GGally)
# first showing the errors that occur in the R Console window
# due to missing values
# previous versions just crashed now you get the graphics
ggpairs(reduced, lower = list(continuous = "smooth"))
###
ggpairs(mydataframe_complete, lower = list(continuous = "smooth"))

######  correlations for comparison
use="pairwise.complete.obs"
round(cor(reduced), 4)
round(cor(reduced,use="pairwise.complete.obs"), 4)
round(cor(mydataframe_complete), 4)
round(cor(mydataframe_complete, use="pairwise.complete.obs"), 4)

####  example from GGally package authors
library(GGally)
ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
ds$sex = as.factor(ifelse(ds$female==1, "female", "male"))
ds$housing = as.factor(ifelse(ds$homeless==1, "homeless", "housed"))
smallds = subset(ds, select=c("housing", "sex", "i1", "cesd"))
ggpairs(smallds, diag=list(continuous="density", discrete="bar"), axisLabels="show")

ggpairs(ds, columns=c("housing", "sex", "i1", "cesd"),
    diag=list(continuous="density",   discrete="bar"), axisLabels="show")

# 57.1.7	Lines for repeated measures
##################################
# repeated measures example from
# http://egret.psychol.cam.ac.uk/statistics/R/graphs2.html
# have downloaded csv files to book2data directory then set it as the working directory
#  data prep
fig30A<- read.csv(("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/MD_fig30A.csv"))
fig30A
fig30B<- read.csv(("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/MD_fig30B.csv"))
fig30B
################  To copy FROM chapter  . . .
mdfig30_long = rbind( data.frame(type="programmed",  delay=mdfig30a$delay, lesion="sham", 
sqrtpresses = mdfig30a$sham, ysem = mdfig30a$SEMsham, xsem = 0 ),
data.frame(type="programmed",  delay=mdfig30a$delay, lesion="AcbC", sqrtpresses = mdfig30a$AcbC, ysem = mdfig30a$SEMAcbC, xsem = 0 ),
data.frame(type="experienced", delay=mdfig30b$shamdelay, lesion="sham", sqrtpresses = mdfig30b$shamresponses, 
ysem = mdfig30b$SEMshamresponses, xsem = mdfig30b$SEMshamdelay ),
data.frame(type="experienced", delay=mdfig30b$AcbCdelay, lesion="AcbC", sqrtpresses = mdfig30b$AcbCresponses, 
ysem = mdfig30b$SEMAcbCresponses, xsem = mdfig30b$SEMAcbCdelay ) ) 

# 57.1.7.2	Developing the plot
## the plot adapted from original code - removed themes aspects
the_plot = ggplot(data = mdfig30_long, aes(x = delay, y = sqrtpresses, fill = lesion) ) +
     facet_grid(. ~ type) +
     geom_line() +
     geom_errorbar(aes(ymin=sqrtpresses - ysem, ymax=sqrtpresses + ysem), width=0.8) +
     geom_errorbarh(aes(xmin=delay - xsem, xmax=delay + xsem), height=0.05) + 
     geom_point(shape=21, colour="black", size=5) +
     scale_fill_manual(values=c("white","black")) +
     scale_x_continuous("Delay from response to reinforcer (s)") +
     scale_y_continuous(expression(sqrt("active lever presses per minute in session 10")), 
	 expand=c(0,0), limits = c(0, 3.5), breaks=0:3)   
the_plot 

###########################################
# original from cambridge website
mdfig30a = read.table("MD_fig30A.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE)
mdfig30b = read.table("MD_fig30B.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE)

mdfig30_long = rbind( data.frame(type="programmed",  delay=mdfig30a$delay, lesion="sham", 
sqrtpresses = mdfig30a$sham, ysem = mdfig30a$SEMsham, xsem = 0 ),
data.frame(type="programmed",  delay=mdfig30a$delay, lesion="AcbC", sqrtpresses = mdfig30a$AcbC, ysem = mdfig30a$SEMAcbC, xsem = 0 ),
data.frame(type="experienced", delay=mdfig30b$shamdelay, lesion="sham", sqrtpresses = mdfig30b$shamresponses, ysem = mdfig30b$SEMshamresponses, xsem = mdfig30b$SEMshamdelay ),
data.frame(type="experienced", delay=mdfig30b$AcbCdelay, lesion="AcbC", sqrtpresses = mdfig30b$AcbCresponses, ysem = mdfig30b$SEMAcbCresponses, xsem = mdfig30b$SEMAcbCdelay ) )
#  

mdfig30_annotations = data.frame(
    x = c(0,10,20,10),
    y = c(3.2,1.7,1.2,1.4),
    type = c("programmed", "programmed", "programmed", "experienced"),
    label = c("*", "*", "***", "###"),
    lesion = c("sham", "sham", "sham", "sham") ) # irrelevant but must not add new levels or they'll go into the legend

f6 = ggplot(data = mdfig30_long, aes(x = delay, y = sqrtpresses, fill = lesion) ) +
     facet_grid(. ~ type) +
     geom_line() +
     geom_errorbar(aes(ymin=sqrtpresses - ysem, ymax=sqrtpresses + ysem), width=0.8) +
     geom_errorbarh(aes(xmin=delay - xsem, xmax=delay + xsem), height=0.05) + # somewhat hard to get it to match other error bar size (different scale!)
     geom_point(shape=21, colour="black", size=5) +
     scale_fill_manual(values=c("white","black")) +
     scale_x_continuous("Delay from response to reinforcer (s)") +
     scale_y_continuous(expression(sqrt("active lever presses per minute in session 10")), expand=c(0,0), limits = c(0, 3.5), breaks=0:3) +
     geom_text(data = mdfig30_annotations, aes(x = x, y = y, label = label) ) +
     theme_bw() +
     opts(plot.margin=unit(c(1,1,1.5,2), "lines"),
          axis.title.x = theme_text(size=12, vjust=-1),
          axis.title.y = theme_text(size=12, angle=90, vjust=0), # likewise
          legend.title = theme_blank(),
          legend.key.size = unit(1.5, "lines"),
          legend.key = theme_blank()
     )
print(f6)
##################################################
# 57.1.8	Kaplan-Meier plots
# see book
##################################################
# 57.1.9	Drawing distributions and adding text
# 57.1.9.1	Creating the distributions and formatted text
xval  <- 1:9
dataf <- data.frame(group="cont", x = xval, y = dbinom(xval, 9, 0.5))
dataf2 <- data.frame(group="ZMapp", x = xval, y = dbinom(xval, 9, 0.7))
totaldata <- rbind(dataf, dataf2)
totaldata
We now obtain the binomial test results we want to show on the plot:
result_ZMapp <- binom.test(7, 9, alternative='two.sided', p=.5, conf.level=.95)
mytext <- paste ("p-value= ", result_ZMapp$p.value)
mytext2 <- paste ("95%CI= ", round(result_ZMapp$conf.int[1],3), "to", round(result_ZMapp$conf.int[2],3 ))
mytext; mytext2
# To add some text with subscripts and the greek letter pi (=?):
text3 <- expression(survival~proportions:~hat(pi)[0.5]~versus~hat(pi)[0.7])

# 57.1.9.2	Plotting the distributions and formatted text
library(ggplot2)
library(plyr) 
plot <- ggplot(totaldata, 
aes(x,y, colour=group, fill=group)) +
geom_bar(subset=.(group=="ZMapp"),stat="identity")+
geom_line(subset=.(group=="cont"),size = 1) +
scale_y_continuous(expand = c(0.1, 0)) + 
scale_x_continuous(breaks=c(3,5,7), labels=c("0.3", "0.5", "0.7")) +
xlab(text3) + 
ylab("Density") + 
labs(title = "Ebola West Africa 1/09/2014 \n Survival rates traditional and ZMapp (n=9)") + 
theme_grey(base_size = 14, base_family = "sans") + 
theme(plot.title = element_text(size = rel(1.2), vjust = 1.5)) +
annotate("text",x=7,y=.32,label=mytext, fontface = 3, size=4) + 
annotate("text",x=7,y=.3,label=mytext2, fontface = 3, size=4)   + 
scale_fill_manual(values=c("red", "blue")) +
scale_colour_manual(values=c("red", "blue"))
print(plot)
###
# To produce a greyscale version of the above requires the editing 
# of the scale_fill_manual() and  scale_colour_manual() lines, show below:
plot <- ggplot(totaldata,
    . . . . .
    . . . . . .
 scale_fill_manual(values=c("gray40", "gray75")) +
 scale_colour_manual(values=c("gray40", "gray75"))
 print(plot)
 ###
 # To  change the labels in the legend you can edit the scale_fill_maual() function, 
 # adding a labels= option:
scale_fill_manual(values=c("gray40", "gray75"), labels=c("trad care","experimental drug"))

# 57.2	Using point and click front ends that generate ggplot2 plots 

# 57.2.1	The R Commander plugin for ggplot2 plots - RcmdrPlugin.KMggplot2
install.packages("tcltk2") # specifically you need the latest version of this 
install.packages("RcmdrPlugin.KMggplot2")
library(Rcmdr)
mydataframe <- read.delim("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/survival_collett_p290.dat", header=TRUE)
names(mydataframe); mydataframe
mydataframe$treat <- as.factor(mydataframe$treat)

# 57.2.2	Deducer
install.packages("Deducer", dependencies=TRUE)
library(Deducer)