################################
###Within this file are scripts for 
#Tosa, MI, Biel M, Graves TA. 2023. Bighorn sheep interactions: understanding tradeoffs of sociality and implications for disease transmission.  Peer J
###There are 7 scripts:
#1_Where are the divisions among subpopulations?
#2_How do dyad characteristics influence contact rate
#3_When do direct contacts occur? Visualization
#4a_When do direct contacts occur? Modelling
#4b_Where are bighorn? General habitat use
#4c_Where do direct contacts occur?
#4d_ General habitat versus direct contacts

##Data are available here:
#Graves, TA, MI Tosa, EP Flesch. K.Keating. 2023.  Glacier Waterton International Peace Park bighorn sheep (Ovis canadensis) data 2002-2012. U.S Geological Survey data release. doi:10.5066/P9ANTBOF
###Predicted SSF, contactRSF, and SSF*contactRSF are available here:
#Tosa, MI, Graves, TA. 2023.  Glacier Waterton International Peace Park bighorn sheep (Ovis canadensis) data 2002-2012. U.S Geological Survey data release. doi:10.5066/P9H078A1



###############################
#1. Where are the divisions among subpopulations?
###############################

#load packages
require(plyr)
require(data.table)
require(igraph)

require(ggplot2)
require(ggpubr)
require(ggrepel)

require(sf)
require(rgdal)
require(raster)


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

#set working directory
path <- "C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper"
setwd(path)

###############################
#load data

gps.data <- read.table("AllGPSDataCoordOnly.txt", header=T)

gps.data$ANIMALID <- as.character(gps.data$ANIMALID.)
gps.data$ANIMALID. <- gsub(gps.data$ANIMALID., pattern="-", replacement="")

gps.data$dt <- as.POSIXct(paste(gps.data$MSDATE., gps.data$MSTIME.), tz="UTC", format="%Y.%m.%d %H:%M:%S")
table(gps.data$ANIMALID.)

sheep <- unique(gps.data[,c("ANIMALID.","SEX.")])
table(sheep$SEX.)

#clean up gps.data since some duplicates in file, started with 168,422 rows, 168,380 rows non-duplicated time stamps
gps.data <- gps.data[!duplicated(gps.data[,c("dt","ANIMALID.")]),]

#########
#calculate centroid of each sheep's home range
centroid <- ddply(gps.data, .(ANIMALID., COLLARID., GROUP., SEX., AGE), 
                  summarize, mean_x=mean(UTME_NAD83, na.rm=T), min_x=min(UTME_NAD83, na.rm=T), max_x=max(UTME_NAD83, na.rm=T),
                  mean_y=mean(UTMN_NAD83, na.rm=T), min_y=min(UTMN_NAD83, na.rm=T), max_y=max(UTMN_NAD83, na.rm=T))

centroid.p <- ggplot(centroid) +
  geom_point(aes(x=mean_x, y=mean_y, col=SEX.)) +
  geom_text_repel(aes(x=mean_x, y=mean_y, col=SEX., label=ANIMALID.)) +
  theme_bw(base_size=20)
centroid.p
#ggsave(centroid.p, filename = "Figures/centroid_labeled_2022-03-07.tiff", height=8, width=8, units="in", dpi=300, compression="lzw")

#########
#calculate distance between each bighorn location

gps.dt <- data.table(gps.data[order(gps.data$dt),])
DT <- gps.dt
coords=c('UTME_NAD83','UTMN_NAD83')

#
indirect <- data.frame()
for(sheep1 in unique(DT$ANIMALID.))
{
  print(sheep1)
  if(sheep1 == "0420")
    next
  # sheep1 <- "0420"
  DT.sheep1 <- DT[DT$ANIMALID. == sheep1,]
  
  for(sheep2 in unique(DT$ANIMALID.))
  {
    if(sheep1 == sheep2)
      next
    if(sheep2 == "0420")
      next
    if(nrow(indirect[indirect$sheep1 == sheep1 & indirect$sheep2 == sheep2,]) > 0) #skip if already calculated the number of indirect contacts for the same dyad
      next
    if(nrow(indirect[indirect$sheep1 == sheep2 & indirect$sheep2 == sheep1,]) > 0) #skip if already calculated the number of indirect contacts for the reciprocal dyad
      next
    
    print(sheep2)
    DT.sheep2 <- DT[DT$ANIMALID. == sheep2,]
    
    #reset distance counters
    n10 <- 0
    n25 <- 0
    n50 <- 0
    n100 <- 0
    for(p1 in 1:nrow(DT.sheep1))
    {
      for(p2 in 1:nrow(DT.sheep2))
      {
        d <- sqrt((DT.sheep1$UTME_NAD83[p1]-DT.sheep2$UTME_NAD83[p2])^2 + (DT.sheep1$UTMN_NAD83[p1]-DT.sheep2$UTMN_NAD83[p2])^2)
        if(d < 10)
          n10 <<- n10 + 1
        if(d < 25)
          n25 <<- n25 + 1
        if(d < 50)
          n50 <<- n50 + 1
        if(d < 100)
          n100 <<- n100 + 1
      }
    }
    indirect <- rbind(indirect, data.frame(sheep1 = sheep1, sheep2 = sheep2, indirect_10 = n10, indirect_25 = n25, indirect_50 = n50, indirect_100 = n100))
    # write.table(indirect, file="indirect_contact_counts_2022-03-07.txt", sep=",", row.names=F)
    rm(DT.sheep2)
    gc()
  }
}

done <- table(indirect[,c("sheep1","sheep2")])
done

indirect <- indirect[!duplicated(indirect),]
# write.table(indirect, file = "indirect_contact_counts_2022-03-07.txt", sep=",", row.names = F)

#########
indirect <- read.table("indirect_contact_counts_2022-03-07.txt", sep=",", header=T, colClasses=c("character","character","numeric","numeric","numeric","numeric")) #includes distance = 10m, 25m, 50m, 100m

#merge with sex
indirect <- merge(indirect, sheep, by.x="sheep1", by.y="ANIMALID.", all.x=T)
indirect <- merge(indirect, sheep, by.x="sheep2", by.y="ANIMALID.", all.x=T, suffixes=c("1","2"))
indirect$d.sex <- paste(indirect$SEX.1, indirect$SEX.2, sep="")

sum(indirect$indirect_10, na.rm=T) #138,469 indirect contacts
sum(indirect[!is.na(indirect$indirect_10),]$indirect_25, na.rm=T) #can't just divide by 2 because some dyads not in table in reciprocal form

sum(indirect$indirect_25, na.rm=T) #1,259,882 indirect contacts w/ 25m

indirect %>% group_by(d.sex) %>% dplyr::summarize(sum10=sum(indirect_10, na.rm=T), sum25=sum(indirect_25, na.rm=T))
#   d.sex sum10  sum25
# 1 FF    47,470 443,426
# 2 FM    31,333 222,043
# 3 MF    18,723 221,913
# 4 MM    40,943 372,500

#    31333 + 18723 = 50,056 MF 10m
# 222043 + 221913 = 443,956 MF 25m

#make sure no duplicate dyads
indirect[duplicated(indirect),]

edges <- indirect
names(edges) <- c("from","to","indirect_25","indirect_50","indirect_100","indirect_10","sex.i","sex.j","d.sex")
edges <- edges[!is.na(edges$indirect_10),]

#create networks
net10 <- graph_from_data_frame(edges[edges$indirect_10 > 0, c("from","to","indirect_10")], vertices=centroid, directed=F)
plot(net10)

net25 <- graph_from_data_frame(edges[edges$indirect_25 > 0, c("from","to","indirect_25")], vertices=centroid, directed=F)
plot(net25)

net50 <- graph_from_data_frame(edges[edges$indirect_50 > 0, c("from","to","indirect_50")], vertices=centroid, directed=F)
plot(net50)

net100 <- graph_from_data_frame(edges[edges$indirect_100 > 0, c("from","to","indirect_100")], vertices=centroid, directed=F)
plot(net100)

#find groups that are connected
gr10 <- groups(components(net10))
gr10

gr25 <- groups(components(net25))
gr25

gr50 <- groups(components(net50))
gr50

gr100 <- groups(components(net100))
gr100

#########
#only plot connected nodes
v.colors <- c(F = "pink", M="steelblue") #color nodes

#indirect contacts within 10 m
#group 1: Many Glacier
net10_1 <- graph_from_data_frame(edges[edges$indirect_10 > 0 & edges$from %in% gr10$`1`, c("from","to","indirect_10")], 
                                 vertices = centroid[centroid$ANIMALID. %in% gr10$`1`,], directed=F)
V(net10_1)$color <- v.colors[V(net10_1)$SEX.]
V(net10_1)$label.color <- c("black")
plot(net10_1, 
     layout=as.matrix(centroid[centroid$ANIMALID. %in% gr10$`1`,c("mean_x","mean_y")]), 
     edge.width=E(net10_1)$indirect_10/100,
     vertex.label.color="black",
     vertex.frame.color=v.colors,
     vertex.color=v.colors,
     vertex.shape="square")

#group 2:
net10_2 <- graph_from_data_frame(edges[edges$indirect_10 > 0 & edges$from %in% gr10$`2`, c("from","to","indirect_10")], 
                                   vertices = centroid[centroid$ANIMALID. %in% gr10$`2`,], directed=F)
V(net10_2)$color <- v.colors[V(net10_2)$SEX.]
V(net10_2)$label.color <- c("black")
plot(net10_2, 
     layout=as.matrix(centroid[centroid$ANIMALID. %in% gr10$`2`,c("mean_x","mean_y")]), 
     edge.width=E(net10_2)$indirect_10/100,
     vertex.label.color="black",
     vertex.frame.color=v.colors,
     vertex.color=v.colors,
     vertex.shape="square")


#indirect contacts within 25 m
#group 1: Many Glacier
net25_1 <- graph_from_data_frame(edges[edges$indirect_25 > 0 & edges$from %in% gr25$`1`, c("from","to","indirect_25")], 
                                 vertices = centroid[centroid$ANIMALID. %in% gr25$`1`,], directed=F)
V(net25_1)$color <- v.colors[V(net25_1)$SEX.]
V(net25_1)$label.color <- c("black")
plot(net25_1, 
     layout=as.matrix(centroid[centroid$ANIMALID. %in% gr25$`1`,c("mean_x","mean_y")]), 
     edge.width=E(net25_1)$indirect_25/100,
     vertex.label.color="black",
     vertex.frame.color=v.colors,
     vertex.color=v.colors,
     vertex.shape="square")

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

