# Analysis of killdeer data illustrating model averaging
# 2019-05-01 CJS Initial version
# 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(RMark)
## This is RMark 2.2.6
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
source(file.path("..","..","RMark.additional.functions.r"))
# 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 you could standardize covariates
kill.ddl <- make.design.data(kill.proc)
kill.ddl$S$Time2 <- (kill.ddl$S$Time-20)^2
kill.ddl
## $S
## par.index model.index group age time Age Time Time2
## 1 1 1 1 0 1 0 0 400
## 2 2 2 1 1 2 1 1 361
## 3 3 3 1 2 3 2 2 324
## 4 4 4 1 3 4 3 3 289
## 5 5 5 1 4 5 4 4 256
## 6 6 6 1 5 6 5 5 225
## 7 7 7 1 6 7 6 6 196
## 8 8 8 1 7 8 7 7 169
## 9 9 9 1 8 9 8 8 144
## 10 10 10 1 9 10 9 9 121
## 11 11 11 1 10 11 10 10 100
## 12 12 12 1 11 12 11 11 81
## 13 13 13 1 12 13 12 12 64
## 14 14 14 1 13 14 13 13 49
## 15 15 15 1 14 15 14 14 36
## 16 16 16 1 15 16 15 15 25
## 17 17 17 1 16 17 16 16 16
## 18 18 18 1 17 18 17 17 9
## 19 19 19 1 18 19 18 18 4
## 20 20 20 1 19 20 19 19 1
## 21 21 21 1 20 21 20 20 0
## 22 22 22 1 21 22 21 21 1
## 23 23 23 1 22 23 22 22 4
## 24 24 24 1 23 24 23 23 9
## 25 25 25 1 24 25 24 24 16
## 26 26 26 1 25 26 25 25 25
## 27 27 27 1 26 27 26 26 36
## 28 28 28 1 27 28 27 27 49
## 29 29 29 1 28 29 28 28 64
## 30 30 30 1 29 30 29 29 81
## 31 31 31 1 30 31 30 30 100
## 32 32 32 1 31 32 31 31 121
## 33 33 33 1 32 33 32 32 144
## 34 34 34 1 33 34 33 33 169
## 35 35 35 1 34 35 34 34 196
## 36 36 36 1 35 36 35 35 225
## 37 37 37 1 36 37 36 36 256
## 38 38 38 1 37 38 37 37 289
## 39 39 39 1 38 39 38 38 324
##
## $pimtypes
## $pimtypes$S
## $pimtypes$S$pim.type
## [1] "all"
# 3. Set up the set of model to fit
model.list.csv <- textConnection(
" S
~1
~Time
~Time+Time2
")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list$model.number <- 1:nrow(model.list)
model.list
## S model.number
## 1 ~1 1
## 2 ~Time 2
## 3 ~Time+Time2 3
# fit the models
myobs <- ls()
myobs <- myobs[ grepl("m...",myobs,fixed=TRUE)]
cat("Removing ", myobs, "\n")
## Removing
rm(list=myobs)
model.fits <- plyr::dlply(model.list, "model.number", function(x,input.data, input.ddl){
cat("\n\n***** Starting ", unlist(x), "\n")
fit <- RMark::mark(input.data, ddl=input.ddl,
model="Nest",
model.parameters=list(
S =list(formula=as.formula(eval(x$S)))
)
#,brief=TRUE,output=FALSE, delete=TRUE
#,invisible=TRUE,output=TRUE # set for debugging
)
mnumber <- paste("m...",formatC(x$model.number, width = 3, format = "d", flag = "0"),sep="")
assign( mnumber, fit, envir=.GlobalEnv)
#browser()
fit
},input.data=kill.proc, input.ddl=kill.ddl)
##
##
## ***** Starting ~1 1
## Warning: Unknown or uninitialised column: 'ch'.
## Warning: Unknown or uninitialised column: 'AgeDay1'.
##
## Output summary for Nest model
## Name : S(~1)
##
## Npar : 1
## -2lnL: 42.51028
## AICc : 44.52951
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 3.557002 0.4141776 2.745214 4.368791
##
##
## Real Parameter S
## 1 2 3 4 5 6 7
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 8 9 10 11 12 13 14
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 15 16 17 18 19 20 21
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 22 23 24 25 26 27 28
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 29 30 31 32 33 34 35
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 36 37 38 39
## 0.9722669 0.9722669 0.9722669 0.9722669
##
##
## ***** Starting ~Time 2
## 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
##
##
## ***** Starting ~Time+Time2 3
## Warning: Unknown or uninitialised column: 'ch'.
## Warning: Unknown or uninitialised column: 'AgeDay1'.
##
## Output summary for Nest model
## Name : S(~Time + Time2)
##
## Npar : 3
## -2lnL: 38.08174
## AICc : 44.19824
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 5.3641150 1.1259280 3.1572960 7.5709339000
## S:Time -0.0474837 0.0382368 -0.1224279 0.0274605000
## S:Time2 -0.0073305 0.0040110 -0.0151920 0.0005309712
##
##
## Real Parameter S
## 1 2 3 4 5 6 7
## 0.9192284 0.9352522 0.9475521 0.9570223 0.9643453 0.9700358 0.9744798
## 8 9 10 11 12 13 14
## 0.9779657 0.9807093 0.9828721 0.9845751 0.9859089 0.9869413 0.9877225
## 15 16 17 18 19 20 21
## 0.988289 0.9886665 0.9888719 0.9889139 0.9887946 0.9885086 0.9880434
## 22 23 24 25 26 27 28
## 0.9873782 0.9864824 0.9853133 0.9838129 0.9819032 0.9794804 0.9764054
## 29 30 31 32 33 34 35
## 0.9724929 0.967495 0.9610801 0.9528053 0.9420809 0.9281289 0.9099389
## 36 37 38 39
## 0.8862323 0.855461 0.8158832 0.7657827
# examine individula model results
model.number <-1
names(model.fits[[model.number]])
## [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"
model.fits[[model.number]]$model.name
## [1] "S(~1)"
summary(model.fits[[model.number]])
## Output summary for Nest model
## Name : S(~1)
##
## Npar : 1
## -2lnL: 42.51028
## AICc : 44.52951
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 3.557002 0.4141776 2.745214 4.368791
##
##
## Real Parameter S
## 1 2 3 4 5 6 7
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 8 9 10 11 12 13 14
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 15 16 17 18 19 20 21
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 22 23 24 25 26 27 28
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 29 30 31 32 33 34 35
## 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669 0.9722669
## 36 37 38 39
## 0.9722669 0.9722669 0.9722669 0.9722669
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## S g1 a0 t1 0.9722669 0.0111679 0.9396425 0.9874919
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## S:(Intercept) 3.557002 0.4141776 2.745214 4.368791
model.fits[[model.number]]$results$derived
## $`S Overall Survival`
## estimate se lcl ucl
## 1 0.3339134 0.1495833 0.1182906 0.6519547
get.real(model.fits[[model.number]], "S", se=TRUE)
## all.diff.index par.index estimate se lcl
## S g1 a0 t1 1 1 0.9722669 0.0111679 0.9396425
## S g1 a1 t2 2 1 0.9722669 0.0111679 0.9396425
## S g1 a2 t3 3 1 0.9722669 0.0111679 0.9396425
## S g1 a3 t4 4 1 0.9722669 0.0111679 0.9396425
## S g1 a4 t5 5 1 0.9722669 0.0111679 0.9396425
## S g1 a5 t6 6 1 0.9722669 0.0111679 0.9396425
## S g1 a6 t7 7 1 0.9722669 0.0111679 0.9396425
## S g1 a7 t8 8 1 0.9722669 0.0111679 0.9396425
## S g1 a8 t9 9 1 0.9722669 0.0111679 0.9396425
## S g1 a9 t10 10 1 0.9722669 0.0111679 0.9396425
## S g1 a10 t11 11 1 0.9722669 0.0111679 0.9396425
## S g1 a11 t12 12 1 0.9722669 0.0111679 0.9396425
## S g1 a12 t13 13 1 0.9722669 0.0111679 0.9396425
## S g1 a13 t14 14 1 0.9722669 0.0111679 0.9396425
## S g1 a14 t15 15 1 0.9722669 0.0111679 0.9396425
## S g1 a15 t16 16 1 0.9722669 0.0111679 0.9396425
## S g1 a16 t17 17 1 0.9722669 0.0111679 0.9396425
## S g1 a17 t18 18 1 0.9722669 0.0111679 0.9396425
## S g1 a18 t19 19 1 0.9722669 0.0111679 0.9396425
## S g1 a19 t20 20 1 0.9722669 0.0111679 0.9396425
## S g1 a20 t21 21 1 0.9722669 0.0111679 0.9396425
## S g1 a21 t22 22 1 0.9722669 0.0111679 0.9396425
## S g1 a22 t23 23 1 0.9722669 0.0111679 0.9396425
## S g1 a23 t24 24 1 0.9722669 0.0111679 0.9396425
## S g1 a24 t25 25 1 0.9722669 0.0111679 0.9396425
## S g1 a25 t26 26 1 0.9722669 0.0111679 0.9396425
## S g1 a26 t27 27 1 0.9722669 0.0111679 0.9396425
## S g1 a27 t28 28 1 0.9722669 0.0111679 0.9396425
## S g1 a28 t29 29 1 0.9722669 0.0111679 0.9396425
## S g1 a29 t30 30 1 0.9722669 0.0111679 0.9396425
## S g1 a30 t31 31 1 0.9722669 0.0111679 0.9396425
## S g1 a31 t32 32 1 0.9722669 0.0111679 0.9396425
## S g1 a32 t33 33 1 0.9722669 0.0111679 0.9396425
## S g1 a33 t34 34 1 0.9722669 0.0111679 0.9396425
## S g1 a34 t35 35 1 0.9722669 0.0111679 0.9396425
## S g1 a35 t36 36 1 0.9722669 0.0111679 0.9396425
## S g1 a36 t37 37 1 0.9722669 0.0111679 0.9396425
## S g1 a37 t38 38 1 0.9722669 0.0111679 0.9396425
## S g1 a38 t39 39 1 0.9722669 0.0111679 0.9396425
## ucl fixed note group age time Age Time
## S g1 a0 t1 0.9874919 1 0 1 0 0
## S g1 a1 t2 0.9874919 1 1 2 1 1
## S g1 a2 t3 0.9874919 1 2 3 2 2
## S g1 a3 t4 0.9874919 1 3 4 3 3
## S g1 a4 t5 0.9874919 1 4 5 4 4
## S g1 a5 t6 0.9874919 1 5 6 5 5
## S g1 a6 t7 0.9874919 1 6 7 6 6
## S g1 a7 t8 0.9874919 1 7 8 7 7
## S g1 a8 t9 0.9874919 1 8 9 8 8
## S g1 a9 t10 0.9874919 1 9 10 9 9
## S g1 a10 t11 0.9874919 1 10 11 10 10
## S g1 a11 t12 0.9874919 1 11 12 11 11
## S g1 a12 t13 0.9874919 1 12 13 12 12
## S g1 a13 t14 0.9874919 1 13 14 13 13
## S g1 a14 t15 0.9874919 1 14 15 14 14
## S g1 a15 t16 0.9874919 1 15 16 15 15
## S g1 a16 t17 0.9874919 1 16 17 16 16
## S g1 a17 t18 0.9874919 1 17 18 17 17
## S g1 a18 t19 0.9874919 1 18 19 18 18
## S g1 a19 t20 0.9874919 1 19 20 19 19
## S g1 a20 t21 0.9874919 1 20 21 20 20
## S g1 a21 t22 0.9874919 1 21 22 21 21
## S g1 a22 t23 0.9874919 1 22 23 22 22
## S g1 a23 t24 0.9874919 1 23 24 23 23
## S g1 a24 t25 0.9874919 1 24 25 24 24
## S g1 a25 t26 0.9874919 1 25 26 25 25
## S g1 a26 t27 0.9874919 1 26 27 26 26
## S g1 a27 t28 0.9874919 1 27 28 27 27
## S g1 a28 t29 0.9874919 1 28 29 28 28
## S g1 a29 t30 0.9874919 1 29 30 29 29
## S g1 a30 t31 0.9874919 1 30 31 30 30
## S g1 a31 t32 0.9874919 1 31 32 31 31
## S g1 a32 t33 0.9874919 1 32 33 32 32
## S g1 a33 t34 0.9874919 1 33 34 33 33
## S g1 a34 t35 0.9874919 1 34 35 34 34
## S g1 a35 t36 0.9874919 1 35 36 35 35
## S g1 a36 t37 0.9874919 1 36 37 36 36
## S g1 a37 t38 0.9874919 1 37 38 37 37
## S g1 a38 t39 0.9874919 1 38 39 38 38
# Model comparision and averaging
# collect models and make AICc table
model.set <- RMark::collect.models( type="Nest")
model.set
## model npar AICc DeltaAICc weight Deviance
## 3 S(~Time + Time2) 3 44.19824 0.0000000 0.3992023 38.08174
## 1 S(~1) 1 44.52951 0.3312659 0.3382669 42.51028
## 2 S(~Time) 2 45.03644 0.8381992 0.2625309 40.97847
names(model.set)
## [1] "m...001" "m...002" "m...003" "model.table"
model.set$model.table
## S model npar AICc DeltaAICc weight
## 3 ~Time + Time2 S(~Time + Time2) 3 44.19824 0.0000000 0.3992023
## 1 ~1 S(~1) 1 44.52951 0.3312659 0.3382669
## 2 ~Time S(~Time) 2 45.03644 0.8381992 0.2625309
## Deviance
## 3 38.08174
## 1 42.51028
## 2 40.97847
# model averaged values
get.real(model.set[[1]], "S", se=TRUE)
## all.diff.index par.index estimate se lcl
## S g1 a0 t1 1 1 0.9722669 0.0111679 0.9396425
## S g1 a1 t2 2 1 0.9722669 0.0111679 0.9396425
## S g1 a2 t3 3 1 0.9722669 0.0111679 0.9396425
## S g1 a3 t4 4 1 0.9722669 0.0111679 0.9396425
## S g1 a4 t5 5 1 0.9722669 0.0111679 0.9396425
## S g1 a5 t6 6 1 0.9722669 0.0111679 0.9396425
## S g1 a6 t7 7 1 0.9722669 0.0111679 0.9396425
## S g1 a7 t8 8 1 0.9722669 0.0111679 0.9396425
## S g1 a8 t9 9 1 0.9722669 0.0111679 0.9396425
## S g1 a9 t10 10 1 0.9722669 0.0111679 0.9396425
## S g1 a10 t11 11 1 0.9722669 0.0111679 0.9396425
## S g1 a11 t12 12 1 0.9722669 0.0111679 0.9396425
## S g1 a12 t13 13 1 0.9722669 0.0111679 0.9396425
## S g1 a13 t14 14 1 0.9722669 0.0111679 0.9396425
## S g1 a14 t15 15 1 0.9722669 0.0111679 0.9396425
## S g1 a15 t16 16 1 0.9722669 0.0111679 0.9396425
## S g1 a16 t17 17 1 0.9722669 0.0111679 0.9396425
## S g1 a17 t18 18 1 0.9722669 0.0111679 0.9396425
## S g1 a18 t19 19 1 0.9722669 0.0111679 0.9396425
## S g1 a19 t20 20 1 0.9722669 0.0111679 0.9396425
## S g1 a20 t21 21 1 0.9722669 0.0111679 0.9396425
## S g1 a21 t22 22 1 0.9722669 0.0111679 0.9396425
## S g1 a22 t23 23 1 0.9722669 0.0111679 0.9396425
## S g1 a23 t24 24 1 0.9722669 0.0111679 0.9396425
## S g1 a24 t25 25 1 0.9722669 0.0111679 0.9396425
## S g1 a25 t26 26 1 0.9722669 0.0111679 0.9396425
## S g1 a26 t27 27 1 0.9722669 0.0111679 0.9396425
## S g1 a27 t28 28 1 0.9722669 0.0111679 0.9396425
## S g1 a28 t29 29 1 0.9722669 0.0111679 0.9396425
## S g1 a29 t30 30 1 0.9722669 0.0111679 0.9396425
## S g1 a30 t31 31 1 0.9722669 0.0111679 0.9396425
## S g1 a31 t32 32 1 0.9722669 0.0111679 0.9396425
## S g1 a32 t33 33 1 0.9722669 0.0111679 0.9396425
## S g1 a33 t34 34 1 0.9722669 0.0111679 0.9396425
## S g1 a34 t35 35 1 0.9722669 0.0111679 0.9396425
## S g1 a35 t36 36 1 0.9722669 0.0111679 0.9396425
## S g1 a36 t37 37 1 0.9722669 0.0111679 0.9396425
## S g1 a37 t38 38 1 0.9722669 0.0111679 0.9396425
## S g1 a38 t39 39 1 0.9722669 0.0111679 0.9396425
## ucl fixed note group age time Age Time
## S g1 a0 t1 0.9874919 1 0 1 0 0
## S g1 a1 t2 0.9874919 1 1 2 1 1
## S g1 a2 t3 0.9874919 1 2 3 2 2
## S g1 a3 t4 0.9874919 1 3 4 3 3
## S g1 a4 t5 0.9874919 1 4 5 4 4
## S g1 a5 t6 0.9874919 1 5 6 5 5
## S g1 a6 t7 0.9874919 1 6 7 6 6
## S g1 a7 t8 0.9874919 1 7 8 7 7
## S g1 a8 t9 0.9874919 1 8 9 8 8
## S g1 a9 t10 0.9874919 1 9 10 9 9
## S g1 a10 t11 0.9874919 1 10 11 10 10
## S g1 a11 t12 0.9874919 1 11 12 11 11
## S g1 a12 t13 0.9874919 1 12 13 12 12
## S g1 a13 t14 0.9874919 1 13 14 13 13
## S g1 a14 t15 0.9874919 1 14 15 14 14
## S g1 a15 t16 0.9874919 1 15 16 15 15
## S g1 a16 t17 0.9874919 1 16 17 16 16
## S g1 a17 t18 0.9874919 1 17 18 17 17
## S g1 a18 t19 0.9874919 1 18 19 18 18
## S g1 a19 t20 0.9874919 1 19 20 19 19
## S g1 a20 t21 0.9874919 1 20 21 20 20
## S g1 a21 t22 0.9874919 1 21 22 21 21
## S g1 a22 t23 0.9874919 1 22 23 22 22
## S g1 a23 t24 0.9874919 1 23 24 23 23
## S g1 a24 t25 0.9874919 1 24 25 24 24
## S g1 a25 t26 0.9874919 1 25 26 25 25
## S g1 a26 t27 0.9874919 1 26 27 26 26
## S g1 a27 t28 0.9874919 1 27 28 27 27
## S g1 a28 t29 0.9874919 1 28 29 28 28
## S g1 a29 t30 0.9874919 1 29 30 29 29
## S g1 a30 t31 0.9874919 1 30 31 30 30
## S g1 a31 t32 0.9874919 1 31 32 31 31
## S g1 a32 t33 0.9874919 1 32 33 32 32
## S g1 a33 t34 0.9874919 1 33 34 33 33
## S g1 a34 t35 0.9874919 1 34 35 34 34
## S g1 a35 t36 0.9874919 1 35 36 35 35
## S g1 a36 t37 0.9874919 1 36 37 36 36
## S g1 a37 t38 0.9874919 1 37 38 37 37
## S g1 a38 t39 0.9874919 1 38 39 38 38
get.real(model.set[[2]], "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
get.real(model.set[[3]], "S", se=TRUE)
## all.diff.index par.index estimate se lcl
## S g1 a0 t1 1 1 0.9192284 0.1102151 0.3828261
## S g1 a1 t2 2 2 0.9352522 0.0815663 0.5075650
## S g1 a2 t3 3 3 0.9475521 0.0606762 0.6227081
## S g1 a3 t4 4 4 0.9570223 0.0455284 0.7178004
## S g1 a4 t5 5 5 0.9643453 0.0345784 0.7902535
## S g1 a5 t6 6 6 0.9700358 0.0266789 0.8426819
## S g1 a6 t7 7 7 0.9744798 0.0209919 0.8795312
## S g1 a7 t8 8 8 0.9779657 0.0169097 0.9050662
## S g1 a8 t9 9 9 0.9807093 0.0139921 0.9226596
## S g1 a9 t10 10 10 0.9828721 0.0119187 0.9347530
## S g1 a10 t11 11 11 0.9845751 0.0104560 0.9430419
## S g1 a11 t12 12 12 0.9859089 0.0094339 0.9486830
## S g1 a12 t13 13 13 0.9869413 0.0087303 0.9524598
## S g1 a13 t14 14 14 0.9877225 0.0082589 0.9549026
## S g1 a14 t15 15 15 0.9882890 0.0079603 0.9563700
## S g1 a15 t16 16 16 0.9886665 0.0077946 0.9571045
## S g1 a16 t17 17 17 0.9888719 0.0077362 0.9572696
## S g1 a17 t18 18 18 0.9889139 0.0077690 0.9569744
## S g1 a18 t19 19 19 0.9887946 0.0078838 0.9562890
## S g1 a19 t20 20 20 0.9885086 0.0080767 0.9552551
## S g1 a20 t21 21 21 0.9880434 0.0083473 0.9538910
## S g1 a21 t22 22 22 0.9873782 0.0086984 0.9521941
## S g1 a22 t23 23 23 0.9864824 0.0091355 0.9501402
## S g1 a23 t24 24 24 0.9853133 0.0096674 0.9476797
## S g1 a24 t25 25 25 0.9838129 0.0103073 0.9447285
## S g1 a25 t26 26 26 0.9819032 0.0110748 0.9411536
## S g1 a26 t27 27 27 0.9794804 0.0120013 0.9367474
## S g1 a27 t28 28 28 0.9764054 0.0131391 0.9311841
## S g1 a28 t29 29 29 0.9724929 0.0145787 0.9239480
## S g1 a29 t30 30 30 0.9674950 0.0164795 0.9142168
## S g1 a30 t31 31 31 0.9610801 0.0191170 0.9006827
## S g1 a31 t32 32 32 0.9528053 0.0229478 0.8813107
## S g1 a32 t33 33 33 0.9420809 0.0286777 0.8530743
## S g1 a33 t34 34 34 0.9281289 0.0373170 0.8118107
## S g1 a34 t35 35 35 0.9099389 0.0502081 0.7525133
## S g1 a35 t36 36 36 0.8862323 0.0690157 0.6706610
## S g1 a36 t37 37 37 0.8554610 0.0956219 0.5652110
## S g1 a37 t38 38 38 0.8158832 0.1317955 0.4425227
## S g1 a38 t39 39 39 0.7657827 0.1784566 0.3174569
## ucl fixed note group age time Age Time
## S g1 a0 t1 0.9952336 1 0 1 0 0
## S g1 a1 t2 0.9950842 1 1 2 1 1
## S g1 a2 t3 0.9949689 1 2 3 2 2
## S g1 a3 t4 0.9948965 1 3 4 3 3
## S g1 a4 t5 0.9948760 1 4 5 4 4
## S g1 a5 t6 0.9949149 1 5 6 5 5
## S g1 a6 t7 0.9950177 1 6 7 6 6
## S g1 a7 t8 0.9951837 1 7 8 7 7
## S g1 a8 t9 0.9954054 1 8 9 8 8
## S g1 a9 t10 0.9956682 1 9 10 9 9
## S g1 a10 t11 0.9959527 1 10 11 10 10
## S g1 a11 t12 0.9962378 1 11 12 11 11
## S g1 a12 t13 0.9965047 1 12 13 12 12
## S g1 a13 t14 0.9967391 1 13 14 13 13
## S g1 a14 t15 0.9969315 1 14 15 14 14
## S g1 a15 t16 0.9970765 1 15 16 15 15
## S g1 a16 t17 0.9971710 1 16 17 16 16
## S g1 a17 t18 0.9972126 1 17 18 17 17
## S g1 a18 t19 0.9971983 1 18 19 18 18
## S g1 a19 t20 0.9971232 1 19 20 19 19
## S g1 a20 t21 0.9969796 1 20 21 20 20
## S g1 a21 t22 0.9967558 1 21 22 21 21
## S g1 a22 t23 0.9964346 1 22 23 22 22
## S g1 a23 t24 0.9959918 1 23 24 23 23
## S g1 a24 t25 0.9953941 1 24 25 24 24
## S g1 a25 t26 0.9945968 1 25 26 25 25
## S g1 a26 t27 0.9935423 1 26 27 26 26
## S g1 a27 t28 0.9921604 1 27 28 27 27
## S g1 a28 t29 0.9903739 1 28 29 28 28
## S g1 a29 t30 0.9881134 1 29 30 29 29
## S g1 a30 t31 0.9853459 1 30 31 30 30
## S g1 a31 t32 0.9821082 1 31 32 31 31
## S g1 a32 t33 0.9785252 1 32 33 32 32
## S g1 a33 t34 0.9747849 1 33 34 33 33
## S g1 a34 t35 0.9710755 1 34 35 34 34
## S g1 a35 t36 0.9675311 1 35 36 35 35
## S g1 a36 t37 0.9642170 1 36 37 36 36
## S g1 a37 t38 0.9611466 1 37 38 37 37
## S g1 a38 t39 0.9583049 1 38 39 38 38
S.ma <- RMark::model.average(model.set, param="S")
S.ma
## par.index estimate se fixed note group age time
## S g1 a0 t1 1 0.9565048 0.07681821 1 0 1
## S g1 a1 t2 2 0.9627822 0.05732738 1 1 2
## S g1 a2 t3 3 0.9675655 0.04313985 1 2 3
## S g1 a3 t4 4 0.9712112 0.03293324 1 3 4
## S g1 a4 t5 5 0.9739912 0.02569978 1 4 5
## S g1 a5 t6 6 0.9761105 0.02068439 1 5 6
## S g1 a6 t7 7 0.9777225 0.01731392 1 6 7
## S g1 a7 t8 8 0.9789420 0.01513857 1 7 8
## S g1 a8 t9 9 0.9798542 0.01379832 1 8 9
## S g1 a9 t10 10 0.9805232 0.01301229 1 9 10
## S g1 a10 t11 11 0.9809964 0.01257501 1 10 11
## S g1 a11 t12 12 0.9813093 0.01234702 1 11 12
## S g1 a12 t13 13 0.9814882 0.01223956 1 12 13
## S g1 a13 t14 14 0.9815523 0.01219879 1 13 14
## S g1 a14 t15 15 0.9815153 0.01219339 1 14 15
## S g1 a15 t16 16 0.9813865 0.01220597 1 15 16
## S g1 a16 t17 17 0.9811717 0.01222771 1 16 17
## S g1 a17 t18 18 0.9808733 0.01225520 1 17 18
## S g1 a18 t19 19 0.9804911 0.01228859 1 18 19
## S g1 a19 t20 20 0.9800219 0.01233069 1 19 20
## S g1 a20 t21 21 0.9794594 0.01238654 1 20 21
## S g1 a21 t22 22 0.9787940 0.01246350 1 21 22
## S g1 a22 t23 23 0.9780123 0.01257164 1 22 23
## S g1 a23 t24 24 0.9770957 0.01272481 1 23 24
## S g1 a24 t25 25 0.9760197 0.01294244 1 24 25
## S g1 a25 t26 26 0.9747517 0.01325286 1 25 26
## S g1 a26 t27 27 0.9732484 0.01369915 1 26 27
## S g1 a27 t28 28 0.9714529 0.01434895 1 27 28
## S g1 a28 t29 29 0.9692893 0.01531045 1 28 29
## S g1 a29 t30 30 0.9666569 0.01675582 1 29 30
## S g1 a30 t31 31 0.9634215 0.01895197 1 30 31
## S g1 a31 t32 32 0.9594041 0.02229424 1 31 32
## S g1 a32 t33 33 0.9543674 0.02733818 1 32 33
## S g1 a33 t34 34 0.9479987 0.03482923 1 33 34
## S g1 a34 t35 35 0.9398924 0.04573329 1 34 35
## S g1 a35 t36 36 0.9295359 0.06125596 1 35 36
## S g1 a36 t37 37 0.9163089 0.08280358 1 36 37
## S g1 a37 t38 38 0.8995135 0.11179910 1 37 38
## S g1 a38 t39 39 0.8784623 0.14924532 1 38 39
## Age Time
## S g1 a0 t1 0 0
## S g1 a1 t2 1 1
## S g1 a2 t3 2 2
## S g1 a3 t4 3 3
## S g1 a4 t5 4 4
## S g1 a5 t6 5 5
## S g1 a6 t7 6 6
## S g1 a7 t8 7 7
## S g1 a8 t9 8 8
## S g1 a9 t10 9 9
## S g1 a10 t11 10 10
## S g1 a11 t12 11 11
## S g1 a12 t13 12 12
## S g1 a13 t14 13 13
## S g1 a14 t15 14 14
## S g1 a15 t16 15 15
## S g1 a16 t17 16 16
## S g1 a17 t18 17 17
## S g1 a18 t19 18 18
## S g1 a19 t20 19 19
## S g1 a20 t21 20 20
## S g1 a21 t22 21 21
## S g1 a22 t23 22 22
## S g1 a23 t24 23 23
## S g1 a24 t25 24 24
## S g1 a25 t26 25 25
## S g1 a26 t27 26 26
## S g1 a27 t28 27 27
## S g1 a28 t29 28 28
## S g1 a29 t30 29 29
## S g1 a30 t31 30 30
## S g1 a31 t32 31 31
## S g1 a32 t33 32 32
## S g1 a33 t34 33 33
## S g1 a34 t35 34 34
## S g1 a35 t36 35 35
## S g1 a36 t37 36 36
## S g1 a37 t38 37 37
## S g1 a38 t39 38 38
# plot the estimates from the model plus the model averaging
length(model.set)
## [1] 4
plotdata.indiv <- plyr::ldply(model.set[-length(model.set)], function(x){
est <- get.real(x, "S", se=TRUE)
est$Model <- x$model.name
est
})
plotdata.ma <- S.ma
plotdata.ma$Model <- "Model avg"
# model average the predicted nest suvival probability
final.fit <- ggplot(data=plotdata.indiv, aes(x=time, y=estimate, color=Model))+
ggtitle("Comparing the model fits")+
geom_line(aes(group=Model))+
geom_line(data=plotdata.ma, color='black', size=2, group=1)+
theme(legend.justification=c(0,0), legend.position=c(0,0))
final.fit
ggsave(final.fit,
file=file.path("..","..","..","..","MyStuff","Images","killdeer-modavg.png"),h=4, w=6, units="in", dpi=300)
# model average the derived parameter of the oveall nest survival
model.fits[[1]]$results$derived$"S Overall Survival"
## estimate se lcl ucl
## 1 0.3339134 0.1495833 0.1182906 0.6519547
model.fits[[2]]$results$derived$"S Overall Survival"
## estimate se lcl ucl
## 1 0.3226129 0.1526737 0.1080227 0.651927
model.fits[[3]]$results$derived$"S Overall Survival"
## estimate se lcl ucl
## 1 0.1684153 0.1505241 0.02404565 0.624727
RMark.model.average.derived(model.set, param="S Overall Survival")
## Loading required package: plyr
## Loading required package: boot
## estimate se lcl ucl
## 1 0.2648794 0.1701028 0.06105668 0.6662856
# cleanup
cleanup(ask=FALSE)