Compute real estimates from beta parameters

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

Compute real estimates from beta parameters

Postby natprice » Wed Oct 17, 2018 3:22 pm

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)

natprice
 
Posts: 5
Joined: Wed Nov 01, 2017 10:50 am

Re: Compute real estimates from beta parameters

Postby natprice » Mon Oct 22, 2018 4:49 pm

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)
natprice
 
Posts: 5
Joined: Wed Nov 01, 2017 10:50 am

Re: Compute real estimates from beta parameters

Postby jlaake » Mon Oct 22, 2018 4:57 pm

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.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: Compute real estimates from beta parameters

Postby egc » Mon Oct 22, 2018 6:32 pm

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
egc
Site Admin
 
Posts: 201
Joined: Thu May 15, 2003 3:25 pm

Re: Compute real estimates from beta parameters

Postby natprice » Mon Oct 22, 2018 9:09 pm

Thanks for your help!
natprice
 
Posts: 5
Joined: Wed Nov 01, 2017 10:50 am


Return to RMark

Who is online

Users browsing this forum: No registered users and 17 guests