# Analysis of the Pied Flycatcher Data: Random Slopes

#### Simon Bonner

#### 2018-01-28

This vignette continues the example of the pied flycatchers discussed in the vignette “pied-flycatchers-1” by adding a random slope to the fixed effects model. I highly recommend that you work through the other vignette first.

# Library

First you have to load the package:

    ## Load package
    library(dalmatian)

# Raw data

The raw data for this example is provided in the package and can be accessed with:

    ## Load pied flycatcher data
    data(pied_flycatchers_1)

As before we will consider load as a response variable rounded to the nearest whole number and define the lower and upper bounds:

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

    ## Warning in log(pfdata$load - 0.049): NaNs produced

    pfdata$upper=log(pfdata$load+.05)

# Model

The model we consider will be the similar to the random effects model in “pied-flycatchers-1”. We will modify this model by including a random, individual slope for the effect of \(log(IVI)\) as well as the random intercept in the mean model. The new model for the mean will be: \[\mu=\beta_0 + \beta_1 \mathrm{log(IVI)_{ij}} + \beta_2 \mathrm{broodsize}_i + \beta_3 \mathrm{sex}_i + \epsilon_{1i} + \epsilon_{2i} \mathrm{log(IVI)_{ij}}\] For convenience, we will also simplify the model of the variance by assuming that the variance is constant across all observations. This can be done by setting: \[\log(\sigma)=\psi_0\].

The objects defining the new mean and variance components are specified as:

    # Random component of mean
    mymean=list(fixed=list(name="alpha",
                    formula=~ log(IVI) + broodsize + sex,
                    priors=list(c("dnorm",0,.001))),
        random=list(name="epsilon",formula=~-1 + indidx + indidx:log(IVI)))

    # Random component of variance
    myvar=list(fixed=list(name="psi",
                   link="log",
                   formula=~1,
                   priors=list(c("dnorm",0,.001))))

## Running the Model with dalmatian

The new model can now be run using `dalmatian` in exactly the same way as above. However, we will make one change. In order to shorten the convergence time, we will base initial values on the results of the fixed effects model and then provide these as input. In particular, we will define initial values for the fixed effects of both the mean and variance components by taking the values from the final iterations of each of the chains run for the fixed effects model. For convenience, we set the random effects variances equal to 1 in all three chains. Note that the chains have been run for a relatively small number of iterations and heavily subsampled to satisfy the size requirements for packages on . Neither of these is recommended in practice.

    ## 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_3_jags.R"),n.adapt=1000)

    ## Define list of arguments for coda.samples()
    cs.args <- list(n.iter=1000)

    ## Run the model using dalmatian
    if(runModels){
      pfresults3 <- dalmatian(df=pfdata,
                              mean.model=mymean,
                              variance.model=myvar,
                              jags.model.args=jm.args,
                              coda.samples.args=cs.args,
                              rounding=TRUE,
                              lower="lower",
                              upper="upper",
                              debug=FALSE)

      save(pfresults2,"pfresults3.RData")
    }
    if(!runModels){
      ## Load output from previously run chain
      load(system.file("Pied_Flycatchers_2","pfresults3.RData",package="dalmatian"))
    }

## Results

Once again, we can use the wrapper functions provided with `dalmatian` to conduct convergence diagnostics and compute summary statistics for the posterior distribution from which we can make inference.

1. Convergence diagnostics

    ## Compute convergence diagnostics
    pfconvergence2 <- convergence(pfresults3,raftery=list(r=.01))

    ## Gelman-Rubin diagnostics
    pfconvergence2$gelman

    ## Potential scale reduction factors:
    ##
    ##                            Point est. Upper C.I.
    ## alpha.(Intercept)                1.01       1.05
    ## alpha.log(IVI)                   1.02       1.07
    ## alpha.broodsize                  1.01       1.06
    ## alpha.sex                        1.01       1.02
    ## psi.(Intercept)                  1.02       1.08
    ## sd.epsilon.indidx                1.21       1.49
    ## sd.epsilon.indidx:log(IVI)       1.23       1.72

    ## Raftery diagnostics
    pfconvergence2$raftery

    ##
    ## Quantile (q) = 0.025
    ## Accuracy (r) = +/- 0.01
    ## Probability (s) = 0.95
    ##
    ## You need a sample size of at least 937 with these values of q, r and s

    ## Effective sample size
    pfconvergence2$effectiveSize

    ##          alpha.(Intercept)             alpha.log(IVI)
    ##                  123.37170                  247.82614
    ##            alpha.broodsize                  alpha.sex
    ##                  111.13724                  105.82065
    ##            psi.(Intercept)          sd.epsilon.indidx
    ##                  300.00000                  175.63886
    ## sd.epsilon.indidx:log(IVI)
    ##                   50.61415

