Here, we’ll run the single-season occupancy models for data on Weta that you analyzed directly in MARK in lab 10. Here, we’re dealing with 2 types of parameters: ‘Psi’ & ‘p’.
library(RMark)## This is RMark 2.2.4
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip# you can ignore warning about final line issue that might appear
weta <- convert.inp("http://www.montana.edu/rotella/documents/502/wetaTwoGroups.inp",
group.df = data.frame(Browse = c("Browsed", "Unbrowsed")),
covariates = c(
"obs11", "obs12", "obs13", "obs14", "obs15",
"obs21", "obs22", "obs23", "obs24", "obs25"))## Warning in readLines(inp.filename): incomplete final line found on 'http://
## www.montana.edu/rotella/documents/502/wetaTwoGroups.inp'# store the number of sites
n.sites <- length(weta$ch)
# store the maximum number of visits
n.visits <- max(nchar(weta$ch))
# make a factor that records whether observer 1, 2, or 3 visited a site
# on each visit. A '9' was used to indicate that no observer visited the site.
# So, if there's a 9 for either observer 1 or 2 on a given date, we
# assign a '.' for that day
obs <- matrix(NA, n.sites, n.visits)
for (i in 1:n.sites) {
obs[i, 1] = ifelse(weta$obs11[i] == 1, 1,
ifelse(weta$obs21[i] == 1, 2,
ifelse(weta$obs11[i] == 9, ".", 3)))
obs[i, 2] = ifelse(weta$obs12[i] == 1, 1,
ifelse(weta$obs22[i] == 1, 2,
ifelse(weta$obs12[i] == 9, ".", 3)))
obs[i, 3] = ifelse(weta$obs13[i] == 1, 1,
ifelse(weta$obs23[i] == 1, 2,
ifelse(weta$obs13[i] == 9, ".", 3)))
obs[i, 4] = ifelse(weta$obs14[i] == 1, 1,
ifelse(weta$obs24[i] == 1, 2,
ifelse(weta$obs14[i] == 9, ".", 3)))
obs[i, 5] = ifelse(weta$obs15[i] == 1, 1,
ifelse(weta$obs25[i] == 1, 2,
ifelse(weta$obs15[i] == 9, ".", 3)))
}
head(weta)## ch freq Browse obs11 obs12 obs13 obs14 obs15 obs21 obs22 obs23
## 1:1 0000. 1 Browsed 1 0 0 0 9 0 0 1
## 1:2 0000. 1 Browsed 1 0 0 0 9 0 0 1
## 1:3 0001. 1 Browsed 1 0 0 0 9 0 0 1
## 1:5 0000. 1 Browsed 1 0 0 0 9 0 0 1
## 1:10 1100. 1 Browsed 1 0 0 0 9 0 0 1
## 1:11 00.10 1 Browsed 1 0 9 0 1 0 0 9
## obs24 obs25
## 1:1 0 9
## 1:2 0 9
## 1:3 0 9
## 1:5 0 9
## 1:10 0 9
## 1:11 0 0head(obs)## [,1] [,2] [,3] [,4] [,5]
## [1,] "1" "3" "2" "3" "."
## [2,] "1" "3" "2" "3" "."
## [3,] "1" "3" "2" "3" "."
## [4,] "1" "3" "2" "3" "."
## [5,] "1" "3" "2" "3" "."
## [6,] "1" "3" "." "3" "1"# can drop obs11, obs12, ..., obs25 from weta since not used
weta <- weta[,-c(4:13)]
# add 'Obs' values for each day to 'weta'
weta$Obs1 = as.factor(obs[, 1])
weta$Obs2 = as.factor(obs[, 2])
weta$Obs3 = as.factor(obs[, 3])
weta$Obs4 = as.factor(obs[, 4])
weta$Obs5 = as.factor(obs[, 5])
head(weta)## ch freq Browse Obs1 Obs2 Obs3 Obs4 Obs5
## 1:1 0000. 1 Browsed 1 3 2 3 .
## 1:2 0000. 1 Browsed 1 3 2 3 .
## 1:3 0001. 1 Browsed 1 3 2 3 .
## 1:5 0000. 1 Browsed 1 3 2 3 .
## 1:10 1100. 1 Browsed 1 3 2 3 .
## 1:11 00.10 1 Browsed 1 3 . 3 1First, the function puts the observer and time covariates into proper form for use in modeling and processes the data. Second, thefunction decleares the structures for ‘p’ and ‘Psi’ that will be used in all-possible combinations of the structures for each parameter type when the function is executed.
# Create function to fit the 18 models in the book
fit.weta.models <- function() {
# use 'make.time.factor' to create time-varying dummy variables Obs1 and Obs2
# observer 3 is used as the intercept to match what we did in MARK
weta = make.time.factor(weta, "Obs", 1:5, intercept = 3)
# Process data
weta.process = process.data(weta, model = "Occupancy", groups = "Browse")
weta.ddl = make.design.data(weta.process)
# time factor variable copied to 'Day' to match name used in the lab
weta.ddl$p$Day = weta.ddl$p$time
# Define p models
p.dot = list(formula = ~ 1)
p.browse = list(formula = ~ Browse)
p.day = list(formula = ~ Day)
p.obs = list(formula = ~ Obs1 + Obs2)
p.day.obs = list(formula = ~ Day + Obs1 + Obs2)
p.day.browse = list(formula = ~ Day + Browse)
p.obs.browse = list(formula = ~ Obs1 + Obs2 + Browse)
p.day.obs.browse = list(formula = ~ Day + Obs1 + Obs2 + Browse)
# Define Psi models
Psi.dot = list(formula = ~ 1)
Psi.browse = list(formula = ~ Browse)
# Create model list
cml = create.model.list("Occupancy")
# Run and return marklist of models
return(mark.wrapper(cml, data = weta.process, ddl = weta.ddl))
}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 when the function written above is called. If you don’t want the output from every model to appear automatically, see the note above the “mark.wrapper” command for how to set the option of suppressing the output. One can also examine model-specific output in other ways as shown below.
weta.models <- fit.weta.models()##
## p.browse.Psi.browse##
## Output summary for Occupancy model
## Name : p(~Browse)Psi(~Browse)
##
## Npar : 4
## -2lnL: 258.2606
## AICc : 266.8576
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6107331 0.3022127 -1.2030700 -0.0183963
## p:BrowseUnbrowsed -0.0296221 0.4859581 -0.9821000 0.9228558
## Psi:(Intercept) 1.1336537 0.6935108 -0.2256276 2.4929350
## Psi:BrowseUnbrowsed -1.1980364 0.8416136 -2.8475991 0.4515262
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3518920 0.3518920 0.3518920 0.3518920 0.3518920
## Group:BrowseUnbrowsed 0.3451663 0.3451663 0.3451663 0.3451663 0.3451663
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.7565125
## Group:BrowseUnbrowsed 0.4839099##
## p.day.Psi.browse##
## Output summary for Occupancy model
## Name : p(~Day)Psi(~Browse)
##
## Npar : 7
## -2lnL: 245.4368
## AICc : 261.1868
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6100131 0.4312907 -1.4553430 0.2353168
## p:Day2 -0.1551354 0.5283107 -1.1906245 0.8803536
## p:Day3 -0.9792832 0.5881659 -2.1320883 0.1735219
## p:Day4 -0.1827222 0.5424144 -1.2458546 0.8804101
## p:Day5 0.9811154 0.5456452 -0.0883492 2.0505800
## Psi:(Intercept) 1.2078314 0.6959862 -0.1563015 2.5719643
## Psi:BrowseUnbrowsed -1.2351491 0.7508901 -2.7068938 0.2365956
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3520562 0.3175295 0.1694829 0.3115816 0.5917253
## Group:BrowseUnbrowsed 0.3520562 0.3175295 0.1694829 0.3115816 0.5917253
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.769915
## Group:BrowseUnbrowsed 0.493171##
## p.day.browse.Psi.browse##
## Output summary for Occupancy model
## Name : p(~Day + Browse)Psi(~Browse)
##
## Npar : 8
## -2lnL: 245.408
## AICc : 263.6937
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.5832511 0.4579072 -1.4807492 0.3142470
## p:Day2 -0.1542413 0.5281840 -1.1894819 0.8809993
## p:Day3 -0.9689625 0.5909318 -2.1271889 0.1892639
## p:Day4 -0.1738073 0.5446955 -1.2414105 0.8937958
## p:Day5 0.9942177 0.5513073 -0.0863446 2.0747801
## p:BrowseUnbrowsed -0.0882306 0.5206885 -1.1087802 0.9323189
## Psi:(Intercept) 1.1623669 0.7146130 -0.2382746 2.5630083
## Psi:BrowseUnbrowsed -1.1516424 0.8802055 -2.8768452 0.5735604
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3581849 0.3235527 0.1747668 0.3192853 0.6013196
## Group:BrowseUnbrowsed 0.3381651 0.3045502 0.1624046 0.3004220 0.5799909
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.7617625
## Group:BrowseUnbrowsed 0.5026811##
## p.day.obs.Psi.browse##
## Output summary for Occupancy model
## Name : p(~Day + Obs1 + Obs2)Psi(~Browse)
##
## Npar : 9
## -2lnL: 239.5994
## AICc : 260.5027
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.2297010 0.4732890 -1.1573474 0.6979454
## p:Day2 -0.1548925 0.5504441 -1.2337630 0.9239780
## p:Day3 -0.9431918 0.5997284 -2.1186595 0.2322759
## p:Day4 -0.0706551 0.5629219 -1.1739820 1.0326718
## p:Day5 1.0354542 0.5639409 -0.0698700 2.1407784
## p:Obs1 -1.0700177 0.4601063 -1.9718261 -0.1682093
## p:Obs2 -0.3435561 0.4357548 -1.1976355 0.5105234
## Psi:(Intercept) 1.1933135 0.6757956 -0.1312459 2.5178729
## Psi:BrowseUnbrowsed -1.1693203 0.7412937 -2.6222560 0.2836154
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3545613 0.3246093 0.1728514 0.3453743 0.626219
## Group:BrowseUnbrowsed 0.3545613 0.3246093 0.1728514 0.3453743 0.626219
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.7673332
## Group:BrowseUnbrowsed 0.5059980##
## p.day.obs.browse.Psi.browse##
## Output summary for Occupancy model
## Name : p(~Day + Obs1 + Obs2 + Browse)Psi(~Browse)
##
## Npar : 10
## -2lnL: 239.5993
## AICc : 263.2059
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.2312156 0.4977480 -1.2068017 0.7443704
## p:Day2 -0.1547383 0.5506911 -1.2340928 0.9246162
## p:Day3 -0.9437800 0.6027615 -2.1251925 0.2376325
## p:Day4 -0.0708173 0.5631909 -1.1746714 1.0330367
## p:Day5 1.0348547 0.5672487 -0.0769527 2.1466621
## p:Obs1 -1.0703414 0.4612909 -1.9744716 -0.1662111
## p:Obs2 -0.3437551 0.4362502 -1.1988055 0.5112954
## p:BrowseUnbrowsed 0.0052088 0.5293598 -1.0323365 1.0427541
## Psi:(Intercept) 1.1957593 0.7216656 -0.2187053 2.6102240
## Psi:BrowseUnbrowsed -1.1740827 0.8861536 -2.9109438 0.5627783
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3541856 0.3242823 0.1725302 0.3449677 0.6256992
## Group:BrowseUnbrowsed 0.3553780 0.3254247 0.1732751 0.3461456 0.6269183
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.7677695
## Group:BrowseUnbrowsed 0.5054189##
## p.dot.Psi.browse##
## Output summary for Occupancy model
## Name : p(~1)Psi(~Browse)
##
## Npar : 3
## -2lnL: 258.2643
## AICc : 264.6172
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6222546 0.2366510 -1.0860906 -0.1584187
## Psi:(Intercept) 1.1492336 0.6557678 -0.1360713 2.4345384
## Psi:BrowseUnbrowsed -1.2252779 0.7205607 -2.6375769 0.1870211
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3492688 0.3492688 0.3492688 0.3492688 0.3492688
## Group:BrowseUnbrowsed 0.3492688 0.3492688 0.3492688 0.3492688 0.3492688
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.7593709
## Group:BrowseUnbrowsed 0.4809981##
## p.obs.Psi.browse##
## Output summary for Occupancy model
## Name : p(~Obs1 + Obs2)Psi(~Browse)
##
## Npar : 5
## -2lnL: 252.0421
## AICc : 262.9512
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.2158421 0.3270212 -0.8568038 0.4251195
## p:Obs1 -1.0294215 0.4324360 -1.8769960 -0.1818469
## p:Obs2 -0.2796125 0.4033747 -1.0702269 0.5110020
## Psi:(Intercept) 1.1178711 0.6269315 -0.1109148 2.3466569
## Psi:BrowseUnbrowsed -1.1792046 0.7017341 -2.5546035 0.1961943
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3629248 0.3689858 0.3590562 0.3701787 0.381149
## Group:BrowseUnbrowsed 0.3629248 0.3689858 0.3590562 0.3701787 0.381149
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.7535936
## Group:BrowseUnbrowsed 0.4846714##
## p.obs.browse.Psi.browse##
## Output summary for Occupancy model
## Name : p(~Obs1 + Obs2 + Browse)Psi(~Browse)
##
## Npar : 6
## -2lnL: 252.0156
## AICc : 265.308
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.2435662 0.3699471 -0.9686625 0.4815301
## p:Obs1 -1.0359711 0.4345995 -1.8877862 -0.1841560
## p:Obs2 -0.2822415 0.4036974 -1.0734883 0.5090054
## p:BrowseUnbrowsed 0.0801247 0.4916892 -0.8835862 1.0438356
## Psi:(Intercept) 1.1560851 0.6937771 -0.2037180 2.5158883
## Psi:BrowseUnbrowsed -1.2478563 0.8329962 -2.8805288 0.3848162
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3560008 0.3620285 0.3521193 0.3632330 0.3741922
## Group:BrowseUnbrowsed 0.3745749 0.3807311 0.3706075 0.3819606 0.3931360
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.7606206
## Group:BrowseUnbrowsed 0.4770733##
## p.browse.Psi.dot##
## Output summary for Occupancy model
## Name : p(~Browse)Psi(~1)
##
## Npar : 3
## -2lnL: 260.4498
## AICc : 266.8027
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.4905233 0.2676562 -1.0151295 0.0340828
## p:BrowseUnbrowsed -0.4891207 0.4133656 -1.2993172 0.3210758
## Psi:(Intercept) 0.6363421 0.4480518 -0.2418394 1.5145237
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3797703 0.3797703 0.3797703 0.3797703 0.3797703
## Group:BrowseUnbrowsed 0.2729624 0.2729624 0.2729624 0.2729624 0.2729624
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.6539261
## Group:BrowseUnbrowsed 0.6539261##
## p.day.Psi.dot##
## Output summary for Occupancy model
## Name : p(~Day)Psi(~1)
##
## Npar : 6
## -2lnL: 248.7948
## AICc : 262.0871
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6211679 0.4312580 -1.4664335 0.2240977
## p:Day2 -0.1466836 0.5280647 -1.1816906 0.8883233
## p:Day3 -0.9798368 0.5873288 -2.1310013 0.1713276
## p:Day4 -0.1642081 0.5426922 -1.2278849 0.8994686
## p:Day5 0.9972785 0.5470389 -0.0749178 2.0694748
## Psi:(Intercept) 0.5311520 0.3961420 -0.2452864 1.3075904
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3495159 0.316944 0.1678412 0.3131624 0.5929347
## Group:BrowseUnbrowsed 0.3495159 0.316944 0.1678412 0.3131624 0.5929347
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.6297518
## Group:BrowseUnbrowsed 0.6297518##
## p.day.browse.Psi.dot##
## Output summary for Occupancy model
## Name : p(~Day + Browse)Psi(~1)
##
## Npar : 7
## -2lnL: 247.1824
## AICc : 262.9324
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.5072781 0.4370940 -1.3639824 0.3494261
## p:Day2 -0.1403519 0.5220271 -1.1635251 0.8828212
## p:Day3 -0.9099775 0.5840456 -2.0547069 0.2347520
## p:Day4 -0.1068076 0.5373022 -1.1599200 0.9463048
## p:Day5 1.0599832 0.5383203 0.0048753 2.1150911
## p:BrowseUnbrowsed -0.5489152 0.4206973 -1.3734820 0.2756516
## Psi:(Intercept) 0.7125470 0.4723259 -0.2132118 1.6383058
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3758318 0.3435238 0.1950922 0.3511278 0.6347630
## Group:BrowseUnbrowsed 0.2580376 0.2320904 0.1228008 0.2381224 0.5009475
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.6709637
## Group:BrowseUnbrowsed 0.6709637##
## p.day.obs.Psi.dot##
## Output summary for Occupancy model
## Name : p(~Day + Obs1 + Obs2)Psi(~1)
##
## Npar : 8
## -2lnL: 242.5466
## AICc : 260.8323
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.2384162 0.4726739 -1.1648570 0.6880246
## p:Day2 -0.1376623 0.5473984 -1.2105633 0.9352387
## p:Day3 -0.9640981 0.5967627 -2.1337530 0.2055568
## p:Day4 -0.0470487 0.5606895 -1.1460001 1.0519027
## p:Day5 1.0484584 0.5618024 -0.0526744 2.1495912
## p:Obs1 -1.1043723 0.4588470 -2.0037124 -0.2050322
## p:Obs2 -0.3562762 0.4349636 -1.2088049 0.4962525
## Psi:(Intercept) 0.5804494 0.4093417 -0.2218603 1.3827591
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3498032 0.3239235 0.1668317 0.3462006 0.6249679
## Group:BrowseUnbrowsed 0.3498032 0.3239235 0.1668317 0.3462006 0.6249679
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.6411708
## Group:BrowseUnbrowsed 0.6411708##
## p.day.obs.browse.Psi.dot##
## Output summary for Occupancy model
## Name : p(~Day + Obs1 + Obs2 + Browse)Psi(~1)
##
## Npar : 9
## -2lnL: 241.4147
## AICc : 262.318
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1481327 0.4801259 -1.0891794 0.7929140
## p:Day2 -0.1651114 0.5413288 -1.2261158 0.8958931
## p:Day3 -0.9041529 0.5926704 -2.0657869 0.2574810
## p:Day4 -0.0471992 0.5539722 -1.1329848 1.0385864
## p:Day5 1.0834174 0.5521178 0.0012665 2.1655682
## p:Obs1 -1.0500604 0.4556463 -1.9431271 -0.1569937
## p:Obs2 -0.3159094 0.4276196 -1.1540439 0.5222250
## p:BrowseUnbrowsed -0.4698723 0.4303053 -1.3132708 0.3735261
## Psi:(Intercept) 0.7375144 0.4801905 -0.2036591 1.6786878
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3758849 0.3432312 0.1928306 0.3719965 0.6581578
## Group:BrowseUnbrowsed 0.2735025 0.2462340 0.1299282 0.2702147 0.5461738
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.6764521
## Group:BrowseUnbrowsed 0.6764521##
## p.dot.Psi.dot##
## Output summary for Occupancy model
## Name : p(~1)Psi(~1)
##
## Npar : 2
## -2lnL: 261.7871
## AICc : 265.9611
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6218311 0.2363737 -1.0851235 -0.1585387
## Psi:(Intercept) 0.4751239 0.3745432 -0.2589808 1.2092287
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3493651 0.3493651 0.3493651 0.3493651 0.3493651
## Group:BrowseUnbrowsed 0.3493651 0.3493651 0.3493651 0.3493651 0.3493651
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.6165958
## Group:BrowseUnbrowsed 0.6165958##
## p.obs.Psi.dot##
## Output summary for Occupancy model
## Name : p(~Obs1 + Obs2)Psi(~1)
##
## Npar : 4
## -2lnL: 255.3578
## AICc : 263.9548
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.2216020 0.3270493 -0.8626185 0.4194146
## p:Obs1 -1.0413814 0.4314191 -1.8869627 -0.1958000
## p:Obs2 -0.2729065 0.4037457 -1.0642481 0.5184351
## Psi:(Intercept) 0.4862374 0.3740421 -0.2468851 1.2193600
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.361047 0.3674223 0.3573965 0.3684661 0.3795116
## Group:BrowseUnbrowsed 0.361047 0.3674223 0.3573965 0.3684661 0.3795116
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.6192197
## Group:BrowseUnbrowsed 0.6192197##
## p.obs.browse.Psi.dot##
## Output summary for Occupancy model
## Name : p(~Obs1 + Obs2 + Browse)Psi(~1)
##
## Npar : 5
## -2lnL: 254.5045
## AICc : 265.4136
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1342856 0.3424409 -0.8054698 0.5368985
## p:Obs1 -0.9958830 0.4292169 -1.8371482 -0.1546178
## p:Obs2 -0.2570374 0.3990736 -1.0392217 0.5251469
## p:BrowseUnbrowsed -0.4044099 0.4288018 -1.2448614 0.4360416
## Psi:(Intercept) 0.6197218 0.4419715 -0.2465424 1.4859860
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3851866 0.3915447 0.3817047 0.3925023 0.4032983
## Group:BrowseUnbrowsed 0.2948381 0.3004336 0.2917853 0.3012788 0.3108495
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.6501553
## Group:BrowseUnbrowsed 0.6501553Once the models are run, we can examine the model-selection table. We can also examine model-specific output.
options(width = 160)
weta.models## model npar AICc DeltaAICc weight Deviance
## 7 p(~Day + Obs1 + Obs2)Psi(~Browse) 9 260.5027 0.0000000 0.196120806 239.5994
## 8 p(~Day + Obs1 + Obs2)Psi(~1) 8 260.8323 0.3296085 0.166322153 242.5466
## 9 p(~Day)Psi(~Browse) 7 261.1868 0.6841742 0.139301929 -162.9434
## 10 p(~Day)Psi(~1) 6 262.0871 1.5844719 0.088809612 -159.5854
## 6 p(~Day + Obs1 + Obs2 + Browse)Psi(~1) 9 262.3180 1.8152900 0.079129506 241.4147
## 4 p(~Day + Browse)Psi(~1) 7 262.9324 2.4297042 0.058199614 -161.1979
## 15 p(~Obs1 + Obs2)Psi(~Browse) 5 262.9512 2.4485251 0.057654499 252.0421
## 5 p(~Day + Obs1 + Obs2 + Browse)Psi(~Browse) 10 263.2059 2.7032416 0.050760071 239.5993
## 3 p(~Day + Browse)Psi(~Browse) 8 263.6937 3.1910585 0.039773529 -162.9723
## 16 p(~Obs1 + Obs2)Psi(~1) 4 263.9548 3.4521591 0.034905739 255.3578
## 11 p(~1)Psi(~Browse) 3 264.6172 4.1145554 0.025064516 -150.1160
## 13 p(~Obs1 + Obs2 + Browse)Psi(~Browse) 6 265.3080 4.8052919 0.017744665 252.0156
## 14 p(~Obs1 + Obs2 + Browse)Psi(~1) 5 265.4136 4.9109451 0.016831603 254.5045
## 12 p(~1)Psi(~1) 2 265.9611 5.4583972 0.012801114 -146.5931
## 2 p(~Browse)Psi(~1) 3 266.8027 6.3000354 0.008404045 -147.9305
## 1 p(~Browse)Psi(~Browse) 4 266.8576 6.3549091 0.008176599 -150.1197Here, we use AIC so that we can compare direclty with results presented in the Occupancy book. You can see that things match (as they should) for most models but that something appears amiss in the book with the Psi(.)p(Obs+Browse) model as it should have 5 parameters but is only listed as having 4 in the book (it is assigned k=5 here).
options(width = 160)
# Store and view a model table that uses AIC rather than AICc
AIC.table = weta.models
AIC.table$model.table = model.table(weta.models, use.AIC = TRUE)
AIC.table## model npar AIC DeltaAIC weight Deviance
## 7 p(~Day + Obs1 + Obs2)Psi(~Browse) 9 257.5994 0.00000 0.275827626 239.5994
## 8 p(~Day + Obs1 + Obs2)Psi(~1) 8 258.5466 0.94712 0.171780264 242.5466
## 6 p(~Day + Obs1 + Obs2 + Browse)Psi(~1) 9 259.4147 1.81529 0.111289078 241.4147
## 9 p(~Day)Psi(~Browse) 7 259.4368 1.83740 0.110065553 -162.9434
## 5 p(~Day + Obs1 + Obs2 + Browse)Psi(~Browse) 10 259.5993 1.99991 0.101475879 239.5993
## 10 p(~Day)Psi(~1) 6 260.7948 3.19539 0.055817148 -159.5854
## 4 p(~Day + Browse)Psi(~1) 7 261.1824 3.58293 0.045984810 -161.1979
## 3 p(~Day + Browse)Psi(~Browse) 8 261.4080 3.80857 0.041078757 -162.9723
## 15 p(~Obs1 + Obs2)Psi(~Browse) 5 262.0421 4.44266 0.029917576 252.0421
## 16 p(~Obs1 + Obs2)Psi(~1) 4 263.3578 5.75837 0.015496143 255.3578
## 13 p(~Obs1 + Obs2 + Browse)Psi(~Browse) 6 264.0156 6.41621 0.011152583 252.0156
## 11 p(~1)Psi(~Browse) 3 264.2643 6.66484 0.009848864 -150.1160
## 14 p(~Obs1 + Obs2 + Browse)Psi(~1) 5 264.5045 6.90508 0.008734111 254.5045
## 12 p(~1)Psi(~1) 2 265.7871 8.18771 0.004599378 -146.5931
## 1 p(~Browse)Psi(~Browse) 4 266.2606 8.66112 0.003629940 -150.1197
## 2 p(~Browse)Psi(~1) 3 266.4498 8.85032 0.003302290 -147.9305# View the model names
names(weta.models)## [1] "p.browse.Psi.browse" "p.browse.Psi.dot" "p.day.browse.Psi.browse" "p.day.browse.Psi.dot" "p.day.obs.browse.Psi.browse"
## [6] "p.day.obs.browse.Psi.dot" "p.day.obs.Psi.browse" "p.day.obs.Psi.dot" "p.day.Psi.browse" "p.day.Psi.dot"
## [11] "p.dot.Psi.browse" "p.dot.Psi.dot" "p.obs.browse.Psi.browse" "p.obs.browse.Psi.dot" "p.obs.Psi.browse"
## [16] "p.obs.Psi.dot" "model.table"# examine the output from top-ranked model (#7) and
# store top model in object with short name
top <- weta.models$p.day.obs.Psi.browse
# look at summary of top model's output
summary(top, showall = FALSE)## Output summary for Occupancy model
## Name : p(~Day + Obs1 + Obs2)Psi(~Browse)
##
## Npar : 9
## -2lnL: 239.5994
## AICc : 260.5027
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.2297010 0.4732890 -1.1573474 0.6979454
## p:Day2 -0.1548925 0.5504441 -1.2337630 0.9239780
## p:Day3 -0.9431918 0.5997284 -2.1186595 0.2322759
## p:Day4 -0.0706551 0.5629219 -1.1739820 1.0326718
## p:Day5 1.0354542 0.5639409 -0.0698700 2.1407784
## p:Obs1 -1.0700177 0.4601063 -1.9718261 -0.1682093
## p:Obs2 -0.3435561 0.4357548 -1.1976355 0.5105234
## Psi:(Intercept) 1.1933135 0.6757956 -0.1312459 2.5178729
## Psi:BrowseUnbrowsed -1.1693203 0.7412937 -2.6222560 0.2836154
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowseBrowsed 0.3545613 0.3246093 0.1728514 0.3453743 0.626219
## Group:BrowseUnbrowsed 0.3545613 0.3246093 0.1728514 0.3453743 0.626219
##
##
## Real Parameter Psi
## 1
## Group:BrowseBrowsed 0.7673332
## Group:BrowseUnbrowsed 0.5059980Here, model-averaging is done for Psi in each habitat setting.
Mod.Avg.Psi <- model.average(weta.models, "Psi")
Mod.Avg.Psi## par.index estimate se fixed note group age time Age Time Browse
## Psi gBrowsed a0 t1 11 0.7102669 0.1268297 Browsed 0 1 0 0 Browsed
## Psi gUnbrowsed a0 t1 12 0.5670989 0.1314801 Unbrowsed 0 1 0 0 UnbrowsedThere is much more that can be done in RMark. And, in the help for the RMark package, the weta example has been worked out in detail. You can see the details by typing “?weta” after loading the RMark package. If you read through the example, you’ll see that much of what’s provided on pages 116-122 of the Occupancy book is worked out in the RMark ‘weta’ example.