###################################################################################################################
# This R script computes correlations between tree radial growth and monthly lagged climate for a variety of climate
# variables, then summarizes the tree-level correlations by species and site. 
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
setwd('')
require(lattice)
require(SPEI)
require(plyr)
require(SDMTools) # for wt.mean and wt.sd

#---------------------------------------------------------------------------------------
# set time periods and variable 
#---------------------------------------------------------------------------------------
strata <- c('wj','pp','gf')
plot <- 1:5
lags <- 1:36
n.lags <- length(lags)
yoi <- 1964:2013
n.yoi <- length(yoi)
yoi.lags <- (min(yoi)-n.lags/12):max(yoi)
n.yoi.lags <- length(yoi.lags)
sp.oi <- c('juoc','pipo','abgr')

#---------------------------------------------------------------------------------------
# Read in RW, climate, and biomass data
#---------------------------------------------------------------------------------------
# read in the master ring width index data
rw <- read.table('data/treerings/treering_master_data.csv', sep = ',', header = T)
rw <- rw[rw$species%in%sp.oi,]
rw <- subset(rw, reliable <= min(yoi)) # take only trees that have atleast 50 years of reliable rw mmnts
rw <- subset(rw, year >=min(yoi))

# read in the prism climate data (tmean, ppt, cmi, pet)
clim <- read.table('output/site_prism_data.csv', sep = ',', header = T)
clim$strata <- tolower(clim$strata)
clim$plot.id <- tolower(clim$plot.id)
clim <- subset(clim, year >= min(yoi.lags)-2)
clim$month <- rep(1:12, nrow(clim)/12)

# clim.pars <- unique(clim$clim.par)
# clim.funs.df <- data.frame(clim.par=clim.pars, funs=c('mean','sum','sum','sum'))
#vpd <- subset(clim, clim.par == 'vpd')
clim.pars <- c('vpd','freeze','cmi','defecit','spei')
clim.funs <- c('mean','sum','sum','sum','NA')
clim.funs.df <- data.frame(clim.par=clim.pars, funs=clim.funs)

# biomass for each species in each subplot; used for biomass-weighted plot mean rw-clim correlations
treeC.by.sp.subplot <- read.csv('data/subplot_sp_treeC.csv', header=T)
treeC.by.sp.subplot$subplot.sp.id <- tolower(paste(treeC.by.sp.subplot$strata, treeC.by.sp.subplot$plot, treeC.by.sp.subplot$subplot, treeC.by.sp.subplot$species, sep='.'))
treeC.by.sp.subplot$subplot.sp.id <- factor(treeC.by.sp.subplot$subplot.sp.id, levels = unique(treeC.by.sp.subplot$subplot.sp.id))

# clim.pars <- c('spei')
# clim.funs.df <- data.frame(clim.par=clim.pars, funs=c('NA'))

#---------------------------------------------------------------------------------------
# Correlate RWI against climate data lagged to varying extents
#---------------------------------------------------------------------------------------
# 
# i='wj'
# j=1
# k='vpd'
# l=1
# m=1
# n=1
# p=1
# strata <- c('pp')
# plot <- 3:5
# clim.pars <- c('vpd')
# lags <- 1:12
# n.lags <- length(lags)

