R/calculateDIC.R In polySegratioMM: Bayesian mixture models for marker dosage in autopolyploids

calculateDIC <-
function(mcmc.mixture, model, priors, seg.ratios, chain=1, print.DIC=FALSE)
{

## calculate DIC_4 from Celeax et al

if (class(mcmc.mixture) != "segratioMCMC")
stop("'mcmc.mixture' must be of class 'segratioMCMC'")
if (class(model) != "modelSegratioMM")
stop("'model' must be of class 'modelSegratioMM'")
if (class(priors) != "priorsSegratioMM")
stop("'priors' must be of class 'priorsSegratioMM'")
if (class(seg.ratios) != "segRatio") {
stop("'seg.ratios' must be of class 'segRatio'")
}

mcmc <- mcmc.mixture$mcmc.list[[chain]] y <- gtools::logit(seg.ratios$seg.ratio)
K <- model$n.components n <- length(y) m <- dim(mcmc)[1] p.ind <- grep("P\\[", colnames(mcmc)) # 1,2,3,4 z.ind <- grep("T\\[", colnames(mcmc)) # 5 ...295 mu.ind <- grep("mu\\[", colnames(mcmc)) # 296 297 298 299 sigma.ind <- grep("sigma\\[", colnames(mcmc)) # 300 301 302 303 equal.variances <- model$equal.variances   # modify if equal vars
if (equal.variances)  {
sigma.ind <- grep("sigma", colnames(mcmc))
}

## need to think about this but it should be straight forward - loops
## to avoid massive memory

## I really need to learn to vectorise more of these calcs!

first <- 0
for (l in 1:m){ # l <- 1
tij <- mcmc[l,p.ind]

ynorm <- NULL
if (equal.variances){
for (j in 1:K){
ynorm <- cbind(ynorm, pnorm(y, mcmc[l,mu.ind[j]],mcmc[l,sigma.ind]))
}
} else {
for (j in 1:K){
ynorm <- cbind(ynorm, pnorm(y, mcmc[l,mu.ind[j]],mcmc[l,sigma.ind[j]]))
}
}
###   for (j in 1:K){
###      ynorm <- cbind(ynorm, pnorm(y, mcmc[l,mu.ind[j]],mcmc[l,sigma.ind[j]]))
###    }

tij <- mcmc[l,p.ind]*ynorm
tij <- tij/rowSums(tij)
## dim(mcmc[m,p.ind] %o% rep(1,n))
first <- first + sum(tij* log( rep(1,n) %o% mcmc[l,p.ind] * ynorm),
na.rm=TRUE)
}

mu <- mcmc[,mu.ind]   # set up the original parameterisation
for (i in length(mu.ind):2) {
mu[,i] <- mu[,i] -  mu[,i-1]
}

## need priors for second sum on page 13 Celeaux et al (2006)
mu.0 <- priors$params$logit.means
var.0 <- 1/priors$params$logit.prec
nu.0 <- 2*priors$params$A
var.nu.0 <- 2* priors$params$B
if (length(priors$params$dirichlet) == 0)
alpha <- rep(1, length(p.ind))

n.j <- rep(1,K)     # imagine that it should be 1 given specified priors

mjl <- t(apply(mcmc[,z.ind], 1, table))
p.jl <- (mjl+1)/(K+n)            # \bar{p}_j^(l)

muhat.jl <- NA*mcmc[,p.ind]# calculate \hat{\mu_j^(l) and \bar{\mu_j^(l)
s2hat.jl <- muhat.jl   # and also \hat{\sigma_j^2(l) and \bar{\sigma_j^2(l)
for (i in 1:m){
muhat.jl[i,] <- tapply(y, mcmc[i,z.ind], mean)
for (j in 1:K){ # j <- 1
s2hat.jl[i,j] <- sum((y[mcmc[i,z.ind]==j] - muhat.jl[i,j])^2)
}
}
mu.jl <- ( rep(1,m) %o% (n.j*mu.0) + muhat.jl*mjl )/( rep(1,m) %o% n.j + mjl)
if (equal.variances){
s2.jl <- ( var.nu.0 + s2hat.jl  +
(rep(1,m) %o% n.j * mjl)/(rep(1,m) %o% n.j + mjl)*
(muhat.jl - rep(1,m)%o%mu.0)^2 /(nu.0 + mjl -2) )
} else {
s2.jl <- ((rep(1,m) %o% var.nu.0) + s2hat.jl  +
(rep(1,m) %o% n.j * mjl)/(rep(1,m) %o% n.j + mjl)*
(muhat.jl - rep(1,m)%o%mu.0)^2 /(rep(1,m)%o%nu.0 + mjl -2) )
}

second <- 0
for (i in 1:m){
second <- second + sum(log(p.jl[i,mcmc[i,z.ind]]) +
pnorm(y, mu.jl[i,mcmc[i,z.ind]],
s2.jl[i,mcmc[i,z.ind]], log.p = TRUE),
na.rm=TRUE)
}

DIC <- - (4/m) * first + (2/m)*second

if (print.DIC)
cat("DIC:",DIC,"\n")

return(DIC)
}


Try the polySegratioMM package in your browser

Any scripts or data that you put into this service are public.

polySegratioMM documentation built on May 2, 2019, 4:41 p.m.