multistate model with unobservable states

questions concerning analysis/theory using the R package 'marked'

multistate model with unobservable states

Postby markmiller » Mon Apr 03, 2023 9:25 am

I created a deterministic data set for a 4-state Cormack-Jolly-Seber model with live recaptures and four capture occasions. I analyzed that data set with the most-general model using MARK, RMark and marked. All three approaches returned essentially the same parameter estimates. Those estimates matched the true values I used to create the data.

Then I imposed a bunch of constraints on the model [creating two unobservable states] and recreated the data set. This time MARK and RMark returned essentially the same parameter estimates. Those estimates matched the true values used to create the data. However, the parameter estimates returned by marked are incorrect.

Could somebody suggest what I am doing incorrectly with 'marked'?

I provide the marked and the RMark code below for the constrained model as well as the parameter estimates returned by MARK, RMark and marked. I also provide the full data set for the constrained model. I posted the original unconstrained data set in the RMark section of this forum a few days ago.

Note that the constrained data set only includes two states when using marked but I added the following two lines when using RMark to create the other two states. I did not know how else to do it.

I would really prefer to run this model with marked rather than with MARK or RMark because marked can do somethings that others cannot.

Code: Select all
0003 1  ;
0004 1  ;


Here is the marked code with automated messages included as comments.

Code: Select all
library(marked)
# Loading required package: lme4
# Loading required package: Matrix
# Loading required package: parallel
# This is marked 1.2.6

library(R2admb)
# Attaching package: ‘R2admb’
#
# The following object is masked from ‘package:marked’:
#
#     setup_admb


library(TMB)
library(ggplot2)


#####@@@@@

# this example requires admb
# mwm: Use ?crm to locate the 'prepare_admb' function.

prepare_admb=function()
{
  Sys.setenv(PATH = paste("C:/ADMB-13.1/bin;C:/ADMB-13.1/utilities;C:/ADMB-13.1/utilities/mingw64/bin",
        Sys.getenv("PATH"), sep = ";"))
    Sys.setenv(ADMB_HOME = "C:/ADMB-13.1")
    invisible()
}
prepare_admb()

#####@@@@@


# Convert a 'MARK' data set into a 'marked' data set:
fourstate.raw.pre <- read.table("fake_2obs_2unobs_state_data_April3_2023.inp",
                header = FALSE, stringsAsFactors = FALSE, na.strings = "NA",
                colClasses = c("character", "numeric", "character"))
head(fourstate.raw.pre)
#     V1       V2 V3
# 1 1100 136.2051  ;
# 2 1101  23.5197  ;
# 3 1102  16.6752  ;
# 4 1110  63.5040  ;
# 5 1111  46.6560  ;
# 6 1112  19.4400  ;


colnames(fourstate.raw.pre) <- c("ch", "freq", "the.end")
fourstate.raw <- fourstate.raw.pre[, !colnames(fourstate.raw.pre) %in% "the.end"]
head(fourstate.raw)
#     ch     freq
# 1 1100 136.2051
# 2 1101  23.5197
# 3 1102  16.6752
# 4 1110  63.5040
# 5 1111  46.6560
# 6 1112  19.4400

tail(fourstate.raw)
#      ch   freq
# 73 0010 490.00
# 74 0011 360.00
# 75 0012 150.00
# 76 0020 595.75
# 77 0021  68.25
# 78 0022 336.00

ms4.processed = process.data(fourstate.raw, model = "Mscjs", strata.labels=c("1","2","3","4"))

ms4.ddl = make.design.data(ms4.processed)


# examine data
head(ms4.ddl$S)
#   id occ stratum time cohort age freq Time Cohort Age order
# 1  1   1       1    1      3   0  490    0      2   0     1
# 2  1   1       2    1      3   0  490    0      2   0     2
# 3  1   1       3    1      3   0  490    0      2   0     3
# 4  1   1       4    1      3   0  490    0      2   0     4
# 5  1   2       1    2      3   0  490    1      2   0     5
# 6  1   2       2    2      3   0  490    1      2   0     6

head(ms4.ddl$p)
#   id occ stratum time cohort age freq Time Cohort Age fix order
# 1  1   2       1    2      3   0  490    0      2   0   0     1
# 2  1   2       2    2      3   0  490    0      2   0   0     2
# 3  1   2       3    2      3   0  490    0      2   0   0     3
# 4  1   2       4    2      3   0  490    0      2   0   0     4
# 5  1   3       1    3      3   0  490    1      2   0   1     5
# 6  1   3       2    3      3   0  490    1      2   0   1     6

head(ms4.ddl$Psi)
#   id occ stratum tostratum time cohort age freq Time Cohort Age fix order
# 1  1   1       1         1    1      3   0  490    0      2   0   1     1
# 2  1   1       1         2    1      3   0  490    0      2   0  NA     2
# 3  1   1       1         3    1      3   0  490    0      2   0  NA     3
# 4  1   1       1         4    1      3   0  490    0      2   0  NA     4
# 5  1   1       2         1    1      3   0  490    0      2   0  NA     5
# 6  1   1       2         2    1      3   0  490    0      2   0   1     6

