
###################### Range-wide indicators of African great ape density distribution
###################### Supporting Information 2

#Install and load packages, functions ## see line 1162 for coding of the self-written functions 
install.packages("polyclip")
install.packages("shapefiles")
install.packages("splancs")
install.packages("maps")
install.packages("gtools")
install.packages("sp")
library(polyclip)
library(shapefiles)
library(splancs)
library(maps)
library(gtools)
library(sp)
setwd("")##need to set the appropriate working directory
source("functions.r")
source("x:/Statistics/R scripts/Roger/no_int_plot.r")
source("X:/Statistics/R scripts/Roger/glmnb2.r")


###Bring in data
abundance.data.old=read.table(file="Y:/primint/Isabel/ape_density_distribution/R/Model_jan2020/abundance_data_utm_fixed_ifl_gdpconst.txt", header=T, sep="\t")
abundance.data.old=subset(abundance.data.old, !is.na(X) & !is.na(Y))
#New column for number of taxa
abundance.data.old$taxa.num=1
abundance.data.old$taxa.num[grepl(x=as.character(abundance.data.old$species), pattern=" and", fixed=T)]=2
abundance.data.old$taxa.num.d=abundance.data.old$taxa.num-1

###Add some new data:                 ## NOTE: all data, including the ones added here, are found in the Supplementary info. 3 dataset
new.ab.data=read.table(file="Y:/primint/Isabel/ape_density_distribution/R/Model_jan2020/update_abund_estimates_jan2020.txt", header=T, sep="\t") 
new.ab.data=new.ab.data[!is.na(new.ab.data$year), ]
xx=as.character(abundance.data.old$year_of_es)
xx[is.na(xx)]=2015
xx=lapply(strsplit(xx, split="-", fixed=T), as.numeric)
xx=lapply(xx, function(x){
	if(length(x)>1){
		return(x[1]:x[length(x)])
	}else{
		return(x)
	}
})
abundance.data.old$year_of_es=unlist(lapply(xx, mean))
abundance.data.old$abund_esti[match(new.ab.data$site, abundance.data.old$site_name)]=new.ab.data$new_abund
abundance.data.old$year_of_es[match(new.ab.data$site, abundance.data.old$site_name)]=new.ab.data$year
###Bring in new sites
abundance.data.new=read.csv(file="Y:/primint/Isabel/ape_density_distribution/R/Model_june2020/abundance_data_new_sites_june.csv", header=T, sep=",")
abundance.data.new=subset(abundance.data.new, !is.na(X) & !is.na(Y))
#65 sites
#New column for number of taxa
abundance.data.new$taxa.num=1
abundance.data.new$taxa.num[grepl(x=as.character(abundance.data.new$species), pattern=" and", fixed=T)]=2
abundance.data.new$taxa.num.d=abundance.data.new$taxa.num-1
#Preparing to assign random year to pseudo-absence sites
xx=as.character(abundance.data.new$year_of_es)
xx=lapply(strsplit(xx, split="-", fixed=T), as.numeric)
xx=lapply(xx, function(x){
	if(length(x)>1){
		return(x[1]:x[length(x)])
	}else{
		return(x)
	}
})
abundance.data.new$year_of_es=unlist(lapply(xx, mean))

###Merge the two datasets
###Create new data frame with the columns we need
abundance.data.new$GDP.2000=NA
abundance.data.new$GDP.2016=NULL
setdiff(names(abundance.data.new), names(abundance.data.old))
setdiff(names(abundance.data.old), names(abundance.data.new))
###Merge
abundance.data<- merge(abundance.data.old, abundance.data.new, all = TRUE) #376 sites

write.table(x=abundance.data, file="Y:/primint/Isabel/ape_density_distribution/R/Model_june2020/abundance_data_merged.txt", row.names=F, col.names=T, sep="\t")
abundance.data=read.csv(file="Y:/primint/Isabel/ape_density_distribution/R/Model_june2020/abundance_data_merged.csv", header=T, sep=",") ###merged dataset with corrected coordinates for 5 sites

###Subset to areas within great ape range:
sec.range=read.shp("Y:/primint/Isabel/ape_density_distribution/env_vars/ga_range/all_apes_Merge.shp")$shp
sec.range=split.shape(shp=sec.range)

io=lapply(1:length(sec.range), function(x){
	print(x)
	x=polyoffset(A=list(list(x=sec.range[[x]]$X, y=sec.range[[x]]$Y)), 1)
	inout(pts=cbind(x=abundance.data$X, y=abundance.data$Y), poly=cbind(x=x[[1]]$x, y=x[[1]]$y), bound=T)
})

io=matrix(unlist(io), ncol=length(io), byrow=F)
io=apply(io, 1, any)

###Checking that the subset is correct
plot.shape2<-function(shape, shape.names=NULL, cex=1, xlim=NULL, ylim=NULL, add=F, col=NA){
	if(!add){
		if(length(xlim)==0){
			xlim=lapply(shape, function(x){range(x$X)})
		}
		if(length(ylim)==0){	
			ylim=lapply(shape, function(x){range(x$Y)})
		}
		xlocs=apply(matrix(unlist(xlim), ncol=2, byrow=T), 1, mean)
		ylocs=apply(matrix(unlist(ylim), ncol=2, byrow=T), 1, mean)
		xlim=range(unlist(xlim))
		ylim=range(unlist(ylim))
		par(mar=c(3, 3, 0.2, 0.2), mgp=c(1.7, 0.4, 0), tcl=-0.25)
		plot(1,1, type="n", asp=1, xlim=xlim, ylim=ylim)
	}
	lapply(shape, polygon, col=col)
	if(length(shape.names)>0){
		text(x=xlocs, y=ylocs, labels=shape.names, cex=cex)
	}
}

plot.shape2(shape=sec.range, shape.names=NULL, cex=1, xlim=NULL, ylim=NULL, col=adjustcolor("blue", alpha.f=0.5), add=F)
points(abundance.data$X, abundance.data$Y, pch=19, col=c("green", "yellow")[as.numeric(io)+1])
sel.abundance.data=subset(abundance.data, io)

###Check buffer, added aug 2018:
plot(abundance.data[, c("X", "Y")], asp=1)
points(sel.abundance.data[, c("X", "Y")], asp=1, pch=19)
for(i in 1:length(sec.range)){
	polygon(sec.range[[i]], col=NA, border="red")
}
#We took sites within 100km around great ape range
points(sel.abundance.data[, c("X", "Y")], asp=1, pch=19, col="red")

###Area size
sel.abundance.data$area.size.km=sel.abundance.data$area.size/(1000^2)

###Survey year:
xx=as.character(sel.abundance.data$year_of_es)
xx[is.na(xx)]=2015
xx=lapply(strsplit(xx, split="-", fixed=T), as.numeric)
xx=lapply(xx, function(x){
	if(length(x)>1){
		return(x[1]:x[length(x)])
	}else{
		return(x)
	}
})
sel.abundance.data$survey.year=unlist(lapply(xx, mean))
sel.abundance.data$survey.year=sel.abundance.data$year_of_es

###Deal with the pseudo-absence sites which all 'took place' in 2015:
fict.areas=as.character(sel.abundance.data$site_name[sel.abundance.data$survey.year==2015 & sel.abundance.data$abund_esti==0])
#158 pseudo-absence areas
###Prepare for matching of the survey year and predictor variable year
#predictors not changing:
constant.vars=c("precip.driest.month", "elevation", "prop.muslims")
#predictors changing by year:
changing.vars=c("GDP.const", "corruption", "human.footprint", "ifl.prop") 
changing.vars.suffixes=list(
	 GDP.const=2000:2016,
	 corruption=1996:2015,
	 human.footprint=c(1993, 2009),
	# poverty=c(1993, 2000),
	 #forest.in.neighborhood=c(1992, 2000),
	 #pop.dens=c(1995, 2000, 2005, 2010, 2015),
	 ifl.prop=c(2000, 2013, 2016)
	# ifl.conn=c(2000.2013, 2016)
	# veg.cov=c(2000, 2005, 2010)
	# prop.forest.loss=2001:2016
)
###Checking distributions	
hist(sel.abundance.data$precip.driest.month)
hist(log(sel.abundance.data$precip.driest.month+1)) #best
hist(sel.abundance.data$elevation)
hist(sqrt(sel.abundance.data$elevation)) #best
hist(sel.abundance.data$prop.muslims)
hist(log(sel.abundance.data$prop.muslims))#best
hist(sel.abundance.data$GDP.const)
hist(log(sel.abundance.data$GDP.const)) #best
hist(sel.abundance.data$human.foot.prt)
hist(log(sel.abundance.data$human.foot.prt+1))#best
hist(sel.abundance.data$corruption)
hist(log(sel.abundance.data$corruption)) #best
hist(sel.abundance.data$ifl.prop.2013) #best
hist(log(sel.abundance.data$ifl.prop.2013+1))
#Transformations
sel.abundance.data$log.precip=log(sel.abundance.data$precip)
sel.abundance.data$log.precip.driest.month=log(sel.abundance.data$precip.driest.month+1)
sel.abundance.data$sqrt.elevation=sqrt(sel.abundance.data$elevation)
#sel.abundance.data$sqrt.human.infl=sqrt(sel.abundance.data$human.infl)
sel.abundance.data$log.prop.muslims=log(sel.abundance.data$prop.muslims)

####Prepare predictor variable data for prediction of density distribution
prediction.data=read.table(file="y:/primint/Isabel/ape_density_distribution/R/Model_jan2020/prediction_data_ifl_rowcheck_RM_gdpconst.txt",header=T)
prediction.data$prop.ifl.cl.yr=prediction.data$prop.ifl.overl
	test.vars<-c("log.human.foot.prt", "sqrt.elevation", "log.GDP.const", "prop.ifl.cl.yr","log.corruption","log.precip.driest.month","log.prop.muslims","survey.year.r") 
	all.vars<-c("log.human.foot.prt", "sqrt.elevation", "log.GDP.const","prop.ifl.cl.yr","log.corruption","log.precip.driest.month","log.prop.muslims", "long", "lat", "country", "survey.year.r", "area.size") 
	prediction.data$survey.year.r=2015
	#names(prediction.data)[names(prediction.data)=="forest.in.neighbourhood"]="forest.in.neighborhood"
	pred.data.cell.width=0.08334
	u.lat=unique(prediction.data$lat)
	xx=unlist(lapply(u.lat, function(x){
		dist.wgs(N1=x, E1=0, N2=x, E2=pred.data.cell.width)
	}))
	xx=xx*pred.data.cell.width*1852*60/(1000^2)
	prediction.data$area.size=xx[match(prediction.data$lat, u.lat)]
	mean(prediction.data$area.size)
	range(prediction.data$area.size)
	plot(prediction.data$long, prediction.data$lat)
