# Traverse angle versus landscape slope: a worked analysis of the table that
# Split Lines into Topography Segments writes, run in R.
#
# Make the segments with Segmentation = "At the vertices only", so that each
# GPS step is one segment, and choose these attributes: Relative Slope,
# Bearing, Average Landscape Slope, Average Landscape Compass Direction,
# Traverse Angle, Planimetric Length and Source Line ID, plus a copied
# animal-ID field. For the step-selection model in step 4, also open the
# tool's "Advanced: alternative steps" group and ask for alternative steps
# (10 or 35 per real step): the output then carries every real step and its
# alternatives together, marked by Seg_Used (1 = real, 0 = alternative) and
# tied by Seg_Step_ID. Export the attribute table (Table To Table, or Export
# Table) as segments.csv, then set the names below to match your fields.
#
# The block marked SIMULATED builds a made-up table of the same shape so the
# script runs on its own; delete it and read your own file instead.

library(survival)   # clogit: conditional logistic regression, ships with R
library(lme4)       # lmer: mixed models (install.packages("lme4"))

## ---- SIMULATED data, stand-in for read.csv("segments.csv") -----------------
## Four animals, 300 real steps each, and 10 alternative steps per real step.
## (The simulation prices every alternative on the same hillside as its real
## step; the tool measures each alternative on the ground it actually
## crosses, which is the point of casting them in GIS.)
set.seed(1)
animals <- c("F01", "F02", "M01", "M02")
K <- 10
d_all <- do.call(rbind, lapply(seq_along(animals), function(ai) {
  a <- animals[ai]; n <- 300
  land_slope <- runif(n, 0, 40)                       # terrain slope, degrees
  k <- c(F01 = 8, F02 = 12, M01 = 20, M02 = 30)[a]     # how quickly each animal
  trav <- 90 * (1 - exp(-land_slope / k)) +           #   turns to the contour
          rnorm(n, 0, 12)
  trav <- pmin(pmax(trav, 0), 90)
  path <- atan(tan(land_slope * pi / 180) * cos(trav * pi / 180)) * 180 / pi
  step_id <- (ai - 1) * n + seq_len(n)
  real <- data.frame(AnimalID = a, Seg_Source_Line_ID = ai,
                     Seg_Used = 1L, Seg_Step_ID = step_id, Seg_Alt_Number = 0L,
                     Seg_Average_Landscape_Slope = land_slope,
                     Seg_Traverse_Angle = trav,
                     Seg_Relative_Slope = path * sample(c(-1, 1), n, TRUE),
                     Seg_Planimetric_Length = exp(rnorm(n, log(150), 0.5)))
  alt <- real[rep(seq_len(n), each = K), ]
  alt$Seg_Used <- 0L
  alt$Seg_Alt_Number <- rep(seq_len(K), n)
  alt$Seg_Traverse_Angle <- runif(nrow(alt), 0, 90)
  alt$Seg_Relative_Slope <- atan(tan(alt$Seg_Average_Landscape_Slope * pi / 180) *
                                 cos(alt$Seg_Traverse_Angle * pi / 180)) * 180 / pi *
                            sample(c(-1, 1), nrow(alt), TRUE)
  rbind(real, alt)
}))
## ---------------------------------------------------------------------------
# Replace "segments.csv" below with the path to your CSV file, then
# uncomment the line (and delete the SIMULATED block above) to run the
# script on your own data.
# d_all <- read.csv("segments.csv")

slope_col <- "Seg_Average_Landscape_Slope"
trav_col  <- "Seg_Traverse_Angle"
path_col  <- "Seg_Relative_Slope"
len_col   <- "Seg_Planimetric_Length"
# The grouping variable: the animal. Use the animal-ID field you copied
# from the source table. If each line IS one animal and you copied no ID
# field, use "Seg_Source_Line_ID" instead; the tool always writes it. If an
# animal has several lines (sessions, seasons), keep AnimalID here and see
# the note at step 3.
group_col <- "AnimalID"   # <-- Change 'AnimalID' here to the appropriate field name

# Real steps only for steps 1-3 (the alternatives, if any, are for step 4).
has_alt <- "Seg_Used" %in% names(d_all)
d <- if (has_alt) d_all[d_all$Seg_Used == 1, ] else d_all
d$group <- factor(d[[group_col]])

# 1. Aspect, and so the traverse angle, means little on near-flat ground:
#    keep steps where the terrain slope is at least 5 degrees.
d <- d[!is.na(d[[trav_col]]) & d[[slope_col]] >= 5, ]
d$path_angle <- abs(d[[path_col]])