# constrain survival probabilities:
# 1. constant and equal in States 1 & 3 and in States 2 & 4
# 2. constant over time within a state.
ms4.ddl$S$mystratum[ms4.ddl$S$stratum %in% c(1,3)] <- 1
ms4.ddl$S$mystratum[ms4.ddl$S$stratum %in% c(2,4)] <- 2
ms4.ddl$S$mystratum <- factor(ms4.ddl$S$mystratum)
head(ms4.ddl$S)
#   id occ stratum time cohort age freq Time Cohort Age order mystratum
# 1  1   1       1    1      3   0  490    0      2   0     1         1
# 2  1   1       2    1      3   0  490    0      2   0     2         2
# 3  1   1       3    1      3   0  490    0      2   0     3         1
# 4  1   1       4    1      3   0  490    0      2   0     4         2
# 5  1   2       1    2      3   0  490    1      2   0     5         1
# 6  1   2       2    2      3   0  490    1      2   0     6         2

# Constrain detection probability
# 1. constant and equal in States 1 & 2
# 2. zero in States 3 & 4
ms4.ddl$p$mystratum[ms4.ddl$p$stratum %in% c(1,2)] <- 1
ms4.ddl$p$mystratum[ms4.ddl$p$stratum %in% c(3,4)] <- 2
ms4.ddl$p$mystratum <- factor(ms4.ddl$p$mystratum)
ms4.ddl$p$fix <- NA
ms4.ddl$p$fix[ms4.ddl$p$mystratum == 2] <- 0
head(ms4.ddl$p)
#   id occ stratum time cohort age freq Time Cohort Age fix order mystratum
# 1  1   2       1    2      3   0  490    0      2   0  NA     1         1
# 2  1   2       2    2      3   0  490    0      2   0  NA     2         1
# 3  1   2       3    2      3   0  490    0      2   0   0     3         2
# 4  1   2       4    2      3   0  490    0      2   0   0     4         2
# 5  1   3       1    3      3   0  490    1      2   0  NA     5         1
# 6  1   3       2    3      3   0  490    1      2   0  NA     6         1

# constrain transition probability = 0 for transitions that are not possible
ms4.ddl$Psi$fix[ms4.ddl$Psi$stratum==1 & ms4.ddl$Psi$tostratum == 4] = 0
ms4.ddl$Psi$fix[ms4.ddl$Psi$stratum==2 & ms4.ddl$Psi$tostratum == 3] = 0
ms4.ddl$Psi$fix[ms4.ddl$Psi$stratum==3 & ms4.ddl$Psi$tostratum == 4] = 0
ms4.ddl$Psi$fix[ms4.ddl$Psi$stratum==4 & ms4.ddl$Psi$tostratum == 3] = 0
head(ms4.ddl$Psi, 20)
#    id occ stratum tostratum time cohort age freq Time Cohort Age fix order
# 1   1   1       1         1    1      3   0  490    0      2   0   1     1
# 2   1   1       1         2    1      3   0  490    0      2   0  NA     2
# 3   1   1       1         3    1      3   0  490    0      2   0  NA     3
# 4   1   1       1         4    1      3   0  490    0      2   0   0     4
# 5   1   1       2         1    1      3   0  490    0      2   0  NA     5
# 6   1   1       2         2    1      3   0  490    0      2   0   1     6
# 7   1   1       2         3    1      3   0  490    0      2   0   0     7
# 8   1   1       2         4    1      3   0  490    0      2   0  NA     8
# 9   1   1       3         1    1      3   0  490    0      2   0  NA     9
# 10  1   1       3         2    1      3   0  490    0      2   0  NA    10
# 11  1   1       3         3    1      3   0  490    0      2   0   1    11
# 12  1   1       3         4    1      3   0  490    0      2   0   0    12
# 13  1   1       4         1    1      3   0  490    0      2   0  NA    13
# 14  1   1       4         2    1      3   0  490    0      2   0  NA    14
# 15  1   1       4         3    1      3   0  490    0      2   0   0    15
# 16  1   1       4         4    1      3   0  490    0      2   0   1    16
# 17  1   2       1         1    2      3   0  490    1      2   0   1    17
# 18  1   2       1         2    2      3   0  490    1      2   0  NA    18
# 19  1   2       1         3    2      3   0  490    1      2   0  NA    19
# 20  1   2       1         4    2      3   0  490    1      2   0   0    20

my.S.ddl <- ms4.ddl$S[ms4.ddl$S$id == 25,]
my.S.desmat = model.matrix(object = ~ -1 + mystratum, data=my.S.ddl)
my.S.desmat
#     mystratum1 mystratum2
# 289          1          0
# 290          0          1
# 291          1          0
# 292          0          1
# 293          1          0
# 294          0          1
# 295          1          0
# 296          0          1
# 297          1          0
# 298          0          1
# 299          1          0
# 300          0          1
# attr(,"assign")
# [1] 1 1
# attr(,"contrasts")
# attr(,"contrasts")$mystratum
# [1] "contr.treatment"

my.p.ddl <- ms4.ddl$p[ms4.ddl$p$id == 25,]
my.p.desmat = model.matrix(object = ~ -1 + mystratum, data=my.p.ddl)
my.p.desmat
#     mystratum1 mystratum2
# 289          1          0
# 290          1          0
# 291          0          1
# 292          0          1
# 293          1          0
# 294          1          0
# 295          0          1
# 296          0          1
# 297          1          0
# 298          1          0
# 299          0          1
# 300          0          1
# attr(,"assign")
# [1] 1 1
# attr(,"contrasts")
# attr(,"contrasts")$mystratum
# [1] "contr.treatment"

