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.
First you have to load the package:
## Load package library(dalmatian)
## 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)
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:
## 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.
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)
.
## 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)
load(system.file("pfresults4.RData",package="dalmatian"))
The following illustrates the functions for processing the results. Again, these are fully described in vignette(pied-flycatchers-1)
.
## Compute convergence diagnostics pfconvergence4 <- convergence(pfmcmc4)
## Gelman-Rubin diagnostics pfconvergence4$gelman ## Raftery diagnostics pfconvergence4$raftery ## Effective sample size pfconvergence4$effectiveSize
## Generate traceplots pftraceplots4 <- traceplots(pfmcmc4,show=FALSE,nthin=10)
## Fixed effects for mean pftraceplots4$meanFixed ## Fixed effects for dispersion pftraceplots4$dispersionFixed
## Compute numerical summaries pfsummary4 <- summary(pfmcmc4)
## Print numerical summaries
pfsummary4
## Generate caterpillar pfcaterpillar4 <- caterpillar(pfmcmc4,show = FALSE)
## Fixed effects for mean pfcaterpillar4$meanFixed ## Fixed effects for dispersion pfcaterpillar4$dispersionFixed
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.