###Checking distributions
hist(prediction.data$precip.driest.month)
hist(log(prediction.data$precip.driest.month+1)) #best
hist(prediction.data$elevation)
hist(sqrt(prediction.data$elevation)) #best
hist(prediction.data$prop.muslims)
hist(log(prediction.data$prop.muslims+1))#best (range from 0)
hist(prediction.data$GDP.const)
hist(log(prediction.data$GDP.const)) #best
hist(prediction.data$human.foot.prt)
hist(log(prediction.data$human.foot.prt+1))#best
hist(prediction.data$corruption)
hist(log(prediction.data$corruption)) #best
hist(prediction.data$prop.ifl.cl.yr) #best
hist(log(prediction.data$prop.ifl.cl.yr+1))	
#Transformations
	prediction.data$sqrt.elevation=sqrt(prediction.data$elevation)
	prediction.data$log.prop.muslims=log(prediction.data$prop.muslims+1)
	prediction.data$log.precip.driest.month=log(prediction.data$precip.driest.month+1)
	prediction.data$log.human.foot.prt=log(prediction.data$human.foot.prt+1)
	prediction.data$log.GDP.const=log(prediction.data$GDP.const)
	prediction.data$log.GDP=log(prediction.data$GDP)
	prediction.data$log.corruption=log(prediction.data$corruption)
	prediction.data=subset(prediction.data, taxa.num>0) ##subsetting to most recent range
##end prepare for prediction

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

###First part: Model 
library(rela)
aic.fun<-function(xsd){
	xac.term=get.1d.ac.from.dist.mat.2(resis=resis, ac.sd=xsd, dist.mat=dist.mat, contr.fac=as.factor(model.data$site_name%in%fict.areas))
	model.data$ac.term=xac.term
	ifull=try(fit.fun(formula="abund_esti~1+z.log.human.foot.prt+z.sqrt.elevation+z.log.GDP.const+I(z.log.GDP.const^2)+z.prop.ifl.cl.yr+z.log.corruption+z.log.precip.driest.month+z.log.prop.muslims+z.survey.year.r+ac.term", 
		data=model.data, offset=model.data$area.size.km, tol=0.00001), silent=T)
	if(class(ifull)!="try-error"){
		return(-ifull$logLik)
	}else{
		return(-full$logLik)
	}
}

n.iter=100
all.res=vector("list", n.iter)
all.years=vector("list", n.iter)
for(i in 1:n.iter){
	print(i)
	xdone=F
	set.seed(i)
	while(!xdone){
		#deal with the pseudo-absence sites which all 'took place' in 2015, assigning random years to them:
		sel.abundance.data$survey.year.r=sel.abundance.data$survey.year
		sel.abundance.data$survey.year.r[sel.abundance.data$site_name%in%fict.areas]=sample(2005:2015, size=sum(sel.abundance.data$site_name%in%fict.areas), replace=T)
		all.years[[i]]=sel.abundance.data$survey.year.r[sel.abundance.data$site_name%in%fict.areas]

		for(j in 1:(length(changing.vars))){
			xx=unlist(lapply(1:nrow(sel.abundance.data), function(x){
				xx=changing.vars.suffixes[[j]][which.min(abs(sel.abundance.data$survey.year.r[x]-changing.vars.suffixes[[j]]))]
				sel.abundance.data[x, paste(c(changing.vars[j], xx), collapse=".")] ###links data from closest year
			}))
				sel.abundance.data[, changing.vars[j]]=xx
		}
		 sel.abundance.data$log.human.foot.prt=log(sel.abundance.data$human.footprint+1)
		 sel.abundance.data$log.corruption=log(sel.abundance.data$corruption)
		 sel.abundance.data$log.GDP.const=log(sel.abundance.data$GDP.const)
		sel.abundance.data$prop.ifl.cl.yr=sel.abundance.data$ifl.prop
		##############################################################################
		##Subset to keep complete cases
		test.vars<-c( "log.human.foot.prt", "sqrt.elevation", "log.GDP.const", "prop.ifl.cl.yr", "log.corruption", "log.precip.driest.month", "X", "Y", "area.size", "utm.zone", "mid.long", "mid.lat", "utm.long", "utm.lat", "country", 
					"country_co", "site_name", "abund_esti", "log.prop.muslims","survey.year.r") 
		model.data=subset(sel.abundance.data, complete.cases(sel.abundance.data[, test.vars]))  #to keep complete cases for test.vars
		nrow(model.data) 
		#z-transformations
		test.vars=c("log.human.foot.prt", "sqrt.elevation", "log.GDP.const", "prop.ifl.cl.yr", "log.corruption","log.precip.driest.month","log.prop.muslims","survey.year.r") 
		xx=data.frame(scale(model.data[, test.vars]))
		names(xx)=paste("z", names(xx), sep=".")
		model.data=data.frame(model.data, xx)
		#prepare for autocorrelation:
		dist.mat=lapply(1:nrow(model.data), function(x){
			dist.wgs.vectorized(N1=model.data$Y[x], E1=model.data$X[x], N2v=model.data$Y, E2v=model.data$X)
		})
		dist.mat=matrix(unlist(dist.mat), ncol=length(dist.mat), byrow=T)
		dist.mat=dist.mat/1000

		#Full null model comp:
		full=try(fit.fun(formula="abund_esti~1+z.log.human.foot.prt+z.sqrt.elevation+z.log.GDP.const+I(z.log.GDP.const^2)+z.prop.ifl.cl.yr+z.log.corruption+z.log.precip.driest.month+z.log.prop.muslims+z.survey.year.r"
			data=model.data, offset=model.data$area.size.km, tol=0.00001), silent=T)
		if(class(full)!="try-error"){
			#get autocorrelation in there:
			resis=full$pearson.resids
			xsd=optimize(aic.fun, lower=1, upper=20)
			xac.term=get.1d.ac.from.dist.mat.2(resis=resis, ac.sd=xsd$minimum, dist.mat=dist.mat, contr.fac=as.factor(model.data$site_name%in%fict.areas))
			model.data$ac.term=xac.term
			model.data$z.ac.term=as.vector(scale(xac.term))
			full=try(fit.fun(formula="abund_esti~1+z.log.human.foot.prt+z.sqrt.elevation+z.log.GDP.const+I(z.log.GDP.const^2)+z.prop.ifl.cl.yr+z.log.corruption+z.log.precip.driest.month+z.log.prop.muslims+z.survey.year.r+z.ac.term"
				data=model.data, offset=model.data$area.size.km, tol=0.00001), silent=T)
			null=try(fit.fun(formula="abund_esti~1+z.sqrt.elevation+z.log.precip.driest.month+z.survey.year.r+z.ac.term",
				data=model.data, offset=model.data$area.size.km, tol=0.00001), silent=T)
			if(class(full)[[1]]!="try-error" & class(null)[[1]]!="try-error"){
				chisq=2*(full$logLik-null$logLik)
				df=nrow(full$coefs)-nrow(null$coefs)
				P=pchisq(chisq, df=df, lower.tail=F)
				#store results:
				all.res[[i]]=list(
					full.coeffs=full$coefs, 
					fn.comp=data.frame(chisq=chisq, df=df, P=P),
					disp.par=full$disp.par,
					n=nrow(model.data),
					model.data=model.data,
					ac.term.opt=xsd,
					resis=resis
				)
				xdone=T
			}
		}
	}
}

###Multi model inference:
library(gtools)
allmodels=all.models(model="z.log.human.foot.prt+z.sqrt.elevation+z.log.GDP.const+I(z.log.GDP.const^2)+z.prop.ifl.cl.yr+z.log.corruption+z.log.precip.driest.month+z.log.prop.muslims+z.survey.year.r")
allmodels=paste(allmodels, "z.ac.term", sep="+")
allmodels=paste("abund_esti", allmodels, sep="~")
library(parallel)
cl <- makeCluster(getOption("cl.cores", detectCores()))
cl=cl[1:(length(cl)-2)]

for(i in 1:length(all.res)){
	print(i)
	model.data=all.res[[i]]$model.data
	clusterExport(cl=cl, varlist=c("model.data", "fit.fun"))
	all.m.res=parLapply(cl=cl, X=allmodels, fun=function(x){
		full=try(fit.fun(formula=x, data=model.data, offset=model.data$area.size.km, tol=0.00001), silent=T)
		if(class(full)[[1]]!="try-error"){
			ests=full$coefs[, "Estimate"]
			names(ests)=rownames(full$coefs)
			logLik=full$logLik
		}else{
			ests=NA
			logLik=NA
		}
		return(list(coef=ests, logLik=logLik))
	})
	save(file=paste(c("Y:/primint/Isabel/ape_density_distribution/R/Model_june2020/MMI_ifl_humaninf_ele_gdpconst_corr_minprec_foodtab_june20/test_mmi_", i, ".RData"), collapse="."), #saving the mmi for the 100 data sets
		list=c("all.m.res", "model.data"))
}
parLapply(cl=cl, X=1:length(cl), fun=function(x){rm(list=ls())})
stopCluster(cl)
length(all.m.res)

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

###Second part: predicted density distribution
print("predictions")
flist=paste("Y:/primint/Isabel/ape_density_distribution/R/Model_june2020/MMI_ifl_humaninf_ele_gdpconst_corr_minprec_foodtab_june20/test_mmi_.", 1:length(all.res), "..RData", sep="")
all.pred=vector("list", length(flist))
test.mean.ab=vector("list", length(flist))
thresh.density=10
kept.and.discarded=vector("list", length(flist))
for(i in 1:length(all.res)){
	print(i)
	load(flist[i])
	#mmi:
	sample.size=nrow(all.res[[i]]$model.data)
	all.ll=unlist(lapply(all.m.res, function(x){x$logLik}))
	all.k=unlist(lapply(all.m.res, function(x){length(x$coef)}))
	all.aicc=-2*all.ll+2*all.k+aic.c.fac(N=sample.size, k=all.k)
	models2keep=!is.na(all.aicc)
	models.conv=allmodels[!is.na(all.aicc)]
	all.aicc=all.aicc[!is.na(all.aicc)]
	c.set=conf.set(aic=all.aicc)
	model.data=all.res[[i]]$model.data
	pred.env.fa=prediction.data
	pred.env.fa=subset(pred.env.fa, complete.cases(pred.env.fa[, c(test.vars)])) 
	z.pred.data=data.frame(scale(pred.env.fa[, test.vars], center=apply(model.data[, test.vars], 2, mean), 
		scale=apply(model.data[, test.vars], 2, sd)))
	colnames(z.pred.data)=paste("z", colnames(z.pred.data), sep=".")
	z.pred.data=model.matrix(object=~z.log.human.foot.prt+z.sqrt.elevation+z.log.GDP.const+I(z.log.GDP.const^2)+z.prop.ifl.cl.yr+z.log.corruption+z.log.precip.driest.month+z.log.prop.muslims+z.survey.year.r, data=as.data.frame(z.pred.data))
	all.preds=lapply(all.m.res, function(x){
		ests=x$coef
		if(length(ests)>1){
			ests=ests[!names(ests)%in%c("z.ac.term", "log.theta")]
			z.pred.data[, names(ests), drop=F]%*%ests
		}else{
			rep(NA, nrow(z.pred.data))
		}
	})
	all.preds=matrix(unlist(all.preds), ncol=length(all.preds), byrow=F)
	to.cut=exp(all.preds)>thresh.density
	to.cut[is.na(to.cut)]=F
	all.preds[to.cut]=log(thresh.density)
	to.cut.ind.mod=table(to.cut)
	#test.mean.ab[[i]]=xx
	all.preds=all.preds[, models2keep]
	all.preds=all.preds%*%c.set$w.aic ####predictions per model for all data sets
	all.preds=exp(all.preds+log(pred.env.fa$area.size))
	to.cut=all.preds/pred.env.fa$area.size>thresh.density
	all.preds[to.cut]=thresh.density*pred.env.fa$area.size[to.cut]
	all.pred[[i]]=all.preds
	kept.and.discarded[[i]]=list(to.cut.ind.mod=to.cut.ind.mod, to.cut.avg=table(to.cut))
}

