Estimation of demographic parameters

class: center, middle, title-slide

# Estimation of demographic parameters - live demos
### Olivier Gimenez for the team

---

class: center, middle
background-size: cover

## Live demo #1: A simple CJS model

---
## Load RMark package

pre {
  font-size: 22px
}
</style>

```r
library(RMark)
```

---
## Read in data (1)

Read in csv file:

```r
dipper_csv <- read_csv("dat/dipper.csv")
```

---
## Read in data (2)

```r
dipper_csv
```

```
# A tibble: 294 × 9
   year_1981 year_1982 year_1983 year_1984 year_1985 year_1986 year_1987 sex  
       <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl> <chr>
 1         1         1         1         1         1         1         0 M    
 2         1         1         1         1         1         0         0 F    
 3         1         1         1         1         0         0         0 M    
 4         1         1         1         1         0         0         0 F    
 5         1         1         0         1         1         1         0 F    
 6         1         1         0         0         0         0         0 M    
 7         1         1         0         0         0         0         0 M    
 8         1         1         0         0         0         0         0 M    
 9         1         1         0         0         0         0         0 M    
10         1         1         0         0         0         0         0 F    
# … with 284 more rows, and 1 more variable: wing_length <dbl>
```

---
## Convert to Mark format (1)

Build a data.frame with encounter histories, sex and (scaled) wing length:

```r
dipper <- data.frame(ch = unite(dipper_csv %>%
                                  select(year_1981:year_1987), 
                                col = "ch", 
                                sep = ""), 
                     sex = dipper_csv$sex,
                     wglength = scale(dipper_csv$wing_length))
```

---
## Convert to Mark format (2)

Inspect first lines:

```r
head(dipper)
```

```
       ch sex wglength
1 1111110   M     0.76
2 1111100   F    -0.87
3 1111000   M     0.53
4 1111000   F    -1.56
5 1101110   F    -1.33
6 1100000   M     1.22
```

---
## Fit model with parameters constant over time

```r
phi.p <- mark(dipper)
```

```

Output summary for CJS model
Name : Phi(~1)p(~1)

Npar :  2
-2lnL:  667
AICc :  671

Beta
                estimate   se   lcl  ucl
Phi:(Intercept)     0.24 0.10 0.042 0.44
p:(Intercept)       2.23 0.33 1.589 2.86

Real Parameter Phi
 
     1    2    3    4    5    6
1 0.56 0.56 0.56 0.56 0.56 0.56
2      0.56 0.56 0.56 0.56 0.56
3           0.56 0.56 0.56 0.56
4                0.56 0.56 0.56
5                     0.56 0.56
6                          0.56

Real Parameter p
 
    2   3   4   5   6   7
1 0.9 0.9 0.9 0.9 0.9 0.9
2     0.9 0.9 0.9 0.9 0.9
3         0.9 0.9 0.9 0.9
4             0.9 0.9 0.9
5                 0.9 0.9
6                     0.9
```

---
## Explore results (1)

We get a list with many components:

```r
str(phi.p, max.level = 0)
```

```
List of 25
 - attr(*, "class")= chr [1:2] "mark" "CJS"
```

---
## Explore results (2)

What are the components of this list?

```r
names(phi.p)
```

```
 [1] "data"             "model"            "title"            "model.name"      
 [5] "links"            "mixtures"         "call"             "parameters"      
 [9] "time.intervals"   "number.of.groups" "group.labels"     "nocc"            
[13] "begin.time"       "covariates"       "fixed"            "design.matrix"   
[17] "pims"             "design.data"      "strata.labels"    "mlogit.list"     
[21] "profile.int"      "simplify"         "model.parameters" "results"         
[25] "output"          
```

---
## Explore results (3)

Get parameters estimates:

```r
phi.p$results$real
```

```
                estimate    se  lcl  ucl fixed note
Phi g1 c1 a0 t1     0.56 0.025 0.51 0.61           
p g1 c1 a1 t2       0.90 0.029 0.83 0.95           
```

---
class: center, middle
background-size: cover

## Variation around the CJS model

---
## Preliminaries

+ Process data (number of capture occasions, time intervals, groups):

```r
dipper.proc <- process.data(data = dipper, groups = "sex")
```

+ Create design matrix with time and group structure:

```r
dipper.ddl <- make.design.data(dipper.proc)
```

+ Define structure for survival and detection probabilities:

```r
phit <- list(formula=~time) # time 
phi <- list(formula=~1)     # constant
pt <- list(formula=~time)   # time
p <- list(formula=~1)       # constant
```

---
## Fit CJS model to dipper data

```r
phit.pt <- mark(dipper.proc,
                dipper.ddl,
                model.parameters = list(Phi = phit, p = pt))
```

---

```

Output summary for CJS model
Name : Phi(~time)p(~time)

Npar :  12  (unadjusted=11)
-2lnL:  657
AICc :  682  (unadjusted=679.58789)

Beta
                estimate   se    lcl    ucl
Phi:(Intercept)    0.935 0.77 -0.571  2.442
Phi:time2         -1.198 0.87 -2.905  0.508
Phi:time3         -1.023 0.80 -2.600  0.555
Phi:time4         -0.420 0.81 -2.006  1.166
Phi:time5         -0.536 0.80 -2.110  1.038
Phi:time6          0.248 0.00  0.248  0.248
p:(Intercept)      0.829 0.78 -0.707  2.365
p:time3            1.656 1.29 -0.875  4.187
p:time4            1.522 1.07 -0.581  3.625
p:time5            1.377 0.99 -0.561  3.314
p:time6            1.795 1.07 -0.300  3.890
p:time7           -0.015 0.00 -0.015 -0.015

Real Parameter Phi
Group:sexF 
     1    2    3    4   5    6
1 0.72 0.43 0.48 0.63 0.6 0.77
2      0.43 0.48 0.63 0.6 0.77
3           0.48 0.63 0.6 0.77
4                0.63 0.6 0.77
5                     0.6 0.77
6                         0.77

Group:sexM 
     1    2    3    4   5    6
1 0.72 0.43 0.48 0.63 0.6 0.77
2      0.43 0.48 0.63 0.6 0.77
3           0.48 0.63 0.6 0.77
4                0.63 0.6 0.77
5                     0.6 0.77
6                         0.77

Real Parameter p
Group:sexF 
    2    3    4   5    6    7
1 0.7 0.92 0.91 0.9 0.93 0.69
2     0.92 0.91 0.9 0.93 0.69
3          0.91 0.9 0.93 0.69
4               0.9 0.93 0.69
5                   0.93 0.69
6                        0.69

Group:sexM 
    2    3    4   5    6    7
1 0.7 0.92 0.91 0.9 0.93 0.69
2     0.92 0.91 0.9 0.93 0.69
3          0.91 0.9 0.93 0.69
4               0.9 0.93 0.69
5                   0.93 0.69
6                        0.69
```

---

```r
phit.pt$results$real
```

