######################################################
#Lecture 3 Script
#Graphing using ggplot2
#Anastasia Ejova 13/10/2014
#R101/PSY532
######################################################

#The workspace is the same as for Lecture/Seminar 2.

#install.packages("ggplot2")
library (ggplot2)
#setwd("~/Masaryk/Teaching/L3 - Graphing") #useful to do at the outset to save graphs

#Basic principle of ggplot2 (see slide 7): "geoms" (reps of data) have "aesthetics" 
#(properties) that can be "mapped" to variables in the dataset

#-----------------------------------------------------
#QUICK PLOTS USING qplot
#-----------------------------------------------------
?qplot

#Scatter plots
qplot(data = SS, x = PreDBC_IOC, y = Post100RoundsGoal, geom = "point") #scatter plot
#--

#Scatter plots with regression lines 
qplot(data = SS, x = PreDBC_IOC, y = Post100RoundsGoal, geom = c("point", "smooth"), 
      span = .2) #scatter plot with smoothing; small span leads to very wiggly line

qplot(data = SS, x = PreDBC_IOC, y = Post100RoundsGoal, geom = c("point", "smooth"),
      span = 1) #less wiggly line

qplot(data = SS, x = PreDBC_IOC, y = Post100RoundsGoal, geom = c("point", "smooth"), 
      method = "lm") #linear model, as in last lecture; default method is loess
      #for full list of methods, see: http://docs.ggplot2.org/current/stat_smooth.html

qplot(data = SS, x = PreDBC_IOC, y = Post100RoundsGoal, geom = c("point", "smooth"), 
      method = "lm", se = FALSE)

#install.packages("MASS")
library(MASS)
qplot(data = SS, x = PreDBC_IOC, y = Post100RoundsGoal, geom = c("point", "smooth"), 
      method = "rlm") #lm with robust fitting algorithm (not as sensitive to outliers)
#--

#Boxplots
qplot(data = SS, x = SeqCond, y= Post100RoundsGoal, geom = "boxplot") #boxplot
#for help, google "geom boxplot ggplot2" - this leads to official documentation
#and examples

#Boxplot for continuous variables
#install.packages("plyr")
library(plyr)
?round_any
qplot(data = SS, x = PreDBC_IOC, y = Post100RoundsGoal, geom = "boxplot", 
      group = round_any(x = PreDBC_IOC, accuracy = 10, 
                        f = floor))
#--

#Jittered plots
qplot(data = SS, x = SeqCond, y= Post100RoundsGoal, geom = "jitter")

qplot(data = SS, x = SeqCond, y= Post100RoundsGoal, geom = "jitter", 
      alpha = I(1/2)) #jittered plot with more transparent points

qplot(data = SS, x = SeqCond, y= Post100RoundsGoal, geom = "jitter", 
      alpha = I(1/10)) #jittered plot with even more transparent points

qplot(data = SS, x = SeqCond, y= Post100RoundsGoal, geom = c("jitter", "boxplot"), 
      alpha = I(1/2)) #boxplot overlaid
#--

#Histograms
qplot(data = SS, x = Post100RoundsGoal, geom = "histogram") #histogram
  #Note: histogram  is the default geom if only x is specified, as in:
  #qplot(data = SS, x = Post100RoundsGoal)

qplot(data = SS, x = Post100RoundsGoal, geom = "histogram", binwidth = 1)

qplot(data = SS, x = Post100RoundsGoal, geom = "histogram", binwidth = 10)

qplot(data = SS, x = Post100RoundsGoal, geom = "histogram", binwidth = 10, 
      fill = SeqCond)

qplot(data = SS, x = PostHowManySingleWins, geom = "histogram", binwidth = 5, 
      fill = PostHowManySingleCaptionType, alpha = .5, position = "identity")

qplot(data = SS, x = PostHowManySingleWins, geom = "histogram", binwidth = 5, 
      fill = PostHowManySingleCaptionType, alpha = .5, position = "dodge")
  #With position = dodge, you probably don't need transparency.