#plot in ggplot
ggedges <- merge(edges, centroid[,c("ANIMALID.","mean_x","mean_y")], by.x="from", by.y="ANIMALID.", all.x=T)
ggedges <- merge(ggedges, centroid[,c("ANIMALID.","mean_x","mean_y")], by.x="to", by.y="ANIMALID.", all.x=T, suffixes=c(".from",".to"))

#min and max x and y coordinates for each section
sheepextents <- data.frame(area=c("waterton","many glacier","two medicine","bellyr","westwaterton"),
                           xmin=c(279000,297500,317500,292500,279000),
                           xmax=c(297500,311500,328000,311500,286000),
                           ymin=c(5427500,5395000,5356000,5395000,5430000),
                           ymax=c(5444000,5425000,5390000,5435000,5444000))

centroid$area <- "manyglacier"
centroid[centroid$mean_x < sheepextents[sheepextents$area == "waterton",]$xmax & centroid$mean_x > sheepextents[sheepextents$area == "waterton",]$xmin & 
           centroid$mean_y < sheepextents[sheepextents$area == "waterton",]$ymax & centroid$mean_y > sheepextents[sheepextents$area == "waterton",]$ymin,]$area <- "waterton"
centroid[centroid$mean_x < sheepextents[sheepextents$area == "two medicine",]$xmax & centroid$mean_x > sheepextents[sheepextents$area == "two medicine",]$xmin & 
           centroid$mean_y < sheepextents[sheepextents$area == "two medicine",]$ymax & centroid$mean_y > sheepextents[sheepextents$area == "two medicine",]$ymin,]$area <- "twomedicine"
centroid[centroid$mean_x < sheepextents[sheepextents$area == "bellyr",]$xmax & centroid$mean_x > sheepextents[sheepextents$area == "bellyr",]$xmin & 
           centroid$mean_y < sheepextents[sheepextents$area == "bellyr",]$ymax & centroid$mean_y > sheepextents[sheepextents$area == "bellyr",]$ymin,]$area <- "bellyr"

centroid.sp <- SpatialPointsDataFrame(coords=centroid[,c("mean_x","mean_y")], data=centroid, proj4string = CRS("+init=epsg:6341")) #UTME_NAD83 zone 12N
#centroid.sp <- spTransform(centroid.sp, CRSobj = CRS("+init=epsg:6340"))
centroid.sf <- st_as_sf(centroid.sp)

#all sheep
i10.plot <- ggplot() +
  #geom_segment(data=ggedges, aes(x=mean_x.from, y=mean_y.from, xend=mean_x.to, yend=mean_y.to, size=indirect_10), col="grey50") +
  geom_segment(data=ggedges, aes(x=mean_x.from, y=mean_y.from, xend=mean_x.to, yend=mean_y.to, size=indirect_25), col="grey50", alpha=0.5) +
  scale_size_continuous(trans="log10", range=c(0.4,2.5), name="number\nindirect\ncontacts") +
  #geom_point(data=centroid, aes(x=mean_x, y=mean_y, col=SEX.), size=5, alpha=0.75) +
  geom_sf(data=centroid.sf, aes(col=SEX.), size=5, alpha=0.75) +
  scale_color_manual(values=c("#F87462","steelblue"), name="Sex") +
  xlab("") + ylab("") +
  theme_bw(base_size=12) +
  theme(panel.grid = element_blank(), axis.ticks = element_blank())
# i10.plot

#waterton
w.plot <- i10.plot +
  #  coord_cartesian(xlim=c(279000, 300000), ylim=c(5425000, 5444000)) +
  # coord_cartesian(xlim=as.numeric(sheepextents[sheepextents$area == "westwaterton",c("xmin","xmax")]),
  #                 ylim=as.numeric(sheepextents[sheepextents$area == "westwaterton",c("ymin","ymax")])) +
  coord_sf(xlim=as.numeric(sheepextents[sheepextents$area == "westwaterton",c("xmin","xmax")]),
           ylim=as.numeric(sheepextents[sheepextents$area == "westwaterton",c("ymin","ymax")])) #+
# geom_text_repel(data=centroid[centroid$area == "waterton",],
#                 aes(x=mean_x, y=mean_y, col=SEX., label=ANIMALID.), size=5, point.padding = 1) +
# w.plot

#many glacier
#no waterton sheep
# mg.plot <- i10.plot +
#   coord_sf(xlim=as.numeric(sheepextents[sheepextents$area == "many glacier",c("xmin","xmax")]),
#            ylim=as.numeric(sheepextents[sheepextents$area == "many glacier",c("ymin","ymax")]))
# geom_text_repel(data=centroid[centroid$area == "manyglacier",], 
#                 aes(x=mean_x, y=mean_y, col=SEX., label=ANIMALID.), size=5)+

#includes east waterton
mg.plot <- i10.plot +
  coord_sf(xlim=as.numeric(sheepextents[sheepextents$area == "bellyr",c("xmin","xmax")]),
           ylim=as.numeric(sheepextents[sheepextents$area == "bellyr",c("ymin","ymax")]))
# geom_text_repel(data=centroid[centroid$area == "bellyr",], 
#                 aes(x=mean_x, y=mean_y, col=SEX., label=ANIMALID.), size=5)+
# mg.plot

#two medicine
tm.plot <- i10.plot +
  coord_sf(xlim=as.numeric(sheepextents[sheepextents$area == "two medicine",c("xmin","xmax")]),
           ylim=as.numeric(sheepextents[sheepextents$area == "two medicine",c("ymin","ymax")])) +
  theme(legend.position = c(0.8, 0.85))
# geom_text_repel(data=centroid[centroid$area == "twomedicine",], 
#                 aes(x=mean_x, y=mean_y, col=SEX., label=ANIMALID.), size=5)+
# tm.plot

all <- ggarrange(w.plot + rremove("axis.text") + rremove("legend"), mg.plot + rremove("axis.text") + rremove("legend"), tm.plot + rremove("axis.text"), nrow=1, align="hv")
all

#####################################################################################################################################################
##########################
#2 How do dyad characteristics influence contact rate?
##########################

#load packages
require(spatsoc)
require(igraph)
require(plyr)
require(data.table)
require(asnipe)

require(ggplot2)
require(ggpubr)

#set working directory
path <- "C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper"
setwd(path)

###############################
#load data

gps.data <- read.table("AllGPSDataCoordOnly.txt", header=T)

gps.data$ANIMALID. <- as.character(gps.data$ANIMALID.)
gps.data$ANIMALID. <- gsub(gps.data$ANIMALID., pattern="-", replacement="")

gps.data$dt <- as.POSIXct(paste(gps.data$MSDATE., gps.data$MSTIME.), tz="UTC", format="%Y.%m.%d %H:%M:%S")
table(gps.data$ANIMALID.)

sheep <- unique(gps.data[,c("ANIMALID.","SEX.")])
table(sheep$SEX.)
#  F  M 
# 43 54 

#clean up gps.data since some duplicates in file, started with 168,422 rows, 168,380 rows non-duplicated time stamps
gps.data <- gps.data[!duplicated(gps.data[,c("dt","ANIMALID.")]),]

#########
#calculate min and max dates for sheep, 97 bighorn
summary.gps <- ddply(gps.data, .(ANIMALID., COLLARID., GROUP., SEX., AGE), summarize, min_dt = min(dt), max_dt = max(dt))

summary.plot <- ggplot(summary.gps) + 
  geom_errorbar(aes(x=ANIMALID., ymin=min_dt, ymax=max_dt, col=SEX.), lwd=1) +
  scale_color_manual(values=c("pink","steelblue")) + 
  coord_flip() +
  ylab("Sheep ID") +
  theme_bw(base_size = 15)
summary.plot
#ggsave(summary.plot, filename = "Figures/GPS_duration_summary.tiff", height=10, width=8, units="in", dpi=300, compression="lzw")

length(unique(summary.gps[summary.gps$GROUP. %in% c("WW","RC","MD","MB","EW","SS","VP"),]$ANIMALID.)) #13, waterton parks
length(unique(summary.gps[summary.gps$GROUP. %in% c("MG","CM","SM","IP","A","GM"),]$ANIMALID.)) #59, north glacier
length(unique(summary.gps[summary.gps$GROUP. %in% c("STM","NTM","SP","DF"),]$ANIMALID.)) #25, two medicine area

#########
#calculate centroid of each sheep's home range
centroid <- ddply(gps.data, .(ANIMALID., COLLARID., GROUP., SEX., AGE), 
                  summarize, mean_x=mean(UTME_NAD83, na.rm=T), min_x=min(UTME_NAD83, na.rm=T), max_x=max(UTME_NAD83, na.rm=T),
                  mean_y=mean(UTMN_NAD83, na.rm=T), min_y=min(UTMN_NAD83, na.rm=T), max_y=max(UTMN_NAD83, na.rm=T))
#########
#groupings

#time
#to use spatsoc functions, need to convert to data.table
threshold.time <- '15 minutes'
# threshold.time <- '1 hour'

gps.dt <- data.table(gps.data[order(gps.data$dt),])
group_times(gps.dt, datetime='dt', threshold=threshold.time)

###############
#group points based on geographic distance, uses chain rule to determine which animals are in the same group

#pick one of these
threshold = 25
#threshold = 50
#threshold = 100

gps.dt <- gps.dt[!duplicated(gps.dt[,c("ANIMALID.","timegroup")]),] #cut down to 168,315 rows, remove rows if timegroup is duplicated for the same sheep
group_pts(gps.dt, threshold=threshold, id='ANIMALID.', coords=c('UTME_NAD83','UTMN_NAD83'), timegroup='timegroup')

pts.in.group <- ddply(gps.dt, .(group), nrow)
pts.in.group <- ddply(gps.dt, .(group, SEX.), nrow)

direct <- pts.in.group[pts.in.group$V1 >= 2,] #find groups with >1 animal in the same time group within threshold distance
direct.locs <- gps.dt[gps.dt$group %in% direct$group,]
nrow(direct.locs)

