noise.model: Estimate model parameters for noise

In BioSSA/BioSSA: A collection of methods for singular spectrum analysis of drosophila embryo (2d and 3d) gene expression

Description

Functions for storing, plotting and model construction for noise

Usage

## Default S3 method:
noise.model(x, trend,
        offset = 0,
        model = "estimate",
        reg.type = c("winsor", "trim"),
        reg.level = 0,
        averaging.type = c("sliding-window", "quantile-break", "equal-break", "none"),
        breaks = 2, window = 51,
        FUN = median,
        FUN.trend = FUN,
        na.rm = TRUE,
        ...)
## S3 method for class 'BioSSA2d'
noise.model(x, groups, ...)
## S3 method for class 'noise.model'
plot(x,
        absolute = FALSE, #TODO Mb, remove it????
        relative = TRUE,
        draw.residuals = TRUE,
        draw.means.fitted = FALSE,
        print.alpha = TRUE,
        ref = TRUE,
        symmetric = !absolute,
        ...,
        dots.residuals = list(),
        dots.means.fitted = list(),
        digits = max(3, getOption("digits") - 3))
## S3 method for class 'noise.model'
summary(object, digits = max(3, getOption("digits") - 3), ...)
## S3 method for class 'noise.model'
print(x, digits = max(3, getOption("digits") - 3), ...)

Arguments

x, object	the input object. This might be ‘BioSSA2d’/‘BioSSA3d’ object for BioSSA2d/BioSSA3d method, or just a numeric vector of residuals for default implementation
trend	numeric vector, trend for noise model estimation; this parameter can be used only for default method
offset	numeric value, trend offset for noise model estimation
model	model name ('additive', 'multiplicative' or 'poisson') or multiplicity power value. Use 'estimate' for perform model estimation
reg.type	regularization type for residuals and trend values averaging
reg.level	numeric value, quantile level for regularization
averaging.type	character, averaging method for power estimation and model visualization
breaks	the number of intervals to trend values breaking. Used for 'equal-break' and 'quantile-break' averaging types
window	length of sliding window. Used for 'sliding-window' averaging type
FUN	averaging function for logarithms of residuals absolute values (like mean or median) should take numeric vector and return one numeric value. Used for all averaging types except 'none'
FUN.trend	averaging function for logarithms of trend values, like previous
na.rm	a logical value indicating whether 'NA' values should be stripped before the computation proceeds
...	additional arguments, passed to inner function calls
groups	numeric vector, component indices in BioSSA decomposition for trend extraction
absolute	logical, whether plot absolute values of residuals instead of residuals itself
relative	logical, whether plot relative residuals instead of absolute ones
draw.residuals	logical, whether plot residuals
draw.means.fitted	logical, whether plot regression line drawn by averaged residuals
print.alpha	logical, whether output multiplicity power estimation
ref	logical, whether plot zero level line
symmetric	logical, whether y-scale should be symmetric by default
dots.residuals	list of additional arguments passed to xyplot function used for plotting of the residuals
dots.means.fitted	same as previous, but for regression line plot (not used, supposed to be excluded)
digits	integer, how many significant digits are to be used for numbers formatting

Details

noise.model applied to ‘BioSSA2d’/‘BioSSA3d’ object calculates the residual and trend and then call the default method passing the further arguments.

plot method plots residuals and smoothed residuals in dependence on the trend.

Data consisting of trend and residual values are ordered by the trend values, then FUN function is applied to residuals and the result is depicted again FUN.trend applied to the ordered trend values.

Value

Object of class ‘noise.model’ for noise.model. Trellis object for plot. Original object (invisibly) for print and summary.

Object of ‘noise.model’ is a list with following fields:

alpha

is numeric value, which is estimated if model = 'estimate', or is equal to the user specified value (multiplicity power)

sigma

is standard deviation of the relative noise, (see ‘Model description’)

sd

is square root of mean square deviation of relative residuals

residuals

is vector of residuals

trend

is vector of trend values

residuals.means

is vector of smoothed (averaged) residuals absolute values

residuals.means.fitted

is vector of fitted values (by lm) of smoothed residuals absolute values

trend.means

is vector of smoothed trend values

offset

is offset value

averaging.type

is character name of used averaging type

call

is object of class ‘call’, constructor call (MB, drop it? This is command, which has created object, ‘call’ means ‘function call’)

Model description

Generalized multiplicative noise model is considered:

res_i = σ (trend_i + offset)^α \cdot ξ_i,

where ξ_i have standard normal distribution.

Value of alpha is estimated as follows: logarithms of absolute values of residuals and ofsetted trend are considered, then they both ordered by the trend values and averaging procedure performed. There are following averaging methods: 'none' means nothing averaging, 'sliding-window' means sliding window averaging (default approach) with window length denoted by window argument, 'equal-break' and 'quantile-break' mean splitting all trend values into buckets ('quantile-break' means buckets with equal quantity of elements in each bucket and 'equal-break' means equal size buckets), correspondingly residuals splitting and averaging logarithmed residuals and trend values in each bucket.

Then linear regression on averaged logarithms is provided. Slope is alpha estimation and intercept is estimation of logarithm of the standard deviation of relative noise, i.e. sigma = exp(Intercept).

Examples

xlim <- c(22, 88)
ylim <- c(32, 68)
L <- c(15, 15)

file <- system.file("extdata/data", "ab16.txt", package = "BioSSA")
df <- read.emb.data(file)

bs <- BioSSA(cad ~ AP + DV, data = df, ylim = ylim, xlim = xlim, L = L)
nm <- noise.model(bs, 1:3, averaging.type = "none")
plot(nm)
summary(nm)

nm <- noise.model(bs, 1:3, averaging.type = "sliding")
plot(nm)
summary(nm)

nm <- noise.model(bs, 1:3, averaging.type = "equal")
plot(nm)
summary(nm)

nm <- noise.model(bs, 1:3, averaging.type = "quantile")
plot(nm)
summary(nm)

nm <- noise.model(bs, 1:3, model = "poisson")
plot(nm)
summary(nm)

nm <- noise.model(bs, 1:3, model = "additive")
plot(nm)
summary(nm)

nm <- noise.model(bs, 1:3, model = "multiplicative")
plot(nm)
summary(nm)

nm <- noise.model(bs, 1:3, model = -1.2)
plot(nm)
summary(nm)



nm.none <- noise.model(bs, 1:3, model = "estimate", averaging.type = "none")
nm <- noise.model(bs, 1:3, model = nm.none$alpha, averaging.type = "sliding")

plot(nm)

BioSSA/BioSSA documentation built on May 5, 2019, 3:47 p.m.