setwd("C:/Users/scheld/Desktop/Derelict Gear/Data Analysis/Final/W&M Digital Archive")

## Read in data ##
DDG_Data<-read.csv("DDG_Data.csv")
#^ Column definitions:"WB"=Water body; "YEAR"=Year; "LBS"=LBS of blue crab harvest; "POTS"=Number of pots; "CRAB"=Estimated abundance (millions of individuals); "RMVS"=Derelict pot removals

## Create dataset for translog model; 1 unit added to all removal observations ##
TLM_Data<-cbind(DDG_Data[,c("WB","YEAR")],log(DDG_Data[,c("LBS","POTS","CRAB")]),log(DDG_Data[,"RMVS"]+1),
                log(DDG_Data[,"POTS"])*log(DDG_Data[,"POTS"]),log(DDG_Data[,"CRAB"])*log(DDG_Data[,"CRAB"]),
                log(DDG_Data[,"RMVS"]+1)*log(DDG_Data[,"RMVS"]+1),log(DDG_Data[,"POTS"])*log(DDG_Data[,"RMVS"]+1),
                log(DDG_Data[,"POTS"])*log(DDG_Data[,"CRAB"]),log(DDG_Data[,"RMVS"]+1)*log(DDG_Data[,"CRAB"]),
                ifelse(DDG_Data$YEAR>=2009,1,0))
colnames(TLM_Data)<-c("WB","YEAR","LNLBS","LNPOTS","LNCRAB","LNRMVS","LNPOTS2","LNCRAB2","LNRMVS2",
                      "LNPOTS_LNRMVS","LNPOTS_LNCRAB","LNRMVS_LNCRAB","DUM_0914")
  
## Load lme4 package and fit random effects translog model ##
require(lme4)
TLM<-lmer(LNLBS~LNPOTS+LNPOTS2+LNCRAB+LNCRAB2+LNRMVS+LNRMVS2+LNPOTS_LNCRAB+LNPOTS_LNRMVS+LNRMVS_LNCRAB+DUM_0914+(1|WB),data=TLM_Data,REML=F)

## Semi-parametric (residual) bootstrap of model parameters ##
ParBS<-function(x){return(fixef(x))}
Parameters<-bootMer(TLM,ParBS,10000,parallel="multicore",ncpus=4,type="semiparametric",use.u=T)

# Harvest predictions with and without removals #
YearIndex<-which(DDG_Data$YEAR>=2009)
XMatA<-model.matrix(TLM)
XMatCF<-model.matrix(TLM)
XMatCF[,c("LNRMVS","LNRMVS2","LNPOTS_LNRMVS","LNRMVS_LNCRAB")]<-0
Actual<-exp(XMatA%*%t(Parameters$t)+ranef(TLM)$WB[match(TLM_Data$WB,rownames(ranef(TLM)$WB)),])
Counterfactual<-exp(XMatCF%*%t(Parameters$t)+ranef(TLM)$WB[match(TLM_Data$WB,rownames(ranef(TLM)$WB)),])
Effects<-Actual-Counterfactual

# Mean and sd of total harvest effects #
mean(apply(Effects[YearIndex,],2,sum))
sd(apply(Effects[YearIndex,],2,sum))


##############################################################
## Additional code available upon request (scheld@vims.edu) ##
##############################################################