write.table(direct.locs, file="direct_contact_locs_2022-03-11.txt", sep=",", row.names=F) #4704
#####################################################
direct.locs <- read.table(file="direct_contact_locs_2022-03-11.txt", sep=",", header=T, colClasses = c("character","numeric","character","character","numeric","numeric","character","character",
                                                                                                       "numeric","numeric","numeric","character","numeric","numeric","numeric","numeric","character","numeric","numeric","numeric"))

sheep <- unique(direct.locs[,c("ANIMALID.","SEX.","GROUP.")])

#turn direct.locs into dyad values
summary.locs <- data.frame()
for(a in unique(direct.locs$ANIMALID.))
{
  a.locs <- direct.locs[direct.locs$ANIMALID. == a,]
  for(b in unique(direct.locs$ANIMALID.))
  {
    if(a == b)
      next
    b.locs <- direct.locs[direct.locs$ANIMALID. == b,]
    summary.locs <- rbind(summary.locs, data.frame(i=a, j=b, d25=nrow(a.locs[a.locs$group %in% b.locs$group,])))
  }
}

# write.table(summary.locs, file="DyadStructure_spatsoc_2022-03-11.txt", sep=",", row.names=F)
summary.locs <- read.table(file="DyadStructure_spatsoc_2022-03-11.txt", sep=",", header=T, colClasses = c("character","character","numeric"))
summary.locs <- merge(summary.locs, sheep[,c("ANIMALID.","SEX.")], by.x="i", by.y="ANIMALID.", all.x=T)
summary.locs <- merge(summary.locs, sheep[,c("ANIMALID.","SEX.")], by.x="j", by.y="ANIMALID.", all.x=T, suffixes=c("i","j"))
summary.locs$d.sex <- paste(summary.locs$SEX.i, summary.locs$SEX.j, sep="")

sum(summary.locs$d25) #5324 direct contacts
summary.locs %>%
  group_by(d.sex) %>%
  dplyr::summarize(sum(d25))
# 1 FF          2098
# 2 FM           255
# 3 MF           255
# 4 MM          2716

summary.locs <- summary.locs[summary.locs$d25 > 0,] #remove dyads with 0 d25 contacts
nrow(summary.locs) #252 dyads > 0 d25 contacts

###################
#censor within group dyads
#remove these dyads from analysis because correlated movement with other bighorn
#during winter: 03-16.03-19, 02-06.02-07, 05-07.05-15, 02-04.02-05, 03-07.03-20
#during lambing: 02-01.03-02, 09-01.09-02, 11-04.11-05, 02-05.02-07, 05-09.05-10, 06-01.06-03, 03-10.03-11
#during fall: 02-04.02-05, 02-05.02-07, 02-04.02-07, 03-16.03-24

group.dyads <- c("02-01.03-02","02-04.02-05","02-04.02-07","02-05.02-04","02-05.02-07","02-06.02-07","02-07.02-04","02-07.02-05",
                 "02-07.02-06","03-02.02-01","03-07.03-20","03-10.03-11","03-11.03-10","03-16.03-19","03-16.03-24","03-19.03-16",
                 "03-20.03-07","03-24.03-16","04-09.04-13","04-13.04-09","05-07.05-15","05-09.05-10","05-10.05-09","05-15.05-07",
                 "06-01.06-03","06-03.06-01","09-01.09-02","09-02.09-01","11-04.11-05","11-05.11-04")
group.dyads <- gsub(group.dyads, pattern="-", replacement="")
group.dyads <- gsub(group.dyads, pattern="[.]", replacement="_")

summary.locs$ij <- paste(summary.locs$i, summary.locs$j, sep="_")

removed.dyads <- summary.locs[summary.locs$ij %in% group.dyads,]
removed.dyads %>% group_by(d.sex) %>% dplyr::summarize(sum(d25), ndyads=length(ij))
#   d.sex `sum(d25)` ndyads
# 1 FF           870     16
# 2 FM             1      1
# 3 MF             1      1
# 4 MM          1492     12

kept.dyads <- summary.locs[!summary.locs$ij %in% group.dyads,]
kept.dyads %>% group_by(d.sex) %>% dplyr::summarize(sum(d25), ndyads=length(ij))

#   d.sex `sum(d25)` ndyads
# 1 FF          1228     54
# 2 FM           254     49
# 3 MF           254     49
# 4 MM          1224     70

#########################
#MRQAP

#add dyad characteristics to frequency of contact
data <- read.table(file="Dyadstructure_spatsoc_2022-03-18.txt", sep=",", header=T) #from spatsoc 222 dyads with contacts, 558 dyads total

data$rate <- data$d25/data$daysactive.num #change to spatsoc numbers
data$QGM2 <- data$QGM^2

#dummy coding sex
data$ff <- 0
data[data$d.sex == "FF",]$ff <- 1
data$mm <- 0
data[data$d.sex == "MM",]$mm <- 1

#convert from long to wide format matrix for each metric
#dependent variable
# freq.m <- as.matrix(dcast(data, i ~ j, value.var="Freq")[,-1]) #frequency matrix
freq.m <- as.matrix(dcast(data, i ~ j, value.var="d25")[,-1]) #frequency matrix
rate.m <- as.matrix(dcast(data, i ~ j, value.var="rate")[,-1]) #rate matrix

#independent variables
ff.m <- as.matrix(dcast(data, i ~ j, value.var="ff")[,-1])
mm.m <- as.matrix(dcast(data, i ~ j, value.var="mm")[,-1])
#d.sex.m <- as.matrix(dcast(data, i ~ j, value.var="d.sex")[,-1])

qgm.m <- as.matrix(dcast(data, i ~ j, value.var="QGM")[,-1])
qgm2.m <- as.matrix(dcast(data, i ~ j, value.var="QGM2")[,-1])

di.m <- as.matrix(dcast(data, i ~ j, value.var="DI")[,-1])
di.w.m <- as.matrix(dcast(data, i ~ j, value.var="W.AvgDI")[,-1])
di.l.m <- as.matrix(dcast(data, i ~ j, value.var="L.AvgDI")[,-1])
di.f.m <- as.matrix(dcast(data, i ~ j, value.var="F.AvgDI")[,-1])

vi.m <- as.matrix(dcast(data, i ~ j, value.var="VI")[,-1])
vi.w.m <- as.matrix(dcast(data, i ~ j, value.var="W.AvgVI")[,-1])
vi.l.m <- as.matrix(dcast(data, i ~ j, value.var="L.AvgVI")[,-1])
vi.f.m <- as.matrix(dcast(data, i ~ j, value.var="F.AvgVI")[,-1])

sex.m <- as.matrix(dcast(data, i ~ j, value.var="sexmatch")[,-1])
age.m <- as.matrix(dcast(data, i ~ j, value.var="ageclassmatch")[,-1])

#############
#run models

###########
#dyad characteristics

#VI vs. DI
mrqap.dsp(vi.l.m ~ di.l.m) #these are highly correlated

#VI
vi <- mrqap.dsp(rate.m ~ vi.m); vi.w <- mrqap.dsp(rate.m ~ vi.w.m); vi.l <- mrqap.dsp(rate.m ~ vi.l.m); vi.f <- mrqap.dsp(rate.m ~ vi.f.m)
vi$coefficients; vi.w$coefficients; vi.l$coefficients; vi.f$coefficients
#vi.l is best

#DI
di <- mrqap.dsp(rate.m ~ di.m); di.w <- mrqap.dsp(rate.m ~ di.w.m); di.l <- mrqap.dsp(rate.m ~ di.l.m); di.f <- mrqap.dsp(rate.m ~ di.f.m)
di$coefficients; di.w$coefficients; di.l$coefficients; di.f$coefficients
#di.l is best

qgm <- mrqap.dsp(rate.m ~ qgm.m)
qgm2 <- mrqap.dsp(rate.m ~ qgm.m + qgm2.m)

#d.sex <- mrqap.dsp(rate.m ~ d.sex.m) #can't do this because these are categorical, need to dummy code
d.sex <- mrqap.dsp(rate.m ~ ff.m + mm.m)

# d.glob <- mrqap.dsp(rate.m ~ ff.m + mm.m + vi.l.m + qgm.m + qgm2.m)
d.glob <- mrqap.dsp(rate.m ~ ff.m + mm.m + vi.l.m + qgm.m + qgm2.m)
d.glob

# MRQAP with Double-Semi-Partialing (DSP)
# 
# Formula:  rate.m ~ ff.m + mm.m + vi.l.m + qgm.m + qgm2.m 
# 
# Coefficients:
#           Estimate    P(β>=r) P(β<=r) P(|β|<=|r|)
# intercept -0.01196813 0.024   0.976   0.024      
# ff.m       0.02039911 1.000   0.000   0.000      
# mm.m       0.01519254 0.995   0.005   0.009      
# vi.l.m     0.11618324 1.000   0.000   0.000      
# qgm.m      0.01195665 0.768   0.232   0.471      
# qgm2.m    -0.10390907 0.087   0.913   0.155      
# 
# Residual standard error: 0.03121 on 163 degrees of freedom
# F-statistic: 32.75 on 5 and 163 degrees of freedom, p-value:     0 
# Multiple R-squared: 0.5011 	Adjusted R-squared: 0.4858 
# AIC: -631.7559 

d.glob.beta <- mrqap.dsp(rate.m ~ ff.m + mm.m + vi.l.m + qgm.m + qgm2.m, test.statistic = "beta") #produces the same results as test.statistic="estimates"

############
#homophily
sex <- mrqap.dsp(rate.m ~ sex.m)
age <- mrqap.dsp(rate.m ~ age.m)
sex.age <- mrqap.dsp(rate.m ~ sex.m + age.m)
sex.age

# MRQAP with Double-Semi-Partialing (DSP)
# 
# Formula:  rate.m ~ sex.m + age.m 
# 
# Coefficients:
#           Estimate    P(β>=r) P(β<=r) P(|β|<=|r|)
# intercept 0.002002357 0.605   0.395   0.817      
# sex.m     0.022476466 1.000   0.000   0.000      
# age.m     0.005952903 0.704   0.296   0.544      
# 
# Residual standard error: 0.0416 on 276 degrees of freedom
# F-statistic:  10.5 on 2 and 276 degrees of freedom, p-value: 4.026e-05 
# Multiple R-squared: 0.07071 	Adjusted R-squared: 0.06398 
# AIC: -582.2491 

#
glob <- mrqap.dsp(rate.m ~ ff.m + mm.m + vi.l.m + qgm.m + qgm2.m + age.m)
glob

