## S3 method for class 'matrix'
fitPrincipalCurve(X, ..., verbose=FALSE)
Arguments
X
An NxK matrix (K>=2) where the columns represent the dimension.
...
Other arguments passed to principal.curve.
verbose
A logical or a Verbose object.
Value
Returns a principal.curve object (which is a list).
See principal.curve for more details.
Missing values
The estimation of the normalization function will only be made
based on complete observations, i.e. observations that contains no NA
values in any of the channels.
Author(s)
Henrik Bengtsson
References
[1] Hastie, T. and Stuetzle, W, Principal Curves, JASA, 1989.
[2] H. Bengtsson, A. Ray, P. Spellman and T.P. Speed, A single-sample method for normalizing and combining full-resolutioncopy numbers from multiple platforms, labs and analysis methods, Bioinformatics, 2009.
See Also
backtransformPrincipalCurve().
principal.curve.
Examples
# Simulate data from the model y <- a + bx + x^c + eps(bx)
J <- 1000
x <- rexp(J)
a <- c(2,15,3)
b <- c(2,3,4)
c <- c(1,2,1/2)
bx <- outer(b,x)
xc <- t(sapply(c, FUN=function(c) x^c))
eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(b), mean=0, sd=0.1*x))
y <- a + bx + xc + eps
y <- t(y)
# Fit principal curve through (y_1, y_2, y_3)
fit <- fitPrincipalCurve(y, verbose=TRUE)
# Flip direction of 'lambda'?
rho <- cor(fit$lambda, y[,1], use="complete.obs")
flip <- (rho < 0)
if (flip) {
fit$lambda <- max(fit$lambda, na.rm=TRUE)-fit$lambda
}
# Backtransform (y_1, y_2, y_3) to be proportional to each other
yN <- backtransformPrincipalCurve(y, fit=fit)
# Same backtransformation dimension by dimension
yN2 <- y
for (cc in 1:ncol(y)) {
yN2[,cc] <- backtransformPrincipalCurve(y, fit=fit, dimensions=cc)
}
stopifnot(identical(yN2, yN))
xlim <- c(0, 1.04*max(x))
ylim <- range(c(y,yN), na.rm=TRUE)
# Pairwise signals vs x before and after transform
layout(matrix(1:4, nrow=2, byrow=TRUE))
par(mar=c(4,4,3,2)+0.1)
for (cc in 1:3) {
ylab <- substitute(y[c], env=list(c=cc))
plot(NA, xlim=xlim, ylim=ylim, xlab="x", ylab=ylab)
abline(h=a[cc], lty=3)
mtext(side=4, at=a[cc], sprintf("a=%g", a[cc]),
cex=0.8, las=2, line=0, adj=1.1, padj=-0.2)
points(x, y[,cc])
points(x, yN[,cc], col="tomato")
legend("topleft", col=c("black", "tomato"), pch=19,
c("orignal", "transformed"), bty="n")
}
title(main="Pairwise signals vs x before and after transform", outer=TRUE, line=-2)
# Pairwise signals before and after transform
layout(matrix(1:4, nrow=2, byrow=TRUE))
par(mar=c(4,4,3,2)+0.1)
for (rr in 3:2) {
ylab <- substitute(y[c], env=list(c=rr))
for (cc in 1:2) {
if (cc == rr) {
plot.new()
next
}
xlab <- substitute(y[c], env=list(c=cc))
plot(NA, xlim=ylim, ylim=ylim, xlab=xlab, ylab=ylab)
abline(a=0, b=1, lty=2)
points(y[,c(cc,rr)])
points(yN[,c(cc,rr)], col="tomato")
legend("topleft", col=c("black", "tomato"), pch=19,
c("orignal", "transformed"), bty="n")
}
}
title(main="Pairwise signals before and after transform", outer=TRUE, line=-2)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(aroma.light)
aroma.light v3.2.0 (2016-01-06) successfully loaded. See ?aroma.light for help.
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/aroma.light/fitPrincipalCurve.Rd_%03d_medium.png", width=480, height=480)
> ### Name: fitPrincipalCurve
> ### Title: Fit a principal curve in K dimensions
> ### Aliases: fitPrincipalCurve fitPrincipalCurve.matrix
> ### Keywords: methods
>
> ### ** Examples
>
>
> # Simulate data from the model y <- a + bx + x^c + eps(bx)
> J <- 1000
> x <- rexp(J)
> a <- c(2,15,3)
> b <- c(2,3,4)
> c <- c(1,2,1/2)
> bx <- outer(b,x)
> xc <- t(sapply(c, FUN=function(c) x^c))
> eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(b), mean=0, sd=0.1*x))
> y <- a + bx + xc + eps
> y <- t(y)
>
> # Fit principal curve through (y_1, y_2, y_3)
> fit <- fitPrincipalCurve(y, verbose=TRUE)
Fitting principal curve...
Data size: 1000x3
Identifying missing values...
Identifying missing values...done
Data size after removing non-finite data points: 1000x3
Calling principal.curve()...
Starting curve---distance^2: 1497703
Iteration 1---distance^2: 332.6797
Iteration 2---distance^2: 332.3741
Converged: TRUE
Number of iterations: 2
Processing time/iteration: 0.0s (0.0s/iteration)
Calling principal.curve()...done
Fitting principal curve...done
>
> # Flip direction of 'lambda'?
> rho <- cor(fit$lambda, y[,1], use="complete.obs")
> flip <- (rho < 0)
> if (flip) {
+ fit$lambda <- max(fit$lambda, na.rm=TRUE)-fit$lambda
+ }
>
>
> # Backtransform (y_1, y_2, y_3) to be proportional to each other
> yN <- backtransformPrincipalCurve(y, fit=fit)
>
> # Same backtransformation dimension by dimension
> yN2 <- y
> for (cc in 1:ncol(y)) {
+ yN2[,cc] <- backtransformPrincipalCurve(y, fit=fit, dimensions=cc)
+ }
> stopifnot(identical(yN2, yN))
>
>
> xlim <- c(0, 1.04*max(x))
> ylim <- range(c(y,yN), na.rm=TRUE)
>
>
> # Pairwise signals vs x before and after transform
> layout(matrix(1:4, nrow=2, byrow=TRUE))
> par(mar=c(4,4,3,2)+0.1)
> for (cc in 1:3) {
+ ylab <- substitute(y[c], env=list(c=cc))
+ plot(NA, xlim=xlim, ylim=ylim, xlab="x", ylab=ylab)
+ abline(h=a[cc], lty=3)
+ mtext(side=4, at=a[cc], sprintf("a=%g", a[cc]),
+ cex=0.8, las=2, line=0, adj=1.1, padj=-0.2)
+ points(x, y[,cc])
+ points(x, yN[,cc], col="tomato")
+ legend("topleft", col=c("black", "tomato"), pch=19,
+ c("orignal", "transformed"), bty="n")
+ }
> title(main="Pairwise signals vs x before and after transform", outer=TRUE, line=-2)
>
>
> # Pairwise signals before and after transform
> layout(matrix(1:4, nrow=2, byrow=TRUE))
> par(mar=c(4,4,3,2)+0.1)
> for (rr in 3:2) {
+ ylab <- substitute(y[c], env=list(c=rr))
+ for (cc in 1:2) {
+ if (cc == rr) {
+ plot.new()
+ next
+ }
+ xlab <- substitute(y[c], env=list(c=cc))
+ plot(NA, xlim=ylim, ylim=ylim, xlab=xlab, ylab=ylab)
+ abline(a=0, b=1, lty=2)
+ points(y[,c(cc,rr)])
+ points(yN[,c(cc,rr)], col="tomato")
+ legend("topleft", col=c("black", "tomato"), pch=19,
+ c("orignal", "transformed"), bty="n")
+ }
+ }
> title(main="Pairwise signals before and after transform", outer=TRUE, line=-2)
>
>
>
>
>
> dev.off()
null device
1
>