# 2. The fall-line fidelity index: path angle regressed on terrain slope
#    through the origin, per animal. 1 = takes the hill as it comes,
#    near 0 = contours. It is run twice, UNWEIGHTED (2a) and WEIGHTED by
#    step length (2b); keep whichever fits your question and delete the
#    other. The weighted version is the unweighted one with the
#    "weights = ..." argument added, and nothing else.
#
#    Unweighted, every step is one vote, whatever its length: the right
#    choice when each step is one decision, as with GPS fixes at a fixed
#    time interval, where a longer step only means the animal moved faster.
#    Weight by planimetric length instead when the question is about
#    distance travelled ("what share of the route follows the contour?")
#    or when the vertices are unevenly spaced, as on a digitized trail.
#    Planimetric length is the better weight: surface length grows with
#    the path angle, so it would favor steep steps.

# 2a. UNWEIGHTED
fidelity <- sapply(split(d, d$group), function(x)
  coef(lm(path_angle ~ 0 + Seg_Average_Landscape_Slope, data = x)))
print(round(fidelity, 3))

# 2b. WEIGHTED by planimetric length (delete this block for an unweighted
#     analysis, or delete the "weights = x[[len_col]]" argument, which
#     turns it back into 2a)
fidelity_w <- sapply(split(d, d$group), function(x)
  coef(lm(path_angle ~ 0 + Seg_Average_Landscape_Slope, data = x,
          weights = x[[len_col]])))
print(round(fidelity_w, 3))

# 3. Does the traverse angle rise with terrain slope? A mixed model with the
#    animal as a random effect, so that the many steps of one animal do not
#    count as independent evidence. Again run UNWEIGHTED (3a) and WEIGHTED
#    (3b); keep one.
#    With several lines per animal (sessions, seasons), add the line as a
#    second grouping nested in the animal, so steps within one line are not
#    treated as independent either:
#      (1 | group) + (1 | Seg_Source_Line_ID)

# 3a. UNWEIGHTED
m <- lmer(Seg_Traverse_Angle ~ Seg_Average_Landscape_Slope + (1 | group),
          data = d)
print(summary(m)$coefficients)

# 3b. WEIGHTED by step length: the same call with a "weights =" argument
#     (lmer treats weights as relative precisions, so use them for the
#     distance-travelled question and read the standard errors with that
#     in mind). Delete this block for an unweighted analysis.
m_w <- lmer(Seg_Traverse_Angle ~ Seg_Average_Landscape_Slope + (1 | group),
            data = d, weights = d[[len_col]] / mean(d[[len_col]]))
print(summary(m_w)$coefficients)

# 4. A step-selection function. Every real step is compared with the
#    alternative steps the tool cast from the same start point (Seg_Used
#    0), each measured on the ground it crosses, and conditional logistic
#    regression, stratified by Seg_Step_ID so that a step is only ever
#    compared with its own alternatives, asks whether steep path angles
#    were avoided, and whether that avoidance strengthens on steeper
#    terrain. Needs the tool's "Alternative steps per real step" set.
#    The cluster() term gives standard errors that allow for the steps of
#    one animal (or line) not being independent of each other; it needs
#    method = "efron", because clogit's default exact method cannot
#    compute a robust variance.
if (has_alt) {
  ssf <- d_all[!is.na(d_all[[path_col]]) & !is.na(d_all[[slope_col]]), ]
  ssf$path_angle <- abs(ssf[[path_col]])
  ssf$land_slope <- ssf[[slope_col]]
  ssf$cluster_id <- ssf[[group_col]]
  # keep only strata that still have the real step AND at least one alternative
  ok <- tapply(ssf$Seg_Used, ssf$Seg_Step_ID, function(u) any(u == 1) && any(u == 0))
  ssf <- ssf[ssf$Seg_Step_ID %in% names(ok)[ok], ]
  fit <- clogit(Seg_Used ~ path_angle + path_angle:land_slope + strata(Seg_Step_ID)
                + cluster(cluster_id), data = ssf, method = "efron")
  print(summary(fit)$coefficients)
  # A negative path_angle coefficient: steep steps are avoided. A negative
  # interaction: the avoidance strengthens as the hillside steepens. Read
  # the "robust se" column.
} else {
  cat("No Seg_Used column: run the tool with alternative steps for step 4.\n")
}

# 5. Figures, in base R (nothing to install). Set out_dir to a folder; four
#    PNG files are written there. They follow the analysis above: 5a and 5b
#    need steps 2-3, 5c and 5d need the alternative steps of step 4.
out_dir <- "."

