Spatial-Mark-Resight-Marginal/NimbleModel SMR Multisession Poisson Dcov Marginal.R at main · benaug/Spatial-Mark-Resight-Marginal · GitHub

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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
NimModel <- nimbleCode({
  #detection function priors - shared across sessions
  lam0.fixed ~ dunif(0,15)
  sigma.fixed ~ dunif(0,10)

  #sample type observation probabilities, shared over sessions.
  alpha.marked[1] <- 1
  alpha.marked[2] <- 1
  alpha.marked[3] <- 1
  alpha.unmarked[1] <- 1
  alpha.unmarked[2] <- 1
  theta.marked.fixed[1:3] ~ ddirch(alpha.marked[1:3])
  theta.unmarked.fixed[1] <- 0
  theta.unmarked.fixed[2:3] ~ ddirch(alpha.unmarked[1:2])

  for(g in 1:N.session){
    #not sharing Dcov parameters over sessions, but can do that
    D0[g] ~ dunif(0,100) #uninformative, diffuse dnorm on log scale can cause neg bias
    D.beta1[g] ~ dnorm(0,sd=10)
    #Density model
    D.intercept[g] <- D0[g]*cellArea[g]
    lambda.cell[g,1:n.cells[g]] <- InSS[g,1:n.cells[g]]*exp(D.beta1[g]*D.cov[g,1:n.cells[g]])
    pi.cell[g,1:n.cells[g]] <- lambda.cell[g,1:n.cells[g]]/pi.denom[g] #expected proportion of total N in cell c
    pi.denom[g] <- sum(lambda.cell[g,1:n.cells[g]])
    lambda.N[g] <- D.intercept[g]*pi.denom[g] #Expected N
    N[g] ~ dpois(lambda.N[g]) #realized N in state space.

    #plug in shared df parameter for each session. Must use lam0[g] and sigma[g] here for custom update.
    #alternatively, can be estimated separately or with random effects.
    lam0[g] <- lam0.fixed
    sigma[g] <- sigma.fixed

    #sample type observation model priors (Dirichlet)
    #If not shared across sessions, use this
    #into g indices here
    # alpha.marked[g,1] <- 1
    # alpha.marked[g,2] <- 1
    # alpha.marked[g,3] <- 1
    # alpha.unmarked[g,1] <- 1
    # alpha.unmarked[g,2] <- 1
    # theta.marked[g,1:3] ~ ddirch(alpha.marked[g,1:3])
    # theta.unmarked[g,1] <- 0
    # theta.unmarked[g,2:3] ~ ddirch(alpha.unmarked[g,1:2])
    #if shared over sessions use this
    theta.marked[g,1:3] <- theta.marked.fixed[1:3]
    theta.unmarked[g,1:3] <- theta.unmarked.fixed[1:3]

    #Marked individuals first
    for(i in 1:n.marked[g]){
      #dunif() here implies uniform distribution within a grid cell
      #also tells nimble s's are in continuous space, not discrete
      s[g,i,1] ~ dunif(xlim[g,1],xlim[g,2])
      s[g,i,2] ~ dunif(ylim[g,1],ylim[g,2])
      #get cell s_i lives in using look-up table
      s.cell[g,i] <- cells[g,trunc(s[g,i,1]/res[g])+1,trunc(s[g,i,2]/res[g])+1]
      #categorical likelihood for this cell, equivalent to zero's trick
      #also disallowing s's in non-habitat
      dummy.data[g,i] ~ dCell(pi.cell[g,s.cell[g,i]])
      lam[g,i,1:J[g]] <- GetDetectionRate(s = s[g,i,1:2], X = X[g,1:J[g],1:2], J=J[g],sigma=sigma[g], lam0=lam0[g], z=z[g,i])
      y.mID[g,i,1:J[g]] ~ dPoissonVector(lam[g,i,1:J[g]]*K1D[g,1:J[g]]*theta.marked[g,1],z=z[g,i]) #marked and identified detections
    }

    #Then unmarked individuals
    for(i in (n.marked[g]+1):M[g]){
      #dunif() here implies uniform distribution within a grid cell
      #also tells nimble s's are in continuous space, not discrete
      s[g,i,1] ~ dunif(xlim[g,1],xlim[g,2])
      s[g,i,2] ~ dunif(ylim[g,1],ylim[g,2])
      #get cell s_i lives in using look-up table
      s.cell[g,i] <- cells[g,trunc(s[g,i,1]/res[g])+1,trunc(s[g,i,2]/res[g])+1]
      #categorical likelihood for this cell, equivalent to zero's trick
      #also disallowing s's in non-habitat
      dummy.data[g,i] ~ dCell(pi.cell[g,s.cell[g,i]])
      lam[g,i,1:J[g]] <- GetDetectionRate(s = s[g,i,1:2], X = X[g,1:J[g],1:2], J=J[g],sigma=sigma[g], lam0=lam0[g], z=z[g,i])
    }

    #Unidentified detections by type
    #1 marked with no ID detections
    bigLam.marked[g,1:J[g]] <- GetbigLam(lam=lam[g,1:n.marked[g],1:J[g]],z=z[g,1:n.marked[g]])
    lam.mnoID[g,1:J[g]] <- bigLam.marked[g,1:J[g]]*K1D[g,1:J[g]]*theta.marked[g,2]
    y.mnoID[g,1:J[g]] ~ dPoissonVector(lam.mnoID[g,1:J[g]],z=1) #plug in z=1 to reuse dPoissonVector

    #2 unmarked detections
    bigLam.unmarked[g,1:J[g]] <- GetbigLam(lam=lam[g,(n.marked[g]+1):M[g],1:J[g]],z=z[g,(n.marked[g]+1):M[g]])
    lam.um[g,1:J[g]] <- bigLam.unmarked[g,1:J[g]]*K1D[g,1:J[g]]*theta.unmarked[g,2]
    y.um[g,1:J[g]] ~ dPoissonVector(lam.um[g,1:J[g]],z=1) #plug in z=1 to reuse dPoissonVector

    #3 unknown marked status
    lam.unk[g,1:J[g]] <- bigLam.marked[g,1:J[g]]*K1D[g,1:J[g]]*theta.marked[g,3] +
      bigLam.unmarked[g,1:J[g]]*K1D[g,1:J[g]]*theta.unmarked[g,3]
    y.unk[g,1:J[g]] ~ dPoissonVector(lam.unk[g,1:J[g]],z=1) #plug in z=1 to reuse dPoissonVector

    #If you have telemetry
    for(i in 1:n.tel.inds[g]){
      for(m in 1:n.locs.ind[g,i]){
        locs[g,i,m,1] ~ dnorm(s[g,tel.inds[g,i],1],sd=sigma[g])
        locs[g,i,m,2] ~ dnorm(s[g,tel.inds[g,i],2],sd=sigma[g])
      }
    }
  }
})# end model