qplot(data = SS, x = Post100RoundsGoal, geom = "density") #density plot - smooth of
  #histogram

qplot(data = SS, x = Post100RoundsGoal, geom = "density", adjust = 10)

qplot(data = SS, x = Post100RoundsGoal, geom = "density", adjust = 10, fill = SeqCond)

qplot(data = SS, x = Post100RoundsGoal, geom = "density", adjust = 10, 
      colour = SeqCond)
#--

#Bar charts
qplot(data = SS, x = Post100RoundsGoal, geom = "bar") #bar chart - same as histogram

qplot(data = SS, x = SeqCond, y = PostHowManySingleWins, geom = "bar", stat = "summary",
      fun.y = "mean") #bar chart with categories on x axis

qplot(data = SS, x = SeqCond, y = PostHowManySingleWins, geom = "bar", stat = "summary",
      fun.y = "mean", width = .5)

qplot(data = SS, x = SeqCond, y = PostHowManySingleWins, 
      geom = "bar", 
      position = "dodge", #without this, would have been a stacked bar chart; now: clustered bar chart
      stat = "summary",
      fun.y = "mean", 
      fill = PostHowManySingleCaptionType)
#--

#Faceting
qplot(data = SS, x = SeqCond, y = PostHowManySingleWins, geom = "bar", 
      stat = "summary",
      fun.y = "mean", facets = . ~ PostHowManySingleCaptionType)

qplot(data = SS, x = SeqCond, y = PostHowManySingleWins, geom = "bar", 
      stat = "summary",
      fun.y = "mean", facets = PostHowManySingleCaptionType ~ .) 

qplot(data = SS, x = SeqCond, y = PostHowManySingleWins, geom = "bar", 
      stat = "summary",
      fun.y = "mean", facets = ~ PostHowManySingleCaptionType) 
  #facet_wrap as opposed to facet_grid (this and the next two examples)

qplot(data = SS, x = SeqCond, y = PostHowManySingleWins, geom = "bar", 
      stat = "summary",
      fun.y = "mean", facets = PreSoccerInterest ~ PostHowManySingleCaptionType) 

qplot(data = SS, x = SeqCond, y = PostNaturalIoC, geom = "bar", 
      stat = "summary",
      fun.y = "mean", facets = ~ PreSoccerInterest + Gender)
#--

#Time series (line) plot
subset <- SS2[,c("PNo", "SeqCond", "PlayerChange_First24",
                 "PlayerChange_Last24", "DirEntropy_First24",
                 "DirEntropy_Last24")]

#install.packages(reshape)
library(reshape)
longsub <- reshape(subset, idvar = c("SeqCond", "PNo"),
               varying = rbind(c("PlayerChange_First24",  
                                 "PlayerChange_Last24"),
                               c("DirEntropy_First24", 
                                 "DirEntropy_Last24")),
               v.names = c("PlayerChange", "DirEntropy"),
               timevar = "Time",
               times = c(1,2),
               direction = "long")

qplot(x = Time, y = mean(PlayerChange), data = longsub, 
      geom = c("point", "path")) #let's improve this!

qplot(x = Time, y = PlayerChange, data = longsub, 
      geom = c("point", "path"), colour = SeqCond, 
      stat = "summary",
      fun.y = "mean")
  #comment out stat = "summary" and see what happens: the means are no longer plotted

qplot(x = Time, y = PlayerChange, data = longsub, 
      geom = c("point", "line"), colour = SeqCond, 
      stat = "summary",
      fun.y = "mean")
  #"path" and "line" are interchangeable in this case...

#...but not in this case
qplot(data = SS, x = PreDBC_IOC, y = Post100RoundsGoal, geom = "path")
qplot(data = SS, x = PreDBC_IOC, y = Post100RoundsGoal, geom =  "line")
  #although I'm not sure you'd ever want to draw these hedgehog plots :)
#--

#More on the notion of 'mapping'
qplot(x = Time, y = PlayerChange, data = longsub, 
      geom = c("point", "path"), 
      linetype = SeqCond, #mapping SeqCond on linetype instead of colour
      stat = "summary",
      fun.y = "mean") 

