Lab 11 - Analyses in RMark
WILD 502 - Jay Rotella
Single-Species Multi-Season Occupancy Models
Here, we’ll run the single-species, multi-season occupancy models for data on northern spotted owls that you analyzed directly in MARK in lab 11. Here, we’re dealing with 5 types of parameters: ‘Psi’, ‘Epsilon’, ‘Gamma’ & ‘p’. There are different parameterizations available. In 5 of the 6 models, we’ll directly model Psi in year 1 as well as Epsilon and Gamma for each year, and p for each sampling occasion. For those models, Psi in years 2-5 will be derived from Psi in year 1 and the relevant Epsilon and Gamma estimates as reviewed in class. In the 6th model, we’ll model Psi in each year directly along with Gamma and p and we’ll derive Epsilon for each year from the estimates of annual estimates of Psi and Gamma.
Bring in the data & process it
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
nso <- convert.inp("http://www.montana.edu/rotella/documents/502/nso.inp")
# Process data and set up the time intervals: 0 for intervals that
# occur between secondary occasions in the same year and 1 for
# intervals that occur between years (8 surveys )
# Here, we use the 'Psi(yr1), Epsilon(t), Gamma(t)' parameterization
nso.pr.EG <- process.data(nso,model = "RDOccupEG",
time.intervals = c(0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0))
nso.ddl = make.design.data(nso.pr.EG)
# add design data for models that set Epsilon to -1 and Gamma to +1
# and that hold Epsilon and Gamma constant across years
nso.ddl$Epsilon$eps = -1
nso.ddl$Gamma$eps = 1
# Make a different version of processed data for the Psi & Gamma
# parameterization that is used in model 6, i.e., 'Psi(t),Gamma(t)'
nso.pr.PG = process.data(nso, model = "RDOccupPG",
time.intervals = c(0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0))
head(nso)## ch freq
## 1 0111....010.....11......011.....0101.... 1
## 2 00......000.....000.....000110..0011.... 1
## 3 1000111.0100....0001....000000..00000... 1
## 4 111.....10......011.....00000100000011.. 1
## 5 000000..0000000.000000..000000..0111.... 1
## 6 1.......000111..0111....01111...011111.. 1nso.ddl## $Psi
## par.index model.index group age time Age Time
## 1 1 1 1 0 1 0 0
##
## $Epsilon
## par.index model.index group age time Age Time eps
## 1 1 2 1 0 1 0 0 -1
## 2 2 3 1 1 2 1 1 -1
## 3 3 4 1 2 3 2 2 -1
## 4 4 5 1 3 4 3 3 -1
##
## $Gamma
## par.index model.index group age time Age Time eps
## 1 1 6 1 0 1 0 0 1
## 2 2 7 1 1 2 1 1 1
## 3 3 8 1 2 3 2 2 1
## 4 4 9 1 3 4 3 3 1
##
## $p
## par.index model.index group time session Time
## 1 1 10 1 1 1 0
## 2 2 11 1 2 1 1
## 3 3 12 1 3 1 2
## 4 4 13 1 4 1 3
## 5 5 14 1 5 1 4
## 6 6 15 1 6 1 5
## 7 7 16 1 7 1 6
## 8 8 17 1 8 1 7
## 9 9 18 1 1 2 0
## 10 10 19 1 2 2 1
## 11 11 20 1 3 2 2
## 12 12 21 1 4 2 3
## 13 13 22 1 5 2 4
## 14 14 23 1 6 2 5
## 15 15 24 1 7 2 6
## 16 16 25 1 8 2 7
## 17 17 26 1 1 3 0
## 18 18 27 1 2 3 1
## 19 19 28 1 3 3 2
## 20 20 29 1 4 3 3
## 21 21 30 1 5 3 4
## 22 22 31 1 6 3 5
## 23 23 32 1 7 3 6
## 24 24 33 1 8 3 7
## 25 25 34 1 1 4 0
## 26 26 35 1 2 4 1
## 27 27 36 1 3 4 2
## 28 28 37 1 4 4 3
## 29 29 38 1 5 4 4
## 30 30 39 1 6 4 5
## 31 31 40 1 7 4 6
## 32 32 41 1 8 4 7
## 33 33 42 1 1 5 0
## 34 34 43 1 2 5 1
## 35 35 44 1 3 5 2
## 36 36 45 1 4 5 3
## 37 37 46 1 5 5 4
## 38 38 47 1 6 5 5
## 39 39 48 1 7 5 6
## 40 40 49 1 8 5 7
##
## $pimtypes
## $pimtypes$Psi
## $pimtypes$Psi$pim.type
## [1] "all"
##
##
## $pimtypes$Epsilon
## $pimtypes$Epsilon$pim.type
## [1] "all"
##
##
## $pimtypes$Gamma
## $pimtypes$Gamma$pim.type
## [1] "all"
##
##
## $pimtypes$p
## $pimtypes$p$pim.type
## [1] "all"Build function to create models
First, the function puts the observer and time covariates into proper form for use in modeling and processes the data. Second, the function declares 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.nso.models <- function() {
# Define formulas for constructing models
# Note that in the design data for 'p', 'session' is a
# factor based on primary occasions
p.year = list(formula = ~ session)
Epsilon.dot = list(formula = ~ 1)
Epsilon.zero = list(formula = ~ 1,
fixed = list(time = 1:8, value = 0))
Epsilon.random.shared = list(formula = ~ -1 + eps,
share = TRUE)
Epsilon.random.yr = list(formula = ~ time:eps,
share = TRUE,
remove.intercept = TRUE)
Epsilon.year = list(formula = ~ time)
Gamma.dot = list(formula = ~ 1)
Gamma.zero = list(formula = ~ 1,
fixed = list(time = 1:8, value = 0))
Gamma.year = list(formula = ~ time)
Psi.dot = list(formula = ~ 1)
# Define the models
# Model 1 - no change in occupancy (epsilon & gamma=0)
mod.1 = mark(nso.pr.EG, nso.ddl,
model.parameters = list(
p = p.year,
Psi = Psi.dot,
Epsilon = Epsilon.zero,
Gamma = Gamma.zero))
# Model 2 - random changes in occupancy (epsilon=1-gamma)
mod.2 = mark(nso.pr.EG, nso.ddl,
model.parameters = list(
p = p.year,
Psi = Psi.dot,
Epsilon = Epsilon.random.shared))
# Model 3 - random changes in occupancy (epsilon=1-gamma)
mod.3 = mark(nso.pr.EG, nso.ddl,
model.parameters = list(
p = p.year,
Psi = Psi.dot,
Epsilon = Epsilon.random.yr))
# Model 4 - annual variation in epsilon & gamma
mod.4 = mark(nso.pr.EG, nso.ddl,
model.parameters = list(
p = p.year,
Psi = Psi.dot,
Epsilon = Epsilon.dot,
Gamma = Gamma.dot))
# Model 5 - annual variation in epsilon & gamma
mod.5=mark(nso.pr.EG, nso.ddl,
model.parameters = list(
p = p.year,
Psi = Psi.dot,
Epsilon = Epsilon.year,
Gamma = Gamma.year))
# Model 6 - Work with different parameterization for which model psi by year.
# For this model, psi & gamma are held constant across years
mod.6 = mark(nso.pr.PG,
model.parameters = list(
p = p.year,
Psi = Psi.dot,
Gamma = Gamma.dot))
# Collect the models to run when the function is called
return(collect.models() )
}Run the models and examine the output
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. One can also examine model-specific output in other ways as shown below.
nso.models <- fit.nso.models()##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~session)
##
## Npar : 6
## -2lnL: 1542.564
## AICc : 1554.878
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 1.6315383 0.3645234 0.9170724 2.3460041
## p:(Intercept) -0.2824360 0.1358812 -0.5487632 -0.0161088
## p:session2 -0.2959648 0.1930479 -0.6743388 0.0824091
## p:session3 -0.7590465 0.2025130 -1.1559720 -0.3621210
## p:session4 -0.7627181 0.1989523 -1.1526646 -0.3727716
## p:session5 -0.3100959 0.1870379 -0.6766901 0.0564983
##
##
## Real Parameter Psi
##
## 1
## 0.8363803
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## NA NA NA NA
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## NA NA NA NA
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.4298567 0.4298567 0.4298567 0.4298567 0.4298567 0.4298567 0.4298567
## 8
## 0.4298567
##
## Session:2
## 1 2 3 4 5 6 7
## 0.3593006 0.3593006 0.3593006 0.3593006 0.3593006 0.3593006 0.3593006
## 8
## 0.3593006
##
## Session:3
## 1 2 3 4 5 6 7 8
## 0.260864 0.260864 0.260864 0.260864 0.260864 0.260864 0.260864 0.260864
##
## Session:4
## 1 2 3 4 5 6 7
## 0.2601567 0.2601567 0.2601567 0.2601567 0.2601567 0.2601567 0.2601567
## 8
## 0.2601567
##
## Session:5
## 1 2 3 4 5 6 7
## 0.3560541 0.3560541 0.3560541 0.3560541 0.3560541 0.3560541 0.3560541
## 8
## 0.3560541
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~-1 + eps)Gamma()p(~session)
##
## Npar : 7
## -2lnL: 1429.532
## AICc : 1443.951
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5103584 0.2902428 -0.0585176 1.0792343
## Epsilon:eps 0.3440392 0.1455373 0.0587860 0.6292924
## p:(Intercept) 0.3528209 0.1661479 0.0271709 0.6784709
## p:session2 -0.2593688 0.2360301 -0.7219877 0.2032502
## p:session3 -0.7435364 0.2533379 -1.2400786 -0.2469941
## p:session4 -0.8259734 0.2447913 -1.3057643 -0.3461824
## p:session5 -0.2099905 0.2295678 -0.6599434 0.2399625
##
##
## Real Parameter Psi
##
## 1
## 0.6248905
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.4148286 0.4148286 0.4148286 0.4148286
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.5851714 0.5851714 0.5851714 0.5851714
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5873015 0.5873015 0.5873015 0.5873015 0.5873015 0.5873015 0.5873015
## 8
## 0.5873015
##
## Session:2
## 1 2 3 4 5 6 7 8
## 0.523346 0.523346 0.523346 0.523346 0.523346 0.523346 0.523346 0.523346
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4035451 0.4035451 0.4035451 0.4035451 0.4035451 0.4035451 0.4035451
## 8
## 0.4035451
##
## Session:4
## 1 2 3 4 5 6 7
## 0.3838704 0.3838704 0.3838704 0.3838704 0.3838704 0.3838704 0.3838704
## 8
## 0.3838704
##
## Session:5
## 1 2 3 4 5 6 7 8
## 0.535647 0.535647 0.535647 0.535647 0.535647 0.535647 0.535647 0.535647
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~time:eps)Gamma()p(~session)
##
## Npar : 10
## -2lnL: 1429.32
## AICc : 1450.153
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5103576 0.2902428 -0.0585184 1.0792335
## Epsilon:time1:eps 0.4099552 0.2866077 -0.1517960 0.9717064
## Epsilon:time2:eps 0.2664923 0.2987904 -0.3191368 0.8521215
## Epsilon:time3:eps 0.4132185 0.3047061 -0.1840056 1.0104425
## Epsilon:time4:eps 0.2888087 0.2773775 -0.2548512 0.8324686
## p:(Intercept) 0.3528209 0.1661480 0.0271708 0.6784709
## p:session2 -0.2626339 0.2366930 -0.7265523 0.2012844
## p:session3 -0.7332369 0.2542976 -1.2316601 -0.2348136
## p:session4 -0.8346149 0.2481969 -1.3210807 -0.3481490
## p:session5 -0.2082451 0.2294520 -0.6579711 0.2414809
##
##
## Real Parameter Psi
##
## 1
## 0.6248903
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.3989229 0.4337684 0.3981406 0.4282955
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.6010771 0.5662316 0.6018594 0.5717045
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5873015 0.5873015 0.5873015 0.5873015 0.5873015 0.5873015 0.5873015
## 8
## 0.5873015
##
## Session:2
## 1 2 3 4 5 6 7
## 0.5225315 0.5225315 0.5225315 0.5225315 0.5225315 0.5225315 0.5225315
## 8
## 0.5225315
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4060266 0.4060266 0.4060266 0.4060266 0.4060266 0.4060266 0.4060266
## 8
## 0.4060266
##
## Session:4
## 1 2 3 4 5 6 7
## 0.3818286 0.3818286 0.3818286 0.3818286 0.3818286 0.3818286 0.3818286
## 8
## 0.3818286
##
## Session:5
## 1 2 3 4 5 6 7
## 0.5360811 0.5360811 0.5360811 0.5360811 0.5360811 0.5360811 0.5360811
## 8
## 0.5360811
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~session)
##
## Npar : 8
## -2lnL: 1337.523
## AICc : 1354.064
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5096833 0.2853277 -0.0495591 1.0689256
## Epsilon:(Intercept) -1.7963669 0.2687696 -2.3231553 -1.2695785
## Gamma:(Intercept) -1.5242751 0.2929070 -2.0983729 -0.9501773
## p:(Intercept) 0.3663497 0.1630641 0.0467440 0.6859554
## p:session2 -0.2764254 0.2294433 -0.7261343 0.1732835
## p:session3 -0.7406938 0.2449777 -1.2208501 -0.2605374
## p:session4 -0.8353613 0.2386106 -1.3030381 -0.3676845
## p:session5 -0.2197445 0.2267269 -0.6641291 0.2246402
##
##
## Real Parameter Psi
##
## 1
## 0.6247322
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.1422939 0.1422939 0.1422939 0.1422939
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.1788328 0.1788328 0.1788328 0.1788328
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5905766 0.5905766 0.5905766 0.5905766 0.5905766 0.5905766 0.5905766
## 8
## 0.5905766
##
## Session:2
## 1 2 3 4 5 6 7
## 0.5224659 0.5224659 0.5224659 0.5224659 0.5224659 0.5224659 0.5224659
## 8
## 0.5224659
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4074918 0.4074918 0.4074918 0.4074918 0.4074918 0.4074918 0.4074918
## 8
## 0.4074918
##
## Session:4
## 1 2 3 4 5 6 7
## 0.3848502 0.3848502 0.3848502 0.3848502 0.3848502 0.3848502 0.3848502
## 8
## 0.3848502
##
## Session:5
## 1 2 3 4 5 6 7
## 0.5365858 0.5365858 0.5365858 0.5365858 0.5365858 0.5365858 0.5365858
## 8
## 0.5365858
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~time)Gamma(~time)p(~session)
##
## Npar : 14
## -2lnL: 1327.639
## AICc : 1357.255
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5301133 0.2856344 -0.0297301 1.0899567
## Epsilon:(Intercept) -2.3314280 0.6232623 -3.5530222 -1.1098338
## Epsilon:time2 0.4654249 0.8563261 -1.2129742 2.1438240
## Epsilon:time3 1.1739848 0.7855196 -0.3656336 2.7136033
## Epsilon:time4 0.3276642 0.8600679 -1.3580690 2.0133974
## Gamma:(Intercept) -2.1158499 0.7697179 -3.6244970 -0.6072028
## Gamma:time2 -0.4832600 1.2829017 -2.9977475 2.0312275
## Gamma:time3 1.6524506 0.8902902 -0.0925183 3.3974195
## Gamma:time4 0.0876465 1.1419653 -2.1506055 2.3258985
## p:(Intercept) 0.3612745 0.1634012 0.0410080 0.6815409
## p:session2 -0.2839538 0.2291279 -0.7330446 0.1651369
## p:session3 -0.7063157 0.2460576 -1.1885887 -0.2240427
## p:session4 -0.8327926 0.2453485 -1.3136756 -0.3519096
## p:session5 -0.2168003 0.2274620 -0.6626258 0.2290252
##
##
## Real Parameter Psi
##
## 1
## 0.6295095
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.0885533 0.1340049 0.2391322 0.1188083
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.1075658 0.0691957 0.3861797 0.1162734
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5893489 0.5893489 0.5893489 0.5893489 0.5893489 0.5893489 0.5893489
## 8
## 0.5893489
##
## Session:2
## 1 2 3 4 5 6 7
## 0.5193205 0.5193205 0.5193205 0.5193205 0.5193205 0.5193205 0.5193205
## 8
## 0.5193205
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4145854 0.4145854 0.4145854 0.4145854 0.4145854 0.4145854 0.4145854
## 8
## 0.4145854
##
## Session:4
## 1 2 3 4 5 6 7 8
## 0.384257 0.384257 0.384257 0.384257 0.384257 0.384257 0.384257 0.384257
##
## Session:5
## 1 2 3 4 5 6 7
## 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558
## 8
## 0.5360558
##
## Output summary for RDOccupPG model
## Name : Psi(~1)Gamma(~1)p(~session)
##
## Npar : 7
## -2lnL: 1337.951
## AICc : 1352.37
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.3848143 0.2087567 -0.0243489 0.7939775
## Gamma:(Intercept) -1.4283492 0.2436051 -1.9058151 -0.9508832
## p:(Intercept) 0.3692626 0.1627455 0.0502814 0.6882439
## p:session2 -0.2790083 0.2292946 -0.7284257 0.1704091
## p:session3 -0.7460380 0.2448595 -1.2259625 -0.2661134
## p:session4 -0.8407514 0.2383948 -1.3080053 -0.3734976
## p:session5 -0.2253229 0.2267040 -0.6696628 0.2190170
##
##
## Real Parameter Psi
##
## 1 2 3 4 5
## 0.5950337 0.5950337 0.5950337 0.5950337 0.5950337
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.193356 0.193356 0.193356 0.193356
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5912808 0.5912808 0.5912808 0.5912808 0.5912808 0.5912808 0.5912808
## 8
## 0.5912808
##
## Session:2
## 1 2 3 4 5 6 7
## 0.5225483 0.5225483 0.5225483 0.5225483 0.5225483 0.5225483 0.5225483
## 8
## 0.5225483
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4069049 0.4069049 0.4069049 0.4069049 0.4069049 0.4069049 0.4069049
## 8
## 0.4069049
##
## Session:4
## 1 2 3 4 5 6 7
## 0.3842639 0.3842639 0.3842639 0.3842639 0.3842639 0.3842639 0.3842639
## 8
## 0.3842639
##
## Session:5
## 1 2 3 4 5 6 7
## 0.5359229 0.5359229 0.5359229 0.5359229 0.5359229 0.5359229 0.5359229
## 8
## 0.5359229## Warning in model.table(x, type, pf = 2, adjust = adjust): Model list contains models of differing typesExamine model-selection table
Once the models are run, we can examine the model-selection table. We can also examine model-specific output as shown below.
options(width = 160)
nso.models## model npar AICc DeltaAICc weight Deviance
## 6 Psi(~1)Gamma(~1)p(~session) 7 1352.370 0.000000 0.65973001 897.1443
## 4 Psi(~1)Epsilon(~1)Gamma(~1)p(~session) 8 1354.064 1.693578 0.28288538 896.7160
## 5 Psi(~1)Epsilon(~time)Gamma(~time)p(~session) 14 1357.254 4.884109 0.05738461 886.8325
## 2 Psi(~1)Epsilon(~-1 + eps)Gamma()p(~session) 7 1443.951 91.580800 0.00000000 988.7251
## 3 Psi(~1)Epsilon(~time:eps)Gamma()p(~session) 10 1450.154 97.783158 0.00000000 988.5135
## 1 Psi(~1)Epsilon(~1)Gamma(~1)p(~session) 6 1554.878 202.507257 0.00000000 1101.7575Look at output from top model
# View the model names; here they are just mod.1, mod.2, ..., mod.6,
# so we don't really need to rename the model like we did in some labs
names(nso.models)## [1] "mod.1" "mod.2" "mod.3" "mod.4" "mod.5" "mod.6" "model.table"# examine the output from top-ranked model (#6) and
nso.models$mod.6$results$real## estimate se lcl ucl fixed note
## Psi g1 a0 t1 0.5950337 0.0503038 0.4939131 0.6886847
## Gamma g1 a0 t1 0.1933560 0.0379950 0.1294517 0.2787072
## p g1 s1 t1 0.5912808 0.0393304 0.5125677 0.6655762
## p g1 s2 t1 0.5225483 0.0404801 0.4433058 0.6006717
## p g1 s3 t1 0.4069049 0.0442077 0.3239210 0.4955650
## p g1 s4 t1 0.3842639 0.0412502 0.3072093 0.4676008
## p g1 s5 t1 0.5359229 0.0392718 0.4587071 0.6114540# for mod.6, estimates of Epsilon are derived
nso.models$mod.6$results$derived## $`epsilon Probabiliy Extinction`
## estimate lcl ucl
## 1 0.1315937 0.07797401 0.1852133
## 2 0.1315937 0.07797401 0.1852133
## 3 0.1315937 0.07797401 0.1852133
## 4 0.1315937 0.07797401 0.1852133
##
## $`lambda Rate of Change`
## estimate lcl ucl
## 1 1 1 1
## 2 1 1 1
## 3 1 1 1
## 4 1 1 1
##
## $`log odds lambda`
## estimate lcl ucl
## 1 1 1 1
## 2 1 1 1
## 3 1 1 1
## 4 1 1 1Model-Specific Estimates of Psi for each year
Here, we look at model-specific estimates of psi in year 1 and for each of the subsequent years. Remember, the latter are derived parameters in models 1-5.
# store & view model-specific estimates of Psi
(psi.mod.1 = nso.models$mod.1$results$derived$`psi Probability Occupied`)## estimate lcl ucl
## 1 0.8363803 0.7386068 0.9341537
## 2 0.8363803 0.7386068 0.9341537
## 3 0.8363803 0.7386068 0.9341537
## 4 0.8363803 0.7386068 0.9341537
## 5 0.8363803 0.7386068 0.9341537(psi.mod.2 = nso.models$mod.2$results$derived$`psi Probability Occupied`)## estimate lcl ucl
## 1 0.6248905 0.4915446 0.7582363
## 2 0.5851714 0.5159273 0.6544154
## 3 0.5851714 0.5159273 0.6544154
## 4 0.5851714 0.5159273 0.6544154
## 5 0.5851714 0.5159273 0.6544154(psi.mod.3 = nso.models$mod.3$results$derived$`psi Probability Occupied`)## estimate lcl ucl
## 1 0.6248903 0.4915444 0.7582362
## 2 0.6010771 0.4663785 0.7357758
## 3 0.5662316 0.4223932 0.7100699
## 4 0.6018594 0.4587497 0.7449690
## 5 0.5717045 0.4385847 0.7048242(psi.mod.4 = nso.models$mod.4$results$derived$`psi Probability Occupied`)## estimate lcl ucl
## 1 0.6247322 0.4936224 0.7558420
## 2 0.6029468 0.5021380 0.7037556
## 3 0.5881573 0.4871583 0.6891564
## 4 0.5781171 0.4664751 0.6897591
## 5 0.5713011 0.4486144 0.6939879(psi.mod.5 = nso.models$mod.5$results$derived$`psi Probability Occupied`)## estimate lcl ucl
## 1 0.6295095 0.4989388 0.7600803
## 2 0.6136165 0.4825076 0.7447253
## 3 0.5581250 0.4206184 0.6956315
## 4 0.5953025 0.4550518 0.7355532
## 5 0.5716312 0.4385636 0.7046987# for mod.6, Psi-hat is held constant across all years, so
# we can just repeat that row 6 times (which takes a bit of work
# to get the results to look like the results for other models)
psi.mod.6 <- do.call("rbind",
replicate(5,
nso.models$mod.6$results$real[1, -c(2, 5:6)],
simplify = FALSE))
row.names(psi.mod.6) <- 1:5
psi.mod.6## estimate lcl ucl
## 1 0.5950337 0.4939131 0.6886847
## 2 0.5950337 0.4939131 0.6886847
## 3 0.5950337 0.4939131 0.6886847
## 4 0.5950337 0.4939131 0.6886847
## 5 0.5950337 0.4939131 0.6886847Graph the estimates of psi for each year from the various models
First, we’ll store the estimates in a data.frame that records which model they came from. Then, we’ll use the ‘ggplot2’ package to plot the estimates.
psi.hats <- rbind(psi.mod.1,
psi.mod.2,
psi.mod.3,
psi.mod.4,
psi.mod.5,
psi.mod.6)
# create a variable with the model number
psi.hats$model = rep(1:6, each = 5)
# create a variable for year
psi.hats$year = rep(1997:2001, 6)
head(psi.hats)## estimate lcl ucl model year
## 1 0.8363803 0.7386068 0.9341537 1 1997
## 2 0.8363803 0.7386068 0.9341537 1 1998
## 3 0.8363803 0.7386068 0.9341537 1 1999
## 4 0.8363803 0.7386068 0.9341537 1 2000
## 5 0.8363803 0.7386068 0.9341537 1 2001
## 6 0.6248905 0.4915446 0.7582363 2 1997# make the graph
library(ggplot2)
ggplot(psi.hats, aes(x = year, y = estimate, group = factor(model))) +
geom_line(size = 1.5, aes(colour = factor(model))) +
geom_ribbon(aes(ymin = lcl, ymax = ucl), alpha = 0.2) +
scale_colour_brewer(palette = "Set1", name = "Model #") +
ylim(0, 1) +
theme(legend.position = c(0.2, 0.02),
legend.justification = c(1, 0)) +
xlab("Year") +
ylab("Estimated Occupancy Rate") 