Introduction

Over the past couple of years, students have been trying to implement SAMS (a type of ecological memory modeling) to find weights in the day to day decisions of birds. Do a few specific days of migration affect their reproductive success more? Is it more important for them to eat more during migration? Right before migration? Mitch has hoped to realize these questions from the SAM. But, all analyses done by students have shown the same thing…the weights are all the same, across the whole migration path.

Is this because ecologically the weights aren’t important? Meaning every day is contributing the same? Do we have too long of a time series? Is the within day time series greater than variation among individuals? I’m hoping to answer these questions by creating simulations and following Ogle et al., 2015 (as her data produced weights).

There’s three different ways in which we can approach this analysis.

Utilize Ogle’s data and figure out why her data works, and ours doesn’t
Simulate PTF and FTI success using a variety of different weights that we CHOOSE and see if the model will output those successfully

Meaning, if we simulate data with weights, will the model show us weights?

In all honesty, I did not keep very good track of my simulations and such, as I more wanted to see if SAMS would work or not…not try and figure out why they weren’t working. Once we kind of determined that it was most likely due to the Bernoulli response variable, I stopped going. I’m happy to keep working on this as needed, but this markdown file has the code needed to contiune figuring out why the SAMS may not be working.

Our data

Before I get into the simulations, I want to first show you what our data looks like. Here is a practice data set of 49 individuals.

METADATA

id: indexing for individual
year2: indexing for year
birdid: band number associated with each bird
year: year of the bird’s life
FTI: Full term incubation (1 meaning they were successful, 0 meaning they were not)
rel_day: relative day in the bird’s migration. For ease of getting this data together, dates for each bird were just filtered so that each bird had a max of 39 days. If a bird was banded that year, they may have less. In the real analysis, only the days that the bird is actively migrating will be considered. Birds have between 28 and 39 days.
max_day: max relative day for each bird
scaled_ODBA: Overall dynamic body acceleration scaled
scaled_PTF: Proportion of time feeding scaled

We can see here that 32 birds failed whereas 17 were successful.

Here is what the SAM model looks like. This model was built off of by many members of the lab, most recently Alec Schindler. I changed it to look at how do the antecedent effects of proportion of time feeding and ODBA (movement) affect the reproductive success of individuals. In the future, I hope to add additional covariates such as age or arrival date on breeding grounds to explain this, but for now, I’ve kept it minimal.

Additionally, I’m trying to code in NIMBLE now, but since last year I ran it in JAGS, I’ve kept it in that code.