#Get range of aicc-values: (corrected for small sample)
xx.aic=vector("list", length(all.res))
for(i in 1:length(all.res)){
	print(i)
	load(flist[i])
	#mmi:
	all.ll=unlist(lapply(all.m.res, function(x){x$logLik}))
	all.k=unlist(lapply(all.m.res, function(x){length(x$coef)}))
	all.aicc=-2*all.ll+2*all.k+aic.c.fac(N=sample.size, k=all.k)
	all.aicc=all.aicc[!is.na(all.aicc)]
	xx.aic[[i]]=all.aicc
}
range(unlist(xx.aic))
xx=unlist(lapply(all.pred, sum, na.rm=T))
paste(round(range(xx)), collapse=" to ")

###Extract coefficients
#Mean across 100 datasets
xx=c()
for(i in 1:length(all.res)){
	yy=data.frame(term=rownames(all.res[[i]]$full.coeffs), all.res[[i]]$full.coeffs)
	xx=rbind(xx, yy)
}
xx=xx[order(xx$term), ]
to.plot=xx
to.plot=droplevels(subset(to.plot, term!="theta"))

###Estimate variation across 100 datasets
xx=mapply(FUN=tapply, X=to.plot[, -1], MoreArgs=list(INDEX=to.plot$term, FUN=mean))
round(xx, 3)

#Range across 100 datasets (model coefficients)
xx=c()
for(i in 1:length(all.res)){
	yy=data.frame(term=rownames(all.res[[i]]$full.coeffs), all.res[[i]]$full.coeffs)
	xx=rbind(xx, yy)
}
xx=xx[order(xx$term), ]
maxs=mapply(FUN=tapply, X=xx[, -1], MoreArgs=list(INDEX=xx[, 1], FUN=max))
mins=mapply(FUN=tapply, X=xx[, -1], MoreArgs=list(INDEX=xx[, 1], FUN=min))
xx=paste(round(mins, 3), round(maxs, 3), sep=" to ")
xx=matrix(unlist(xx), ncol=4, byrow=F)
colnames(xx)=colnames(maxs)
rownames(xx)=rownames(maxs)
round(xx, 3)

#Sample size
xx=unlist(lapply(all.res, function(x){nrow(x$model.data)}))
paste(range(xx), collapse=" to ")
#List all countries 
all.countries=lapply(all.res, function(x){
	x$model.data$country
})
all.countries=sort(unique(unlist(all.countries)))

#Predictions per country, mean and range across 100 datasets
pred.env.fa=droplevels(pred.env.fa)
xx=lapply(all.pred, function(x){tapply(x, pred.env.fa$country, sum, na.rm=T)})
pred.per.country=matrix(unlist(xx), ncol=length(all.pred), byrow=F)
rownames(pred.per.country)=levels(pred.env.fa$country)
round(apply(pred.per.country, 1, mean, na.rm=T))
round(apply(pred.per.country, 1, mean, na.rm=T))
ound(apply(pred.per.country, 1, range, na.rm=T))

##### select dataset to.show
#Get nationwide predictions for comparison with published estimates:
xx=lapply(all.pred, function(x){tapply(x, pred.env.fa$country, sum, na.rm=T)})
pred.per.country=matrix(unlist(xx), ncol=length(all.pred), byrow=F)
rownames(pred.per.country)=levels(pred.env.fa$country)
apes.per.country=read.table(file="Y:/primint/Isabel/ape_density_distribution/R/Model_april2020/apes_per_country.txt", header=T, sep="\t") ### info. included in Supplementary 1, Table S3.4
setdiff(levels(apes.per.country$country), rownames(pred.per.country)) 
setdiff(rownames(pred.per.country), levels(apes.per.country$country))
rownames(pred.per.country)
levels(apes.per.country$country)
rownames(pred.per.country)[rownames(pred.per.country)=="Angola"]="AGO"
rownames(pred.per.country)[rownames(pred.per.country)=="EquatorialGuinea"]="GNQ"
rownames(pred.per.country)[rownames(pred.per.country)=="Burundi"]="BDI"
rownames(pred.per.country)[rownames(pred.per.country)=="GuineaBissau"]="GIB"
rownames(pred.per.country)[rownames(pred.per.country)=="MLI"]="MAL"
rownames(pred.per.country)[rownames(pred.per.country)=="SouthSudan"]="SSD"
intersect(rownames(pred.per.country), levels(apes.per.country$country))
countries.2.keep=intersect(levels(apes.per.country$country), rownames(pred.per.country))
apes.per.country=droplevels(subset(apes.per.country, country%in%countries.2.keep))
pred.per.country=pred.per.country[rownames(pred.per.country)%in%countries.2.keep, ]

###Select data set.
rel.change=(pred.per.country[as.character(apes.per.country$country), ]-apes.per.country$estimate)/apes.per.country$estimate 
to.show=which.min(apply(abs(rel.change), 2, sum, na.rm=T)) #relative error (gives more weight to larger relative deviations)
round(sum(all.pred[[to.show]], na.rm=T))

###Extract names of sites actually used:
all.sites.used=lapply(all.res, function(x){
	x$model.data$site_name[!x$model.data$site_name%in%fict.areas]
})
all.sites.used=sort(unique(as.character(unlist(all.sites.used))))
#Difference between all.sites.used and all real survey sites
sites.not.used=setdiff(levels(abundance.data$site_name), all.sites.used)
xx=as.character(abundance.data$site_name[is.na(abundance.data$year_of_es) & abundance.data$abund_esti==0])
sites.not.used=setdiff(sites.not.used, xx)
####Pseudo-absence sites and real survey sites
all.sites=lapply(all.res, function(x){
	x$model.data$site_name
})
all.sites=sort(unique(as.character(unlist(all.sites))))
fict.sites.used=lapply(all.res, function(x){
	x$model.data$site_name[x$model.data$site_name%in%fict.areas] 
})
fict.sites.used=sort(unique((unlist(fict.sites.used))))
length(all.sites.used)#156
length(all.sites)#290
length(fict.sites.used)#134
all.sites.used

#Predictions per country for selected dataset
round(pred.per.country[, to.show])
round(pred.per.country[, to.show])

#Extract full null model comp. for selected dataset:
round(all.res[[to.show]]$fn.comp, 3)

#Plot prediction sum against average sampling year:
xx=unlist(lapply(all.pred, sum, na.rm=T))
yy=unlist(lapply(all.years, mean))
plot(yy, xx, xlab="mean sampling year", ylab="abundance estimate")
points(yy[to.show], xx[to.show], col="red", pch=19)
sum(all.pred[[to.show]], na.rm=T)

###Bootstrap predicted abundance for the selected random assignment (to.show):
#Prepare prediction data:
model.data=all.res[[to.show]]$model.data
pred.env.fa=prediction.data
pred.env.fa=subset(pred.env.fa, complete.cases(pred.env.fa[, c(test.vars)]))
z.pred.data=data.frame(scale(pred.env.fa[, test.vars], center=apply(model.data[, test.vars], 2, mean), 
	scale=apply(model.data[, test.vars], 2, sd)))
colnames(z.pred.data)=paste("z", colnames(z.pred.data), sep=".")
z.pred.data=model.matrix(object=~z.log.human.foot.prt+z.sqrt.elevation+z.log.GDP.const+I(z.log.GDP.const^2)+z.prop.ifl.cl.yr+z.log.corruption+
	z.log.prop.muslims+z.log.precip.driest.month+z.survey.year.r, data=as.data.frame(z.pred.data))
	
#Prepare parallelization:
library(parallel)
cl <- makeCluster(getOption("cl.cores", detectCores()))
cl=cl[1:(length(cl)-1)]

##Get model data
model.data=all.res[[to.show]]$model.data
sample.size=nrow(all.res[[i]]$model.data)
to.sel.from=list((1:nrow(model.data))[model.data$site_name%in%fict.areas], (1:nrow(model.data))[!model.data$site_name%in%fict.areas])

print("bootstrap abundance per cell")
n.boots=1000
all.pred.boot=vector("list", n.boots)
for(i in 1:n.boots){
	set.seed(i)
	print(i)
	xx=c(sample(x=to.sel.from[[1]], size=length(to.sel.from[[1]]), replace=T), sample(x=to.sel.from[[2]], size=length(to.sel.from[[2]]), replace=T))
	boot.data=model.data[xx, ]#
	clusterExport(cl=cl, varlist=c("boot.data", "fit.fun"))
	all.m.res=parLapply(cl=cl, X=allmodels, fun=function(x){
		full=try(fit.fun(formula=x, data=boot.data, offset=boot.data$area.size.km, tol=0.00001), silent=T)
		if(class(full)[[1]]!="try-error"){
			ests=full$coefs[, "Estimate"]
			names(ests)=rownames(full$coefs)
			logLik=full$logLik
			k=length(full$coefs[, "Estimate"])
		}else{
			ests=NA
			logLik=NA
			k=NA
		}
		return(list(coef=ests, logLik=logLik, k=k))
	})
	#Get confidence set:
	all.ll=unlist(lapply(all.m.res, function(x){x$logLik}))
	all.k=unlist(lapply(all.m.res, function(x){length(x$coef)}))
	all.aicc=-2*all.ll+2*all.k+aic.c.fac(N=sample.size, k=all.k)
	models2keep=!is.na(all.aicc)
	models.conv=allmodels[!is.na(all.aicc)]
	all.m.res=all.m.res[!is.na(all.aicc)]
	all.aicc=all.aicc[!is.na(all.aicc)]
	c.set=conf.set(aic=all.aicc)
	#Get predictions for individual models:
	i.all.preds=lapply(all.m.res, function(x){
		ests=x$coef
		if(length(ests)>1){
			ests=ests[!names(ests)%in%c("log.theta", "z.ac.term")]
			z.pred.data[, names(ests), drop=F]%*%ests
		}else{
			rep(NA, nrow(z.pred.data))
		}
	})
	i.all.preds=matrix(unlist(i.all.preds), ncol=length(i.all.preds), byrow=F)
	i.all.preds=exp(i.all.preds+log(pred.env.fa$area.size))
	i.all.preds=i.all.preds%*%c.set$w.aic
	all.pred.boot[[i]]=i.all.preds	
}
xx=unlist(lapply(all.pred.boot, sum, na.rm=T))
which.max(xx) #262
plot(xx, log="y")
round(quantile(xx, prob=c(0.025, 0.975)))

par(mfrow=c(2, 1))
plot(table(boot.data$abund_esti))
plot(table(model.data$abund_esti))
exp(mean(log(model.data$abund_esti+1)))-1
exp(mean(log(boot.data$abund_esti+1)))-1
min(model.data$abund_esti[!model.data$site_name%in%fict.areas])
model.data$site_name[model.data$abund_esti==0 & !model.data$site_name%in%fict.areas]