my.psi.ddl <- ms4.ddl$Psi[ms4.ddl$Psi$id == 25,]
my.psi.desmat = model.matrix(object = ~ -1 + stratum:tostratum, data=my.psi.ddl)
colnames(my.psi.desmat)
# [1] "stratum1:tostratum1" "stratum2:tostratum1" "stratum3:tostratum1" "stratum4:tostratum1"
#     "stratum1:tostratum2" "stratum2:tostratum2" "stratum3:tostratum2" "stratum4:tostratum2"
#     "stratum1:tostratum3" "stratum2:tostratum3" "stratum3:tostratum3" "stratum4:tostratum3"
#     "stratum1:tostratum4" "stratum2:tostratum4" "stratum3:tostratum4" "stratum4:tostratum4"

table(colSums(my.psi.desmat))
#
#  3
# 16


# fit model
S.stratum   = list(formula =~ -1 + mystratum)
p.stratum   = list(formula =~ -1 + mystratum)
Psi.stratum = list(formula =~ -1 + stratum:tostratum)

mod = crm(ms4.processed, ms4.ddl,
          model.parameters = list(S   = S.stratum,
                                  p   = p.stratum,
                                  Psi = Psi.stratum),
          hessian=TRUE, use.admb = FALSE, strata.labels=c("1","2","3","4"))
# compiling with args: '  -s ' ...
# compile output:
#   *** Parse: multistate.tpl xxglobal.tmp xxhtop.tmp
#          header.tmp xxalloc.tmp xxtopm.tmp
# 1 file(s) copied. tpl2cpp   multistate 
#   *** Compile: multistate.cpp g++ -c -std=c++17 -O2
#          -D_FILE_OFFSET_BITS=64 -DUSE_ADMB_CONTRIBS
#          -D_USE_MATH_DEFINES -I. -I"C:\ADMB-13.1\include"
#          -I"C:\ADMB-13.1\include\contrib" -o multistate.obj multistate.cpp
#   *** Linking: multistate.obj  g++ -static  -o multistate.exe multistate.obj
#   "C:\ADMB-13.1\lib\libadmb-contrib-mingw64-g++12.a"
#   Successfully built 'multistate.exe'.
# compile log:
#
#
# Elapsed time in minutes:  0.2932
#
# Warning message:
# In read_pars(fn) :
#   std file missing: some problem with fit,
#   but retrieving parameter estimates anyway

mod$results
#
# crm Model Summary
#
# Npar :  11
# -2lnL:  19318.68
# AIC  :  19340.68
#
# Beta
#                           Estimate           se           lcl           ucl
# S.mystratum1             16.3726193 5.914354e+03 -1.157576e+04  1.160851e+04
# S.mystratum2             17.5368972 1.394114e+04 -2.730709e+04  2.734217e+04
# p.mystratum1              0.8145821 7.018123e-02  6.770269e-01  9.521373e-01
# Psi.stratum2:tostratum1  -1.5906032 6.232242e-02 -1.712755e+00 -1.468451e+00
# Psi.stratum3:tostratum1  -3.3521341 5.624320e-01 -4.454501e+00 -2.249767e+00
# Psi.stratum4:tostratum1  -4.1004661 4.844378e-01 -5.049964e+00 -3.150968e+00
# Psi.stratum1:tostratum2  -0.8666090 4.539042e-02 -9.555742e-01 -7.776438e-01
# Psi.stratum3:tostratum2 -78.8418429 1.000000e+05 -1.960788e+05  1.959212e+05
# Psi.stratum4:tostratum2  -3.9666574 5.964215e-01 -5.135644e+00 -2.797671e+00
# Psi.stratum1:tostratum3  -0.6721677 7.699075e-02 -8.230696e-01 -5.212659e-01
# Psi.stratum2:tostratum4  -0.1440854 5.802741e-02 -2.578191e-01 -3.035166e-02


# from RMark
# S.sitetime.p.sitetime.psi.sitetime$results$beta
#                           estimate        se         lcl        ucl
# S:mystratum1             1.3864154 0.9894669  -0.5529398  3.3257706
# S:mystratum2             0.8474320 1.1155115  -1.3389706  3.0338345
# p:mystratum1             1.0984849 0.9659791  -0.7948342  2.9918039
# Psi:stratum2:tostratum1 -1.5939329 0.0640032  -1.7193791 -1.4684866
# Psi:stratum3:tostratum1 -0.8028675 4.1493275  -8.9355494  7.3298145
# Psi:stratum4:tostratum1 -1.7496550 3.6444092  -8.8926972  5.3933872
# Psi:stratum1:tostratum2 -0.8754691 0.0474696  -0.9685096 -0.7824286
# Psi:stratum3:tostratum2 -1.2884317 4.5637118 -10.2333070  7.6564435
# Psi:stratum4:tostratum2 -1.2902143 4.3111878  -9.7401426  7.1597140
# Psi:stratum1:tostratum3 -1.3863451 0.6027196  -2.5676756 -0.2050146
# Psi:stratum2:tostratum4 -1.0233519 0.6688294  -2.3342575  0.2875538