for (i in strata){
  for (j in plot){
    # SKIP WJ-4 AND PP-4
    if (i=='wj' & j == 4){
      next()
    } else if (i=='pp' & j==4){
      next()
    } else if (i=='gf' & j==3){
      next()
    }
    
    # get the rw data for the plot
    plot.rw.data <- subset(rw, strata==i & plot==j)
    plot.rw.series <- unique(plot.rw.data$id) 
    if (length(plot.rw.series)==0){
      next()
    }
    plot.n.rw.series <- length(plot.rw.series)
    for (k in clim.pars){
      #------------------------------------------------------------
      # build data frame to store correlation results
      rw.clim.cor.df <- as.data.frame(matrix(NA, nrow=0, ncol=10))
      colnames(rw.clim.cor.df) <- c('strata','plot','tree.id','species','clim.par','month','lag','cor.r','cor.p','cor.cat')
      
      #----------------------------------------------------------
      # calculate RW-SPEI correlations 
      if (k=='spei'){
        clim.par.allmonths <- subset(clim, strata == i & plot == j & clim.par=='cmi')
        for (m in lags){
          spei.lag <- spei(clim.par.allmonths$values,m)  
          clim.par.allmonths$spei <- as.vector(spei.lag$fitted)
          for (l in 1:12){
            spei.month <- subset(clim.par.allmonths, month==l & year >= min(yoi)-1)
            # detrend and pre-whiten spei data
            spei.lm <- lm(spei~year, spei.month)
            spei.month$spei.dt <- mean(spei.month$spei)+as.vector(spei.lm$residuals)
            #spei.month$spei.ar <- ar(spei.month$spei.dt)$resid
            spei.month <- spei.month[-1,]
            for (p in 1:plot.n.rw.series){
              rwi.series <- subset(plot.rw.data, id==plot.rw.series[p])
              cor.results <- cor.test(rev(rwi.series$rwi.ar),na.omit(spei.month$spei)) # reverse order of the RW data
              # categorize rw-spei correlations (pos/ns/neg)
              cor.cat <- if(cor.results$p.value<0.05){
                if (cor.results$estimate > 0){
                  'pos'
                } else if (cor.results$estimate < 0){
                  'neg'
                }} else {'ns'}
              
              # store correlation results in data frame
              rw.clim.cor.df <- rbind(rw.clim.cor.df, 
                                      data.frame(strata=i,plot=j,tree.id=plot.rw.series[p],species=unique(rwi.series$species),clim.par=k,
                                                 month=l,lag=m,cor.r=round(as.numeric(cor.results$estimate),3),
                                                 cor.p=round(as.numeric(cor.results$p.value),3),cor.cat=cor.cat))
            }
          }
        }
        outname <- paste('output/rw_clim_relations/tmp/', i,'.',j,'.',k,'.climcor.csv', sep='')
        write.table(rw.clim.cor.df, outname, sep=',', row.names = F, col.names = T)
      }
      #----------------------------------------------------------
      # calcute other RW-climate correlations
      else {
        clim.par.allmonths <- subset(clim, strata == i & plot == j & clim.par==k)
        clim.par.fun <- as.character(clim.funs.df$funs[which(clim.funs.df$clim.par%in%k)])
        for (l in 1:12){
          for (m in lags){
            # build data frame to temporarly hold the lagged climate data 
            clim.par.lag <- data.frame(matrix(NA, (n.yoi+1)*m, ncol=3)) # take one extra year to account for year lost after removing autocor.
            colnames(clim.par.lag) <- c('year','lag','value')
            clim.par.lag$year <- sort(rep(c(min(yoi)-1,yoi), m))
            clim.par.lag$lag <- rep(lags[m], n.yoi+1)
            clim.par.allmonths.lags <- subset(clim.par.allmonths, year>=min(yoi)-(1+ceiling((m+1)/12)))
            endpt.indx <- which(clim.par.allmonths.lags$month%in%l) # ID the endpoint indices for each year
            endpt.indx <- endpt.indx[-c(1:ceiling((m+1)/12))] # Trim the first m indices from start of series
            lag.indx <- c()
            for (n in 1:length(endpt.indx)){
              lag.indx <- c(lag.indx, seq((endpt.indx[n]-m+1), endpt.indx[n], 1))  
            }
            clim.par.lag$value <- clim.par.allmonths.lags$values[lag.indx]
            clim.par.agg <- aggregate(value~year, clim.par.lag, FUN = clim.par.fun) #aggregate the climate data across lags
            # detrend, anomalized, and prewhiten
            clim.par.agg$zscore <- with(clim.par.agg, (value-mean(value, na.rm = T))/sd(value, na.rm = T))
            if (all(is.na(clim.par.agg$zscore))){ # to account for months where there are no freezing days
              clim.par.agg$zscore <- 0
              clim.par.agg$clim.dt <- 0
              clim.par.agg$clim.ar <- 0
            } else {
                clim.lm <- lm(zscore~year, clim.par.agg)
                clim.par.agg$clim.dt <- mean(clim.par.agg$zscore)+as.vector(clim.lm$residuals)
                clim.par.agg$clim.ar <- ar(clim.par.agg$clim.dt)$resid
            }  
            clim.par.agg <- subset(clim.par.agg, year>=min(yoi))
            for (p in 1:plot.n.rw.series){
              rwi.series <- subset(plot.rw.data, id==plot.rw.series[p])
              cor.results <- cor.test(rev(rwi.series$rwi.ar),clim.par.agg$clim.ar) # reverse order of the RW data to match clim data order (oldest to newest)
              if (all(clim.par.agg$clim.ar==0)){ # correlation won't run for summer freeze days since all = 0
                cor.results$p.value <- 1
                cor.results$estimate <- 0
              }
              
              cor.cat <- if(cor.results$p.value<0.05){
                if (cor.results$estimate > 0){
                  'pos'
                } else if (cor.results$estimate < 0){
                  'neg'
                }} else {'ns'}
  
              rw.clim.cor.df <- rbind(rw.clim.cor.df, data.frame(strata=i,plot=j,tree.id=plot.rw.series[p],species=unique(rwi.series$species),
                                                                 clim.par=k,month=l,lag=m,cor.r=round(as.numeric(cor.results$estimate),3),
                                                                cor.p=round(as.numeric(cor.results$p.value),3),cor.cat=cor.cat))
            }
          }
        }
      }
      outname <- paste('output/rw_clim_relations/tmp/', i,'.',j,'.',k,'.climcor.csv', sep='')
      write.table(rw.clim.cor.df, outname, sep=',', row.names = F, col.names = T)
      print(paste('finished:',i,j,k, sep=' '))
    }
  }
}

