options("rgdal_show_exportToProj4_warnings"="none")
library(inlabru) 
library(INLA)
library(rgdal)
library(RColorBrewer)
library(ggplot2)
library(sp)
library(raster)
library(ggpubr)
library(ggpolypath)
library(dplyr)
library(readxl)
library(rgeos)

sessionInfo()
source("Source_covariates.R") 


##colour scale
colsc <- function(...) {
  scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(11,"RdYlBu")),
                       limits = range(..., na.rm=TRUE))
}


##projection
crs_km <- sp::CRS("+proj=utm +zone=18 +south +ellps=GRS80 +units=km +no_defs") 

aacf <- readRDS("new_aac_spatialpixelsdfv2.rds")
roostsf <-read_excel("senasa_fieldt_km_MODELS.xlsx")
coords <- cbind(roostsf$x, roostsf$y)
colnames(coords) <- c("x", "y")
vbr <- SpatialPointsDataFrame(coords=coords, data=roostsf, proj4string = crs_km)
expval <- readRDS("boundaries.rds")
crs(expval) <- crs(vbr)


####MESH 
bound1 <- inla.sp2segment(expval)
mesh1 <- inla.mesh.2d(loc = coords, boundary = bound1, cutoff=2, max.edge=c(4,40), min.angle=c(25, 21), crs=crs_km)

plot(mesh1)  
plot(vbr, col="red", cex = 1, pch=20, add=TRUE) 


pcmatern1 <- inla.spde2.pcmatern(mesh1, alpha=2,      
                                 prior.sigma = c(0.2, 0.001),  
                                 prior.range = c(50, 0.01)) 


df <- pixels(mesh1, mask=expval)

effort <- readRDS("grid_eff_region.rds")
plot(effort)


######################Variables to test
# The ones I used in the model, but one by one - repeat the same procedure with the new resolution ones - slope, TRI and ebin 
# The new covariates - min temp, temp seasonality and precipitation seasonality 


##new covs - slope4, elev4, TRI4, tseas4, pseas4, tmin4, ebin4

comp_sp1 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + slope(slope4, model="linear") 

fit_sp1 <- lgcp(comp_sp1,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + slope, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))

INLA:::summary.inla(fit_sp1) 
sum(fit_sp1$cpo$cpo, na.rm=TRUE)





comp_sp2 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + ebinar(ebin4, model="factor_contrast") 

fit_sp2 <- lgcp(comp_sp2,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + ebinar, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))


INLA:::summary.inla(fit_sp2) 
summary(fit_sp2$summary.random$ebinar)
sum(fit_sp2$cpo$cpo, na.rm=TRUE)

comp_sp3 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + Tseas(tseas4, model="linear") 

fit_sp3 <- lgcp(comp_sp3,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + Tseas, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))


INLA:::summary.inla(fit_sp31) 
sum(fit_sp31$cpo$cpo, na.rm=TRUE)

############################try with 1km resolution 
#tseask1 <- as(tseas1k1, "SpatialPixelsDataFrame") 
comp_sp3 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + Tseas1k(tseask1, model="linear") 

fit_sp3 <- lgcp(comp_sp3,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + Tseas1k, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))


INLA:::summary.inla(fit_sp3) 
sum(fit_sp3$cpo$cpo, na.rm=TRUE)






comp_sp4 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + Pseas(pseas4, model="linear") 

fit_sp4 <- lgcp(comp_sp4,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + Pseas, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))


INLA:::summary.inla(fit_sp41) 
sum(fit_sp41$cpo$cpo)

#################1km 
#pseask1 <- as(pseas1k1, "SpatialPixelsDataFrame") 

comp_sp4 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + Pseas1k(pseask1, model="linear") 

fit_sp4 <- lgcp(comp_sp4,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + Pseas1k, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))


INLA:::summary.inla(fit_sp4) 
sum(fit_sp4$cpo$cpo)





comp_sp5 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + tminc(tmin4, model="linear") 

fit_sp5 <- lgcp(comp_sp5,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + tminc, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))


INLA:::summary.inla(fit_sp51) 
sum(fit_sp51$cpo$cpo)

#################1km 
#tmink1 <- as(tmin1k1, "SpatialPixelsDataFrame")
comp_sp5 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + Tminc1k(tmink1, model="linear") 

fit_sp5 <- lgcp(comp_sp5,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + Tminc1k, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))


INLA:::summary.inla(fit_sp5) 
sum(fit_sp5$cpo$cpo)


#######################For meanT
comp_sp61 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + 
  meantemp(mtemp, model="linear")

fit_sp61 <- lgcp(comp_sp61,
                 data = vbr,
                 samplers = expval,
                 domain = list(coordinates = mesh1),
                 formula = coordinates ~ Intercept + mysmooth + log(eff) + meantemp, 
                 options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                                control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                                verbose = TRUE))


INLA:::summary.inla(fit_sp61) 
sum(fit_sp61$cpo$cpo)


#tmeank1 <- as(mtemp1k1, "SpatialPixelsDataFrame")

comp_sp6 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + Meantemp1k(tmeank1, model="linear")

fit_sp6 <- lgcp(comp_sp6,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + Meantemp1k, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))


INLA:::summary.inla(fit_sp6) 
sum(fit_sp6$cpo$cpo)

################




comp_sp6 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + 
  meantemp(mtemp, model="linear") + mtemp2(mtemp, model="linear")

