################################################################################
################################################################################
# 
# MRM2 model including Canadian studies (32 observations)
#
################################################################################
################################################################################

# Note: to run the entire script, select everything, then click on "Run"
# otherwise "Run" will execute only the selected line
  
# Preliminary housekeeping
################################################################################
# set working directory to keep file paths short
setwd('C:/Users/55338/ICF/WOTUS Reconsideration (ICF Internal Site) - General/Wetland Meta-Analysis/Wetland MRM - R Code/')

rm(list = ls(all = TRUE)) #first, clear R's workspace
options(prompt = "R> ", digits = 4) # show 4 digits by default
options(continue=" ") #remove continuation prompt
#
set.seed(37)  #sets the random number generator so we can reproduce results
#
tic<-proc.time() #start stop watch
#
# load packages we'll need
############################
library("xtable")
library("MASS")   #for draws from the multivariate normal
library("stats")  #for draws from the gamma
library("mvtnorm")#for draws from the multivariate t
library("MCMCpack") #for draws from the inverse gamma

################################################################################
#Load and Prepare Data
################################################################################
#
# read tab-delimited text data (easiest way to load)
# header=TRUE also loads variable names
#data<- read.table('./data/wetlanddataFeb2022.txt', sep="\t", header=TRUE)
data<- read.table('./data/wetlanddataJune2022_corrected.txt', sep="\t", header=TRUE)
#
# save in R-format, for future loading
#save(data, file = "./data/wetlanddataFeb2022.rda")
save(data, file = "./data/wetlanddataJune2022_corrected.rda")

# % contents of data
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# % 1   studyid         original unique study id
# % 2   obsid           original unique observation id
# % 3   lnwtp           log of wtp2021
# % 4   wlfresh         1=freshwater wetland, 0 = salt/coastal
# % 5   lnyear          log of yr of data collection - 1988
# % 6   lninc           log of inc2021
# % 7   sagulf          1 = S-Atlantic/Gulf, study involves AL,GA,SC,LA
# % 8   nema            1 = ME/mid-Atlantic, study involves DE,MD,NJ,PA,RI
# % 9   nmw             N/Mid-West, study involves KY,MI,NE,OH,WI
# % 10  CAN             1 = Canadian study
# % 11  local           f1= local scale (target pop at sub-state level)
# % 12  prov            1=provisioning function affected
# % 13  reg             1=regulating function affected
# % 14  cult            1=cultural function affected
# % 15  forest          1=forested wetland            
# % 16  q0              baseline acreage
# % 17  q1              policy acreage
# % 18  volunt          1= pay mech = voluntary contribution
# % 19  lumpsum         1=single payment
# % 20  ce              choice experiment
# % 21  nrev            1=not peer-reviewed
# % 22  median          1=orig. wtp based on econometric medians
# % 23  q0alt            alternative definition of baseline acres - see notes
# % 24  q1alt            alternative definition of policy acres - see notes

# Drop Rudd et al observations
################################
#data<-subset(data,data[,1]!=128)

# Drop Hindsley and Yoskowitz observations (remove comment to drop)
################################
#data<-subset(data,data[,1]!=130)

# Drop Nijhum thesis observations (remove comment to drop)
################################
#data<-subset(data,data[,1]!=132)

# 'attach' makes it possible to refer to the variables in the data frame by their 
# names alone, rather than as components of the data frame (e.g., height rather 
# than women$height).
attach(data)

# keep freshwater only
###############################
data<-subset(data,data[,4]==1)

# define y and X
###############################
q0<-data[,16]/1000
q1<-data[,17]/1000

# takes third column of data (lnwtp) to do log(wtp) - log(q1-q0) which = log(wtp/(q1-q0)). Log(WTP/change in acres)
y<-as.matrix(data[,3]-log(q1-q0))
n<-length(y)

# rep(1,n) creates the constant term (vector of 1's). Matrix combines this with variables in columns 1-15, (q0+q1)/2, 18-19
X<-as.matrix(cbind(rep(1,n),data[,5:15],(q0+q1)/2, data[18:19]))
k<-ncol(X)

# % contents of X
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# % 1   const           constant
# % 2   lnyear          log of yr of data collection - 1988
# % 3   lninc           log of inc2017
# % 4   sagulf          1 = S-Atlantic/Gulf, study involves AL,GA,SC,LA
# % 5   nema            1 = ME/mid-Atlantic, study involves DE,MD,NJ,PA,RI
# % 6   nmw             N/Mid-West, study involves KY,MI,NE,OH,WI
# % 8   local           f1= local scale (target pop at sub-state level)
# % 9   prov            1=provisioning function affected
# % 10   reg             1=regulating function affected
# % 11  cult            1=cultural function affected
# % 12  forest          1=forested wetland
# % 13  (q0+q1)/2
# % 14  volunt          1= pay mech = voluntary contribution
# % 15  lumpsum         1=single payment


################################################################################
# Starting values, priors, and tuners
################################################################################

r1<-100000  #burn-ins
r2<-100000  #keepers
R<-r1+r2
  
# generic OLS to get starting values
####################################

# solves for the betas or parameters using OLS
bols<-solve((t(X)) %*% X) %*% (t(X) %*% y)

# error between actual and predicted values.
e<-y-X%*%bols 

# sum of squared residuals
SSR<-(t(e)%*%e)

# SSR/(n-k)
s2<-(t(e)%*%e)/(n-k) 
  
#PRIORS:
###################
#for beta:
# initial settings for the normal distribution
mu0<-rep(0,k)
V0<-10*diag(k)

#for sig2:
# initial settings for the inverse gamma distribution (sig2 has inverse gamma and v0 has gamma). 
# v0 is parameter for weighting.
v0<-1/2
tau0<-1/2


# STARTING VALUES
######################
betadraw<-bols
sig2draw<-as.vector(s2)

################################################################################
# call Gibbs Sampler function and run sampler
################################################################################

source("./functions/gs_normal_independent.R")  #call the function

int<-gs_normal_independent(X,y,n,k,r1,r2,mu0,V0,betadraw,v0,tau0,sig2draw)
# run the sampler - you should see a flurry of counters in the console,
# in steps of 1000's

betamat<-int[[1]]
sig2mat<-int[[2]]

################################################################################
# get convergence diagnostics
################################################################################

M<-rbind(betamat,sig2mat)

source("./functions/kdiagnostics_greater0.R")  #call the function
diagnostics=kdiagnostics_greater0(M) #apply the function

################################################################################
# output table
################################################################################
# coefficients, followed by sig2
metaout_df<-data.frame(diagnostics) 
colnames(metaout_df)<-c("post. mean","post.std","p(>0)","nse","IEF","M*","CD")
rownames(metaout_df)<- c("const","lnyear","lninc","sagulf","nema","nmw", "CAN","local",
                  "prov","reg","cult","forest","(q0+q1/2)","volunt","lumpsum","sig2")

# you can now open the table in the "Global Environment" window and copy/paste it
# to Excel etc.
# Below file will write to output folder (./ represents current working directory).
#write.csv(metaout_df,file="./output using Feb22 data/MRM2_CAN_32obs.csv",row.names=TRUE)
write.csv(metaout_df,file="./output using June22 data/MRM2_CAN_32obs.csv",row.names=TRUE)

################################################################################
# save your output for further processing
################################################################################
#save(X,y,betamat,sig2mat,metaout_df,file = "./output using Feb22 data/MRM2_CAN_32obs.rda")
save(X,y,betamat,sig2mat,metaout_df,file = "./output using June22 data/MRM2_CAN_32obs.rda")

proc.time()-tic #capture run time

