aroma.core (version 3.3.1)

fitGenotypeCone.matrix: 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 matrix
fitGenotypeCone(y, flavor=c("sfit", "expectile"), ...)

Value

Returns a named list structure.

Arguments

y

A numeric Nx2 matrix with one column for each allele and where N is the number of data points.

flavor

A character string specifying what model/algorithm should be used to fit the genotype cone.

...

Additional arguments passed to the internal fit function.

Author

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)