```
                estimate    se  lcl  ucl fixed note
Phi gF c1 a0 t1     0.72 0.156 0.36 0.92           
Phi gF c1 a1 t2     0.43 0.069 0.31 0.57           
Phi gF c1 a2 t3     0.48 0.060 0.36 0.59           
Phi gF c1 a3 t4     0.63 0.059 0.50 0.73           
Phi gF c1 a4 t5     0.60 0.056 0.49 0.70           
Phi gF c1 a5 t6     0.77 0.000 0.77 0.77           
p gF c1 a1 t2       0.70 0.166 0.33 0.91           
p gF c1 a2 t3       0.92 0.073 0.62 0.99           
p gF c1 a3 t4       0.91 0.058 0.71 0.98           
p gF c1 a4 t5       0.90 0.054 0.74 0.97           
p gF c1 a5 t6       0.93 0.046 0.77 0.98           
p gF c1 a6 t7       0.69 0.000 0.69 0.69           
```

---
## Fit model with time-varying survival and constant detection

```r
phit.p <- mark(dipper.proc,
                dipper.ddl,
                model.parameters = list(Phi = phit, p = p))
```

---

```

Output summary for CJS model
Name : Phi(~time)p(~1)

Npar :  7
-2lnL:  660
AICc :  674

Beta
                estimate   se   lcl  ucl
Phi:(Intercept)   0.5144 0.48 -0.42 1.45
Phi:time2        -0.6981 0.55 -1.78 0.39
Phi:time3        -0.6009 0.53 -1.64 0.44
Phi:time4        -0.0061 0.53 -1.05 1.04
Phi:time5        -0.0757 0.53 -1.11 0.96
Phi:time6        -0.1781 0.53 -1.21 0.85
p:(Intercept)     2.2204 0.33  1.58 2.87

Real Parameter Phi
Group:sexF 
     1    2    3    4    5    6
1 0.63 0.45 0.48 0.62 0.61 0.58
2      0.45 0.48 0.62 0.61 0.58
3           0.48 0.62 0.61 0.58
4                0.62 0.61 0.58
5                     0.61 0.58
6                          0.58

Group:sexM 
     1    2    3    4    5    6
1 0.63 0.45 0.48 0.62 0.61 0.58
2      0.45 0.48 0.62 0.61 0.58
3           0.48 0.62 0.61 0.58
4                0.62 0.61 0.58
5                     0.61 0.58
6                          0.58

Real Parameter p
Group:sexF 
    2   3   4   5   6   7
1 0.9 0.9 0.9 0.9 0.9 0.9
2     0.9 0.9 0.9 0.9 0.9
3         0.9 0.9 0.9 0.9
4             0.9 0.9 0.9
5                 0.9 0.9
6                     0.9

Group:sexM 
    2   3   4   5   6   7
1 0.9 0.9 0.9 0.9 0.9 0.9
2     0.9 0.9 0.9 0.9 0.9
3         0.9 0.9 0.9 0.9
4             0.9 0.9 0.9
5                 0.9 0.9
6                     0.9
```

---
## Fit model with constant survival and time-varying detection

```r
phi.pt <- mark(dipper.proc,
                dipper.ddl,
                model.parameters = list(Phi = phi, p = pt))
```

---

```

Output summary for CJS model
Name : Phi(~1)p(~time)

Npar :  7
-2lnL:  664
AICc :  679

Beta
                estimate   se     lcl  ucl
Phi:(Intercept)     0.21 0.11 -0.0066 0.43
p:(Intercept)       1.30 0.74 -0.1622 2.75
p:time3             0.80 1.16 -1.4800 3.08
p:time4             0.65 1.00 -1.3124 2.61
p:time5             1.00 0.95 -0.8554 2.85
p:time6             1.47 1.03 -0.5537 3.49
p:time7             1.99 3.06 -4.0157 8.00

Real Parameter Phi
Group:sexF 
     1    2    3    4    5    6
1 0.55 0.55 0.55 0.55 0.55 0.55
2      0.55 0.55 0.55 0.55 0.55
3           0.55 0.55 0.55 0.55
4                0.55 0.55 0.55
5                     0.55 0.55
6                          0.55

Group:sexM 
     1    2    3    4    5    6
1 0.55 0.55 0.55 0.55 0.55 0.55
2      0.55 0.55 0.55 0.55 0.55
3           0.55 0.55 0.55 0.55
4                0.55 0.55 0.55
5                     0.55 0.55
6                          0.55

Real Parameter p
Group:sexF 
     2    3    4    5    6    7
1 0.79 0.89 0.88 0.91 0.94 0.96
2      0.89 0.88 0.91 0.94 0.96
3           0.88 0.91 0.94 0.96
4                0.91 0.94 0.96
5                     0.94 0.96
6                          0.96

Group:sexM 
     2    3    4    5    6    7
1 0.79 0.89 0.88 0.91 0.94 0.96
2      0.89 0.88 0.91 0.94 0.96
3           0.88 0.91 0.94 0.96
4                0.91 0.94 0.96
5                     0.94 0.96
6                          0.96
```

---
## Model ranking with AIC

```r
collect.models()
```

```
               model npar AICc DeltaAICc weight Deviance
1       Phi(~1)p(~1)    2  671       0.0 0.8112       58
3    Phi(~time)p(~1)    7  674       3.1 0.1694       77
2    Phi(~1)p(~time)    7  679       7.9 0.0158       82
4 Phi(~time)p(~time)   12  682      10.8 0.0036       74
```

+ Best model has constant parameters

---
class: center, middle
background-size: cover

## Live demo #2: Group effect

---
## Define sex-specific survival

```r
phisex <- list(formula=~sex)
```

---
## Fit model with sex-specific survival, and constant detection

```r
phisex.p <- mark(dipper.proc,
                 dipper.ddl,
                 model.parameters = list(Phi = phisex, p = p))
```

---

```

Output summary for CJS model
Name : Phi(~sex)p(~1)

Npar :  3
-2lnL:  667
AICc :  673

Beta
                estimate   se   lcl  ucl
Phi:(Intercept)    0.204 0.14 -0.07 0.48
Phi:sexM           0.079 0.20 -0.31 0.47
p:(Intercept)      2.227 0.33  1.59 2.86

Group:sexM 
     1    2    3    4    5    6
1 0.57 0.57 0.57 0.57 0.57 0.57
2      0.57 0.57 0.57 0.57 0.57
3           0.57 0.57 0.57 0.57
4                0.57 0.57 0.57
5                     0.57 0.57
6                          0.57

Group:sexM 
    2   3   4   5   6   7
1 0.9 0.9 0.9 0.9 0.9 0.9
2     0.9 0.9 0.9 0.9 0.9
3         0.9 0.9 0.9 0.9
4             0.9 0.9 0.9
5                 0.9 0.9
6                     0.9
```

---
## Parameter estimates (1)

+ On the logit scale:

```r
phisex.p$results$beta
```