# from RMark
# S.sitetime.p.sitetime.psi.sitetime$results$real
#                            estimate        se          lcl       ucl fixed note
# S s1 g1 c1 a0 o1 t1       0.8000194 0.1583032 0.3651826000 0.9653024         
# S s2 g1 c1 a0 o1 t1       0.7000282 0.2342448 0.2076794000 0.9540795         
# p s1 g1 c1 a1 o1 t2       0.7499761 0.1811326 0.3111316000 0.9522025         
# p s3 g1 c1 a1 o1 t2       0.0000000 0.0000000 0.0000000000 0.0000000 Fixed   
# Psi s1 to2 g1 c1 a0 o1 t1 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c1 a1 o2 t2 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c1 a2 o3 t3 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c2 a0 o2 t2 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c2 a1 o3 t3 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c3 a0 o3 t3 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to3 g1 c1 a0 o1 t1 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c1 a1 o2 t2 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c1 a2 o3 t3 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c2 a0 o2 t2 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c2 a1 o3 t3 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c3 a0 o3 t3 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s2 to1 g1 c1 a0 o1 t1 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c1 a1 o2 t2 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c1 a2 o3 t3 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c2 a0 o2 t2 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c2 a1 o3 t3 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c3 a0 o3 t3 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to4 g1 c1 a0 o1 t1 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c1 a1 o2 t2 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c1 a2 o3 t3 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c2 a0 o2 t2 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c2 a1 o3 t3 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c3 a0 o3 t3 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s3 to1 g1 c1 a0 o1 t1 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c1 a1 o2 t2 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c1 a2 o3 t3 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c2 a0 o2 t2 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c2 a1 o3 t3 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c3 a0 o3 t3 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to2 g1 c1 a0 o1 t1 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c1 a1 o2 t2 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c1 a2 o3 t3 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c2 a0 o2 t2 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c2 a1 o3 t3 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c3 a0 o3 t3 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s4 to1 g1 c1 a0 o1 t1 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c1 a1 o2 t2 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c1 a2 o3 t3 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c2 a0 o2 t2 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c2 a1 o3 t3 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c3 a0 o3 t3 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to2 g1 c1 a0 o1 t1 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c1 a1 o2 t2 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c1 a2 o3 t3 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c2 a0 o2 t2 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c2 a1 o3 t3 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c3 a0 o3 t3 0.1899262 0.5808033 0.0001434280 0.9973972


Here is the RMark code containing results from RMark and MARK:

Code: Select all
#
# RMark code
#
library(RMark)

# Read the fake two observable and two unobservable state deterministic data
ms <- convert.inp("fake_2obs_2unobs_state_datai.inp")
head(ms)

# Process data
mstrata.processed = process.data(ms, begin.time = 1,
     model = "Multistrata")

# Create default design data
mstrata.ddl = make.design.data(mstrata.processed)

# examine data
head(mstrata.ddl$S)
head(mstrata.ddl$p)
head(mstrata.ddl$Psi)


# constrain survival probabilities:
# 1. equal in States 1 & 3 and in States 2 & 4
# 2. constant over time within a state.
mstrata.ddl$S$mystratum[mstrata.ddl$S$stratum %in% c(1,3)] <- 1
mstrata.ddl$S$mystratum[mstrata.ddl$S$stratum %in% c(2,4)] <- 2
mstrata.ddl$S$mystratum <- factor(mstrata.ddl$S$mystratum)
mstrata.ddl$S


mstrata.ddl$p$mytime1 <- NA
mstrata.ddl$p$mytime1[mstrata.ddl$p$time == 1] = 1
mstrata.ddl$p$mytime1[mstrata.ddl$p$time == 2] = 2
mstrata.ddl$p$mytime1[mstrata.ddl$p$time %in% c(3,4)] = 3
mstrata.ddl$p$mytime1 <- factor(mstrata.ddl$p$mytime1)

# Constrain detection probability
# 1. constant and equal in States 1 & 2
# 2. zero in States 3 & 4
mstrata.ddl$p$mystratum[mstrata.ddl$p$stratum %in% c(1,2)] <- 1
mstrata.ddl$p$mystratum[mstrata.ddl$p$stratum %in% c(3,4)] <- 2
mstrata.ddl$p$mystratum <- factor(mstrata.ddl$p$mystratum)
mstrata.ddl$p$fix <- NA
mstrata.ddl$p$fix[mstrata.ddl$p$mystratum == 2] <- 0
mstrata.ddl$p

# constrain transition probability = 0 for transitions that are not possible
mstrata.ddl$Psi$fix=NA
mstrata.ddl$Psi$fix[mstrata.ddl$Psi$stratum==1 & mstrata.ddl$Psi$tostratum == 4] = 0
mstrata.ddl$Psi$fix[mstrata.ddl$Psi$stratum==2 & mstrata.ddl$Psi$tostratum == 3] = 0
mstrata.ddl$Psi$fix[mstrata.ddl$Psi$stratum==3 & mstrata.ddl$Psi$tostratum == 4] = 0
mstrata.ddl$Psi$fix[mstrata.ddl$Psi$stratum==4 & mstrata.ddl$Psi$tostratum == 3] = 0

# survival probability
S.site.time = list(formula = ~ -1 + mystratum)

my.S.mstrata.ddl <- mstrata.ddl$S[mstrata.ddl$S$cohort == 1,]
my.S.desmat = model.matrix(object = ~ -1 + mystratum, data=my.S.mstrata.ddl)
my.S.desmat
#    mystratum1 mystratum2
# 1           1          0
# 2           1          0
# 3           1          0
# 7           0          1
# 8           0          1
# 9           0          1
# 13          1          0
# 14          1          0
# 15          1          0
# 19          0          1
# 20          0          1
# 21          0          1

# detection probability
p.site.time = list(formula = ~ -1 + mystratum)

my.p.mstrata.ddl <- mstrata.ddl$p[mstrata.ddl$p$cohort == 1,]
my.p.desmat = model.matrix(object = ~ -1 + mystratum, data=my.p.mstrata.ddl)
my.p.desmat
#    mystratum1 mystratum2
# 1           1          0
# 2           1          0
# 3           1          0
# 7           1          0
# 8           1          0
# 9           1          0
# 13          0          1
# 14          0          1
# 15          0          1
# 19          0          1
# 20          0          1
# 21          0          1