2. Traceplots

    ## Generate traceplots
    pftraceplots2 <- traceplots(pfresults3,plot = FALSE)

    ## Fixed effects for mean
    pftraceplots2$meanFixed

    ## Fixed effects for variance
    pftraceplots2$varianceFixed

    ## Random effects variances for mean
    pftraceplots2$meanRandom

    ## Random effects variances for variances
    pftraceplots2$varianceRandom

    ## NULL

3. Numerical summaries

    ## Compute numerical summaries
    summary(pfresults3)

    ##
    ## Iterations = 1001:1991
    ## Thinning interval = 10
    ## Number of chains = 3
    ## Sample size per chain = 100
    ##
    ## Posterior Summary Statistics for Each Model Component
    ##
    ## Mean Model: Fixed Effects
    ##              Mean Median   SD Lower 95% Lower 50% Upper 50% Upper 95%
    ## (Intercept) -3.12  -3.12 0.17     -3.44     -3.18     -2.96     -2.78
    ## log(IVI)     0.12   0.12 0.01      0.09      0.11      0.12      0.14
    ## broodsize    0.12   0.12 0.07     -0.03      0.08      0.17      0.26
    ## sex         -0.08  -0.08 0.07     -0.24     -0.12     -0.02      0.05
    ##
    ## Mean Model: Random Effects
    ##                 Mean Median   SD Lower 95% Lower 50% Upper 50% Upper 95%
    ## indidx          0.22   0.21 0.06      0.14      0.20      0.24      0.29
    ## indidx:log(IVI) 0.02   0.02 0.01      0.00      0.02      0.03      0.04
    ##
    ## Variance Model: Fixed Effects
    ##              Mean Median   SD Lower 95% Lower 50% Upper 50% Upper 95%
    ## (Intercept) -0.37  -0.37 0.01      -0.4     -0.38     -0.37     -0.34

4. Graphical summaries

    ## Generate caterpillar
    pfcaterpillar2 <- caterpillar(pfresults3,plot = FALSE)

    ## Fixed effects for mean
    pfcaterpillar2$meanFixed

    ## Fixed effects for variance
    pfcaterpillar2$varianceFixed

The convergence diagnostics and traceplots suggest that the chains should be run for longer. In particular, the effective sample size for the variance of the random slopes is very small, but we will proceed with inference for illustration. In the numerical summaries we see that the estimated variance of the random slopes is very small, .02 with 95% credible interval (.00,.04). This suggests that there is very little inter-individual variation in the effect of `log(IVI)`.

Estimates of the individual random effects (predictions for frequentists) can be obtained from the function `ranef()`. This function computes posterior summary statistics for the individual random effects. In this case, the random effects for the mean include both the intercept and the effect of `log(IVI)`.

The following code computes the summary statistics and plots the predicted values of the random slopes with HPD 95% confidence intervals.

    ## Compute summary statistics for random effects
    ranef2 <- ranef(pfresults3)

    ## Load ggplot2
    library(ggplot2)

    ## Identify number of individuals
    nind <- nlevels(pfdata$indidx)

    ## Plot predicted random slopes
    ggplot(data=as.data.frame(ranef2$mean[nind+(1:nind),]),aes(x=1:nind,y=Mean)) +
      geom_point() +
      geom_errorbar(aes(ymin=`Lower 95%`,ymax=`Upper 95%`)) +
      geom_abline(intercept=0,slope=0)

Note that the 95% crebile intervals for all of the values include 0. This is further evidence that the effect of `log(IVI)` does not vary significantly between the individuals.