#Import and analyze data from paper

library(boot)
library(MASS)
rm(list=ls())


#-------------------------LOAD DATA--------------------------
load('data.RData')

#----------------------DEFINE FUNCTION-----------------------
mean2<-function(x,i) mean(x[i])

#-----------FLUCTUATING ENVIRONMENTS EXPERIMENT---------------
attach(data.fluctuating.environments)

mean=mean(delta.t.g); mean
stdErr=sd(delta.t.g)/sqrt(length(delta.t.g)); stdErr
mean(delta.t.q)
stdErr=sd(delta.t.q)/sqrt(length(delta.t.q)); stdErr

x<-delta.t.g-delta.t.q;
test.pooled<-boot(x,statistic=mean2,100000,stype='i')
mean.diff<-test.pooled$t0
q<-quantile(test.pooled$t,c(0.025,0.975));q
e<-ecdf(test.pooled$t)
p<-e(0);p

x<-delta.t.g[treatment==0]-delta.t.q[treatment==0];
test.low<-boot(x,statistic=mean2,100000,stype='i')
mean.diff<-test.low$t0
q<-quantile(test.low$t,c(0.025,0.975));q
e<-ecdf(test.low$t)
p<-e(0);p

x<-delta.t.g[treatment==1]-delta.t.q[treatment==1];
test.med<-boot(x,statistic=mean2,100000,stype='i')
mean.diff<-test.med$t0
q<-quantile(test.med$t,c(0.025,0.975));q
e<-ecdf(test.med$t)
p<-e(0);p

x<-delta.t.g[treatment==2]-delta.t.q[treatment==2];
test.high<-boot(x,statistic=mean2,100000,stype='i')
mean.diff<-test.high$t0
q<-quantile(test.high$t,c(0.025,0.975));q
e<-ecdf(test.high$t)
p<-e(0);p

x11(height=4,width=5); par(mfrow=c(1,2))
truehist(delta.t.g144,h=2,col='grey',xlab='duration of growth (days)',xlim=c(0,20))
truehist(delta.t.q,h=2,col='grey',xlab='duration of decline (days)',xlim=c(0,20))


#-----------------HABITAT SIZE EXPERIMENT--------------------
attach(data.habitat.size)

mean=mean(delta.t.g); mean
stdErr=sd(delta.t.g)/sqrt(length(delta.t.g)); stdErr
mean(delta.t.q)
stdErr=sd(delta.t.q)/sqrt(length(delta.t.q)); stdErr

x<-delta.t.g-delta.t.q
test<-boot(x,statistic=mean2,100000,stype='i')
mean.diff<-test$t0
q<-quantile(test$t,c(0.025,0.975));q
e<-ecdf(test$t)
p<-e(0)

model<-aov(log(delta.t.q+1)~log(lambda+1)+size+food); summary(model)

x11(height=4,width=5); par(mfrow=c(1,2))
truehist(delta.t.g,h=2,col='grey',xlab='duration of growth (wks)',xlim=c(0,25))
truehist(delta.t.q,h=2,col='grey',xlab='duration of decline (wks)',xlim=c(0,25))


x11(height=4,width=4); plot(log(lambda[size=='Large']+1),log(delta.t.q[size=='Large']+1),
	xlab=expression(log(lambda+1)),ylab='log(E+1)')
points(log(lambda[size=='Small']+1),log(delta.t.q[size=='Small']+1),pch=19)