```
                estimate   se   lcl  ucl
Phi:(Intercept)    0.204 0.14 -0.07 0.48
Phi:sexM           0.079 0.20 -0.31 0.47
p:(Intercept)      2.227 0.33  1.59 2.86
```

---
## Parameter estimates (2)

+ Female survival:

```r
logitphiF <- phisex.p$results$beta[1,1]
plogis(logitphiF)
```

```
[1] 0.55
```

Male survival:

```r
logitphiM <- logitphiF + phisex.p$results$beta[2,1]
plogis(logitphiM)
```

```
[1] 0.57
```

---
## Parameter estimates (3)

All at once:

```r
phisex.p$results$real
```

```
                estimate    se  lcl  ucl fixed note
Phi gF c1 a0 t1     0.55 0.035 0.48 0.62           
Phi gM c1 a0 t1     0.57 0.035 0.50 0.64           
p gF c1 a1 t2       0.90 0.029 0.83 0.95           
```

---
## Model ranking with AIC

```r
collect.models()
```

```
               model npar AICc DeltaAICc weight Deviance
1       Phi(~1)p(~1)    2  671       0.0 0.6150       58
3     Phi(~sex)p(~1)    3  673       1.9 0.2418       84
4    Phi(~time)p(~1)    7  674       3.1 0.1285       77
2    Phi(~1)p(~time)    7  679       7.9 0.0119       82
5 Phi(~time)p(~time)   12  682      10.8 0.0027       74
```

---
## Model with additive effect of sex and time on survival and constant detection

```r
phisexpt <- list(formula=~sex+time)

phisexpt.p <- mark(dipper.proc,
                   dipper.ddl,
                   model.parameters = list(Phi = phisexpt, 
                                           p = p))
```

---

```

Output summary for CJS model
Name : Phi(~sex + time)p(~1)

Npar :  8
-2lnL:  660
AICc :  676

Beta
                estimate   se   lcl  ucl
Phi:(Intercept)   0.4829 0.49 -0.48 1.44
Phi:sexM          0.0569 0.20 -0.34 0.45
Phi:time2        -0.6916 0.55 -1.78 0.39
Phi:time3        -0.5962 0.53 -1.64 0.44
Phi:time4        -0.0023 0.53 -1.05 1.04
Phi:time5        -0.0729 0.53 -1.11 0.96
Phi:time6        -0.1754 0.53 -1.21 0.86
p:(Intercept)     2.2216 0.33  1.58 2.87

Real Parameter Phi
Group:sexF 
     1    2    3    4   5    6
1 0.62 0.45 0.47 0.62 0.6 0.58
2      0.45 0.47 0.62 0.6 0.58
3           0.47 0.62 0.6 0.58
4                0.62 0.6 0.58
5                     0.6 0.58
6                         0.58

Group:sexM 
     1    2    3    4    5    6
1 0.63 0.46 0.49 0.63 0.61 0.59
2      0.46 0.49 0.63 0.61 0.59
3           0.49 0.63 0.61 0.59
4                0.63 0.61 0.59
5                     0.61 0.59
6                          0.59

Group:sexM 
    2   3   4   5   6   7
1 0.9 0.9 0.9 0.9 0.9 0.9
2     0.9 0.9 0.9 0.9 0.9
3         0.9 0.9 0.9 0.9
4             0.9 0.9 0.9
5                 0.9 0.9
6                     0.9
```

---

```r
phisexpt.p$results$real
```

```
                estimate    se  lcl  ucl fixed note
Phi gF c1 a0 t1     0.62 0.115 0.38 0.81           
Phi gF c1 a1 t2     0.45 0.070 0.32 0.59           
Phi gF c1 a2 t3     0.47 0.063 0.35 0.59           
Phi gF c1 a3 t4     0.62 0.062 0.49 0.73           
Phi gF c1 a4 t5     0.60 0.060 0.48 0.71           
Phi gF c1 a5 t6     0.58 0.062 0.45 0.69           
Phi gM c1 a0 t1     0.63 0.113 0.40 0.82           
Phi gM c1 a1 t2     0.46 0.072 0.33 0.60           
Phi gM c1 a2 t3     0.49 0.064 0.36 0.61           
Phi gM c1 a3 t4     0.63 0.062 0.50 0.74           
Phi gM c1 a4 t5     0.61 0.059 0.49 0.72           
Phi gM c1 a5 t6     0.59 0.062 0.47 0.70           
p gF c1 a1 t2       0.90 0.029 0.83 0.95           
```

---

---

```r
# build dataset
df <- data.frame(year = c(1981:1986, 1981:1986),
                 survival = c(phisexpt.p$results$real[1:6,1],
                              phisexpt.p$results$real[7:12,1]),
                 sex = c(rep("female",6), rep("male", 6)))

# visualization
df %>%
  ggplot() + 
  aes(x = year, y = survival, color = sex) +
  geom_line() + 
  labs(color = NULL) + 
  theme_light(base_size = 18) + 
  scale_color_viridis_d(begin = 0, end = 0.7)
```

---
## Model ranking with AIC

```r
collect.models()
```

```
                  model npar AICc DeltaAICc weight Deviance
1          Phi(~1)p(~1)    2  671       0.0 0.5872       58
3        Phi(~sex)p(~1)    3  673       1.9 0.2309       84
5       Phi(~time)p(~1)    7  674       3.1 0.1227       77
4 Phi(~sex + time)p(~1)    8  676       5.1 0.0452       77
2       Phi(~1)p(~time)    7  679       7.9 0.0114       82
6    Phi(~time)p(~time)   12  682      10.8 0.0026       74
```

---
class: center, middle
background-size: cover

## Live demo #2: Continuous individual covariate

---
## Fit model with linear effect of wing length on survival

```r
phiwl <- list(formula=~wglength)
```

```r
phiwl.p <- mark(dipper.proc,
                dipper.ddl,
                model.parameters = list(Phi = phiwl, 
                                        p = p))
```

---
## Parameter estimates (1)

```

Output summary for CJS model
Name : Phi(~wglength)p(~1)

Npar :  3
-2lnL:  667
AICc :  673

Beta
                estimate    se    lcl  ucl
Phi:(Intercept)    0.242 0.102  0.042 0.44
Phi:wglength      -0.018 0.094 -0.202 0.17
p:(Intercept)      2.226 0.325  1.589 2.86

Real Parameter Phi
Group:sexF 
     1    2    3    4    5    6
1 0.56 0.56 0.56 0.56 0.56 0.56
2      0.56 0.56 0.56 0.56 0.56
3           0.56 0.56 0.56 0.56
4                0.56 0.56 0.56
5                     0.56 0.56
6                          0.56

Group:sexM 
     1    2    3    4    5    6
1 0.56 0.56 0.56 0.56 0.56 0.56
2      0.56 0.56 0.56 0.56 0.56
3           0.56 0.56 0.56 0.56
4                0.56 0.56 0.56
5                     0.56 0.56
6                          0.56

Group:sexM 
    2   3   4   5   6   7
1 0.9 0.9 0.9 0.9 0.9 0.9
2     0.9 0.9 0.9 0.9 0.9
3         0.9 0.9 0.9 0.9
4             0.9 0.9 0.9
5                 0.9 0.9
6                     0.9
```