qplot(x = Time, y = PlayerChange, data = longsub, 
      geom = c("point", "path"), 
      shape = SeqCond, #mapping SeqCond on point shape
      stat = "summary",
      fun.y = "mean") 

#Saving plots: demonstration, as per lecture slide


#-----------------------------------------------------
#MORE FLEXIBLE PLOTS USING ggplot2
#-----------------------------------------------------
#Added advantages: Layering and easy changes (or overriding) to settings.
#Each plot treated as a variable.

#Layering
  p1 <- ggplot(data = SS, mapping = aes(x= PreDBC_IOC, y=Post100RoundsGoal)) 
  p1 #displays nothing... until we add a layer
  
  p1 <- p1 + geom_point()
  p1
  
  p2 <- p1 + geom_smooth(method = "lm")
  p2
  
  p3 <- p1 + stat_smooth(method = "lm")
  p3 #all stats have a default geom and vice versa
  
  p4 <- p1 + geom_smooth(method = "lm", colour = "red", fill = "red") #adding settings
  p4 
  
  p5 <- p1 + layer(geom = "smooth", 
                   geom_params = list(colour = "red", fill = "red"), 
                   stat_params = list(method = "lm")) 
    #alternative specification of geom_smooth(); resulting graph is the same
  p5
  
  #Advantages of layering
  #Overlaid regression lines: least squares (lm) and robust
  library(MASS)
  p6 <- p5 + geom_smooth(colour = "blue", fill = "blue", method = "rlm")
  p6

  #Overlaid histograms
  #Normal distribution
  library(psych)
  describe(SS$Post100RoundsGoal)
  #mean = 21.26, sd = 16.02

  histo <- ggplot(data=SS, aes(x=Post100RoundsGoal)) + 
  geom_histogram(aes(y = ..density..)) +
  stat_function(fun = dnorm, colour="red", 
                args=list(mean=21.26, sd = 16.02))
  histo
  #The y-axis of the histogram does not come from the original data, but
  #from the statistical transformation that counts the number of observations in
  #each bin when stat bin (the default for histograms) is applied. 
  #Using ..density.. tells ggplot2 to  map the density to the y-axis instead of using
  #the default count also calculated by stat bin.
#--

#Adding error bars
#The trick is to make use of a 'stat' layer instead of a 'geom' layer at one point.
#Stat_summary has a default associated geom "pointrange", and we have to change
#that default to "errorbar". This function also automatically calculates the min
#and max values for CIs around relevant means in the dataframe.
#Details: http://rstudio-pubs-static.s3.amazonaws.com/10674_02d3c0dbc5ff4cab8c444015d3c735d1.html
#You'd want to compare your visible CIs with CIs calculated using a descriptive stats 
#function (e.g., group.CI in Rmisc package).

#Another useful reference if happy to write own functions:
#http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)

#install.packages("Rmisc")
library(Rmisc)
group.CI(x = PlayerChange~SeqCond + Time, longsub, ci = 0.95)

#cip: ci plot
cip <- ggplot(data = longsub, aes(x=Time, y=PlayerChange, colour=SeqCond)) + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min, #min and max are calculated by stat_summary
               ymax = max, 
               width = .05) + #width refers to size of horiz ends of error bars
  geom_point(stat = "summary", stat_params = list(fun.y = mean)) +
  geom_line(stat = "summary", stat_params = list(fun.y = mean))

cip #note: "mean_cl_boot" also available as an option

#To fix overlapping error bars and x axis
pd <- position_dodge(width = .1) # ensure a distance of .1 (on x-axis) between 
#points/bars; http://docs.ggplot2.org/0.9.3.1/position_dodge.html

cipd <- ggplot(data = longsub, aes(x=Time, y=PlayerChange, colour=SeqCond)) + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .05,
               position = pd) + #aadd this line to previous graph
  geom_point(stat = "summary", stat_params = list(fun.y = mean), position = pd) +
  geom_line(stat = "summary", stat_params = list(fun.y = mean), position = pd) +
  scale_x_continuous(breaks = 1:2, labels = c("Trials 1-24", "Trials 25-48"))