############
#individual characteristics
data$i.sex <- 0
data[data$SEX.i == "M",]$i.sex <- 1

i.sex.m <- as.matrix(dcast(data, i ~ j, value.var="i.sex")[,-1])
i.age.m <- as.matrix(dcast(data, i ~ j, value.var="AGE.i")[,-1])

sex.i <- mrqap.dsp(rate.m ~ i.sex.m)
age.i <- mrqap.dsp(rate.m ~ i.age.m)

sex.age.i <- mrqap.dsp(rate.m ~ i.sex.m + i.age.m)
sex.age.i

# MRQAP with Double-Semi-Partialing (DSP)
# 
# Formula:  rate.m ~ i.sex.m + i.age.m 
# 
# Coefficients:
#           Estimate    P(β>=r) P(β<=r) P(|β|<=|r|)
# intercept 0.010658040 0.942   0.058   0.084      
# i.sex.m   0.002455240 0.675   0.325   0.659      
# i.age.m   0.001252548 0.881   0.119   0.212      
# 
# Residual standard error: 0.04299 on 276 degrees of freedom
# F-statistic: 1.038 on 2 and 276 degrees of freedom, p-value: 0.3555 
# Multiple R-squared: 0.007465 	Adjusted R-squared: 0.0002731 
# AIC: -791.9164 


##################################################################################################################################################################
#########################
#3 when do direct contacts occur? Visualization
#########################

########
#load packages
########
require(lubridate)

path <- "C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper"
setwd(path)

########
#read in data
########

direct <- read.table(file="direct_contact_locs_2022-03-11.txt", sep=",", header=T, colClasses = c("character","numeric","character","character","numeric","numeric","character","character",
                                                                                                  "numeric","numeric","numeric","character","numeric","numeric","numeric","numeric","character","numeric","numeric","numeric"))
direct$MSDATE. <- as.POSIXct(direct$MSDATE., format="%Y.%m.%d", tz="UTC")
direct$yr <- year(direct$MSDATE.)
direct$month <- month(direct$MSDATE.)
nrow(direct) #5104 direct contact

sheep <- unique(direct[,c("ANIMALID.","SEX.","GROUP.")])

######
#condense by timegroup and sex type
direct.ddply <- ddply(direct[order(direct$ANIMALID.),], .(timegroup,DOY,month,yr), summarize, num.sheep=length(timegroup), ij=paste(ANIMALID., collapse="_"), sex.ij=paste(SEX., collapse=""))
nrow(direct.ddply) #2315

unique(direct.ddply[nchar(direct.ddply$sex) > 2,]$sex)
direct.ddply$sex <- direct.ddply$sex.ij
direct.ddply[direct.ddply$sex %in% c("FFF","FFFF"),]$sex <- "FF"
direct.ddply[direct.ddply$sex %in% c("MMM","MMMM"),]$sex <- "MM"
direct.ddply[nchar(direct.ddply$sex) > 2,]$sex <- "FM"
direct.ddply[direct.ddply$sex == "MF",]$sex <- "FM"

group.dyads <- c("02-01.03-02","02-04.02-05","02-04.02-07","02-05.02-04","02-05.02-07","02-06.02-07","02-07.02-04","02-07.02-05",
                 "02-07.02-06","03-02.02-01","03-07.03-20","03-10.03-11","03-11.03-10","03-16.03-19","03-16.03-24","03-19.03-16",
                 "03-20.03-07","03-24.03-16","04-09.04-13","04-13.04-09","05-07.05-15","05-09.05-10","05-10.05-09","05-15.05-07",
                 "06-01.06-03","06-03.06-01","09-01.09-02","09-02.09-01","11-04.11-05","11-05.11-04")
group.dyads <- gsub(group.dyads, pattern="-", replacement="")
group.dyads <- gsub(group.dyads, pattern="[.]", replacement="_")

direct.ddply <- direct.ddply[!direct.ddply$ij %in% group.dyads,] #remove within group contacts

#sum direct contacts by month and dyad sex
direct.bymonth <- ddply(direct.ddply, .(month, yr, sex), summarize, numcontacts=length(timegroup), ndyads=length(unique(ij)))
direct.bymonth$perdyad <- direct.bymonth$numcontacts/direct.bymonth$ndyads 

direct.bymonth$month <- factor(direct.bymonth$month)
direct.bymonth$sex <- factor(direct.bymonth$sex, levels=c("FF","MM","FM"))

#total contacts per month
ggplot(data=direct.bymonth, aes(x=month, y=numcontacts, fill=sex, col=sex)) + geom_boxplot(lwd=0.5) + #facet_wrap(~sex, ncol=1, scales="free_y") + 
  # ggplot(data=direct.bymonth, aes(x=month, y=numcontacts, group=month, col=sex)) + geom_boxplot(lwd=0.75) + facet_wrap(~sex, ncol=1, scales="free_y") + 
  scale_color_manual(values=c("#B80C09","#084C61","black"), name="dyad sex") +
  scale_fill_manual(values=c("#F9B5AC","#B4D2E7","grey50"), name="dyad sex") +
  ylab("number of direct contacts") +
  theme_bw(base_size=20) + theme(panel.grid=element_blank(), legend.position = c(0.9, 0.8)) #, legend.position = "topright")
ggsave(filename = "Figs_Final/Fig3.png", height=6, width=10, units="in", dpi=400)

#contacts per month per dyad
ggplot(data=direct.bymonth, aes(x=month, y=perdyad, fill=sex, col=sex)) + geom_boxplot(lwd=0.5) + #facet_wrap(~sex, ncol=1, scales="free_y") + 
  # ggplot(data=direct.bymonth, aes(x=month, y=numcontacts, group=month, col=sex)) + geom_boxplot(lwd=0.75) + facet_wrap(~sex, ncol=1, scales="free_y") + 
  scale_color_manual(values=c("#B80C09","#084C61","black"), name="dyad sex") +
  scale_fill_manual(values=c("#F9B5AC","#B4D2E7","grey50"), name="dyad sex") +
  ylab("number of direct contacts per dyad") + coord_cartesian(ylim=c(NA,10)) +
  theme_bw(base_size=20) + theme(panel.grid=element_blank(), legend.position = c(0.9, 0.8)) #, legend.position = "topright")
ggsave(filename = "Figs_Final/Fig3_2.png", height=6, width=10, units="in", dpi=400)

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

#################################################
#4a when do contacts occur? Modelling
#################################################

########
#load packages
########
require(lme4)
require(AICcmodavg)
require(dplyr)

path <- "C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper"
setwd(path)

########
#load data
########
#data 
data.contact <- read.table("C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper/ContactsRSF_spatsoc_2022-03-23.txt", sep=",", header=T)
data.contact %>% group_by(d.sex) %>% dplyr::summarize(length(i), sum(used))
1011 + 403 + 956 #2370 contacts

data.contact[data.contact$used == 1,] %>% group_by(d.sex) %>% dplyr::summarize(length(unique(ij)))
52+93+67 #212 dyads

data.contact %>% group_by(d.sex) %>% dplyr::summarize(length(unique(ij)))

########
#need to reduce number of available points
########
data <- data.contact
data <- data[data$Prox < 500,] #another sheep within 500 m? 500 m = 3rd quantile of step lengths

data %>% group_by(d.sex) %>% dplyr::summarize(length(i), sum(used))
data %>% group_by(d.sex) %>% dplyr::summarize(length(unique(ij)))

data$ij <- factor(data$ij)
data$month <- factor(data$month)
data$HOD <- factor(data$HOD)

data.ff <- data[data$d.sex == "FF",]
data.mm <- data[data$d.sex == "MM",]
data.mf <- data[data$d.sex == "MF",]

########
#when do contacts occur?
########

data <- data.ff
data <- data.mm
data <- data.mf

family = binomial(link="logit")
null <- glmer(used ~ (1|ij), data=data, family=family) #random effect = ij (dyad)

#univariate models
around <- glmer(used ~ Freq400m + (1|ij), data=data, family=family)
hod <- glmer(used ~ HOD + (1|ij), data=data, family=family)
  ss <- getME(hod, c("theta","fixef"))
  hod <- update(hod, start=ss, control=glmerControl(optCtrl=list(maxfun=2e6)))
light <- glmer(used ~ light + (1|ij), data=data, family=family)
season <- glmer(used ~ season + (1|ij), data=data, family=family)
speed <- glmer(used ~ scale(speed_mhr) + (1|ij), data=data, family=family)
month <- glmer(used ~ month  + (1|ij), data=data, family=family)
  ss <- getME(month, c("theta","fixef"))
  month <- update(month, start=ss, control=glmerControl(optCtrl=list(maxfun=2e6)))
  
#global model
global <- glmer(used ~ Freq400m + scale(speed_mhr) + light + season + (1|ij), data=data, family=family)
summary(global)

#rank by aicc
objs =  mget(ls())
modname <- names(Filter(function(i) inherits(i, "glmerMod"), objs))
aictab(lapply(modname, function(x) eval(parse(text=x))), modname)

#coefficients from global model
ff.coef <- data.frame(summary(global)$coefficients)
ff.coef$d.sex <- "FF"
ff.coef$coefficient <- row.names(ff.coef)

mm.coef <- data.frame(summary(global)$coefficients)
mm.coef$d.sex <- "MM"
mm.coef$coefficient <- row.names(mm.coef)

mf.coef <- data.frame(summary(global)$coefficients)
mf.coef$d.sex <- "MF"
mf.coef$coefficient <- row.names(mf.coef)

coef <- rbind(ff.coef, mm.coef, mf.coef)
coef$coefficient <- factor(coef$coefficient, labels=c("Number collars nearby","Diurnal","Nocturnal","Speed","Summer","Winter","(Intercept)"), levels=c("Freq400m","lightday","lightnight","scale(speed_mhr)","seasonsummer","seasonwinter","(Intercept)"))

coef.plot <- ggplot(data=coef[!coef$coefficient %in% c("(Intercept)"),], aes(x=Estimate, y=coefficient, col=d.sex)) + 
  geom_vline(aes(xintercept=0), lwd=0.5, lty="dashed", col="black") + 
  geom_point(size=3) +
  geom_errorbar(aes(xmin=Estimate - Std..Error, xmax=Estimate + Std..Error), lwd=1, width=0.25) +
  scale_color_manual(values=c("pink","black","lightblue")) +
  theme_bw(base_size=15) + theme(panel.grid=element_blank(), legend.position="none")
