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)