# specify model
sink("SAM_bern.txt")
cat("
    model {
    
    ### Priors

    # intercept
    alpha ~ dnorm(0, 0.01)
      
    # betas
    beta_PTF ~ dnorm(0, 0.01)
    beta_ODBA ~ dnorm(0, 0.01)

    # random year effect
    for(t in 1:nyears){
      eps_year[t] ~ dnorm(0, tau_year)
    }
    tau_year ~ dgamma(0.01, 0.01)
    sd_year <- 1/sqrt(tau_year)
    
    for(k in 1:nobs1) {
      # Dirichlet prior for daily PTF weights
      delta_PTF[ids1[k],years1[k],days1[k]] ~ dgamma(1,1)
      # Dirichlet prior for daily ODBA weights
      delta_ODBA[ids1[k],years1[k],days1[k]] ~ dgamma(1,1) # (rate, shape)
      
      # normalize PTF weights to 1
      weight_PTF[ids1[k],years1[k],days1[k]] <- delta_PTF[ids1[k],years1[k],days1[k]]/sum_D1[ids1[k],years1[k]]
      # normalize ODBA weights to 1
      weight_ODBA[ids1[k],years1[k],days1[k]] <- delta_ODBA[ids1[k],years1[k],days1[k]]/sum_D2[ids1[k],years1[k]]
      
      # For each time period into the past, compute weighted PTF variable
      weighted_PTF[ids1[k],years1[k],days1[k]] <- weight_PTF[ids1[k],years1[k],days1[k]]*PTF[ids1[k],years1[k],max_days1[k]-days1[k]+1]
      # For each time period into the past, compute weighted ODBA variable
      weighted_ODBA[ids1[k],years1[k],days1[k]] <- weight_ODBA[ids1[k],years1[k],days1[k]]*ODBA[ids1[k],years1[k],max_days1[k]-days1[k]+1]
      
      # Reorder weights from recent to older
      weight_ordered_PTF[ids1[k],years1[k],max_days1[k]-days1[k]+1] <- weight_PTF[ids1[k],years1[k],days1[k]]
      weight_ordered_ODBA[ids1[k],years1[k],max_days1[k]-days1[k]+1] <- weight_ODBA[ids1[k],years1[k],days1[k]]
      
      # Compute cumulative daily weights
      cum_weight_PTF[ids1[k],years1[k],days1[k]] <- sum(weight_ordered_PTF[ids1[k],years1[k],1:days1[k]])
      cum_weight_ODBA[ids1[k],years1[k],days1[k]] <- sum(weight_ordered_ODBA[ids1[k],years1[k],1:days1[k]])
    }
    
    for(l in 1:nobs2){
      # Compute sum of deltas
      sum_D1[ids2[l],years2[l]] <- sum(delta_PTF[ids2[l],years2[l],1:max_days2[l]])
      sum_D2[ids2[l],years2[l]] <- sum(delta_ODBA[ids2[l],years2[l],1:max_days2[l]])
      
      # calculate antecedent PTF
      ant_PTF[ids2[l],years2[l]] <- sum(weighted_PTF[ids2[l],years2[l],1:max_days2[l]])
      # calculate antecedent ODBA
      ant_ODBA[ids2[l],years2[l]] <- sum(weighted_ODBA[ids2[l],years2[l],1:max_days2[l]])
    
    ### Likelihood
      BS[l] ~ dbern(p[ids2[l],years2[l]])
      logit(p[ids2[l],years2[l]]) <- alpha + beta_PTF * ant_PTF[ids2[l],years2[l]] + beta_ODBA * ant_ODBA[ids2[l],years2[l]] +  eps_year[years2[l]]
    }
    
    }", fill=TRUE)
sink()

# Prepare the data list for SAM
jags.data_SAM <- list(
  nyears = length(unique(totaldataforSAM$year2)),
  nobs1 = nrow(totaldataforSAM),
  ids1 = totaldataforSAM$id,
  years1 = totaldataforSAM$year2,
  days1 = totaldataforSAM$rel_day,
  max_days1 = totaldataforSAM$max_day,
  PTF = PTF.data,# using the predicted values from the DLM
  ODBA = ODBA.data,
  nobs2 = nrow(totaldataforSAM.dat2),
  ids2 = totaldataforSAM.dat2$id,
  years2 = totaldataforSAM.dat2$year2,
  max_days2 = totaldataforSAM.dat2$max_day,
  BS = totaldataforSAM.dat2$FTI
)


# Initial values for SAM
inits_SAM <- function() {
  list(
    alpha = rnorm(1),
    beta_PTF = rnorm(1),
    beta_ODBA = rnorm(1),
    tau_year = rgamma(1, 0.01, 0.01)
  )
}

# Parameters monitored in SAM
params_SAM <- c("alpha", "beta_PTF", "beta_ODBA", "p", "eps_year", "tau_year", "sd_year",
                "weight_PTF", "weight_ODBA", "cum_weight_PTF", "cum_weight_ODBA")

# MCMC settings
ni <- 10000
nb <- 2000
nt <- 10
nc <- 3

# Run SAM model
library(jagsUI)
# out_SAM <- jags(jags.data_SAM, inits_SAM, params_SAM, "SAM_bern.txt", n.thin=nt, n.chains=nc, n.burnin=nb, n.iter=ni, parallel=T)

I ran this model last summer yet I’m not sure what I saved the data as. I was more interested in getting the code to run, not the results. Instead, I ran it again now. Here’s a head of the results output. When I ran this last year, I’m almost positive beta_ODBA didn’t overlap zero but it is now. Weights may not show because our parameters aren’t “significant”, but I don’t think that’s the only reason why.

To show the current issue with our weights, I’m going to present some of the individuals and their weights. You can see here that the confidence intervals are quite large and no weight deviates from zero. This is what all members of the lab have experienced.

Looking at Ogle’s data

To determine what the issue was, I brought in Ogle’s data and replicated her results. Please have Ogle et al., 2015 to reference.

Here is the data:

### data

# dataset #1

block = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
          18, 19, 20, 21, 22, 23, 24, 25, 25, 26, 26, 27, 27, 28, 28,
          29, 29, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 34,
          34, 35, 35, 35, 36, 36, 36, 37, 37, 37, 38,
          38, 38)
