www.phidot.org

For the abundance parameters in POPAN and Jolly-Seber Lambda -- Burnham models, the relationship between the beta parameter estimates and the real parameter estimates is not given by the inverse link function. It is not clear to me what the functional relationship is between the model coefficients and the real estimates.

Is there some offset term in the formulation of the linear model in MARK?
How can I compute the real estimates of abundance from the beta parameter estimates using RMark?

EDIT 1: Added reproducible example
EDIT 2: Added Jolly model example. Rephrased post title and question.

Code: Select all


# Load data ---------------------------------------------------------------
data(dipper)

# POPAN model -------------------------------------------------------------
# Define formulas
Phi.dot = list(formula =  ~ 1)
p.dot = list(formula =  ~ 1)
pent.dot = list(formula =  ~ 1)
N.dot = list(formula =  ~ 1)

popan.processed <- process.data(dipper, model="POPAN")
popan.ddl <- make.design.data(popan.processed)

# Run model
popanMod <-
  mark(
    popan.processed,
    popan.ddl,
    model.parameters = list(
      Phi = Phi.dot,
      p = p.dot,
      pent = pent.dot,
      N = N.dot
    )
  )

# MARK results
popanReal <- popanMod$results$real
N1.real.mark <- popanReal[grepl("N ", row.names(popanReal)), 1]

# Compute results
popanCompute <- compute.real(model = popanMod, design = popanMod$design.matrix)
N1.real.compute <- popanCompute[grepl("N ", row.names(popanReal)), 1]

# Manual results
beta_popan <- popanMod$results$beta
Beta_N1 <- beta_popan[grepl("N:", row.names(beta_popan)), 1]
N1.real.manual <- exp(Beta_N1)

# RMark / MARK results do not match compute results
all.equal(N1.real.mark, N1.real.compute, scale = 1)

# Compute results match my calculation
all.equal(N1.real.compute, N1.real.manual, scale = 1)

# Jolly model -------------------------------------------------------------
# Define formulas
Phi.dot = list(formula =  ~ 1)
p.dot = list(formula =  ~ 1)
Lambda.dot = list(formula =  ~ 1)

jolly.processed <- process.data(dipper, model="Jolly")
jolly.ddl <- make.design.data(jolly.processed)

# Run model
jollyMod <-
  mark(jolly.processed,
       jolly.ddl,
       model.parameters = list(Phi = Phi.dot, p = p.dot, Lambda = Lambda.dot),
       options= 'SIMANNEAL')

# MARK results
jollyReal <- jollyMod$results$real
N2.real.mark <- jollyReal[grepl("N ", row.names(jollyReal)), 1]

# Compute results
jollyCompute <- compute.real(model = jollyMod, design = jollyMod$design.matrix)
N2.real.compute <- jollyCompute[grepl("N ", row.names(jollyReal)), 1]

# Manual results
beta_jolly <- jollyMod$results$beta
Beta_N2 <- beta_jolly[grepl("N:", row.names(beta_jolly)), 1]
N2.real.manual <- exp(Beta_N2)

# RMark / MARK results do not match compute results
all.equal(N2.real.mark, N2.real.compute, scale = 1)

# Compute results match my calculation
all.equal(N2.real.compute, N2.real.manual, scale = 1)

# Rerun POPAN model with identity link function ---------------------------

# Get default link functions
dummy <-
  make.mark.model(popan.processed,
                  popan.ddl,
                  simplify = FALSE,
                  parm.specific = TRUE)
input.links <- dummy$links

# Get model index for N
log.indices <- popan.ddl$N$model.index

# Assign identity link function
input.links[log.indices] <- "Identity"

# Run model
popanMod <-
  mark(
    popan.processed,
    popan.ddl,
    input.links = input.links,
    model.parameters = list(
      Phi = Phi.dot,
      p = p.dot,
      pent = pent.dot,
      N = N.dot
    )
  )

# MARK results
popanReal <- popanMod$results$real
N1.real.mark <- popanReal[grepl("N ", row.names(popanReal)), 1]

