fitGenotypeConeBySfit.matrix: Fits an affine transformation to allele A and allele B data

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

Fits an affine transformation to allele A and allele B data

Description

Fits an affine transformation to allele A and allele B data using robust estimators.

Usage

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

Arguments

y	A numeric Nx2 matrix with one column for each allele 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() of the sfit package.

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 backtransformGenotypeCone(). Internally, the cfit() function the sfit package is used.

Examples

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Fit genotype cone based on available methods (==packages)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
flavors <- c("expectile", "sfit")
## The 'expectile' package does not do proper S3 registration
if (getRversion() >= "3.6.0") flavors <- setdiff(flavors, "expectile")
keep <- sapply(flavors, FUN=require, character.only=TRUE)
flavors <- flavors[keep]
cat("Available fitting flavors:", paste(flavors, collapse=", "), "\n")
hasSfit <- is.element("sfit", flavors)


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Simulate data (taken from the cfit.matrix() example of 'sfit')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#set.seed(0xbeef)

N <- 1000

# Simulate genotypes
g <- sample(c("AA", "AB", "AB", "BB"), size=N, replace=TRUE)

# Simulate concentrations of allele A and allele B
X <- matrix(rexp(N), nrow=N, ncol=2)
colnames(X) <- c("A", "B")
X[g == "AA", "B"] <- 0
X[g == "BB", "A"] <- 0
X[g == "AB",] <- X[g == "AB",] / 2

# Transform noisy X
xi <- matrix(rnorm(2*N, mean=0, sd=0.05), ncol=2)
a0 <- c(0,0)+0.3
A <- matrix(c(0.9, 0.1, 0.1, 0.8), nrow=2, byrow=TRUE)
A <- apply(A, MARGIN=2, FUN=function(u) u / sqrt(sum(u^2)))
Z <- t(a0 + A %*% t(X + xi))

# Add noise to Y
eps <- matrix(rnorm(2*N, mean=0, sd=0.05), ncol=2)
Y <- Z + eps


lim <- c(-1/2,6)
xlab <- "Allele A"
ylab <- "Allele B"
plot(Y, xlab=xlab, ylab=ylab, xlim=lim, ylim=lim)
lines(x=c(0,0,2*lim[2]), y=c(2*lim[2],0,0), col="#aaaaaa", lwd=3, lty=3)

# Different alpha sequences to illustrate the impact
alphas <- c(0.075, 0.05, 0.01, 0.03, 0.002, 0.001)

cols <- seq(from=2, to=length(alphas)+1)
legend("topright", sprintf("%.3f", alphas), col=cols, lwd=4, title="Alphas")

for (kk in seq_along(alphas)) {
  for (flavor in flavors) {
    fit <- fitGenotypeCone(Y, alpha=alphas[kk], flavor=flavor)
    YN <- backtransformGenotypeCone(Y, fit=fit)
    if (hasSfit) {
      radials(fit$M, lwd=3, col=cols[kk], lty=ifelse(flavor == "sfit", 1,2))
      drawApex(fit$M, col=cols[kk], pch=19, cex=2)
    }
    points(YN, col=cols[kk])
  }
}
lines(x=c(0,0,2*lim[2]), y=c(2*lim[2],0,0), col="#aaaaaa", lwd=3, lty=3)

aroma.core documentation built on June 25, 2024, 1:15 a.m.