#----------------------------------------------------------------------
# Compile the individual files into one large data frame
#----------------------------------------------------------------------
f <- list.files('output/rw_clim_relations/tmp/')
for (i in 1:length(f)){
  if (i==1){
    inname <- paste('output/rw_clim_relations/tmp/', f[i], sep = '')
    rw.clim.cor.df <- read.table(inname, sep = ',', header = T)
  } else {
    inname <- paste('output/rw_clim_relations/tmp/', f[i], sep = '')
    df.new <- read.table(inname, sep = ',', header = T)[,1:10]
    rw.clim.cor.df <- rbind(rw.clim.cor.df, df.new)
  }
}

rw.clim.cor.df$subplot <- factor(substr(rw.clim.cor.df$tree.id, 6,6))
rw.clim.cor.df$subplot.sp.id <- paste(rw.clim.cor.df$strata, rw.clim.cor.df$plot, rw.clim.cor.df$subplot, rw.clim.cor.df$species, sep='.')

# Add spcies biomass weight to data frames
rw.clim.cor.df$sp.frac <- rep(NA, nrow(rw.clim.cor.df))
subplot.sp.id <- unique(rw.clim.cor.df$subplot.sp.id)
for (i in subplot.sp.id){
  sp.frac <- subset(treeC.by.sp.subplot, subplot.sp.id==i)
  index <- which(rw.clim.cor.df$subplot.sp.id%in%i)
  rw.clim.cor.df$sp.frac[index] <- rep(sp.frac$sp.frac, length(index)) 
}

write.table(rw.clim.cor.df, 'output/rw_clim_relations/rw_clim_cors_all_species.csv', sep = ',', row.names = F, col.names = T)
head(rw.clim.cor.df)
dim(rw.clim.cor.df)

#----------------------------------------------------------------------
# Mean RW-climate correlation by species and plot
#----------------------------------------------------------------------
#-------------- by species --------------
rw.clim.species.avg <- ddply(rw.clim.cor.df, .(species,month,lag,clim.par), plyr::summarize, 
                             species=unique(species), n.trees=length(tree.id), cor.r.avg=mean(cor.r), cor.r.sd=sd(cor.r),  cor.p.avg=mean(cor.p))