block <- matrix(block, 5, 12)


N <- 52
Nlag <- 5
Nyrs <- 91
Nblocks <- 38

# dataset #2
dataset2 <- read.csv("test_data.csv") %>%
  mutate(Event1 = as.numeric(Event1),
         Event2 = as.numeric(Event2),
         Event3 = as.numeric(Event3),
         Event4 = as.numeric(Event4),
         NPP = as.numeric(NPP))

Event <- dataset2[c("Event1", "Event2", "Event3", "Event4")]
YearID <- dataset2$YearID
NPP <- dataset2$NPP


# dataset #3
dataset3 <- read.csv("test_data2.csv")

ppt <- dataset3[c("ppt1", "ppt2", "ppt3", "ppt4", "ppt5", "ppt6", "ppt7", "ppt8",
             "ppt9", "ppt10", "ppt11", "ppt12")]


numeric_vector <- as.numeric(as.matrix(ppt))

Here is the model

# sink("OgleSAM.txt")
# cat("
# 
#       model{
#     for(i in 1:N){
#       # Data model (or likelihood) for the observed NPP data:
#       NPP[i] ~ dnorm(mu[i], tau)
#       # Generate “replicated data” to evaluate model fit.
#       NPP.rep[i] ~ dnorm(mu[i],tau)
#       # Define model for latent (mean) NPP; Event[,k] represents the amount
#       # of precipitation received in different size classes, where k indexes
#       # the even size class (k=1 for < 5 mm; k=2 for 5-15 mm; k=3 for 15-
#       # 30 mm; k=4 for >30 mm); convert antecedent precipitation (antX) from
#       # inches to mm.
#       mu[i] <- a[1] + a[2]*antX[YearID[i]]*25.4 + a[3]*Event[i,1] +
#         a[4]*Event[i,2] + a[5]*Event[i,3] + a[6]*Event[i,4]
#       # Some of the precipitation event data are missing, so specify a simple
#       # data model for the Event data for the purpose of estimating the
#       # missing data:
#       for(k in 1:4){
#         Event[i,k] ~ dnorm(mu.ev[k], tau.ev[k])
#     }
#   }
#   # Dirichlet prior for monthly precipitation weights (due to restrictions
#   # on when the built-in dirichlet distribution can be used, we are required
#   # to use the relationship between the gamma distribution and the dirichlet
#   # to assign the dirichlet prior. For each time block into the past, assign
#   # the unnormalized weight (deltaX) a gamma(1,1) prior:
#   for(j in 1:Nblocks){
#     deltaX[j] ~ dgamma(1,1)
#   }
#   for(t in 1:Nlag){
#     # Compute the yearly weights:
#     yr.w[t] <- sum(weight[,t])
#     alphad[t] <- 1
#     for(m in 1:12){
#       # Redefine the unnormalized monthly weights to account for post-ANPP # harvest period; i.e., 2nd part involving equals and step functions # sets weight = 0 if in year 1 and in Oct, Nov, or Dec (i.e., post- # ANPP harvest).
#       delta[m,t] <- (deltaX[block[t,m]])*(1-equals(t,1)*step(m-9.5))
#       # Compute normalized monthly weights, which will be between 0 and 1,
#       # and will some to one.
#       weight[m,t] <- delta[m,t]/sumD
#       # Reorder the weights to go from most recent month (Dec of current
#       # year) to “oldest” month (Jan at past year = Nlag).
#       weightOrdered[(t-1)*12 + (12-m+1)] <- weight[m,t]
#       # For each time into the past, compute the weighted precipitation
#       # variable.
#       for(i in Nlag:Nyrs){
#         antX1[i,m,t] <- weight[m,t]*ppt[i-t+1,m]
#       }
#     } }
#   
#   
#   # Compute sum of deltas (unnormalized weights), to be used to compute
#   # the normalized antecedent weights:
#   for(t in 1:Nlag){
#     sumD1[t] <- sum(delta[,t])
#   }
#   sumD <- sum(sumD1[])
#   # Compute the cumulative monthly weights:
#   for(t in 1:(12*Nlag)){
#     cum.weight[t] <- sum(weightOrdered[1:t])
#   }
#   # Compute the month within year weights (alpha’s = wP,m in Box 1 in main
#   # text); that is, these weights sum to 1 within each past year
#   for(m in 1:12){
#     for(t in 1:Nlag){
#       alpha[m,t] <- delta[m,t]/sum(delta[,t])
#     }
#   }
#   # Compute antecedent precipitation by summing the weighted precipitation
#   # variable over months and past years:
#   for(i in Nlag:Nyrs){
#     for(t in 1:Nlag){
#       ant.sum1[i,t] <- sum(antX1[i,,t])
#     }
#     antX[i] <- sum(ant.sum1[i,])
#   }
#   # Assign priors to the ANPP regression parameters (covariate effects):
#   for(k in 1:6){
#     a[k] ~ dnorm(0,0.0000001)
#   }
#   # Prior for residual (observation) standard deviation, and compute
#   # associated precision:
#   sig ~ dunif(0,100)
#   tau <- pow(sig,-2)
#   # Priors for parameters in the Event missing data model:
#   for(k in 1:4){
#     mu.ev[k] ~ dunif(0,500)
#     sig.ev[k] ~ dunif(0,500)
#     tau.ev[k] <- pow(sig.ev[k],-2)
#   } 
# }", fill=TRUE)
# sink()
# 
# 
# data = list(
#   'block'=block, 
#   'YearID'=YearID, 
#   'Event'=Event, 
#   'ppt'=ppt,
#   'Nlag'=Nlag, 
#   'N'=N, 
#   'Nyrs'=Nyrs, 
#   'Nblocks'=Nblocks,
#   'NPP'=NPP)
# 
# 
# 
# params <- c("a", "mu", "weight", "yr.w", "antX", "Event", "cum.weight")
# 
# 
# 
# # MCMC settings
# ni <- 10000
# nb <- 2000
# nt <- 10
# nc <- 3
# 
# 
# # Run SAM model
# library(jagsUI)

