I have a mix of PIT-tagged and telemetered individuals with 22 capture occasions. While all individuals were PIT-tagged for the duration of the study, some were also telemetered at varying times for varying durations. I am using Burnham’s joint recapture/recovery model and would like to fix p=1 for the telemetered individuals during the intervals in which they were telemetered. It seems like a time-varying individual covariate (say “cov”) would be the way to accomplish this, where “cov” was 1 if the individual had a functioning transmitter and 0 otherwise. For starters, I tried fitting a model where p was modeled as a function of “Telem” [i.e., S(.)p(cov)r(.)F(.)].
However, when I fit this model I had 21 (k-1) real parameter estimates for p. This wasn’t what I expected so I wanted to see if my results are indeed legitimate or I have set my data up incorrectly.
My data consists of a 23-column data frame, the first column containing the capture histories (44 values in LD format for 22 occasions) and 22 columns denoting the “cov” status (0/1) for each individual on that occasion. The 22 “cov” columns are labeled “cov1”, “cov2”, …, “cov22”.
- Code: Select all
> str(CH)
'data.frame': 99 obs. of 23 variables:
$ ch : chr "10000000000000000000000000000000000000000000" "10000000000000000000000000000000000000000000" "11000000000000000000000000000000000000000000" "10000000000000000000000000000000000000000000" ...
$ cov1 : num 0 0 0 0 0 0 0 0 0 0 ...
$ cov2 : num 0 0 0 0 0 0 0 0 0 0 ...
$ cov3 : num 0 0 0 0 0 0 0 0 0 0 ...
$ cov4 : num 0 0 0 0 0 0 0 0 0 0 ...
$ cov5 : num 0 0 0 0 0 0 0 0 0 0 ...
$ cov6 : num 0 0 0 0 0 0 0 0 0 0 ...
$ cov7 : num 0 0 0 0 0 0 0 0 0 0 ...
$ cov8 : num 0 0 0 0 0 0 0 0 0 0 ...
$ cov9 : num 0 0 0 0 0 0 0 0 0 0 ...
$ cov10: num 0 0 0 0 0 0 0 0 0 0 ...
$ cov11: num 0 0 0 0 0 0 0 0 0 0 ...
$ cov12: num 0 0 0 0 0 0 0 0 0 0 ...
$ cov13: num 0 0 0 0 1 0 0 0 0 0 ...
$ cov14: num 0 0 0 0 1 0 0 0 0 0 ...
$ cov15: num 0 0 0 0 1 0 0 0 0 0 ...
$ cov16: num 0 0 0 0 1 0 0 0 0 0 ...
$ cov17: num 0 0 0 0 1 0 0 0 0 0 ...
$ cov18: num 0 0 0 0 1 0 0 0 0 0 ...
$ cov19: num 0 0 0 0 0 1 0 0 0 0 ...
$ cov20: num 0 0 0 0 0 1 0 0 0 0 ...
$ cov21: num 0 0 0 0 0 0 0 0 0 0 ...
$ cov22: num 0 0 0 0 0 0 0 0 0 0 ...
t1_process <- process.data(data = CH,model="Burnham")
t1.ddl <- make.design.data(t1_process)
> head(t1_process$data)
ch cov1 cov2 cov3 cov4 cov5 cov6 cov7 cov8 cov9 cov10 cov11 cov12
1 10000000000000000000000000000000000000000000 0 0 0 0 0 0 0 0 0 0 0 0
2 10000000000000000000000000000000000000000000 0 0 0 0 0 0 0 0 0 0 0 0
3 11000000000000000000000000000000000000000000 0 0 0 0 0 0 0 0 0 0 0 0
4 10000000000000000000000000000000000000000000 0 0 0 0 0 0 0 0 0 0 0 0
5 10000010101010100000100010101010101000000000 0 0 0 0 0 0 0 0 0 0 0 0
6 10101010101010000000000000000000000010100100 0 0 0 0 0 0 0 0 0 0 0 0
cov13 cov14 cov15 cov16 cov17 cov18 cov19 cov20 cov21 cov22 freq
1 0 0 0 0 0 0 0 0 0 0 1
2 0 0 0 0 0 0 0 0 0 0 1
3 0 0 0 0 0 0 0 0 0 0 1
4 0 0 0 0 0 0 0 0 0 0 1
5 1 1 1 1 1 1 0 0 0 0 1
6 0 0 0 0 0 0 1 1 0 0 1
t1 <- mark(data = t1_process,
ddl = t1.ddl,
model.parameters = list(S = list(formula=~1), p = list(formula=~cov),
F = list(formula=~1,fixed=1), r = list(formula=~1)),
model = "Burnham")
> t1$results$beta
estimate se lcl ucl
S:(Intercept) 1.830187 0.1464575 1.5431304 2.1172437
p:(Intercept) -0.438305 0.1407934 -0.7142601 -0.1623500
p:cov 20.501717 2454.5176000 -4790.3529000 4831.3564000
r:(Intercept) -1.023930 0.2801911 -1.5731046 -0.4747556
> t1$results$real
estimate se lcl ucl fixed note
S g1 c1 a0 t1 0.8617840 0.0174449 8.239193e-01 0.8925679
p g1 c1 a1 t2 0.4424519 6.1162633 6.236393e-22 1.0000000
p g1 c1 a2 t3 0.5456126 18.4401140 5.839770e-64 1.0000000
p g1 c1 a3 t4 0.3921449 0.0335605 3.286582e-01 0.4595014
p g1 c1 a4 t5 0.4424519 6.1162633 6.236393e-22 1.0000000
p g1 c1 a5 t6 0.3921449 0.0335605 3.286582e-01 0.4595014
p g1 c1 a6 t7 0.4424519 6.1162633 6.236393e-22 1.0000000
p g1 c1 a7 t8 0.4424519 6.1162633 6.236393e-22 1.0000000
p g1 c1 a8 t9 0.4424519 6.1162633 6.236393e-22 1.0000000
p g1 c1 a9 t10 0.4424519 6.1162633 6.236393e-22 1.0000000
p g1 c1 a10 t11 0.4424519 6.1162633 6.236393e-22 1.0000000
p g1 c1 a11 t12 0.4424519 6.1162633 6.236393e-22 1.0000000
p g1 c1 a12 t13 0.6908772 31.7697870 5.288171e-127 1.0000000
p g1 c1 a13 t14 0.6908772 31.7697870 5.288171e-127 1.0000000
p g1 c1 a14 t15 0.7332767 33.9436050 5.116009e-148 1.0000000
p g1 c1 a15 t16 0.6908772 31.7697870 5.288171e-127 1.0000000
p g1 c1 a16 t17 0.7717817 34.9354160 4.949429e-169 1.0000000
p g1 c1 a17 t18 0.8061972 34.8638080 4.788258e-190 1.0000000
p g1 c1 a18 t19 0.9561522 17.6707520 2.150015e-303 1.0000000
p g1 c1 a19 t20 0.9561522 17.6707520 2.150015e-303 1.0000000
p g1 c1 a20 t21 0.9561522 17.6707520 2.150015e-303 1.0000000
p g1 c1 a21 t22 0.9640592 15.4630380 2.644714e-303 1.0000000
r g1 c1 a0 t1 0.2642626 0.0544769 1.717743e-01 0.3834913
F g1 c1 a0 t1 1.0000000 0.0000000 1.000000e+00 1.0000000 Fixed
> t1$design.matrix
S:(Intercept) p:(Intercept) p:cov r:(Intercept)
S g1 c1 a0 t1 "1" "0" "0" "0"
p g1 c1 a1 t2 "0" "1" "cov2" "0"
p g1 c1 a2 t3 "0" "1" "cov3" "0"
p g1 c1 a3 t4 "0" "1" "cov4" "0"
p g1 c1 a4 t5 "0" "1" "cov5" "0"
p g1 c1 a5 t6 "0" "1" "cov6" "0"
p g1 c1 a6 t7 "0" "1" "cov7" "0"
p g1 c1 a7 t8 "0" "1" "cov8" "0"
p g1 c1 a8 t9 "0" "1" "cov9" "0"
p g1 c1 a9 t10 "0" "1" "cov10" "0"
p g1 c1 a10 t11 "0" "1" "cov11" "0"
p g1 c1 a11 t12 "0" "1" "cov12" "0"
p g1 c1 a12 t13 "0" "1" "cov13" "0"
p g1 c1 a13 t14 "0" "1" "cov14" "0"
p g1 c1 a14 t15 "0" "1" "cov15" "0"
p g1 c1 a15 t16 "0" "1" "cov16" "0"
p g1 c1 a16 t17 "0" "1" "cov17" "0"
p g1 c1 a17 t18 "0" "1" "cov18" "0"
p g1 c1 a18 t19 "0" "1" "cov19" "0"
p g1 c1 a19 t20 "0" "1" "cov20" "0"
p g1 c1 a20 t21 "0" "1" "cov21" "0"
p g1 c1 a21 t22 "0" "1" "cov22" "0"
r g1 c1 a0 t1 "0" "0" "0" "1"
F g1 c1 a0 t1 "0" "0" "0" "0"
It seems (and I would like) to have two real parameter estimates, one when “cov”=0 and the other when “cov”=1 but I haven’t been able to do this, even though I feel like this must be a simple fix. Any guidance would be much appreciated!
Also, if I do get this to work correctly, would the following syntax be used to fix p=1 when “cov=1”? It threw an error message when I tried it with my current set up, again suggesting something is wrong with my data formatting(?).
- Code: Select all
> t1 <- mark(data = t1_process,
+ ddl = t1.ddl,
+ model.parameters = list(S = list(formula=~1),
+ p = list(formula=~cov,fixed=list(cov=1,value=1)),
+ F = list(formula=~1,fixed=1), r = list(formula=~1)),
+ model = "Burnham")
Error in make.mark.model(data.proc, title = title, parameters = model.parameters, :
Unrecognized structure for fixed parameters = 1
Unrecognized structure for fixed parameters = 1
Error in mark(data = t1_process, ddl = t1.ddl, model.parameters = list(S = list(formula = ~1), :
Misspecification of model or internal error in code
Thanks so much,
Javan