# Analysis of killdeer data illustrating basic RMark features
# Fit a model with DSR differing in first and second half of the study
# 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.
# Here we create a variable for first/second half of the study
kill.ddl <- make.design.data(kill.proc)
kill.ddl$S$studyhalf <- car::recode(kill.ddl$S$Time,
" lo:20='1st'; 21:hi='2nd'")
str(kill.ddl)
## List of 2
## $ S :'data.frame': 39 obs. of 8 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 ...
## ..$ studyhalf : chr [1:39] "1st" "1st" "1st" "1st" ...
## $ pimtypes:List of 1
## ..$ S:List of 1
## .. ..$ pim.type: chr "all"
kill.ddl
## $S
## par.index model.index group age time Age Time studyhalf
## 1 1 1 1 0 1 0 0 1st
## 2 2 2 1 1 2 1 1 1st
## 3 3 3 1 2 3 2 2 1st
## 4 4 4 1 3 4 3 3 1st
## 5 5 5 1 4 5 4 4 1st
## 6 6 6 1 5 6 5 5 1st
## 7 7 7 1 6 7 6 6 1st
## 8 8 8 1 7 8 7 7 1st
## 9 9 9 1 8 9 8 8 1st
## 10 10 10 1 9 10 9 9 1st
## 11 11 11 1 10 11 10 10 1st
## 12 12 12 1 11 12 11 11 1st
## 13 13 13 1 12 13 12 12 1st
## 14 14 14 1 13 14 13 13 1st
## 15 15 15 1 14 15 14 14 1st
## 16 16 16 1 15 16 15 15 1st
## 17 17 17 1 16 17 16 16 1st
## 18 18 18 1 17 18 17 17 1st
## 19 19 19 1 18 19 18 18 1st
## 20 20 20 1 19 20 19 19 1st
## 21 21 21 1 20 21 20 20 1st
## 22 22 22 1 21 22 21 21 2nd
## 23 23 23 1 22 23 22 22 2nd
## 24 24 24 1 23 24 23 23 2nd
## 25 25 25 1 24 25 24 24 2nd
## 26 26 26 1 25 26 25 25 2nd
## 27 27 27 1 26 27 26 26 2nd
## 28 28 28 1 27 28 27 27 2nd
## 29 29 29 1 28 29 28 28 2nd
## 30 30 30 1 29 30 29 29 2nd
## 31 31 31 1 30 31 30 30 2nd
## 32 32 32 1 31 32 31 31 2nd
## 33 33 33 1 32 33 32 32 2nd
## 34 34 34 1 33 34 33 33 2nd
## 35 35 35 1 34 35 34 34 2nd
## 36 36 36 1 35 36 35 35 2nd
## 37 37 37 1 36 37 36 36 2nd
## 38 38 38 1 37 38 37 37 2nd
## 39 39 39 1 38 39 38 38 2nd
##
## $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=~studyhalf)
)
)
## Warning: Unknown or uninitialised column: 'ch'.
## Warning: Unknown or uninitialised column: 'AgeDay1'.
##
## Output summary for Nest model
## Name : S(~studyhalf)
##
## Npar : 2
## -2lnL: 41.92815
## AICc : 45.98612
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 3.9413016 0.7140400 2.541783 5.340820
## S:studyhalf2nd -0.6514382 0.8772159 -2.370781 1.067905
##
##
## Real Parameter S
## 1 2 3 4 5 6 7
## 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471
## 8 9 10 11 12 13 14
## 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471
## 15 16 17 18 19 20 21
## 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471
## 22 23 24 25 26 27 28
## 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794
## 29 30 31 32 33 34 35
## 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794
## 36 37 38 39
## 0.9640794 0.9640794 0.9640794 0.9640794
summary(mod.res)
## Output summary for Nest model
## Name : S(~studyhalf)
##
## Npar : 2
## -2lnL: 41.92815
## AICc : 45.98612
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 3.9413016 0.7140400 2.541783 5.340820
## S:studyhalf2nd -0.6514382 0.8772159 -2.370781 1.067905
##
##
## Real Parameter S
## 1 2 3 4 5 6 7
## 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471
## 8 9 10 11 12 13 14
## 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471
## 15 16 17 18 19 20 21
## 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471 0.9809471
## 22 23 24 25 26 27 28
## 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794
## 29 30 31 32 33 34 35
## 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794 0.9640794
## 36 37 38 39
## 0.9640794 0.9640794 0.9640794 0.9640794
# 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) 3.9413016 0.7140400 2.541783 5.340820
## S:studyhalf2nd -0.6514382 0.8772159 -2.370781 1.067905
mod.res$results$real# on the regular 0-1 scale for each site
## estimate se lcl ucl fixed note
## S g1 a0 t1 0.9809471 0.0133453 0.9270196 0.9952309
## S g1 a21 t22 0.9640794 0.0176463 0.9081389 0.9864618
# 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.3456117 0.1507159 0.1251511 0.6610013
# alternatively
get.real(mod.res, "S", se=TRUE)
## all.diff.index par.index estimate se lcl
## S g1 a0 t1 1 1 0.9809471 0.0133453 0.9270196
## S g1 a1 t2 2 1 0.9809471 0.0133453 0.9270196
## S g1 a2 t3 3 1 0.9809471 0.0133453 0.9270196
## S g1 a3 t4 4 1 0.9809471 0.0133453 0.9270196
## S g1 a4 t5 5 1 0.9809471 0.0133453 0.9270196
## S g1 a5 t6 6 1 0.9809471 0.0133453 0.9270196
## S g1 a6 t7 7 1 0.9809471 0.0133453 0.9270196
## S g1 a7 t8 8 1 0.9809471 0.0133453 0.9270196
## S g1 a8 t9 9 1 0.9809471 0.0133453 0.9270196
## S g1 a9 t10 10 1 0.9809471 0.0133453 0.9270196
## S g1 a10 t11 11 1 0.9809471 0.0133453 0.9270196
## S g1 a11 t12 12 1 0.9809471 0.0133453 0.9270196
## S g1 a12 t13 13 1 0.9809471 0.0133453 0.9270196
## S g1 a13 t14 14 1 0.9809471 0.0133453 0.9270196
## S g1 a14 t15 15 1 0.9809471 0.0133453 0.9270196
## S g1 a15 t16 16 1 0.9809471 0.0133453 0.9270196
## S g1 a16 t17 17 1 0.9809471 0.0133453 0.9270196
## S g1 a17 t18 18 1 0.9809471 0.0133453 0.9270196
## S g1 a18 t19 19 1 0.9809471 0.0133453 0.9270196
## S g1 a19 t20 20 1 0.9809471 0.0133453 0.9270196
## S g1 a20 t21 21 1 0.9809471 0.0133453 0.9270196
## S g1 a21 t22 22 2 0.9640794 0.0176463 0.9081389
## S g1 a22 t23 23 2 0.9640794 0.0176463 0.9081389
## S g1 a23 t24 24 2 0.9640794 0.0176463 0.9081389
## S g1 a24 t25 25 2 0.9640794 0.0176463 0.9081389
## S g1 a25 t26 26 2 0.9640794 0.0176463 0.9081389
## S g1 a26 t27 27 2 0.9640794 0.0176463 0.9081389
## S g1 a27 t28 28 2 0.9640794 0.0176463 0.9081389
## S g1 a28 t29 29 2 0.9640794 0.0176463 0.9081389
## S g1 a29 t30 30 2 0.9640794 0.0176463 0.9081389
## S g1 a30 t31 31 2 0.9640794 0.0176463 0.9081389
## S g1 a31 t32 32 2 0.9640794 0.0176463 0.9081389
## S g1 a32 t33 33 2 0.9640794 0.0176463 0.9081389
## S g1 a33 t34 34 2 0.9640794 0.0176463 0.9081389
## S g1 a34 t35 35 2 0.9640794 0.0176463 0.9081389
## S g1 a35 t36 36 2 0.9640794 0.0176463 0.9081389
## S g1 a36 t37 37 2 0.9640794 0.0176463 0.9081389
## S g1 a37 t38 38 2 0.9640794 0.0176463 0.9081389
## S g1 a38 t39 39 2 0.9640794 0.0176463 0.9081389
## ucl fixed note group age time Age Time
## S g1 a0 t1 0.9952309 1 0 1 0 0
## S g1 a1 t2 0.9952309 1 1 2 1 1
## S g1 a2 t3 0.9952309 1 2 3 2 2
## S g1 a3 t4 0.9952309 1 3 4 3 3
## S g1 a4 t5 0.9952309 1 4 5 4 4
## S g1 a5 t6 0.9952309 1 5 6 5 5
## S g1 a6 t7 0.9952309 1 6 7 6 6
## S g1 a7 t8 0.9952309 1 7 8 7 7
## S g1 a8 t9 0.9952309 1 8 9 8 8
## S g1 a9 t10 0.9952309 1 9 10 9 9
## S g1 a10 t11 0.9952309 1 10 11 10 10
## S g1 a11 t12 0.9952309 1 11 12 11 11
## S g1 a12 t13 0.9952309 1 12 13 12 12
## S g1 a13 t14 0.9952309 1 13 14 13 13
## S g1 a14 t15 0.9952309 1 14 15 14 14
## S g1 a15 t16 0.9952309 1 15 16 15 15
## S g1 a16 t17 0.9952309 1 16 17 16 16
## S g1 a17 t18 0.9952309 1 17 18 17 17
## S g1 a18 t19 0.9952309 1 18 19 18 18
## S g1 a19 t20 0.9952309 1 19 20 19 19
## S g1 a20 t21 0.9952309 1 20 21 20 20
## S g1 a21 t22 0.9864618 1 21 22 21 21
## S g1 a22 t23 0.9864618 1 22 23 22 22
## S g1 a23 t24 0.9864618 1 23 24 23 23
## S g1 a24 t25 0.9864618 1 24 25 24 24
## S g1 a25 t26 0.9864618 1 25 26 25 25
## S g1 a26 t27 0.9864618 1 26 27 26 26
## S g1 a27 t28 0.9864618 1 27 28 27 27
## S g1 a28 t29 0.9864618 1 28 29 28 28
## S g1 a29 t30 0.9864618 1 29 30 29 29
## S g1 a30 t31 0.9864618 1 30 31 30 30
## S g1 a31 t32 0.9864618 1 31 32 31 31
## S g1 a32 t33 0.9864618 1 32 33 32 32
## S g1 a33 t34 0.9864618 1 33 34 33 33
## S g1 a34 t35 0.9864618 1 34 35 34 34
## S g1 a35 t36 0.9864618 1 35 36 35 35
## S g1 a36 t37 0.9864618 1 36 37 36 36
## S g1 a37 t38 0.9864618 1 37 38 37 37
## S g1 a38 t39 0.9864618 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.3456112
mod.res$design.matrix
## S:(Intercept) S:studyhalf2nd
## S g1 a0 t1 1 0
## S g1 a21 t22 1 1
# 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 two groups of S")+
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-2groups-S.png"), h=4, w=6, units="in", dpi=300)