# 5a. Traverse angle against terrain slope: every real step, the mean in
#     each 5-degree band of slope with its 95% confidence interval, and the
#     mixed-model line.
png(file.path(out_dir, "fig_traverse_vs_slope.png"), width = 1600, height = 1100, res = 200)
plot(d[[slope_col]], d[[trav_col]], pch = 16, cex = 0.5, col = rgb(0, 0, 0, 0.25),
     xlab = "Terrain slope under the step (degrees)", ylab = "Traverse angle (degrees)",
     main = "Real steps: traverse angle against terrain slope")
b <- fixef(m)
abline(a = b[1], b = b[2], lwd = 2, col = "firebrick")
bins <- cut(d[[slope_col]], breaks = seq(0, ceiling(max(d[[slope_col]]) / 5) * 5, by = 5))
mids <- tapply(d[[slope_col]], bins, mean)
mn   <- tapply(d[[trav_col]], bins, mean)
se   <- tapply(d[[trav_col]], bins, function(v) sd(v) / sqrt(length(v)))
ok5  <- !is.na(mn) & tapply(d[[trav_col]], bins, length) >= 5
arrows(mids[ok5], (mn - 1.96 * se)[ok5], mids[ok5], (mn + 1.96 * se)[ok5],
       angle = 90, code = 3, length = 0.03, col = "navy")
points(mids[ok5], mn[ok5], pch = 15, col = "navy")
legend("topleft", c("one real step", "mean in a 5-degree band (95% CI)", "mixed-model line"),
       pch = c(16, 15, NA), lty = c(NA, NA, 1), col = c("grey40", "navy", "firebrick"), bty = "n")
dev.off()

# 5b. The fall-line fidelity index by group, with the pooled value.
pooled <- coef(lm(path_angle ~ 0 + Seg_Average_Landscape_Slope, data = d))
f <- sort(fidelity)
names(f) <- sub("[.]Seg_Average_Landscape_Slope$", "", names(f))
png(file.path(out_dir, "fig_fidelity_by_group.png"), width = 1400, height = 1300, res = 200)
dotchart(f, pch = 16, xlim = c(0, 1), xlab = "Fall-line fidelity index (path angle per degree of terrain slope)",
         main = "Fidelity index by group")
abline(v = pooled, lty = 2, col = "firebrick")
abline(v = c(0, 1), col = "grey75")
mtext("0 = follows the contour; 1 = takes the hill as it comes; dashed = all groups pooled", side = 1, line = 4, cex = 0.8)
dev.off()

if (has_alt) {
  # 5c. Where the real step ranks among its alternatives: how many of them
  #     were steeper than the step the animal took. With no preference the
  #     bars would be level (dashed line); a pile-up at the left is selection.
  steeper <- tapply(seq_len(nrow(ssf)), ssf$Seg_Step_ID, function(idx) {
    u <- ssf$path_angle[idx][ssf$Seg_Used[idx] == 1][1]
    sum(ssf$path_angle[idx][ssf$Seg_Used[idx] == 0] > u)
  })
  K_alt <- max(tapply(ssf$Seg_Used == 0, ssf$Seg_Step_ID, sum))
  png(file.path(out_dir, "fig_rank_among_alternatives.png"), width = 1600, height = 1100, res = 200)
  barplot(table(factor(steeper, levels = 0:K_alt)), col = "steelblue",
          xlab = "Number of the step's alternatives that were steeper than the step taken",
          ylab = "Real steps", main = "The step taken, ranked against its alternatives")
  abline(h = length(steeper) / (K_alt + 1), lty = 2)
  dev.off()

  # 5d. The model's own curve: the odds of a candidate step being the one
  #     taken, relative to a level step, against its path angle, on a level
  #     start and on a 30-degree hillside (the interaction).
  bb <- coef(fit); pa <- seq(0, 40, by = 0.5)
  png(file.path(out_dir, "fig_relative_odds.png"), width = 1600, height = 1100, res = 200)
  plot(pa, exp(bb[1] * pa), type = "l", lwd = 2, col = "firebrick", ylim = c(0, 1),
       xlab = "Path angle of a candidate step (degrees)",
       ylab = "Odds of being the step taken, relative to a level step",
       main = "Step-selection model: the price of a steep step")
  lines(pa, exp(bb[1] * pa + bb[2] * pa * 30), lwd = 2, lty = 2, col = "navy")
  legend("topright", c("on a level start (terrain slope 0)", "on a 30-degree hillside"),
         lty = c(1, 2), col = c("firebrick", "navy"), bty = "n")
  dev.off()
}