# Modeling rock glacier stability in relation to terrain characteristics
# Marcer et al, 2018. Inferring the degradation susceptibility of mountain permafrost in the French Alps using an inventory of destabilized rock glaciers 
#
# Code created the 13/03/2018
#
# WARNING : THIS CODE WAS BUILT UNDER R VERSION 3.2.2
# PACKAGE sperrorest may not work with the present syntax in later versions
##########################################################################################################

# Instructions :

# 1 . Requires the packages "sperrorest", "rgdal", "mgcv". Use install.packages()
# 2 . Set working directory to PATH. Use "/" in the path string

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

# Step 0 - Set up working station

  PATH <- "C:/Users/.../Code"
	setwd(PATH)

	#Path to database
	VECPATH <- paste(PATH,"data",sep='/')

	#Load required libraries
	library(rgdal)
	library(sperrorest)
	library(mgcv)

# Model specifications 

	nsp<-5 	# How many points per rock glacier to be used in the model?

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

# Step 1 - Prepare Modeling Database

# 1.1 Load shapefile with sampling points and initialize database
	P<-readOGR(dsn=VECPATH,layer="sampling_points") # load sampling points in rock glaciers. Sampling already performed

	# Define which degree of destabilization use in the modeling 
	P1f<-subset(P,ID.1==0 & Dest_Index=="3") # unstable areas of rock glaciers (ID.1 = 0) with destabilization index = 3
	P0f<-subset(P,is.na(ID.1)==TRUE & strtoi(Dest_Index)<2 & Moving==strtoi(1)) # Stable areas of observed active rock glaciers with destabilization index = 0 and 1

	# Assign dependent variables
	P1f$p<-1
	P0f$p<-0

	# Initialize Database
	DB<-rbind(P1f,P0f)

# 1.2 Get IDs of rock glaciers (There are multiple points per rock glacier)

	idx<-unique(as.numeric(as.character(DB$ID_complet)))	  
	nP<-length(idx)

# 1.3 Obtain a random subsample of nsp points per rock glacier. All rock glaciers are used
	
	# Initialize modeling database for modeling based on model specifications
	NDB<-array(NA,dim=c(nsp*nP,ncol(DB)+2))
	colnames(NDB)<-names(as.data.frame(DB))

	# Create database for modeling based on model specifications
	for (i in 1:nP){
		temp<-subset(DB,ID_complet==idx[[i]])								# Get points in the ith rock glaciers
		if (length(temp)>nsp){mysample <- temp[sample(1:nrow(temp),nsp,replace=FALSE),]}	# If in the ith rock glaciers there are more than nsp points, than random subset
		else{mysample <- temp[sample(1:nrow(temp),nsp,replace=TRUE),]}				# Otherwise, random subset with replacement
		NDB[(nsp*i-nsp+1):(nsp*i),]<-as.numeric(as.matrix(as.data.frame(mysample,stringsAsFactors = FALSE)))
		}

# 1.4 Final arrangements to the ,modeling database

	destab<-data.frame(NDB)
	
	destab$x<-destab$X
	destab$y<-destab$Y
	destab$P<-factor(destab$p)

	colnames(destab)<-colnames(NDB) 

	keeps <- c("X","Y","p","DEM25","PISR25","TWI","PTP","Slope25","Downslope_C","ID_complet")
	destab<-destab[,keeps]
	
	destab$unstable <- factor(destab$p == 1)

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

# Step 2 - Sperrorest