rw.clim.species.avg$cor.r.se <- with(rw.clim.species.avg, cor.r.sd/sqrt(n.trees))
write.table(rw.clim.species.avg, 'output/rw_clim_relations/rw_clim_sp_avg_cors.csv', sep = ',', row.names = F, col.names = T)

# summarize top 5% of strongest rw-clim relations for each climate parameter
rw.clim.sp.top.cors <- as.data.frame(matrix(NA, nrow=0, ncol=11))
colnames(rw.clim.sp.top.cors) <- c('clim.par','species','n.trees','month.avg','month.sd','lag.avg','lag.sd','cor.avg','cor.sd','cor.se','cor.p')

# i='tmean'
#  j='juoc'
# i='defecit'
for (i in clim.pars){
  for (j in sp.oi){
    sp.dat <- subset(rw.clim.species.avg, clim.par==i & species == j)
    cutoff <- quantile(abs(sp.dat$cor.r.avg), 0.95)
    sp.dat.top <- sp.dat[which(abs(sp.dat$cor.r.avg)>cutoff),]
    sp.dat.top.summary <- ddply(sp.dat.top, .(species,clim.par), plyr::summarize, species=unique(species), n.tree=unique(n.trees), 
                                month.avg=round(mean(month),1), month.sd=round(sd(month),1), 
                                lag.avg=round(mean(lag),1), lag.sd=round(sd(month),1), 
                                cor.avg=round(mean(cor.r.avg),2), cor.sd=round(mean(cor.r.sd),2), cor.se=round(mean(cor.r.se),2),
                                cor.p=round(mean(cor.p.avg),3))
    rw.clim.sp.top.cors <- rbind(rw.clim.sp.top.cors, sp.dat.top.summary)
  }
}

write.table(rw.clim.sp.top.cors, 'output/rw_clim_relations/rw_clim_sp_top_cors.csv', sep = ',', row.names = F, col.names = T)

#-------------- by plot, with biomass-weighted averaging by species in each subplot --------------

rw.clim.subplot.wavg <- ddply(rw.clim.cor.df, .(strata,plot,subplot,month,lag,clim.par), plyr::summarize, 
                              n.trees=length(cor.r), cor.r.wavg=wt.mean(cor.r, sp.frac), cor.r.wsd=wt.sd(cor.r,sp.frac), cor.p.wavg=wt.mean(cor.p,sp.frac))
rw.clim.plot.avg <- ddply(rw.clim.subplot.wavg, .(strata,plot,month,lag,clim.par), plyr::summarize, cor.r.avg=mean(cor.r.wavg), cor.r.sd=sd(cor.r.wavg),cor.p.avg=mean(cor.p.wavg))
rw.clim.plot.avg$cor.r.se <- with(rw.clim.plot.avg, cor.r.sd/sqrt(4))
write.table(rw.clim.plot.avg, 'output/rw_clim_relations/rw_clim_plot_avg_cors.csv', sep = ',', row.names = F, col.names = T)

rw.clim.plot.top.cors <- as.data.frame(matrix(NA, nrow=0, ncol=12))
colnames(rw.clim.plot.top.cors) <- c('clim.par','strata','plot','n.trees','month.avg','month.sd','lag.avg','lag.sd','cor.avg','cor.sd','cor.se','cor.p')