---
## Parameter estimates (2)

```r
phiwl.p$results$beta
```

```
                estimate    se    lcl  ucl
Phi:(Intercept)    0.242 0.102  0.042 0.44
Phi:wglength      -0.018 0.094 -0.202 0.17
p:(Intercept)      2.226 0.325  1.589 2.86
```

---
## Graphical representation (1)

```r
# standardize wing length
scaled_wl <- scale(dipper_csv$wing_length)
# build a grid of 100 values between min and max of wing length
grid_wl <- seq(min(scaled_wl), max(scaled_wl), length = 100)
# predict survival on logit scale
logit_predicted_survival <- phiwl.p$results$beta[1,1] + 
                    phiwl.p$results$beta[2,1] * grid_wl
# back-transform to get predicted survival
predicted_survival <- plogis(logit_predicted_survival)
# create data.frame
df <- data.frame(winglength = grid_wl,
                 survival = predicted_survival)
```

Note: The two steps prediction and back-transformation can be done at once with `RMark` function `covariate.predictions()`.

---
## Graphical representation (2)

```r
# visualise
df %>% 
  ggplot() +
  aes(x = winglength, y = survival) + 
  geom_line()
```

---

---
## Quadratic relationship

```r
# square standardised covariate and add it to dipper data
dipper$wglengthsq <- scale(dipper_csv$wing_length)^2
# process data
dipper.proc <- process.data(data = dipper, 
                            groups = ("sex"))
# build design matrices
dipper.ddl <- make.design.data(dipper.proc)
```

---
## Model with quadratic effect of wing length on survival

```r
phiwl2 <- list(formula=~wglength + wglengthsq)
```

```r
phiwl2.p <- mark(dipper.proc,
                 dipper.ddl,
                 model.parameters = list(Phi = phiwl2, p = p))
```

---

```

Output summary for CJS model
Name : Phi(~wglength + wglengthsq)p(~1)

Npar :  4
-2lnL:  664
AICc :  672

Beta
                estimate    se    lcl  ucl
Phi:(Intercept)    0.017 0.162 -0.301 0.33
Phi:wglength      -0.026 0.096 -0.213 0.16
Phi:wglengthsq     0.206 0.117 -0.023 0.43
p:(Intercept)      2.231 0.325  1.594 2.87

Group:sexM 
    2   3   4   5   6   7
1 0.9 0.9 0.9 0.9 0.9 0.9
2     0.9 0.9 0.9 0.9 0.9
3         0.9 0.9 0.9 0.9
4             0.9 0.9 0.9
5                 0.9 0.9
6                     0.9
```

---

```r
phiwl2.p$results$beta
```

```
                estimate    se    lcl  ucl
Phi:(Intercept)    0.017 0.162 -0.301 0.33
Phi:wglength      -0.026 0.096 -0.213 0.16
Phi:wglengthsq     0.206 0.117 -0.023 0.43
p:(Intercept)      2.231 0.325  1.594 2.87
```

---

---

```r
# Build dataset
scaled_wl <- scale(dipper_csv$wing_length)
grid_wl <- seq(min(scaled_wl), max(scaled_wl), length = 100)
logit_predicted_survival <- phiwl2.p$results$beta[1,1] + phiwl2.p$results$beta[2,1] * grid_wl + phiwl2.p$results$beta[2,1] * grid_wl^2
predicted_survival <- plogis(logit_predicted_survival)
df <- data.frame(winglength = grid_wl,
                 survival = predicted_survival)

# Visualize
df %>% 
  ggplot() +
  aes(x = winglength, y = survival) + 
  geom_line() + 
  theme_light(base_size = 18)
```

---
## Model ranking with AIC

```r
collect.models()
```

```
                             model npar AICc DeltaAICc weight Deviance
1                     Phi(~1)p(~1)    2  671      0.00 0.3665       58
8 Phi(~wglength + wglengthsq)p(~1)    4  672      0.84 0.2404      664
3                   Phi(~sex)p(~1)    3  673      1.87 0.1441       84
7              Phi(~wglength)p(~1)    3  673      1.99 0.1355      667
5                  Phi(~time)p(~1)    7  674      3.13 0.0766       77
4            Phi(~sex + time)p(~1)    8  676      5.13 0.0282       77
2                  Phi(~1)p(~time)    7  679      7.88 0.0071       82
6               Phi(~time)p(~time)   12  682     10.84 0.0016       74
```

---
class: center, middle
background-size: cover

## Live demo #3: Environmental covariate

---
## Flood effect on survival?

+ A flood took place toward the end of the second survival interval and continue into the beginning of the third survival interval

+ Question is whether survival was lower in the two flood years

---
## Explore further model structure

```r
head(dipper.ddl$Phi)
```

```
  par.index model.index group cohort age time occ.cohort Cohort Age Time sex
1         1           1     F      1   0    1          1      0   0    0   F
2         2           2     F      1   1    2          1      0   1    1   F
3         3           3     F      1   2    3          1      0   2    2   F
4         4           4     F      1   3    4          1      0   3    3   F
5         5           5     F      1   4    5          1      0   4    4   F
6         6           6     F      1   5    6          1      0   5    5   F
```

---
## Create dummy Flood variable

```r
dipper.ddl$Phi$Flood <- 0
head(dipper.ddl$Phi)
```

```
  par.index model.index group cohort age time occ.cohort Cohort Age Time sex
1         1           1     F      1   0    1          1      0   0    0   F
2         2           2     F      1   1    2          1      0   1    1   F
3         3           3     F      1   2    3          1      0   2    2   F
4         4           4     F      1   3    4          1      0   3    3   F
5         5           5     F      1   4    5          1      0   4    4   F
6         6           6     F      1   5    6          1      0   5    5   F
  Flood
1     0
2     0
3     0
4     0
5     0
6     0
```

---
## Add flood effect

```r
dipper.ddl$Phi$Flood[dipper.ddl$Phi$time==2 | 
                       dipper.ddl$Phi$time==3] <- 1
head(dipper.ddl$Phi)
```

```
  par.index model.index group cohort age time occ.cohort Cohort Age Time sex
1         1           1     F      1   0    1          1      0   0    0   F
2         2           2     F      1   1    2          1      0   1    1   F
3         3           3     F      1   2    3          1      0   2    2   F
4         4           4     F      1   3    4          1      0   3    3   F
5         5           5     F      1   4    5          1      0   4    4   F
6         6           6     F      1   5    6          1      0   5    5   F
  Flood
1     0
2     1
3     1
4     0
5     0
6     0
```

---
## Fit model with flood effect

```r
phiFlood <- list(formula=~Flood)
```

```r
phiFlood.p <- mark(dipper.proc,
                   dipper.ddl,
                   model.parameters = list(Phi = phiFlood, 
                                           p = p))
```

---

