############ Length prelim analyses ############
library(ggplot2)
library(glmmTMB)
library(lme4)
library(lmerTest)
library(DHARMa)
library(MuMIn)
library(optimx)
library(mefa)
library(nlme)
library(ape)
library(AICcmodavg)
library(mapproj)

### read in data ##############
#set an appropriate working directory
#read data
diveDat = read.csv("Puget Sound_RCA_DataForAnalysis.csv")

# data already skew corrected

diveDat$Date <- as.POSIXct(diveDat$Date, format = "%d-%m-%y")
diveDat$Year <- year(diveDat$Date)

# average complexity
diveDat$Complexity = apply(diveDat[,c("Complex1","Complex2","Complex3")],1,function(x)mean(x,na.rm=T))


#######  AIC tables ##########
diveDatSub = diveDat
#filter by sample size
SampleSizeAllSpModel=data.frame(tapply(diveDatSub$Multiplier,list(diveDatSub$InRCA,diveDatSub$RCA_Name),function(x)sum(x)))

#output a sample size table like this
SampleSizeAllSpModel2=data.frame(tapply(diveDatSub$Multiplier,list(diveDatSub$Species,diveDatSub$RCA_Name,diveDatSub$InRCA),function(x)sum(x)))

##### model aic selection ##########################

set.seed(1)
AlbersCoords = mapproject(diveDatSub$Lat, diveDatSub$Long, "albers", param=c(51, 53))
diveDatSub$X = AlbersCoords$x
diveDatSub$Y = AlbersCoords$y

diveDatSub$X = diveDatSub$X+runif(length(diveDatSub$X),min = -0.00001,max = 0.00001)
diveDatSub$Y = diveDatSub$Y+runif(length(diveDatSub$Y),min = -0.00001,max = 0.00001)


#add a tiny bit of spatial jitter
#set.seed(1)
#diveDatSub$Lat = diveDatSub$Lat+runif(length(diveDatSub$Lat),min = -0.000000001,max = 0.000000001)
#diveDatSub$Long = diveDatSub$Long+runif(length(diveDatSub$Lat),min = -0.000000001,max = 0.000000001)


print("All Species: Global Model")

start_time <- Sys.time()
mk = gls(data= diveDatSub,
         TLnormalized~InRCA*RCA_Name+Complexity+Depth,
         correlation=corGaus(form=~X+Y,nugget=TRUE),
         weights=varFixed(~Multiplier),
         method="ML")
end_time <- Sys.time()
end_time - start_time



print("dredge")

options(warn=1)

Allk=dredge(mk,fixed = c("Complexity","Depth"),trace=T,rank="AICc")

### best models

Allk2=get.models(Allk,subset=delta==0)[[1]]
Allk2=update(Allk2,method="REML")

# plot model parameter values ###################

library(MuMIn)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggplot2)
library(snakecase)

RVI = data.frame(variable = c("Complexity","Depth","InRCA","RCA_Name","InRCA-RCA_Name"),
                 SumWts = c(1,1,0.4667,0.2751,0.111))

# sink("puget sound RVI.txt")
# print(RVI)
# sink()


jpeg("20200630_variable plot PugetSound_highRes.jpeg", width = 1360, height = 1360,
     pointsize = 12, res=300)

set_theme(base = theme_classic(),
          axis.title.size = 0.9,
          axis.textsize.x = 0.9,
          axis.textsize.y = 0.9,
          title.align = "r")

plot_model(Allk2,  axis.title = "Standardized Coefficient", 
           title = "Predictor                                                                  RVI", 
           color = c("black", "black"),
           wrap.title = 150,wrap.labels = 50)


dev.off()