# transition probs
Psi.site.time = list(formula = ~ -1 + stratum:tostratum)


S.sitetime.p.sitetime.psi.sitetime = mark(mstrata.processed,
     mstrata.ddl,
     model.parameters = list(S = S.site.time,
     p = p.site.time, Psi = Psi.site.time),
     output=FALSE, silent=TRUE, delete=TRUE, mlogit0=TRUE)

# examine the design matrix for transition probabilites in the above model
# Just use Cohort 1
my.mstrata.ddl <- mstrata.ddl$Psi[mstrata.ddl$Psi$cohort == 1,]
my.desmat = model.matrix(object = ~ -1 + stratum:tostratum,
                         data=my.mstrata.ddl)

colnames(my.desmat)
table(colSums(my.desmat))


# examine output betas
S.sitetime.p.sitetime.psi.sitetime$results$beta
#                           estimate        se         lcl        ucl
# S:mystratum1             1.3864154 0.9894669  -0.5529398  3.3257706
# S:mystratum2             0.8474320 1.1155115  -1.3389706  3.0338345
# p:mystratum1             1.0984849 0.9659791  -0.7948342  2.9918039
# Psi:stratum2:tostratum1 -1.5939329 0.0640032  -1.7193791 -1.4684866
# Psi:stratum3:tostratum1 -0.8028675 4.1493275  -8.9355494  7.3298145
# Psi:stratum4:tostratum1 -1.7496550 3.6444092  -8.8926972  5.3933872
# Psi:stratum1:tostratum2 -0.8754691 0.0474696  -0.9685096 -0.7824286
# Psi:stratum3:tostratum2 -1.2884317 4.5637118 -10.2333070  7.6564435
# Psi:stratum4:tostratum2 -1.2902143 4.3111878  -9.7401426  7.1597140
# Psi:stratum1:tostratum3 -1.3863451 0.6027196  -2.5676756 -0.2050146
# Psi:stratum2:tostratum4 -1.0233519 0.6688294  -2.3342575  0.2875538

# from MARK
# PARM-SPECIFIC Link Function Parameters of {phi(g) p(.) phi(g) DM}
#     Parameter       Beta      Standard Error
#  -------------- ------------  --------------
#     1:          1.3861429       0.9905899     phi(state 1 & 3)
#     2:          0.8471931       1.1179774     phi(state 2 & 4)
#     3:          1.0987339       0.9700386     p(state 1 & 2)
#     4:          0.0000000       0.0000000     p(state 3 & 4)
#     5:         -0.8754678       0.0474701     Psi(1 to 2)
#     6:         -1.3862928       0.6049720     Psi(1 to 3)
#     7:          0.0000000       0.0000000     Psi(1 to 4)
#     8:         -1.5939344       0.0640055     Psi(2 to 1)
#     9:          0.0000000       0.0000000     Psi(2 to 3)
#    10:         -1.0233975       0.6684302     Psi(2 to 4)
#    11:         -0.8017422       4.1581095     Psi(3 to 1)
#    12:         -1.2872080       4.5739741     Psi(3 to 2)
#    13:          0.0000000       0.0000000     Psi(3 to 4)
#    14:         -1.7488127       3.6549281     Psi(4 to 1)
#    15:         -1.2891991       4.3241670     Psi(4 to 2)
#    16:          0.0000000       0.0000000     Psi(4 to 3)


# examine output reals
S.sitetime.p.sitetime.psi.sitetime$results$real
#                            estimate        se          lcl       ucl fixed note
# S s1 g1 c1 a0 o1 t1       0.8000194 0.1583032 0.3651826000 0.9653024         
# S s2 g1 c1 a0 o1 t1       0.7000282 0.2342448 0.2076794000 0.9540795         
# p s1 g1 c1 a1 o1 t2       0.7499761 0.1811326 0.3111316000 0.9522025         
# p s3 g1 c1 a1 o1 t2       0.0000000 0.0000000 0.0000000000 0.0000000 Fixed   
# Psi s1 to2 g1 c1 a0 o1 t1 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c1 a1 o2 t2 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c1 a2 o3 t3 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c2 a0 o2 t2 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c2 a1 o3 t3 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to2 g1 c3 a0 o3 t3 0.2500018 0.0239926 0.2059679000 0.2998943         
# Psi s1 to3 g1 c1 a0 o1 t1 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c1 a1 o2 t2 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c1 a2 o3 t3 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c2 a0 o2 t2 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c2 a1 o3 t3 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s1 to3 g1 c3 a0 o3 t3 0.1499935 0.0768010 0.0514014000 0.3649408         
# Psi s2 to1 g1 c1 a0 o1 t1 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c1 a1 o2 t2 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c1 a2 o3 t3 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c2 a0 o2 t2 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c2 a1 o3 t3 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to1 g1 c3 a0 o3 t3 0.1299990 0.0207751 0.0944049000 0.1763992         
# Psi s2 to4 g1 c1 a0 o1 t1 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c1 a1 o2 t2 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c1 a2 o3 t3 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c2 a0 o2 t2 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c2 a1 o3 t3 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s2 to4 g1 c3 a0 o3 t3 0.2300065 0.1183304 0.0746174000 0.5252997         
# Psi s3 to1 g1 c1 a0 o1 t1 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c1 a1 o2 t2 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c1 a2 o3 t3 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c2 a0 o2 t2 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c2 a1 o3 t3 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to1 g1 c3 a0 o3 t3 0.2599238 0.6095555 0.0007046141 0.9943161         
# Psi s3 to2 g1 c1 a0 o1 t1 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c1 a1 o2 t2 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c1 a2 o3 t3 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c2 a0 o2 t2 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c2 a1 o3 t3 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s3 to2 g1 c3 a0 o3 t3 0.1599441 0.4417536 0.0003026264 0.9917185         
# Psi s4 to1 g1 c1 a0 o1 t1 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c1 a1 o2 t2 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c1 a2 o3 t3 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c2 a0 o2 t2 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c2 a1 o3 t3 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to1 g1 c3 a0 o3 t3 0.1199644 0.2872790 0.0006576324 0.9657980         
# Psi s4 to2 g1 c1 a0 o1 t1 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c1 a1 o2 t2 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c1 a2 o3 t3 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c2 a0 o2 t2 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c2 a1 o3 t3 0.1899262 0.5808033 0.0001434280 0.9973972         
# Psi s4 to2 g1 c3 a0 o3 t3 0.1899262 0.5808033 0.0001434280 0.9973972