#Evaluate bootstraps per cell in the prediction grid:
pred.env.fa
xx=matrix(unlist(all.pred.boot), ncol=length(all.pred.boot), byrow=F)
ci.per.cell=apply(xx, 1, quantile, prob=c(0.025, 0.975), na.rm=T)
ci.per.cell=t(ci.per.cell)
plot(ci.per.cell[, 2])
table(round(ci.per.cell[, 2]))
mean(ci.per.cell[, 2]>100)

all.abes=rep(NA, n.boots)
for(i in 1:n.boots){
	set.seed(i)
	xx=c(sample(x=to.sel.from[[1]], size=length(to.sel.from[[1]]), replace=T), sample(x=to.sel.from[[2]], size=length(to.sel.from[[2]]), replace=T))
	boot.data=model.data[xx, ]#sample(x=to.sel.from, size=length(to.sel.from), replace=T), ]
	all.abes[i]=mean(boot.data$survey.year.r)
}
xx=unlist(lapply(all.pred.boot, sum, na.rm=T))
cor(all.abes, xx, method="spearman") #
plot(all.abes, xx, log="y")
abline(v=mean(model.data$survey.year.r))

###Confidence intervals for nationwide estimates
cl.country=as.character(pred.env.fa$country)
lowcl=as.numeric(ci.per.cell[1:32428])
lowcl=as.numeric(ci.per.cell[,1])
upcl=as.numeric(ci.per.cell[32429:64856])
upcl=as.numeric(ci.per.cell[,2])
x<-pred.env.fa$lat
y<-pred.env.fa$long
cl.data=data.frame(lat=x, long=y, lowcl, upcl, cl.country=as.character(pred.env.fa$country),stringsAsFactors=F)
aggregate(cl.data$lowcl, by=list(country=cl.data$cl.country), FUN=sum)
group_by(cl.data,cl.country)%>%summary(sum_lowcl=sum(lowcl,na.rm=T),sum_upcl=sum(upcl,na.rm=T))

##Evaluation: plotting variation of predictors in model data vs. predictor data
xx=allmodels[length(allmodels)]
xx=gsub(xx, pattern="abund_esti~1+", replacement="", fixed=T)
xx=unlist(strsplit(xx, split="+", fixed=T))
xx=xx[!grepl(xx, pattern="(", fixed=T)]
par(mfrow=c(2, 4), mar=c(3, 3, 1.5, 0.2))
for(i in 1:length(xx)){
	d1=density(model.data[, xx[i]])
	d2=density(z.pred.data[, xx[i]])
	plot(x=d1$x, y=d1$y/max(d1$y), type="l", xlim=range(c(d1$x, d2$x)), col="red")
	lines(x=d2$x, y=d2$y/max(d2$y), col="blue")
	mtext(text=xx[i], side=3, line=0.2)
}
legend("topright", legend=c("model data", "prediction data"), text.col=c("red", "blue"), bty="n")

###Bootstrap model coefficients:
library(parallel)
cl <- makeCluster(getOption("cl.cores", detectCores()))
print("bootstrap model coefficients")
n.boots=1000
all.coef.boots=vector("list", length(all.res))
for(i in 1:length(all.res)){
	print(i)
	model.data=all.res[[i]]$model.data
	to.sel.from=1:nrow(model.data)
	clusterExport(cl=cl, varlist=c("model.data", "to.sel.from", "fit.fun"))
	iall.boots=parLapply(cl=cl, X=1:n.boots, fun=function(x){
	#iall.boots=lapply(X=1:n.boots, FUN=function(x){
		#print(x)
		xdone=FALSE
		while(!xdone){
			boot.data=model.data[sample(x=to.sel.from, size=length(to.sel.from), replace=T), ]
			ifull=try(fit.fun(formula="abund_esti~1+z.log.human.foot.prt+z.sqrt.elevation+z.log.GDP.const+I(z.log.GDP.const^2)+z.prop.ifl.cl.yr+z.log.corruption+z.log.prop.muslims+z.log.precip.driest.month+z.survey.year.r+z.ac.term",
				data=boot.data, offset=model.data$area.size.km, tol=0.00001), silent=T)
			if(class(ifull)[1]!="try-error"){
				xdone=TRUE
				return(ifull$coefs)
			}
		}
	})
	#iall.boots=matrix(unlist(iall.boots), nrow=length(iall.boots), byrow=T)
	xx=lapply(iall.boots, function(x){
		return(x[, "Estimate"])
	})
	xx=matrix(unlist(xx), nrow=length(xx), byrow=T)
	colnames(xx)=rownames(iall.boots[[1]])
	all.coef.boots[[i]]=apply(xx, 2, quantile, probs=c(0.025, 0.975))
}
xx=lapply(all.coef.boots, function(x){x["2.5%", ]})
xx=matrix(unlist(xx), nrow=length(xx), byrow=T)
colnames(xx)=colnames(all.coef.boots[[1]])
all.ci=apply(xx, 2, mean)
all.ci=cbind(all.ci, t(apply(xx, 2, range)))
colnames(all.ci)=c("lwr.mean", "lwr.min", "lwr.max")
xx=lapply(all.coef.boots, function(x){x["97.5%", ]})
xx=matrix(unlist(xx), nrow=length(xx), byrow=T)
colnames(xx)=colnames(all.coef.boots[[1]])
all.ci=cbind(all.ci, "upr.mean"=apply(xx, 2, mean))
all.ci=cbind(all.ci, "upr.min"=apply(xx, 2, min))
all.ci=cbind(all.ci, "upr.max"=apply(xx, 2, max))
round(all.ci, 3)

###################################################################################################################################################
#### Get results from models (all.res) i.e. across 100 data sets
#Sample size:
xx=unlist(lapply(all.res, function(x){nrow(x$model.data)}))
median(xx)
paste(range(xx), collapse=" to ")
#Dispersion parameters:
round(mean(unlist(lapply(all.res, function(x){x$disp.par}))), 3)
paste(round(range(unlist(lapply(all.res, function(x){x$disp.par}))), 3), collapse=" to ")
#Sample size across data sets
unlist(lapply(all.res, function(x){x$n}))
#Full null model comparison:
xx=matrix(unlist(lapply(all.res, function(x){x$fn.comp})), nrow=length(all.res), byrow=T)
colnames(xx)=names(all.res[[1]]$fn.comp)
paste(round(apply(xx, 2, mean), 3), collapse=", ")#28.963, 6, 0
(apply(xx, 2, max))[3] #9.95882593532986e-05

###Plot showing range of estimates:
dev.off()
windows(height=5)
par(mfrow=c(1,1), mar=c(3, 7, 0.2, 0.5), mgp=c(1.7, 0.3, 0), tcl=-0.15)
plot(x=to.plot$Estimate, y=as.numeric(to.plot$term), yaxt="n", ylab="", xlab="estimate", pch=19, col=grey(alpha=0.2, level=0.2))
to.add=levels(to.plot$term)
to.add=gsub(to.add, pattern="log.", replacement="", fixed=T)
to.add=gsub(to.add, pattern="sqrt.", replacement="", fixed=T)
to.add=gsub(to.add, pattern="z.", replacement="", fixed=T)
to.add=gsub(to.add, pattern="I(", replacement="", fixed=T)
to.add=gsub(to.add, pattern=")", replacement="", fixed=T)
to.add=gsub(to.add, pattern="(", replacement="", fixed=T)
to.add=gsub(to.add, pattern="GDP.const^2", replacement="GDP", fixed=T)
to.add[to.add=="GDP.const"]="GDP"
to.add=gsub(to.add, pattern="survey.year.r", replacement="survey year", fixed=T)
to.add=gsub(to.add, pattern="prop.muslims", replacement="food taboos", fixed=T)
to.add[to.add=="precip.driest.month"]="min. precipitation"
to.add[to.add=="prop.ifl.cl.yr"]="proportion IFL"
#to.add=gsub(to.add, pattern="prop.ifl.cl.yr", replacement="proportion IFL", fixed=T)
to.add=gsub(to.add, pattern="human.foot.prt", replacement="human footprint", fixed=T)
#to.add=gsub(to.add, pattern="(Intercept", replacement="intercept",  fixed=T)
to.add[to.add=="Intercept"]="intercept"
to.add[to.add=="ac.term"]="autocorrelation"
#to.add=tolower(to.add) #
#to.add[to.add=="GDP.const"]="GDP"
mtext(text=bquote(.(tolower(to.add[12]))), side=2, at=12, las=1, line=0.2, cex=0.9)
mtext(text=to.add[3:12], side=2, at=3:12, las=1, line=0.2, cex=0.9)
#mtext(text=to.add[3:length(to.add)], side=2, at=3:length(levels(to.plot$term)), las=1, line=0.2, cex=0.9)
mtext(text=bquote(.(to.add[2])^2), side=2, at=2, las=1, line=0.2, cex=0.9)
mtext(text=bquote(.(to.add[1])), side=2, at=1, las=1, line=0.2, cex=0.9)
abline(v=0, lty=3)
       
###Estimate variation
xx=mapply(FUN=tapply, X=xx[, -1], MoreArgs=list(INDEX=xx[, 1], FUN=mean))
round(xx, 3)
xx=c()
for(i in 1:length(all.res)){
	yy=data.frame(term=rownames(all.res[[i]]$full.coeffs), all.res[[i]]$full.coeffs)
	xx=rbind(xx, yy)
}
xx=xx[order(xx$term), ]
mins=mapply(FUN=tapply, X=xx[, -1], MoreArgs=list(INDEX=xx[, 1], FUN=min))
maxs=mapply(FUN=tapply, X=xx[, -1], MoreArgs=list(INDEX=xx[, 1], FUN=max))
xx=paste(round(mins, 3), round(maxs, 3), sep=" to ")
xx=matrix(unlist(xx), ncol=4, byrow=F)
colnames(xx)=colnames(maxs)
rownames(xx)=rownames(maxs)
xx

###Model stability:
print("model stability")
all.m.stab=vector("list", length(all.res))
for(i in 1:length(all.res)){
	print(i)
	model.data=all.res[[i]]$model.data
	i.stab=matrix(NA, nrow=nrow(model.data), ncol=nrow(all.res[[i]]$full.coeffs))
	colnames(i.stab)=rownames(all.res[[i]]$full.coeffs)
	for(j in 1:nrow(model.data)){
		ifull=try(fit.fun(formula="abund_esti~1+z.log.human.foot.prt+z.sqrt.elevation+z.log.GDP.const+I(z.log.GDP.const^2)+z.prop.ifl.cl.yr+z.log.corruption+z.log.prop.muslims+z.log.precip.driest.month+z.survey.year.r+z.ac.term",
			data=model.data[-j, ], offset=model.data[-j, ]$area.size.km, tol=0.00001), silent=T)
		if(class(ifull)!="try-error"){
			i.stab[j, ]=ifull$coefs[, "Estimate"]
			colnames(i.stab)=rownames(ifull$coefs)
		}
	}
	n=sum(!is.na(i.stab[, 1]))
	#colnames(i.stab)=rownames(
	i.stab=cbind(all.res[[i]]$full.coeffs[, "Estimate"], t(apply(i.stab, 2, range, na.rm=T))[rownames(all.res[[i]]$full.coeffs), ])
	colnames(i.stab)=c("orig", "min", "max")
	all.m.stab[[i]]=list(m.stab=i.stab, n=n)
}
xx=c()
for(i in 1:length(all.m.stab)){
	x=all.m.stab[[i]]$m.stab
	x=data.frame(term=rownames(x), x)
	xx=rbind(xx, x)
}
#xx=mapply(FUN=tapply, X=xx[, -1], MoreArgs=list(INDEX=xx[, 1], FUN=mean))
all.stabs=aggregate(xx[, "orig"], list(term=xx[, "term"]), mean)
names(all.stabs)[2]="orig"
all.stabs$min=aggregate(xx[, "min"], list(term=xx[, "term"]), min)$x
all.stabs$max=aggregate(xx[, "max"], list(term=xx[, "term"]), max)$x
all.stabs
rownames(all.stabs)=all.stabs$term
par(mfrow=c(1, 1))
m.stab.plot(all.stabs[, -1])