coef.plot
ggsave(coef.plot, filename="Figs_Final/directcontact_timing.tiff", height=4, width=8, units="in", compression="lzw", dpi=350)

coef %>% mutate_if(is.numeric, round, 3)

#                   Estimate Std..Error z.value Pr...z.. d.sex           coefficient
# (Intercept)         -3.233      0.153 -21.167    0.000    FF           (Intercept)
# Freq400m             0.658      0.063  10.492    0.000    FF Number collars nearby***
# scale(speed_mhr)     0.069      0.034   2.020    0.043    FF                 Speed
# lightday             0.007      0.085   0.078    0.938    FF               Diurnal
# lightnight           0.081      0.086   0.943    0.346    FF             Nocturnal
# seasonsummer         0.244      0.109   2.235    0.025    FF                Summer
# seasonwinter         0.581      0.092   6.303    0.000    FF                Winter
# (Intercept)1        -2.852      0.143 -19.925    0.000    MM           (Intercept)
# Freq400m1            0.814      0.069  11.826    0.000    MM Number collars nearby***
# scale(speed_mhr)1   -0.144      0.046  -3.133    0.002    MM                 Speed
# lightday1           -0.165      0.086  -1.922    0.055    MM               Diurnal
# lightnight1         -0.167      0.090  -1.868    0.062    MM             Nocturnal
# seasonsummer1        0.035      0.110   0.318    0.750    MM                Summer
# seasonwinter1        0.359      0.114   3.146    0.002    MM                Winter
# (Intercept)2        -3.749      0.212 -17.654    0.000    MF           (Intercept)
# Freq400m2            0.691      0.076   9.133    0.000    MF Number collars nearby***
# scale(speed_mhr)2   -0.229      0.079  -2.920    0.004    MF                 Speed
# lightday2           -0.443      0.142  -3.112    0.002    MF               Diurnal
# lightnight2          0.219      0.125   1.748    0.080    MF             Nocturnal
# seasonsummer2       -0.224      0.281  -0.798    0.425    MF                Summer
# seasonwinter2        0.379      0.165   2.300    0.021    MF                Winter

##############################################################################################################################################################################
#################################################
#4b where are bighorn? general habitat use
#step selection function
#################################################

#load packages
require(lme4)
require(mclogit) #Thurfjell et al. 2014 recommended

require(AICcmodavg)
require(reshape2) #for melt

require(lubridate)
require(ggplot2)
require(ggpubr)

require(dplyr)
require(tidyr) #for separate

#install.packages("amt")
require(amt) #to plot predicted probabilities of use from fit_clogit objects
require(ggpubr)

#############
#load data
path.local <- "C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper"
setwd(path.local)

d <- read.table(file="SSF_data_2021-01-14.txt", sep=",", header=T)
d$step_id_ <- as.integer(d$set)

table(d[,c("SEX.","used")])
# SEX.     0      1
#F    261033  52348
#M    270103  54131

############
#variable categories
#linear variables
vars <- names(data2) %>% grep(pattern="\\.s", value=T)
#vars <- c(vars, "cc","vrm","sri","ndvi","irg","mgoat")

vars <- data.frame(s=vars)
vars <- separate(vars, col="s", into=c("var"), remove=F, sep="\\.")

pvars <- data.frame(p=names(data2) %>% grep(pattern="\\.p", value=T))
pvars <- separate(pvars, col="p", into="var", remove=F, sep="\\.")

qvars <- data.frame(q=names(data2) %>% grep(pattern="\\.q", value=T))
qvars <- separate(qvars, col="q", into="var", remove=F, sep="\\.")

v <- merge(vars, pvars, by="var", all.x=T)
v <- merge(v, qvars, by="var", all.x=T)
v$cat <- c("survival","survival","resource","resource","resource",
           "disturbance","resource","disturbance","disturbance","resource",
           "resource","disturbance","resource","survival","resource",
           "resource","survival")

##################
#vars to keep in simplified model
#e.p, change sri.q to sri.s since doesn't make sense otherwise, dh2o_per.q, dmlick.p, irg.q, ndvi.q, exclude swe.p unless sign switch in global model
#desc.s, cc.s, vrm.p exclude slope.q since correlated with desc
#exclude mgoat.p, dtr.s, drd.s, dheli.s unless sign switch in global model

########
#compare figures of global and global.s
#global, no interactions
g <- d %>% fit_clogit(used ~ 
                        e.p+e.s+sri.q+sri.s+dh2o_per.q+dh2o_per.s+
                        dmlick.p+dmlick.s+irg.q+irg.s+ndvi.q+ndvi.s+swe.p+swe.s+
                        vrm.p+vrm.s+desc.p+desc.s+cc.q+cc.s+
                        #dtr.p+dtr.s+dheli.p+dheli.s+drd.p+drd.s +
                        strata(step_id_), model=T)

#global simplified, for m+f, f only, m only
g.s <- d %>% fit_clogit(used ~
                          e.p+e.s+sri.s+dh2o_per.q+dh2o_per.s+
                          dmlick.p+dmlick.s+irg.q+irg.s+ndvi.s+swe.s+ #check swe effect #removed swe.p since doesnt make sense
                          vrm.p+vrm.s+desc.s+cc.s+
                          strata(step_id_), model=T)
g.s.f <- d %>% filter(SEX. == "F") %>%
  fit_clogit(used ~ e.p+e.s+sri.s+dh2o_per.q+dh2o_per.s+
               dmlick.p+dmlick.s+irg.q+irg.s+ndvi.s+swe.s+ #check swe effect
               vrm.p+vrm.s+desc.s+cc.s+
               strata(step_id_), model=T)
g.s.m <- d %>% filter(SEX. == "M") %>%
  fit_clogit(used ~ e.p+e.s+sri.s+dh2o_per.q+dh2o_per.s+
               dmlick.p+dmlick.s+irg.q+irg.s+ndvi.s+swe.s+ #check swe effect
               vrm.p+vrm.s+desc.s+cc.s+
               strata(step_id_), model=T)

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

#function to plot coefficients from clogit models
plotclogit <- function(b, prefix, mod1, mod2, mod3=NULL, col1, col2, col3=NULL)
  #for(b in unique(gsub(v, pattern=".[spq]$", replacement="")))
{
  print(b)
  #x1 <- data.frame(array(data=0, dim=c(nrow(d), length(f$model$coefficients))))
  x1 <- data.frame(array(data=0, dim=c(5000, length(mod1$model$coefficients))))
  names(x1) <- names(mod1$model$coefficients)
  bvars <- grep(names(mod1$model$coefficients), pattern=paste("^", b, sep=""), value=T)
  bvars <- grep(bvars, pattern=":", invert=T, value=T)
  x1[,bvars] <- d[sample(nrow(d), size=5000, replace=F),][,bvars]
  #x1[,c("e.p","e.s")] <- d[sample(nrow(d), size=5000, replace=F),][,c("e.p","e.s")]
  x2 <- data.frame(array(data=0, dim=c(1,length(names(x1)))))
  names(x2) <- names(x1)
  logRSS1 <- log_rss(object=mod1, x1=x1, x2=x2, ci="se")
  logRSS2 <- log_rss(object=mod2, x1=x1, x2=x2, ci="se")
  logRSS3 <- NULL
  logRSS3 <- try(log_rss(object=mod3, x1=x1, x2=x2, ci="se"))
  
  
  p <- ggplot() + 
    #mod1
    geom_ribbon(data=logRSS1$df, aes(x=eval(parse(text=paste(b,".s_x1",sep=""))), ymin = lwr, ymax = upr), fill = col1, alpha=0.25) +
    geom_line(data=logRSS1$df, aes(x=eval(parse(text=paste(b,".s_x1",sep=""))), y=log_rss), lwd=1, lty="dashed", col=col1) + 
    #mod2
    geom_ribbon(data=logRSS2$df, aes(x=eval(parse(text=paste(b,".s_x1",sep=""))), ymin = lwr, ymax = upr), fill = col2, alpha=0.25) +
    geom_line(data=logRSS2$df, aes(x=eval(parse(text=paste(b,".s_x1",sep=""))), y=log_rss), lwd=1, lty="dashed", col=col2) + 
    xlab(paste(b,".s", sep="")) +
    theme_bw()
  if(!is.null(logRSS3))
  {
    p <- p +
      #mod3
      geom_ribbon(data=logRSS3$df, aes(x=eval(parse(text=paste(b,".s_x1",sep=""))), ymin = lwr, ymax = upr), fill = col3, alpha=0.25) +
      geom_line(data=logRSS3$df, aes(x=eval(parse(text=paste(b,".s_x1",sep=""))), y=log_rss), lwd=1, lty="dashed", col=col3)
  }
  
  ggsave(p, filename=paste("C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper/Figures_SSF/", prefix, "_", b, ".tiff", sep=""), height=5, width=5, units="in", dpi=300, compression="lzw")
}

