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 ;