# from MARK
# Real Function Parameters of {phi(g) p(.) phi(g) DM}
#      Parameter            Estimate       Standard Error
#  --------------------  --------------  --------------
#      1:S 1:State 1       0.7999758       0.1585088
#      2:S 1:State 1       0.7999758       0.1585088
#      3:S 1:State 1       0.7999758       0.1585088
#      4:S 2:State 2       0.6999780       0.2347851
#      5:S 2:State 2       0.6999780       0.2347851
#      6:S 2:State 2       0.6999780       0.2347851
#      7:S 3:State 3       0.7999758       0.1585088
#      8:S 3:State 3       0.7999758       0.1585088
#      9:S 3:State 3       0.7999758       0.1585088
#     10:S 4:State 4       0.6999780       0.2347851
#     11:S 4:State 4       0.6999780       0.2347851
#     12:S 4:State 4       0.6999780       0.2347851
#     13:p 1:State 1       0.7500228       0.1818712
#     14:p 1:State 1       0.7500228       0.1818712
#     15:p 1:State 1       0.7500228       0.1818712
#     16:p 2:State 2       0.7500228       0.1818712
#     17:p 2:State 2       0.7500228       0.1818712
#     18:p 2:State 2       0.7500228       0.1818712
#     19:p 3:State 3       0.0000000       0.0000000   Fixed
#     20:p 3:State 3       0.0000000       0.0000000   Fixed
#     21:p 3:State 3       0.0000000       0.0000000   Fixed
#     22:p 4:State 4       0.0000000       0.0000000   Fixed
#     23:p 4:State 4       0.0000000       0.0000000   Fixed
#     24:p 4:State 4       0.0000000       0.0000000   Fixed
#     25:Psi 1 to 2        0.2500001       0.0240701
#     26:Psi 1 to 2        0.2500001       0.0240701
#     27:Psi 1 to 2        0.2500001       0.0240701
#     28:Psi 1 to 3        0.1500002       0.0770905
#     29:Psi 1 to 3        0.1500002       0.0770905
#     30:Psi 1 to 3        0.1500002       0.0770905
#     31:Psi 1 to 4        0.0000000       0.0000000   Fixed
#     32:Psi 1 to 4        0.0000000       0.0000000   Fixed
#     33:Psi 1 to 4        0.0000000       0.0000000   Fixed
#     34:Psi 2 to 1        0.1300002       0.0207629
#     35:Psi 2 to 1        0.1300002       0.0207629
#     36:Psi 2 to 1        0.1300002       0.0207629
#     37:Psi 2 to 3        0.0000000       0.0000000   Fixed
#     38:Psi 2 to 3        0.0000000       0.0000000   Fixed
#     39:Psi 2 to 3        0.0000000       0.0000000   Fixed
#     40:Psi 2 to 4        0.2299985       0.1182567
#     41:Psi 2 to 4        0.2299985       0.1182567
#     42:Psi 2 to 4        0.2299985       0.1182567
#     43:Psi 3 to 1        0.2600894       0.6108861
#     44:Psi 3 to 1        0.2600894       0.6108861
#     45:Psi 3 to 1        0.2600894       0.6108861
#     46:Psi 3 to 2        0.1600617       0.4428943
#     47:Psi 3 to 2        0.1600617       0.4428943
#     48:Psi 3 to 2        0.1600617       0.4428943
#     49:Psi 3 to 4        0.0000000       0.0000000   Fixed
#     50:Psi 3 to 4        0.0000000       0.0000000   Fixed
#     51:Psi 3 to 4        0.0000000       0.0000000   Fixed
#     52:Psi 4 to 1        0.1200302       0.2881501
#     53:Psi 4 to 1        0.1200302       0.2881501
#     54:Psi 4 to 1        0.1200302       0.2881501
#     55:Psi 4 to 2        0.1900633       0.5828210
#     56:Psi 4 to 2        0.1900633       0.5828210
#     57:Psi 4 to 2        0.1900633       0.5828210
#     58:Psi 4 to 3        0.0000000       0.0000000   Fixed
#     59:Psi 4 to 3        0.0000000       0.0000000   Fixed
#     60:Psi 4 to 3        0.0000000       0.0000000   Fixed

