# Expected deviance Model selection criterion based on posterior predictive distributions and approximations of the expected deviance. ```r expected.deviance(target, postsumstat, kernel = "gaussian", subset=NULL, print=TRUE) ``` ## Arguments - `target`: a vector of the observed summary statistics. - `postsumstat`: a vector, matrix or data frame of summary statistics simulated a posteriori. - `kernel`: a character string specifying the kernel to be used when. Defaults to `"gaussian"`. See `density` for details. - `subset`: a logical expression indicating elements or rows to keep. Missing values in `postsumstat` are taken as `FALSE`. - `print`: prints out what percent of the distances have been zero. ## Details This function implements an approximation for the expected deviance based on simulation performed a posteriori. Thus, after the posterior distribution of parameters or the posterior model probabilities have been determined, users need to re-simulate data using the posterior. The Monte-Carlo estimate of the expected deviance is computed from the simulated data as follows: c("$D=-\\frac{2}{n}\\sum_{j=1}^{n}\\log(K_\\epsilon(\\parallel\n$", "$ s^j-s_0\\parallel))$"), where n is number of simulations, $K$ is the statistical kernel, $\epsilon$ is the error, i.e. difference between the observed and simulated summary statistics below which simualtions were accepted in the original call to `postpr`, the $s^j$'s are the summary statistics obtained from the posterior predictive simualtions, and $s_0$ are the observed values of the summary statistics. The expected devaince averaged over the posterior distribution to compute a deviance information criterion (DIC). ## References Francois O, Laval G (2011) Deviance information criteria for model selection in approximate Bayesian computation **arXiv**:0240377 . ## Returns A list with the following components: - **expected.deviance**: The approximate expected deviance. - **dist**: The Euclidean distances for summary statistics simulated a posteriori. ## Examples ```r ## Function definitions skewness <- function(x) { sk <- mean((x-mean(x))^3)/(sd(x)^3) return(sk) } kurtosis <- function(x) { k <- mean((x-mean(x))^4)/(sd(x)^4) - 3 return(k) } ## Observed summary statistics obs.sumstat <- c(2.004821, 3.110915, -0.7831861, 0.1440266) ## Model 1 (Gaussian) ## ################## ## Simulate data theta <- rnorm(10000, 2, 10) zeta <- 1/rexp(10000, 1) param <- cbind(theta, zeta) y <- matrix(rnorm(200000, rep(theta, each = 20), sd = rep(sqrt(zeta), each = 20)), nrow = 20, ncol = 10000) ## Calculate summary statistics s <- cbind(apply(y, 2, mean), apply(y, 2, sd), apply(y, 2, skewness), apply(y, 2, kurtosis)) ## ABC inference gaus <- abc(target=obs.sumstat, param = param, sumstat=s, tol=.1, hcorr = FALSE, method = "loclinear") param.post <- gaus$adj.values ## Posterior predictive simulations postpred.gaus <- matrix(rnorm(20000, rep(param.post[,1], each = 20), sd = rep(sqrt(param.post[,2]), each = 20)), nrow = 20, ncol = 1000) statpost.gaus <- cbind(apply(postpred.gaus, 2, mean),apply(postpred.gaus, 2, sd),apply(postpred.gaus, 2,skewness),apply(postpred.gaus, 2,kurtosis)) # Computation of the expected deviance expected.deviance(obs.sumstat, statpost.gaus)$expected.deviance expected.deviance(obs.sumstat, statpost.gaus, kernel = "epanechnikov")$expected.deviance ## Modele 2 (Laplace) ## ################## ## Simulate data zeta <- rexp(10000) param <- cbind(theta, zeta) y <- matrix(theta + sample(c(-1,1),200000, replace = TRUE)*rexp(200000, rep(zeta, each = 20)), nrow = 20, ncol = 10000) ## Calculate summary statistics s <- cbind( apply(y, 2, mean), apply(y, 2, sd), apply(y, 2, skewness), apply(y, 2, kurtosis)) ## ABC inference lapl <- abc(target=obs.sumstat, param = param, sumstat=s, tol=.1, hcorr = FALSE, method = "loclinear") param.post <- lapl$adj.values ## Posterior predictive simulations postpred.lapl <- matrix(param.post[,1] + sample(c(-1,1),20000, replace = TRUE)*rexp(20000, rep(param.post[,2], each = 20)), nrow = 20, ncol = 1000) statpost.lapl <- cbind(apply(postpred.lapl, 2, mean),apply(postpred.lapl, 2, sd),apply(postpred.lapl, 2,skewness),apply(postpred.lapl, 2,kurtosis)) ## Computation of the expected deviance expected.deviance(obs.sumstat, statpost.lapl)$expected.deviance expected.deviance(obs.sumstat, statpost.lapl, kernel = "epanechnikov")$expected.deviance ```

expected.deviance function