##
#
#             HMSC year model - simple script
#
##--------------------------------------------------------------------------

# load prepared data
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/2.ObsInfo")
load('prep_centCLC_J23.RDATA')


## ++++++++++++++++++++++++++++++++++++++++++++++ Modeling with HMSC
library(Hmsc)

#####        with random effects
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## fitting param
nChains=4
thin=100
samples=250
transient=ceiling(0.5*samples*thin)
verbose=ceiling(0.5*samples*thin)


# Env info
XData = data.frame(prec=Env.br$precip, tempWo= Env.br$GDD5w, tempW= Env.br$CDD5w, tempS= Env.br$GDD5s, snow=Env.br$snow, clumpy=Env.br$clumpy,
                   ent = Env.br$ent, enn=Env.br$enn, BLforest=Env.br$BLforest, Cforest=Env.br$Cforest, Mforest=Env.br$Mforest,
                   SNHherb=Env.br$SNHherb, Effort=log(Env.br$Effort))
dim(XData)

# Community info
Ystart = sp40Y_spread[, -c(1,2)]
dim(Ystart)

Ystart[is.na(Ystart)] = 0

# Trait
TrData = data.frame(Species = spT$Species, WingspanF = spT$WingspanF, HostPlantForm = spT$HostPlantForm, Diet = spT$Diet)
# re order
TrData = TrData[ order(match(TrData$Species, colnames(Ystart))),]
TrData = TrData[, -1]
rownames(TrData) = colnames(Ystart)

###  Define spatio temporal design
sample.id = Env.br$SiteID
year = Env.br$YearStart
region = Env.br$BR


# Define study design
studyDesign=data.frame(sample=as.factor(sample.id), year=as.factor(year), region=as.factor(region))

# define random effect
rLsample=HmscRandomLevel(units = sample.id)
rLregion=HmscRandomLevel(units = region)
rLyear=HmscRandomLevel(units = year)


# Formula
Fallpoly2 =~ poly(prec, degree=2, raw=TRUE) +  poly(tempW, degree=2, raw=TRUE)  + 
  poly(tempS, degree=2, raw=TRUE) + poly(snow, degree=2, raw=TRUE) + 
  clumpy + ent +  #enn + 
  poly(BLforest, degree=2, raw=TRUE) + poly(Cforest, degree=2, raw=TRUE) + 
  poly(Mforest, degree=2, raw=TRUE) + poly(SNHherb, degree=2, raw=TRUE) + Effort

formTr1 = ~ HostPlantForm + WingspanF


##                                                        poisson
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

models=list()

for (i in 1) {
  
  if(i==1){ Y=1*(Ystart>0)}
  if(i==2){ 
    Y = Ystart
    Y[Y==0] = NA
    Y =log(Y)
  }
  
  tmp = list()
  #for (j in 1) {
  m = Hmsc(Y=Y, XData=XData, XFormula= Fallpoly2, distr = switch(i, 'probit', 'normal'),
           TrData = TrData, TrFormula = formTr1,
           studyDesign=studyDesign, ranLevels = list('sample'=rLsample, 'year'=rLyear,
                                                     'region'=rLregion),
           YScale=TRUE)
  
  # tmp[[j]]=m
  #}
  models[[i]]=m
}

#start_time <- Sys.time()
for (i in 1) {
  models[[i]] = sampleMcmc(models[[i]], thin=thin, samples=samples, transient=transient, nChains=nChains,
                           nParallel = nChains,
                           verbose=verbose, adaptNf=rep(ceiling(0.4*samples*thin),models[[i]]$nr))
}

#models[[1]]
W1 = computeWAIC(models[[1]])
#W2 = computeWAIC(models[[2]])
#WAIC= c(W1, W2)


save(models, W1, file='test3107cCLC_yearPA.RDATA')
