fitMultiDimensionalCone.matrix: Fits an affine transformation to multi-dimensional data

In HenrikBengtsson/aroma.core: Core Methods and Classes Used by 'aroma.*' Packages Part of the Aroma Framework

Fits an affine transformation to multi-dimensional data

Description

Fits an affine transformation to multi-dimensional data using robust estimators.

Usage

## S3 method for class 'matrix'
fitMultiDimensionalCone(y, alpha=c(0.1, 0.075, 0.05, 0.03, 0.01), q=2, Q=98, ...,
  flavor=c("sfit", "expectile"))

Arguments

y	A numeric NxK matrix with one column for each dimension and where N is the number of data points.
alpha	A numeric vector of decreasing values in (0,1). This parameter "determines how far we are willing to press the boundary of the [genotype cone]". Lowering alpha expand the cone. When alpha goes to zero, all data points will be on or inside the cone.
q,Q	Percentiles in [0,100] for which data points that are below (above) will be assigned zero weight in the fitting of the parameters.
...	Additional arguments passed to the cfit() function of the sfit package.
flavor	A character string specifying what model/algorithm should be used to fit the genotype cone.

Value

Returns the parameter estimates as a named list with elements:

M	An estimate of the three vertices defining the genotype triangle. These three vertices are describes as an 2x3 matrix with column origin, AA, and BB.
Minv	The inverse of M.
origin	The estimate of the shift.
W	The estimate of shear/rotation matrix with columns AA and BB.
Winv	The inverse of W.
params	The parameters used for the fit, i.e. alpha, q, Q, and those passed in ....
dimData	The dimension of the input data.

Author(s)

Henrik Bengtsson

See Also

To backtransform data fitted using this method, see *backtransformMultiDimensionalCone(). Internally, the cfit() function the sfit package is used.

Examples

if (require("sfit")) {
 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Simulate data (taken from the cfit.matrix() example of 'sfit')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
N <- 1000

# Simulate sequences
nucleotides <- c("A", "C", "G", "T")
g <- sample(nucleotides, size=N, replace=TRUE)
ndim <- length(nucleotides)

# Simulate concentrations of allele A and allele B
X <- matrix(rexp(N), nrow=N, ncol=ndim)
colnames(X) <- nucleotides
for (nucleotide in nucleotides) {
  cc <- match(nucleotide, nucleotides)
  X[g == nucleotide, -cc] <- 0
}

# The true offset
a0 <- rep(0.3, times=ndim)

# The crosstalk matrix
A <- matrix(c(
  0.9, 0.3, 0.2, 0.1,
  0.1, 0.8, 0.1, 0.1,
  0.3, 0.4, 0.7, 0.1,
  0.1, 0.1, 0.6, 0.9
), nrow=ndim, byrow=TRUE)
A <- apply(A, MARGIN=2, FUN=function(u) u / sqrt(sum(u^2)))

# Simulate random errors on the input
xi <- matrix(rnorm(ndim*N, mean=0, sd=0.05), ncol=ndim)

# Generate the noisy crosstalk affected input data
Z <- t(a0 + A %*% t(X + xi))

# Generate the noisy observations of the latter
eps <- matrix(rnorm(ndim*N, mean=0, sd=0.05), ncol=ndim)
Y <- Z + eps

# Fit crosstalk model and calibrate data accordingly
fit <- fitMultiDimensionalCone(Y, flavor="sfit")
Yc <- backtransformMultiDimensionalCone(Y, fit=fit)

lim <- c(-0.5,6)
layout(matrix(c(1,2,3,0,4,5,0,0,6), nrow=3, ncol=3, byrow=TRUE))
par(mar=c(5,4,1,1)+0.1)
for (ii in 1:(ndim-1)) {
  for (jj in (ii+1):ndim) {
    cc <- c(jj,ii)
    labs <- nucleotides[cc]
    plot(Y[,cc], cex=0.8, xlim=lim, ylim=lim, xlab=labs[1], ylab=labs[2])
    points(Yc[,cc], cex=0.8, col="blue")
    Mcc <- fit$M[c(1,1+cc),cc]
    class(Mcc) <- class(fit$M)
    lines(Mcc, lwd=2, col="red")
  }
}

}

HenrikBengtsson/aroma.core documentation built on Feb. 20, 2024, 9:17 p.m.