###Extract akaike weights per predictor, and other stuff:
flist=paste("Y:/primint/Isabel/ape_density_distribution/R/Model_june2020/MMI_ifl_humaninf_ele_gdpconst_corr_minprec_foodtab_june20/test_mmi_.", 1:100, "..RData", sep="") 
all.mmi.res=vector("list", length(flist))
for(i in 1:length(flist)){
	load(flist[i])
	all.ll=unlist(lapply(all.m.res, function(x){x$logLik}))
	all.k=unlist(lapply(all.m.res, function(x){length(x$coef)}))
	all.aicc=-2*all.ll+2*all.k+aic.c.fac(N=sample.size, k=all.k)
	models.conv=allmodels[!is.na(all.aicc)]
	all.aicc=all.aicc[!is.na(all.aicc)]
	c.set=conf.set(aic=all.aicc)
	x=c.set$w.aic
	names(x)=gsub(x=models.conv, pattern="abund_esti~", replacement="", fixed=T)
	all.mmi.res[[i]]=list(
		#ak.wgt.per.pv=ak.wgt.per.pv,
		n.in.cset=sum(c.set$c.set), 
		null.in.cset=c.set$c.set[1]==1 & models.conv[1]==allmodels[1],
		n.conv=sum(!is.na(all.ll))
		#exp.waic.per.pv=exp.waic.per.pv
	)
}

###MMI, converged models, mean and range, etc
round(mean(unlist(lapply(all.mmi.res, function(x){x$n.conv}))), 1)
paste(range(unlist(lapply(all.mmi.res, function(x){x$n.conv}))), collapse=" to ")
length(allmodels)
#Number models in c.set:
round(mean(unlist(lapply(all.mmi.res, function(x){x$n.in.cset}))), 1)
round(mean(100*unlist(lapply(all.mmi.res, function(x){x$n.in.cset}))/unlist(lapply(all.mmi.res, function(x){x$n.conv}))), 1)
paste(range(unlist(lapply(all.mmi.res, function(x){x$n.in.cset}))), collapse=" to ")
#Range number models in conf. set:
paste(round(range(100*unlist(lapply(all.mmi.res, function(x){x$n.in.cset}))/
	unlist(lapply(all.mmi.res, function(x){x$n.conv}))), 1), collapse=" to ")
table(unlist(lapply(all.mmi.res, function(x){x$null.in.cset})))
#=> Null model never in conf. set across the mmis of the 100 datasets

##VIF:
library(car)
all.vif=lapply(all.res, function(x){
	xx=lm(abund_esti~1+z.log.human.foot.prt+z.sqrt.elevation+z.log.GDP.const+I(z.log.GDP.const^2)+z.prop.ifl.cl.yr+z.log.corruption+z.log.prop.muslims+z.log.precip.driest.month+z.survey.year.r, data=x$model.data)
	vif(xx)
})
max(unlist(all.vif))

#Means and SDs model data
pvs=c("sqrt.elevation", "log.human.foot.prt", "log.GDP.const", "prop.ifl.cl.yr", "log.prop.muslims", "log.corruption", 
		"log.precip.driest.month", "survey.year.r")
all.means=c()
all.sds=c()
for(i in 1:length(all.res)){
	model.data=all.res[[i]]$model.data
	all.means=rbind(all.means, apply(model.data[, pvs], 2, mean))
	all.sds=rbind(all.sds, apply(model.data[, pvs], 2, sd))
}
round(cbind(apply(all.means, 2, mean), apply(all.sds, 2, mean)),3)
round(cbind(apply(all.means, 2, mean), apply(all.sds, 2, mean)),3)

#Means and sds for prediction data
pvs.pred=c("sqrt.elevation", "log.human.foot.prt", "log.GDP.const", "prop.ifl.cl.yr", "log.prop.muslims", "log.corruption", 
		"log.precip.driest.month", "survey.year.r")
apply(pred.env.fa[, pvs.pred], 2, mean)
apply(pred.env.fa[, pvs.pred], 2, sd)

#Mean, median, min, max nationwide predicted estimate across 100 datasets
round(apply(pred.per.country, 1, mean, na.rm=T))
round(apply(pred.per.country, 1, median, na.rm=T))
round(apply(pred.per.country, 1, min, na.rm=T))
round(apply(pred.per.country, 1, max, na.rm=T))
#Mean and median survey year, all sites, pesudo-absence sites, real sites
round(mean(sel.abundance.data$survey.year.r[as.character(sel.abundance.data$site_name)%in%all.sites.used]))
round(mean(sel.abundance.data$survey.year.r[as.character(sel.abundance.data$site_name)%in%fict.sites.used]))
round(mean(sel.abundance.data$survey.year.r))
round(median(sel.abundance.data$survey.year.r[as.character(sel.abundance.data$site_name)%in%all.sites.used]))
round(median(sel.abundance.data$survey.year.r[as.character(sel.abundance.data$site_name)%in%fict.sites.used]))

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

####Extract results for to.show (selected dataset)
#Model coeffs
all.res[[to.show]]$full.coeffs
round(all.res[[to.show]]$full.coeffs, 3)
#Conf int for model coeffs
t(all.coef.boots[[to.show]])
#Model stability
round(all.m.stab[[to.show]]$m.stab,3)
#Likelihood ratio test from full-null model comparison
all.res[[to.show]]$fn.comp
#Dispersion parameter
round(all.res[[to.show]]$disp.par, 3)
#Sample size
all.res[[to.show]]$n
#MMI models
range(all.mmi.res[[to.show]])
#Collinearity
max(all.vif[[to.show]]) 
#Models in c.set:
all.mmi.res[[to.show]]$n.in.cset
#Range of predictions per km2
xx=as.vector(all.pred[[to.show]])/pred.env.fa$area.size
round(range(xx, na.rm=T),3)

#Plotting human footprint for each country
pred.env.fa$country_names=pred.env.fa$country
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="CAF", replacement="Central African Rep.", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="CIV", replacement="Côte d'Ivoire", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="CMR", replacement="Cameroon", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="COD", replacement="DRC", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="COG", replacement="Rep. of Congo", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="GAB", replacement="Gabon", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="GHA", replacement="Ghana", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="GIN", replacement="Guinea", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="LBR", replacement="Liberia", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="MLI", replacement="Mali", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="NGA", replacement="Nigeria", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="RWA", replacement="Rwanda", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="SEN", replacement="Senegal", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="SLE", replacement="Sierra Leone", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="TZA", replacement="Tanzania", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="UGA", replacement="Uganda", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="GuineaBissau", replacement="Guinea Bissau", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="EquatorialGuinea", replacement="Equatorial Guinea", fixed=T)
pred.env.fa$country_names=gsub(x=pred.env.fa$country_names, pattern="SouthSudan", replacement="South Sudan", fixed=T)
pred.env.fa$country_names=as.factor(pred.env.fa$country_names)
cc=unique(pred.env.fa$country_names)
countries=as.factor(cc)
par(mar=c(6.6, 2.5, 0.2, 0.2), mgp=c(2.5, 0.2, 0), tcl=-0.15)
plot(x=pred.env.fa$country_names, y=pred.env.fa$log.human.foot.prt, xlab="", ylab="", cex=0.8, 
 asp=1, xaxt="n")
mtext(text="Human footprint", side=2, line=1.5, cex=0.8, las=0)	
mtext(text=levels(countries), side=1, line=0.2, at=1:21, cex=0.8, las=2)	
range(pred.env.fa$log.human.foot.prt, na.rm=T)

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

###Plotting
#Significant vars, Human Footprint and GDP
par(mfrow=c(1, 2))
	### Human footprint
	layout(matrix(c(3, 1:2), nrow=1), widths=c(0.2, rep(5, 3)))
	layout.show(3)
	icoefs=all.res[[to.show]]$full.coeffs[, "Estimate"]
	names(icoefs)=rownames(all.res[[to.show]]$full.coeffs)
	range(all.res[[to.show]]$model.data$human.footprint)
	x.labels=2^(0:6)
	x.at=(log(x.labels)-mean(log(all.res[[to.show]]$model.data$human.footprint)))/sd(log(all.res[[to.show]]$model.data$human.footprint))
	xx=no.int.plot(
		plot.data=all.res[[to.show]]$model.data, covariate="z.log.human.foot.prt", coefs=icoefs, 
		response=all.res[[to.show]]$model.data$abund_esti/all.res[[to.show]]$model.data$area.size.km, 
		link="log",
		grid.resol=20, 
		xlab="Human Footprint", ylab="",
		x.at=x.at, x.labels=x.labels, y.at=c(0,1,2), y.labels=c(0,1,2), xlim=NULL, ylim=c(0,2.5),
		cex.lab=1.2, cex.axis=1.2, size.fac=1, col="black", pch=1, lty=2,
		my.par=list(mar=c(3.4, 2.8, 0.4, 0.8), mgp=c(2, 0.8, 0), tcl=-0.15, las=1), quiet=F, reset.par=F, mean.fun=median)
	mtext(text=bquote(paste("Great ape density ", group("(", ind/km^2, ")"))), side=2, line=1.8, las=0, cex=0.8)
	#mtext(text=bquote(paste("great ape density ", group("[", frac(individuals, km^2), "]"))), side=2, line=2.7, las=0, cex=0.8)
	range(xx$N)
	legend(x="topleft", legend="A", bty="n", cex=1)
	### GDP:
	range(all.res[[to.show]]$model.data$GDP.const)
	x.labels=2^(0:9)
	x.at=(log(x.labels*1000000000)-mean(log(all.res[[to.show]]$model.data$GDP.const)))/sd(log(all.res[[to.show]]$model.data$GDP.const))
	xx=no.int.plot(
		plot.data=all.res[[to.show]]$model.data, covariate="z.log.GDP.const", coefs=icoefs, 
		response=all.res[[to.show]]$model.data$abund_esti/all.res[[to.show]]$model.data$area.size.km, 
		link="log",
		grid.resol=20, 
		xlab="Gross Domestic Product (billion USD)", ylab="",
		x.at=x.at, x.labels=x.labels, y.at=NULL, y.labels=NULL, xlim=NULL, ylim=c(0, 0.5),
		cex.lab=1.2, cex.axis=1.2, size.fac=1, col="black", pch=1,lty=2,
		my.par=list(tcl=-0.15, las=1), quiet=F, reset.par=F, mean.fun=median)
		range(xx$N) # 1-69
		legend(x="topleft", legend="B", bty="n", cex=1)