# bothout <- jags(data, params, inits = NULL, "OgleSAM.txt", n.thin=nt, n.chains=nc, n.burnin=nb, n.iter=ni, parallel=T)

Here is the recreation of Figure 2a from the paper. You can see that it matches up quite nicely. Weights deviated a little most likely due to the prior. The code is able to replicate her results.

Now, what would happen if her response variable was a Bernoulii distribution like ours? To do this, I extracted the beta coefficients of the predictors and calculated what the NPP value would be if it were a binary variable (0 or 1).

See the code used to do that.

model_summary <- read.csv("Ogle_data.csv", row.names = 1)

# Extract the means of the a[1] to a[6] parameters
a_means <- model_summary[grep("^a\\[", rownames(model_summary)), "mean"]

# Pull all antX values
antX_full <- model_summary[grep("^antX\\[", rownames(model_summary)), "mean"]

# Align to 52 observed individuals using YearID (which has values from 5 to 91)
antX_post <- antX_full[YearID]

# Use the estimated coefficients
intercept <- a_means[1]
beta_antX <- a_means[2]
beta_event <- a_means[3:6]

# Calculate linear predictor
lp <- intercept + 
      beta_antX * antX_post + 
      beta_event[1] * Event$Event1 +
      beta_event[2] * Event$Event2 +
      beta_event[3] * Event$Event3 +
      beta_event[4] * Event$Event4

# Convert to probability using inverse logit
p <- plogis(lp)

# Simulate binary response
set.seed(2025)  # for reproducibility
NPP_binary <- rbinom(length(p), 1, prob = p)

dataset2$NPP_binary <- NPP_binary

NPP_binary

##  [1] NA  0 NA NA  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## [26]  0  0  0  1  0 NA  0  0  0  0  0  0  0  0  0  1  0  0  1  0  0  0  0 NA NA
## [51] NA NA

The, I reran the model with a Bernoulli distribution for the response variable and produced the following graph. You can see that this lines up with the prior distribution (lines up with the gray credible intervals in Figure 2a). So, the estimated weights that Ogle got in her model from the paper just return the prior in the model with a Bernoulli response variable.

This was what got me originally thinking that a Bernoulli response variable was our “problem”.

Manual Simulations

