---
title: "Pied Flycatchers 3: Joint Effects"
author: "Simon Bonner"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Pied Flycatchers 3: Joint Effects}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

Joint effects are fixed or random effects that take the same value in linear predictors describing variation in both the mean and the dispersion parameter. As an example, we consider the simplest model of the pied flycatcher data with only fixed effects on the mean and dispersion, except that we assume that the effect of broodsize is the same in both linear predictors. 

# Library

First you have to load the package:
```{r}
## Load package
library(dalmatian)
```

# Data

```{r}
## Load pied flycatcher data
data(pied_flycatchers_1)

## Create variables bounding the true load
pfdata$lower=ifelse(pfdata$load==0,log(.001),log(pfdata$load-.049))
pfdata$upper=log(pfdata$load+.05)
```

# Defining the Model
Our model will assume that the logarithm of the load returned on the $j$th trip by adult $i$ is normally distributed with mean $\mu_{ij}$ and variance $\phi$ where:
$$\mu=\alpha_0 + \alpha_1 \mathrm{log(IVI)}_{ij} + \beta_2 \mathrm{sex}_i + \gamma_1 \mathrm{broodsize}_i$$
and
$$\log(\phi)=\psi_0 + \psi_1\mathrm{sex}_i$ + \gamma_1 \mathrm{broodsize}_i$.
Here $\gamma_1$ is a fixed effect that is common to both the linear predictors describing the variation in both the mean and dispersion. This model can be defined with the following three components:
```{r}
## Mean model
mymean <- list(fixed=list(name="alpha",
                          formula=~ log(IVI) + sex,
                          priors=list(c("dnorm",0,.001))))

## Dispersion model
mydisp=list(fixed=list(name="psi",
                       formula=~sex,
                       priors=list(c("dnorm",0,.001))),
            link="log")

## Joint components
myjoint <- list(fixed = list(name = "gamma",
                            formula = ~-1 + broodsize,
                            priors = list(c("dnorm",0,.001))))

```
These three objects will now be used to generate the JAGS code, data, and initial values for running the model. Note that a link function cannot be specified as part of the joint component. 

# Running the Model 

The model can then be fit with `dalmatian` by defining the arguments that control the construction and updating of the JAGS model. Full details are included in `vignette(pied-flycatchers-1)`. 
```{r, eval = FALSE}

## Set working directory
## By default uses a system temp directory. You probably want to change this.
workingDir <- tempdir()

## Define list of arguments for jags.model()
jm.args <- list(file=file.path(workingDir,"pied_flycatcher_joint_1_jags.R"),n.adapt=1000)

## Define list of arguments for coda.samples()
cs.args <- list(n.iter=5000,thin=20)

## Run the model using dalmatian
pfmcmc4 <- dalmatian(df=pfdata,
                     mean.model=mymean,
                     dispersion.model=mydisp,
                     joint.model=myjoint,
                     jags.model.args=jm.args,
                     coda.samples.args=cs.args,
                     rounding=TRUE,
                     lower="lower",
                     upper="upper",
                     residuals = FALSE,
                     debug=FALSE)
```

```{r, echo = FALSE}
load(system.file("pfresults4.RData",package="dalmatian"))
```

# Results

The following illustrates the functions for processing the results. Again, these are fully described in `vignette(pied-flycatchers-1)`.

## Numerical Convergence Diagnostics
```{r, eval = FALSE}
## Compute convergence diagnostics
pfconvergence4 <- convergence(pfmcmc4)

```

```{r}
## Gelman-Rubin diagnostics
pfconvergence4$gelman

## Raftery diagnostics
pfconvergence4$raftery

## Effective sample size
pfconvergence4$effectiveSize
```

## Graphical Convergence Diagnostics
```{r, eval = FALSE}
## Generate traceplots
pftraceplots4 <- traceplots(pfmcmc4,show=FALSE,nthin=10)

```

```{r,fig.width=6,fig.align="center"}
## Fixed effects for mean
pftraceplots4$meanFixed

## Fixed effects for dispersion
pftraceplots4$dispersionFixed

```

## Numerical Summaries
```{r, eval = FALSE}
## Compute numerical summaries
pfsummary4 <- summary(pfmcmc4)
```

```{r}
## Print numerical summaries
pfsummary4
```

## Graphical Summaries
```{r, eval = FALSE}
## Generate caterpillar
pfcaterpillar4 <- caterpillar(pfmcmc4,show = FALSE)

```

```{r,fig.width=6,fig.align="center"}
## Fixed effects for mean
pfcaterpillar4$meanFixed

## Fixed effects for dispersion
pfcaterpillar4$dispersionFixed
```
