###########################################
##                                       ##
##  Bar charts as a means of presenting  ##
## multiple samples of quantitative data ##
##                                       ##
###########################################

###########################################################################

##      BAR CHARTS    ##
##  WITH PIG LITTERS  ##
##        DATA        ##

piglets <- read.table(file.choose(), header = T, sep = ",") # load the data in from the Excel CSV. file
                                                            # 'pig_litters' and name it 'piglets'
View(piglets) # view the data

cleanpiglets <- lapply(piglets, function(col)col[!is.na(col)]) # remove NA values
View(cleanpiglets) # view the cleaned data

lit1<-cleanpiglets$litter1 # give each sample to plot a convenient name
lit2<-cleanpiglets$litter2
lit3<-cleanpiglets$litter3
lit4<-cleanpiglets$litter4
lit5<-cleanpiglets$litter5
lit6<-cleanpiglets$litter6
lit7<-cleanpiglets$litter7

n <- c(length(lit1),length(lit2),length(lit3),length(lit4),
       length(lit5),length(lit6),length(lit7)) # create list of sample sizes

n # view n

### Creating a simple bar chart          

y.means <- c(mean(lit1),mean(lit2),mean(lit3),mean(lit4),
             mean(lit5),mean(lit6),mean(lit7)) # create a list of the mean birth weights

mybar<-barplot(y.means, # produce a simple bar chart of the means called 'mybar'
               names.arg=c(1:7), # set x-axis tick labels
               ylab="Mean birth weight (lb)", xlab="Litter", # set axes labels
               col="cadetblue1", # choose colour
               ylim=c(0,5),yaxs = "i") # set y-axis limits
abline(h=0) # add an x-axis line

y.sd <-  c(sd(lit1),sd(lit2),sd(lit3),sd(lit4),
           sd(lit5),sd(lit6),sd(lit7)) # create a list of the standard deviations

arrows(mybar, # use mybar x-coordinates for x0
       y.means-y.sd, # set y0 lowest error bar limits
       mybar, # use mybar x-coordinates for x1
       y.means+y.sd, # set y1 highest error bar limits
       length=0.1, # size of arrowhead in inches
       angle=90, # angle of arrowhead relative to shaft
       code=3) # use error bar type arrowheads

text(mybar, 4.5, # set x- and y-coordinates for text
     labels=paste("n=",n,sep=""), # add list of sample sizes as text
     col="royalblue3") # choose text colour

### Refining our bar chart: adding sample sizes as text at various heights

altbar<-barplot(y.means, names.arg=c(1:7), 
                ylab="Mean birth weight (lb)", xlab="Litter", 
                col="cadetblue1", ylim=c(0,5),yaxs = "i") # draw base bar chart with error bars again,
                                                          # calling it 'altbar'
abline(h=0) # add an x-axis line
arrows(altbar, y.means-y.sd, altbar, y.means+y.sd, 
       length=0.1, angle=90, code=3) # add error bars

text(altbar, y.means+y.sd+1.0, # add sample sizes 1 above bars+standard deviations
     labels=paste("n=",n,sep=""), # add list of sample sizes as text
     col="royalblue3")

###########################################################################

##   GROUPED BAR CHART   ##
## WITH SEED GERMINATION ##
##         DATA          ##

seeds <- read.table(file.choose(), header = T, sep = ",") # load the data in from the Excel CSV. file
                                                          # 'seed_data' and name it 'seeds'
View(seeds) # view the data

uncovered <- subset(seeds, treatment=="uncovered") # subset 'uncovered' data
covered <- subset(seeds, treatment=="covered") # subset 'covered' data

tapply(uncovered$germinated, uncovered$water, mean) # return the means for each treatment
                                                    # for the uncovered sample

uncovmeans <- c(24.25, 46.00, 66.75, 78.00, 72.75) # create a list of the mean 
                                                   # germination values
                        
tapply(uncovered$germinated, uncovered$water, sd) # return the standard deviations 
                                                  # for each treatment for the 
                                                  # uncovered sample

uncovsd <- c(2.217356, 9.273618, 11.586630,  5.830952,  4.856267) # create a list of 
                                                                  # the standard deviations

tapply(covered$germinated, covered$water, mean) # return the means for each treatment
                                                # for the covered sample

covmeans <- c(42.75000, 75.25000, 76.25000, 52.00000, 37.33333) # create a list of the mean 
                                                                # germination values

tapply(covered$germinated, covered$water, sd) # return the standard deviations 
                                              # for each treatment for the covered sample

covsd <- c(1.707825, 6.946222, 3.304038, 9.201449, 7.094599) # create a list of 
                                                             # the standard deviations

seedmeans <- rbind(uncovmeans,covmeans) # join the means from the two subsetted 
                                        # treatments into one dataset

seedmeans # look at the dataset

seedbar<-barplot(seedmeans, beside=TRUE, # produce a simple grouped bar chart
                                         # of the means called 'seedbar'
                 names.arg=c(1:5), # set x-axis tick labels
                 ylab="Seeds germinating per box", xlab="Amount of water", # set axes labels
                 col=c("gray","yellow"), # choose colours
                 ylim=c(0,100),yaxs = "i") # set y-axis limits
abline(h=0) # add an x-axis line

seedsd <- rbind(uncovsd,covsd) # join the standard deviations from the two 
                               # subsetted treatments into one dataset

arrows(seedbar, seedmeans-seedsd, seedbar, seedmeans+seedsd, # add error bars
       length=0.06, # select shorter length of arrowhead
       angle=90, code=3) 

legend("topleft", # add legend to the top left of plot
       bty="n", # do not encase legend in a box
       legend=c("Uncovered","Covered"), # set names of samples
       fill=c("gray","yellow")) # set colours of samples