fit_sp6 <- lgcp(comp_sp6,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + meantemp + mtemp2, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))



INLA:::summary.inla(fit_sp6) 
sum(fit_sp6$cpo$cpo)



comp_sp7 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") +
  precipitation(prec, model="linear")


fit_sp7 <- lgcp(comp_sp7,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + precipitation, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))



INLA:::summary.inla(fit_sp7) 
sum(fit_sp7$cpo$cpo)


comp_sp8 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") +
  TRIx(TRI4, model="linear")


fit_sp8 <- lgcp(comp_sp8,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + TRIx, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))



INLA:::summary.inla(fit_sp8) 
sum(fit_sp8$cpo$cpo)


comp_sp9 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") +
  treecover(treecov, model="linear")


fit_sp9 <- lgcp(comp_sp9,
                data = vbr,
                samplers = expval,
                domain = list(coordinates = mesh1),
                formula = coordinates ~ Intercept + mysmooth + log(eff) + treecover, 
                options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                               control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                               verbose = TRUE))



INLA:::summary.inla(fit_sp9) 
sum(fit_sp9$cpo$cpo)


comp_sp10 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") +
  cropcover(crop, model="linear")


fit_sp10 <- lgcp(comp_sp10,
                 data = vbr,
                 samplers = expval,
                 domain = list(coordinates = mesh1),
                 formula = coordinates ~ Intercept + mysmooth + log(eff) + cropcover, 
                 options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                                control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                                verbose = TRUE))



INLA:::summary.inla(fit_sp10) 
sum(fit_sp10$cpo$cpo)


comp_sp11 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + 
  centrospoblados(centros2k, model="linear")


fit_sp11 <- lgcp(comp_sp11,
                 data = vbr,
                 samplers = expval,
                 domain = list(coordinates = mesh1),
                 formula = coordinates ~ Intercept + mysmooth + log(eff) + centrospoblados, 
                 options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                                control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                                verbose = TRUE))



INLA:::summary.inla(fit_sp11) 
sum(fit_sp11$cpo$cpo)

##centros poblados 5km 
comp_sp11 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + 
  centrospoblados5(centros5k, model="linear")


fit_sp11 <- lgcp(comp_sp11,
                 data = vbr,
                 samplers = expval,
                 domain = list(coordinates = mesh1),
                 formula = coordinates ~ Intercept + mysmooth + log(eff) + centrospoblados5, 
                 options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                                control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                                verbose = TRUE))



INLA:::summary.inla(fit_sp11) 
sum(fit_sp11$cpo$cpo)


####Cattle density at 5km 

comp_spc5 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + 
  cattle(cattle5k, model="linear")


fit_spc5 <- lgcp(comp_spc5,
                 data = vbr,
                 samplers = expval,
                 domain = list(coordinates = mesh1),
                 formula = coordinates ~ Intercept + mysmooth + log(eff) + cattle, 
                 options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                                control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                                verbose = TRUE))



INLA:::summary.inla(fit_spc5) 
sum(fit_spc5$cpo$cpo)


#####Cattle density at 10km 
comp_spc10 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + 
  cattle10(cattle10k, model="linear")


fit_spc10 <- lgcp(comp_spc10,
                  data = vbr,
                  samplers = expval,
                  domain = list(coordinates = mesh1),
                  formula = coordinates ~ Intercept + mysmooth + log(eff) + cattle10, 
                  options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                                 control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                                 verbose = TRUE))



INLA:::summary.inla(fit_spc10) 
sum(fit_spc10$cpo$cpo)





###################################Model selection by WAIC - add previous variables 1 by 1, except variables high high VIF
comp_full1 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + centrospoblados(centros2k, model="linear") + 
  ebinar(ebin4, model="factor_contrast") + tminc(tmin4, model="linear") + precipitation(prec, model="linear") + cropcover(crop, model="linear")

fit_full1 <- lgcp(comp_full1,
              data = vbr,
              samplers = expval,
              domain = list(coordinates = mesh1),
              formula = coordinates ~ Intercept + mysmooth + log(eff) + centrospoblados + ebinar + tminc + precipitation + cropcover, 
              options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                             control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                             verbose = TRUE))


INLA:::summary.inla(fit_full1)
sum(fit_full1$cpo$cpo)
summary(fit_full1$summary.random$ebinar)



comp_full2 <- coordinates ~ Intercept(1) + mysmooth(coordinates, model=pcmatern1) + eff(effort, model="offset") + 
  ebinar(ebin4, model="factor_contrast") + meantemp(mtemp, model="linear") + mtemp2(mtemp, model="linear") + 
  centrospoblados(centros2k, model="linear") + precipitation(prec, model="linear") + cropcover(crop, model="linear")


fit_full2 <- lgcp(comp_full2,
                 data = vbr,
                 samplers = expval,
                 domain = list(coordinates = mesh1),
                 formula = coordinates ~ Intercept + mysmooth + log(eff) + meantemp + mtemp2 + centrospoblados + precipitation + cropcover + ebinar, 
                 options = list(control.compute = list(dic = TRUE, waic = TRUE, config = TRUE, cpo = TRUE),
                                control.inla = list(strategy = 'simplified.laplace', int.strategy='eb'),
                                verbose = TRUE))


INLA:::summary.inla(fit_full2) 
summary(fit_full2$summary.random$ebinar)
sum(fit_full2$cpo$cpo)