# Sperrorest using GAM and partition.kmeans

  # Note that I am limiting the flexibility of the GAM
	# manually using s(...,k=4) to avoid weird oscillating
	# behaviour in the transformation function

	fo<- unstable ~ DEM25 + Slope25 + PISR25 + Downslope_C
 
	my_gam <- function(formula, data, ...) {
  		response <- all.vars(formula)[1]
  		predictors <- all.vars(formula)[-1]
  		s_pred <- paste("s(", predictors, ",k=4)", sep="")
  		formula <- paste( response, "~",
                    paste(s_pred, collapse="+") )
  		formula <- as.formula(formula)
  		fit <- mgcv::gam(formula, data, ...)
  		return(fit)
	}

	# check that it works:
	testfit <- my_gam(fo, destab, family="binomial")
	summary(testfit)
	plot(testfit)
	
	out <- sperrorest(data = destab, formula = fo,
	                  coords = c("X","Y"),
	                  model_fun = my_gam,
	                  model_args = list(family = "binomial"),
	                  pred_fun = predict,
	                  pred_args = list(type = "response"),
	                  smp_fun = partition_kmeans,
	                  smp_args = list(nfold=2, repetition=1:3, seed1=123), #use 100 repetitions in final run
	                  err_fun = err_default,
                   	 par_args = list(par_mode="sequential"),
                  	 importance = TRUE, 
	                  imp_permutations = 10 # use 100 in your final run
	                  # --> may take several hours
				)
	
	# Training set AUROC, and spatial cross-validation estimate of AUROC:	
	#summary(out$error_rep)
	round(summary(out$error_rep)[c("train_auroc","test_auroc"),c("mean","sd")],3)
	
	imp <- summary(out$importance)
	nms <- rownames(imp)
	imp <- imp[,"mean.auroc"]
	names(imp) <- nms
	# mean AUROC reduction based on 
	# permutation of each predictor:
	sort(imp, decr=TRUE) 
	
	out.gam <- out
	imp.gam <- imp
	
	save(out.gam, imp.gam, file="results_gam_dem.Rdata")
	
 #####################################################################################################################

# Step 3 - Transformation Plots

	point0<-subset(destab,p==0)
	point1<-subset(destab,p==1)

	# Some graphic parameters
	ep=0.2
	M=4	
	m=-7

	# GAM:
	fit <- my_gam(fo, destab, family="binomial")
	summary(fit)

	plot(fit,scheme = 1,unconditional = TRUE,pages=0,select=1,ylim=c(-8,5),cex.axis=2)
	points(point1$DEM25,runif(length(point1$DEM25),M-ep,M+ep),
	pch=19,col=rgb(0.5,0.5,0.5,alpha=0.3),cex=0.7)
	points(point0$DEM25,runif(length(point0$DEM25),m-ep,m+ep),
	pch=19,col=rgb(0.5,0.5,0.5,alpha=0.3),cex=0.7)

	plot(fit,scheme = 1,unconditional = TRUE,pages=0,select=2,ylim=c(-8,5),cex.axis=2)
	points(point1$Slope25,runif(length(point1$Slope25),M-ep,M+ep),
	pch=19,col=rgb(0.5,0.5,0.5,alpha=0.3),cex=0.7)
	points(point0$Slope25,runif(length(point0$Slope25),m-ep,m+ep),
	pch=19,col=rgb(0.5,0.5,0.5,alpha=0.3),cex=0.7)

	plot(fit,scheme = 1,unconditional = TRUE,pages=0,select=3,ylim=c(-8,5),cex.axis=2)
	points(point1$PISR25,runif(length(point1$PISR25),M-ep,M+ep),
	pch=19,col=rgb(0.5,0.5,0.5,alpha=0.3),cex=0.7)
	points(point0$PISR25,runif(length(point0$PISR25),m-ep,m+ep),
	pch=19,col=rgb(0.5,0.5,0.5,alpha=0.3),cex=0.7)

	plot(fit,scheme = 1,unconditional = TRUE,pages=0,select=4,ylim=c(-8,5),cex.axis=2)
	points(point1$Downslope_C,runif(length(point1$PISR25),M-ep,M+ep),
	pch=19,col=rgb(0.5,0.5,0.5,alpha=0.3),cex=0.7)
	points(point0$Downslope_C,runif(length(point0$PISR25),m-ep,m+ep),
	pch=19,col=rgb(0.5,0.5,0.5,alpha=0.3),cex=0.7)