Alec Schindler was kind enough to write up this quick simulation code that I used to test different response variables and variations. The code forces weights into the code in certain patterns (constant, ascending, descending, triangle, early, late, and random). If the model works, we would expect to see weights that follow that specific distribution. See the code below for the Bernoulli response variable.

I ran the following variations

Smaller PTF variation, 30 days migration, Bernoulli response variable, alpha = -3, beta = 6.5
Larger PTF variation, 30 days migration, Bernoulli response variable, alpha = -3, beta = 6.5
Smaller PTF variation, 30 days migration, Poisson response variable, alpha = -3, beta = 11

In Ogle’s code, the SD of the response variable was around 25. When I make the beta 11, the SD of the response variable is also around 25.

Multiplied PTF by 10, 30 days migration, Poisson response variable, alpha = -3, beta = 1.1

Since our PTF data has boundaries (0,1) I wondered if that impacted anything. I also lowered beta down to again, get the response variable to have a SD around 25.

Multiplied PTF by 10, 30 days migration, Bernoulli response variable, alpha = -3, beta = 0.6

Since our PTF data has boundaries (0,1) I wondered if that impacted anything. I also lowered beta down to make sure the response variable had 0s and 1s.

Here is the code used to develop the simulations. Feel free to play around with it

# mean PTF (proportion of time feeding) for different groups
# (can manipulate these values to test how large differences need to be between groups
#  or how large differences need to be between days to detect weights and effects)
lower_ptf1 <- 0.5 # more likely to succeed
upper_ptf1 <- 0.75
lower_ptf2 <- .25 # more likely to fail
upper_ptf2 <- 0.4

# # another example (larger variation among days)
# lower_ptf1 <- 0.4  # more likely to succeed
# upper_ptf1 <- 0.85
# lower_ptf2 <- 0.1  # more likely to fail
# upper_ptf2 <- 0.55

# length of migration (can change to match ecology of the species)
max_day <- 30

# function to add different weight structures
calc_weights <- function(max_day){
  # constant weights: all weights have the same value
  constant <- rep(1 / max_day, max_day)
  # ascending weights: weights increase over time
  ascending <- rdirichlet(1, 1:max_day)
  # descending weights: weights decrease over time
  descending <- rdirichlet(1, max_day:1)
  # early weights: descending weights early in the time series & 0 later
  early <- c(rdirichlet(1, floor(max_day * 0.25):1), rep(0, ceiling(max_day * 0.75)))
  # late weights: ascending weights late in the time series & 0 earlier
  late <- c(rep(0, ceiling(max_day * 0.75)), rdirichlet(1, 1:floor(max_day * 0.25)))
  # triangle weights: weights increase and then decrease (highest in the middle of the time series)
  triangle <- rdirichlet(1, c(1:floor(median(1:max_day) - 0.5), floor(median(1:max_day)):1))
  # combine weights in a data frame
  data.frame(rel_day = 1:max_day,
             constant = constant, 
             ascending = t(ascending), 
             descending = t(descending), 
             early = early, 
             late = late,
             triangle = t(triangle))
}

# create weights
weights <- calc_weights(max_day)

# number of individuals in simulated data set
nind <- 50
nind1 <- ceiling(nind * 0.4)
nind2 <- floor(nind * 0.6)

# simulate data (for group of birds with higher PTF - aka better reproductive success)
sim_dat1 <- data.frame(id = rep(1:nind1, each = max_day),
                       rel_day = rep(1:max_day, times = nind1),
                       PTF = runif(nind1 * max_day, lower_ptf1, upper_ptf1))

# simulate data (for group of birds with lower PTF)
sim_dat2 <- data.frame(id = rep((nind1 + 1):(nind1 + nind2), each = max_day),
                       rel_day = rep(1:max_day, times = nind2),
                       PTF = runif(nind2 * max_day, lower_ptf2, upper_ptf2))

# combine simulated data and add weights
sim_dat <- bind_rows(sim_dat1, sim_dat2) %>% 
  left_join(weights, by = "rel_day")

# add weights for random individual heterogeneity
sim_dat <- sim_dat %>% 
  mutate(random = rbeta(nind * max_day, 1, 1)) %>% 
  group_by(id) %>% 
  mutate(random = random / sum(random)) %>% 
  ungroup()

