setwd("~/WILD4580/data/exercise_dat") # set working directory
leks<-read.csv("gsg_leks.csv") #import dataset

AnnualLekSummary <- function(lek_data, year) {
  leks <- read.csv(lek_data) # read in lek data
  # create empty array to store data by year
  leks_year = array(0, dim=c(length(which(leks$year==year)),5))
  # set column names
  dimnames(leks_year)[[2]]<-list("year", "disturbance", "lek_id", "complex", "tot_male")
  # convert output array to data frame
  leks_year<- data.frame(leks_year)
  j = 0
  # loop through all rows and pull write those of the given year to the output df
  for (i in 1:length(leks$year)){
      if (leks$year[i]==year){ j=j+1; leks_year[j,]<-leks[i,]} }
  # calculate mean by complex in given year and specify column names
  mean <- aggregate(leks_year$tot_male,by=list(leks_year$complex),FUN=mean)
  names(mean)<-c("complex","mean_tot_male"); mean
  # calculate standard deviation and two times sd by complex in given year
  sd <- aggregate(leks_year$tot_male,by=list(leks_year$complex),FUN=sd)
  sd2 <- sd*2
  names(sd2) <- c("complex","2sd_tot_male"); sd2
  # specify n for each complex in given year
  n = aggregate(leks_year$tot_male, by=list(leks_year$complex), FUN=length)
  # calculate 90% confidence intervals per complex in given year
  lwr_CI = array(0, dim=c(length(mean$complex),2))
  lwr_CI[,1] <- mean[,1]
  dimnames(lwr_CI)[[2]]<-list("complex", "lwr_90%_CI")
  for (i in 1:length(mean$complex)){
    lwr_CI[i,2] <- mean[i,2] - qnorm(0.95)*sd[i,2]/sqrt(n[i,2]) }
  upr_CI = array(0, dim=c(length(mean$complex),2))
  upr_CI[,1] <- mean[,1]
  dimnames(upr_CI)[[2]]<-list("complex", "upr_90%_CI")
  for (i in 1:length(mean$complex)){
    upr_CI[i,2] <- mean[i,2] + qnorm(0.95)*sd[i,2]/sqrt(n[i,2]) }
  # Merge these into single data object, and export as a .csv file
  lek_stats = cbind(mean, sd2[,2], n[,2], lwr_CI[,2], upr_CI[,2])
  # set column names
  names(lek_stats)<-c("complex", "mean_tot_male", "2*sd", "n", "lwr_90_CI", "upr_90_CI")
  # specify output filename by year
  fn = paste(year, "lek_statistics.csv", sep="_")
  # export statistics data frame as .csv file
  write.csv(lek_stats, file = fn)}

AnnualLekSummary(lek_data = "C:/Users/Maggie/Documents/WILD4580/data/exercise_dat/gsg_leks.csv", year = 2007); read.csv("2007_lek_statistics.csv")

AnnualLekSummary(lek_data = "C:/Users/Maggie/Documents/WILD4580/data/exercise_dat/gsg_leks.csv", year = 1996); read.csv("1996_lek_statistics.csv")