#1. download R; download and open R studio
# some basic things:
# R can work just like a calculator
2+2
# R can also name and save "objects" which can be used later on
four <- 2+2
four
six <- 2 + four
six
#####################################
#2. import data - working directories, etc...
# set your working directory using the setwd() function
setwd("C:/Users/Patrick/Desktop/Data Expeditions/")
# you can also click in the toolbar above to Session --> Set Working Directory --> Choose Directory

# load your data using the read.csv() function. 
# loading your data requires you to name it as an object. Here we use an arrow to tell R to name our loaded data "data"
data <- read.csv("Dataset for DDE.csv")
# if your data does not load, make sure it is in your current working directory and is a .csv format
#####################################
#3. inspect data - str() function tells you the structure of the data
str(data)
# the head() function lets you look at the first 6 rows of your data
head(data)
#####################################
#4. subset and index data
# R organizes data using indexing, which works using square brackets [] and is organized as: [row, column]
# so, to see what is in the first row and third column of your dataset:
data[1,3]
# if you want all the entries (columns) in one row, you can leave the columns field of your indexing blank:
data[1,]
# if you want to see only the data that match certain criteria, you can also use indexing
# in this case, indexing works both using brackets and the dollar sign. The dollar sign subsets your dataframe (data) to a specific column
# say you want to see only the individuals where the species name is Ancylomenes pedersoni
Ap<-data[data$Species=="Ancylomenes_pedersoni",]
#note that we saved this as a new object, Ap
str(Ap)
# what if you want only Ancylomenes pedersoni, AND only individual #3?
Ap3<-data[data$Species=="Ancylomenes_pedersoni"&data$Individual==3,]
str(Ap3)
## spend more time fiddling with indexing, because this is really important!
# e.g., what line of code would you use to get only Lysmata amboinensis? Lysmata amboinensis individual #2?
# quiz: pull out the data for Lysmata amboinensis individual #4, with wavelength value = 380?
#####################################
#5. plot data for one individual of one species - par() and other functions
# plotting the data lets you get a basic idea of what it looks like
# we can plot wavelength of light as the x-axis, and response as the y-axis
plot(x=Ap3$Wavelength,y=Ap3$Response)
# we can also make this look prettier - we can make the dots any shape or color we want, change the axes labels, etc...
# the par() function lets us do all of this
# use R studio to get help to learn about the par() function
?par()
# par() tells us that we can change the color of the plot using col
# we do this by adding a "col" argument into our plot() function
plot(x=Ap3$Wavelength,y=Ap3$Response,col='red')
#we can also change the type of plot; for example, showing both points and lines connecting points
plot(x=Ap3$Wavelength,y=Ap3$Response,col='red',type='b')
# we can change the point type using the "pch" argument - ask R for help to choose your favorite point type
?pch()
plot(x=Ap3$Wavelength,y=Ap3$Response,col='red',pch=18,type='b')
#add whatever you want! mess with par() more - change axes labels, etc....
#let's change the names of the x- and y-axes
plot(x=Ap3$Wavelength,y=Ap3$Response,col='red',pch=18,type='b',xlab='Wavelength',ylab='Response')
#this is great so far!
#but, there is lots of noise in this dataset.
#fit a smoothed curve to the data using the lowess() function instead of x and y
points(lowess(x=Ap3$Wavelength,y=Ap3$Response),type='l',col='red',lwd=2)
#we use lwd=2 to make the line width larger (easier to see)
# now we have a cleaner way of visualizing the response for this individual
#####################################
#7. we can use the points() function to add data from more individuals to the same plot
# set up a new plot, with just the smoothed curve from individual 3
plot(lowess(x=Ap3$Wavelength,y=Ap3$Response),type='l',col='red',lwd=2,xlab='Wavelength',ylab='Response')
# instead of making a new dataset for each new individual, we can use indexing to accomplish this task!
points(lowess(x=Ap$Wavelength[Ap$Individual==1],y=Ap$Response[Ap$Individual==1]),col='blue',type='l',lwd=2)
points(lowess(x=Ap$Wavelength[Ap$Individual==2],y=Ap$Response[Ap$Individual==2]),col='dark green',type='l',lwd=2)
points(lowess(x=Ap$Wavelength[Ap$Individual==4],y=Ap$Response[Ap$Individual==4]),col='orange',type='l',lwd=2)
points(lowess(x=Ap$Wavelength[Ap$Individual==5],y=Ap$Response[Ap$Individual==5]),col='purple',type='l',lwd=2)
#OH NO! our original plot didn't have the right axes!!!
#go back and change the axes limits to make sure the y-axis goes all the way to 1
plot(lowess(x=Ap3$Wavelength,y=Ap3$Response),type='l',col='red',lwd=2,xlab='Wavelength',ylab='Response',xlim=c(300,700),ylim=c(0,1.0))
#no go up and re-run the code (lines 75-78) to add the other individuals
#finally, add the last individual in our dataset 
points(lowess(x=Ap$Wavelength[Ap$Individual==6],y=Ap$Response[Ap$Individual==6]),col='black',lwd=2,type='l')
#####################################
#8. - average curve for a species; plot all three species together
# we want to create an average plot for each species, and then plot those three lines together
# the aggregate function can pull out the mean value of all individuals for each species at a given wavelength
speciesavg<-aggregate(x=data$Response,by=list(data$Species,data$Wavelength),FUN=mean)
head(speciesavg)
#rename the columns
colnames(speciesavg)<-c('species','wavelength','response')
head(speciesavg)
#now, let's make a separate dataset for each species
avgAp<-speciesavg[speciesavg$species=="Ancylomenes_pedersoni",]
avgLa<-speciesavg[speciesavg$species=="Lysmata_amboinensis",]
avgUa<-speciesavg[speciesavg$species=="Urocaridella_antonbruunii",]
#now, let's plot the lowess curve for each species
plot(lowess(x=avgAp$wavelength,y=avgAp$response),type='l',col='red',lwd=2,xlab='Wavelength',ylab='Response',xlim=c(300,700),ylim=c(0,1))
points(lowess(x=avgLa$wavelength,y=avgLa$response),type='l',col='purple',lwd=2,lty=2)
points(lowess(x=avgUa$wavelength,y=avgUa$response),type='l',col='blue',lwd=2,lty=3)
#let's make a legend so we can tell who is who
legend("topleft",c('Ancylomenes pedersoni','Lysmata amboinensis','Urocaridella antonbruunii'),col=c('red','purple','blue'),lty=c(1,2,3),lwd=c(2,2,2))
#####################################
#9. let's pull numbers out of these data - the mean and standard deviation of the peak sensitivity for each species
#use the which() function to find the row for any case where response = 1 (maximum response)
maxrow<-which(data$Response==1)
maxdata<-data[maxrow,]
#within a species, get the average a standard deviation across wavelengths with max responses
Apmean<-mean(maxdata$Wavelength[maxdata$Species=="Ancylomenes_pedersoni"])
Apsd<-sd(maxdata$Wavelength[maxdata$Species=="Ancylomenes_pedersoni"])
Lamean<-mean(maxdata$Wavelength[maxdata$Species=="Lysmata_amboinensis"])
Lasd<-sd(maxdata$Wavelength[maxdata$Species=="Lysmata_amboinensis"])
Uamean<-mean(maxdata$Wavelength[maxdata$Species=="Urocaridella_antonbruunii"])
Uasd<-sd(maxdata$Wavelength[maxdata$Species=="Urocaridella_antonbruunii"])
#make a dataframe of these data!
summarydata<-cbind(c(Apmean,Apsd),c(Lamean,Lasd),c(Uamean,Uasd))
colnames(summarydata)<-c("A. pedersoni","L. amboinensis","U. antonbruunii")
rownames(summarydata)<-c("mean","sd")
summarydata
#export the data as a csv file
write.csv(summarydata,file="C:/Users/Patrick/Desktop/Data Expeditions/summaryERGdata.csv")