for (i in clim.pars){
  for (j in strata){
    for (k in plot){
      if (i=='wj' & j == 4){
        next()
      } else if (i=='pp' & j==4){
        next()
      } else if (i=='gf' & j==3){
        next()
      }
      plot.dat <- subset(rw.clim.plot.avg, clim.par==i & strata == j & plot == k)
      cutoff <- quantile(abs(plot.dat$cor.r.avg), 0.95)
      plot.dat.top <- plot.dat[which(abs(plot.dat$cor.r.avg)>cutoff),]
      plot.dat.top.summary <- ddply(plot.dat.top, .(clim.par,strata,plot), plyr::summarize, strata=unique(strata), plot=unique(plot),
                                    month.avg=round(mean(month),1), month.sd=round(sd(month),1), 
                                    lag.avg=round(mean(lag),1), lag.sd=round(sd(month),1), 
                                    cor.avg=round(mean(cor.r.avg),2), cor.sd=round(mean(cor.r.sd),2), cor.se=round(mean(cor.r.se),2),
                                    cor.p=round(mean(cor.p.avg),3))
      rw.clim.plot.top.cors <- rbind(rw.clim.plot.top.cors, plot.dat.top.summary)
    }
  }
}
write.table(rw.clim.plot.top.cors, 'output/rw_clim_relations/rw_clim_plot_top_cors.csv', sep = ',', row.names = F, col.names = T)

#----------------------------------------------------------------------
# Fractional response by species and plot
#----------------------------------------------------------------------

#-------------- by species -------------------
# count total number of correlations 
rw.clim.species.cor.sum <- ddply(rw.clim.cor.df, .(clim.par,species,month,lag), summarise, cor.cnt=length(cor.cat), .drop=F)
rw.clim.species.cor.sum.expanded <- rw.clim.species.cor.sum[sort(rep(seq.int(1,nrow(rw.clim.species.cor.sum)),3)),1:5] # expand df, repping each row 3 times

# count number of correlation in each response category
rw.clim.species.cor.cnt <- ddply(rw.clim.cor.df, .(clim.par,species,month,lag,cor.cat), summarise, cor.cnt=length(cor.cat), .drop=F)

# calculate fraction of trees in each category
rw.clim.species.cor.cnt$cor.frac <- rw.clim.species.cor.cnt$cor.cnt/rw.clim.species.cor.sum.expanded$cor.cnt
rw.clim.species.cor.cnt$N <- rw.clim.species.cor.sum.expanded$cor.cnt
rw.clim.species.cor.cnt$cor.frac.sd <- with(rw.clim.species.cor.cnt, sqrt((cor.frac*(1-cor.frac))/N))
rw.clim.species.cor.cnt$cor.pcnt.avg <- round(rw.clim.species.cor.cnt$cor.frac*100)
rw.clim.species.cor.cnt$cor.pcnt.sd <- round(rw.clim.species.cor.cnt$cor.frac.sd*100)

write.table(rw.clim.species.cor.cnt, 'output/rw_clim_relations/rw_clim_sp_frac_resp.csv', sep = ',', row.names = F, col.names = T)

# identify the strongest positive  population response 
rw.clim.sp.frac.resp.pos <- subset(rw.clim.species.cor.cnt, cor.cat=='pos')
rw.clim.sp.pcnt.resp.pos.max <- as.data.frame(matrix(NA, nrow=0, ncol=9))
colnames(rw.clim.sp.pcnt.resp.pos.max) <- c('clim.par','strata','month.avg','month.sd','lag.avg','lag.sd','pos.cor.pcnt.avg','pos.cor.pcnt.sd')

for (i in clim.pars){
  for (j in sp.oi){
    sp.dat <- subset(rw.clim.sp.frac.resp.pos, clim.par==i & species == j)
    cutoff <- quantile(sp.dat$cor.frac, 0.95)
    sp.dat.top <- sp.dat[which(sp.dat$cor.frac>cutoff),]
    sp.dat.top.summary <- ddply(sp.dat.top, .(species,clim.par), plyr::summarize, species=unique(species), month.avg=round(mean(month),1), month.sd=round(sd(month),1), 
                                lag.avg=round(mean(lag),1), lag.sd=round(sd(month),1), 
                                pos.cor.pcnt.avg=round(mean(cor.pcnt.avg),2), pos.cor.pcnt.sd=round(mean(cor.pcnt.sd),2))
    rw.clim.sp.pcnt.resp.pos.max <- rbind(rw.clim.sp.pcnt.resp.pos.max, sp.dat.top.summary)
  }
}

write.table(rw.clim.sp.pcnt.resp.pos.max, 'output/rw_clim_relations/rw_clim_sp_frac_pos_resp_max.csv', sep = ',', row.names = F, col.names = T)