```

Output summary for CJS model
Name : Phi(~Flood)p(~1)

Npar :  3
-2lnL:  660
AICc :  666

Beta
                estimate   se   lcl   ucl
Phi:(Intercept)     0.44 0.13  0.18  0.69
Phi:Flood          -0.56 0.22 -0.98 -0.14
p:(Intercept)       2.19 0.32  1.56  2.83

Real Parameter Phi
Group:sexF 
     1    2    3    4    5    6
1 0.61 0.47 0.47 0.61 0.61 0.61
2      0.47 0.47 0.61 0.61 0.61
3           0.47 0.61 0.61 0.61
4                0.61 0.61 0.61
5                     0.61 0.61
6                          0.61

Group:sexM 
     1    2    3    4    5    6
1 0.61 0.47 0.47 0.61 0.61 0.61
2      0.47 0.47 0.61 0.61 0.61
3           0.47 0.61 0.61 0.61
4                0.61 0.61 0.61
5                     0.61 0.61
6                          0.61

Group:sexM 
    2   3   4   5   6   7
1 0.9 0.9 0.9 0.9 0.9 0.9
2     0.9 0.9 0.9 0.9 0.9
3         0.9 0.9 0.9 0.9
4             0.9 0.9 0.9
5                 0.9 0.9
6                     0.9
```

---

```r
phiFlood.p$results$real
```

```
                estimate    se  lcl  ucl fixed note
Phi gF c1 a0 t1     0.61 0.031 0.55 0.67           
Phi gF c1 a1 t2     0.47 0.043 0.39 0.55           
p gF c1 a1 t2       0.90 0.029 0.83 0.94           
```

+ Clearly, survival in flood years was much lower than survival in non-flood years

---
## Model ranking with AIC

```r
collect.models()
```

```
                             model npar AICc DeltaAICc  weight Deviance
3                 Phi(~Flood)p(~1)    3  666       0.0 0.79405       78
1                     Phi(~1)p(~1)    2  671       4.7 0.07549       58
9 Phi(~wglength + wglengthsq)p(~1)    4  672       5.6 0.04951      664
4                   Phi(~sex)p(~1)    3  673       6.6 0.02968       84
8              Phi(~wglength)p(~1)    3  673       6.7 0.02790      667
6                  Phi(~time)p(~1)    7  674       7.8 0.01577       77
5            Phi(~sex + time)p(~1)    8  676       9.8 0.00581       77
2                  Phi(~1)p(~time)    7  679      12.6 0.00147       82
7               Phi(~time)p(~time)   12  682      15.5 0.00033       74
```

---
class: center, middle
background-size: cover

## Live demo #4: Age in capture-recapture models

---
## Build and fit a model with age-dependent survival

```r
phi.age <- list(formula =~ age) # already coded in RMark
```

```r
phiage.p <- mark(dipper.proc,
                 dipper.ddl,
                 model.parameters = list(Phi = phi.age, 
                                         p = p))
```

---

```

Output summary for CJS model
Name : Phi(~age)p(~1)

Npar :  7  (unadjusted=6)
-2lnL:  663
AICc :  677  (unadjusted=675.02107)

Beta
                estimate      se       lcl     ucl
Phi:(Intercept)    0.205    0.14    -0.068    0.48
Phi:age1           0.021    0.25    -0.476    0.52
Phi:age2           0.359    0.37    -0.365    1.08
Phi:age3          -0.058    0.54    -1.113    1.00
Phi:age4           0.240    0.89    -1.498    1.98
Phi:age5         -16.372 2027.27 -3989.814 3957.07
p:(Intercept)      2.255    0.33     1.611    2.90

Real Parameter Phi
Group:sexF 
     1    2    3    4    5       6
1 0.55 0.56 0.64 0.54 0.61 9.5e-08
2      0.55 0.56 0.64 0.54 6.1e-01
3           0.55 0.56 0.64 5.4e-01
4                0.55 0.56 6.4e-01
5                     0.55 5.6e-01
6                          5.5e-01

Group:sexM 
     1    2    3    4    5       6
1 0.55 0.56 0.64 0.54 0.61 9.5e-08
2      0.55 0.56 0.64 0.54 6.1e-01
3           0.55 0.56 0.64 5.4e-01
4                0.55 0.56 6.4e-01
5                     0.55 5.6e-01
6                          5.5e-01

Real Parameter p
Group:sexF 
     2    3    4    5    6    7
1 0.91 0.91 0.91 0.91 0.91 0.91
2      0.91 0.91 0.91 0.91 0.91
3           0.91 0.91 0.91 0.91
4                0.91 0.91 0.91
5                     0.91 0.91
6                          0.91

Group:sexM 
     2    3    4    5    6    7
1 0.91 0.91 0.91 0.91 0.91 0.91
2      0.91 0.91 0.91 0.91 0.91
3           0.91 0.91 0.91 0.91
4                0.91 0.91 0.91
5                     0.91 0.91
6                          0.91
```

---
## Parameter estimates

```r
phiage.p$results$real
```

```
                estimate      se      lcl  ucl fixed note
Phi gF c1 a0 t1  5.5e-01 0.03447  4.8e-01 0.62           
Phi gF c1 a1 t2  5.6e-01 0.05027  4.6e-01 0.65           
Phi gF c1 a2 t3  6.4e-01 0.07948  4.7e-01 0.78           
Phi gF c1 a3 t4  5.4e-01 0.12937  2.9e-01 0.76           
Phi gF c1 a4 t5  6.1e-01 0.20896  2.2e-01 0.90           
Phi gF c1 a5 t6  9.5e-08 0.00019 5.3e-316 1.00           
p gF c1 a1 t2    9.1e-01 0.02825  8.3e-01 0.95           
```

---
## What if we'd like only 2 age classes?

+ Create a (0, 1+) age variable:

```r
dipper.ddl <- add.design.data(dipper.proc,
                              dipper.ddl, 
                              "Phi",
                              type = "age",
                              bins = c(0, 1, 7), # (0, 1+)
                              name = "ageclass") # var name
```

---
## Fit model with 2 age classes on survival

```r
phi.age2 <- list(formula=~ageclass) # age effect on survival
```

```r
phiage2.p <- mark(dipper.proc,
                  dipper.ddl,
                  model.parameters = list(Phi = phi.age2, 
                                          p = p))
```

---

```

Output summary for CJS model
Name : Phi(~ageclass)p(~1)

Npar :  3
-2lnL:  667
AICc :  673

Beta
                  estimate   se     lcl  ucl
Phi:(Intercept)       0.22 0.11 -0.0038 0.44
Phi:ageclass(1,7]     0.15 0.28 -0.3957 0.70
p:(Intercept)         2.24 0.33  1.6001 2.87

Real Parameter Phi
Group:sexF 
     1    2    3    4    5    6
1 0.55 0.55 0.59 0.59 0.59 0.59
2      0.55 0.55 0.59 0.59 0.59
3           0.55 0.55 0.59 0.59
4                0.55 0.55 0.59
5                     0.55 0.55
6                          0.55

Group:sexM 
     1    2    3    4    5    6
1 0.55 0.55 0.59 0.59 0.59 0.59
2      0.55 0.55 0.59 0.59 0.59
3           0.55 0.55 0.59 0.59
4                0.55 0.55 0.59
5                     0.55 0.55
6                          0.55

Group:sexM 
    2   3   4   5   6   7
1 0.9 0.9 0.9 0.9 0.9 0.9
2     0.9 0.9 0.9 0.9 0.9
3         0.9 0.9 0.9 0.9
4             0.9 0.9 0.9
5                 0.9 0.9
6                     0.9
```

