############ Length prelim analyses ############
library(ggplot2)
library(MuMIn)
library(optimx)
library(nlme)
library(AICcmodavg)
library(mapproj)

### read in data ##############
diveDat = read.csv("Puget Sound_RCA_DataForAnalysis.csv")
# data already skew corrected

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

#######  AIC tables ##########
diveDatSub = diveDat

# changing coordinates from Lat Long
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)

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")

sink("20200407_Puget Sound TL model summary_MLselection.txt")
print(summary(Allk2))
sink()

sink("20200407_Puget Sound TL Dive Length AIC table_MLselection.txt")
print(Allk)
sink()

### prediction tables
newDat = data.frame(Depth = rep(rep(c(5,10,15,20,30),each=32),7),
                    Complexity = rep(rep(c(1,2,3,4),each=4),70))

test = predictSE(Allk2,newDat,se.fit=T)

newDat$predictions = test[[1]]
newDat$SE = test[[2]]

write.csv(newDat,"PugetSoundPredictionTable_MLselection_DepthComplex.csv")