
library(Hmsc)
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/3.ModelRun")


##...................................................................................................................
### Results from HMSC run:

load('test10908cent_yearPA.RDATA')
modelyear = list()
modelyear[[1]] = models[[1]]

load('test10908cent_yearAB.RDATA')
modelyear[[2]] = models[[2]]

##...................................................................................................................
### Diagnostics

# check convergence
for (i in 1) {
  mpost = convertToCodaObject(modelyear[[i]])
  png(filename = )
  plot(mpost$Beta) 
}
dev.off()


ess.beta = gd.beta = ess.omega = gd.omega =  ess.Gamma = gd.Gamma = rep(list(list()), 2)
for (i in 1:2) {
  mpost = convertToCodaObject(modelyear[[i]])
  
  ess.beta[[i]] <- effectiveSize(mpost$Beta)
  gd.beta[[i]] <- gelman.diag(mpost$Beta, multivariate = FALSE)$psrf
  
  ess.omega[[i]] <- effectiveSize(mpost$Omega[[1]])
  gd.omega[[i]] <- gelman.diag(mpost$Omega[[1]], multivariate = FALSE)$psrf
  
  ess.Gamma[[i]] <- effectiveSize(mpost$Gamma)
  gd.Gamma[[i]] <- gelman.diag(mpost$Gamma, multivariate = FALSE)$psrf
  
}


#summary(gelman.diag(mpost$Beta))



par(mfrow=c(2,3), oma = c(0, 0, 0, 0),
    mar = c(4,4,3,0))
main = c("PA - hurdle P1", "AB - hurdle P2" )
#ylab = c('0.05km', '0.1km', '1.5km', '1.5km', '1.5km', '6km')
for (i in 1) {
  hist(ess.beta[[i]], main =main[i], xlab='Eff size beta')
}

for (i in 1) {
  hist(ess.Gamma[[i]], main =main[i], xlab='Eff size gamma')
}

for (i in 1) {
  hist(ess.omega[[i]], main =main[i], xlab='Eff size omega')
}


for (i in 1) {
  hist(gd.beta[[i]], main =main[i], xlab='potential scale reduction factor beta')
}

for (i in 1) {
  hist(gd.Gamma[[i]], main =main[i], xlab='potential scale reduction factor gamma')
}

for (i in 1) {
  hist(gd.omega[[i]], main =main[i], xlab='potential scale reduction factor omega')
}


summary(ess.beta[[1]])
summary(ess.omega[[1]])
summary(ess.Gamma[[1]])

summary(gd.beta[[1]])
summary(gd.omega[[1]])
summary(gd.Gamma[[1]])



summary(ess.beta[[2]])
summary(ess.omega[[2]])
summary(ess.Gamma[[2]])

summary(gd.beta[[2]])
summary(gd.omega[[2]])
summary(gd.Gamma[[2]])


# Model performance

# perf# perf
preds = MF = TabMF =  list()
for (i in 1) {
  preds[[i]] = computePredictedValues(modelyear[[i]])
  MF[[i]]= evaluateModelFit(hM=modelyear[[i]], predY = preds[[i]])
  
  TabMF[[i]] = cbind(as.data.frame(MF[[i]]), species = modelyear[[i]]$spNames)
}


#TabMF

mat.MF =  cbind(TabMF[[1]]$AUC, TabMF[[1]]$TjurR2, TabMF[[1]]$species)
colnames(mat.MF) = c('AUC', 'TjurR2', 'species')


# the first column
AUC <- as.data.frame(mat.MF[,1])
AUC=apply(AUC, 2, as.numeric)

TjurR2 <- as.data.frame(mat.MF[,2])
TjurR2=apply(TjurR2, 2, as.numeric)

par(mfrow=c(1,1), xpd=F)
plot(TjurR2~AUC, xlim=c(0,1), ylim=c(0,1))
abline(coef = c(0,1), col='red')

par(mfrow=c(2,1))

boxplot(AUC, ylab='AUC')
points(apply(AUC, 2, mean), col='red')

boxplot(TjurR2, ylab='Tjur R2')
points(apply(TjurR2, 2, mean), col='red')



library(ggplot2)
temp.auc <- cbind.data.frame(auc = as.numeric(TabMF[[1]]$AUC), 
                             species = as.character(modelyear[[1]]$spNames))
ggplot(data = temp.auc, aes(y=auc, x=species))+
  geom_point()+ theme(axis.text.x = element_text(angle = 90, vjust = 1.5, hjust=1))


temp.TjurR2 <- cbind.data.frame(TjurR2 = as.numeric(TabMF[[1]]$TjurR2), 
                                species = as.character(modelyear[[1]]$spNames))

ggplot() + geom_point(data=temp.TjurR2, aes(y=TjurR2, x=species)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 1.5, hjust=1))


##...................................................................................................................
## predictive power based on cross-validation  # 
partition = createPartition(modelyear[[1]], nfolds = 2, column = 'sample')


preds.part = computePredictedValues(modelyear[[1]], partition = partition)
MF.part = evaluateModelFit(hM = modelyear[[1]], predY = preds.part)

mat.MF.part =  cbind(MF.part$AUC, MF.part$TjurR2)
colnames(mat.MF.part) = c('AUC', 'TjurR2')

AUC.part <- as.data.frame(mat.MF.part[,1])
AUC.part=apply(AUC.part, 2, as.numeric)

TjurR2.part <- as.data.frame(mat.MF.part[,2])
TjurR2.part = apply(TjurR2.part, 2, as.numeric)

par(mfrow=c(2,1))
boxplot(cbind(AUC, AUC.part), ylab='AUC', col=c('gray28', 'gray65'), xaxt='n')
axis(1, at=c(1,2), labels=c('expl', 'pred'))
points(c(apply(AUC, 2, mean),apply(AUC.part, 2, mean)), col='red')

boxplot(cbind(TjurR2, TjurR2.part), ylab='Tjur R2', col=c('gray28', 'gray65'), xaxt='n')
axis(1, at=c(1,2), labels=c('expl', 'pred'))
points(c(apply(TjurR2, 2, mean),apply(TjurR2.part, 2, mean)), col='red')

