########################################################################################################################
# ggplot2 tutorial
# created by jdegen on 02/12/2011
#
# requires functions "myCenter", "beehive", and "se", can be downloaded as myFunctions.R here: 
# http://wiki.bcs.rochester.edu:2525/HlpLab/StatsCourses/Feb2011Regression?action=AttachFile&do=get&target=myFunctions.R
# 
# for ggplot2 documentation, see http://had.co.nz/ggplot2/
########################################################################################################################

library(ggplot2)

data(lexdec)
lexdec$realRT <- exp(lexdec$RT)
clex <- myCenter(lexdec)
lexdec <- cbind(lexdec,clex)

###############################
# HISTOGRAMS
###############################

# create simple histogram of RTs
ggplot(lexdec, aes(x=RT)) +
	geom_histogram()

# adjust bin size
ggplot(lexdec, aes(x=RT)) +
	geom_histogram(binwidth=0.25)	

# density instead of histogram. 	
ggplot(lexdec, aes(x=RT)) +
	geom_density()
	
# density and histogram		
ggplot(lexdec, aes(x=RT)) +
	geom_histogram(aes(y = ..density..)) + 
	geom_density()	

# to plot the distributions for each subjects: facet (like conditionalizing. plots distribution for each subject subset of the data). 
ggplot(lexdec, aes(x=RT)) +
	geom_density() +
	facet_wrap( ~ Subject)

####################################
# BAR CHARTS
####################################		

## bar chart of mean reaction times by native language ##
	# using stat_summary
	ggplot(lexdec, aes(x=NativeLanguage,y=RT,fill=Class)) +
		stat_summary(fun.y=mean, geom="bar", position="dodge") + # mean RTs by condition as bars
		scale_fill_manual(values=c(rgb(1,0,0),"#483D8B")) + # change bar colors
		coord_cartesian(ylim=c(5.5,7)) + # zoom in to y range
		scale_y_continuous(breaks=seq(5.5,7,by=.25)) + # control y axis tick marks
		xlab("Language") # change x axis label
	
	# same thing, but aggregate data first
	agr <- with(lexdec, aggregate(RT, by=list(NativeLanguage, Class), FUN=mean)) # create data.frame with mean RTs by NativeLanguage and Class
	names(agr) <- c("NativeLanguage","Class","Mean") # assign meaningful column names

	# same plot as before 
	ggplot(agr, aes(x=NativeLanguage, y=Mean, fill=Class)) +
		geom_bar(position="dodge") +
		scale_fill_manual(values=c(rgb(1,0,0),"#483D8B")) + 
		coord_cartesian(ylim=c(5.5,7)) + 
		scale_y_continuous(breaks=seq(5.5,7,by=.25)) + 
		xlab("Language") 
				
## bar plot of lmer model coefficients by effect significance ##
	# run the model
	mod <- lmer(RT ~ cFrequency*cNativeLanguage*cClass + (1|Subject), data=lexdec)	
	# get coefficient estimates and 95% confidence intervals
	pvls <- pvals.fnc(mod)
	o <- as.vector(as.numeric(pvls[["fixed"]]$MCMCmean))[-1]
	lo <- as.vector(as.numeric(pvls[["fixed"]]$HPD95lower))[-1]
	hi <- as.vector(as.numeric(pvls[["fixed"]]$HPD95upper))[-1]
	n <- names(mod@fixef)[2:length(names(mod@fixef))]

	# create data.frame with data to plot & extra variable coding whether effect is significant
	vls <- data.frame(Mean=o, Low=lo, High=hi, Effect=n)
	vls$Significant <- ifelse(vls$Mean < 0, ifelse(vls$High < 0, "significant","nonsignificant"), ifelse(vls$Low >0, "significant", "nonsignificant"))
	vls$Significant <- as.factor(as.character(vls$Significant))
	vls$Effect <- as.factor(as.character(vls$Effect))

	# plot
	ggplot(vls, aes(x=Effect, y=Mean, fill=Significant)) + 
		geom_bar(position="dodge") + # bar layer
		geom_errorbar(aes(ymax=vls$High, ymin=vls$Low), width=0.25) + # error bar layer
		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
		scale_y_continuous(name="Coefficient") + # change y axis label
		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
		
	
####################################
# POINTS AND LINES
####################################

# scatterplot of mean RTs by frequency
ggplot(lexdec, aes(x=Frequency, y=RT)) +
	stat_summary(fun.y=mean, geom="point") +
	stat_summary(fun.y=mean, geom="line")

# connect observations with line
ggplot(lexdec, aes(x=Frequency, y=RT)) +
	stat_summary(fun.y=mean, geom="point") +
	stat_summary(fun.y=mean, geom="line")

# fit smoothed lines through data cloud by NativeLanguage condition
ggplot(lexdec, aes(x=Frequency, y=RT, colour=NativeLanguage)) +
	stat_summary(fun.y=mean, geom="point") +
	stat_smooth(method="lm", aes(group=NativeLanguage), colour="black")

	
######################################
# PLOT MODEL PREDICTIONS OF MIXED LOGIT MODELS 
# 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:
http://hlplab.wordpress.com/2009/01/19/plotting-effects-for-glmer-familybimomial-models/
######################################

mod <- lmer(Correct ~ cFrequency + (1|Subject),data=lexdec,family="binomial")

# beehive plot in log odds space
beehive(mod, "cFrequency", predictor=lexdec$Frequency, name.outcome="correct answer", xlab="Word frequency")

# beehive plot in probability space
beehive(mod, "cFrequency", predictor=lexdec$Frequency, name.outcome="correct answer", xlab="Word frequency",fun=plogis)
	
	
##########################################
# PIE CHARTS
##########################################

pie <- ggplot(lexdec, aes(x=factor(1), fill=Class)) +
	geom_bar(width=1) +
	coord_polar(theta="y")