`getting-started.Rmd`

To install `moultmcmc`

please follow the instructions in the README file.

For demonstration purposes we use the `sanderlings`

dataset from Underhill and Zucchini (1988). It contains observations of sampling dates `Day`

and moult indices `MIndex`

for 164 adult Sanderlings (*Calidris alba*) trapped on 11 days in the South Africa in the austral summer of 1978/79. The sample contains 85 pre-moult individuals, 66 in active moult and 13 post-moult individuals. The data are augmented with a column describing categorical moult status `MCat`

, so we can demonstrate fitting Type 1 models when data are coded categorically.

```
data("sanderlings")
sanderlings$MCat <- case_when(sanderlings$MIndex == 0 ~ 1,
sanderlings$MIndex == 1 ~3,
TRUE ~ 2)
#plot the data
ggplot(sanderlings, aes(y=MIndex, x=Day, col = factor(MCat))) + geom_point() +theme_classic()
```

The central function for model fitting is the `moultmcmc()`

function. Model variants of the standard Underhill-Zucchini models are specified using the `type`

argument to designate the relevant input data type. Additionally the user has to designate the input data columns containing moult records, and dates, respectively, as well as linear predictor structures for the main moult parameters - the start date, duration and the standard deviation of the start date. To fit an intercept-only type 1 model we therefore use the following specification:

```
mcmc1 = moultmcmc(moult_column = "MCat",
date_column = "Day",
start_formula = ~1,
duration_formula = ~1,
sigma_formula = ~1,
data = sanderlings,
type=1,
init = "auto",
chains = 2,
refresh=1000)
#>
#> SAMPLING FOR MODEL 'uz1_linpred' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000134 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.34 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.788 seconds (Warm-up)
#> Chain 1: 0.346 seconds (Sampling)
#> Chain 1: 1.134 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'uz1_linpred' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 5.3e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.53 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.702 seconds (Warm-up)
#> Chain 2: 0.339 seconds (Sampling)
#> Chain 2: 1.041 seconds (Total)
#> Chain 2:
```

Most MCMC parameters, such as the number of chains (here 2), iterations, or cores are set in the same way as for `rstan::sampling`

, however, for initial values `moultmcmc`

defaults to an additional option `"auto"`

which will provide initial values for the MCMC chains based on the input data. Custom initial values can be specified instead as detailed in the documentation of `rstan::sampling`

.

Once the MCMC has completed, a summary table of parameter estimates can be displayed with

```
summary_table(mcmc1)
#> # A tibble: 5 × 11
#> parame…¹ estim…² se_mean sd lci uci n_eff Rhat prob method type
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 mean_(I… 134. 0.0881 3.19 128. 140. 1314. 1.00 0.95 MCMC 1
#> 2 duratio… 97.6 0.293 10.8 77.4 120. 1352. 1.00 0.95 MCMC 1
#> 3 log_sd_… 3.33 0.00651 0.228 2.86 3.76 1227. 1.00 0.95 MCMC 1
#> 4 sd_(Int… 28.7 0.180 6.45 17.5 42.8 1278. 1.00 0.95 MCMC 1
#> 5 lp__ -104. 0.0363 1.14 -107. -103. 987. 1.00 0.95 MCMC 1
#> # … with abbreviated variable names ¹parameter, ²estimate
```

and standard model assessment can be done as for any posterior sample from Stan, by directly accessing the `stanfit`

slot of the returned S3 object. For example we can look at the traceplots for the model using `rstan::stan_trace`

`stan_trace(mcmc1$stanfit)`

We can also plot the model fit against the data using the `moult_plot`

function

`moult_plot(mcmc1)`

For comparison we can fit the same model using maximum likelihood estimation via the `moult`

package (Erni et al. 2013) and compare parameter estimates visually with the `compare_plot`

function.

```
ml1 = moult(MIndex ~ Day,data = sanderlings, type = 1)
compare_plot(ml1, mcmc1)
```

The comparison plot shows very similar estimates for both methods.

Similarly we can fit the corresponding Type 2 model, that makes use of the full information contained in the moult indices, by setting `type=2`

. All `moultmcmc`

models default to an intercept only model, so we do not need to explicitly specify the individual linear predictors here.

```
#fit using MCMC
mcmc2 = moultmcmc(moult_column = "MIndex",
date_column = "Day",
type=2,
data = sanderlings,
chains = 2,
refresh = 1000)
#>
#> SAMPLING FOR MODEL 'uz2_linpred' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 6.3e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.63 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.693 seconds (Warm-up)
#> Chain 1: 0.248 seconds (Sampling)
#> Chain 1: 0.941 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'uz2_linpred' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 4.3e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.501 seconds (Warm-up)
#> Chain 2: 0.243 seconds (Sampling)
#> Chain 2: 0.744 seconds (Total)
#> Chain 2:
#fit using ML
ml2 = moult(MIndex ~ Day,data = sanderlings, type = 2)
#compare parameter estimates
compare_plot(ml2, mcmc2)
```

As before the comparison plot shows very similar estimates for both methods. Perhaps more interesting is a comparison of the two model types, which illustrates the gain in precision achieved by exploiting the information contained in the moult indices. Note that models can be given custom names for visualisation purposes using the `names`

argument to `compare_plot`

`compare_plot(mcmc1, mcmc2, names = c('UZ1','UZ2'))`

Erni, Birgit, Bo T. Bonnevie, Hans-Dieter Oschadleus, Res Altwegg, and Les G. Underhill. 2013. “Moult: An R Package to Analyze Moult in Birds.” *Journal of Statistical Software* 52. https://doi.org/10.18637/jss.v052.i08.

Underhill, Les G., and Walter Zucchini. 1988. “A Model for Avian Primary Moult.” *Ibis* 130: 358–72. https://doi.org/10.1111/j.1474-919x.1988.tb00993.x.