# 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 best done via visual display of the results, separately for each parameters of interest.

### Usage

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | ```
dctable(x, ...)
## Default S3 method:
dctable(x, ...)
## S3 method for class 'dctable'
plot(x, which = 1:length(x),
type = c("all", "var", "log.var"),
position = "topright", box.cex = 0.75, box.bg, ...)
extractdctable(x, ...)
## Default S3 method:
extractdctable(x, ...)
dcdiag(x, ...)
## Default S3 method:
dcdiag(x, ...)
## S3 method for class 'dcdiag'
plot(x, which = c("all", "lambda.max",
"ms.error", "r.squared", "log.lambda.max"),
position = "topright", ...)
extractdcdiag(x, ...)
## Default S3 method:
extractdcdiag(x, ...)
``` |

### Arguments

`x` |
An MCMC or a 'dctable' object. |

`...` |
Optionally more fitted model objects for function |

`which` |
What to plot. For |

`type` |
Type of plot to be drawn. See Details. |

`position` |
Position for the legend, as for |

`box.cex` |
Scaling factor for the interquartile boxes. |

`box.bg` |
Background color 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 for `dctable`

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"`

results 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 = "log.var"`

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.

Lele et al. (2010) introduced diagnostic measures
for checking the convergence of the data cloning algorithm
which are based on the joint posterior distribution
and not only on single parameters. These
include to calculate the largest eigenvalue of the posterior
variance covariance matrix (`lambda.max`

as returned by
`lambdamax.diag`

),
or to calculate mean squared error (`ms.error`

) and another
correlation-like fit statistic (`r.squared`

) based on a
Chi-squared approximation
(as returned by `chisq.diag`

). The maximum
eigenvalue reflects the degenerateness of the
posterior distribution, while the two fit measures reflect
if the Normal approximation is adequate. All
three statistics should converge to zero as the number of clones
increases. If this happens, different prior specifications are no
longer influencing the results (Lele et al., 2007, 2010).
These are conveniently collected by the `dcdiag`

function.

**IMPORTANT!**
Have you checked if different prior specifications
lead to the same results?

### 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.

`dcdiag`

returns a data frame with convergence diagnostics.

The `plot`

methods produce graphs as side effect.

### Author(s)

Peter Solymos, 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.

Lele, S. R., K. Nadeem and B. Schmuland, 2010.
Estimability and likelihood inference for generalized
linear mixed models using data cloning.
*Journal of the American Statistical Association*
**105**, 1617–1625.

Solymos, P., 2010. dclone: Data Cloning in R.
*The R Journal* **2(2)**, 29–37.
URL: http://journal.r-project.org/archive/2010-2/RJournal_2010-2_Solymos.pdf

### See Also

Data cloning: `dclone`

Model fitting: `jags.fit`

, `bugs.fit`

,
`dc.fit`

### Examples

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 | ```
## 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 <- crossprod(t(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)
``` |