rsurvival-analysisrecoverystate-space

State-SpaceMark-Recovery Model for Survival and Recovery by Age


Can anyone help me with my issue? I simulated using a state-space model to estiamte recovery(p) and survival (phi) by age class.

My estimates are inaccurate.

I get the following:


    beta [1] = 0.331
    beta [2] = 0.634 
    delta [1] = 0.469
    delta [2] = 0.529

According to my simulated data (or what I think should happen)


    beta [1] = 0.8
    beta [2] = 0.8
    delta [1] = 0.2
    delta [2] = 0.1

Any help would be extremely appreciated! Also let me know if you can't read my code.


    ### Finding first capture year
    get.first <- function(x)min(which(x!=0))
    
    known.state.mr <- function (mr) {
      state <- matrix (NA, nrow=dim(mr)[1], ncol= dim(mr)[2])
      rec <- which(rowSums(mr)==2)
      for (i in 1: length(rec)){
        n1<- min(which(mr[rec[i],]==1)) #n1 time of capture
        n2<- max(which(mr[rec[i],]==1)) # Recovery time
        state[rec[i],n1:n2]<-1 
        state[rec[i],n1]<-NA 
        state[rec[i],n2:dim(mr)[2]]<- 0 
      }
      return(state)
    }
    
    
    mr.init.z<-function(mr){
      ch<-matrix(NA, nrow=dim(mr)[1],ncol=dim(mr)[2])
      rec<- which(rowSums(mr)==1) 
      for(i in 1:length(rec)){ 
        n1<-which(mr[rec[i],]==1)
        ch[rec[i], n1:dim(mr)[2]]<-0
        ch[rec[i],n1] <- NA 
      }
      return (ch)
    }
    
    # Define function to simulate mark-recovery data
    simul.mr <- function(S, R, marked){
      n.occasions <- dim(S)[2]
      MR <- matrix(NA, ncol = n.occasions+1, nrow = sum(marked))
      # Define a vector with the occasion of marking
      mark.occ <- rep(1:n.occasions, marked)
      # Fill the CH matrix
      for (i in 1:sum(marked)){
        MR[i, mark.occ[i]] <- 1    
        for (t in mark.occ[i]:n.occasions){
          # Bernoulli trial: has individual survived occasion? 
          sur <- rbinom(1, 1, S[i,t])
          if (sur==1) next     
           
          rp <- rbinom(1, 1, R[i,t])
          if (rp==0){
            MR[i,t+1] <- 0
            break
          }
          if (rp==1){
            MR[i,t+1] <- 1
            break
          }
        } #t
      } #i
      # Replace the NA in the file by 0
      MR[which(is.na(MR))] <- 0
      return(MR)
    }
    
    
    ########### Simulate Data #########################
    
    # Juvenile Matrix
    n.occasions <- 14                 # Number of release occasions
    marked <- rep(50, n.occasions)    # Annual number of newly 
    marked individuals
    as <- rep(0.8, n.occasions)
    ar <- rep(0.2, n.occasions)
    
    # Define matrices with survival and recovery probabilities
    aS <- matrix(as, ncol = n.occasions, nrow = sum(marked))
    aR <- matrix(ar, ncol = n.occasions, nrow = sum(marked))
    
    
    # Execute function
    MR.A <- simul.mr(aS, aR, marked)
    
    
    ### Juvenile Matrix #####
    js <- rep(0.8, n.occasions)
    jr <- rep(0.1, n.occasions)
    
    # Define matrices with survival and recovery probabilities
    jS <- matrix(js, ncol = n.occasions, nrow = sum(marked))
    jR <- matrix(jr, ncol = n.occasions, nrow = sum(marked))
    
    # Execute function
    MR.J <- simul.mr(jS, jR, marked)
    
    get.first <- function(x) min(which(x!=0))
    f.j <- apply(MR.J, 1, get.first)
    f.a <- apply(MR.A, 1, get.first)
    
    # Create matrices X indicating age classes
    x.j <- matrix(NA, ncol = dim(MR.J)[2]-1, nrow = dim(MR.J)[1])
    x.a <- matrix(NA, ncol = dim(MR.A)[2]-1, nrow = dim(MR.A)[1])
    
    for (i in 1:dim(MR.J)[1]){
      for (t in f.j[i]:(dim(MR.J)[2]-1)){
        x.j[i,t] <- 2
        x.j[i,f.j[i]] <- 1   
      } #t
    } #i
    
    for (i in 1:dim(MR.A)[1]){
      for (t in f.a[i]:(dim(MR.A)[2]-1)){
        x.a[i,t] <- 2
      } #t
    } #i
    
    MR <- rbind(MR.J, MR.A)
    f <- c(f.j, f.a)
    x <- rbind(x.j, x.a)
    
    
    # Specify model in BUGS language
    sink("Model_MR")
    cat("
    model {
    # Priors and constraints
    for (i in 1:nind){
       for (t in f[i]:(n.occasions-1)){
          
          phi[i,t] <- beta[x[i,t]]
          p[i,t] <- delta[x[i,t]]
          
          } #t
       } #i
       
    for (u in 1:2){
       
       beta[u] ~ dunif(0,1) # Priors for age-specific survival
       delta[u] ~ dunif(0,1)
       
    }
       
    mean.p ~ dunif(0, 1)                  # Prior for mean recapture
    # Likelihood 
    for (i in 1:nind){
       
       # Define latent state at first capture
       
       z[i,f[i]] <- 1
       for (t in (f[i]+1):n.occasions){
          
          # State process
          
          z[i,t] ~ dbern(mu1[i,t])
          mu1[i,t] <- phi[i,t-1] * z[i,t-1]
          
          # Observation process
          
          y[i,t] ~ dbern(mu2[i,t])
          mu2[i,t] <- p[i,t-1] * (z[i,t-1]-z[i,t])
          
          } #t
       } #i
    }
    ",fill = TRUE)
    sink()
    
    # Bundle data
    bugs.data <- list(y = MR, 
                      f = f, 
                      nind = dim(MR)[1],
                      n.occasions = dim(MR)[2],
                      z = known.state.mr(MR), 
                      x = x)
    
    # Initial values
    inits <- function(){list(z = mr.init.z(MR), 
                             beta = runif(2, 0, 1), 
                             mean.p = runif(1, 0, 1))}  
    
    # Parameters monitored
    parameters <- c("beta", "mean.p","delta")
    
    # MCMC settings
    ni <- 2000
    nt <- 3
    nb <- 1000
    nc <- 3
    
    # Call WinBUGS from R (BRT 3 min)
    cjs.age <- jags(bugs.data, inits, parameters, "Model_MR", 
    n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
    
    print(cjs.age, digits = 3


Solution

  • I figured out my issue... I was simulating the data wrong. Please don't come after me! Here's the corrected code.

    The issue was I didn't simulate transition probability for when juveniles become adults.

    Now you know!

    
        # Juvenile Matrix
        n.occasions <- 10                 # Number of release occasions
        marked.j <- rep(50, n.occasions)    # Annual number of newly marked individuals
        marked.a <- rep(50, n.occasions)
        marked.j <- rep(200, n.occasions-1) # Annual number of newly marked juveniles
        marked.a <- rep(30, n.occasions-1)  # Annual number of newly marked adults
        phi <- rep(0.8, n.occasions-1)        # Recapture
        p.juv <- 0.3                      # Juvenile annual survival
        p.ad <- 0.65 
        p.j <- c(p.juv, rep(p.ad, n.occasions-2))
        p.a <- rep(p.ad, n.occasions-1)
        
        P.J <- matrix(0, ncol = n.occasions-1, nrow = sum(marked.j))
        
        for (i in 1:length(marked.j)){
          P.J[(sum(marked.j[1:i])-marked.j[i]+1):sum(marked.j[1:i]),i:(n.occasions-1)] <- matrix(rep(p.j[1:(n.occasions-i)],marked.j[i]), ncol = n.occasions-i, byrow = TRUE)
        }
        
        PHI.J <- matrix(rep(phi, sum(marked.j)), ncol = n.occasions-1, nrow = sum(marked.j), byrow = TRUE)
        P.A <- matrix(rep(p.a, sum(marked.a)), ncol = n.occasions-1, nrow = sum(marked.a), byrow = TRUE)
        PHI.A <- matrix(rep(phi, sum(marked.a)), ncol = n.occasions-1, nrow = sum(marked.a), byrow = TRUE)
        
        # Execute function
        MR.A <- simul.mr(PHI.A, P.A, marked.a)
        
        # Execute function
        MR.J <- simul.mr(PHI.J, P.J, marked.j)
        
        get.first <- function(x) min(which(x!=0))
        f.j <- apply(MR.J, 1, get.first)
        f.a <- apply(MR.A, 1, get.first)
        
        # Create matrices X indicating age classes
        x.j <- matrix(NA, ncol = dim(MR.J)[2]-1, nrow = dim(MR.J)[1])
        x.a <- matrix(NA, ncol = dim(MR.A)[2]-1, nrow = dim(MR.A)[1])
        
        for (i in 1:dim(MR.J)[1]){
          for (t in f.j[i]:(dim(MR.J)[2]-1)){
            x.j[i,t] <- 2
            x.j[i,f.j[i]] <- 1   
          } #t
        } #i
        
        for (i in 1:dim(MR.A)[1]){
          for (t in f.a[i]:(dim(MR.A)[2]-1)){
            x.a[i,t] <- 2
          } #t
        } #i
        
        MR <- rbind(MR.J, MR.A)
        f <- c(f.j, f.a)
        x <- rbind(x.j, x.a)