RJaCGH(y, Chrom = NULL, Pos = NULL, Dist=NULL, model = "genome",
max.dist=NULL, normal.reference=0, normal.ref.percentile=0.95,
burnin = 10000, TOT =10000, k.max = 6, stat = NULL, mu.alfa = NULL,
mu.beta = NULL, ka = NULL,
g = NULL, prob.k = NULL, jump.parameters=list(), start.k = NULL,
RJ=TRUE, auto.label=0.60)length(y)-1. Note that when Chrom is not NULL,
every last value of every Chromosome is not used.model="genome", the same model is fitted for
the whole genome. If model="Chrom", a different model is
fitted for each chromosome.max.dist, they
are considered independent. That is, the state of that spot
does not affect the state of the other. If NULL0.k.max. If NULL, it is assumed a
uniform distribution for every model.prob.k is chosen.RJaCGH.array is returned, with components:model is "genome", an object of class RJaCGH.genome
is returned, with components:model is "Chrom", an object of class RJaCGH.Chrom is
returned, with the following components:RJaCGH (See below).model was specified and no Chrom was given, an
object of class RJaCGH is returned, with components k,
prob.b, prob.d, prob.s, prob.c, y,
Pos, x as described before, plus a list with as many
components of number of max hidden states fitted.
For every hidden markov model fitted, a list is returned with
components:mu in the
Metropolis-Hastings step.sigma.2 in the
Metropolis-Hastings step.beta in the
Metropolis-Hastings step.prob.k. If NULL, a uniform distribution between 1
and k.max is used.
The hidden states follow a normal distribution which mean (mu) follows
itself a normal distribution with mean
mu.alfa and stdev mu.beta. If NULL, these are the
median of the data and the range. The variance (sigma.2)of the hidden states
follows an Inverse Gamma distribution with parameters ka and
g. If NULL, these are 2 and diff(range(y))^2 / 50,
respectively.
The model for the transition matrix is based on a random matrix
beta whose diagonal is zero. The transition matrix, Q,
has the form:
Q[i,j] = exp(-beta[i,j] + beta[i,j]*x) / sum(i,.) {exp(-beta[i,.] +
beta[i,.]*x}
The prior distribution for beta is gamma with parameters 1, 1.
The x are the distances between positions, normalized to lay
between zero and 1 (x=diff(Pos) / max(diff(Pos)))
RJaCGH performs Markov Chain MonteCarlo with Reversible Jump to sample for the posterior distribution. Every sweep has 3 steps:
1.- A Metropolis-Hastings move is used to update, for a fixed number
of hidden states, mu, sigma.2 and beta. A
symmetric proposal with a normal distribution and standard deviation
sigma.tau.mu, sigma.tau.sigma.2 and
sigma.tau.beta is sampled.
2.- A transdimensional move is chosen, between birth (a new hidden state is sampled from the prior) or death (an existing hidden state is erased).
3.- Another transdimensional move is performed; an split move (divide
an existing state in two) or a combine move (join two adjacent
states). The length of the split is sampled from a normal distribution
with standard deviation tau.split.mu for the mu and
tau.split.beta for beta.
jump.parameters must be a list with the parameters for the
moves. It must have components sigma.tau.mu,
sigma.tau.sigma.2, sigma.tau.beta These are vectors of
length k.max. tau.split.mu, tau.split.beta are vectors of
length 1. If any of them is NULL, a call to the internal function
get.jump() is made to find 'good' values.
A relabelling of hidden states is performed to match biological
states. The states whom have the normal.reference value
inside a normal.ref.percentile% confidence
based on a normal distribution with means the median of mu
and sd the square root of the median of sigma.2 are labelled as
'Normal'. If no state is labelled as 'Normal', the closest state to
normal.reference is chosen.
If auto.label is not null, closest states to 'Normal' are also
labelled as 'Normal' until a proportion of auto.label is
reached. Please note that the default value is 0.60, so at least the
60% of the observations will be labelled as 'Normal'.
summary.RJaCGH,
states, model.averaging,
plot.RJaCGH, trace.plot,
gelman.brooks.plot, collapseChainy <- c(rnorm(100, 0, 1), rnorm(10, -3, 1), rnorm(20, 3, 1),
rnorm(100,0, 1))
Pos <- sample(x=1:500, size=230, replace=TRUE)
Pos <- cumsum(Pos)
Chrom <- rep(1:23, rep(10, 23))
jp <- list(sigma.tau.mu=rep(0.05, 4), sigma.tau.sigma.2=rep(0.03, 4),
sigma.tau.beta=rep(0.07, 4), tau.split.mu=0.1, tau.split.beta=0.1)
fit.chrom <- RJaCGH(y=y, Pos=Pos, Chrom=Chrom, model="Chrom",
burnin=10, TOT=1000, k.max = 4,
jump.parameters=jp)
##RJ results for chromosome 5
table(fit.chrom[[5]]$k)
fit.genome <- RJaCGH(y=y, Pos=Pos, Chrom=Chrom, model="genome",
burnin=100, TOT=1000, jump.parameters=jp, k.max = 4)
## Results for the model with 3 states:
apply(fit.genome[[3]]$mu, 2, summary)
apply(fit.genome[[3]]$sigma.2, 2, summary)
apply(fit.genome[[3]]$beta, c(1,2), summary)Run the code above in your browser using DataLab