# Examine the estimated transition matrix
Psilist = get.real(S.sitetime.p.sitetime.psi.sitetime,"Psi",vcv=TRUE)
Psivalues = Psilist$estimates

# For the 1st-to-2nd capture occasion
TransitionMatrix(Psivalues[Psivalues$time==1,])



Here is the full constrained data set:

Code: Select all
1100 136.2051  ;
1101 23.5197  ;
1102 16.6752  ;
1110 63.504  ;
1111 46.656  ;
1112 19.44  ;
1120 32.1705  ;
1121 3.6855  ;
1122 18.144  ;
1200 76.9003125  ;
1201 3.89655  ;
1202 8.5656375  ;
1210 5.016375  ;
1211 3.6855  ;
1212 1.535625  ;
1220 30.0258  ;
1221 3.4398  ;
1222 16.9344  ;
1000 342.0290575  ;
1001 20.33331  ;
1002 15.9851325  ;
1010 32.012925  ;
1011 23.5197  ;
1012 9.799875  ;
1020 27.59514  ;
1021 3.16134  ;
1022 15.56352  ;
2100 25.822216875  ;
2101 4.458943125  ;
2102 3.16134  ;
2110 12.0393  ;
2111 8.8452  ;
2112 3.6855  ;
2120 6.098990625  ;
2121 0.698709375  ;
2122 3.4398  ;
2200 172.2567  ;
2201 8.728272  ;
2202 19.187028  ;
2210 11.23668  ;
2211 8.25552  ;
2212 3.4398  ;
2220 67.257792  ;
2221 7.705152  ;
2222 37.933056  ;
2000 484.2792660625  ;
2001 10.8773266875  ;
2002 17.51215725  ;
2010 12.72873  ;
2011 9.35172  ;
2012 3.89655  ;
2020 34.0198569375  ;
2021 3.8973650625  ;
2022 19.187028  ;
0100 378.3475  ;
0101 65.3325  ;
0102 46.32  ;
0110 176.4  ;
0111 129.6  ;
0112 54  ;
0120 89.3625  ;
0121 10.2375  ;
0122 50.4  ;
0200 512.66875  ;
0201 25.977  ;
0202 57.10425  ;
0210 33.4425  ;
0211 24.57  ;
0212 10.2375  ;
0220 200.172  ;
0221 22.932  ;
0222 112.896  ;
0010 490  ;
0011 360  ;
0012 150  ;
0020 595.75  ;
0021 68.25  ;
0022 336  ;
markmiller
 
Posts: 49
Joined: Fri Nov 08, 2013 6:23 pm

Re: multistate model with unobservable states

Postby markmiller » Tue Apr 04, 2023 6:00 am

I have been able to get the R package 'marked' to return good estimates matching those from MARK when using two observable states (States 1 & 2) and one unobservable state (State 3). Here is the full R code and full deterministic data set. I also include the beta estimates from 'marked' and 'MARK'. I can keep trying to get good estimates with two observable states and two unobservable states:

Here is the R code:

Code: Select all
library(marked)
library(R2admb)
library(TMB)
library(ggplot2)

#####@@@@@

# mwm: Use ?crm to locate the 'prepare_admb' function.

prepare_admb=function()
{
  Sys.setenv(PATH = paste("C:/ADMB-13.1/bin;C:/ADMB-13.1/utilities;C:/ADMB-13.1/utilities/mingw64/bin",
        Sys.getenv("PATH"), sep = ";"))
    Sys.setenv(ADMB_HOME = "C:/ADMB-13.1")
    invisible()
}
prepare_admb()

#####@@@@@

three.raw.pre <- read.table("fake_2obs_1unobs_state_in_marked_Apr4_2023.inp",
                header = FALSE, stringsAsFactors = FALSE, na.strings = "NA",
                colClasses = c("character", "numeric", "character"))
head(three.raw.pre)

colnames(three.raw.pre) <- c("ch", "freq", "the.end")
three.raw <- three.raw.pre[, !colnames(three.raw.pre) %in% "the.end"]
head(three.raw)
tail(three.raw)

ms3.processed = process.data(three.raw, model = "Mscjs", strata.labels=c("1","2","3"))
ms3.ddl = make.design.data(ms3.processed)

# constrain survival probabilities:
# 1. constant and equal in States 1 & 3
# 2. constant over time within a state.
ms3.ddl$S$mystratum <- NA
ms3.ddl$S$mystratum[ms3.ddl$S$stratum %in% c(1,3)] <- 1
ms3.ddl$S$mystratum[ms3.ddl$S$stratum %in% c(2)] <- 2
ms3.ddl$S$mystratum <- factor(ms3.ddl$S$mystratum)
head(ms3.ddl$S)

# Constrain detection probability
# 1. constant and equal in States 1 & 2
# 2. zero in States 3 & 4
ms3.ddl$p$mystratum <- NA
ms3.ddl$p$mystratum[ms3.ddl$p$stratum %in% c(1,2)] <- 1
ms3.ddl$p$mystratum[ms3.ddl$p$stratum %in% c(3)] <- 2
ms3.ddl$p$mystratum <- factor(ms3.ddl$p$mystratum)
ms3.ddl$p$fix <- NA
ms3.ddl$p$fix[ms3.ddl$p$mystratum == 2] <- 0

# fit model
S.stratum.3   = list(formula =~ -1 + mystratum)
p.dot.3       = list(formula =~ -1 + mystratum)
Psi.stratum.3 = list(formula =~ -1 + stratum:tostratum)

