library(lme4)
library(plyr)
library(ggplot2)
library(tidyr)
library(gridExtra)
set.seed(2016)

################################################################################################
########### simulation to create continuous covariate 
################################################################################################
#set number of individuals and number of GPS locs per individual as desired 
indv <- 50
locs <- 1000
#set min and max of mean available across individuals
min <- 0.25 
max <- 0.6  
#variance in the simulated data; lower to make data better behaved, raise to simulate messy data
std <-.2

#make sequence of mean values for different scenarios using min and max above;
#these values will determine the mean of the distribution used to simulate
#the data (so the 'constant' scenario draws random numbers from a normal
#distribution with a mean of 0.5 and a standard deviation
# of 0.2).  To change any of the sequences, you can change the min and max above,
# or the constants  in the code below
avail <- seq(min,max,length=indv)
prop.use <- seq(min,max,length=indv)
constant <- rep(0.5,indv)
additive <- seq(min,max,length=indv)+0.2
specialist <- seq(min,max,length=indv)+seq(0.3,0,length=indv)
increasing <- seq(min,max,length=indv)+seq(0.10,0.25,length=indv)
trade <- seq(min,max,length=indv)+seq(0.1,-0.1,length=indv)
decreasing <- seq(min,max,length=indv)+ seq(-0.09,-0.18,length=indv)
relaxed <- seq(min,max,length=indv)+seq(-0.5,0,length=indv)
use.list <- list(prop.use,constant,additive,specialist,increasing,trade,decreasing,relaxed)

#create empty data set for simulated values, size will change depending on parameters above
data <- matrix(nrow=indv*locs,ncol=(length(use.list)+2))
#loop that populates values; change sd parameter for more or less variation within individuals
for (i in 1:indv) {
	r <- 1:locs+(locs*i)-locs
	data[r,1] <- rep(i,locs)
	data[r,2] <- rnorm(n=locs,sd=std,mean=avail[i])
	for (j in 1:length(use.list)){
		data[r,j+2] <- rnorm(n=locs,sd=std,mean=use.list[[j]][i])
}
}

#convert matrix into dataframe
output <- data.frame(data)
names(output) <- c("id","avail","prop.use","constant", "additive", "specialist","increasing","trade","decreasing","relaxed")

#change any values above or below 1 and 0 to 1 or 0
output[,2:10][output[,2:10]>1]  <- 1
output[,2:10][output[,2:10]<0]  <- 0

#convert data from wide to long for use in regressions
all.data <- output%>%gather(type,covariate,c(avail,prop.use,constant,additive,specialist,increasing,trade,decreasing,relaxed))
all.data$type <- as.factor(all.data$type)
#add y variable indicating use (1) or availability (0)
all.data$y[all.data$type=="avail"] <- 0
all.data$y[all.data$type!="avail"] <- 1


#####################################################################################
###USE SIMULATED DATA ABOVE TO RUN IN 4 ANALYSIS APPROACHES#################
#####################################################################################

# show different scenario names in dataset
types <- levels(all.data$type)

# select which scenario to analyze; change name in quotes to any of the scenarios given in 'types'
i <- "decreasing"

#create reduced dataset of only available locations and used locations of the type selected above
data <- all.data[all.data$type=="avail"|all.data$type==i,]
data <- droplevels(data)
summary(data)

################################# use/avail approach #################################
#create summary of mean use and available for each individual
used <- ddply(data[data$type!="avail",], .(id), summarize, used.mean = mean(covariate)) 
avail <- ddply(data[data$type=="avail",], .(id), summarize, avail = mean(covariate))
#make into new dataframe
dat <- as.data.frame(cbind(used,avail.mean=avail[,2])) 
#run linear model of mean use to mean available
mod1 <- lm(used.mean ~ avail.mean, data=dat)
#predict values from model, with confidence intervals
PredUsed <- data.frame(predict(mod1, interval="confidence", level = 0.95, type="response"))
#model output to test presence of functional response
coef(mod1)
confint(mod1)
summary(mod1)$r.squared

############################### log use/log avail approach ##############################
#run linear model using log of mean use and available, dataset from above
mod2 <- lm(log(used.mean) ~ log(avail.mean), data=dat)
#predict values from model, with confidence interval
PredUsedlog <- data.frame(predict(mod2, interval="confidence", level = 0.95, type="response"))
#make log values of avail and use, for graphing
log.avail <- log(dat$avail.mean)
log.use <- log(dat$used.mean)
#model output to test presence of functional response
coef(mod2)
confint(mod2)
summary(mod2)$r.squared

################################## random slope approach ###############################
#run general linear mixed model with individual slopes and intercept
mod3 <- glmer(as.factor(y)~ covariate + (1|id)+(0+covariate|id), family="binomial", data=data)
#get modeled estimate of slope for each individual
dat$slopes <- coef(mod3)$id[,"covariate"]
#run linear model of individual slopes as a function of log availability
mod3a <- lm(slopes ~ log(avail.mean), data=dat)
#predict values from linear model, with confidence intervals
PredSlope <- data.frame(predict(mod3a, interval="confidence", level = 0.95, type="response"))
#get model summary
summary(mod3a)
confint(mod3a, level=0.95)
summary(mod3a)$r.squared