###Import great ape range polygon for the density distribution map
ga.range=read.shp("Y:/primint/Isabel/ape_density_distribution/R/new_ga_range/all_apes_with_new_ptverus.shp")$shp
ga.range=split.shape(shp=ga.range)
io=lapply(1:length(ga.range), function(x){
	print(x)
	x=polyoffset(A=list(list(x=ga.range[[x]]$X, y=ga.range[[x]]$Y)), 1)
	inout(pts=cbind(x=abundance.data$X, y=abundance.data$Y), poly=cbind(x=x[[1]]$x, y=x[[1]]$y), bound=T)
})
io=matrix(unlist(io), ncol=length(io), byrow=F)
io=apply(io, 1, any)

###Dissolved great ape range (no borders inside)
ga.range=read.shp("Y:/primint/Isabel/ape_density_distribution/R/new_ga_range/dissolved_range_new_ptverus.shp")$shp
ga.range=split.shape(shp=ga.range)
io=lapply(1:length(ga.range), function(x){
	print(x)
	x=polyoffset(A=list(list(x=ga.range[[x]]$X, y=ga.range[[x]]$Y)), 1)
	inout(pts=cbind(x=abundance.data$X, y=abundance.data$Y), poly=cbind(x=x[[1]]$x, y=x[[1]]$y), bound=T)
})
io=matrix(unlist(io), ncol=length(io), byrow=F)
io=apply(io, 1, any)

####Plot range-wide prediction
par(mfrow=c(1, 1))
library(maps)
#Attach predicted dens to pred.env.fa dataframe
pred.env.fa.1=pred.env.fa
pred.env.fa.1$predicted.ab=all.pred[[to.show]]
#Subset pred.env.fa (exclude DRC, LBR, SSD, CAR)
pred.env.fa.sub=droplevels(subset(pred.env.fa.1, country%in%c("Angola", "Burundi",  "CMR", "COG","EquatorialGuinea",
 "GAB",  "GHA",  "GIN",   "GuineaBissau",  "MLI",  "NGA", "RWA", "SEN", "SLE",  "TZA",  "UGA")))
ab.to.plot=as.vector(pred.env.fa.sub$predicted.ab/pred.env.fa.sub$area.size)
par(mar=c(2.4, 2, 0.2, 0.2), mgp=c(1.7, 0.3, 0), tcl=-0.15, las=1)
plot(x=pred.env.fa.sub$long, y=pred.env.fa.sub$lat, asp=1, type="n", xlab="", ylab="")
mtext(text="longitude", side=1, line=1.2, cex=0.8)
mtext(text="latitude", side=2, line=1, cex=0.8, las=0)	
hrw=0.5*pred.data.cell.width #width of cell in degrees
map(add=T, fill=T, col="snow3", boundary=F)
for(i in 1:length(ga.range)){
	polygon(ga.range[[i]], col="white", border=NA, lwd=2) ###make countries with no pred to show blank
}
col=(ab.to.plot-min(ab.to.plot, na.rm=T))/diff(range(ab.to.plot, na.rm=T))
ramp=colorRamp(colors=c("azure2","gold", "yellowgreen", "forestgreen", "darkgreen", "deepskyblue4")) 
col.val=rgb(red=ramp(col[!is.na(ab.to.plot)]), maxColorValue=255)
rect(xleft=pred.env.fa.sub$long[!is.na(ab.to.plot)]-hrw, xright=pred.env.fa.sub$long[!is.na(ab.to.plot)]+hrw, 
	ybottom=pred.env.fa.sub$lat[!is.na(ab.to.plot)]-hrw, ytop=pred.env.fa.sub$lat[!is.na(ab.to.plot)]+hrw, border=col.val, col=col.val)	
map(add=T, boundary=T, col="gray31")
#Legend
xx=as.vector(pred.env.fa.sub$predicted.ab/pred.env.fa.sub$area.size)
range(xx, na.rm=T)
leg.vals=0:6
leg.n=400
segments(x0=min(pred.env.fa.sub$long), x1=min(pred.env.fa.sub$long)+diff(range(pred.env.fa.sub$long))*0.05,
	y0=seq(-8, 2, length.out=leg.n),
	y1=seq(-8, 2, length.out=leg.n),
	col=rgb(red=ramp(seq(0, 1, length.out=leg.n)), maxColorValue=255), lend=1)
text(labels=round(leg.vals, 6), x=min(pred.env.fa.sub$long)+diff(range(pred.env.fa.sub$long))*0.05,
	y=-8+10*(leg.vals-min(ab.to.plot, na.rm=T))/diff(range(ab.to.plot, na.rm=T)),
	pos=4, cex=0.8)
	text(labels=bquote(frac(individuals, km^2)), x=min(pred.env.fa.sub$long)+diff(range(pred.env.fa.sub$long))*0.025, y=2, pos=3, cex=0.8)
	for(i in 1:length(ga.range)){
	polygon(ga.range[[i]], col=NA, border="firebrick", lwd=2) ########### add new range mediumpurple3
}
box(lty="solid", col="black")
map.scale(x=-5, y=-7, relwidth = 0.1, metric = TRUE, ratio = FALSE, cex=0.7)
	
################Plot our estimates vs. estimates from other studies
par(mfrow=c(1, 1))
###country comparison
esti.data=read.table(file="Y:/primint/Isabel/ape_density_distribution/R/Model_june2020/country_comparison.txt", header=T, sep="\t") ## info. included in Supplementary 1, Table S3.4
str(esti.data)
esti.data=esti.data[complete.cases(esti.data), ]
hist(esti.data$our_study)
hist(log(esti.data$our_study))
hist(esti.data$other_studies)
hist(log(esti.data$other_studies))
esti.data$log.other.studies=log(esti.data$other_studies)
esti.data$log.our.study=log(esti.data$our_study)

plot(x=esti.data$other_studies, y=esti.data$our_study)
par(mar=c(3, 3.3, 0.2, 1.2), mgp=c(2.5, 0.2, 0), tcl=-0.15, las=1)
options(scipen=8)
plot(x=esti.data$log.other.studies, y=esti.data$log.our.study, pch=ifelse(esti.data$country=="Total", 19, 16), col=ifelse(esti.data$country=="Total", "red", "black"), xlab="", 
	ylab="", asp=1, xaxt="n", yaxt="n", cex=ifelse(esti.data$country=="Total", 1.4, 1),
	ylim=range(c(esti.data$log.other.studies, esti.data$log.our.study), na.rm=T), 
	xlim=range(c(esti.data$log.other.studies, esti.data$log.our.study), na.rm=T))
mtext(text="predicted nationwide estimates", side=2, line=2.4, cex=0.8, las=0)	
mtext(text="estimates from other studies", side=1, line=1.7, cex=0.8)
### estimates from other studies
pos=rep(4, nrow(esti.data))
pos[esti.data$country=="Democratic Republic of Congo"]=2
pos[esti.data$country=="Republic of Congo"]=2
pos[esti.data$country=="Central African Republic"]=2
pos[esti.data$country=="Burundi"]=1
pos[esti.data$country=="Ivory Coast"]=2
pos[esti.data$country=="Ghana"]=3
pos[esti.data$country=="Gabon"]=2
pos[esti.data$country=="Mali"]=2
pos[esti.data$country=="Tanzania"]=1
pos[esti.data$country=="Uganda"]=4
pos[esti.data$country=="Cameroon"]=2
pos[esti.data$country=="Senegal"]=2
pos[esti.data$country=="Rwanda"]=4
pos[esti.data$country=="Liberia"]=2
pos[esti.data$country=="Sierra Leone"]=2
pos[esti.data$country=="Guinea"]=2
pos[esti.data$country=="Total"]=2
text(x=esti.data$log.other.studies, y=esti.data$log.our.study, labels=esti.data$country, pos=pos, cex=0.8)
model=lm(esti.data$log.our.study~0+esti.data$log.other.studies)
#abline(model, lty=3)
abline(a=0, b=1, lty=2) #to compare line with actual data points and how much they diverge from ideal scenario
pretty(esti.data$our_study, na.rm=T)
range(esti.data$our_study)
xvals=c(sort(c(1*10^(1:5), 3*10^(0:5)))) #5*10^(1:5)
yvals=xvals
axis(side=1, at=log(xvals), labels=xvals, las=1, cex.axis=0.8)
axis(side=2, at=log(yvals), labels=yvals, las=1, cex.axis=0.8)

###Save dens dist layer to work on it in qgis
x<-pred.env.fa$lat
y<-pred.env.fa$long
z<-all.pred[[to.show]]
mt=cbind(x, y, z)
map.data=as.data.frame(mt)		 
names(map.data)[names(map.data) == "x"] <- "lat"
names(map.data)[names(map.data) == "y"] <- "long"	 
names(map.data)[names(map.data) == "z"] <- "ab"
write.table(x=map.data, file="Y:/primint/Isabel/ape_density_distribution/R/Model_june2020/ddm.txt", row.names=F, col.names=T, sep="\t")

###Plot pseudo-absence sites and real sites
library(maps)
library(rgdal)
par(mfrow=c(1, 1))
ab.to.plot=as.vector(all.pred[[to.show]])/pred.env.fa$area.size
par(mar=c(3, 2.8, 0.2, 0.2), mgp=c(1.7, 0.3, 0), tcl=-0.15, las=1)
plot(x=pred.env.fa$long, y=pred.env.fa$lat, asp=1, type="n", xlab="longitude", ylab="latitude")
map(add=T, fill=T, col="gray86")
map(add=T, col="black")	
for(i in 1:length(ga.range)){
	polygon(ga.range[[i]], col=NA, border="deepskyblue2", lwd=2)
}
to.add=(1:nrow(sel.abundance.data))[as.character(sel.abundance.data$site_name)%in%all.sites.used]
for(row.number in to.add){
	layer.name<-unlist(strsplit(as.character(sel.abundance.data$Path_polygons[row.number]), split="\\", fixed=T))
	layer.name<-layer.name[length(layer.name)]
	layer.name<-gsub(layer.name, pattern=".shp", replacement="", fixed=T)
	fname=sel.abundance.data$Path_polygons[row.number]
	#fname=gsub(x=fname, pattern="Y:", replacement="/home/roger_mundry/mnt/primatology", fixed=T)
	fname=gsub(x=fname, pattern="\\", replacement="/", fixed=T)
	test<-readOGR(dsn=fname, layer=layer.name) 
	new.shp=mk.list.ogr.shp(shp=test)
	#new.shp=new.shp[shapes.to.keep[[row.number]]]
	for(j in 1:length(new.shp)){
		polygon(new.shp[[j]], lwd=1, col="chartreuse2")
	}
}
to.add=(1:nrow(sel.abundance.data))[as.character(sel.abundance.data$site_name)%in%fict.sites.used]
for(row.number in to.add){
	layer.name<-unlist(strsplit(as.character(sel.abundance.data$Path_polygons[row.number]), split="\\", fixed=T))
	layer.name<-layer.name[length(layer.name)]
	layer.name<-gsub(layer.name, pattern=".shp", replacement="", fixed=T)
	fname=sel.abundance.data$Path_polygons[row.number]
	#fname=gsub(x=fname, pattern="Y:", replacement="/home/roger_mundry/mnt/primatology", fixed=T)
	fname=gsub(x=fname, pattern="\\", replacement="/", fixed=T)
	test<-readOGR(dsn=fname, layer=layer.name) 
	new.shp=mk.list.ogr.shp(shp=test)
	#new.shp=new.shp[shapes.to.keep[[row.number]]]
	for(j in 1:length(new.shp)){
		polygon(new.shp[[j]], lwd=1, col="darkgreen")
	}
}
box(lty="solid", col="black")
sites.legend=as.factor(c("real survey sites (n=156)", "fictitious sites (n=134)"))
legend(x="bottomleft", legend=levels(sites.legend), pch=16, cex=1, col=c("darkgreen", "chartreuse2"), bty = "n")