# identify the strongest negitive  population response 
rw.clim.sp.frac.resp.neg <- subset(rw.clim.species.cor.cnt, cor.cat=='neg')
rw.clim.sp.pcnt.resp.neg.max <- as.data.frame(matrix(NA, nrow=0, ncol=9))
colnames(rw.clim.sp.pcnt.resp.neg.max) <- c('clim.par','strata','month.avg','month.sd','lag.avg','lag.sd','neg.cor.pcnt.avg','neg.cor.pcnt.sd')

for (i in clim.pars){
  for (j in sp.oi){
    sp.dat <- subset(rw.clim.sp.frac.resp.neg, clim.par==i & species == j)
    cutoff <- quantile(sp.dat$cor.frac, 0.95)
    sp.dat.top <- sp.dat[which(sp.dat$cor.frac>cutoff),]
    sp.dat.top.summary <- ddply(sp.dat.top, .(species,clim.par), plyr::summarize, species=unique(species), month.avg=round(mean(month),1), month.sd=round(sd(month),1), 
                                lag.avg=round(mean(lag),1), lag.sd=round(sd(month),1), 
                                neg.cor.pcnt.avg=round(mean(cor.pcnt.avg),2), neg.cor.pcnt.sd=round(mean(cor.pcnt.sd),2))
    rw.clim.sp.pcnt.resp.neg.max <- rbind(rw.clim.sp.pcnt.resp.neg.max, sp.dat.top.summary)
  }
}

write.table(rw.clim.sp.pcnt.resp.neg.max, 'output/rw_clim_relations/rw_clim_sp_frac_neg_resp_max.csv', sep = ',', row.names = F, col.names = T)

#------------- by plot -----------------------
rw.clim.subplot.cor.sum <- ddply(rw.clim.cor.df, .(clim.par,strata,plot,subplot,month,lag), summarise, cor.cnt=length(cor.cat), .drop=F)
rw.clim.subplot.cor.sum.expanded <- rw.clim.subplot.cor.sum[sort(rep(seq.int(1,nrow(rw.clim.subplot.cor.sum)),3)),1:7] # expand df, repping each row 3 times

# count number of correlation in each response category
rw.clim.subplot.cor.cnt <- ddply(rw.clim.cor.df, .(clim.par,strata,plot,subplot,month,lag,cor.cat), summarise, cor.cnt=length(cor.cat), .drop=F)

# calculate fraction of trees in each category
rw.clim.subplot.cor.cnt$cor.frac <- rw.clim.subplot.cor.cnt$cor.cnt/rw.clim.subplot.cor.sum.expanded$cor.cnt
rw.clim.subplot.cor.cnt$N <- rw.clim.subplot.cor.sum.expanded$cor.cnt
rw.clim.subplot.cor.cnt$cor.frac.sd <- with(rw.clim.subplot.cor.cnt, sqrt((cor.frac*(1-cor.frac))/N))
rw.clim.plot.cor.cnt <- ddply(rw.clim.subplot.cor.cnt, .(clim.par,strata,plot,month,lag,cor.cat), summarise, 
                              cor.frac.avg=mean(cor.frac), cor.frac.sd=mean(cor.frac.sd), .drop=F)
rw.clim.plot.cor.cnt$cor.pcnt.avg <- round(rw.clim.plot.cor.cnt$cor.frac.avg*100)
rw.clim.plot.cor.cnt$cor.pcnt.sd <- round(rw.clim.plot.cor.cnt$cor.frac.sd*100)
write.table(rw.clim.plot.cor.cnt, 'output/rw_clim_relations/rw_clim_plot_frac_resp.csv', sep = ',', row.names = F, col.names = T)

# identify the strongest positive population response 
rw.clim.plot.cor.cnt.pos <- subset(rw.clim.plot.cor.cnt, cor.cat=='pos')