cipd

#sep: standard error plot - almost the same code
sep <- ggplot(data = longsub, aes(x=Time, y=PlayerChange, colour=SeqCond)) + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .05,
               mult = 1, #add this to calculate one se from mean; http://stackoverflow.com/questions/19258460/standard-error-bars-using-stat-summary
               position = pd) + 
  geom_point(stat = "summary", stat_params = list(fun.y = mean), position = pd) +
  geom_line(stat = "summary", stat_params = list(fun.y = mean), position = pd) +
  scale_x_continuous(breaks = 1:2, labels = c("Trials 1-24", "Trials 25-48"))

sep  

#to keep y-axis range fixed 
sep <- sep + coord_cartesian(ylim=c(0,3)) #do not use scale_y_continuous; otherwise, 
#data outside limit range will be omitted
#i.e., do not use sep + scale_y_continuous(limits = c(0,3))
sep 


cipd <- cipd + coord_cartesian(ylim=c(0,3))
cipd


#bse: bar plot with standard errors
bse <- ggplot(data = SS, aes(x=SeqCond, y=PostHowManySingleWins, 
                              fill=PostHowManySingleCaptionType)) +
  stat_summary(fun.y = "mean", geom = "bar", position = "dodge") + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .5,
               mult = 1,
               position = "dodge")

bse

#error bars don't dodge to the same degree as bars, so experiment with different
#widths for dodging (http://docs.ggplot2.org/0.9.3.1/position_dodge.html)
pd2 <- position_dodge(width = .9) 

bse <- ggplot(data = SS2, aes(x=SeqCond, y=PostHowManySingleWins, 
                              fill=PostHowManySingleCaptionType)) +
  stat_summary(fun.y = "mean", geom = "bar", position = "dodge") + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .5,
               mult = 1, 
               position = pd2)

bse
#--

#Faceting
bsef1 <- ggplot(data = SS, aes(x=SeqCond, y=PostHowManySingleWins)) + 
  facet_grid(PostHowManySingleCaptionType ~ .) + #Facet grid: column
  stat_summary(fun.y = "mean", geom = "bar") + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .5,
               mult = 1,
               position = pd2)
bsef1

bsef2 <- ggplot(data = SS, aes(x=SeqCond, y=PostHowManySingleWins)) + 
  facet_wrap(~ PostHowManySingleCaptionType) + #Facet wrap
  stat_summary(fun.y = "mean", geom = "bar") + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .5,
               mult = 1,
               position = pd2)
bsef2

bsef3 <- ggplot(data = SS, aes(x=SeqCond, y=PostHowManySingleWins)) +
  facet_grid(. ~ PostHowManySingleCaptionType) + #Facet grid: row
  stat_summary(fun.y = "mean", geom = "bar") + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .5,
               mult = 1,
               position = pd2)
bsef3

bsef4 <- ggplot(data = SS, aes(x=SeqCond, y=PostHowManySingleWins)) +
  facet_grid(. ~ PostHowManySingleCaptionType, 
             margins = c("PostHowManySingleCaptionType", 
                         "grand_row")) + #adding a "grand row" and/or "grand column"
  stat_summary(fun.y = "mean", geom = "bar") + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .5,
               mult = 1,
               position = pd2)
bsef4
#--

#Faceting with free scales
#Useful when you want to display two trends side-by-side, but the scales of the y
#axes differ (e.g., one trend is visible over millions of dollars
#while the other is visible in single dollars). 
#In our example below, player changes differ on a scale of 1-5 while direction
#entropy is measured on a logarithmic scale, ranging from 1 to 2. We will see in the
#graph that, despite the differnces in scales, there's a downwards trend over time
#in the Ascending condition.

#First, need to combine all DVs into one variable
longsub2 <- reshape(subset, idvar = c("SeqCond", "PNo"),
               varying = rbind(c("PlayerChange_First24",
                                 "DirEntropy_First24",
                                 "PlayerChange_Last24",
                                 "DirEntropy_Last24")),
               v.names = c("Value"),
               timevar = c("Time", "Measure"),
               times = c(c(1,1), c(2,2)),
               direction = "long")