# calculate antecedent PTF
sim_dat_ant <- sim_dat %>% 
  mutate(across(constant:random, .names = "PTF_{.col}") * PTF) %>% 
  group_by(id) %>% 
  summarise(across(starts_with("PTF_"), ~ sum(.x), .names = "ant_{.col}")) %>% 
  ungroup()

# function to simulate reproductive success
sim_rs <- function(antPTF, alpha, beta){
  # logistic regression
  logit_p <- alpha + beta * antPTF
  # convert to probability scale
  p <- 1 / (1 + exp(-logit_p))
  # simulate reproductive success
  rbinom(length(antPTF), 1, p)
}

# function to simulate reproductive success POISSON DISTRIBUTION

# sim_rs_poisson <- function(antPTF, alpha, beta){
#   # linear predictor
#   eta <- alpha + beta * antPTF
#   # convert to Poisson rate (lambda) using log link
#   lambda <- exp(eta)
#   # simulate counts
#   rpois(length(antPTF), lambda)
# }
# 

# simulate reproductive success
sim_repro_success <- sim_dat_ant %>%
  mutate(across(starts_with("ant"), 
                # can manipulate alpha/beta based on expected relationship
                ~ sim_rs(.x, alpha = -3, beta = 6.5),  
                .names = "success_{.col}"))

## the likelihoods

# BS[l] ~ dbern(p[ids2[l]])
# logit(p[ids2[l]]) <- alpha + beta_PTF * ant_PTF[ids2[l]]

# BS[l] ~ dpois(lambda[ids2[l]])
# log(lambda[ids2[l]]) <- alpha + beta_PTF * ant_PTF[ids2[l]]

#  specify model
sink("SAM_bern.txt")
cat("
    model {
    
    ### Priors

    # intercept
    alpha ~ dnorm(0, 0.01)
      
    # betas
    beta_PTF ~ dnorm(0, 0.01)
    
    for(k in 1:nobs1) {
      # Dirichlet prior for daily PTF weights
      delta_PTF[ids1[k], days1[k]] ~ dgamma(1,1)
      
      # normalize PTF weights to 1
      weight_PTF[ids1[k], days1[k]] <- 
        delta_PTF[ids1[k], days1[k]] / sum_D1[ids1[k]]
      
      # For each time period into the past, compute weighted PTF variable
      weighted_PTF[ids1[k], days1[k]] <- 
        weight_PTF[ids1[k], days1[k]] * PTF[ids1[k], max_days1[k] - days1[k] + 1]
      
      # Reorder weights from recent to older
      weight_ordered_PTF[ids1[k], max_days1[k] - days1[k] + 1] <- 
        weight_PTF[ids1[k], days1[k]]
      
      # Compute cumulative daily weights
      cum_weight_PTF[ids1[k], days1[k]] <- 
        sum(weight_ordered_PTF[ids1[k], 1:days1[k]])
    }
    
    for(l in 1:nobs2){
      # Compute sum of deltas
      sum_D1[ids2[l]] <- sum(delta_PTF[ids2[l], 1:max_days2[l]])
      
      # calculate antecedent PTF
      ant_PTF[ids2[l]] <- sum(weighted_PTF[ids2[l], 1:max_days2[l]])
    
    ### Likelihood
      BS[l] ~ dpois(lambda[ids2[l]])
      log(lambda[ids2[l]]) <- alpha + beta_PTF * ant_PTF[ids2[l]]

    }
    
  }", fill=TRUE)
sink()

# bundle data
jags_data <- list(
  # k loop
  nobs1 = nrow(sim_dat), 
  ids1 = sim_dat$id, 
  days1 = sim_dat$rel_day,
  max_days1 = rep(max_day, nrow(sim_dat)), 
  PTF = reshape2::acast(sim_dat, id ~ rel_day, value.var = "PTF", na.rm = FALSE),
  
  # l loop
  nobs2 = nrow(sim_repro_success), 
  ids2 = sim_repro_success$id, 
  max_days2 = rep(max_day, nrow(sim_repro_success)))

# add different simulated reproductive success data
jags_data_constant <- jags_data_ascending <- jags_data_descending <- 
  jags_data_early <- jags_data_late <- jags_data_triangle <- jags_data_random <- jags_data

jags_data_constant$BS <- sim_repro_success$success_ant_PTF_constant
jags_data_ascending$BS <- sim_repro_success$success_ant_PTF_ascending
jags_data_descending$BS <- sim_repro_success$success_ant_PTF_descending
jags_data_early$BS <- sim_repro_success$success_ant_PTF_early
jags_data_late$BS <- sim_repro_success$success_ant_PTF_late
jags_data_triangle$BS <- sim_repro_success$success_ant_PTF_triangle
jags_data_random$BS <- sim_repro_success$success_ant_PTF_random

# initial values
inits <- function() {list(alpha = rnorm(1), 
                          beta_PTF = rnorm(1))}

# parameters monitored
params <- c("alpha", 
            "beta_PTF", 
            "p", 
            "weight_PTF", 
            "cum_weight_PTF")

# MCMC settings
ni <- 9000
nb <- 5000
nt <- 1
nc <- 3

1) Smaller PTF variation, 30 days migration, Bernoulli response variable, alpha = -3, beta = 6.5