############################################## FUNCTIONS #############################################################

##Functions used for Ordaz-Nemeth et al.
##all functions written by Roger Mundry
#functions determining the grand cricle distance between two points of surface the earth
#arguments are the coordinates (N, E) of the two points (1, 2)
#coordinates need to be in wgs (longitude, latitude)
#both return the distance in meters
dist.wgs<-function(N1,E1,N2,E2){
	XX = sin(N1 * pi / 180) * sin(N2 * pi / 180) + cos(N1 * pi / 180) * cos(N2 * pi / 180) * cos(abs(E1 - E2) * pi / 180)
	if ((N1 != N2) | (E1 != E2)) {XX = atan(-XX / sqrt(-XX * XX + 1)) + 1.5707963267949} else {XX = 0}
	XX = XX * 180 / pi
	XX = XX * 60
	return(1852 * XX)
}
dist.wgs.vectorized<-function(N1, E1, N2v, E2v){
	XX = sin(N1 * pi / 180) * sin(N2v * pi / 180) + cos(N1 * pi / 180) * cos(N2v * pi / 180) * cos(abs(E1 - E2v) * pi / 180)
	sel=(N1 != N2v) | (E1 != E2v)
	XX[sel]=atan(-XX[sel] / sqrt(-XX[sel] * XX[sel] + 1)) + 1.5707963267949
	XX[!sel] = 0
	XX = XX * 180 / pi
	XX = XX * 60
	return(1852 * XX)
}

################################################################################
#function turning a shape object into a list with one data frame for each disjunkt part of a shape (makes dealing with shape 
	#files easier in which individual shapes consist of parts (e.g., borders of countries which comprise islands)
split.shape<-function(shp){
	xxx=c()
	lapply(shp, function(x){
		if(x$num.parts==1){
			xxx<<-c(xxx, list(x$points))
		}else{
			parts=rep(0, nrow(x$points))
			parts[x$parts+1]=1
			parts=cumsum(parts)
			lapply(unique(parts), function(ip){
				xxx<<-c(xxx, list(x$points[parts==ip,]))
			})
		}
	})
	return(xxx)
}

################################################################################
#function plotting estimates together with 'error bars' (e.g., CI or range of estimates after dropping cases,
	#one at a time)
m.stab.plot<-function(est, lower=NULL, upper=NULL, xnames=NULL, col="black", center.at.null=F, reset.par=T, pch=18, xlim=NULL){
	#version from Aug 12 2018
	#function to plot model stability or CIs (not really supposed to be nice, but to give a quick overview/rapid diagnostic)
	#input:
		#either a three columns data frame or matrix (with rownames) handed over to argument est (first column needs to comprise the original estimate)
		#or several vectors handed over as follows:
			#est: numeric; estimated coefficients of the model
			#lower: numeric; lower limits of the estimates (either from model stability of from bootstrap)
			#upper: numeric; upper limits of the estimates (either from model stability of from bootstrap)
			#xnames: character; names to be depicted besides the error bars
		#col: character; name of the color with which data should be depicted
	if(ncol(est)==3){
		lower=est[, 2]
		upper=est[, 3]
		xnames=rownames(est)
		est=est[, 1]
		x.at=est
	}
	if(is.null(xlim)){
		xlim=range(c(0, lower, upper))
	}
  old.par = par(no.readonly = TRUE)
	par(mar=c(3, 0.5, 0.5, 0.5), mgp=c(1, 0.4, 0), tcl=-0.2)
	plot(x=est, y=1:length(est), pch=pch, xlab="estimate", ylab="", yaxt="n", xlim=xlim, type="n", ylim=c(1, length(est)+1))
	abline(v=0, lty=3)
	if(center.at.null){
		x.at=rep(0, length(est))
		text(labels=xnames, x=x.at, y=(1:length(est))+0.3, cex=0.8)
	}else{
		xx=range(est)#c(lower, upper))
		for(i in 1:length(est)){
			text(labels=xnames[i], x=x.at[i], y=i+0.3, cex=0.8, adj=c((est[i]-min(xx))/diff(xx), 0))#pos=c(2, 4)[1+as.numeric(est<0)])
		}
	}
	points(x=est, y=1:length(est), pch=pch, col=col)
	segments(x0=lower, x1=upper, y0=1:length(est), y1=1:length(est), col=col)
	if(reset.par){par(old.par)}
}

################################################################################
#function determining the correction factor for AIC to be added in order to get AICc
aic.c.fac<-function(N, k){
  return(2*k*(k+1)/(N-k-1))
}

################################################################################
#function determining Akaike weights for a set of AIC-values and also the 95% best model confidence set
conf.set<-function(aic){
  o=order(aic)
  #aic=aic[o]
  d.aic=aic-min(aic)
  e.aic=exp(-d.aic/2)
  w.aic= e.aic/sum(e.aic)
  cum=cumsum(w.aic[o])[order(o)]
  c.set=cum<=0.95
  c.set[which.min(aic)]=1
  if (max(cum[c.set==T])<0.95){
    c.set[cum==min(cum[c.set==F])]=T
  }
  m.rank=rank(aic)
  result=data.frame(aic, d.aic, w.aic, cum, c.set, m.rank)
  #o=order(o)
  #result=result[o,]
  return(result)
}

################################################################################
#function creating all subset models for a given model
library(gtools)
all.models<-function(model){
  #browser()
  model=attr(terms(as.formula(paste(c("~", model), collapse=""))), "term.labels")
  count.char<-Vectorize(function(xtext, char){
    sum(unlist(strsplit(xtext, split=""))==char)
  })
  #browser()
  model=model[order(count.char(model, "^"), count.char(model, ":"))]
  #build models with main effects only:
  to.treat=model[!grepl(model, pattern=":", fixed=T) & !grepl(model, pattern="^", fixed=T)]
  all.models=cbind(T, permutations(n=2, r=length(to.treat), v=c(F, T), repeats.allowed=T))
  to.treat=c(1, to.treat)
  all.models=lapply(1:nrow(all.models), function(x){to.treat[all.models[x,]]})
  #include squared terms (not involved in interactions):
  to.treat=model[!grepl(model, pattern=":", fixed=T) & grepl(model, pattern="^", fixed=T)]
  term.names=to.treat
  term.names=gsub(term.names, pattern="I(", replacement="", fixed=T)
  term.names=gsub(term.names, pattern="^2)", replacement="", fixed=T)
  if(length(to.treat)>0){
    for(i in 1:length(term.names)){
      all.models=c(all.models, lapply(all.models[unlist(lapply(all.models, function(x){length(intersect(term.names[i], x))==1}))], function(y){c(y, to.treat[i])}))
    }
  }
  #include interactions:
  to.treat=model[grepl(model, pattern=":", fixed=T)]
  if(length(to.treat)>0){
    for(i in 1:length(to.treat)){
      xx=unlist(strsplit(to.treat[i], split=":", fixed=T))
      yy=permutations(n=2, r=length(xx), v=c(F, T), repeats.allowed=T)
      yy=yy[-c(1, nrow(yy)),]
      yy=apply(yy, 1, function(x){paste(xx[x], collapse=":")})
      xx=yy[grepl(yy, pattern="^", fixed=T)]
      if(length(xx)>0){
        xx=c(xx, to.treat[i])
        xx=gsub(xx, pattern="I(", replacement="", fixed=T)
        xx=gsub(xx, pattern="^2)", replacement="", fixed=T)
      }
      yy=c(yy, xx)
      yy=unlist(lapply(all.models, function(x){length(intersect(x, yy))==length(yy)}))
      all.models=c(all.models, lapply(all.models[yy], function(x){c(x, to.treat[i])}))
    }
  }
  all.models=unlist(lapply(all.models, paste, collapse="+"))
  return(all.models)
}

################################################################################
#functions averaging for each data point the residuals of all others (potentially of only those which belong to the same level
	#of a grouping variable); averages are weighted by distance
get.2d.ac<-function(resis, ac.sd, t.ac.sd, lat, long=NULL, xtime, contr.fac=NULL){
  if(length(contr.fac)==0){contr.fac=rep(1, length(resis))}
  contr.fac=as.factor(contr.fac)
  if(length(long)==0){long=rep(0, length(lat))}
  all.ac=rep(NA, length(resis))
  c.fac.lev=levels(contr.fac)
  for (i in 1:length(c.fac.lev)){
    ind.resis=resis[contr.fac==c.fac.lev[i]]
    ind.lat=lat[contr.fac==c.fac.lev[i]]
    ind.long=long[contr.fac==c.fac.lev[i]]
    ind.time=xtime[contr.fac==c.fac.lev[i]]
    xx=unlist(lapply(1:length(ind.resis), function(cell){
      ac.weights=dnorm(x=sqrt((ind.lat[cell]-ind.lat[-cell])^2+ (ind.long[cell]-ind.long[-cell])^2), mean=0, sd=ac.sd)
      t.ac.weights=dnorm(x=abs(ind.time[cell]-ind.time[-cell]),mean=0, sd=t.ac.sd)
      ac.weights[is.nan(ac.weights)]=0
      t.ac.weights[is.nan(t.ac.weights)]=0
      weighted.mean(ind.resis[-cell], w=ac.weights*t.ac.weights)
    }))
    all.ac[contr.fac==c.fac.lev[i]]=xx
  }
  all.ac[is.nan(all.ac)]=0
  all.ac[is.na(all.ac)]=0
  return(all.ac)
}

get.1d.ac.from.dist.mat<-function(resis, ac.sd, dist.mat, excl.contr.fac=NULL){
  #if(length(contr.fac)==0){contr.fac=rep(1, length(resis))}
  #contr.fac=as.factor(contr.fac)
  #c.fac.lev=levels(contr.fac)
	if(is.null(excl.contr.fac)){excl.contr.fac=1:length(resis)}
  excl.contr.fac=as.factor(excl.contr.fac)
  excl.contr.fac.levels=levels(excl.contr.fac)
  all.ac=rep(NA, length(resis))
  for(i in 1:length(excl.contr.fac.levels)){
    ind.resis=resis[excl.contr.fac!=excl.contr.fac.levels[i]]
    ind.dist.mat=dist.mat[excl.contr.fac!=excl.contr.fac.levels[i], ]
    xx=unlist(lapply(1:length(ind.resis), function(cell){
      ac.weights=dnorm(x=as.vector(ind.dist.mat[, cell]), mean=0, sd=ac.sd)
      ac.weights[is.nan(ac.weights)]=0
      #browser()
      weighted.mean(ind.resis, w=ac.weights)
    }))
    all.ac[excl.contr.fac!=excl.contr.fac.levels[i]]=xx
  }
  all.ac[is.nan(all.ac)]=0
  all.ac[is.na(all.ac)]=0
  return(all.ac)
}

