SCR - NA or NAN values for z[]
21 views
Skip to first unread message
Erika Anderson
Jun 14, 2024, 2:09:23 PMJun 14
to nimble-users
Hello Nimble users,
I have successfully ran an SCR model with a continuous state space and am now working on a discrete space. I am running into an issue where the model is returning NA values for z at the nimbleModel() phase. I've tried a number of different iterations (e.g. different initial values, simplifying the model by removing alpha and beta coefficients) but in the final output each model results in either NA for all results or the same value for each parameter for all posterior samples (i.e. N=147 for all 6000 samples, b0=4.033 for all 600 samples, etc). I have included my code below, and some screen shots of what my data looks like in R (note there are couple of extra elements that I have removed from the model code - nimble tells me they are ignored when running nimbleModel(). I'm not sure if this is an issue with my data format, initial values, or something else. Any thoughts are appreciated!
Erika
#################################################
marten <- nimbleCode ({
# Priors and constraints
#psi ~ dunif(0,1) # Probability of inclusion (aka augmented individuals are real)
#alpha0 ~ dunif(-20,20) # intercept for detection probability
#alpha1 ~ dnorm(0,0.01) # Effect of previous detection on capture probability
beta0 ~ dunif(-20,20) # intercept for density
beta1 ~ dnorm(0,0.01) # Effect of ddi on density
#beta2 ~ dnorm(0,0.01) # Effect of snow
#beta3 ~ dnorm(0,0.01) # Effect of ladder fuels
#beta4 ~ dnorm(0,0.01) # Effect of fisher
#beta5 ~ dnorm(0,0.01) # Effect of bobcat
sigma ~ dunif(0,2400)
p0 ~ dunif(0,5) # lambda = baseline trap encounter rate
for(g in 1:ngrid) { # For each grid cell
mu[g] <- exp(beta0 + beta1 * ddi[g] + beta2 * snow[g] + beta3 * ladder[g] + beta4 * fisher[g] + beta5 * bobcat[g]) * pixArea
# Priors and constraints
#psi ~ dunif(0,1) # Probability of inclusion (aka augmented individuals are real)
#alpha0 ~ dunif(-20,20) # intercept for detection probability
#alpha1 ~ dnorm(0,0.01) # Effect of previous detection on capture probability
beta0 ~ dunif(-20,20) # intercept for density
beta1 ~ dnorm(0,0.01) # Effect of ddi on density
#beta2 ~ dnorm(0,0.01) # Effect of snow
#beta3 ~ dnorm(0,0.01) # Effect of ladder fuels
#beta4 ~ dnorm(0,0.01) # Effect of fisher
#beta5 ~ dnorm(0,0.01) # Effect of bobcat
sigma ~ dunif(0,2400)
p0 ~ dunif(0,5) # lambda = baseline trap encounter rate
for(g in 1:ngrid) { # For each grid cell
mu[g] <- exp(beta0 + beta1 * ddi[g] + beta2 * snow[g] + beta3 * ladder[g] + beta4 * fisher[g] + beta5 * bobcat[g]) * pixArea
probs[g] <- mu[g] / EN # probability of an activity center being located in grid cell g
} #g, grid cells
EN <- sum(mu[1:ngrid])
psi <- EN/M
for (i in 1:M){ # For each individual (n+nz=M)
# Make activity centers
z[i] ~ dbern(psi)
s[i] ~ dcat(probs[1:ngrid])
x0g[i] <- grid[s[i],1]
y0g[i] <- grid[s[i],2]
#logit(p0[i]) <- alpha0 + alpha1 * sex[i]
for (j in 1:J){ # For each hair snare
Dist2[i,j] <- sqrt((x0g[i] - trapmatrix[j,1])^2 + (y0g[i] - trapmatrix[j,2])^2)
p[i,j] <- z[i] * p0 * exp(-(1 / (2 * sigma * sigma)) * Dist2[i,j] * Dist2[i,j])
y[i,j] ~ dbern(p[i,j]) # p = probability of detection/encounter
} #j, hair snares
} #i, individuals
N <- sum(z[]) # population size (N)
})
#################################################
load("Data/NimbleData_inits.RData") #V1 inits
load("Data/NimbleData_constants.RData")
load("Data/NimbleData_data_V2.RData")
inits2 <- list (z=inits$z,
psi=0.5,
beta0 =runif(1, -5, 5),
s= inits$s,
sigma= 870,
p0=0.1) #V2 inits
inits2$z[is.na(inits2$z)] <- 1 #V3 inits
# parameters
parameters_fish <- c("N", "EN", "psi", "mu", "sigma", "z", "s", "beta0", "beta1", "beta2", "beta3", "beta4", "beta5")
} #g, grid cells
EN <- sum(mu[1:ngrid])
psi <- EN/M
for (i in 1:M){ # For each individual (n+nz=M)
# Make activity centers
z[i] ~ dbern(psi)
s[i] ~ dcat(probs[1:ngrid])
x0g[i] <- grid[s[i],1]
y0g[i] <- grid[s[i],2]
#logit(p0[i]) <- alpha0 + alpha1 * sex[i]
for (j in 1:J){ # For each hair snare
Dist2[i,j] <- sqrt((x0g[i] - trapmatrix[j,1])^2 + (y0g[i] - trapmatrix[j,2])^2)
p[i,j] <- z[i] * p0 * exp(-(1 / (2 * sigma * sigma)) * Dist2[i,j] * Dist2[i,j])
y[i,j] ~ dbern(p[i,j]) # p = probability of detection/encounter
} #j, hair snares
} #i, individuals
N <- sum(z[]) # population size (N)
})
#################################################
load("Data/NimbleData_inits.RData") #V1 inits
load("Data/NimbleData_constants.RData")
load("Data/NimbleData_data_V2.RData")
inits2 <- list (z=inits$z,
psi=0.5,
beta0 =runif(1, -5, 5),
s= inits$s,
sigma= 870,
p0=0.1) #V2 inits
inits2$z[is.na(inits2$z)] <- 1 #V3 inits
# parameters
parameters_fish <- c("N", "EN", "psi", "mu", "sigma", "z", "s", "beta0", "beta1", "beta2", "beta3", "beta4", "beta5")
Erika Anderson
Jun 14, 2024, 2:58:59 PMJun 14
to nimble-users
Follow up: I realize my code above wasn't the correct version. Here it is fixed.
#################################################
marten <- nimbleCode ({
# Priors and constraints
# Priors and constraints
beta0 ~ dunif(-20,20) # intercept for density
beta1 ~ dnorm(0,0.01) # Effect of ddi on density
beta2 ~ dnorm(0,0.01) # Effect of snow
beta3 ~ dnorm(0,0.01) # Effect of ladder fuels
beta4 ~ dnorm(0,0.01) # Effect of fisher
beta5 ~ dnorm(0,0.01) # Effect of bobcat
sigma ~ dunif(0,2400)
p0 ~ dunif(0,5) # lambda = baseline trap encounter rate
for(g in 1:ngrid) { # For each grid cell
mu[g] <- exp(beta0 + beta1 * ddi[g]
sigma ~ dunif(0,2400)
p0 ~ dunif(0,5) # lambda = baseline trap encounter rate
for(g in 1:ngrid) { # For each grid cell
mu[g] <- exp(beta0 + beta1 * ddi[g]
+ beta2 * snow[g] + beta3 * ladder[g] + beta4 * fisher[g] + beta5 * bobcat[g]) * pixArea
probs[g] <- mu[g] / EN # probability of an activity center being located in grid cell g
} #g, grid cells
EN <- sum(mu[1:ngrid])
psi <- EN/M
for (i in 1:M){ # For each individual (n+nz=M)
# Make activity centers
z[i] ~ dbern(psi)
s[i] ~ dcat(probs[1:ngrid])
x0g[i] <- grid[s[i],1]
y0g[i] <- grid[s[i],2]
#logit(p0[i]) <- alpha0 + alpha1 * sex[i]
for (j in 1:J){ # For each hair snare
Dist2[i,j] <- sqrt((x0g[i] - trapmatrix[j,1])^2 + (y0g[i] - trapmatrix[j,2])^2)
p[i,j] <- z[i] * p0 * exp(-(1 / (2 * sigma * sigma)) * Dist2[i,j] * Dist2[i,j])
y[i,j] ~ dbern(p[i,j]) # p = probability of detection/encounter
} #j, hair snares
} #i, individuals
N <- sum(z[]) # population size (N)
})
} #g, grid cells
EN <- sum(mu[1:ngrid])
psi <- EN/M
for (i in 1:M){ # For each individual (n+nz=M)
# Make activity centers
z[i] ~ dbern(psi)
s[i] ~ dcat(probs[1:ngrid])
x0g[i] <- grid[s[i],1]
y0g[i] <- grid[s[i],2]
#logit(p0[i]) <- alpha0 + alpha1 * sex[i]
for (j in 1:J){ # For each hair snare
Dist2[i,j] <- sqrt((x0g[i] - trapmatrix[j,1])^2 + (y0g[i] - trapmatrix[j,2])^2)
p[i,j] <- z[i] * p0 * exp(-(1 / (2 * sigma * sigma)) * Dist2[i,j] * Dist2[i,j])
y[i,j] ~ dbern(p[i,j]) # p = probability of detection/encounter
} #j, hair snares
} #i, individuals
N <- sum(z[]) # population size (N)
})
Heather Gaya
Jun 14, 2024, 3:06:35 PMJun 14
to Erika Anderson, nimble-users
Hi Erika,
It's hard to tell without having your data in front of me, but I suspect the problem is how psi is initialized. Since psi is just sum(mu)/M, it's possible that sum(mu) is > M. Try increasing M or initializing all your betas as much lower. An easy to do this is often to initialize all but one of them 0 and then randomize the intercept and the remaining beta.
I'm not sure how you're running your model, but if you do:
prepnim <- nimbleModel(code = yourmodel, constants = nc.Foo,
data = nd.Foo, inits = ni.Foo)
prepnim$calculate()
and don't see an NA or a -INF or something, then you know you've fixed the problem. You can also check by running:
data = nd.Foo, inits = ni.Foo)
prepnim$calculate()
and don't see an NA or a -INF or something, then you know you've fixed the problem. You can also check by running:
prepnim$initializeInfo()
prepnim$psi
prepnim$z
to get a better idea of what the model is initializing these values at.
Hope that helps!
Best,
Heather Gaya
Postdoctoral Researcher
University of Georgia
--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/34146171-5eb0-4809-b0f3-1e6a80ab2a99n%40googlegroups.com.
Erika Anderson
Jun 14, 2024, 7:59:07 PMJun 14
to nimble-users
Hi Heather,
Thanks so much for the quick reply. I took another look at the coefficients using techniques you suggested (e.g. prepnim$psi) and saw that not only was psi returning NA, but mu and all beta (1-5) coefficients were also coming back as NA. I tried a simple solution of setting initial values for the beta coefficients and the model successfully ran with results that make sense!
Reply all
Reply to author
Forward
0 new messages