Nothing
Pied Flycatchers 1: Initial Models
In dalmatian: Automating the Fitting of Double Linear Mixed Models in 'JAGS' and 'nimble'
As a first example of the functions in this package we will fit a model to the pied flycatcher data. We will start by fitting a model with fixed effects on the mean and dispersion components. Specifically, we will model the load returned on each trip by the adult flycatchers a independent normal random variables with mean given by a linear combination of the logarithm of the inter-visit interval (IVI), the brood size manipulation (broodsize), and the adult's sex, and standard deviation given by a linear combination of broodsize and sex, on the log scale. We will then incorporate individual random effects into this model.
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)
This data contains records on 5795 different food deliveries recorded during the pied flycatcher experiment. The response variable, load, records the load rounded to the nearest .1 gram. We want to model the logarithm of the load, and we also want to account for the rounding. To do this we will create two new variables in the data frame, lower and upper, which bound the logarithm of the true load. Not that we add .001 when the observed value is zero to avoid problems with the logarithm.
## Create variables bounding the true load
pfdata$lower=ifelse(pfdata$load==0,log(.001),log(pfdata$load-.049))
pfdata$upper=log(pfdata$load+.05)
Model 1: Fixed Effects
Our initial model for the logarithm of the load returned on the $j$th trip by adult $i$ will be $N(\mu_{ij},\sigma_{ij}^2)$ where:
$$\mu=\beta_0 + \beta_1 \mathrm{log(IVI)}{ij} + \beta_2 \mathrm{broodsize}_i + \beta_3 \mathrm{sex}_i$$
and
$$\log(\sigma{ij}^2)=\psi_0 + \psi_1 \mathrm{broodsize}i + \psi_2\mathrm{sex}_i.$$
In the framework of generalized linear models, the variance, $\sigma^2_i$, is equal to the dispersion parameter, which we denote generically with $\phi$. Hence, the model can then be written as:
$Y{ij} \sim N(\mu_{ij},\phi_{ij})$ where:
$$\mu=\beta_0 + \beta_1 \mathrm{log(IVI)}{ij} + \beta_2 \mathrm{broodsize}_i + \beta_3 \mathrm{sex}_i$$
and
$$\log(\phi{ij})=\psi_0 + \psi_1 \mathrm{broodsize}_i + \psi_2\mathrm{sex}_i.$$
Model
The first step is to define the models for the mean and dispersion components. These models are simply lists with (up to) two named components, fixed and random, which provide the details on the fixed and random effects. Both lists must contain a variable name specifying the basename for the coefficients and a variable formula specifying the effects. The fixed list should also contain an object called priors which specifies the priors for coefficients in the fixed effects portion of the model. Random effects are currently assumed to be normal with mean zero and unknown variances which are assigned the half $t$-distribution with ? degrees of freedom and scale ?. The optional parameter link can also be used to specify a link function for either component.
The components of the model for the pied flycatcher data specified above would be generated with:
## Mean model
mymean <- list(fixed=list(name="alpha",
formula=~ log(IVI) + broodsize + sex,
priors=list(c("dnorm",0,.001))))
## Dispersion model
mydisp <- list(fixed=list(name="psi",
formula=~broodsize + sex,
priors=list(c("dnorm",0,.001))),
link="log")
These two objects will now be used to generate the JAGS code, data, and initial values for running the model.
Running the Model with dalmatian
The primary function in the package dalmatian automates the creation of the JAGS code, data, and initial values and then passes these to JAGS via the functions in the rjags package. The dalmatian function takes requires several arguments including the data frame for the analysis, the model objects created above, and the name of the JAGS model file. It also requires two lists containing named arguments for the functions jags.model and coda.samples from the rjags package. Descriptions of the arguments for these two functions are available in their own help pages. Any arguments that are not specified in these lists will take the default values given by rjags. The two exceptions are the file argument of jags.model and the n.iter argument of coda.samples which do not have default values and must be specified. The number of parallel chains will be taken from the jags.model list. If this value is not specified then three chains will be run in order that convergence diagnostics can be computed. 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 \texttt{CRAN}. 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_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
pfmcmc1 <- dalmatian(df=pfdata,
mean.model=mymean,
dispersion.model=mydisp,
jags.model.args=jm.args,
coda.samples.args=cs.args,
rounding=TRUE,
lower="lower",
upper="upper",
residuals = FALSE,
debug=FALSE)
}
load(system.file("pfresults1.RData",package="dalmatian"))
Results
The function dalmatian returns list containing three objects:
1. jags.model.args
The full list of arguments passed to jags.model. This includes the variables in the input list along with the name of the model file, the formatted data, and the initial values.
2. coda.samples.args
The full list of arguments passed to coda.samples. This includes the variables in the input list along with the complete object returned by jags.model.args and the character vector including the names of the monitored parameters.
3. coda
An object of class mcmc.list and length n.chains containing the samples generated by JAGS.
Results
Inference is based on the output from the MCMC chain which is stored in the coda object. Summaries of the MCMC sampler and the posterior distribution can be constructed from coda with a variety of tools designed to work with objects of the mcmc.list class. For example, we can use the tools in the package coda to generate traceplots, density plots and histograms, numerical summaries of the posterior, and convergence diagnostics. The ggmcmc package can also be used to create fancier looking traceplots, density plots, and plots summarizing the posterior distribution in terms of posterior means and quantiles (credible intervals). However, several wrapper functions exist within the package to automate some of these tasks.
The first thing we should do is to check the convergence of the chains. Summary statistics computed from the output in coda are only valid if the chain has converged and is sampling from the correct distribution. The wrapper function convergence() applies three of the convergence diagnostics from the coda package.
## Compute convergence diagnostics
pfconvergence1 <- convergence(pfmcmc1)
## Gelman-Rubin diagnostics
pfconvergence1$gelman
## Raftery diagnostics
pfconvergence1$raftery
## Effective sample size
pfconvergence1$effectiveSize
Details of these measures can be found in the help files of the coda package. Not surprisingly, these all suggest that the chains need to be run for longer.
Convergence and mixing of the chains can also be examined visually through traceplots. The wrapper function traceplots generates traces using the functions from the ggmcmc package and separates the variables by the model components. One of the advantages of using ggmcmc which is based on ggplot2 is that the plots can be saved and easily reproduced. Here I have thinned the chains further prior to plotting to reduce the size of the vignette. Note that thinning is conducted relative to the original chain, not the already thinned chain. I.e., thinning by 20 selects the output from every 20th iteration of the original chain (every 2nd iteration from the retained samples).
## Generate traceplots
pfraceplots1 <- traceplots(pfmcmc1,show=FALSE,nthin=20)
## Fixed effects for mean
pftraceplots1$meanFixed
## Fixed effects for dispersion
pftraceplots1$dispersionFixed
Numerical summaries of the posterior distribution computed from the output in coda can be obtained with the summary function.
## Compute numerical summaries
pfsummary1 <- summary(pfmcmc1)
## Print numerical summaries
pfsummary1
For each component of the model the tables provide statistics including the estimated posterior mean, median, and standard deviation and the end points of the 50% and 95% highest posterior density credible intervals. Note that these summaries are only valid if the previous diagnostics have suggested that the chain has converged and the effective sample size is large enough. Graphical summaries can also be computed with the caterpillar which, like traceplots, is a simple wrapper for the functions in the ggmcmc package.
## Generate caterpillar
pfcaterpillar1 <- caterpillar(pfmcmc1,show = FALSE)
## Fixed effects for mean
pfcaterpillar1$meanFixed
## Fixed effects for dispersion
pfcaterpillar1$dispersionFixed
Model 2: Random Effects
We will now consider adding individual (bird) random effects to both the fixed and random components of our model. This can be done by simply adding to the model components created above and then re-runnig the model. The new model will have mean and variance
$$\mu_{ij}=\beta_0 + \beta_1 \mathrm{log(IVI)}{ij} + \beta_2 \mathrm{broodsize}_i + \beta_3 \mathrm{sex}_i + \epsilon{i}}$$
and
$$\log(\phi_{ij})=\psi_0 + \psi_1 \mathrm{broodsize}i + \psi_2\mathrm{sex}_i + \xi_j$$
where $\epsilon_i \sim N(0,\tau^2\mu)$ and $\xi_i \sim N(0,\tau^2_\phi)$. By default the variances $\tau^2_\epsilon$ and $\tau^2_\xi$ are assigned half-$t$ prior distributions.
Model
Instead of creating new model objects we can simply add to the old objects by creating the random fields. As with the fixed effects, these fields are lists which must include a basename for the random effects and formula specifying the random effects model.
# Random component of mean
mymean$random <- list(name="epsilon",formula=~-1 + indidx)
# Random component of dispersion
mydisp$random <- list(name="xi",formula=~-1 + indidx)
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 dispersion components by taking the values from the final iterations of each of the chains run for the fixed effects model. I have also found that it is best to set initial values for the response variable, y, if the response is rounded. Otherwise, JAGS often fails to generate initial values that fall between the rounding limits. 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 \texttt{CRAN}. Neither of these is recommended in practice.
## Define initial values
inits <- lapply(1:3,function(i){
setJAGSInits(mymean,
mydisp,
y = runif(nrow(pfdata),pfdata$lower,pfdata$upper),
fixed.mean = tail(pfmcmc1$coda[[i]],1)[1:4],
fixed.dispersion = tail(pfmcmc1$coda[[i]],1)[5:7],
sd.mean = 1,
sd.dispersion=1)
})
## Define list of arguments for jags.model()
jm.args <- list(file=file.path(workingDir,"pied_flycatcher_2_jags.R"),inits=inits,n.adapt=1000)
## Define list of arguments for coda.samples()
cs.args <- list(n.iter=5000,thin=10)
## Run the model using dalmatian
pfmcmc2 <- dalmatian(df=pfdata,
mean.model=mymean,
dispersion.model=mydisp,
jags.model.args=jm.args,
coda.samples.args=cs.args,
rounding=TRUE,
lower="lower",
upper="upper",
residuals=FALSE,
debug=FALSE)
load(system.file("pfresults2.RData",package="dalmatian"))
Results
Once again, we can use the wrapper functions within the package to examine the convergence and compute summaries of the posterior distribution. Note that the chains have been thinned in creating the traceplots to reduce the size of the package.
- Convergence diagnostics
## Compute convergence diagnostics
pfconvergence2 <- convergence(pfmcmc2)
## Gelman-Rubin diagnostics
pfconvergence2$gelman
## Raftery diagnostics
pfconvergence2$raftery
## Effective sample size
pfconvergence2$effectiveSize
- Traceplots
## Generate traceplots
pftraceplots2 <- traceplots(pfmcmc2, show = FALSE, nthin=20)
## Fixed effects for mean
pftraceplots2$meanFixed
## Fixed effects for dispersion
pftraceplots2$dispersionFixed
## Random effects variances for mean
pftraceplots2$meanRandom
## Random effects variances for dispersion
pftraceplots2$dispersionRandom
- Numerical summaries
## Compute numerical summaries
pfsummary2 <- summary(pfmcmc2)
## Print numerical summaries
pfsummary2
- Graphical summaries
## Generate caterpillar
pfcaterpillar2 <- caterpillar(pfmcmc2,show = FALSE)
## Fixed effects for mean
pfcaterpillar2$meanFixed
## Fixed effects for dispersion
pfcaterpillar2$dispersionFixed
Predictions
By using 'prediction' function, predictions for the mean and dispersion models can be calculated. The 'prediction' function should take a data frame including all variables required and their names defined for the usage of 'dalmatian' function. For the Pied Flycatcher data, previously 'log(IVI)' was used to define the mean and dispersion models. Hence, it should be included in the data frame as following:
## Add 'log(IVI)' variable in pfdata
pfdata$'log(IVI)' <- log(pfdata$IVI)
As an example we computed the fitted means and dispersion parameters from both models for the first 5 observations.
## Compute predictions for model 1
pffitted1 <- predict(pfmcmc,
newdata = pfdata[1:5,])
pffitted1$mean
pffitted1$dispersion
## Compute predictions for model 2
pffitted2 <- predict(pfmcmc2,
newdata = pfdata[1:5,]))
pffitted2$mean
pffitted2$dispersion
Try the dalmatian package in your browser
Any scripts or data that you put into this service are public.
dalmatian documentation built on Nov. 23, 2021, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
As a first example of the functions in this package we will fit a model to the pied flycatcher data. We will start by fitting a model with fixed effects on the mean and dispersion components. Specifically, we will model the load returned on each trip by the adult flycatchers a independent normal random variables with mean given by a linear combination of the logarithm of the inter-visit interval (IVI), the brood size manipulation (broodsize), and the adult's sex, and standard deviation given by a linear combination of broodsize and sex, on the log scale. We will then incorporate individual random effects into this model.
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)
This data contains records on 5795 different food deliveries recorded during the pied flycatcher experiment. The response variable, load, records the load rounded to the nearest .1 gram. We want to model the logarithm of the load, and we also want to account for the rounding. To do this we will create two new variables in the data frame, lower and upper, which bound the logarithm of the true load. Not that we add .001 when the observed value is zero to avoid problems with the logarithm.
## Create variables bounding the true load pfdata$lower=ifelse(pfdata$load==0,log(.001),log(pfdata$load-.049)) pfdata$upper=log(pfdata$load+.05)
Model 1: Fixed Effects
Our initial model for the logarithm of the load returned on the $j$th trip by adult $i$ will be $N(\mu_{ij},\sigma_{ij}^2)$ where: $$\mu=\beta_0 + \beta_1 \mathrm{log(IVI)}{ij} + \beta_2 \mathrm{broodsize}_i + \beta_3 \mathrm{sex}_i$$ and $$\log(\sigma{ij}^2)=\psi_0 + \psi_1 \mathrm{broodsize}i + \psi_2\mathrm{sex}_i.$$ In the framework of generalized linear models, the variance, $\sigma^2_i$, is equal to the dispersion parameter, which we denote generically with $\phi$. Hence, the model can then be written as: $Y{ij} \sim N(\mu_{ij},\phi_{ij})$ where: $$\mu=\beta_0 + \beta_1 \mathrm{log(IVI)}{ij} + \beta_2 \mathrm{broodsize}_i + \beta_3 \mathrm{sex}_i$$ and $$\log(\phi{ij})=\psi_0 + \psi_1 \mathrm{broodsize}_i + \psi_2\mathrm{sex}_i.$$
Model
The first step is to define the models for the mean and dispersion components. These models are simply lists with (up to) two named components, fixed and random, which provide the details on the fixed and random effects. Both lists must contain a variable name specifying the basename for the coefficients and a variable formula specifying the effects. The fixed list should also contain an object called priors which specifies the priors for coefficients in the fixed effects portion of the model. Random effects are currently assumed to be normal with mean zero and unknown variances which are assigned the half $t$-distribution with ? degrees of freedom and scale ?. The optional parameter link can also be used to specify a link function for either component.
The components of the model for the pied flycatcher data specified above would be generated with:
## Mean model mymean <- list(fixed=list(name="alpha", formula=~ log(IVI) + broodsize + sex, priors=list(c("dnorm",0,.001)))) ## Dispersion model mydisp <- list(fixed=list(name="psi", formula=~broodsize + sex, priors=list(c("dnorm",0,.001))), link="log")
These two objects will now be used to generate the JAGS code, data, and initial values for running the model.
Running the Model with dalmatian
The primary function in the package dalmatian automates the creation of the JAGS code, data, and initial values and then passes these to JAGS via the functions in the rjags package. The dalmatian function takes requires several arguments including the data frame for the analysis, the model objects created above, and the name of the JAGS model file. It also requires two lists containing named arguments for the functions jags.model and coda.samples from the rjags package. Descriptions of the arguments for these two functions are available in their own help pages. Any arguments that are not specified in these lists will take the default values given by rjags. The two exceptions are the file argument of jags.model and the n.iter argument of coda.samples which do not have default values and must be specified. The number of parallel chains will be taken from the jags.model list. If this value is not specified then three chains will be run in order that convergence diagnostics can be computed. 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 \texttt{CRAN}. 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_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 pfmcmc1 <- dalmatian(df=pfdata, mean.model=mymean, dispersion.model=mydisp, jags.model.args=jm.args, coda.samples.args=cs.args, rounding=TRUE, lower="lower", upper="upper", residuals = FALSE, debug=FALSE) }
load(system.file("pfresults1.RData",package="dalmatian"))
Results
The function dalmatian returns list containing three objects:
1. jags.model.args
The full list of arguments passed to jags.model. This includes the variables in the input list along with the name of the model file, the formatted data, and the initial values.
2. coda.samples.args
The full list of arguments passed to coda.samples. This includes the variables in the input list along with the complete object returned by jags.model.args and the character vector including the names of the monitored parameters.
3. coda
An object of class mcmc.list and length n.chains containing the samples generated by JAGS.
Results
Inference is based on the output from the MCMC chain which is stored in the coda object. Summaries of the MCMC sampler and the posterior distribution can be constructed from coda with a variety of tools designed to work with objects of the mcmc.list class. For example, we can use the tools in the package coda to generate traceplots, density plots and histograms, numerical summaries of the posterior, and convergence diagnostics. The ggmcmc package can also be used to create fancier looking traceplots, density plots, and plots summarizing the posterior distribution in terms of posterior means and quantiles (credible intervals). However, several wrapper functions exist within the package to automate some of these tasks.
The first thing we should do is to check the convergence of the chains. Summary statistics computed from the output in coda are only valid if the chain has converged and is sampling from the correct distribution. The wrapper function convergence() applies three of the convergence diagnostics from the coda package.
## Compute convergence diagnostics pfconvergence1 <- convergence(pfmcmc1)
## Gelman-Rubin diagnostics pfconvergence1$gelman ## Raftery diagnostics pfconvergence1$raftery ## Effective sample size pfconvergence1$effectiveSize
Details of these measures can be found in the help files of the coda package. Not surprisingly, these all suggest that the chains need to be run for longer.
Convergence and mixing of the chains can also be examined visually through traceplots. The wrapper function traceplots generates traces using the functions from the ggmcmc package and separates the variables by the model components. One of the advantages of using ggmcmc which is based on ggplot2 is that the plots can be saved and easily reproduced. Here I have thinned the chains further prior to plotting to reduce the size of the vignette. Note that thinning is conducted relative to the original chain, not the already thinned chain. I.e., thinning by 20 selects the output from every 20th iteration of the original chain (every 2nd iteration from the retained samples).
## Generate traceplots pfraceplots1 <- traceplots(pfmcmc1,show=FALSE,nthin=20)
## Fixed effects for mean pftraceplots1$meanFixed ## Fixed effects for dispersion pftraceplots1$dispersionFixed
Numerical summaries of the posterior distribution computed from the output in coda can be obtained with the summary function.
## Compute numerical summaries pfsummary1 <- summary(pfmcmc1)
## Print numerical summaries
pfsummary1
For each component of the model the tables provide statistics including the estimated posterior mean, median, and standard deviation and the end points of the 50% and 95% highest posterior density credible intervals. Note that these summaries are only valid if the previous diagnostics have suggested that the chain has converged and the effective sample size is large enough. Graphical summaries can also be computed with the caterpillar which, like traceplots, is a simple wrapper for the functions in the ggmcmc package.
## Generate caterpillar pfcaterpillar1 <- caterpillar(pfmcmc1,show = FALSE)
## Fixed effects for mean pfcaterpillar1$meanFixed ## Fixed effects for dispersion pfcaterpillar1$dispersionFixed
Model 2: Random Effects
We will now consider adding individual (bird) random effects to both the fixed and random components of our model. This can be done by simply adding to the model components created above and then re-runnig the model. The new model will have mean and variance $$\mu_{ij}=\beta_0 + \beta_1 \mathrm{log(IVI)}{ij} + \beta_2 \mathrm{broodsize}_i + \beta_3 \mathrm{sex}_i + \epsilon{i}}$$ and $$\log(\phi_{ij})=\psi_0 + \psi_1 \mathrm{broodsize}i + \psi_2\mathrm{sex}_i + \xi_j$$ where $\epsilon_i \sim N(0,\tau^2\mu)$ and $\xi_i \sim N(0,\tau^2_\phi)$. By default the variances $\tau^2_\epsilon$ and $\tau^2_\xi$ are assigned half-$t$ prior distributions.
Model
Instead of creating new model objects we can simply add to the old objects by creating the random fields. As with the fixed effects, these fields are lists which must include a basename for the random effects and formula specifying the random effects model.
# Random component of mean mymean$random <- list(name="epsilon",formula=~-1 + indidx) # Random component of dispersion mydisp$random <- list(name="xi",formula=~-1 + indidx)
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 dispersion components by taking the values from the final iterations of each of the chains run for the fixed effects model. I have also found that it is best to set initial values for the response variable, y, if the response is rounded. Otherwise, JAGS often fails to generate initial values that fall between the rounding limits. 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 \texttt{CRAN}. Neither of these is recommended in practice.
## Define initial values inits <- lapply(1:3,function(i){ setJAGSInits(mymean, mydisp, y = runif(nrow(pfdata),pfdata$lower,pfdata$upper), fixed.mean = tail(pfmcmc1$coda[[i]],1)[1:4], fixed.dispersion = tail(pfmcmc1$coda[[i]],1)[5:7], sd.mean = 1, sd.dispersion=1) }) ## Define list of arguments for jags.model() jm.args <- list(file=file.path(workingDir,"pied_flycatcher_2_jags.R"),inits=inits,n.adapt=1000) ## Define list of arguments for coda.samples() cs.args <- list(n.iter=5000,thin=10) ## Run the model using dalmatian pfmcmc2 <- dalmatian(df=pfdata, mean.model=mymean, dispersion.model=mydisp, jags.model.args=jm.args, coda.samples.args=cs.args, rounding=TRUE, lower="lower", upper="upper", residuals=FALSE, debug=FALSE)
load(system.file("pfresults2.RData",package="dalmatian"))
Results
Once again, we can use the wrapper functions within the package to examine the convergence and compute summaries of the posterior distribution. Note that the chains have been thinned in creating the traceplots to reduce the size of the package.
- Convergence diagnostics
## Compute convergence diagnostics pfconvergence2 <- convergence(pfmcmc2)
## Gelman-Rubin diagnostics pfconvergence2$gelman ## Raftery diagnostics pfconvergence2$raftery ## Effective sample size pfconvergence2$effectiveSize
- Traceplots
## Generate traceplots pftraceplots2 <- traceplots(pfmcmc2, show = FALSE, nthin=20)
## Fixed effects for mean pftraceplots2$meanFixed ## Fixed effects for dispersion pftraceplots2$dispersionFixed ## Random effects variances for mean pftraceplots2$meanRandom ## Random effects variances for dispersion pftraceplots2$dispersionRandom
- Numerical summaries
## Compute numerical summaries pfsummary2 <- summary(pfmcmc2)
## Print numerical summaries
pfsummary2
- Graphical summaries
## Generate caterpillar pfcaterpillar2 <- caterpillar(pfmcmc2,show = FALSE)
## Fixed effects for mean pfcaterpillar2$meanFixed ## Fixed effects for dispersion pfcaterpillar2$dispersionFixed
Predictions
By using 'prediction' function, predictions for the mean and dispersion models can be calculated. The 'prediction' function should take a data frame including all variables required and their names defined for the usage of 'dalmatian' function. For the Pied Flycatcher data, previously 'log(IVI)' was used to define the mean and dispersion models. Hence, it should be included in the data frame as following:
## Add 'log(IVI)' variable in pfdata pfdata$'log(IVI)' <- log(pfdata$IVI)
As an example we computed the fitted means and dispersion parameters from both models for the first 5 observations.
## Compute predictions for model 1 pffitted1 <- predict(pfmcmc, newdata = pfdata[1:5,])
pffitted1$mean pffitted1$dispersion
## Compute predictions for model 2 pffitted2 <- predict(pfmcmc2, newdata = pfdata[1:5,]))
pffitted2$mean pffitted2$dispersion
Try the dalmatian package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.