This is data that is quite similar to mine. The beta and alpha values produced what I thought were reasonable responses. I ran this with 200 birds before I decided that was taking too long.

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights -

late Weights -

Triangle Weights -

Random Weights -

As you can see from all of these, the weights don’t change even when hard coded. I’m not sure why that is.

2) Larger PTF variation, 30 days migration, Bernoulli response variable, alpha = -3, beta = 6.5

This time, I made the variation larger in PTF to see if that would help.

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights -

late Weights -

Triangle Weights -

Random Weights -

As you can see, making the PTF have more variation did not change anything (although variation in ascending is a little greater).

3) Smaller PTF variation, 30 days migration, Poisson response variable, alpha = -3, beta = 11

Next, I simulated a poisson response variable. Hopefully I modifed the function correctly!

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights -

late Weights -

Triangle Weights -

Random Weights -

For this one, we see a bit more variation in the weights, but still not as drastic as the Ogle paper. These are now in a Poisson and have more variation in the response than just a classic 0 vs 1 binary variable. It’s still odd that the weights are doing what we coded them to do.

4) Multiplied PTF by 10, 30 days migration, Poisson response variable, alpha = -3, beta = 1.1

These were the bounds I used for this model. In this case, I made the PTF values greater to see if that would help.

lower_ptf1: 5 upper_ptf1: 7.5 lower_ptf2: 2.5 upper_ptf2: 4

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights -

late Weights -

Triangle Weights -

Random Weights -

We still see some variability that we saw in the other poisson model, just not a lot like in the Ogle paper. Seems like making PTF bigger helped?

5) Multiplied PTF by 10, 30 days migration, Bernoulli response variable, alpha = -3, beta = 0.6

Lastly, I turned back to the Bernoulli response variable once last time with the greater PTF values.

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights - model didn’t work for some reason

late Weights -

Triangle Weights -

Random Weights -

Even though we increased the magnitude of PTF, bernouli still wasn’t able to work.

Conclusions

The main issue is that our data has a Bernoulli response variable. When we simulated a 0 / 1 response variable from Ogle’s own data, the weights just returned the priors.

In terms of why else the model doesn’t work, I’m not sure. The biggest confusion that I (nor Alec) looked into was why the simulation code for the weight wouldn’t return what we hard coded. This would be really interesting to look into.

SAM simulations

Cassidy Waldrep

2025-07-15

Introduction

Our data

METADATA

Looking at Ogle’s data

Manual Simulations

1) Smaller PTF variation, 30 days migration, Bernoulli response variable, alpha = -3, beta = 6.5

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights -

late Weights -

Triangle Weights -

Random Weights -

2) Larger PTF variation, 30 days migration, Bernoulli response variable, alpha = -3, beta = 6.5

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights -

late Weights -

Triangle Weights -

Random Weights -

3) Smaller PTF variation, 30 days migration, Poisson response variable, alpha = -3, beta = 11

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights -

late Weights -

Triangle Weights -

Random Weights -

4) Multiplied PTF by 10, 30 days migration, Poisson response variable, alpha = -3, beta = 1.1

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights -

late Weights -

Triangle Weights -

Random Weights -

5) Multiplied PTF by 10, 30 days migration, Bernoulli response variable, alpha = -3, beta = 0.6

Constant Weights -

Ascending Weights -

Descending Weights -

Early Weights - model didn’t work for some reason

late Weights -

Triangle Weights -

Random Weights -

Conclusions