---
## Parameter estimates

```r
phiage2.p$results$real
```

```
                estimate    se  lcl  ucl fixed note
Phi gF c1 a0 t1     0.55 0.028 0.50 0.61           
Phi gF c1 a2 t3     0.59 0.062 0.47 0.70           
p gF c1 a1 t2       0.90 0.028 0.83 0.95           
```

---
class: center, middle
background-size: cover

## Live demo #5: Multisite capture-recapture models

---

## Read in data (1)

Read in csv file:

```r
geese_csv <- read_csv2("dat/allgeese.csv")
```

---
## Read in data (2)

```r
geese_csv
```

```
# A tibble: 623 × 7
   year_1984 year_1985 year_1986 year_1987 year_1988 year_1989     N
       <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl> <dbl>
 1         0         0         0         0         0         1   158
 2         0         0         0         0         0         2   352
 3         0         0         0         0         0         3   271
 4         0         0         0         0         1         0   317
 5         0         0         0         0         1         1    62
 6         0         0         0         0         1         2    31
 7         0         0         0         0         2         0   748
 8         0         0         0         0         2         1    56
 9         0         0         0         0         2         2   234
10         0         0         0         0         2         3     7
# … with 613 more rows
```

---
## Convert to Mark format (1)

Build a data.frame with encounter histories, and number of individuals:

```r
geese <- data.frame(ch = unite(geese_csv %>% 
                                 select(year_1984:year_1989), 
                               col = "ch", 
                               sep = ""), 
                    freq = geese_csv$N)
```

---
## Convert to Mark format (2)

Inspect first lines:

```r
head(geese)
```

```
      ch freq
1 000001  158
2 000002  352
3 000003  271
4 000010  317
5 000011   62
6 000012   31
```

---
## Preliminaries

```r
geese.proc <- process.data(data = geese,
                           model = "Multistrata")
```

```r
geese.ddl <- make.design.data(geese.proc)
```

```r
phi <- list(formula=~1) # survival constant 
phi.site <- list(formula=~stratum) # site-specific survival 
p <- list(formula=~1) # constant detection
psi.site <- list(formula=~-1+stratum:tostratum,
                 link='mlogit') # movements
```

---
## Run model with site effect on survival:

```r
phis.psi.p <- mark(geese.proc,
                   geese.ddl,
                   model.parameters=list(S = phi.site,
                                         p = p,
                                         Psi = psi.site))
```

---

```

Output summary for Multistrata model
Name : S(~stratum)p(~1)Psi(~-1 + stratum:tostratum)

Npar :  10
-2lnL:  73743
AICc :  73763

Beta
                        estimate    se    lcl    ucl
S:(Intercept)              0.720 0.031  0.659  0.782
S:stratum2                 0.049 0.038 -0.025  0.123
S:stratum3                -0.153 0.049 -0.249 -0.058
p:(Intercept)             -0.336 0.020 -0.374 -0.298
Psi:stratum2:tostratum1   -1.988 0.035 -2.057 -1.920
Psi:stratum3:tostratum1   -2.501 0.101 -2.699 -2.302
Psi:stratum1:tostratum2   -1.141 0.037 -1.213 -1.069
Psi:stratum3:tostratum2   -0.873 0.049 -0.969 -0.777
Psi:stratum1:tostratum3   -4.950 0.234 -5.408 -4.492
Psi:stratum2:tostratum3   -3.666 0.078 -3.818 -3.514

Real Parameter S
 Stratum:1 
     1    2    3    4    5
1 0.67 0.67 0.67 0.67 0.67
2      0.67 0.67 0.67 0.67
3           0.67 0.67 0.67
4                0.67 0.67
5                     0.67

Stratum:2 
     1    2    3    4    5
1 0.68 0.68 0.68 0.68 0.68
2      0.68 0.68 0.68 0.68
3           0.68 0.68 0.68
4                0.68 0.68
5                     0.68

Stratum:3 
     1    2    3    4    5
1 0.64 0.64 0.64 0.64 0.64
2      0.64 0.64 0.64 0.64
3           0.64 0.64 0.64
4                0.64 0.64
5                     0.64

Real Parameter p
 Stratum:1 
     2    3    4    5    6
1 0.42 0.42 0.42 0.42 0.42
2      0.42 0.42 0.42 0.42
3           0.42 0.42 0.42
4                0.42 0.42
5                     0.42

Stratum:2 
     2    3    4    5    6
1 0.42 0.42 0.42 0.42 0.42
2      0.42 0.42 0.42 0.42
3           0.42 0.42 0.42
4                0.42 0.42
5                     0.42

Stratum:3 
     2    3    4    5    6
1 0.42 0.42 0.42 0.42 0.42
2      0.42 0.42 0.42 0.42
3           0.42 0.42 0.42
4                0.42 0.42
5                     0.42

Real Parameter Psi
 Stratum:1 To:2 
     1    2    3    4    5
1 0.24 0.24 0.24 0.24 0.24
2      0.24 0.24 0.24 0.24
3           0.24 0.24 0.24
4                0.24 0.24
5                     0.24

Stratum:1 To:3 
       1      2      3      4      5
1 0.0053 0.0053 0.0053 0.0053 0.0053
2        0.0053 0.0053 0.0053 0.0053
3               0.0053 0.0053 0.0053
4                      0.0053 0.0053
5                             0.0053

Stratum:2 To:1 
     1    2    3    4    5
1 0.12 0.12 0.12 0.12 0.12
2      0.12 0.12 0.12 0.12
3           0.12 0.12 0.12
4                0.12 0.12
5                     0.12

Stratum:2 To:3 
      1     2     3     4     5
1 0.022 0.022 0.022 0.022 0.022
2       0.022 0.022 0.022 0.022
3             0.022 0.022 0.022
4                   0.022 0.022
5                         0.022

Stratum:3 To:1 
      1     2     3     4     5
1 0.055 0.055 0.055 0.055 0.055
2       0.055 0.055 0.055 0.055
3             0.055 0.055 0.055
4                   0.055 0.055
5                         0.055

Stratum:3 To:2 
     1    2    3    4    5
1 0.28 0.28 0.28 0.28 0.28
2      0.28 0.28 0.28 0.28
3           0.28 0.28 0.28
4                0.28 0.28
5                     0.28
```

---

```r
phis.psi.p$results$real
```