plotclogitint <- function(b1, b2,prefix, mod)
  #for(b in unique(gsub(v, pattern=".[spq]$", replacement="")))
{
  print(b1); print(b2)
  x1 <- data.frame(array(data=0, dim=c(5000, length(mod$model$coefficients))))
  names(x1) <- names(mod$model$coefficients)
  b1vars <- grep(names(mod$model$coefficients), pattern=paste("^", b1, sep=""), value = T)
  b1vars <- grep(b1vars, pattern=":", invert=T, value=T)
  b2vars <- grep(names(mod$model$coefficients), pattern=paste("^", b2, sep=""), value = T)
  newdata <- d[sample(nrow(d), size=5000, replace=F),]
  x1[,b1vars] <- newdata[,b1vars]
  #x1[,c("e.p","e.s")] <- d[sample(nrow(d), size=5000, replace=F),][,c("e.p","e.s")]
  x2 <- data.frame(array(data=0, dim=c(1,length(names(x1)))))
  names(x2) <- names(x1)
  
  x1[,b2vars] <- -1
  logRSS1 <- log_rss(object=mod, x1=x1, x2=x2, ci="se")
  logRSS1$df$int <- "-1 SD"
  x1[,b2vars] <- 0
  logRSS2 <- log_rss(object=mod, x1=x1, x2=x2, ci="se")
  logRSS2$df$int <- "average"
  x1[,b2vars] <- 1
  logRSS3 <- log_rss(object=mod, x1=x1, x2=x2, ci="se")
  logRSS3$df$int <- "+1 SD"
  
  
  p <- ggplot() + 
    #logRSS1
    geom_ribbon(data=logRSS1$df, aes(x=eval(parse(text=paste(b1,".s_x1",sep=""))), ymin = lwr, ymax = upr), fill = "pink", alpha=0.25) +
    geom_line(data=logRSS1$df, aes(x=eval(parse(text=paste(b1,".s_x1",sep=""))), y=log_rss, col=int), lwd=1, lty="dashed") + 
    #logRSS2
    geom_ribbon(data=logRSS2$df, aes(x=eval(parse(text=paste(b1,".s_x1",sep=""))), ymin = lwr, ymax = upr), fill = "grey80", alpha=0.25) +
    geom_line(data=logRSS2$df, aes(x=eval(parse(text=paste(b1,".s_x1",sep=""))), y=log_rss, col=int), lwd=1, lty="dashed") +
    #logRSS3
    geom_ribbon(data=logRSS3$df, aes(x=eval(parse(text=paste(b1,".s_x1",sep=""))), ymin = lwr, ymax = upr), fill = "steelblue", alpha=0.25) +
    geom_line(data=logRSS3$df, aes(x=eval(parse(text=paste(b1,".s_x1",sep=""))), y=log_rss, col=int), lwd=1, lty="dashed") +
    scale_color_manual(values=c("pink","steelblue","black"), name=b2) +
    xlab(paste(b1,".s", sep="")) +
    #labs(fill = b2) +
    theme_bw()
  ggsave(p, filename=paste("Figures_SSF/", prefix, "_", b, ".tiff", sep=""), height=5, width=5, units="in", dpi=300, compression="lzw")
}

v1 <- unique(gsub(names(g.s$model$coefficients), pattern=".[spq]$", replacement=""))
v1 <- grep(v1, pattern=":", invert=T, value=T)
map(v1, plotclogit, prefix="g_s", mod1=g.s, mod2=g.s.f, mod3=g.s.m, col1="black",col2="steelblue",col3="pink")

v1 <- names(g$model$coefficients)
map(unique(gsub(v1, pattern=".[spq]$", replacement="")), plotclogit, prefix="test", mod1=g, mod2=g.s)


############
#plot betas of global int s models
############
#globalbeta <- as.data.frame(coefficients(summary(g.int.s)))
globalbeta <- as.data.frame(coefficients(summary(g.s)))
globalbeta$sex <- "all"
globalbeta$variable <- row.names(globalbeta)

#globalbeta.m <- as.data.frame(coefficients(summary(global.int.m)))
globalbeta.m <- as.data.frame(coefficients(summary(g.s.m)))
globalbeta.m$sex <- "M"
globalbeta.m$variable <- row.names(globalbeta.m)

#globalbeta.f <- as.data.frame(coefficients(summary(global.int.f)))
globalbeta.f <- as.data.frame(coefficients(summary(g.s.f)))
globalbeta.f$sex <- "F"
globalbeta.f$variable <- row.names(globalbeta.f)

globalbeta <- rbind(globalbeta, globalbeta.m, globalbeta.f) #combine data from all 3 models

globalbeta <- separate(globalbeta, col="variable", into=c("var","type"), sep="\\.", remove=F)
#globalbeta[is.na(globalbeta$type),]$type <- "s"
#globalbeta[grep(globalbeta$variable, pattern=":"),]$var <- "interaction"
#globalbeta[grep(globalbeta$variable, pattern=":"),]$type <- "i"
globalbeta$var2 <- gsub(globalbeta$variable, pattern="\\.s", replacement="")
#names(globalbeta) <- c("Estimate","StdError","zvalue","Pr(>|z|)","sex","variable","var","type","var2") #for mclogit models
names(globalbeta) <- c("Estimate","exp_coef","StdError","zvalue","Pr(>|z|)","sex","variable","var","type","var2") #for fit_clogit models
#globalbeta$modeltype <- "global.int.s"
globalbeta$modeltype <- "global.noint.s"
globalbeta <- merge(globalbeta, v[,c("var","cat")], by="var", all.x=T)
#globalbeta[is.na(globalbeta$cat),]$cat <- c("resource","resource","survival","survival","resource","resource","resource","resource","survival") #assign categories to interaction terms
#globalbeta[is.na(globalbeta$cat),]$cat <- c("resource","survival","resource","resource","resource","resource","resource","survival","survival")

betaorder <- unique(globalbeta[,c("cat","var","var2","type")])
betaorder$cat <- factor(betaorder$cat, levels=c("resource","survival","disturbance"))
betaorder$type <- factor(betaorder$type, levels=c("s","q","p","i"))

betaorder <- betaorder %>%
  #group_by(cat) %>%
  arrange(desc(cat, var, type, var2))

#simple global with interactions
# globalbeta$var2 <- factor(globalbeta$var2, levels=c("dheli",
#                                                     "cc","desc","vrm.p","vrm","desc:cc",
#                                                     "dh2o_per.q","dh2o_per","dmlick.p","dmlick",
#                                                     "e.p","e","irg.q","irg","ndvi","sri","swe","e.p:sri","e:sri"))
#simple global, no interactions
globalbeta$var2 <- factor(globalbeta$var2, levels=c("dheli",
                                                    "cc","desc","vrm.p","vrm",
                                                    "dh2o_per.q","dh2o_per","dmlick.p","dmlick",
                                                    "e.p","e","irg.q","irg","ndvi","sri","swe"))


# write.table(globalbeta, file="C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper/coefficients_globalsbetas_2021-01-14.txt", sep=",")

# g <- ggplot(data=globalbeta, aes(x=Estimate, y=var2, col=sex)) + 
#   geom_vline(aes(xintercept=0), lty="dashed") + 
#   geom_point(aes(shape=type), size=3, position=position_dodge2(width=0.5)) + 
#   geom_errorbar(aes(xmin=Estimate-StdError, xmax=Estimate+StdError), width=0.5, position=position_dodge2(width=0.5)) +
#   scale_shape_manual(values=c(18,16,17,15)) +
#   scale_color_manual(values=c("black","pink","lightblue")) +
#   ylab("variable") + xlab("beta estimate") + 
#   facet_wrap(~cat, ncol=1, scales="free") +
#   theme_bw(base_size = 15)
# g
# ggsave(g, filename="C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper/coefficients_global_2020-12-07.tiff", compression="lzw", height=10, width=8, units="in", dpi=300)

#resources
r <- ggplot(data=globalbeta[globalbeta$cat=="resource",], aes(x=Estimate, y=var2, col=sex)) + 
  geom_vline(aes(xintercept=0), lty="dashed") + 
  geom_point(aes(shape=type), size=3, position=position_dodge2(width=0.5)) + 
  geom_errorbar(aes(xmin=Estimate-StdError, xmax=Estimate+StdError), width=0.5, position=position_dodge2(width=0.5)) +
  scale_shape_manual(values=c(18,16,17,15)) +
  scale_color_manual(values=c("black","pink","lightblue")) +
  ylab("variable") + 
  theme_bw(base_size = 15) #, axis.x=element_blank(), plot.margin = unit(c(0,0.1,0.1,0), "cm")
#survival
s <- ggplot(data=globalbeta[globalbeta$cat=="survival",], aes(x=Estimate, y=var2, col=sex)) + 
  geom_vline(aes(xintercept=0), lty="dashed") + 
  geom_point(aes(shape=type), size=3, position=position_dodge2(width=0.5)) + 
  geom_errorbar(aes(xmin=Estimate-StdError, xmax=Estimate+StdError), width=0.5, position=position_dodge2(width=0.5)) +
  scale_shape_manual(values=c(18,16,17,15)) +
  scale_color_manual(values=c("black","pink","lightblue")) +
  ylab("variable") + xlab("") + 
  theme_bw(base_size = 15)
#disturbance
# d <- ggplot(data=globalbeta[globalbeta$cat=="disturbance",], aes(x=Estimate, y=var2, col=sex)) + 
#   geom_vline(aes(xintercept=0), lty="dashed") + 
#   geom_point(aes(shape=type), size=3, position=position_dodge2(width=0.5)) + 
#   geom_errorbar(aes(xmin=Estimate-StdError, xmax=Estimate+StdError), width=0.5, position=position_dodge2(width=0.5)) +
#   scale_shape_manual(values=c(18,16,17,15)) +
#   scale_color_manual(values=c("black","pink","lightblue")) +
#   ylab("variable") + xlab("beta estimate") + 
#   theme_bw(base_size = 15)

date <- Sys.Date()
date <- "2021-01-14"
# ggsave(ggarrange(r,s,d, ncol=1, align="v", heights=c(5,2,1), common.legend=T), #, labels="AUTO"), 
ggsave(ggarrange(r,s, ncol=1, align="v", heights=c(5,2.5), common.legend=T), #, labels="AUTO"), 
       filename=paste("C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper/coefficients_global_s_", date,".tiff", sep=""),
       compression="lzw", height=10, width=8, units="in", dpi=300)
	   
	   
##########################################################################################################################################################

#######################################
#4c where do direct contact occur
#######################################

########
#load packages
########
require(lme4)
require(AICcmodavg)
require(dplyr)

path <- "C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper"
setwd(path)

########
#load data
########
#data 
data.contact <- read.table("C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper/ContactsRSF_spatsoc_2022-03-23.txt", sep=",", header=T)
data.contact %>% group_by(d.sex) %>% dplyr::summarize(length(i), sum(used))

# d.sex `length(i)` `sum(used)`
# 1 FF          25580        1011
# 2 MF          36827         403
# 3 MM          36553         956

# ggsave(filename="Figures/DOY_vs_NDVI.tiff", ggplot(data=data, aes(x=DOY, y=NDVI, col=d.sex)) + geom_point() + theme_bw(base_size=20) + facet_wrap(~d.sex, nrow=3), height=8, width=8, units="in", dpi=300, compression="lzw")
# ggsave(filename="Figures/DOY_vs_IRG.tiff", ggplot(data=data, aes(x=DOY, y=IRG, col=d.sex)) + geom_point() + theme_bw(base_size=20) + facet_wrap(~d.sex, nrow=3), height=8, width=8, units="in", dpi=300, compression="lzw")
# ggsave(filename="Figures/DOY_vs_SWE.tiff", ggplot(data=data, aes(x=DOY, y=SWE, col=d.sex)) + geom_point() + theme_bw(base_size=20) + facet_wrap(~d.sex, nrow=3), height=8, width=8, units="in", dpi=300, compression="lzw")
# ggsave(filename="Figures/DOY_vs_dmlicks.tiff", ggplot(data=data, aes(x=DOY, y=d_mlicks, col=d.sex)) + geom_point() + theme_bw(base_size=20) + facet_wrap(~d.sex, nrow=3), height=8, width=8, units="in", dpi=300, compression="lzw")

