1 ########################################################################################################################
   2 # ggplot2 tutorial
   3 # created by jdegen on 02/12/2011
   4 #
   5 # requires functions "myCenter", "beehive", and "se", can be downloaded as myFunctions.R here: 
   6 # http://wiki.bcs.rochester.edu:2525/HlpLab/StatsCourses/Feb2011Regression?action=AttachFile&do=get&target=myFunctions.R
   7 # 
   8 # for ggplot2 documentation, see http://had.co.nz/ggplot2/
   9 ########################################################################################################################
  10 
  11 library(ggplot2)
  12 
  13 data(lexdec)
  14 lexdec$realRT <- exp(lexdec$RT)
  15 clex <- myCenter(lexdec)
  16 lexdec <- cbind(lexdec,clex)
  17 
  18 ###############################
  19 # HISTOGRAMS
  20 ###############################
  21 
  22 # create simple histogram of RTs
  23 ggplot(lexdec, aes(x=RT)) +
  24 	geom_histogram()
  25 
  26 # adjust bin size
  27 ggplot(lexdec, aes(x=RT)) +
  28 	geom_histogram(binwidth=0.25)	
  29 
  30 # density instead of histogram. 	
  31 ggplot(lexdec, aes(x=RT)) +
  32 	geom_density()
  33 	
  34 # density and histogram		
  35 ggplot(lexdec, aes(x=RT)) +
  36 	geom_histogram(aes(y = ..density..)) + 
  37 	geom_density()	
  38 
  39 # to plot the distributions for each subjects: facet (like conditionalizing. plots distribution for each subject subset of the data). 
  40 ggplot(lexdec, aes(x=RT)) +
  41 	geom_density() +
  42 	facet_wrap( ~ Subject)
  43 
  44 ####################################
  45 # BAR CHARTS
  46 ####################################		
  47 
  48 ## bar chart of mean reaction times by native language ##
  49 	# using stat_summary
  50 	ggplot(lexdec, aes(x=NativeLanguage,y=RT,fill=Class)) +
  51 		stat_summary(fun.y=mean, geom="bar", position="dodge") + # mean RTs by condition as bars
  52 		scale_fill_manual(values=c(rgb(1,0,0),"#483D8B")) + # change bar colors
  53 		coord_cartesian(ylim=c(5.5,7)) + # zoom in to y range
  54 		scale_y_continuous(breaks=seq(5.5,7,by=.25)) + # control y axis tick marks
  55 		xlab("Language") # change x axis label
  56 	
  57 	# same thing, but aggregate data first
  58 	agr <- with(lexdec, aggregate(RT, by=list(NativeLanguage, Class), FUN=mean)) # create data.frame with mean RTs by NativeLanguage and Class
  59 	names(agr) <- c("NativeLanguage","Class","Mean") # assign meaningful column names
  60 
  61 	# same plot as before 
  62 	ggplot(agr, aes(x=NativeLanguage, y=Mean, fill=Class)) +
  63 		geom_bar(position="dodge") +
  64 		scale_fill_manual(values=c(rgb(1,0,0),"#483D8B")) + 
  65 		coord_cartesian(ylim=c(5.5,7)) + 
  66 		scale_y_continuous(breaks=seq(5.5,7,by=.25)) + 
  67 		xlab("Language") 
  68 				
  69 ## bar plot of lmer model coefficients by effect significance ##
  70 	# run the model
  71 	mod <- lmer(RT ~ cFrequency*cNativeLanguage*cClass + (1|Subject), data=lexdec)	
  72 	# get coefficient estimates and 95% confidence intervals
  73 	pvls <- pvals.fnc(mod)
  74 	o <- as.vector(as.numeric(pvls[["fixed"]]$MCMCmean))[-1]
  75 	lo <- as.vector(as.numeric(pvls[["fixed"]]$HPD95lower))[-1]
  76 	hi <- as.vector(as.numeric(pvls[["fixed"]]$HPD95upper))[-1]
  77 	n <- names(mod@fixef)[2:length(names(mod@fixef))]
  78 
  79 	# create data.frame with data to plot & extra variable coding whether effect is significant
  80 	vls <- data.frame(Mean=o, Low=lo, High=hi, Effect=n)
  81 	vls$Significant <- ifelse(vls$Mean < 0, ifelse(vls$High < 0, "significant","nonsignificant"), ifelse(vls$Low >0, "significant", "nonsignificant"))
  82 	vls$Significant <- as.factor(as.character(vls$Significant))
  83 	vls$Effect <- as.factor(as.character(vls$Effect))
  84 
  85 	# plot
  86 	ggplot(vls, aes(x=Effect, y=Mean, fill=Significant)) + 
  87 		geom_bar(position="dodge") + # bar layer
  88 		geom_errorbar(aes(ymax=vls$High, ymin=vls$Low), width=0.25) + # error bar layer
  89 		scale_x_discrete(name="Predictor",breaks=levels(vls$Effect),labels=c("Frequency","Frequency x Native", "F x N x P", "Frequency x Previous Type", "Native", "Native x Previous Type", "Previous Type")) + # change x axis labels
  90 		scale_y_continuous(name="Coefficient") + # change y axis label
  91 		opts(axis.text.x=theme_text(size=15,vjust=0,hjust=1,angle=45),legend.position="none") # adjust orientation and size of x axis labels, remove legend
  92 		
  93 	
  94 ####################################
  95 # POINTS AND LINES
  96 ####################################
  97 
  98 # scatterplot of mean RTs by frequency
  99 ggplot(lexdec, aes(x=Frequency, y=RT)) +
 100 	stat_summary(fun.y=mean, geom="point") +
 101 	stat_summary(fun.y=mean, geom="line")
 102 
 103 # connect observations with line
 104 ggplot(lexdec, aes(x=Frequency, y=RT)) +
 105 	stat_summary(fun.y=mean, geom="point") +
 106 	stat_summary(fun.y=mean, geom="line")
 107 
 108 # fit smoothed lines through data cloud by NativeLanguage condition
 109 ggplot(lexdec, aes(x=Frequency, y=RT, colour=NativeLanguage)) +
 110 	stat_summary(fun.y=mean, geom="point") +
 111 	stat_smooth(method="lm", aes(group=NativeLanguage), colour="black")
 112 
 113 	
 114 ######################################
 115 # PLOT MODEL PREDICTIONS OF MIXED LOGIT MODELS 
 116 # use the beehive function you downloaded from the wiki, which is basically the same as the my.glmerplot function that florian wrote (but has some different defaults). documentation can be found here:
 117 http://hlplab.wordpress.com/2009/01/19/plotting-effects-for-glmer-familybimomial-models/
 118 ######################################
 119 
 120 mod <- lmer(Correct ~ cFrequency + (1|Subject),data=lexdec,family="binomial")
 121 
 122 # beehive plot in log odds space
 123 beehive(mod, "cFrequency", predictor=lexdec$Frequency, name.outcome="correct answer", xlab="Word frequency")
 124 
 125 # beehive plot in probability space
 126 beehive(mod, "cFrequency", predictor=lexdec$Frequency, name.outcome="correct answer", xlab="Word frequency",fun=plogis)
 127 	
 128 	
 129 ##########################################
 130 # PIE CHARTS
 131 ##########################################
 132 
 133 pie <- ggplot(lexdec, aes(x=factor(1), fill=Class)) +
 134 	geom_bar(width=1) +
 135 	coord_polar(theta="y")