dctable {dclone}R Documentation

Retrieve descriptive statistics from fitted objects to evaluate convergence

Description

The function is used to retrieve descriptive statistics from fitted objects on order to evaluate convergence of the data cloning algorithm. This is done via visual display of the results, separately for each parameters of interest.

Usage

dctable(x, ...)
## Default S3 method:
dctable(x, ...)
## S3 method for class 'dctable':
plot(x, which = 1:length(x), 
type = c("all", "var", "logvar"),
position = "topright", box.cex = 0.75, ...)
extractdctable(x, ...)
## Default S3 method:
extractdctable(x, ...)

Arguments

x An MCMC or a 'dctable' object.
... Optionally more fitted model objects for function dctable.
which What to plot, character names or integer indices of the estimated parameters are accepted.
type Type of plot to be drawn. See Details.
position Position for the legend, as for legend.
box.cex Scaling factor for the interquartile boxes.

Details

dctable returns the "dctable" attribute of the MCMC object, or if it is NULL, calculates the dctable summaries. If more than one fitted objects are provided, summaries are calculated for all objects, and results are ordered by the number of clones.

The plot method helps in graphical representation of the descriptive statistics. type = "all" results in plotting means, standard deviations and quantiles against the number of clones as boxplot. type = "var" reults in plotting the scaled variances against the number of clones. In this case variances are divided by the variance of the model with smallest number of clones, min(n.clones). type = "logvar" is the same as "var", but on the log scale. Along with the values, the min(n.clones) / n.clones line is plotted for reference.

Value

An object of class 'dctable'. It is a list, and contains as many data frames as the number of parameters in the fitted object. Each data frame contains descriptives as the function of the number of clones.
The plot method produces graphs as side effect.

Author(s)

P\'eter S\'olymos, solymos@ualberta.ca, implementation is based on many discussions with Khurram Nadeem and Subhash Lele.

References

Lele, S.R., B. Dennis and F. Lutscher, 2007. Data cloning: easy maximum likelihood estimation for complex ecological models using Bayesian Markov chain Monte Carlo methods. Ecology Letters 10, 551–563.

See Also

Data cloning: dclone

Model fitting: jags.fit

Examples

## Not run: 
## simulation for Poisson GLMM
set.seed(1234)
n <- 20
beta <- c(2, -1)
sigma <- 0.1
alpha <- rnorm(n, 0, sigma)
x <- runif(n)
X <- model.matrix(~x)
linpred <- X %*% beta + alpha
Y <- rpois(n, exp(linpred))
## JAGS model as a function
jfun1 <- function() {
    for (i in 1:n) {
        Y[i] ~ dpois(lambda[i])
        log(lambda[i]) <- alpha[i] + inprod(X[i,], beta[1,])
        alpha[i] ~ dnorm(0, 1/sigma^2)
    }
    for (j in 1:np) {
        beta[1,j] ~ dnorm(0, 0.001)
    }
    sigma ~ dlnorm(0, 0.001)
}
## data
jdata <- list(n = n, Y = Y, X = X, np = NCOL(X))
## number of clones to be used, etc.
## iteartive fitting
jmod <- dc.fit(jdata, c("beta", "sigma"), jfun1, 
    n.clones = 1:5, multiply = "n", unchanged = "np")
## summary with DC SE and R hat
summary(jmod)
dct <- dctable(jmod)
plot(dct)
## How to use estimates to make priors more informative?
glmm.model.up <- function() {
    for (i in 1:n) {
        Y[i] ~ dpois(lambda[i])
        log(lambda[i]) <- alpha[i] + inprod(X[i,], beta[1,])
        alpha[i] ~ dnorm(0, 1/sigma^2)
    }
    for (j in 1:p) {
        beta[1,j] ~ dnorm(priors[j,1], priors[j,2])
    }
    sigma ~ dgamma(priors[(p+1),2], priors[(p+1),1])
}
## function for updating, x is an MCMC object
upfun <- function(x) {
    if (missing(x)) {
        p <- ncol(X)
        return(cbind(c(rep(0, p), 0.001), rep(0.001, p+1)))
    } else {
        par <- coef(x)
        return(cbind(par, rep(0.01, length(par))))
    }
}
updat <- list(n = n, Y = Y, X = X, p = ncol(X), priors = upfun())
dcmod <- dc.fit(updat, c("beta", "sigma"), glmm.model.up,
    n.clones = 1:5, multiply = "n", unchanged = "p",
    update = "priors", updatefun = upfun)
summary(dcmod)
dct <- dctable(dcmod)
plot(dct)
plot(dct, type = "var")
## End(Not run)

[Package dclone version 1.0-0 Index]