##function fitting a GLM negative binomial error distribution (today one could use glmmTMB instead):
##written by Roger Mundry
fit.fun<-function(formula, data, offset=NULL, tol=0.0001){#browser
	object=as.formula(paste(c("~", unlist(strsplit(formula, split="~", fixed=T))[2]), collapse=""))
	m.mat=model.matrix(object=object, data=data)
	if(is.null(offset)){offset=rep(1, nrow(data))}
	poiss.m=glm(as.formula(paste(c(formula, "offset(log(offset))"), collapse="+")), data, family=poisson)
	ests=coefficients(poiss.m)
	ests=c(ests, log.theta=log(10))
	rv=data[, gsub(x=unlist(strsplit(x=formula, split="~", fixed=T))[1], pattern=" ", fixed=T, replacement="")]
	opt.fun<-function(ests, m.mat, rv, offset){
		mu=m.mat[, names(ests)[-length(ests)], drop=F]%*%ests[-length(ests)]+log(offset)
		return(-sum(dnbinom(x=rv, mu=exp(mu), size=exp(ests["log.theta"]), log=T)))
	}
	old.ll=as.vector(logLik(poiss.m))
	change=T
	while(change){
		opt.res=optim(par=ests, f=opt.fun, m.mat=m.mat, rv=rv, offset=offset, hessian=T, control=list(maxit=10000), method="BFGS")#, control=list(fnscale=-1)
		if(abs(opt.res$value-old.ll)<tol){change=F}
		old.ll=opt.res$value
		ests=opt.res$par
	}
	#ests=opt.res$par
	fitted.vals=as.vector(exp(m.mat[, names(ests)[-length(ests)], drop=F]%*%ests[-length(ests)]+log(offset)))
	fitted.var= fitted.vals + fitted.vals^2/exp(ests["log.theta"])
	resids=(rv-fitted.vals)#/sqrt(fitted.vals)
	disp.par=sum((resids/sqrt(fitted.var))^2)/(length(resids)-length(opt.res$par))
	coefs=data.frame(Estimate=opt.res$par, SE=sqrt(diag(solve(opt.res$hessian))))
	coefs$z=coefs$Estimate/coefs$SE
	coefs$P=2*pnorm(-abs(coefs$z))
	return(list(coefs=coefs, logLik=-opt.res$value, disp.par=disp.par, fitted.vals=fitted.vals, resids=resids, pearson.resids=resids/sqrt(fitted.var), opt.res=opt.res))
}

##an attempt to fit non-linear models (not completed)
glm.nb.nl<-function(formula, non.linear=NULL, data, offset=NULL, tol=0.0001){#browser
	object=as.formula(paste(c("~", unlist(strsplit(formula, split="~", fixed=T))[2]), collapse=""))
	m.mat=model.matrix(object=object, data=data)
	if(is.null(offset)){offset=rep(1, nrow(data))}
	poiss.m=glm(as.formula(paste(c(formula, "offset(log(offset))"), collapse="+")), data, family=poisson)
	ests=coefficients(poiss.m)
	ests=c(ests, theta=10)
	rv=data[, gsub(x=unlist(strsplit(x=formula, split="~", fixed=T))[1], pattern=" ", fixed=T, replacement="")]
	opt.fun<-function(ests, m.mat, rv, offset){
		mu=m.mat[, names(ests)[-length(ests)], drop=F]%*%ests[-length(ests)]+log(offset)
		return(-sum(dnbinom(x=rv, mu=exp(mu), size=ests["theta"], log=T)))
	}
	old.ll=as.vector(logLik(poiss.m))
	change=T
	while(change){
		opt.res=optim(par=ests, f=opt.fun, m.mat=m.mat, rv=rv, offset=offset, hessian=T, control=list(maxit=10000), method="BFGS")#, control=list(fnscale=-1)
		if(abs(opt.res$value-old.ll)<tol){change=F}
		old.ll=opt.res$value
		ests=opt.res$par
	}
	#ests=opt.res$par
	fitted.vals=as.vector(exp(m.mat[, names(ests)[-length(ests)], drop=F]%*%ests[-length(ests)]+log(offset)))
	fitted.var= fitted.vals + fitted.vals^2/ests["theta"]
	resids=(rv-fitted.vals)#/sqrt(fitted.vals)
	disp.par=sum((resids/sqrt(fitted.var))^2)/(length(resids)-length(opt.res$par))
	coefs=data.frame(Estimate=opt.res$par, SE=sqrt(diag(solve(opt.res$hessian))))
	coefs$z=coefs$Estimate/coefs$SE
	coefs$P=2*pnorm(-abs(coefs$z))
	return(list(coefs=coefs, logLik=-opt.res$value, disp.par=disp.par, fitted.vals=fitted.vals, resids=resids, pearson.resids=resids/sqrt(fitted.var), opt.res=opt.res))
}

##function plotting a quantitative response against a covariate and adding the model line
##written by Roger Mundry
no.int.plot<-function(
	plot.data, covariate, coefs=NULL, response, link=c("identity", "logit", "log"),
  grid.resol=11, 
	xlab="", ylab="", x.at=NULL, x.labels=NULL, y.at=NULL, y.labels=NULL, xlim=NULL, ylim=NULL,
	cex.lab=1, cex.axis=1, size.fac=NA, col="black", pch=1, conf.int=NULL, ci.type=c("lines", "area"), lty=2, lwd=1, pt.lwd=1,
	my.par=list(mar=c(3, 3, 0.2, 0.2), mgp=c(1.7, 0.3, 0), tcl=-0.15, las=1), quiet=T, reset.par=F, mean.fun=mean,
	weights=NULL){
	#last modified: 2018 AUG 20, works also with model without intercept
  old.par = par(no.readonly = TRUE)
  #version from april 15th 2016
	if(length(ci.type)>1){
		ci.type=ci.type[1]
	}
	#browser()
  link=link[1]
	#extract coefficients needed:
	if(!is.null(coefs)){
		coefs=coefs[names(coefs)%in%c("(Intercept)", covariate, paste(c("I(", covariate, "^2)"), collapse=""))]
		#coefs=coefs[c(which(names(coefs)=="(Intercept)"), grep(x=names(coefs), pattern=covariate))]
		#coefs=coefs[!grepl(x=names(coefs), pattern=":", fixed=T)]
	}
	if(length(dim(response))==2){
		response=response[, 1]/apply(response, 1, sum)
	}
  if(is.null(weights)){
		weights=rep(1, nrow(plot.data))
	}
	sel=!is.na(weights) & !is.na(response) & !is.na(plot.data[, covariate])
	weights=weights[sel]
	response=response[sel]
	plot.data=subset(plot.data, sel)
  if(!is.na(grid.resol)){
		xvar=seq(min(plot.data[, covariate]), max(plot.data[, covariate]), length.out=grid.resol)
		bin.x=cut(x=plot.data[, covariate[1]], breaks=xvar, labels=F, include.lowest=T)
		bin.x=min(xvar)+diff(xvar)[1]/2+((bin.x-1)*diff(xvar)[1])
  }else{
    bin.x=plot.data[, covariate]
  }

  if(!is.na(size.fac) & !is.na(grid.resol)){
		obs.mean=aggregate(response, list(bin.x), mean.fun)
    #get sample size per 'cell':
    N=aggregate(weights, by=list(bin.x), sum)
    #... and turn the result into a data frame:
		N=N$x
  }else if(!is.na(size.fac) & is.na(grid.resol)){
		obs.mean=aggregate(response, list(bin.x, response), mean.fun)
    N=aggregate(weights, by=list(bin.x, response), sum)
		N=N$x
  }else{
		obs.mean=aggregate(response, list(bin.x), mean.fun)
    N=rep(1, nrow(obs.mean))
    size.fac=1
  }
  names(obs.mean)[1]=covariate
	names(obs.mean)[2]="Freq"
  obs.mean=data.frame(obs.mean, N)
  obs.mean=subset(obs.mean, N>0)
	#get predicted values:
	xvals=seq(min(plot.data[, covariate]), max(plot.data[, covariate]), length.out=100)
  if(any(names(coefs)=="(Intercept)")){
		x.pred.data=data.frame(1, xvals)
		names(x.pred.data)=c("(Intercept)", covariate)
	}else{
		x.pred.data=data.frame(xvals)
		names(x.pred.data)=covariate
	}
	
	if(!is.null(coefs)){
		model.terms=names(coefs)
		if(length(intersect(model.terms, "(Intercept)"))>0){model.terms[model.terms=="(Intercept)"]="1"}
		rr=runif(n=nrow(x.pred.data))
		pvs=model.matrix(object=as.formula(paste(c("rr", paste(model.terms, collapse="+")), collapse="~")), data=x.pred.data)
		yvals=apply(t(coefs*t(pvs[, names(coefs)])), 1, sum)
		#transform according to link function
		if(link=="log"){yvals=exp(yvals)}
		if(link=="logit"){yvals=exp(yvals)/(1+exp(yvals))}
	}

	#prepare axis ticks and labels
  if(length(xlim)==0){
    xlim=range(plot.data[,covariate], na.rm=T)
  }
  if(length(ylim)==0){
		if(!is.null(coefs)){
			ylim=range(c(response, yvals))
		}else{
			ylim=range(response)
		}
	}
  if(length(x.at)==0){
    x.at=pretty(plot.data[, covariate])
    #x.at=x.at[x.at>=min(plot.data[, covariate])]
    #x.at=x.at[x.at<=max(plot.data[, covariate])]
    x.labels=x.at
  }
  if(length(y.at)==0){
    y.at=pretty(ylim)
    #y.at=y.at[y.at>=min(ylim)]
    #y.at=y.at[y.at<=max(ylim)]
    y.labels=as.character(y.at)
  }

  #plot:
  par(my.par)
	plot(obs.mean[, covariate], obs.mean$Freq, cex=size.fac*sqrt(obs.mean$N), type="n", 
		xlim=xlim, ylim=ylim, xaxt="n", yaxt="n", xlab=xlab, ylab=ylab, tcl=-0.25, cex.lab=cex.lab)
	#add axes:
	axis(at=x.at, labels=x.labels, tcl=-0.2, side=1, tcl=-0.25, cex.axis=cex.axis)
	#mtext(text=xaxt.labels, at=xaxt.at, side=1, line=0.3, cex=0.7)
	axis(at=y.at, tcl=-0.2, side=2, las=1, labels=y.labels, tcl=-0.25, cex.axis=cex.axis)
	#add CI:
	if(length(conf.int)>0){
		if(ci.type=="lines"){
			ci.yvals=conf.int[, "lwr"]
			lines(x=conf.int[, covariate], y=conf.int[, "lwr"], lty=3, lwd=lwd)
			lines(x=conf.int[, covariate], y=conf.int[, "upr"], lty=3, lwd=lwd)
		}else if(ci.type=="area"){
			polygon(x=c(conf.int[, covariate], rev(conf.int[, covariate])), 
				y=c(conf.int[, "lwr"], rev(conf.int[, "upr"])), border=NA, col="grey")
		}
	}
	points(obs.mean[, covariate], obs.mean$Freq, cex=size.fac*sqrt(obs.mean$N), pch=pch, col=col, lwd=pt.lwd)
	if(!is.null(coefs)){
		lines(xvals, yvals, lty=lty, lwd=lwd)
	}
	if(reset.par){par(old.par)}
	if(!quiet){return(obs.mean)}
}