########
#need to reduce number of available points
########
data <- data.contact
data <- data[data$Prox < 500,] #another sheep within 500 m? 500 m = 3rd quantile of step lengths

censor <- ddply(data, .(ij, used, d.sex), summarize, Freq=length(FixStartDateTime))
censor <- censor[censor$used == 1,]

table(data[,c("d.sex","used")])

data$ij <- factor(data$ij)
data$month <- factor(data$month)
data$HOD <- factor(data$HOD)

#rename so same as RSF data
data.n <- names(data)
data.n[data.n == "d_h2o_25km"] <- "dh2o"
data.n[data.n == "d_h2o_per"] <- "dh2o_per"
data.n[data.n == "d_mlicks"] <- "dmlick"
data.n[data.n == "desc27"] <- "desc"
data.n[data.n == "slope"] <- "slope"
data.n[data.n == "DEM"] <- "e"
data.n[data.n == "VRM"] <- "vrm"
data.n[data.n == "SRI"] <- "sri"
data.n[data.n == "SWE"] <- "swe"
data.n[data.n == "d_roads"] <- "drd"
data.n[data.n == "d_trails"] <- "dtr"
data.n[data.n == "d_heli"] <- "dheli"
data.n[data.n == "CrownClosure"] <- "cc"
data.n[data.n == "IRG"] <- "irg"
data.n[data.n == "NDVI"] <- "ndvi"

names(data) <- data.n #rename columns
data$cc <- data$cc/100 #make percentage
data$LC <- as.factor(data$LC)

#########
#scale and transform variables
data.s <- data[,c("i","j","ij")] #values that cant be scaled
data.s.1 <- data[,c(18:25, 27:32)] #numeric values

data.s.2 <- as.data.frame(scale(data.s.1, scale=T)) #scale numerical values
names(data.s.2) <- paste(names(data.s.2), ".s", sep="")

data.s.3 <- data.s.1^2 #create quadratic variables
#data.s.3 <- data.s.2^2 #create quadratic variables
data.s.3 <- as.data.frame(scale(data.s.3)) #scale after transformation
#names(data.s.3) <- gsub(names(data.s.3), pattern="\\.s",replacement=".q") #unscaled variables
names(data.s.3) <- paste(names(data.s.3), ".q", sep="")

#data.s.4 <- data.s.1[,c("cc","vrm","sri","irg","mgoat")]^2 #scaled variables
#names(data.s.4) <- paste(names(data.s.4), ".q", sep="")

data.s.5 <- log(data.s.1[,!names(data.s.1) %in% c("sri")] + 0.0001) #pseudo threshold
data.s.5 <- as.data.frame(scale(data.s.5)) #scale after transformation
names(data.s.5) <- paste(names(data.s.5), ".p", sep="")

data2 <- cbind(data, data.s.2, data.s.3, data.s.5)

#########
#linear variables
vars <- names(data2) %>% grep(pattern="\\.s$", value=T)
#vars <- c(vars, "cc","vrm","sri","ndvi","irg","mgoat")

vars <- data.frame(s=vars)
vars <- separate(vars, col="s", into=c("var"), remove=F, sep="\\.")

pvars <- data.frame(p=names(data2) %>% grep(pattern="\\.p", value=T))
pvars <- separate(pvars, col="p", into="var", remove=F, sep="\\.")

qvars <- data.frame(q=names(data2) %>% grep(pattern="\\.q", value=T))
qvars <- separate(qvars, col="q", into="var", remove=F, sep="\\.")

v <- merge(vars, pvars, by="var", all.x=T)
v <- merge(v, qvars, by="var", all.x=T)
v$cat <- c("survival","survival","resource","disturbance","resource","disturbance","disturbance","resource","resource","resource","survival",
           "resource","resource","survival")

########
data.ff <- data[data$d.sex == "FF",]
data.mm <- data[data$d.sex == "MM",]
data.mf <- data[data$d.sex == "MF",]

########
#general habitat use
########

data <- data.ff
data <- data.mm
data <- data.mf

family = binomial(link="logit")
null <- glmer(used ~ (1|ij), data=data, family=family) #random effect = ij (dyad)

#female-female model
glob.ff <- glmer(used ~ 
                   e.p+e.s+sri.s+dh2o_per.q+dh2o_per.s+
                   dmlick.p+dmlick.s+irg.q+irg.s+ndvi.s+swe.s+ #check swe effect
                   #e.p*sri.s + e.s*sri.s +
                   vrm.p+vrm.s+desc.s+cc.s+
                   #desc.s*cc.s+
                   #dheli.s+
                   (1|ij),
                 data=data.ff, family=family)
  ss <- getME(glob.ff, c("theta","fixef"))
  glob.ff <- update(glob.ff, start=ss, control=glmerControl(optCtrl=list(maxfun=2e6)))
summary(glob.ff)

lc.ff <- glmer(used ~ factor(LC) + (1|i), data=data2, family=family) 
summary(lc.ff)


#male-male model
glob.mm <- glmer(used ~ e.p+e.s+sri.s+dh2o_per.q+dh2o_per.s+
                   dmlick.p+dmlick.s+irg.q+irg.s+ndvi.s+swe.s+ #check swe effect
                   #e.p*sri.s + e.s*sri.s +
                   vrm.p+vrm.s+desc.s+cc.s+
                   #desc.s*cc.s+
                   #dheli.s+
                   (1|ij),
                 data=data2, family=family)
ss <- getME(glob.mm, c("theta","fixef"))
glob.mm <- update(glob.mm, start=ss, control=glmerControl(optCtrl=list(maxfun=2e6)))
summary(glob.mm)

lc.mm <- glmer(used ~ factor(LC) + (1|i), data=data2, family=family) 
summary(lc.mm)

#male-female model
glob.mf <- glmer(used ~ e.p+e.s+sri.s+dh2o_per.q+dh2o_per.s+
                   dmlick.p+dmlick.s+irg.q+irg.s+ndvi.s+swe.s+ #check swe effect
                   #e.p*sri.s + e.s*sri.s +
                   vrm.p+vrm.s+desc.s+cc.s+
                   #desc.s*cc.s+
                   #dheli.s+
                   (1|ij),
                 data=data2, family=family)
ss <- getME(glob.mf, c("theta","fixef"))
glob.mf <- update(glob.mf, start=ss, control=glmerControl(optCtrl=list(maxfun=2e6)))
summary(glob.mf)

lc.mf <- glmer(used ~ factor(LC) + (1|i), data=data2, family=family) 
summary(lc.mf)

##############
#combine coefficients from models and plot

globalbeta.mf <- as.data.frame(coefficients(summary(glob.mf)))
globalbeta.mf$sex <- "MF"
globalbeta.mf$variable <- row.names(globalbeta.mf)

globalbeta.mm <- as.data.frame(coefficients(summary(glob.mm)))
globalbeta.mm$sex <- "MM"
globalbeta.mm$variable <- row.names(globalbeta.mm)

globalbeta.ff <- as.data.frame(coefficients(summary(glob.ff)))
globalbeta.ff$sex <- "FF"
globalbeta.ff$variable <- row.names(globalbeta.ff)

globalbeta <- rbind(globalbeta.mf, globalbeta.mm, globalbeta.ff)

globalbeta <- separate(globalbeta, col="variable", into=c("var","type"), sep="\\.", remove=F)
globalbeta[is.na(globalbeta$type),]$type <- "s"
#globalbeta[grep(globalbeta$variable, pattern=":"),]$var <- "interaction"
#globalbeta[grep(globalbeta$variable, pattern=":"),]$type <- "i"
globalbeta$var2 <- gsub(globalbeta$variable, pattern="\\.s", replacement="")
names(globalbeta) <- c("Estimate","StdError","zvalue","Pr(>|z|)","sex","variable","var","type","var2")
globalbeta$modeltype <- "global"
globalbeta <- merge(globalbeta, v[,c("var","cat")], by="var", all.x=T)
globalbeta[is.na(globalbeta$cat),]$cat <- "intercept"
#globalbeta[is.na(globalbeta$cat),]$cat <- c("intercept","resource","survival","disturbance","resource","survival","disturbance","resource","survival","disturbance")

betaorder <- unique(globalbeta[,c("cat","var2","type")])
#betaorder$cat <- factor(betaorder$cat, levels=c("resource","survival","disturbance"))

# globalbeta$var2 <- factor(globalbeta$var2, levels=c("sri:ndvi","dh2o_per","dmlick.p","dmlick","e.p","e","irg.q","irg","ndvi","sri",
#                                                     "desc:cc","cc","desc.p","desc","vrm",
#                                                     "dtr:dheli","dheli","drd","dtr")) 
globalbeta$var2 <- factor(globalbeta$var2, levels=c("dh2o_per.q","dh2o_per","dmlick.p","dmlick","e.p","e","irg.q","irg","ndvi","sri","swe",
                                                    "cc","desc","vrm.p","vrm",
                                                    "dheli",
                                                    "(Intercept)"))

g <- ggplot(data=globalbeta, aes(x=Estimate, y=var2, col=sex)) + 
  geom_vline(aes(xintercept=0), lty="dashed") + 
  geom_point(aes(shape=type), size=3, position=position_dodge2(width=0.5)) + 
  geom_errorbar(aes(xmin=Estimate-StdError, xmax=Estimate+StdError), width=0.5, position=position_dodge2(width=0.5)) +
  scale_shape_manual(values=c(18,16,17,15)) +
  scale_color_manual(values=c("pink","black","lightblue")) +
  ylab("variable") + xlab("beta estimate") + 
  theme_bw(base_size = 15)
g
ggsave(g, filename="C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper/coefficients_global_contacts_2020-11-17.tiff",
       compression="lzw", height=10, width=12, units="in", dpi=300)