longsub2[c(1:176, 353:528),"Measure"] <- "PlayerChange"
longsub2[c(177:352, 529:704),"Measure"] <- "DirEntropy"

#Then, it's quite easy to draw the graph, keeping the y scale free and otherwise 
#using the code from sep
pd <- position_dodge(width = .1)

sepf <- ggplot(data = longsub2, aes(x=Time, y=Value, colour=SeqCond)) + 
  facet_grid(Measure ~ ., scales = "free_y") + #y scale is free
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .05,
               mult = 1,
               position = pd) + 
  geom_point(stat = "summary", stat_params = list(fun.y = mean), position = pd) +
  geom_line(stat = "summary", stat_params = list(fun.y = mean), position = pd) +
  scale_x_continuous(breaks = 1:2, labels = c("Trials 1-24", "Trials 25-48"))

sepf  

#-----------------------------------------------------
#FINAL FORMATTING USING ggplot
#-----------------------------------------------------
#Formatting for presentations and journal papers.

#Themes
#Can be set globally
theme_set(theme_bw()) #theme_grey() with grey background is the default
sepf #grey background gone and border has appeared as required by APA
besf1 #likewise

#... or locally
theme_set(theme_grey())
sepf2 <- sepf + theme_bw()
sepf2
sepf #still in globally set grey theme

#but notice that, even though the theme is b+w, the colour has stayed...
#--

#Colour scale (i.e., colour scheme)
sepf2 + scale_colour_brewer() #a different preset colour palette
sepf2 + scale_colour_grey() #a grayscale palette
sepf2 + scale_colour_hue() #default palette

#A darker grayscale palette is possible...
sepf2 + scale_colour_grey(start = 0.1, end = 0.4) #default start is start .2, end .8

#But, to accommodate b+w it's best to avoid mappings of variables on colour
#Let's map onto linetype instead
theme_set(theme_bw())

sepf3 <- ggplot(data = longsub2, aes(x=Time, y=Value, 
                                group = SeqCond)) + #group replaces colour mapping
  geom_line(mapping = aes(linetype = SeqCond), #mapping to linetype here
            stat = "summary", 
            stat_params = list(fun.y = mean), position = pd) + #moved to top as well
  facet_grid(Measure ~ ., scales = "free_y") + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .05,
               mult = 1, 
               position = pd) + 
  geom_point(stat = "summary", stat_params = list(fun.y = mean), position = pd) +
  scale_x_continuous(breaks = 1:2, labels = c("Trials 1-24", "Trials 25-48"))

sepf3

#Changing the appearance of points, errorbars and lines
#shape options for points: http://sape.inf.usi.ch/quick-reference/ggplot2/shape
#colours: http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
sepf3 <- ggplot(data = longsub2, aes(x=Time, y=Value, group = SeqCond)) +
  geom_line(mapping = aes(linetype = SeqCond), stat = "summary", 
            stat_params = list(fun.y = mean), position = pd, 
            geom_params = list(size = 1)) + #making line thicker
  facet_grid(Measure ~ ., scales = "free_y") + 
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", 
               ymin = min,
               ymax = max, 
               width = .05,
               mult = 1, 
               position = pd,
               colour = "gray41") + #changing error bar colour; no need for geom_params call as for geom_lines
  geom_point(stat = "summary", stat_params = list(fun.y = mean), position = pd,
             geom_params = list(fill = "white", shape = 21, size = 3)) + #changing point settings
  scale_x_continuous(breaks = 1:2, labels = c("Trials 1-24", "Trials 25-48"))

sepf3

sepf3 <- sepf3 + scale_linetype_manual(values=c("dotted", "solid", "longdash", 
                                                "dotdash"))
sepf3
#adjust linetype scale as per http://www.cookbook-r.com/Graphs/Shapes_and_line_types/

#Next week:
  #axis labels
  #text direction
  #legends - text and position
  #titles