########################### habitat selection with interactions ######################
#create new dataset, assigning each individual's mean availability to all its observations
dat.merge <- merge(data, avail, by="id")
#run general linear mixed model predicting use or available (y) given the interaction between the covariate and mean availability
glmer1.1 <- glmer(as.factor(y)~ covariate + covariate:avail + (1|id), family="binomial", data=dat.merge)
#second model run for comparison with no interaction
glmer1.2 = glmer(as.factor(y)~ covariate + (1|id), family="binomial", data=dat.merge)
#difference in AIC and BIC between the two models, to test for a functional response
(diff.aic <- AIC(glmer1.2)-AIC(glmer1.1))
(diff.bic <- BIC(glmer1.2)-BIC(glmer1.1))
summary(glmer1.1)
confint(glmer1.1,parm="beta_",method="Wald", level=0.95)
#used in graph to indicate significance
if(diff.aic>=2){x="*"}else{x=""}

#set range of availabilities for prediction
av = data.frame(seq(min,max,by = 0.01))

#predict probability of selection using model with interactions across a range of availabilties
rsf1 = exp(glmer1.1@beta[2]*mean(dat.merge$covariate) + glmer1.1@beta[3]*(mean(dat.merge$covariate)*av)) #range = 0.05-0.80
interact1 = cbind(av,rsf1)
names(interact1) <- c("Avail","RSF")

#####################################create graphs for each approach##########################################

#compile all results into single dataframe for graphing

results <-as.data.frame(cbind(dat,PredUsed,PredSlope,PredUsedlog,log.avail,log.use))
names(results)<-c("ID","used.mean","avail.mean","slopes","pred.use","use.lwr","use.upr",
	"pred.slope","slope.lwr","slope.upr","pred.log","log.lwr","log.upr","log.avail","log.use")

#function to set decimal places to make graphs better looking
fmt.dec <- function(decimals=0){
	function(x) as.character(round(x,decimals))
}

#create graphs
used <- ggplot(results, aes(x=avail.mean, y=used.mean)) +
	geom_point(shape=19, cex=2) +
	geom_line(data=results, aes(x=avail.mean,y=pred.use), lwd = 1) +
	geom_ribbon(data=results, aes(ymin=use.lwr, ymax=use.upr), alpha=0.3) +
	geom_abline(aes(intercept=0,slope=1), colour = "black", linetype="dashed", lwd = 1)+
	xlab("Available") + ylab("Use") + theme_bw(base_size = 10) +
	scale_x_continuous(limits=c(min(results$use.lwr,results$avail.mean),max(results$use.upr,results$avail.mean)),label=fmt.dec(1))+
	scale_y_continuous(limits=c(min(results$use.lwr,results$avail.mean,results$used.mean),max(results$use.upr,results$avail.mean)),label=fmt.dec(1))+
	theme(axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
	axis.line.y = element_line(colour = "black", size=0.5, linetype='solid'),
	axis.text.x = element_text(colour="black", size=8), axis.text.y = element_text(colour="black", size=8)) +
	theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())

log <- ggplot(results,aes(x=log.avail,y=log.use))+
	geom_point(shape=19, cex=2) +
	geom_line(data=results, aes(x=log.avail,y=pred.log), lwd = 1) +
	geom_ribbon(data=results, aes(ymin=log.lwr, ymax=log.upr), alpha=0.3) +
	geom_abline(aes(intercept=0,slope=1), colour = "black", linetype="dashed", lwd = 1)+
	xlab("Log(Available)") + ylab("Log(Use)") + theme_bw(base_size = 10) +
	scale_x_continuous(limits=c(min(results$log.avail,results$log.use,results$log.lwr),max(results$log.avail,results$log.use,results$log.upr)), label=fmt.dec(1))+
	scale_y_continuous(limits=c(min(results$log.avail,results$log.use,results$log.lwr),max(results$log.avail,results$log.use,results$log.upr)))+
	theme(axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
	axis.line.y = element_line(colour = "black", size=0.5, linetype='solid'),
	axis.text.x = element_text(colour="black", size=8), axis.text.y = element_text(colour="black", size=8)) +
	theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())

slopes <- ggplot(results, aes(x=avail.mean, y=slopes))  +
	geom_point(shape=19, cex=2) +
	geom_line(data=results, aes(x=avail.mean,y=pred.slope), lwd = 1) +
	geom_ribbon(data=results, aes(ymin=slope.lwr, ymax=slope.upr), alpha=0.3) +
	geom_hline(aes(yintercept=0), colour = "black", linetype="dashed", lwd = 1)+
	#scale_x_continuous(labels=fmt_dcimals(1))+
	scale_y_continuous(labels=fmt.dec(1),limits=c(-6,0))+
	xlab("Available") + ylab("Selection Coefficient") + theme_bw(base_size = 10) +
	theme(axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
	axis.line.y = element_line(colour = "black", size=0.5, linetype='solid'),
	axis.text.x = element_text(colour="black", size=8), axis.text.y = element_text(colour="black", size=8)) +
	theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())

interact <- ggplot(interact1, aes(x=Avail, y=RSF)) + 
	geom_line(data=interact1, aes(x=Avail,y=RSF), lwd = 1) +
	geom_hline(aes(yintercept=1),color="black",linetype="dashed",lwd=1)+
	xlab("Available") + ylab("w(x)") + theme_bw(base_size = 10) +
	annotate("text",x=0.55,y=.95,vjust=1,label=x,size=8)+ #max(interact1$Avail) max(interact1$RSF)
	#scale_x_continuous(limits=c(interact1$)+
	#scale_y_continuous(labels=fmt.dec(1))+
	theme(axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
	axis.line.y = element_line(colour = "black", size=0.5, linetype='solid'),
	axis.text.x = element_text(colour="black", size=8), axis.text.y = element_text(colour="black", size=8)) +
	theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())


#compile graphs into a single panel and save

setwd("XXX")

tiff(filename="XXX.tif",width=190,height=40,units="mm",res=600)
grid.arrange(used,log,slopes,interact,nrow=1)
dev.off()