cqmc: Mode and/or MCMC sample of Cq-dPCR posterior

Description Usage Arguments Details Value Author(s) Examples

View source: R/cqmc.R

Description

Estimates the mode and produces an MCMC sample of the posterior distribution of the parameters of the Cq-dPCR model.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
cqmc(data, mc.rep = 10^4, h = NULL, 
  n0 = NULL, n1 = NULL, 
  nt0 = 0, nt1 = 0, n.lo = 0, n.hi = 0, 
  pars0 = NULL, report = 1, probreport = 0.1, 
  extra = c("trendx", "trendy", "disp", "E1"), 
  c0 = 6, maxn0 = 7, 
  mod.method = "Nelder-Mead", mod.rep = 2000, 
  burnin = 0, nreport = 10, cq.xlim = NA, tune = 1, 
  E1.init = 0.8, E.init = 0.9, d.init = 1, 
  prior = TRUE, 
  mu.fun = function(x) dgamma(x, 1.5, 1.5), 
  A.fun = function(x) x^-1, 
  E.fun = function(x) dbeta(x, 60, 5), 
  E1.fun = E.fun, 
  trendx.fun = function(x) 1, trendy.fun = function(x) 1, 
  disp.fun = function(x) dgamma(x, 10, 10))

Arguments

data

cqdat object or data frame

mc.rep

number of MCMC samples. If 0 then MCMC not performed.

h

Threshold value

n0, n1, nt0, nt1, n.lo, n.hi

integer counts of the negative partitions (n0), positive partitions (n1), number of trimmed negative partitions (nt0), trimmed positive partitions (nt1), low outliers (n.lo) and high outliers (n.hi). Not needed if data is cqdat object.

pars0

initial model parameter values

report

if postive tracing information is produced. High values may produce more information.

probreport

if report>0, the probability that "." is printed when a valid log-likelihood is calculated. Useful for showing progress in finding mode or simulating MCMC chain.

extra

vector of names of parameters to include in model in additon to "mu", "E", and "A". Defaults to all possibitilies: c("E1","trendx","trendy","disp").

c0

number of cycles for which exact probabilities are caclulated. This can have a significant impact on speed of computation.

maxn0

the maximum number of initial molecules used in computation.

mod.method

character string, name of optimsation method used by optim to find posterior mode.

mod.rep

maximum number of iterations used by optim to find mode, which is used to initialise MCMC process if mc.rep>0. If 0 then mode not found.

burnin

proportion of MCMC samples discarded as burn-in.

nreport

number of times the function value and acceptance rate of the MCMC sample are printed.

cq.xlim

if NA then no rows trimmed, if vector then cq.xlim is the start and end rows for which Cq data is retained, with positive and negative counts only kept for remaining trimmed rows.

tune

the tuning parameter for the Metropolis sampling. If a vector, then the same length as the parameter vector.

E.init, E1.init

vectors of initial parameter values (probabilities) for E and E_1.

d.init

initial parameter value (positive) for dispersion (ν).

prior

logical. if FALSE use constant priors, in which case the mode is equivalent to the MLE. Otherwise use priors specified by mu.fun, A.fun, E.fun, E1.fun, trendx.fun, trendy.fun and disp.fun.

mu.fun, A.fun, E.fun, E1.fun, trendx.fun, trendy.fun, disp.fun

priors as single parameter functions.

Details

This function can be used to find the posterior mode using link{optim} and/or simulate a posterior MCMC sample of the Enhanced dPCR model using MCMCmetrop1R. If prior is FALSE then the former is equivalent to finding the Maximum Likelihood estimate (MLE).

The mode found for a particular set of initial paramemeters may not be the global maximum. For this reason it is recommended to search from a number combinations of E and E1 set through the vectors E.init and E1.init.

If report>0 then "." is printed when the posterior is calculated as 0 (due to computational limitations). If report>0 and probreport>0 then "'" is printed with probability probreport when a non-zero posterior is calculated. If optim fails to find non-zero posteriors, and thus fails to work then a better initial value may be required.

Value

A list with components:

cqdata

the cqdata object containing the data.

counts

a vector of the count data.

pars0

initial parameters.

logval0

log-posterior at pars0.

h

threshold value.

nx, ny

column and row numbers.

Also included if mod is TRUE:

pars.mod

parameters at posterior mode.

logval.mod

log-posterior at pars.mod.

mod.res

matrix of parameters and log-posteriors found from running optim at different starting values of E and E1.

Also included if mc is TRUE:

pars.mc

sample mean of MCMC sample of posterior.

logval.mc

log-posterior at pars.mc

pars.sum

summary of MCMC sample of posterior.

pars.chain

MCMC sample of posterior.

mc.vals

log-posteriors associated with MCMC sample of posterior

mc.acc

acceptance rate.

mcmod.pars

posterior mode from MCMC posterior sample.

mcmod.val

log-posterior at mcmod.pars.

prior.funs

list of prior distributions.

Author(s)

Philip Wilson

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
## Not run: 
# These examples take some time to run.

dat<-fetch(Exp37a,panel=1)
# Estimate posterior mode for single intitial parameter vector
res1<-cqmc(dat,mc.rep=0,E.init=.9,E1.init=.8,probreport=.1)

res2<-cqmc(dat,mod.rep=10^3,mc.rep=10^3,E.init=.9,E1.init=c(.5,.7,.9),probreport=.1)
# Estimate posterior mode from 9 initial parameter vectors. 
# Simulate MCMC sample from posterior distribution 
#  using estimated mode as initial value.
res3<-cqmc(dat,mc.rep=10^4,
  E.init=rep(c(.9,.7),each=3),
  E1.init=rep(c(.9,.7,.5),2),probreport=.1) 

## End(Not run)

edpcr documentation built on May 2, 2019, 5:22 p.m.