```
                          estimate     se    lcl    ucl fixed note
S s1 g1 c1 a0 o1 t1         0.6727 0.0069 0.6590 0.6860           
S s2 g1 c1 a0 o1 t1         0.6834 0.0051 0.6734 0.6932           
S s3 g1 c1 a0 o1 t1         0.6381 0.0095 0.6193 0.6564           
p s1 g1 c1 a1 o1 t2         0.4167 0.0047 0.4075 0.4261           
Psi s1 to2 g1 c1 a0 o1 t1   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c1 a1 o2 t2   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c1 a2 o3 t3   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c1 a3 o4 t4   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c1 a4 o5 t5   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c2 a0 o2 t2   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c2 a1 o3 t3   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c2 a2 o4 t4   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c2 a3 o5 t5   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c3 a0 o3 t3   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c3 a1 o4 t4   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c3 a2 o5 t5   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c4 a0 o4 t4   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c4 a1 o5 t5   0.2409 0.0067 0.2279 0.2543           
Psi s1 to2 g1 c5 a0 o5 t5   0.2409 0.0067 0.2279 0.2543           
Psi s1 to3 g1 c1 a0 o1 t1   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c1 a1 o2 t2   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c1 a2 o3 t3   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c1 a3 o4 t4   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c1 a4 o5 t5   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c2 a0 o2 t2   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c2 a1 o3 t3   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c2 a2 o4 t4   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c2 a3 o5 t5   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c3 a0 o3 t3   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c3 a1 o4 t4   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c3 a2 o5 t5   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c4 a0 o4 t4   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c4 a1 o5 t5   0.0053 0.0012 0.0034 0.0084           
Psi s1 to3 g1 c5 a0 o5 t5   0.0053 0.0012 0.0034 0.0084           
Psi s2 to1 g1 c1 a0 o1 t1   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c1 a1 o2 t2   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c1 a2 o3 t3   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c1 a3 o4 t4   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c1 a4 o5 t5   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c2 a0 o2 t2   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c2 a1 o3 t3   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c2 a2 o4 t4   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c2 a3 o5 t5   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c3 a0 o3 t3   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c3 a1 o4 t4   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c3 a2 o5 t5   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c4 a0 o4 t4   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c4 a1 o5 t5   0.1178 0.0036 0.1108 0.1251           
Psi s2 to1 g1 c5 a0 o5 t5   0.1178 0.0036 0.1108 0.1251           
Psi s2 to3 g1 c1 a0 o1 t1   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c1 a1 o2 t2   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c1 a2 o3 t3   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c1 a3 o4 t4   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c1 a4 o5 t5   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c2 a0 o2 t2   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c2 a1 o3 t3   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c2 a2 o4 t4   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c2 a3 o5 t5   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c3 a0 o3 t3   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c3 a1 o4 t4   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c3 a2 o5 t5   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c4 a0 o4 t4   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c4 a1 o5 t5   0.0220 0.0017 0.0190 0.0255           
Psi s2 to3 g1 c5 a0 o5 t5   0.0220 0.0017 0.0190 0.0255           
Psi s3 to1 g1 c1 a0 o1 t1   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c1 a1 o2 t2   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c1 a2 o3 t3   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c1 a3 o4 t4   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c1 a4 o5 t5   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c2 a0 o2 t2   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c2 a1 o3 t3   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c2 a2 o4 t4   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c2 a3 o5 t5   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c3 a0 o3 t3   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c3 a1 o4 t4   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c3 a2 o5 t5   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c4 a0 o4 t4   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c4 a1 o5 t5   0.0547 0.0052 0.0453 0.0659           
Psi s3 to1 g1 c5 a0 o5 t5   0.0547 0.0052 0.0453 0.0659           
Psi s3 to2 g1 c1 a0 o1 t1   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c1 a1 o2 t2   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c1 a2 o3 t3   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c1 a3 o4 t4   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c1 a4 o5 t5   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c2 a0 o2 t2   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c2 a1 o3 t3   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c2 a2 o4 t4   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c2 a3 o5 t5   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c3 a0 o3 t3   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c3 a1 o4 t4   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c3 a2 o5 t5   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c4 a0 o4 t4   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c4 a1 o5 t5   0.2786 0.0098 0.2597 0.2982           
Psi s3 to2 g1 c5 a0 o5 t5   0.2786 0.0098 0.2597 0.2982           
```

---

```r
Psilist <- get.real(phis.psi.p, "Psi", vcv = TRUE)
Psi.values = Psilist$estimates
top.psi = TransitionMatrix(Psi.values[Psi.values$time == 1, ], 
                           vcv.real = Psilist$vcv.real)
top.psi
```

```
$TransitionMat
      1    2      3
1 0.754 0.24 0.0053
2 0.118 0.86 0.0219
3 0.055 0.28 0.6662

$se.TransitionMat
       1      2      3
1 0.0067 0.0067 0.0012
2 0.0036 0.0039 0.0017
3 0.0052 0.0098 0.0101

$lcl.TransitionMat
      1    2      3
1 0.740 0.23 0.0034
2 0.111 0.85 0.0189
3 0.045 0.26 0.6461

$ucl.TransitionMat
      1    2      3
1 0.767 0.25 0.0084
2 0.125 0.87 0.0254
3 0.066 0.30 0.6857
```

---
## Run model with state effect on survival

```r
phi.psi.p <- mark(geese.proc,
                  geese.ddl,
                  model.parameters=list(S = phi,
                                        p = p,
                                        Psi = psi.site))
```

---

