# Analysis of killdeer data illustrating basic RMark features
# Fit a linear trend to the data
# 2019-05-01 CJS Code produced.
# This is the killdeer data that ships with RMark
# which has been saved as a *.csv file
library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
library(readxl)
library(RMark)
## This is RMark 2.2.6
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
# The dataframe must contain the following fields with the following names
#
# FirstFound: day the nest was first found
# LastPresent: last day that a chick was present in the nest
# LastChecked: last day the nest was checked
# Fate: fate of the nest; 0=hatch an
# Freq: number of nests with this data
#
# In this example, the multiple visits to a nest have been collapsed
# to a single record for each nest.
# In more complex examples, you may have multple records per nest
# as shown in the mallard example.
#
killdata <- readxl::read_excel(file.path("..","Killdeer.xlsx"),
sheet="killdeer")
head(killdata)
## # A tibble: 6 x 6
## id FirstFound LastPresent LastChecked Fate Freq
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 /*A*/ 1 9 9 0 1
## 2 /*B*/ 5 5 9 1 1
## 3 /*C*/ 5 40 40 0 1
## 4 /*D*/ 9 32 32 0 1
## 5 /*E*/ 7 8 8 0 1
## 6 /*F*/ 3 15 15 0 1
# what are the parameters of the model
# There is only one parameter, the daily survival probality (S)
setup.parameters("Nest", check=TRUE)
## [1] "S"
# 1. Process the data.
# The nocc variable is the data at which hatching occurs
kill.proc <- process.data(killdata, model="Nest", nocc=40)
kill.proc
## $data
## # A tibble: 18 x 6
## id FirstFound LastPresent LastChecked Fate freq
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 /*A*/ 1 9 9 0 1
## 2 /*B*/ 5 5 9 1 1
## 3 /*C*/ 5 40 40 0 1
## 4 /*D*/ 9 32 32 0 1
## 5 /*E*/ 7 8 8 0 1
## 6 /*F*/ 3 15 15 0 1
## 7 /*G*/ 8 32 32 0 1
## 8 /*H*/ 14 14 16 1 1
## 9 /*I*/ 8 14 14 0 1
## 10 /*J*/ 13 14 14 0 1
## 11 /*K*/ 14 33 33 0 1
## 12 /*L*/ 15 37 37 0 1
## 13 /*M*/ 16 37 40 1 1
## 14 /*N*/ 16 28 32 1 1
## 15 /*O*/ 16 17 17 0 1
## 16 /*P*/ 21 28 33 1 1
## 17 /*Q*/ 23 33 34 1 1
## 18 /*R*/ 27 29 29 0 1
##
## $model
## [1] "Nest"
##
## $mixtures
## [1] 1
##
## $freq
## group1
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1
## 7 1
## 8 1
## 9 1
## 10 1
## 11 1
## 12 1
## 13 1
## 14 1
## 15 1
## 16 1
## 17 1
## 18 1
##
## $nocc
## [1] 40
##
## $nocc.secondary
## NULL
##
## $time.intervals
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1
##
## $begin.time
## [1] 1
##
## $age.unit
## [1] 1
##
## $initial.ages
## [1] 0
##
## $group.covariates
## NULL
##
## $nstrata
## [1] 1
##
## $strata.labels
## [1] ""
##
## $counts
## NULL
##
## $reverse
## [1] FALSE
##
## $areas
## NULL
##
## $events
## NULL
# 2. Examine and/or modify the ddl. (Not done here)
kill.ddl <- make.design.data(kill.proc)
str(kill.ddl)
## List of 2
## $ S :'data.frame': 39 obs. of 7 variables:
## ..$ par.index : int [1:39] 1 2 3 4 5 6 7 8 9 10 ...
## ..$ model.index: num [1:39] 1 2 3 4 5 6 7 8 9 10 ...
## ..$ group : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## ..$ age : Factor w/ 39 levels "0","1","2","3",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..$ time : Factor w/ 39 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..$ Age : num [1:39] 0 1 2 3 4 5 6 7 8 9 ...
## ..$ Time : num [1:39] 0 1 2 3 4 5 6 7 8 9 ...
## $ pimtypes:List of 1
## ..$ S:List of 1
## .. ..$ pim.type: chr "all"
kill.ddl
## $S
## par.index model.index group age time Age Time
## 1 1 1 1 0 1 0 0
## 2 2 2 1 1 2 1 1
## 3 3 3 1 2 3 2 2
## 4 4 4 1 3 4 3 3
## 5 5 5 1 4 5 4 4
## 6 6 6 1 5 6 5 5
## 7 7 7 1 6 7 6 6
## 8 8 8 1 7 8 7 7
## 9 9 9 1 8 9 8 8
## 10 10 10 1 9 10 9 9
## 11 11 11 1 10 11 10 10
## 12 12 12 1 11 12 11 11
## 13 13 13 1 12 13 12 12
## 14 14 14 1 13 14 13 13
## 15 15 15 1 14 15 14 14
## 16 16 16 1 15 16 15 15
## 17 17 17 1 16 17 16 16
## 18 18 18 1 17 18 17 17
## 19 19 19 1 18 19 18 18
## 20 20 20 1 19 20 19 19
## 21 21 21 1 20 21 20 20
## 22 22 22 1 21 22 21 21
## 23 23 23 1 22 23 22 22
## 24 24 24 1 23 24 23 23
## 25 25 25 1 24 25 24 24
## 26 26 26 1 25 26 25 25
## 27 27 27 1 26 27 26 26
## 28 28 28 1 27 28 27 27
## 29 29 29 1 28 29 28 28
## 30 30 30 1 29 30 29 29
## 31 31 31 1 30 31 30 30
## 32 32 32 1 31 32 31 31
## 33 33 33 1 32 33 32 32
## 34 34 34 1 33 34 33 33
## 35 35 35 1 34 35 34 34
## 36 36 36 1 35 36 35 35
## 37 37 37 1 36 37 36 36
## 38 38 38 1 37 38 37 37
## 39 39 39 1 38 39 38 38
##
## $pimtypes
## $pimtypes$S
## $pimtypes$S$pim.type
## [1] "all"
# 3. Fit a particular model
# This is a model with S linear over time.
# Notice the use of Time vs time.
mod.res <- RMark::mark(kill.proc, ddl=kill.ddl,
model="Nest",
model.parameters=list(
S =list(formula=~Time)
)
)
## Warning: Unknown or uninitialised column: 'ch'.
## Warning: Unknown or uninitialised column: 'AgeDay1'.
##
## Output summary for Nest model
## Name : S(~Time)
##
## Npar : 2
## -2lnL: 40.97847
## AICc : 45.03644
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 4.9372980 1.3402857 2.3103379 7.5642581
## S:Time -0.0622855 0.0526986 -0.1655748 0.0410037
##
##
## Real Parameter S
## 1 2 3 4 5 6 7
## 0.9928771 0.9924229 0.9919398 0.9914263 0.9908803 0.9902999 0.9896829
## 8 9 10 11 12 13 14
## 0.9890272 0.9883302 0.9875895 0.9868025 0.9859663 0.9850778 0.984134
## 15 16 17 18 19 20 21
## 0.9831316 0.982067 0.9809365 0.9797361 0.9784619 0.9771094 0.9756741
## 22 23 24 25 26 27 28
## 0.9741512 0.9725356 0.9708221 0.969005 0.9670787 0.9650369 0.9628734
## 29 30 31 32 33 34 35
## 0.9605815 0.9581542 0.9555844 0.9528645 0.9499868 0.9469432 0.9437253
## 36 37 38 39
## 0.9403246 0.9367321 0.9329388 0.9289353
summary(mod.res)
## Output summary for Nest model
## Name : S(~Time)
##
## Npar : 2
## -2lnL: 40.97847
## AICc : 45.03644
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 4.9372980 1.3402857 2.3103379 7.5642581
## S:Time -0.0622855 0.0526986 -0.1655748 0.0410037
##
##
## Real Parameter S
## 1 2 3 4 5 6 7
## 0.9928771 0.9924229 0.9919398 0.9914263 0.9908803 0.9902999 0.9896829
## 8 9 10 11 12 13 14
## 0.9890272 0.9883302 0.9875895 0.9868025 0.9859663 0.9850778 0.984134
## 15 16 17 18 19 20 21
## 0.9831316 0.982067 0.9809365 0.9797361 0.9784619 0.9771094 0.9756741
## 22 23 24 25 26 27 28
## 0.9741512 0.9725356 0.9708221 0.969005 0.9670787 0.9650369 0.9628734
## 29 30 31 32 33 34 35
## 0.9605815 0.9581542 0.9555844 0.9528645 0.9499868 0.9469432 0.9437253
## 36 37 38 39
## 0.9403246 0.9367321 0.9329388 0.9289353
# Look the objects returned in more details
names(mod.res)
## [1] "data" "model" "title"
## [4] "model.name" "links" "mixtures"
## [7] "call" "parameters" "time.intervals"
## [10] "number.of.groups" "group.labels" "nocc"
## [13] "begin.time" "covariates" "fixed"
## [16] "design.matrix" "pims" "design.data"
## [19] "strata.labels" "mlogit.list" "profile.int"
## [22] "simplify" "model.parameters" "results"
## [25] "output"
names(mod.res$results)
## [1] "lnl" "deviance" "deviance.df"
## [4] "npar" "n" "AICc"
## [7] "beta" "real" "beta.vcv"
## [10] "derived" "derived.vcv" "covariate.values"
## [13] "singular" "real.vcv"
# look at estimates on beta and original scale
mod.res$results$beta # on the logit scale
## estimate se lcl ucl
## S:(Intercept) 4.9372980 1.3402857 2.3103379 7.5642581
## S:Time -0.0622855 0.0526986 -0.1655748 0.0410037
mod.res$results$real# on the regular 0-1 scale for each site
## estimate se lcl ucl fixed note
## S g1 a0 t1 0.9928771 0.0094787 0.9097296 0.9994816
## S g1 a1 t2 0.9924229 0.0097027 0.9126190 0.9993915
## S g1 a2 t3 0.9919398 0.0099182 0.9153905 0.9992862
## S g1 a3 t4 0.9914263 0.0101238 0.9180448 0.9991630
## S g1 a4 t5 0.9908803 0.0103177 0.9205823 0.9990191
## S g1 a5 t6 0.9902999 0.0104980 0.9230028 0.9988512
## S g1 a6 t7 0.9896829 0.0106632 0.9253056 0.9986556
## S g1 a7 t8 0.9890272 0.0108112 0.9274892 0.9984280
## S g1 a8 t9 0.9883302 0.0109404 0.9295514 0.9981638
## S g1 a9 t10 0.9875895 0.0110492 0.9314887 0.9978576
## S g1 a10 t11 0.9868025 0.0111361 0.9332964 0.9975036
## S g1 a11 t12 0.9859663 0.0112002 0.9349680 0.9970958
## S g1 a12 t13 0.9850778 0.0112410 0.9364948 0.9966275
## S g1 a13 t14 0.9841340 0.0112589 0.9378652 0.9960922
## S g1 a14 t15 0.9831316 0.0112554 0.9390639 0.9954837
## S g1 a15 t16 0.9820670 0.0112338 0.9400711 0.9947966
## S g1 a16 t17 0.9809365 0.0111994 0.9408604 0.9940273
## S g1 a17 t18 0.9797361 0.0111608 0.9413980 0.9931748
## S g1 a18 t19 0.9784619 0.0111303 0.9416397 0.9922427
## S g1 a19 t20 0.9771094 0.0111251 0.9415288 0.9912401
## S g1 a20 t21 0.9756741 0.0111683 0.9409938 0.9901840
## S g1 a21 t22 0.9741512 0.0112892 0.9399461 0.9890999
## S g1 a22 t23 0.9725356 0.0115233 0.9382794 0.9880216
## S g1 a23 t24 0.9708221 0.0119108 0.9358712 0.9869892
## S g1 a24 t25 0.9690050 0.0124937 0.9325866 0.9860437
## S g1 a25 t26 0.9670787 0.0133122 0.9282842 0.9852215
## S g1 a26 t27 0.9650369 0.0144012 0.9228208 0.9845480
## S g1 a27 t28 0.9628734 0.0157890 0.9160546 0.9840351
## S g1 a28 t29 0.9605815 0.0174969 0.9078453 0.9836815
## S g1 a29 t30 0.9581542 0.0195412 0.8980516 0.9834759
## S g1 a30 t31 0.9555844 0.0219351 0.8865300 0.9834012
## S g1 a31 t32 0.9528645 0.0246908 0.8731342 0.9834378
## S g1 a32 t33 0.9499868 0.0278211 0.8577175 0.9835665
## S g1 a33 t34 0.9469432 0.0313400 0.8401384 0.9837694
## S g1 a34 t35 0.9437253 0.0352632 0.8202682 0.9840311
## S g1 a35 t36 0.9403246 0.0396080 0.7980009 0.9843385
## S g1 a36 t37 0.9367321 0.0443937 0.7732654 0.9846806
## S g1 a37 t38 0.9329388 0.0496408 0.7460378 0.9850485
## S g1 a38 t39 0.9289353 0.0553713 0.7163544 0.9854348
# derived variabldes is the nest survival probability over the 40 (nocc) days
names(mod.res$results$derived)
## [1] "S Overall Survival"
mod.res$results$derived$"S Overall Survival"
## estimate se lcl ucl
## 1 0.3226129 0.1526737 0.1080227 0.651927
# alternatively
get.real(mod.res, "S", se=TRUE)
## all.diff.index par.index estimate se lcl
## S g1 a0 t1 1 1 0.9928771 0.0094787 0.9097296
## S g1 a1 t2 2 2 0.9924229 0.0097027 0.9126190
## S g1 a2 t3 3 3 0.9919398 0.0099182 0.9153905
## S g1 a3 t4 4 4 0.9914263 0.0101238 0.9180448
## S g1 a4 t5 5 5 0.9908803 0.0103177 0.9205823
## S g1 a5 t6 6 6 0.9902999 0.0104980 0.9230028
## S g1 a6 t7 7 7 0.9896829 0.0106632 0.9253056
## S g1 a7 t8 8 8 0.9890272 0.0108112 0.9274892
## S g1 a8 t9 9 9 0.9883302 0.0109404 0.9295514
## S g1 a9 t10 10 10 0.9875895 0.0110492 0.9314887
## S g1 a10 t11 11 11 0.9868025 0.0111361 0.9332964
## S g1 a11 t12 12 12 0.9859663 0.0112002 0.9349680
## S g1 a12 t13 13 13 0.9850778 0.0112410 0.9364948
## S g1 a13 t14 14 14 0.9841340 0.0112589 0.9378652
## S g1 a14 t15 15 15 0.9831316 0.0112554 0.9390639
## S g1 a15 t16 16 16 0.9820670 0.0112338 0.9400711
## S g1 a16 t17 17 17 0.9809365 0.0111994 0.9408604
## S g1 a17 t18 18 18 0.9797361 0.0111608 0.9413980
## S g1 a18 t19 19 19 0.9784619 0.0111303 0.9416397
## S g1 a19 t20 20 20 0.9771094 0.0111251 0.9415288
## S g1 a20 t21 21 21 0.9756741 0.0111683 0.9409938
## S g1 a21 t22 22 22 0.9741512 0.0112892 0.9399461
## S g1 a22 t23 23 23 0.9725356 0.0115233 0.9382794
## S g1 a23 t24 24 24 0.9708221 0.0119108 0.9358712
## S g1 a24 t25 25 25 0.9690050 0.0124937 0.9325866
## S g1 a25 t26 26 26 0.9670787 0.0133122 0.9282842
## S g1 a26 t27 27 27 0.9650369 0.0144012 0.9228208
## S g1 a27 t28 28 28 0.9628734 0.0157890 0.9160546
## S g1 a28 t29 29 29 0.9605815 0.0174969 0.9078453
## S g1 a29 t30 30 30 0.9581542 0.0195412 0.8980516
## S g1 a30 t31 31 31 0.9555844 0.0219351 0.8865300
## S g1 a31 t32 32 32 0.9528645 0.0246908 0.8731342
## S g1 a32 t33 33 33 0.9499868 0.0278211 0.8577175
## S g1 a33 t34 34 34 0.9469432 0.0313400 0.8401384
## S g1 a34 t35 35 35 0.9437253 0.0352632 0.8202682
## S g1 a35 t36 36 36 0.9403246 0.0396080 0.7980009
## S g1 a36 t37 37 37 0.9367321 0.0443937 0.7732654
## S g1 a37 t38 38 38 0.9329388 0.0496408 0.7460378
## S g1 a38 t39 39 39 0.9289353 0.0553713 0.7163544
## ucl fixed note group age time Age Time
## S g1 a0 t1 0.9994816 1 0 1 0 0
## S g1 a1 t2 0.9993915 1 1 2 1 1
## S g1 a2 t3 0.9992862 1 2 3 2 2
## S g1 a3 t4 0.9991630 1 3 4 3 3
## S g1 a4 t5 0.9990191 1 4 5 4 4
## S g1 a5 t6 0.9988512 1 5 6 5 5
## S g1 a6 t7 0.9986556 1 6 7 6 6
## S g1 a7 t8 0.9984280 1 7 8 7 7
## S g1 a8 t9 0.9981638 1 8 9 8 8
## S g1 a9 t10 0.9978576 1 9 10 9 9
## S g1 a10 t11 0.9975036 1 10 11 10 10
## S g1 a11 t12 0.9970958 1 11 12 11 11
## S g1 a12 t13 0.9966275 1 12 13 12 12
## S g1 a13 t14 0.9960922 1 13 14 13 13
## S g1 a14 t15 0.9954837 1 14 15 14 14
## S g1 a15 t16 0.9947966 1 15 16 15 15
## S g1 a16 t17 0.9940273 1 16 17 16 16
## S g1 a17 t18 0.9931748 1 17 18 17 17
## S g1 a18 t19 0.9922427 1 18 19 18 18
## S g1 a19 t20 0.9912401 1 19 20 19 19
## S g1 a20 t21 0.9901840 1 20 21 20 20
## S g1 a21 t22 0.9890999 1 21 22 21 21
## S g1 a22 t23 0.9880216 1 22 23 22 22
## S g1 a23 t24 0.9869892 1 23 24 23 23
## S g1 a24 t25 0.9860437 1 24 25 24 24
## S g1 a25 t26 0.9852215 1 25 26 25 25
## S g1 a26 t27 0.9845480 1 26 27 26 26
## S g1 a27 t28 0.9840351 1 27 28 27 27
## S g1 a28 t29 0.9836815 1 28 29 28 28
## S g1 a29 t30 0.9834759 1 29 30 29 29
## S g1 a30 t31 0.9834012 1 30 31 30 30
## S g1 a31 t32 0.9834378 1 31 32 31 31
## S g1 a32 t33 0.9835665 1 32 33 32 32
## S g1 a33 t34 0.9837694 1 33 34 33 33
## S g1 a34 t35 0.9840311 1 34 35 34 34
## S g1 a35 t36 0.9843385 1 35 36 35 35
## S g1 a36 t37 0.9846806 1 36 37 36 36
## S g1 a37 t38 0.9850485 1 37 38 37 37
## S g1 a38 t39 0.9854348 1 38 39 38 38
# Notice that the nest survival is the product of the individual S
prod(get.real(mod.res, "S", se=TRUE)$estimate)
## [1] 0.3226129
mod.res$design.matrix
## S:(Intercept) S:Time
## S g1 a0 t1 1 0
## S g1 a1 t2 1 1
## S g1 a2 t3 1 2
## S g1 a3 t4 1 3
## S g1 a4 t5 1 4
## S g1 a5 t6 1 5
## S g1 a6 t7 1 6
## S g1 a7 t8 1 7
## S g1 a8 t9 1 8
## S g1 a9 t10 1 9
## S g1 a10 t11 1 10
## S g1 a11 t12 1 11
## S g1 a12 t13 1 12
## S g1 a13 t14 1 13
## S g1 a14 t15 1 14
## S g1 a15 t16 1 15
## S g1 a16 t17 1 16
## S g1 a17 t18 1 17
## S g1 a18 t19 1 18
## S g1 a19 t20 1 19
## S g1 a20 t21 1 20
## S g1 a21 t22 1 21
## S g1 a22 t23 1 22
## S g1 a23 t24 1 23
## S g1 a24 t25 1 24
## S g1 a25 t26 1 25
## S g1 a26 t27 1 26
## S g1 a27 t28 1 27
## S g1 a28 t29 1 28
## S g1 a29 t30 1 29
## S g1 a30 t31 1 30
## S g1 a31 t32 1 31
## S g1 a32 t33 1 32
## S g1 a33 t34 1 33
## S g1 a34 t35 1 34
## S g1 a35 t36 1 35
## S g1 a36 t37 1 36
## S g1 a37 t38 1 37
## S g1 a38 t39 1 38
# plot the trend over time
plotdata <- get.real(mod.res, "S", se=TRUE)
plot.res <- ggplot(data=plotdata, aes(x=time, y=estimate, group=1))+
ggtitle("Model with linear trend in S on logit scale")+
geom_line()+
geom_ribbon(aes(x=time, ymin=lcl, ymax=ucl), alpha=0.2)+
#geom_errorbar(aes(ymin=lcl, ymax=ucl), width=.1)+
ylim(0,1)
plot.res
ggsave(plot.res,
file=file.path("..","..","..","..","MyStuff","Images","killdeer-linear-S.png"), h=4, w=6, units="in", dpi=300)