Here, you’ll find code for running each of the models used in lab this week with the data for swifts found in the “swift.inp” data set. You’ll also find the model-selection table, model-specific results, and results based on some follow-up work with the output from the top model.
library(RMark)## Loading required package: snowfall
## Loading required package: snow
## This is RMark 2.1.13swift=convert.inp("http://www.montana.edu/rotella/documents/502/AA.inp",
group.df=data.frame(colony=c("Poor","Good")), covariates=NULL)
head(swift)## ch freq colony
## 1:1 11100000 1 Poor
## 1:2 11100000 1 Poor
## 1:3 11001000 1 Poor
## 1:4 11000000 1 Poor
## 1:5 10000000 1 Poor
## 1:6 10000000 1 Pooroptions(width = 150)
# Process data
swift.pr=process.data(swift, begin.time = 1, model = "CJS", groups=("colony"))
#
# Create default design data
swift.ddl=make.design.data(swift.pr)
# Examine the design data. Because the output is so long,
# only the head of each piece of design data is shown
head(swift.ddl$Phi)## par.index model.index group cohort age time occ.cohort Cohort Age Time colony
## 1 1 1 Good 1 0 1 1 0 0 0 Good
## 2 2 2 Good 1 1 2 1 0 1 1 Good
## 3 3 3 Good 1 2 3 1 0 2 2 Good
## 4 4 4 Good 1 3 4 1 0 3 3 Good
## 5 5 5 Good 1 4 5 1 0 4 4 Good
## 6 6 6 Good 1 5 6 1 0 5 5 Goodhead(swift.ddl$p)## par.index model.index group cohort age time occ.cohort Cohort Age Time colony
## 1 1 57 Good 1 1 2 1 0 1 0 Good
## 2 2 58 Good 1 2 3 1 0 2 1 Good
## 3 3 59 Good 1 3 4 1 0 3 2 Good
## 4 4 60 Good 1 4 5 1 0 4 3 Good
## 5 5 61 Good 1 5 6 1 0 5 4 Good
## 6 6 62 Good 1 6 7 1 0 6 5 GoodHere, we set up a function that contains the structures for ‘Phi’ and ‘p’. It will run all combinations of the structures for each parameter type when executed.
run.swift=function()
{
#
# Process data
swift.pr=process.data(swift, begin.time = 1, model = "CJS", groups=("colony"))
#
# Create default design data
#
swift.ddl=make.design.data(swift.pr)
#
# Define range of models for Phi
#
Phi.dot=list(formula=~1)
Phi.time=list(formula=~time)
Phi.colony=list(formula=~colony)
Phi.colonyXtime=list(formula=~colony*time)
#
# Define range of models for p
#
p.dot=list(formula=~1)
p.time=list(formula=~time)
#
# Run all pairings of models for phi & p.
# NOTE: if you do not want to see the output for each model, you can add the text
# ', output=FALSE' after 'ddl=swift.ddl' below. Here, I don't do that so you can
# see the output for each model, but this might not be desired if you have a
# lot of models!
# I added 'adjust=FALSE' below to prevent RMark from adjusting the number of
# parameters, which would be problematic for models like phi(t),p(t) or
# phi(c*t),p(t) in which not all parameters can be estimated.
swift.model.list=create.model.list("CJS")
swift.results=mark.wrapper(swift.model.list,
data=swift.pr,ddl=swift.ddl, adjust=FALSE)
#
# Return model table and list of models
#
return(swift.results)
}It’s very simple to then run the function and each of the models. As implemented here, the output from each of the models appears in the console where it can be reviewed. One can also examine model-specific output in other ways as shown below. And, if you don’t want to see the output from every model appear in the output here, see the note above for the use of the “mark.wrapper” command where one has the option of suppressing the output.
swift.results=run.swift()##
## Phi.colony.p.dot##
## Output summary for CJS model
## Name : Phi(~colony)p(~1)
##
## Npar : 3
## -2lnL: 367.4051
## AICc : 373.5263
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 1.0804050 0.2122218 0.6644503 1.4963597
## Phi:colonyPoor -0.8582301 0.3646689 -1.5729812 -0.1434791
## p:(Intercept) 0.8815436 0.2388674 0.4133634 1.3497238
##
##
## Real Parameter Phi
## Group:colonyGood
## 1 2 3 4 5 6 7
## 1 0.7465706 0.7465706 0.7465706 0.7465706 0.7465706 0.7465706 0.7465706
## 2 0.7465706 0.7465706 0.7465706 0.7465706 0.7465706 0.7465706
## 3 0.7465706 0.7465706 0.7465706 0.7465706 0.7465706
## 4 0.7465706 0.7465706 0.7465706 0.7465706
## 5 0.7465706 0.7465706 0.7465706
## 6 0.7465706 0.7465706
## 7 0.7465706
##
## Group:colonyPoor
## 1 2 3 4 5 6 7
## 1 0.5553164 0.5553164 0.5553164 0.5553164 0.5553164 0.5553164 0.5553164
## 2 0.5553164 0.5553164 0.5553164 0.5553164 0.5553164 0.5553164
## 3 0.5553164 0.5553164 0.5553164 0.5553164 0.5553164
## 4 0.5553164 0.5553164 0.5553164 0.5553164
## 5 0.5553164 0.5553164 0.5553164
## 6 0.5553164 0.5553164
## 7 0.5553164
##
##
## Real Parameter p
## Group:colonyGood
## 2 3 4 5 6 7 8
## 1 0.707142 0.707142 0.707142 0.707142 0.707142 0.707142 0.707142
## 2 0.707142 0.707142 0.707142 0.707142 0.707142 0.707142
## 3 0.707142 0.707142 0.707142 0.707142 0.707142
## 4 0.707142 0.707142 0.707142 0.707142
## 5 0.707142 0.707142 0.707142
## 6 0.707142 0.707142
## 7 0.707142
##
## Group:colonyPoor
## 2 3 4 5 6 7 8
## 1 0.707142 0.707142 0.707142 0.707142 0.707142 0.707142 0.707142
## 2 0.707142 0.707142 0.707142 0.707142 0.707142 0.707142
## 3 0.707142 0.707142 0.707142 0.707142 0.707142
## 4 0.707142 0.707142 0.707142 0.707142
## 5 0.707142 0.707142 0.707142
## 6 0.707142 0.707142
## 7 0.707142##
## Phi.colonyXtime.p.dot##
## Output summary for CJS model
## Name : Phi(~colony * time)p(~1)
##
## Npar : 15
## -2lnL: 353.9155
## AICc : 386.4962
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 4.021170 6.2821380 -8.2918207 16.334161
## Phi:colonyPoor -2.814173 6.3878992 -15.3344560 9.706110
## Phi:time2 -3.023332 6.4380764 -15.6419620 9.595298
## Phi:time3 -3.034088 6.3075930 -15.3969710 9.328794
## Phi:time4 -2.475614 6.3132448 -14.8495740 9.898346
## Phi:time5 -2.892450 6.3035606 -15.2474290 9.462529
## Phi:time6 -2.804021 6.3015697 -15.1550980 9.547056
## Phi:time7 -4.090684 6.2802760 -16.4000250 8.218657
## Phi:colonyPoor:time2 2.045401 6.6067252 -10.9037810 14.994582
## Phi:colonyPoor:time3 1.844732 6.4828492 -10.8616530 14.551117
## Phi:colonyPoor:time4 1.060909 6.4772950 -11.6345890 13.756408
## Phi:colonyPoor:time5 2.101086 6.4926144 -10.6244380 14.826611
## Phi:colonyPoor:time6 3.717328 7.2038465 -10.4022110 17.836867
## Phi:colonyPoor:time7 1.779563 6.4538965 -10.8700750 14.429200
## p:(Intercept) 1.070237 0.2539357 0.5725229 1.567951
##
##
## Real Parameter Phi
## Group:colonyGood
## 1 2 3 4 5 6 7
## 1 0.9823839 0.7306333 0.7285111 0.824271 0.7556027 0.7715614 0.4826285
## 2 0.7306333 0.7285111 0.824271 0.7556027 0.7715614 0.4826285
## 3 0.7285111 0.824271 0.7556027 0.7715614 0.4826285
## 4 0.824271 0.7556027 0.7715614 0.4826285
## 5 0.7556027 0.7715614 0.4826285
## 6 0.7715614 0.4826285
## 7 0.4826285
##
## Group:colonyPoor
## 1 2 3 4 5 6 7
## 1 0.7697672 0.5570173 0.5044101 0.448259 0.6024379 0.892861 0.2489679
## 2 0.5570173 0.5044101 0.448259 0.6024379 0.892861 0.2489679
## 3 0.5044101 0.448259 0.6024379 0.892861 0.2489679
## 4 0.448259 0.6024379 0.892861 0.2489679
## 5 0.6024379 0.892861 0.2489679
## 6 0.892861 0.2489679
## 7 0.2489679
##
##
## Real Parameter p
## Group:colonyGood
## 2 3 4 5 6 7 8
## 1 0.744642 0.744642 0.744642 0.744642 0.744642 0.744642 0.744642
## 2 0.744642 0.744642 0.744642 0.744642 0.744642 0.744642
## 3 0.744642 0.744642 0.744642 0.744642 0.744642
## 4 0.744642 0.744642 0.744642 0.744642
## 5 0.744642 0.744642 0.744642
## 6 0.744642 0.744642
## 7 0.744642
##
## Group:colonyPoor
## 2 3 4 5 6 7 8
## 1 0.744642 0.744642 0.744642 0.744642 0.744642 0.744642 0.744642
## 2 0.744642 0.744642 0.744642 0.744642 0.744642 0.744642
## 3 0.744642 0.744642 0.744642 0.744642 0.744642
## 4 0.744642 0.744642 0.744642 0.744642
## 5 0.744642 0.744642 0.744642
## 6 0.744642 0.744642
## 7 0.744642##
## Phi.dot.p.dot##
## Output summary for CJS model
## Name : Phi(~1)p(~1)
##
## Npar : 2
## -2lnL: 372.8533
## AICc : 376.9136
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 0.8524384 0.1753794 0.5086948 1.196182
## p:(Intercept) 0.8881232 0.2391869 0.4193169 1.356929
##
##
## Real Parameter Phi
## Group:colonyGood
## 1 2 3 4 5 6 7
## 1 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784
## 2 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784
## 3 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784
## 4 0.7010784 0.7010784 0.7010784 0.7010784
## 5 0.7010784 0.7010784 0.7010784
## 6 0.7010784 0.7010784
## 7 0.7010784
##
## Group:colonyPoor
## 1 2 3 4 5 6 7
## 1 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784
## 2 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784
## 3 0.7010784 0.7010784 0.7010784 0.7010784 0.7010784
## 4 0.7010784 0.7010784 0.7010784 0.7010784
## 5 0.7010784 0.7010784 0.7010784
## 6 0.7010784 0.7010784
## 7 0.7010784
##
##
## Real Parameter p
## Group:colonyGood
## 2 3 4 5 6 7 8
## 1 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027
## 2 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027
## 3 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027
## 4 0.7085027 0.7085027 0.7085027 0.7085027
## 5 0.7085027 0.7085027 0.7085027
## 6 0.7085027 0.7085027
## 7 0.7085027
##
## Group:colonyPoor
## 2 3 4 5 6 7 8
## 1 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027
## 2 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027
## 3 0.7085027 0.7085027 0.7085027 0.7085027 0.7085027
## 4 0.7085027 0.7085027 0.7085027 0.7085027
## 5 0.7085027 0.7085027 0.7085027
## 6 0.7085027 0.7085027
## 7 0.7085027##
## Phi.time.p.dot##
## Output summary for CJS model
## Name : Phi(~time)p(~1)
##
## Npar : 8
## -2lnL: 362.2662
## AICc : 379.0123
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 2.2494541 1.3060918 -0.3104860 4.8093941
## Phi:time2 -1.5542402 1.4417055 -4.3799830 1.2715027
## Phi:time3 -1.4773075 1.3928559 -4.2073051 1.2526901
## Phi:time4 -1.2296931 1.3935547 -3.9610603 1.5016741
## Phi:time5 -1.2532870 1.3949524 -3.9873937 1.4808196
## Phi:time6 -0.9241435 1.4512521 -3.7685977 1.9203107
## Phi:time7 -2.5025972 1.3406122 -5.1301971 0.1250028
## p:(Intercept) 1.0638310 0.2537878 0.5664069 1.5612551
##
##
## Real Parameter Phi
## Group:colonyGood
## 1 2 3 4 5 6 7
## 1 0.9046034 0.6671258 0.683985 0.734926 0.7303043 0.7900639 0.43705
## 2 0.6671258 0.683985 0.734926 0.7303043 0.7900639 0.43705
## 3 0.683985 0.734926 0.7303043 0.7900639 0.43705
## 4 0.734926 0.7303043 0.7900639 0.43705
## 5 0.7303043 0.7900639 0.43705
## 6 0.7900639 0.43705
## 7 0.43705
##
## Group:colonyPoor
## 1 2 3 4 5 6 7
## 1 0.9046034 0.6671258 0.683985 0.734926 0.7303043 0.7900639 0.43705
## 2 0.6671258 0.683985 0.734926 0.7303043 0.7900639 0.43705
## 3 0.683985 0.734926 0.7303043 0.7900639 0.43705
## 4 0.734926 0.7303043 0.7900639 0.43705
## 5 0.7303043 0.7900639 0.43705
## 6 0.7900639 0.43705
## 7 0.43705
##
##
## Real Parameter p
## Group:colonyGood
## 2 3 4 5 6 7 8
## 1 0.743422 0.743422 0.743422 0.743422 0.743422 0.743422 0.743422
## 2 0.743422 0.743422 0.743422 0.743422 0.743422 0.743422
## 3 0.743422 0.743422 0.743422 0.743422 0.743422
## 4 0.743422 0.743422 0.743422 0.743422
## 5 0.743422 0.743422 0.743422
## 6 0.743422 0.743422
## 7 0.743422
##
## Group:colonyPoor
## 2 3 4 5 6 7 8
## 1 0.743422 0.743422 0.743422 0.743422 0.743422 0.743422 0.743422
## 2 0.743422 0.743422 0.743422 0.743422 0.743422 0.743422
## 3 0.743422 0.743422 0.743422 0.743422 0.743422
## 4 0.743422 0.743422 0.743422 0.743422
## 5 0.743422 0.743422 0.743422
## 6 0.743422 0.743422
## 7 0.743422##
## Phi.colony.p.time##
## Output summary for CJS model
## Name : Phi(~colony)p(~time)
##
## Npar : 9
## -2lnL: 350.8705
## AICc : 369.808
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 1.2080437 0.2197772 0.7772804 1.6388071
## Phi:colonyPoor -0.8973285 0.3759843 -1.6342577 -0.1603994
## p:(Intercept) 2.3000713 1.0331517 0.2750940 4.3250486
## p:time3 -1.3106886 1.1495110 -3.5637301 0.9423529
## p:time4 -2.1525277 1.1281043 -4.3636122 0.0585567
## p:time5 -1.4604238 1.1436593 -3.7019960 0.7811485
## p:time6 -0.4981633 1.2596760 -2.9671283 1.9708017
## p:time7 -0.4818736 1.3516499 -3.1311075 2.1673603
## p:time8 -2.4471003 1.0999952 -4.6030910 -0.2911097
##
##
## Real Parameter Phi
## Group:colonyGood
## 1 2 3 4 5 6 7
## 1 0.7699526 0.7699526 0.7699526 0.7699526 0.7699526 0.7699526 0.7699526
## 2 0.7699526 0.7699526 0.7699526 0.7699526 0.7699526 0.7699526
## 3 0.7699526 0.7699526 0.7699526 0.7699526 0.7699526
## 4 0.7699526 0.7699526 0.7699526 0.7699526
## 5 0.7699526 0.7699526 0.7699526
## 6 0.7699526 0.7699526
## 7 0.7699526
##
## Group:colonyPoor
## 1 2 3 4 5 6 7
## 1 0.5770598 0.5770598 0.5770598 0.5770598 0.5770598 0.5770598 0.5770598
## 2 0.5770598 0.5770598 0.5770598 0.5770598 0.5770598 0.5770598
## 3 0.5770598 0.5770598 0.5770598 0.5770598 0.5770598
## 4 0.5770598 0.5770598 0.5770598 0.5770598
## 5 0.5770598 0.5770598 0.5770598
## 6 0.5770598 0.5770598
## 7 0.5770598
##
##
## Real Parameter p
## Group:colonyGood
## 2 3 4 5 6 7 8
## 1 0.9088829 0.728966 0.5368191 0.698391 0.858381 0.8603497 0.4633088
## 2 0.728966 0.5368191 0.698391 0.858381 0.8603497 0.4633088
## 3 0.5368191 0.698391 0.858381 0.8603497 0.4633088
## 4 0.698391 0.858381 0.8603497 0.4633088
## 5 0.858381 0.8603497 0.4633088
## 6 0.8603497 0.4633088
## 7 0.4633088
##
## Group:colonyPoor
## 2 3 4 5 6 7 8
## 1 0.9088829 0.728966 0.5368191 0.698391 0.858381 0.8603497 0.4633088
## 2 0.728966 0.5368191 0.698391 0.858381 0.8603497 0.4633088
## 3 0.5368191 0.698391 0.858381 0.8603497 0.4633088
## 4 0.698391 0.858381 0.8603497 0.4633088
## 5 0.858381 0.8603497 0.4633088
## 6 0.8603497 0.4633088
## 7 0.4633088##
## Phi.colonyXtime.p.time##
## Output summary for CJS model
## Name : Phi(~colony * time)p(~time)
##
## Npar : 20
## -2lnL: 346.7694
## AICc : 391.4103
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 2.6603497 1.718601 -0.7081085 6.0288079
## Phi:colonyPoor -1.7682770 1.937003 -5.5648027 2.0282487
## Phi:time2 -1.5645991 1.898831 -5.2863087 2.1571104
## Phi:time3 -1.3615838 1.904035 -5.0934934 2.3703257
## Phi:time4 -1.2862678 1.891571 -4.9937478 2.4212123
## Phi:time5 -1.7862068 1.794655 -5.3037311 1.7313175
## Phi:time6 -1.6786662 1.819408 -5.2447054 1.8873730
## Phi:time7 -1.5511971 0.000000 -1.5511971 -1.5511971
## Phi:colonyPoor:time2 0.9922982 2.215302 -3.3496940 5.3342903
## Phi:colonyPoor:time3 0.8966518 2.420498 -3.8475242 5.6408278
## Phi:colonyPoor:time4 0.0736330 2.246614 -4.3297299 4.4769958
## Phi:colonyPoor:time5 1.1602933 2.217577 -3.1861587 5.5067453
## Phi:colonyPoor:time6 2.2518374 2.616936 -2.8773579 7.3810328
## Phi:colonyPoor:time7 0.1978668 0.000000 0.1978668 0.1978668
## p:(Intercept) 2.0134253 1.049826 -0.0442329 4.0710836
## p:time3 -1.0096340 1.198442 -3.3585806 1.3393125
## p:time4 -1.9040432 1.162949 -4.1834231 0.3753368
## p:time5 -1.2090470 1.175206 -3.5124513 1.0943572
## p:time6 -0.0983438 1.291917 -2.6305013 2.4338137
## p:time7 -0.0687085 1.480364 -2.9702222 2.8328052
## p:time8 -2.0119669 0.000000 -2.0119669 -2.0119669
##
##
## Real Parameter Phi
## Group:colonyGood
## 1 2 3 4 5 6 7
## 1 0.934646 0.749463 0.7856272 0.7980388 0.705607 0.7274421 0.7519711
## 2 0.749463 0.7856272 0.7980388 0.705607 0.7274421 0.7519711
## 3 0.7856272 0.7980388 0.705607 0.7274421 0.7519711
## 4 0.7980388 0.705607 0.7274421 0.7519711
## 5 0.705607 0.7274421 0.7519711
## 6 0.7274421 0.7519711
## 7 0.7519711
##
## Group:colonyPoor
## 1 2 3 4 5 6 7
## 1 0.7093177 0.5792686 0.6051907 0.4205388 0.5661498 0.8123334 0.3866875
## 2 0.5792686 0.6051907 0.4205388 0.5661498 0.8123334 0.3866875
## 3 0.6051907 0.4205388 0.5661498 0.8123334 0.3866875
## 4 0.4205388 0.5661498 0.8123334 0.3866875
## 5 0.5661498 0.8123334 0.3866875
## 6 0.8123334 0.3866875
## 7 0.3866875
##
##
## Real Parameter p
## Group:colonyGood
## 2 3 4 5 6 7 8
## 1 0.8821995 0.7318033 0.5273183 0.6909103 0.871589 0.8748694 0.5003646
## 2 0.7318033 0.5273183 0.6909103 0.871589 0.8748694 0.5003646
## 3 0.5273183 0.6909103 0.871589 0.8748694 0.5003646
## 4 0.6909103 0.871589 0.8748694 0.5003646
## 5 0.871589 0.8748694 0.5003646
## 6 0.8748694 0.5003646
## 7 0.5003646
##
## Group:colonyPoor
## 2 3 4 5 6 7 8
## 1 0.8821995 0.7318033 0.5273183 0.6909103 0.871589 0.8748694 0.5003646
## 2 0.7318033 0.5273183 0.6909103 0.871589 0.8748694 0.5003646
## 3 0.5273183 0.6909103 0.871589 0.8748694 0.5003646
## 4 0.6909103 0.871589 0.8748694 0.5003646
## 5 0.871589 0.8748694 0.5003646
## 6 0.8748694 0.5003646
## 7 0.5003646##
## Phi.dot.p.time##
## Output summary for CJS model
## Name : Phi(~1)p(~time)
##
## Npar : 8
## -2lnL: 356.4878
## AICc : 373.2339
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 0.9609298 0.1778847 0.6122757 1.3095839
## p:(Intercept) 2.2693089 1.0302722 0.2499753 4.2886425
## p:time3 -1.3175779 1.1430662 -3.5579878 0.9228319
## p:time4 -2.1404247 1.1245544 -4.3445513 0.0637019
## p:time5 -1.4506994 1.1401739 -3.6854402 0.7840414
## p:time6 -0.4538378 1.2575182 -2.9185736 2.0108980
## p:time7 -0.2969835 1.3771912 -2.9962783 2.4023113
## p:time8 -2.3835359 1.0981279 -4.5358666 -0.2312051
##
##
## Real Parameter Phi
## Group:colonyGood
## 1 2 3 4 5 6 7
## 1 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079
## 2 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079
## 3 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079
## 4 0.7233079 0.7233079 0.7233079 0.7233079
## 5 0.7233079 0.7233079 0.7233079
## 6 0.7233079 0.7233079
## 7 0.7233079
##
## Group:colonyPoor
## 1 2 3 4 5 6 7
## 1 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079
## 2 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079
## 3 0.7233079 0.7233079 0.7233079 0.7233079 0.7233079
## 4 0.7233079 0.7233079 0.7233079 0.7233079
## 5 0.7233079 0.7233079 0.7233079
## 6 0.7233079 0.7233079
## 7 0.7233079
##
##
## Real Parameter p
## Group:colonyGood
## 2 3 4 5 6 7 8
## 1 0.9063031 0.7214632 0.5321765 0.6939411 0.8600218 0.8778607 0.4714743
## 2 0.7214632 0.5321765 0.6939411 0.8600218 0.8778607 0.4714743
## 3 0.5321765 0.6939411 0.8600218 0.8778607 0.4714743
## 4 0.6939411 0.8600218 0.8778607 0.4714743
## 5 0.8600218 0.8778607 0.4714743
## 6 0.8778607 0.4714743
## 7 0.4714743
##
## Group:colonyPoor
## 2 3 4 5 6 7 8
## 1 0.9063031 0.7214632 0.5321765 0.6939411 0.8600218 0.8778607 0.4714743
## 2 0.7214632 0.5321765 0.6939411 0.8600218 0.8778607 0.4714743
## 3 0.5321765 0.6939411 0.8600218 0.8778607 0.4714743
## 4 0.6939411 0.8600218 0.8778607 0.4714743
## 5 0.8600218 0.8778607 0.4714743
## 6 0.8778607 0.4714743
## 7 0.4714743##
## Phi.time.p.time##
## Output summary for CJS model
## Name : Phi(~time)p(~time)
##
## Npar : 13
## -2lnL: 354.9445
## AICc : 382.8807
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 1.7439701 0.8654835 0.0476224 3.4403179
## Phi:time2 -0.9670002 1.0306456 -2.9870656 1.0530653
## Phi:time3 -0.5738986 1.1624636 -2.8523273 1.7045301
## Phi:time4 -0.8957171 1.0338528 -2.9220687 1.1306344
## Phi:time5 -0.9809819 0.9802253 -2.9022235 0.9402596
## Phi:time6 -0.6912527 1.0551054 -2.7592593 1.3767538
## Phi:time7 -1.8256772 0.0000000 -1.8256772 -1.8256772
## p:(Intercept) 2.0030687 1.0495306 -0.0540113 4.0601486
## p:time3 -0.9689947 1.1966902 -3.3145075 1.3765182
## p:time4 -1.9340756 1.1630582 -4.2136697 0.3455185
## p:time5 -1.2041766 1.1750316 -3.5072387 1.0988854
## p:time6 -0.0882494 1.2916736 -2.6199297 2.4434309
## p:time7 -0.0861456 1.4799681 -2.9868832 2.8145920
## p:time8 -1.1127934 0.0000000 -1.1127934 -1.1127934
##
##
## Real Parameter Phi
## Group:colonyGood
## 1 2 3 4 5 6 7
## 1 0.8511906 0.6850267 0.7631579 0.7002005 0.6820022 0.7412964 0.4795846
## 2 0.6850267 0.7631579 0.7002005 0.6820022 0.7412964 0.4795846
## 3 0.7631579 0.7002005 0.6820022 0.7412964 0.4795846
## 4 0.7002005 0.6820022 0.7412964 0.4795846
## 5 0.6820022 0.7412964 0.4795846
## 6 0.7412964 0.4795846
## 7 0.4795846
##
## Group:colonyPoor
## 1 2 3 4 5 6 7
## 1 0.8511906 0.6850267 0.7631579 0.7002005 0.6820022 0.7412964 0.4795846
## 2 0.6850267 0.7631579 0.7002005 0.6820022 0.7412964 0.4795846
## 3 0.7631579 0.7002005 0.6820022 0.7412964 0.4795846
## 4 0.7002005 0.6820022 0.7412964 0.4795846
## 5 0.6820022 0.7412964 0.4795846
## 6 0.7412964 0.4795846
## 7 0.4795846
##
##
## Real Parameter p
## Group:colonyGood
## 2 3 4 5 6 7 8
## 1 0.8811189 0.737705 0.5172414 0.6897374 0.8715596 0.8717949 0.708947
## 2 0.737705 0.5172414 0.6897374 0.8715596 0.8717949 0.708947
## 3 0.5172414 0.6897374 0.8715596 0.8717949 0.708947
## 4 0.6897374 0.8715596 0.8717949 0.708947
## 5 0.8715596 0.8717949 0.708947
## 6 0.8717949 0.708947
## 7 0.708947
##
## Group:colonyPoor
## 2 3 4 5 6 7 8
## 1 0.8811189 0.737705 0.5172414 0.6897374 0.8715596 0.8717949 0.708947
## 2 0.737705 0.5172414 0.6897374 0.8715596 0.8717949 0.708947
## 3 0.5172414 0.6897374 0.8715596 0.8717949 0.708947
## 4 0.6897374 0.8715596 0.8717949 0.708947
## 5 0.8715596 0.8717949 0.708947
## 6 0.8717949 0.708947
## 7 0.708947Once the models are run, we can examine the model-selection table. We can also examine model-specific output.
options(width = 150)
swift.results## model npar AICc DeltaAICc weight Deviance
## 2 Phi(~colony)p(~time) 9 369.8080 0.000000 7.264694e-01 111.6644
## 6 Phi(~1)p(~time) 8 373.2339 3.425894 1.310068e-01 117.2816
## 1 Phi(~colony)p(~1) 3 373.5263 3.718302 1.131874e-01 128.1989
## 5 Phi(~1)p(~1) 2 376.9136 7.105612 2.080910e-02 133.6472
## 7 Phi(~time)p(~1) 8 379.0123 9.204334 7.286544e-03 123.0601
## 8 Phi(~time)p(~time) 13 382.8807 13.072740 1.053193e-03 115.7384
## 3 Phi(~colony * time)p(~1) 15 386.4962 16.688205 1.727506e-04 114.7094
## 4 Phi(~colony * time)p(~time) 20 391.4103 21.602314 1.480248e-05 107.5633names(swift.results)## [1] "Phi.colony.p.dot" "Phi.colony.p.time" "Phi.colonyXtime.p.dot" "Phi.colonyXtime.p.time" "Phi.dot.p.dot"
## [6] "Phi.dot.p.time" "Phi.time.p.dot" "Phi.time.p.time" "model.table"# examine the output from top-ranked model (model #2) and
# store top model in object with short name
top=swift.results$Phi.colony.p.time
# look at output
top$results$beta## estimate se lcl ucl
## Phi:(Intercept) 1.2080437 0.2197772 0.7772804 1.6388071
## Phi:colonyPoor -0.8973285 0.3759843 -1.6342577 -0.1603994
## p:(Intercept) 2.3000713 1.0331517 0.2750940 4.3250486
## p:time3 -1.3106886 1.1495110 -3.5637301 0.9423529
## p:time4 -2.1525277 1.1281043 -4.3636122 0.0585567
## p:time5 -1.4604238 1.1436593 -3.7019960 0.7811485
## p:time6 -0.4981633 1.2596760 -2.9671283 1.9708017
## p:time7 -0.4818736 1.3516499 -3.1311075 2.1673603
## p:time8 -2.4471003 1.0999952 -4.6030910 -0.2911097top$results$real[c(1,2),]## estimate se lcl ucl fixed note
## Phi gGood c1 a0 t1 0.7699526 0.0389282 0.6850937 0.8373726
## Phi gPoor c1 a0 t1 0.5770598 0.0771525 0.4233888 0.7171377Here, we adjust using an estimate of 1.035414. To estimate this quantity, you’d need to work in MARK directly. Notice that these SEs have been inflated by the square root of 1.035414.
options(width = 150)
QAICc.table=adjust.chat(1.035414, swift.results)
QAICc.table## model npar QAICc DeltaQAICc weight QDeviance chat
## 2 Phi(~colony)p(~time) 9 357.8072 0.000000 6.844790e-01 107.8451 1.035414
## 1 Phi(~colony)p(~1) 3 360.9600 3.152774 1.414960e-01 123.8142 1.035414
## 6 Phi(~1)p(~time) 8 361.0410 3.233768 1.358803e-01 113.2703 1.035414
## 5 Phi(~1)p(~1) 2 364.1610 6.353739 2.855376e-02 129.0761 1.035414
## 7 Phi(~time)p(~1) 8 366.6218 8.814569 8.342588e-03 118.8511 1.035414
## 8 Phi(~time)p(~time) 13 370.7406 12.933396 1.063920e-03 111.7799 1.035414
## 3 Phi(~colony * time)p(~1) 15 374.3913 16.584056 1.714661e-04 110.7860 1.035414
## 4 Phi(~colony * time)p(~time) 20 379.5498 21.742582 1.300225e-05 103.8843 1.035414# store top model in object with short name
top.qaic=QAICc.table$Phi.colony.p.time
# look at estimates of real parameters and notice that the SEs have been adjusted
phi.hat.qaic=get.real(top.qaic, "Phi", se=TRUE)
# only need to print out rows 1 and 29 to get estimate for each colony
phi.hat.qaic[c(1,29),c(1:6)]## all.diff.index par.index estimate se lcl ucl
## Phi gGood c1 a0 t1 1 1 0.7699526 0.03961147 0.6834601 0.8383996
## Phi gPoor c1 a0 t1 29 2 0.5770598 0.07850670 0.4207360 0.7193386# compare with results obtained from AICc
top$results$real[c(1,2),1:6]## estimate se lcl ucl fixed note
## Phi gGood c1 a0 t1 0.7699526 0.0389282 0.6850937 0.8373726
## Phi gPoor c1 a0 t1 0.5770598 0.0771525 0.4233888 0.7171377In lab, we were interested in estimating the difference in survival rates for the 2 colonies and putting a confidence interval on that estimated difference. Here, we use the output from the model and the delta method to develop the SE and 95% CI for expressing the uncertainty for the estimated difference. First, the code uses the ‘deltamethod’ function in the ‘msm’ package to calculate the SE for the estimated difference. If you don’t have the package ‘msm’ and you’re trying to use the code below, you’ll need to install the package and then load it before executing the code below. The code also shows the actual (& algebraically simple) calculation for estimating that SE based on the variance and covariance of the estimated survival rates for the 2 colonies. For the variance of a difference, we sum the variances and subtract 2 times the covariance.
# store var-cov for real estimates
rr=get.real(top.qaic, "Phi", se=TRUE,vcv=TRUE)
sigma=rr$vcv.real # store var-cov in "sigma"
# load package 'msm' that implements the delta method
# then run delta method on difference in survival
library(msm)
s.good = rr$estimates$estimate[1] # first row w/ output for good colony
s.poor = rr$estimates$estimate[29] # first row w/ output for poor colony
Diff= s.good - s.poor
seDiff=deltamethod(~x1-x2,c(s.good,s.poor),sigma)
lclDiff=Diff-1.96*seDiff
uclDiff=Diff+1.96*seDiff
round(c(Diff,seDiff,lclDiff,uclDiff),3)## [1] 0.193 0.086 0.024 0.362# compare SE to that obtained by direct calculation
# where: variance of difference = var(A)+var(B)-2*cov(A,B)
(se.diff.direct=sqrt(sum(diag(sigma))-2*sigma[1,2]))## [1] 0.08616798# cleanup files created by MARK with the 2 lines below
#rm(list=ls(all=TRUE))
#cleanup(ask = FALSE)