# Compute results
popanCompute <- compute.real(model = popanMod, design = popanMod$design.matrix)
N1.real.compute <- popanCompute[grepl("N ", row.names(popanReal)), 1]

# Manual results
beta_popan <- popanMod$results$beta
Beta_N1 <- beta_popan[grepl("N:", row.names(beta_popan)), 1]
N1.real.manual <- Beta_N1

# RMark / MARK results do not match compute results
all.equal(N1.real.mark, N1.real.compute, scale = 1)

# Compute results match my calculation
all.equal(N1.real.compute, N1.real.manual, scale = 1)



I cross-posted this: viewtopic.php?f=1&t=3716

Here is the solution:

Code: Select all

library(RMark)
library(dplyr)
# Load data ---------------------------------------------------------------
data(dipper)
dipper <- dipper[dipper$sex == "Male",]

# POPAN model -------------------------------------------------------------
# Define formulas
Phi.dot = list(formula =  ~ time)
p.dot = list(formula =  ~ 1)
pent.dot = list(formula =  ~ time)
N.dot = list(formula =  ~ 1)

popan.processed <- process.data(dipper, model="POPAN")
popan.ddl <- make.design.data(popan.processed)

# Run model
popanMod <-
  mark(
    popan.processed,
    popan.ddl,
    model.parameters = list(
      Phi = Phi.dot,
      p = p.dot,
      pent = pent.dot,
      N = N.dot
    )
  )

# MARK results
popanReal <- popanMod$results$real
N1.real.mark <- popanReal[grepl("N ", row.names(popanReal)), 1]

# Compute results
popanCompute <- compute.real(model = popanMod, design = popanMod$design.matrix)
N1.real.compute <- nrow(dipper) + popanCompute[grepl("N ", row.names(popanReal)), 1]

# Manual results
beta_popan <- popanMod$results$beta
Beta_N1 <- beta_popan[grepl("N:", row.names(beta_popan)), 1]
N1.real.manual <- nrow(dipper) +exp(Beta_N1)

# RMark / MARK results do not match compute results
all.equal(N1.real.mark, N1.real.compute, scale = 1, tolerance = 1e-5)

# Compute results match my calculation
all.equal(N1.real.compute, N1.real.manual, scale = 1, tolerance = 1e-5)

# Jolly model -------------------------------------------------------------
# Define formulas
Phi.dot = list(formula =  ~ time)
p.dot = list(formula =  ~ 1)
Lambda.dot = list(formula =  ~ time)

jolly.processed <- process.data(dipper, model="Jolly")
jolly.ddl <- make.design.data(jolly.processed)

# Run model
jollyMod <-
  mark(jolly.processed,
       jolly.ddl,
       model.parameters = list(Phi = Phi.dot, p = p.dot, Lambda = Lambda.dot),
       options= 'SIMANNEAL')

# MARK results
jollyReal <- jollyMod$results$real
N2.real.mark <- jollyReal[grepl("N ", row.names(jollyReal)), 1]

# Compute results
jollyCompute <- compute.real(model = jollyMod, design = jollyMod$design.matrix)
N2.real.compute <- nrow(dipper) + jollyCompute[grepl("N ", row.names(jollyReal)), 1]

# Manual results
beta_jolly <- jollyMod$results$beta
Beta_N2 <- beta_jolly[grepl("N:", row.names(beta_jolly)), 1]
N2.real.manual <- nrow(dipper) + exp(Beta_N2)

# RMark / MARK results do not match compute results
all.equal(N2.real.mark, N2.real.compute, scale = 1, tolerance = 1e-5)

# Compute results match my calculation
all.equal(N2.real.compute, N2.real.manual, scale = 1, tolerance = 1e-5)


Also for the POPAN model there is a function in RMark called popan.derived that derives a set of different abundance estimates. Also, you can get some of those from model$results$derived which are the derived parameters from MARK.

natprice wrote:I cross-posted this: viewtopic.php?f=1&t=3716

Indeed - please don't do that going forward. One thread, one forum.

Management

Thanks for your help!