rw.clim.plot.pcnt.pos.resp.max <- as.data.frame(matrix(NA, nrow=0, ncol=9))
colnames(rw.clim.plot.pcnt.pos.resp.max) <- c('clim.par','strata','plot','month.avg','month.sd','lag.avg','lag.sd','pos.cor.pcnt.avg','pos.cor.pcnt.sd')

for (i in clim.pars){
  for (j in strata){
    for (k in plot){
      if (j=='wj' & k == 4){
        next()
      } else if (j=='pp' & k==4){
        next()
      } else if (j=='gf' & k==3){
        next()
      }
      plot.dat <- subset(rw.clim.plot.cor.cnt.pos, clim.par==i & strata == j & plot == k)
      cutoff <- quantile(plot.dat$cor.pcnt.avg, 0.95)
      plot.dat.top <- plot.dat[which(plot.dat$cor.pcnt.avg>=cutoff),]
      plot.dat.top.summary <- ddply(plot.dat.top, .(clim.par,strata,plot), plyr::summarize, strata=unique(strata), plot=unique(plot), 
                                    month.avg=round(mean(month),1), month.sd=round(sd(month),1), 
                                    lag.avg=round(mean(lag),1), lag.sd=round(sd(month),1), 
                                    pos.cor.pcnt.avg=round(mean(cor.pcnt.avg, na.rm = T),2), pos.cor.pcnt.sd=round(mean(cor.pcnt.sd, na.rm = T),2))
                                  #.drop=F
      rw.clim.plot.pcnt.pos.resp.max <- rbind(rw.clim.plot.pcnt.pos.resp.max, plot.dat.top.summary)
    }
  }
}
rw.clim.plot.pcnt.pos.resp.max
write.table(rw.clim.plot.pcnt.pos.resp.max, 'output/rw_clim_relations/rw_clim_plot_frac_pos_resp_max.csv', sep = ',', row.names = F, col.names = T)

# identify the strongest negative population response 
rw.clim.plot.cor.cnt.neg <- subset(rw.clim.plot.cor.cnt, cor.cat=='neg')

rw.clim.plot.pcnt.neg.resp.max <- as.data.frame(matrix(NA, nrow=0, ncol=9))
colnames(rw.clim.plot.pcnt.neg.resp.max) <- c('clim.par','strata','plot','month.avg','month.sd','lag.avg','lag.sd','neg.cor.pcnt.avg','neg.cor.pcnt.sd')

for (i in clim.pars){
  for (j in strata){
    for (k in plot){
      if (j=='wj' & k == 4){
        next()
      } else if (j=='pp' & k==4){
        next()
      } else if (j=='gf' & k==3){
        next()
      }
      plot.dat <- subset(rw.clim.plot.cor.cnt.neg, clim.par==i & strata == j & plot == k)
      cutoff <- quantile(plot.dat$cor.pcnt.avg, 0.95)
      plot.dat.top <- plot.dat[which(plot.dat$cor.pcnt.avg>=cutoff),]
      plot.dat.top.summary <- ddply(plot.dat.top, .(clim.par,strata,plot), plyr::summarize, strata=unique(strata), plot=unique(plot), 
                                    month.avg=round(mean(month),1), month.sd=round(sd(month),1), 
                                    lag.avg=round(mean(lag),1), lag.sd=round(sd(month),1), 
                                    neg.cor.pcnt.avg=round(mean(cor.pcnt.avg, na.rm = T),2), neg.cor.pcnt.sd=round(mean(cor.pcnt.sd, na.rm = T),2))
      #.drop=F
      rw.clim.plot.pcnt.neg.resp.max <- rbind(rw.clim.plot.pcnt.neg.resp.max, plot.dat.top.summary)
    }
  }
}
rw.clim.plot.pcnt.neg.resp.max
write.table(rw.clim.plot.pcnt.neg.resp.max, 'output/rw_clim_relations/rw_clim_plot_frac_neg_resp_max.csv', sep = ',', row.names = F, col.names = T)



#----------------------------------------------------------------------
save.image('workspaces/rw_climate_cors.RDATA')
#load('workspaces/rw_climate_cors.RDATA')
#----------------------------------------------------------------------