```

Output summary for Multistrata model
Name : S(~1)p(~1)Psi(~-1 + stratum:tostratum)

Npar :  8
-2lnL:  73762
AICc :  73778

Beta
                        estimate    se   lcl   ucl
S:(Intercept)               0.73 0.018  0.69  0.76
p:(Intercept)              -0.34 0.019 -0.38 -0.30
Psi:stratum2:tostratum1    -1.99 0.035 -2.06 -1.92
Psi:stratum3:tostratum1    -2.49 0.101 -2.68 -2.29
Psi:stratum1:tostratum2    -1.14 0.037 -1.21 -1.07
Psi:stratum3:tostratum2    -0.86 0.049 -0.95 -0.76
Psi:stratum1:tostratum3    -4.96 0.233 -5.41 -4.50
Psi:stratum2:tostratum3    -3.68 0.078 -3.83 -3.53

Real Parameter S
 Stratum:1 
     1    2    3    4    5
1 0.67 0.67 0.67 0.67 0.67
2      0.67 0.67 0.67 0.67
3           0.67 0.67 0.67
4                0.67 0.67
5                     0.67

Stratum:2 
     1    2    3    4    5
1 0.67 0.67 0.67 0.67 0.67
2      0.67 0.67 0.67 0.67
3           0.67 0.67 0.67
4                0.67 0.67
5                     0.67

Stratum:3 
     1    2    3    4    5
1 0.67 0.67 0.67 0.67 0.67
2      0.67 0.67 0.67 0.67
3           0.67 0.67 0.67
4                0.67 0.67
5                     0.67

Real Parameter p
 Stratum:1 
     2    3    4    5    6
1 0.42 0.42 0.42 0.42 0.42
2      0.42 0.42 0.42 0.42
3           0.42 0.42 0.42
4                0.42 0.42
5                     0.42

Stratum:2 
     2    3    4    5    6
1 0.42 0.42 0.42 0.42 0.42
2      0.42 0.42 0.42 0.42
3           0.42 0.42 0.42
4                0.42 0.42
5                     0.42

Stratum:3 
     2    3    4    5    6
1 0.42 0.42 0.42 0.42 0.42
2      0.42 0.42 0.42 0.42
3           0.42 0.42 0.42
4                0.42 0.42
5                     0.42

Real Parameter Psi
 Stratum:1 To:2 
     1    2    3    4    5
1 0.24 0.24 0.24 0.24 0.24
2      0.24 0.24 0.24 0.24
3           0.24 0.24 0.24
4                0.24 0.24
5                     0.24

Stratum:2 To:1 
     1    2    3    4    5
1 0.12 0.12 0.12 0.12 0.12
2      0.12 0.12 0.12 0.12
3           0.12 0.12 0.12
4                0.12 0.12
5                     0.12

Stratum:3 To:2 
     1    2    3    4    5
1 0.28 0.28 0.28 0.28 0.28
2      0.28 0.28 0.28 0.28
3           0.28 0.28 0.28
4                0.28 0.28
5                     0.28
```

---

```r
phi.psi.p$results$real
```

```
                          estimate     se    lcl    ucl fixed note
S s1 g1 c1 a0 o1 t1         0.6745 0.0040 0.6665 0.6823           
p s1 g1 c1 a1 o1 t2         0.4156 0.0047 0.4063 0.4249           
Psi s1 to2 g1 c1 a0 o1 t1   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c1 a1 o2 t2   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c1 a2 o3 t3   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c1 a3 o4 t4   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c1 a4 o5 t5   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c2 a0 o2 t2   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c2 a1 o3 t3   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c2 a2 o4 t4   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c2 a3 o5 t5   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c3 a0 o3 t3   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c3 a1 o4 t4   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c3 a2 o5 t5   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c4 a0 o4 t4   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c4 a1 o5 t5   0.2416 0.0067 0.2287 0.2550           
Psi s1 to2 g1 c5 a0 o5 t5   0.2416 0.0067 0.2287 0.2550           
Psi s1 to3 g1 c1 a0 o1 t1   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c1 a1 o2 t2   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c1 a2 o3 t3   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c1 a3 o4 t4   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c1 a4 o5 t5   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c2 a0 o2 t2   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c2 a1 o3 t3   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c2 a2 o4 t4   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c2 a3 o5 t5   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c3 a0 o3 t3   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c3 a1 o4 t4   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c3 a2 o5 t5   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c4 a0 o4 t4   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c4 a1 o5 t5   0.0053 0.0012 0.0034 0.0083           
Psi s1 to3 g1 c5 a0 o5 t5   0.0053 0.0012 0.0034 0.0083           
Psi s2 to1 g1 c1 a0 o1 t1   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c1 a1 o2 t2   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c1 a2 o3 t3   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c1 a3 o4 t4   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c1 a4 o5 t5   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c2 a0 o2 t2   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c2 a1 o3 t3   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c2 a2 o4 t4   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c2 a3 o5 t5   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c3 a0 o3 t3   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c3 a1 o4 t4   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c3 a2 o5 t5   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c4 a0 o4 t4   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c4 a1 o5 t5   0.1175 0.0036 0.1106 0.1248           
Psi s2 to1 g1 c5 a0 o5 t5   0.1175 0.0036 0.1106 0.1248           
Psi s2 to3 g1 c1 a0 o1 t1   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c1 a1 o2 t2   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c1 a2 o3 t3   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c1 a3 o4 t4   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c1 a4 o5 t5   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c2 a0 o2 t2   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c2 a1 o3 t3   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c2 a2 o4 t4   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c2 a3 o5 t5   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c3 a0 o3 t3   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c3 a1 o4 t4   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c3 a2 o5 t5   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c4 a0 o4 t4   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c4 a1 o5 t5   0.0217 0.0016 0.0187 0.0252           
Psi s2 to3 g1 c5 a0 o5 t5   0.0217 0.0016 0.0187 0.0252           
Psi s3 to1 g1 c1 a0 o1 t1   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c1 a1 o2 t2   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c1 a2 o3 t3   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c1 a3 o4 t4   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c1 a4 o5 t5   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c2 a0 o2 t2   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c2 a1 o3 t3   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c2 a2 o4 t4   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c2 a3 o5 t5   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c3 a0 o3 t3   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c3 a1 o4 t4   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c3 a2 o5 t5   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c4 a0 o4 t4   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c4 a1 o5 t5   0.0552 0.0053 0.0458 0.0665           
Psi s3 to1 g1 c5 a0 o5 t5   0.0552 0.0053 0.0458 0.0665           
Psi s3 to2 g1 c1 a0 o1 t1   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c1 a1 o2 t2   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c1 a2 o3 t3   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c1 a3 o4 t4   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c1 a4 o5 t5   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c2 a0 o2 t2   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c2 a1 o3 t3   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c2 a2 o4 t4   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c2 a3 o5 t5   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c3 a0 o3 t3   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c3 a1 o4 t4   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c3 a2 o5 t5   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c4 a0 o4 t4   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c4 a1 o5 t5   0.2817 0.0099 0.2628 0.3014           
Psi s3 to2 g1 c5 a0 o5 t5   0.2817 0.0099 0.2628 0.3014           
```

---

```r
Psilist <- get.real(phi.psi.p, "Psi", vcv = TRUE)
Psi.values = Psilist$estimates
top.psi = TransitionMatrix(Psi.values[Psi.values$time == 1, ], 
                           vcv.real = Psilist$vcv.real)
top.psi
```

```
$TransitionMat
      1    2      3
1 0.754 0.24 0.0053
2 0.118 0.86 0.0217
3 0.055 0.28 0.6640

$se.TransitionMat
       1      2      3
1 0.0067 0.0067 0.0012
2 0.0036 0.0039 0.0016
3 0.0052 0.0098 0.0101

$lcl.TransitionMat
      1    2      3
1 0.740 0.23 0.0034
2 0.111 0.85 0.0187
3 0.046 0.26 0.6440

$ucl.TransitionMat
      1    2      3
1 0.767 0.25 0.0083
2 0.125 0.87 0.0251
3 0.066 0.30 0.6834
```

---

```r
collect.models()
```

```
                                         model npar  AICc DeltaAICc  weight
2 S(~stratum)p(~1)Psi(~-1 + stratum:tostratum)   10 73763         0 0.99952
1       S(~1)p(~1)Psi(~-1 + stratum:tostratum)    8 73778        15 0.00048
  Deviance
2     2587
1     2607
```