r <- ggplot(data=globalbeta[globalbeta$cat == "resource",], aes(x=Estimate, y=var2, col=sex)) + 
  geom_vline(aes(xintercept=0), lty="dashed") + 
  geom_point(aes(shape=type), size=3, position=position_dodge2(width=0.5)) + 
  geom_errorbar(aes(xmin=Estimate-StdError, xmax=Estimate+StdError), width=0.5, position=position_dodge2(width=0.5)) +
  scale_shape_manual(values=c(18,16,17,15)) +
  scale_color_manual(values=c("pink","black","lightblue")) +
  ylab("variable") + xlab("beta estimate") + 
  theme_bw(base_size = 15)
s <- ggplot(data=globalbeta[globalbeta$cat == "survival",], aes(x=Estimate, y=var2, col=sex)) + 
  geom_vline(aes(xintercept=0), lty="dashed") + 
  geom_point(aes(shape=type), size=3, position=position_dodge2(width=0.5)) + 
  geom_errorbar(aes(xmin=Estimate-StdError, xmax=Estimate+StdError), width=0.5, position=position_dodge2(width=0.5)) +
  scale_shape_manual(values=c(18,16,17,15)) +
  scale_color_manual(values=c("pink","black","lightblue")) +
  ylab("variable") + xlab("beta estimate") + 
  theme_bw(base_size = 15)
# d <- ggplot(data=globalbeta[globalbeta$cat == "disturbance",], aes(x=Estimate, y=var2, col=sex)) + 
#   geom_vline(aes(xintercept=0), lty="dashed") + 
#   geom_point(aes(shape=type), size=3, position=position_dodge2(width=0.5)) + 
#   geom_errorbar(aes(xmin=Estimate-StdError, xmax=Estimate+StdError), width=0.5, position=position_dodge2(width=0.5)) +
#   scale_shape_manual(values=c(18,16,17,15)) +
#   scale_color_manual(values=c("pink","black","lightblue")) +
#   ylab("variable") + xlab("beta estimate") + 
#   theme_bw(base_size = 15)


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

###########################
#4d. comparison between general habitat and direct contact probabilities
###########################

require(raster)
require(ggplot2)
require(dplyr)
require(viridis)
require(RColorBrewer)
library(rasterVis)
require(ggpubr)
require(gridExtra)
require(rgdal)

path.local <- "C:/Users/tosam/Documents/1_Post-Master's_2015-present/2016 USGS GNP bighorn sheep/social contacts paper"
setwd(path.local)

boundaries <- readOGR("CCE_Boundary_2008_UTM11NAD83", dsn="CCE Boundary")
waterton <- readOGR(dsn="Park Boundaries", "Waterton Lakes National Park")
boundaries <- readOGR("Waterton_Glacier_PeacePark", dsn="Park Boundaries")
tribal <- readOGR(".", "Indian_Reservations")
tribal <- spTransform(tribal, projection(boundaries))

blackfeet <- tribal[tribal$IND_NAME %in% "Blackfeet Reservation",]

#########################
########## not classes, raw probabilities
#SSF rasters (general use), not classes
f.ssf <- raster("model_rasters/ssf_f.tif")
m.ssf <- raster("model_rasters/ssf_m.tif")
mf.ssf <- raster("model_rasters/ssf_mf.tif")

#contacts RSF, not classes
ff.rsf <- raster("model_rasters/contact_rsf_ff.tif")
mm.rsf <- raster("model_rasters/contact_rsf_mm.tif")
mf.rsf <- raster("model_rasters/contact_rsf_mf.tif")

#calculate Pr(contact) = Pr(general use) * Pr(contact|general use)
f <- f.ssf*ff.rsf
m <- m.ssf*mm.rsf
mf <- mf.ssf*mf.rsf

f <- f/max(values(f), na.rm=T) #rescale from 0-1
m <- m/max(values(m), na.rm=T)
mf <- mf/max(values(mf), na.rm=T)

all <- stack(f, m, mf)
all.plot <- levelplot(all, layout=c(3,1),
                      colorkey=list(space='bottom', axis.line=list(col='black'), width=1),
                      #par.settings=list(axis.line=list(col='transparent')), #remove axes
                      #par.settings=list(strip.border=list(col='transparent'), strip.background=list(col='transparent'), axis.line=list(col='transparent')),
                      #par.strip.text=list(cex=1.5, lines=2, fontface='bold'),
                      scales=list(draw=F), #remove values of x and y axes
                      col.regions=viridis, names.attr=rep("", 3)) + layer(sp.polygons(waterton,lwd=2)) + layer(sp.polygons(boundaries,lwd=2)) + layer(sp.polygons(blackfeet, lwd=2))
all.plot

######################
########## classes
#SSF rasters, deciles
f.ssf <- raster("model_rasters/ssf_class_f.tif")
m.ssf <- raster("model_rasters/ssf_class_m.tif")
mf.ssf <- raster("model_rasters/ssf_class_mf.tif")

#contacts RSF rasters
ff.rsf <- raster("model_rasters/contact_rsf_class_ff.tif")
mm.rsf <- raster("model_rasters/contact_rsf_class_mm.tif")
mf.rsf <- raster("model_rasters/contact_rsf_class_mf.tif")

ssf <- stack(f.ssf, m.ssf, mf.ssf)
rsf <- stack(ff.rsf, mm.rsf, mf.rsf)

comp.stack <- stack(f.ssf, m.ssf, mf.ssf, ff.rsf, mm.rsf, mf.rsf)
comp.vals <- data.frame(values(comp.stack))
names(comp.vals) <- c("ssf_class_f","ssf_class_m","ssf_class_mf","contact_rsf_class_ff","contact_rsf_class_mm","contact_rsf_class_mf")
comp.vals <- comp.vals[!is.na(comp.vals$ssf_class_f),]
comp.vals <- comp.vals[!is.na(comp.vals$ssf_class_m),]
comp.vals <- comp.vals[!is.na(comp.vals$ssf_class_mf),]

########
#spearman rank test
#ssf
cor.test(comp.vals$ssf_class_f, comp.vals$ssf_class_m, method="spearman")

# Spearman's rank correlation rho
# data:  comp.vals$ssf_class_f and comp.vals$ssf_class_m
# S = 5.0873e+15, p-value < 2.2e-16
# alternative hypothesis: true rho is not equal to 0
# sample estimates:
#       rho 
# 0.9853978 

#ssf vs. contact rsf
cor.test(comp.vals$ssf_class_f, comp.vals$contact_rsf_class_ff, method="spearman")
# data:  comp.vals$ssf_class_f and comp.vals$contact_rsf_class_ff
# S = 5.6233e+17, p-value < 2.2e-16
# alternative hypothesis: true rho is not equal to 0
# sample estimates:
#   rho 
# -0.6140703 
cor.test(comp.vals$ssf_class_m, comp.vals$contact_rsf_class_mm, method="spearman")
# data:  comp.vals$ssf_class_m and comp.vals$contact_rsf_class_mm
# S = 6.1415e+17, p-value < 2.2e-16
# alternative hypothesis: true rho is not equal to 0
# sample estimates:
#   rho 
# -0.7627884 
cor.test(comp.vals$ssf_class_mf, comp.vals$contact_rsf_class_mf, method="spearman")
# data:  comp.vals$ssf_class_mf and comp.vals$contact_rsf_class_mf
# S = 5.4448e+17, p-value < 2.2e-16
# alternative hypothesis: true rho is not equal to 0
# sample estimates:
#   rho 
# -0.5628252 

###
#contact vs. contact
###
cor.test(comp.vals$contact_rsf_class_mm, comp.vals$contact_rsf_class_ff, method="spearman")
# data:  comp.vals$contact_rsf_class_mm and comp.vals$contact_rsf_class_ff
# S = 1.186e+17, p-value < 2.2e-16
# alternative hypothesis: true rho is not equal to 0
# sample estimates:
#   rho 
# 0.6595887 

cor.test(comp.vals$contact_rsf_class_mm, comp.vals$contact_rsf_class_mf, method="spearman")
# data:  comp.vals$contact_rsf_class_mm and comp.vals$contact_rsf_class_mf
# S = 1.5692e+17, p-value < 2.2e-16
# alternative hypothesis: true rho is not equal to 0
# sample estimates:
#   rho 
# 0.5495995 

cor.test(comp.vals$contact_rsf_class_ff, comp.vals$contact_rsf_class_mf, method="spearman")
# data:  comp.vals$contact_rsf_class_ff and comp.vals$contact_rsf_class_mf
# S = 4.0932e+16, p-value < 2.2e-16
# alternative hypothesis: true rho is not equal to 0
# sample estimates:
#   rho 
# 0.8825114 

########
#graph figures
comp.vals <- mutate_all(comp.vals, as.factor)

#summarize data
table(comp.vals[,c("ssf_class_f","contact_rsf_class_ff")])
table(comp.vals[,c("ssf_class_m","contact_rsf_class_mm")])
table(comp.vals[,c("ssf_class_mf","contact_rsf_class_mf")])

#make bar chart comparing ssf to contact.rsf
f <- ggplot(comp.vals[!is.na(comp.vals$contact_rsf_class_ff),]) +
  geom_bar(position="fill", aes(x=ssf_class_f, group=contact_rsf_class_ff, fill=contact_rsf_class_ff), stat="count") +
  scale_fill_manual(values=brewer.pal(10, "PRGn")) +
  #scale_fill_manual(values=brewer.pal(10, "PuOr")) +
  xlab("Habitat rank") + ylab("Proportion of habitat rank") +
  labs(fill="Contact rank") +
  theme_bw()
m <- ggplot(comp.vals[!is.na(comp.vals$contact_rsf_class_mm),]) +
  geom_bar(position="fill", aes(x=ssf_class_m, group=contact_rsf_class_mm, fill=contact_rsf_class_mm), stat="count") +
  scale_fill_manual(values=brewer.pal(10, "PRGn")) +
  xlab("Habitat rank") + ylab("Proportion of habitat rank") +
  labs(fill="Contact rank") +
  theme_bw()
mf <- ggplot(comp.vals[!is.na(comp.vals$contact_rsf_class_mf),]) +
  geom_bar(position="fill", aes(x=ssf_class_mf, group=contact_rsf_class_mf, fill=contact_rsf_class_mf), stat="count") +
  scale_fill_manual(values=brewer.pal(10, "PRGn")) +
  xlab("Habitat rank") + ylab("Proportion of habitat rank") +
  labs(fill="Contact rank") +
  theme_bw()