mod3a = crm(ms3.processed, ms3.ddl,
          model.parameters = list(S   = S.stratum.3,
                                  p   = p.dot.3,
                                  Psi = Psi.stratum.3),
          hessian=TRUE, use.admb = FALSE, strata.labels=c("1","2","3"))
mod3a$results
#
# crm Model Summary
#
# Npar :  9
# -2lnL:  19259.54
# AIC  :  19277.54
#
# Beta
#                           Estimate         se        lcl        ucl
# S.mystratum1             1.3862959 0.31764883  0.7637042  2.0088876
# S.mystratum2             0.4054657 0.11608349  0.1779421  0.6329893
# p.mystratum1             0.8472956 0.31948116  0.2211125  1.4734787
# Psi.stratum2:tostratum1 -0.6931472 0.05348020 -0.7979684 -0.5883260
# Psi.stratum3:tostratum1 -0.9163019 1.73460866 -4.3161349  2.4835311
# Psi.stratum1:tostratum2 -1.1786548 0.05192712 -1.2804320 -1.0768777
# Psi.stratum3:tostratum2 -0.5108313 1.18669908 -2.8367615  1.8150989
# Psi.stratum1:tostratum3 -1.4663395 0.49436485 -2.4352946 -0.4973844
# Psi.stratum2:tostratum3 -1.7917642 0.82312392 -3.4050871 -0.1784413
#
#
# From MARK:
#
#      PARM-SPECIFIC Link Function Parameters of {constant2}
#     Parameter                  Beta         Standard Error
#  -------------------------  --------------  --------------
#     1:S 1:State 1            1.3862733       0.3175979
#     2:S 2:State 2            0.4054573       0.1160707
#     3:p 1:State 1            0.8473072       0.3194839
#     4:p 3:State 3            0.0000000       0.0000000
#     5:Psi 1 to 2            -1.1786556       0.0519271
#     6:Psi 1 to 3            -1.4663452       0.4943975
#     7:Psi 2 to 1            -0.6931470       0.0534803
#     8:Psi 2 to 3            -1.7917661       0.8231651
#     9:Psi 3 to 1            -0.9161791       1.7344097
#    10:Psi 3 to 2            -0.5107389       1.1865752
#


Here is the data set:

Code: Select all
1100 144.872  ;
1101 27.763008  ;
1102 18.100992  ;
1110 69.427904  ;
1111 48.228544  ;
1112 14.839552  ;
1120 25.357696  ;
1121 5.136768  ;
1122 10.273536  ;
1200 60.331264  ;
1201 4.478208  ;
1202 4.854528  ;
1210 7.394688  ;
1211 5.136768  ;
1212 1.580544  ;
1220 17.555328  ;
1221 3.556224  ;
1222 7.112448  ;
1000 354.896576  ;
1001 22.378944  ;
1002 20.72448  ;
1010 39.966528  ;
1011 27.763008  ;
1012 8.542464  ;
1020 30.930816  ;
1021 6.265728  ;
1022 12.531456  ;
2100 50.148  ;
2101 9.610272  ;
2102 6.265728  ;
2110 24.032736  ;
2111 16.694496  ;
2112 5.136768  ;
2120 8.777664  ;
2121 1.778112  ;
2122 3.556224  ;
2200 135.745344  ;
2201 10.075968  ;
2202 10.922688  ;
2210 16.638048  ;
2211 11.557728  ;
2212 3.556224  ;
2220 39.499488  ;
2221 8.001504  ;
2222 16.003008  ;
2000 513.144736  ;
2001 12.7176  ;
2002 12.809664  ;
2010 20.951616  ;
2011 14.554176  ;
2012 4.478208  ;
2020 26.959968  ;
2021 5.461344  ;
2022 10.922688  ;
0100 398  ;
0101 76.272  ;
0102 49.728  ;
0110 190.736  ;
0111 132.496  ;
0112 40.768  ;
0120 69.664  ;
0121 14.112  ;
0122 28.224  ;
0200 538.672  ;
0201 39.984  ;
0202 43.344  ;
0210 66.024  ;
0211 45.864  ;
0212 14.112  ;
0220 156.744  ;
0221 31.752  ;
0222 63.504  ;
0010 524  ;
0011 364  ;
0012 112  ;
0020 622  ;
0021 126  ;
0022 252  ;
markmiller
 
Posts: 49
Joined: Fri Nov 08, 2013 6:23 pm

Re: multistate model with unobservable states

Postby markmiller » Thu Apr 06, 2023 3:36 am

Jeff Laake helped me get the model with two observable states and two unobservable states to run using his R package 'markedTMB'. That package returned essentially the same estimates as RMark.
markmiller
 
Posts: 49
Joined: Fri Nov 08, 2013 6:23 pm

Re: multistate model with unobservable states

Postby jlaake » Fri Jun 09, 2023 2:16 pm

The only difference between marked and markedTMB for this analysis is that the latter uses method="nlminb" by default whereas marked uses method="BFGS" by default. I have found that nlminb is a more reliable optimization function. When we set method="nlminb" with use.tmb=TRUE then the correct estimates were produced with marked. Note that Mark was not able to get ADMB to produce the correct estimates even using nlminb, if I understood him correctly. ADMB is no longer going to be developed or supported by the ADMB foundation and I will pull it from a future version of marked. I recommend using use.tmb=TRUE for any applications you use currently.
jlaake
 
Posts: 1410
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA


Return to analysis help

Who is online

Users browsing this forum: No registered users and 1 guest

cron