| Back to Multiple platform build/check report for BioC 3.17: simplified long |
|
This page was generated on 2023-10-16 11:35:11 -0400 (Mon, 16 Oct 2023).
| Hostname | OS | Arch (*) | R version | Installed pkgs |
|---|---|---|---|---|
| nebbiolo1 | Linux (Ubuntu 22.04.2 LTS) | x86_64 | 4.3.1 (2023-06-16) -- "Beagle Scouts" | 4626 |
| palomino3 | Windows Server 2022 Datacenter | x64 | 4.3.1 (2023-06-16 ucrt) -- "Beagle Scouts" | 4379 |
| merida1 | macOS 12.6.4 Monterey | x86_64 | 4.3.1 (2023-06-16) -- "Beagle Scouts" | 4395 |
| Click on any hostname to see more info about the system (e.g. compilers) (*) as reported by 'uname -p', except on Windows and Mac OS X | ||||
| Package 87/2230 | Hostname | OS / Arch | INSTALL | BUILD | CHECK | BUILD BIN | ||||||||
| aroma.light 3.30.0 (landing page) Henrik Bengtsson
| nebbiolo1 | Linux (Ubuntu 22.04.2 LTS) / x86_64 | OK | OK | OK |  | ||||||||
| palomino3 | Windows Server 2022 Datacenter / x64 | OK | OK | OK | OK | | ||||||||
| merida1 | macOS 12.6.4 Monterey / x86_64 | OK | OK | OK | OK | | ||||||||
| kjohnson2 | macOS 12.6.1 Monterey / arm64 | see weekly results here | ||||||||||||
|
To the developers/maintainers of the aroma.light package: - Allow up to 24 hours (and sometimes 48 hours) for your latest push to git@git.bioconductor.org:packages/aroma.light.git to reflect on this report. See Troubleshooting Build Report for more information. - Use the following Renviron settings to reproduce errors and warnings. - If 'R CMD check' started to fail recently on the Linux builder(s) over a missing dependency, add the missing dependency to 'Suggests:' in your DESCRIPTION file. See Renviron.bioc for more information. |
| Package: aroma.light |
| Version: 3.30.0 |
| Command: /home/biocbuild/bbs-3.17-bioc/R/bin/R CMD check --install=check:aroma.light.install-out.txt --library=/home/biocbuild/bbs-3.17-bioc/R/site-library --timings aroma.light_3.30.0.tar.gz |
| StartedAt: 2023-10-15 19:13:50 -0400 (Sun, 15 Oct 2023) |
| EndedAt: 2023-10-15 19:14:56 -0400 (Sun, 15 Oct 2023) |
| EllapsedTime: 66.0 seconds |
| RetCode: 0 |
| Status: OK |
| CheckDir: aroma.light.Rcheck |
| Warnings: 0 |
##############################################################################
##############################################################################
###
### Running command:
###
### /home/biocbuild/bbs-3.17-bioc/R/bin/R CMD check --install=check:aroma.light.install-out.txt --library=/home/biocbuild/bbs-3.17-bioc/R/site-library --timings aroma.light_3.30.0.tar.gz
###
##############################################################################
##############################################################################
* using log directory ‘/home/biocbuild/bbs-3.17-bioc/meat/aroma.light.Rcheck’
* using R version 4.3.1 (2023-06-16)
* using platform: x86_64-pc-linux-gnu (64-bit)
* R was compiled by
gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0
GNU Fortran (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0
* running under: Ubuntu 22.04.3 LTS
* using session charset: UTF-8
* checking for file ‘aroma.light/DESCRIPTION’ ... OK
* this is package ‘aroma.light’ version ‘3.30.0’
* package encoding: latin1
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for hidden files and directories ... NOTE
Found the following hidden files and directories:
inst/rsp/.rspPlugins
These were most likely included in error. See section ‘Package
structure’ in the ‘Writing R Extensions’ manual.
* checking for portable file names ... OK
* checking for sufficient/correct file permissions ... OK
* checking whether package ‘aroma.light’ can be installed ... OK
* checking installed package size ... OK
* checking package directory ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* checking whether the package can be loaded ... OK
* checking whether the package can be loaded with stated dependencies ... OK
* checking whether the package can be unloaded cleanly ... OK
* checking whether the namespace can be loaded with stated dependencies ... OK
* checking whether the namespace can be unloaded cleanly ... OK
* checking loading without being on the library search path ... OK
* checking startup messages can be suppressed ... OK
* checking dependencies in R code ... OK
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... OK
* checking Rd files ... NOTE
checkRd: (-1) averageQuantile.Rd:49: Escaped LaTeX specials: \&
checkRd: (-1) normalizeQuantileRank.Rd:97: Escaped LaTeX specials: \&
checkRd: (-1) normalizeQuantileRank.matrix.Rd:83: Escaped LaTeX specials: \&
* checking Rd metadata ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking examples ... OK
Examples with CPU (user + system) or elapsed time > 5s
user system elapsed
normalizeCurveFit 5.890 0.016 5.907
normalizeAffine 5.653 0.029 5.684
* checking for unstated dependencies in ‘tests’ ... OK
* checking tests ...
Running ‘backtransformAffine.matrix.R’
Running ‘backtransformPrincipalCurve.matrix.R’
Running ‘callNaiveGenotypes.R’
Running ‘distanceBetweenLines.R’
Running ‘findPeaksAndValleys.R’
Running ‘fitPrincipalCurve.matrix.R’
Running ‘fitXYCurve.matrix.R’
Running ‘iwpca.matrix.R’
Running ‘likelihood.smooth.spline.R’
Running ‘medianPolish.matrix.R’
Running ‘normalizeAffine.matrix.R’
Running ‘normalizeAverage.list.R’
Running ‘normalizeAverage.matrix.R’
Running ‘normalizeCurveFit.matrix.R’
Running ‘normalizeDifferencesToAverage.R’
Running ‘normalizeFragmentLength-ex1.R’
Running ‘normalizeFragmentLength-ex2.R’
Running ‘normalizeQuantileRank.list.R’
Running ‘normalizeQuantileRank.matrix.R’
Running ‘normalizeQuantileSpline.matrix.R’
Running ‘normalizeTumorBoost,flavors.R’
Running ‘normalizeTumorBoost.R’
Running ‘robustSmoothSpline.R’
Running ‘rowAverages.matrix.R’
Running ‘sampleCorrelations.matrix.R’
Running ‘sampleTuples.R’
Running ‘wpca.matrix.R’
Running ‘wpca2.matrix.R’
OK
* checking PDF version of manual ... OK
* DONE
Status: 2 NOTEs
See
‘/home/biocbuild/bbs-3.17-bioc/meat/aroma.light.Rcheck/00check.log’
for details.
aroma.light.Rcheck/00install.out
############################################################################## ############################################################################## ### ### Running command: ### ### /home/biocbuild/bbs-3.17-bioc/R/bin/R CMD INSTALL aroma.light ### ############################################################################## ############################################################################## * installing to library ‘/home/biocbuild/bbs-3.17-bioc/R/site-library’ * installing *source* package ‘aroma.light’ ... ** using staged installation ** R ** inst ** byte-compile and prepare package for lazy loading ** help *** installing help indices ** building package indices ** testing if installed package can be loaded from temporary location ** testing if installed package can be loaded from final location ** testing if installed package keeps a record of temporary installation path * DONE (aroma.light)
aroma.light.Rcheck/tests/backtransformAffine.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> X <- matrix(1:8, nrow=4, ncol=2)
> X[2,2] <- NA_integer_
>
> print(X)
[,1] [,2]
[1,] 1 5
[2,] 2 NA
[3,] 3 7
[4,] 4 8
>
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=c(1,5)))
[,1] [,2]
[1,] 0 0
[2,] 1 NA
[3,] 2 2
[4,] 3 3
>
> # Returns a 4x2 matrix
> print(backtransformAffine(X, b=c(1,1/2)))
[,1] [,2]
[1,] 1 10
[2,] 2 NA
[3,] 3 14
[4,] 4 16
>
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=matrix(1:4,ncol=1)))
[,1] [,2]
[1,] 0 4
[2,] 0 NA
[3,] 0 4
[4,] 0 4
>
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=matrix(1:3,ncol=1)))
[,1] [,2]
[1,] 0 4
[2,] 0 NA
[3,] 0 4
[4,] 3 7
>
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=matrix(1:2,ncol=1), b=c(1,2)))
[,1] [,2]
[1,] 0 2
[2,] 0 NA
[3,] 2 3
[4,] 2 3
>
> # Returns a 4x1 matrix
> print(backtransformAffine(X, b=c(1,1/2), project=TRUE))
[,1]
[1,] 2.8
[2,] 1.6
[3,] 5.2
[4,] 6.4
>
> # If the columns of X are identical, and a identity
> # backtransformation is applied and projected, the
> # same matrix is returned.
> X <- matrix(1:4, nrow=4, ncol=3)
> Y <- backtransformAffine(X, b=c(1,1,1), project=TRUE)
> print(X)
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 2 2
[3,] 3 3 3
[4,] 4 4 4
> print(Y)
[,1]
[1,] 1
[2,] 2
[3,] 3
[4,] 4
> stopifnot(sum(X[,1]-Y) <= .Machine$double.eps)
>
>
> # If the columns of X are identical, and a identity
> # backtransformation is applied and projected, the
> # same matrix is returned.
> X <- matrix(1:4, nrow=4, ncol=3)
> X[,2] <- X[,2]*2; X[,3] <- X[,3]*3
> print(X)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 2 4 6
[3,] 3 6 9
[4,] 4 8 12
> Y <- backtransformAffine(X, b=c(1,2,3))
> print(Y)
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 2 2
[3,] 3 3 3
[4,] 4 4 4
> Y <- backtransformAffine(X, b=c(1,2,3), project=TRUE)
> print(Y)
[,1]
[1,] 1
[2,] 2
[3,] 3
[4,] 4
> stopifnot(sum(X[,1]-Y) <= .Machine$double.eps)
>
> proc.time()
user system elapsed
0.231 0.034 0.254
aroma.light.Rcheck/tests/backtransformPrincipalCurve.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Consider the case where K=4 measurements have been done
> # for the same underlying signals 'x'. The different measurements
> # have different systematic variation
> #
> # y_k = f(x_k) + eps_k; k = 1,...,K.
> #
> # In this example, we assume non-linear measurement functions
> #
> # f(x) = a + b*x + x^c + eps(b*x)
> #
> # where 'a' is an offset, 'b' a scale factor, and 'c' an exponential.
> # We also assume heteroscedastic zero-mean noise with standard
> # deviation proportional to the rescaled underlying signal 'x'.
> #
> # Furthermore, we assume that measurements k=2 and k=3 undergo the
> # same transformation, which may illustrate that the come from
> # the same batch. However, when *fitting* the model below we
> # will assume they are independent.
>
> # Transforms
> a <- c(2, 15, 15, 3)
> b <- c(2, 3, 3, 4)
> c <- c(1, 2, 2, 1/2)
> K <- length(a)
>
> # The true signal
> N <- 1000
> x <- rexp(N)
>
> # The noise
> bX <- outer(b,x)
> E <- apply(bX, MARGIN=2, FUN=function(x) rnorm(K, mean=0, sd=0.1*x))
>
> # The transformed signals with noise
> Xc <- t(sapply(c, FUN=function(c) x^c))
> Y <- a + bX + Xc + E
> Y <- t(Y)
>
>
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Fit principal curve
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Fit principal curve through Y = (y_1, y_2, ..., y_K)
> fit <- fitPrincipalCurve(Y)
>
> # 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
+ }
>
> L <- ncol(fit$s)
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Backtransform data according to model fit
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Backtransform toward the principal curve (the "common scale")
> YN1 <- backtransformPrincipalCurve(Y, fit=fit)
> stopifnot(ncol(YN1) == K)
>
>
> # Backtransform toward the first dimension
> YN2 <- backtransformPrincipalCurve(Y, fit=fit, targetDimension=1)
> stopifnot(ncol(YN2) == K)
>
>
> # Backtransform toward the last (fitted) dimension
> YN3 <- backtransformPrincipalCurve(Y, fit=fit, targetDimension=L)
> stopifnot(ncol(YN3) == K)
>
>
> # Backtransform toward the third dimension (dimension by dimension)
> # Note, this assumes that K == L.
> YN4 <- Y
> for (cc in 1:L) {
+ YN4[,cc] <- backtransformPrincipalCurve(Y, fit=fit,
+ targetDimension=1, dimensions=cc)
+ }
> stopifnot(identical(YN4, YN2))
>
>
> # Backtransform a subset toward the first dimension
> # Note, this assumes that K == L.
> YN5 <- backtransformPrincipalCurve(Y, fit=fit,
+ targetDimension=1, dimensions=2:3)
> stopifnot(identical(YN5, YN2[,2:3]))
> stopifnot(ncol(YN5) == 2)
>
>
> # Extract signals from measurement #2 and backtransform according
> # its model fit. Signals are standardized to target dimension 1.
> y6 <- Y[,2,drop=FALSE]
> yN6 <- backtransformPrincipalCurve(y6, fit=fit, dimensions=2,
+ targetDimension=1)
> stopifnot(identical(yN6, YN2[,2,drop=FALSE]))
> stopifnot(ncol(yN6) == 1)
>
>
> # Extract signals from measurement #2 and backtransform according
> # the the model fit of measurement #3 (because we believe these
> # two have undergone very similar transformations.
> # Signals are standardized to target dimension 1.
> y7 <- Y[,2,drop=FALSE]
> yN7 <- backtransformPrincipalCurve(y7, fit=fit, dimensions=3,
+ targetDimension=1)
> stopifnot(ncol(yN7) == 1)
>
> rho <- cor(yN7, yN6)
> print(rho)
[,1]
[1,] 0.9999915
> stopifnot(rho > 0.999)
>
> proc.time()
user system elapsed
0.654 0.066 0.709
aroma.light.Rcheck/tests/callNaiveGenotypes.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> layout(matrix(1:3, ncol=1))
> par(mar=c(2,4,4,1)+0.1)
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A bimodal distribution
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> xAA <- rnorm(n=10000, mean=0, sd=0.1)
> xBB <- rnorm(n=10000, mean=1, sd=0.1)
> x <- c(xAA,xBB)
> fit <- findPeaksAndValleys(x)
> print(fit)
type x density
1 peak 0.0004193677 1.6795936316
2 valley 0.5040406223 0.0004671231
3 peak 0.9947485114 1.6946873067
> calls <- callNaiveGenotypes(x, cn=rep(1,length(x)), verbose=-20)
Calling genotypes from allele B fractions (BAFs)...
Fitting naive genotype model...
Fitting naive genotype model from normal allele B fractions (BAFs)...
Flavor: density
Censoring BAFs...
Before:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.421177 0.001068 0.528332 0.499810 0.997694 1.399127
[1] 20000
After:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-Inf 0.001068 0.528332 0.997694 Inf
[1] 16834
Censoring BAFs...done
Copy number level #1 (C=1) of 1...
Identified extreme points in density of BAF:
type x density
1 peak 0.01501775 1.638625130
2 valley 0.49828621 0.004035251
3 peak 0.97812724 1.645265821
Local minimas ("valleys") in BAF:
type x density
2 valley 0.4982862 0.004035251
Copy number level #1 (C=1) of 1...done
Fitting naive genotype model from normal allele B fractions (BAFs)...done
[[1]]
[[1]]$flavor
[1] "density"
[[1]]$cn
[1] 1
[[1]]$nbrOfGenotypeGroups
[1] 2
[[1]]$tau
[1] 0.4982862
[[1]]$n
[1] 16834
[[1]]$fit
type x density
1 peak 0.01501775 1.638625130
2 valley 0.49828621 0.004035251
3 peak 0.97812724 1.645265821
[[1]]$fitValleys
type x density
2 valley 0.4982862 0.004035251
attr(,"class")
[1] "NaiveGenotypeModelFit" "list"
Fitting naive genotype model...done
Copy number level #1 (C=1) of 1...
Model fit:
$flavor
[1] "density"
$cn
[1] 1
$nbrOfGenotypeGroups
[1] 2
$tau
[1] 0.4982862
$n
[1] 16834
$fit
type x density
1 peak 0.01501775 1.638625130
2 valley 0.49828621 0.004035251
3 peak 0.97812724 1.645265821
$fitValleys
type x density
2 valley 0.4982862 0.004035251
Genotype threshholds [1]: 0.498286209491153
TCN=1 => BAF in {0,1}.
Call regions: A = (-Inf,0.498], B = (0.498,+Inf)
Copy number level #1 (C=1) of 1...done
Calling genotypes from allele B fractions (BAFs)...done
> xc <- split(x, calls)
> print(table(calls))
calls
0 1
10000 10000
> xx <- c(list(x),xc)
> plotDensity(xx, adjust=1.5, lwd=2, col=seq_along(xx), main="(AA,BB)")
> abline(v=fit$x)
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution with missing values
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> xAB <- rnorm(n=10000, mean=1/2, sd=0.1)
> x <- c(xAA,xAB,xBB)
> x[sample(length(x), size=0.05*length(x))] <- NA_real_
> x[sample(length(x), size=0.01*length(x))] <- -Inf
> x[sample(length(x), size=0.01*length(x))] <- +Inf
> fit <- findPeaksAndValleys(x)
> print(fit)
type x density
1 peak 0.0005042883 1.1607776
2 valley 0.2461620955 0.1866657
3 peak 0.4958470800 1.1696998
4 valley 0.7455320644 0.1896356
5 peak 0.9911898717 1.1795212
> calls <- callNaiveGenotypes(x)
> xc <- split(x, calls)
> print(table(calls))
calls
0 0.5 1
9575 9325 9635
> xx <- c(list(x),xc)
> plotDensity(xx, adjust=1.5, lwd=2, col=seq_along(xx), main="(AA,AB,BB)")
> abline(v=fit$x)
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution with clear separation
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> xAA <- rnorm(n=10000, mean=0, sd=0.02)
> xAB <- rnorm(n=10000, mean=1/2, sd=0.02)
> xBB <- rnorm(n=10000, mean=1, sd=0.02)
> x <- c(xAA,xAB,xBB)
> fit <- findPeaksAndValleys(x)
> print(fit)
type x density
1 peak -0.002548063 2.610723e+00
2 valley 0.248126022 3.136140e-05
3 peak 0.495983544 2.602111e+00
4 valley 0.746657629 3.427514e-05
5 peak 0.997331714 2.605358e+00
> calls <- callNaiveGenotypes(x)
> xc <- split(x, calls)
> print(table(calls))
calls
0 0.5 1
10000 10000 10000
> xx <- c(list(x),xc)
> plotDensity(xx, adjust=1.5, lwd=2, col=seq_along(xx), main="(AA',AB',BB')")
> abline(v=fit$x)
>
> proc.time()
user system elapsed
0.516 0.053 0.558
aroma.light.Rcheck/tests/distanceBetweenLines.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> for (zzz in 0) {
+
+ # This example requires plot3d() in R.basic [http://www.braju.com/R/]
+ if (!require(pkgName <- "R.basic", character.only=TRUE)) break
+
+ layout(matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
+
+ ############################################################
+ # Lines in two-dimensions
+ ############################################################
+ x <- list(a=c(1,0), b=c(1,2))
+ y <- list(a=c(0,2), b=c(1,1))
+ fit <- distanceBetweenLines(ax=x$a, bx=x$b, ay=y$a, by=y$b)
+
+ xlim <- ylim <- c(-1,8)
+ plot(NA, xlab="", ylab="", xlim=ylim, ylim=ylim)
+
+ # Highlight the offset coordinates for both lines
+ points(t(x$a), pch="+", col="red")
+ text(t(x$a), label=expression(a[x]), adj=c(-1,0.5))
+ points(t(y$a), pch="+", col="blue")
+ text(t(y$a), label=expression(a[y]), adj=c(-1,0.5))
+
+ v <- c(-1,1)*10
+ xv <- list(x=x$a[1]+x$b[1]*v, y=x$a[2]+x$b[2]*v)
+ yv <- list(x=y$a[1]+y$b[1]*v, y=y$a[2]+y$b[2]*v)
+
+ lines(xv, col="red")
+ lines(yv, col="blue")
+
+ points(t(fit$xs), cex=2.0, col="red")
+ text(t(fit$xs), label=expression(x(s)), adj=c(+2,0.5))
+ points(t(fit$yt), cex=1.5, col="blue")
+ text(t(fit$yt), label=expression(y(t)), adj=c(-1,0.5))
+ print(fit)
+
+
+ ############################################################
+ # Lines in three-dimensions
+ ############################################################
+ x <- list(a=c(0,0,0), b=c(1,1,1)) # The 'diagonal'
+ y <- list(a=c(2,1,2), b=c(2,1,3)) # A 'fitted' line
+ fit <- distanceBetweenLines(ax=x$a, bx=x$b, ay=y$a, by=y$b)
+
+ xlim <- ylim <- zlim <- c(-1,3)
+ dummy <- t(c(1,1,1))*100
+
+ # Coordinates for the lines in 3d
+ v <- seq(-10,10, by=1)
+ xv <- list(x=x$a[1]+x$b[1]*v, y=x$a[2]+x$b[2]*v, z=x$a[3]+x$b[3]*v)
+ yv <- list(x=y$a[1]+y$b[1]*v, y=y$a[2]+y$b[2]*v, z=y$a[3]+y$b[3]*v)
+
+ for (theta in seq(30,140,length.out=3)) {
+ plot3d(dummy, theta=theta, phi=30, xlab="", ylab="", zlab="",
+ xlim=ylim, ylim=ylim, zlim=zlim)
+
+ # Highlight the offset coordinates for both lines
+ points3d(t(x$a), pch="+", col="red")
+ text3d(t(x$a), label=expression(a[x]), adj=c(-1,0.5))
+ points3d(t(y$a), pch="+", col="blue")
+ text3d(t(y$a), label=expression(a[y]), adj=c(-1,0.5))
+
+ # Draw the lines
+ lines3d(xv, col="red")
+ lines3d(yv, col="blue")
+
+ # Draw the two points that are closest to each other
+ points3d(t(fit$xs), cex=2.0, col="red")
+ text3d(t(fit$xs), label=expression(x(s)), adj=c(+2,0.5))
+ points3d(t(fit$yt), cex=1.5, col="blue")
+ text3d(t(fit$yt), label=expression(y(t)), adj=c(-1,0.5))
+
+ # Draw the distance between the two points
+ lines3d(rbind(fit$xs,fit$yt), col="purple", lwd=2)
+ }
+
+ print(fit)
+
+ } # for (zzz in 0)
Loading required package: R.basic
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, :
there is no package called 'R.basic'
> rm(zzz)
>
> proc.time()
user system elapsed
0.306 0.048 0.343
aroma.light.Rcheck/tests/findPeaksAndValleys.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> layout(matrix(1:3, ncol=1))
> par(mar=c(2,4,4,1)+0.1)
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A unimodal distribution
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> x1 <- rnorm(n=10000, mean=0, sd=1)
> x <- x1
> fit <- findPeaksAndValleys(x)
> print(fit)
type x density
1 peak 0.06086266 0.3938079173
2 valley 3.63444591 0.0004188846
3 peak 3.83390172 0.0005647564
> plot(density(x), lwd=2, main="x1")
> abline(v=fit$x)
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> x2 <- rnorm(n=10000, mean=4, sd=1)
> x3 <- rnorm(n=10000, mean=8, sd=1)
> x <- c(x1,x2,x3)
> fit <- findPeaksAndValleys(x)
> print(fit)
type x density
1 peak -0.0610495 0.12363290
2 valley 1.9485093 0.04486347
3 peak 3.9580680 0.12154833
4 valley 6.0028822 0.04289612
5 peak 8.0124410 0.12453937
> plot(density(x), lwd=2, main="c(x1,x2,x3)")
> abline(v=fit$x)
>
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution with clear separation
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> x1b <- rnorm(n=10000, mean=0, sd=0.1)
> x2b <- rnorm(n=10000, mean=4, sd=0.1)
> x3b <- rnorm(n=10000, mean=8, sd=0.1)
> x <- c(x1b,x2b,x3b)
>
> # Illustrating explicit usage of density()
> d <- density(x)
> fit <- findPeaksAndValleys(d, tol=0)
> print(fit)
type x density
1 peak -0.01148341 3.428411e-01
2 valley 1.98167564 1.187180e-06
3 peak 3.97483470 3.422214e-01
4 valley 5.98942558 1.176476e-06
5 peak 7.98258464 3.427964e-01
> plot(d, lwd=2, main="c(x1b,x2b,x3b)")
> abline(v=fit$x)
>
> proc.time()
user system elapsed
0.285 0.056 0.330
aroma.light.Rcheck/tests/fitPrincipalCurve.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # 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: 4222294
Iteration 1---distance^2: 358.2358
Iteration 2---distance^2: 356.1203
Iteration 3---distance^2: 356.146
Converged: TRUE
Number of iterations: 3
Processing time/iteration: 0.1s (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)
>
> proc.time()
user system elapsed
0.928 0.089 1.007
aroma.light.Rcheck/tests/fitXYCurve.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Simulate data from the model y <- a + bx + x^c + eps(bx)
> x <- rexp(1000)
> a <- c(2,15)
> b <- c(2,1)
> c <- c(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(x), mean=0, sd=0.1*x))
> Y <- a + bx + xc + eps
> Y <- t(Y)
>
> lim <- c(0,70)
> plot(Y, xlim=lim, ylim=lim)
>
> # Fit principal curve through a subset of (y_1, y_2)
> subset <- sample(nrow(Y), size=0.3*nrow(Y))
> fit <- fitXYCurve(Y[subset,], bandwidth=0.2)
>
> lines(fit, col="red", lwd=2)
>
> # Backtransform (y_1, y_2) keeping y_1 unchanged
> YN <- backtransformXYCurve(Y, fit=fit)
> points(YN, col="blue")
> abline(a=0, b=1, col="red", lwd=2)
>
> proc.time()
user system elapsed
0.367 0.040 0.396
aroma.light.Rcheck/tests/iwpca.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> for (zzz in 0) {
+
+ # This example requires plot3d() in R.basic [http://www.braju.com/R/]
+ if (!require(pkgName <- "R.basic", character.only=TRUE)) break
+
+ # Simulate data from the model y <- a + bx + eps(bx)
+ x <- rexp(1000)
+ a <- c(2,15,3)
+ b <- c(2,3,4)
+ bx <- outer(b,x)
+ eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(x), mean=0, sd=0.1*x))
+ y <- a + bx + eps
+ y <- t(y)
+
+ # Add some outliers by permuting the dimensions for 1/10 of the observations
+ idx <- sample(1:nrow(y), size=1/10*nrow(y))
+ y[idx,] <- y[idx,c(2,3,1)]
+
+ # Plot the data with fitted lines at four different view points
+ opar <- par(mar=c(1,1,1,1)+0.1)
+ N <- 4
+ layout(matrix(1:N, nrow=2, byrow=TRUE))
+ theta <- seq(0,270,length.out=N)
+ phi <- rep(20, length.out=N)
+ xlim <- ylim <- zlim <- c(0,45)
+ persp <- list()
+ for (kk in seq_along(theta)) {
+ # Plot the data
+ persp[[kk]] <- plot3d(y, theta=theta[kk], phi=phi[kk], xlim=xlim, ylim=ylim, zlim=zlim)
+ }
+
+ # Weights on the observations
+ # Example a: Equal weights
+ w <- NULL
+ # Example b: More weight on the outliers (uncomment to test)
+ w <- rep(1, length(x)); w[idx] <- 0.8
+
+ # ...and show all iterations too with different colors.
+ maxIter <- c(seq(1,20,length.out=10),Inf)
+ col <- topo.colors(length(maxIter))
+ # Show the fitted value for every iteration
+ for (ii in seq_along(maxIter)) {
+ # Fit a line using IWPCA through data
+ fit <- iwpca(y, w=w, maxIter=maxIter[ii], swapDirections=TRUE)
+
+ ymid <- fit$xMean
+ d0 <- apply(y, MARGIN=2, FUN=min) - ymid
+ d1 <- apply(y, MARGIN=2, FUN=max) - ymid
+ b <- fit$vt[1,]
+ y0 <- -b * max(abs(d0))
+ y1 <- b * max(abs(d1))
+ yline <- matrix(c(y0,y1), nrow=length(b), ncol=2)
+ yline <- yline + ymid
+
+ for (kk in seq_along(theta)) {
+ # Set pane to draw in
+ par(mfg=c((kk-1) %/% 2, (kk-1) %% 2) + 1)
+ # Set the viewpoint of the pane
+ options(persp.matrix=persp[[kk]])
+
+ # Get the first principal component
+ points3d(t(ymid), col=col[ii])
+ lines3d(t(yline), col=col[ii])
+
+ # Highlight the last one
+ if (ii == length(maxIter))
+ lines3d(t(yline), col="red", lwd=3)
+ }
+ }
+
+ par(opar)
+
+ } # for (zzz in 0)
Loading required package: R.basic
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, :
there is no package called 'R.basic'
> rm(zzz)
>
> proc.time()
user system elapsed
0.299 0.041 0.329
aroma.light.Rcheck/tests/likelihood.smooth.spline.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Define f(x)
> f <- expression(0.1*x^4 + 1*x^3 + 2*x^2 + x + 10*sin(2*x))
>
> # Simulate data from this function in the range [a,b]
> a <- -2; b <- 5
> x <- seq(a, b, length.out=3000)
> y <- eval(f)
>
> # Add some noise to the data
> y <- y + rnorm(length(y), 0, 10)
>
> # Plot the function and its second derivative
> plot(x,y, type="l", lwd=4)
>
> # Fit a cubic smoothing spline and plot it
> g <- smooth.spline(x,y, df=16)
> lines(g, col="yellow", lwd=2, lty=2)
>
> # Calculating the (log) likelihood of the fitted spline
> l <- likelihood(g)
>
> cat("Log likelihood with unique x values:\n")
Log likelihood with unique x values:
> print(l)
Likelihood of smoothing spline: -298438.8
Log base: 2.718282
Weighted residuals sum of square: 298439
Penalty: -0.1347938
Smoothing parameter lambda: 0.0009257147
Roughness score: 145.6105
>
> # Note that this is not the same as the log likelihood of the
> # data on the fitted spline iff the x values are non-unique
> x[1:5] <- x[1] # Non-unique x values
> g <- smooth.spline(x,y, df=16)
> l <- likelihood(g)
>
> cat("\nLog likelihood of the *spline* data set:\n")
Log likelihood of the *spline* data set:
> print(l)
Likelihood of smoothing spline: -298011
Log base: 2.718282
Weighted residuals sum of square: 298011.1
Penalty: -0.1349871
Smoothing parameter lambda: 0.0009261969
Roughness score: 145.7435
>
> # In cases with non unique x values one has to proceed as
> # below if one want to get the log likelihood for the original
> # data.
> l <- likelihood(g, x=x, y=y)
> cat("\nLog likelihood of the *original* data set:\n")
Log likelihood of the *original* data set:
> print(l)
Likelihood of smoothing spline: -298423.4
Log base: 2.718282
Weighted residuals sum of square: 298423.5
Penalty: -0.1349873
Smoothing parameter lambda: 0.0009261969
Roughness score: 145.7436
>
>
>
>
>
>
> proc.time()
user system elapsed
0.324 0.062 0.376
aroma.light.Rcheck/tests/medianPolish.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Deaths from sport parachuting; from ABC of EDA, p.224:
> deaths <- matrix(c(14,15,14, 7,4,7, 8,2,10, 15,9,10, 0,2,0), ncol=3, byrow=TRUE)
> rownames(deaths) <- c("1-24", "25-74", "75-199", "200++", "NA")
> colnames(deaths) <- 1973:1975
>
> print(deaths)
1973 1974 1975
1-24 14 15 14
25-74 7 4 7
75-199 8 2 10
200++ 15 9 10
NA 0 2 0
>
> mp <- medianPolish(deaths)
> mp1 <- medpolish(deaths, trace=FALSE)
> print(mp)
Median Polish Results (Dataset: "deaths")
Overall: 8
Row Effects:
1-24 25-74 75-199 200++ NA
6 -1 0 2 -8
Column Effects:
1973 1974 1975
0 -1 0
Residuals:
1973 1974 1975
1-24 0 2 0
25-74 0 -2 0
75-199 0 -5 2
200++ 5 0 0
NA 0 3 0
>
> ff <- c("overall", "row", "col", "residuals")
> stopifnot(all.equal(mp[ff], mp1[ff]))
>
> # Validate decomposition:
> stopifnot(all.equal(deaths, mp$overall+outer(mp$row,mp$col,"+")+mp$resid))
>
> proc.time()
user system elapsed
0.246 0.027 0.263
aroma.light.Rcheck/tests/normalizeAffine.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> pathname <- system.file("data-ex", "PMT-RGData.dat", package="aroma.light")
> rg <- read.table(pathname, header=TRUE, sep="\t")
> nbrOfScans <- max(rg$slide)
>
> rg <- as.list(rg)
> for (field in c("R", "G"))
+ rg[[field]] <- matrix(as.double(rg[[field]]), ncol=nbrOfScans)
> rg$slide <- rg$spot <- NULL
> rg <- as.matrix(as.data.frame(rg))
> colnames(rg) <- rep(c("R", "G"), each=nbrOfScans)
>
> rgC <- rg
>
> layout(matrix(c(1,2,0,3,4,0,5,6,7), ncol=3, byrow=TRUE))
>
> for (channel in c("R", "G")) {
+ sidx <- which(colnames(rg) == channel)
+ channelColor <- switch(channel, R="red", G="green")
+
+ # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ # The raw data
+ # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ plotMvsAPairs(rg, channel=channel)
+ title(main=paste("Observed", channel))
+ box(col=channelColor)
+
+ # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ # The calibrated data
+ # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ rgC[,sidx] <- calibrateMultiscan(rg[,sidx], average=NULL)
+
+ plotMvsAPairs(rgC, channel=channel)
+ title(main=paste("Calibrated", channel))
+ box(col=channelColor)
+ } # for (channel ...)
There were 50 or more warnings (use warnings() to see the first 50)
>
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # The average calibrated data
> #
> # Note how the red signals are weaker than the green. The reason
> # for this can be that the scale factor in the green channel is
> # greater than in the red channel, but it can also be that there
> # is a remaining relative difference in bias between the green
> # and the red channel, a bias that precedes the scanning.
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> rgCA <- matrix(NA_real_, nrow=nrow(rg), ncol=2)
> colnames(rgCA) <- c("R", "G")
> for (channel in c("R", "G")) {
+ sidx <- which(colnames(rg) == channel)
+ rgCA[,channel] <- calibrateMultiscan(rg[,sidx])
+ }
>
> plotMvsA(rgCA)
> title(main="Average calibrated")
>
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # The affine normalized average calibrated data
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Create a matrix where the columns represent the channels
> # to be normalized.
> rgCAN <- rgCA
> # Affine normalization of channels
> rgCAN <- normalizeAffine(rgCAN)
>
> plotMvsA(rgCAN)
> title(main="Affine normalized A.C.")
>
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # It is always ok to rescale the affine normalized data if its
> # done on (R,G); not on (A,M)! However, this is only needed for
> # esthetic purposes.
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> rgCAN <- rgCAN * 2^5
> plotMvsA(rgCAN)
> title(main="Rescaled normalized")
>
>
>
> proc.time()
user system elapsed
2.090 0.050 2.129
aroma.light.Rcheck/tests/normalizeAverage.list.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Simulate ten samples of different lengths
> N <- 10000
> X <- list()
> for (kk in 1:8) {
+ rfcn <- list(rnorm, rgamma)[[sample(2, size=1)]]
+ size <- runif(1, min=0.3, max=1)
+ a <- rgamma(1, shape=20, rate=10)
+ b <- rgamma(1, shape=10, rate=10)
+ values <- rfcn(size*N, a, b)
+
+ # "Censor" values
+ values[values < 0 | values > 8] <- NA_real_
+
+ X[[kk]] <- values
+ }
>
> # Add 20% missing values
> X <- lapply(X, FUN=function(x) {
+ x[sample(length(x), size=0.20*length(x))] <- NA_real_
+ x
+ })
>
> # Normalize quantiles
> Xn <- normalizeAverage(X, na.rm=TRUE, targetAvg=median(unlist(X), na.rm=TRUE))
>
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, Xn, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The normalized distributions")
>
> proc.time()
user system elapsed
0.321 0.063 0.375
aroma.light.Rcheck/tests/normalizeAverage.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Simulate three samples with on average 20% missing values
> N <- 10000
> X <- cbind(rnorm(N, mean=3, sd=1),
+ rnorm(N, mean=4, sd=2),
+ rgamma(N, shape=2, rate=1))
> X[sample(3*N, size=0.20*3*N)] <- NA_real_
>
> # Normalize quantiles
> Xn <- normalizeAverage(X, na.rm=TRUE, targetAvg=median(X, na.rm=TRUE))
>
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, Xn, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The three original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The three normalized distributions")
>
> proc.time()
user system elapsed
0.307 0.024 0.322
aroma.light.Rcheck/tests/normalizeCurveFit.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> pathname <- system.file("data-ex", "PMT-RGData.dat", package="aroma.light")
> rg <- read.table(pathname, header=TRUE, sep="\t")
> nbrOfScans <- max(rg$slide)
>
> rg <- as.list(rg)
> for (field in c("R", "G"))
+ rg[[field]] <- matrix(as.double(rg[[field]]), ncol=nbrOfScans)
> rg$slide <- rg$spot <- NULL
> rg <- as.matrix(as.data.frame(rg))
> colnames(rg) <- rep(c("R", "G"), each=nbrOfScans)
>
> layout(matrix(c(1,2,0,3,4,0,5,6,7), ncol=3, byrow=TRUE))
>
> rgC <- rg
> for (channel in c("R", "G")) {
+ sidx <- which(colnames(rg) == channel)
+ channelColor <- switch(channel, R="red", G="green")
+
+ # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ # The raw data
+ # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ plotMvsAPairs(rg[,sidx])
+ title(main=paste("Observed", channel))
+ box(col=channelColor)
+
+ # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ # The calibrated data
+ # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ rgC[,sidx] <- calibrateMultiscan(rg[,sidx], average=NULL)
+
+ plotMvsAPairs(rgC[,sidx])
+ title(main=paste("Calibrated", channel))
+ box(col=channelColor)
+ } # for (channel ...)
>
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # The average calibrated data
> #
> # Note how the red signals are weaker than the green. The reason
> # for this can be that the scale factor in the green channel is
> # greater than in the red channel, but it can also be that there
> # is a remaining relative difference in bias between the green
> # and the red channel, a bias that precedes the scanning.
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> rgCA <- rg
> for (channel in c("R", "G")) {
+ sidx <- which(colnames(rg) == channel)
+ rgCA[,sidx] <- calibrateMultiscan(rg[,sidx])
+ }
>
> rgCAavg <- matrix(NA_real_, nrow=nrow(rgCA), ncol=2)
> colnames(rgCAavg) <- c("R", "G")
> for (channel in c("R", "G")) {
+ sidx <- which(colnames(rg) == channel)
+ rgCAavg[,channel] <- apply(rgCA[,sidx], MARGIN=1, FUN=median, na.rm=TRUE)
+ }
>
> # Add some "fake" outliers
> outliers <- 1:600
> rgCAavg[outliers,"G"] <- 50000
>
> plotMvsA(rgCAavg)
> title(main="Average calibrated (AC)")
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Normalize data
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Weight-down outliers when normalizing
> weights <- rep(1, nrow(rgCAavg))
> weights[outliers] <- 0.001
>
> # Affine normalization of channels
> rgCANa <- normalizeAffine(rgCAavg, weights=weights)
> # It is always ok to rescale the affine normalized data if its
> # done on (R,G); not on (A,M)! However, this is only needed for
> # esthetic purposes.
> rgCANa <- rgCANa *2^1.4
> plotMvsA(rgCANa)
> title(main="Normalized AC")
>
> # Curve-fit (lowess) normalization
> rgCANlw <- normalizeLowess(rgCAavg, weights=weights)
Warning message:
In normalizeCurveFit.matrix(X, method = "lowess", ...) :
Weights were rounded to {0,1} since 'lowess' normalization supports only zero-one weights.
> plotMvsA(rgCANlw, col="orange", add=TRUE)
>
> # Curve-fit (loess) normalization
> rgCANl <- normalizeLoess(rgCAavg, weights=weights)
> plotMvsA(rgCANl, col="red", add=TRUE)
>
> # Curve-fit (robust spline) normalization
> rgCANrs <- normalizeRobustSpline(rgCAavg, weights=weights)
> plotMvsA(rgCANrs, col="blue", add=TRUE)
>
> legend(x=0,y=16, legend=c("affine", "lowess", "loess", "r. spline"), pch=19,
+ col=c("black", "orange", "red", "blue"), ncol=2, x.intersp=0.3, bty="n")
>
>
> plotMvsMPairs(cbind(rgCANa, rgCANlw), col="orange", xlab=expression(M[affine]))
> title(main="Normalized AC")
> plotMvsMPairs(cbind(rgCANa, rgCANl), col="red", add=TRUE)
> plotMvsMPairs(cbind(rgCANa, rgCANrs), col="blue", add=TRUE)
> abline(a=0, b=1, lty=2)
> legend(x=-6,y=6, legend=c("lowess", "loess", "r. spline"), pch=19,
+ col=c("orange", "red", "blue"), ncol=2, x.intersp=0.3, bty="n")
>
>
> proc.time()
user system elapsed
6.200 0.072 6.260
aroma.light.Rcheck/tests/normalizeDifferencesToAverage.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Simulate three shifted tracks of different lengths with same profiles
> ns <- c(A=2, B=1, C=0.25)*1000
> xx <- lapply(ns, FUN=function(n) { seq(from=1, to=max(ns), length.out=n) })
> zz <- mapply(seq_along(ns), ns, FUN=function(z,n) rep(z,n))
>
> yy <- list(
+ A = rnorm(ns["A"], mean=0, sd=0.5),
+ B = rnorm(ns["B"], mean=5, sd=0.4),
+ C = rnorm(ns["C"], mean=-5, sd=1.1)
+ )
> yy <- lapply(yy, FUN=function(y) {
+ n <- length(y)
+ y[1:(n/2)] <- y[1:(n/2)] + 2
+ y[1:(n/4)] <- y[1:(n/4)] - 4
+ y
+ })
>
> # Shift all tracks toward the first track
> yyN <- normalizeDifferencesToAverage(yy, baseline=1)
>
> # The baseline channel is not changed
> stopifnot(identical(yy[[1]], yyN[[1]]))
>
> # Get the estimated parameters
> fit <- attr(yyN, "fit")
>
> # Plot the tracks
> layout(matrix(1:2, ncol=1))
> x <- unlist(xx)
> col <- unlist(zz)
> y <- unlist(yy)
> yN <- unlist(yyN)
> plot(x, y, col=col, ylim=c(-10,10))
> plot(x, yN, col=col, ylim=c(-10,10))
>
> proc.time()
user system elapsed
0.410 0.045 0.443
aroma.light.Rcheck/tests/normalizeFragmentLength-ex1.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Example 1: Single-enzyme fragment-length normalization of 6 arrays
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Number samples
> I <- 9
>
> # Number of loci
> J <- 1000
>
> # Fragment lengths
> fl <- seq(from=100, to=1000, length.out=J)
>
> # Simulate data points with unknown fragment lengths
> hasUnknownFL <- seq(from=1, to=J, by=50)
> fl[hasUnknownFL] <- NA_real_
>
> # Simulate data
> y <- matrix(0, nrow=J, ncol=I)
> maxY <- 12
> for (kk in 1:I) {
+ k <- runif(n=1, min=3, max=5)
+ mu <- function(fl) {
+ mu <- rep(maxY, length(fl))
+ ok <- !is.na(fl)
+ mu[ok] <- mu[ok] - fl[ok]^{1/k}
+ mu
+ }
+ eps <- rnorm(J, mean=0, sd=1)
+ y[,kk] <- mu(fl) + eps
+ }
>
> # Normalize data (to a zero baseline)
> yN <- apply(y, MARGIN=2, FUN=function(y) {
+ normalizeFragmentLength(y, fragmentLengths=fl, onMissing="median")
+ })
>
> # The correction factors
> rho <- y-yN
> print(summary(rho))
V1 V2 V3 V4
Min. :6.650 Min. :7.058 Min. :6.142 Min. :8.007
1st Qu.:7.077 1st Qu.:7.557 1st Qu.:6.505 1st Qu.:8.106
Median :7.440 Median :7.983 Median :6.963 Median :8.299
Mean :7.572 Mean :8.007 Mean :7.078 Mean :8.459
3rd Qu.:8.028 3rd Qu.:8.437 3rd Qu.:7.603 3rd Qu.:8.776
Max. :8.935 Max. :9.080 Max. :8.445 Max. :9.398
V5 V6 V7 V8
Min. :2.418 Min. :2.536 Min. :6.659 Min. :6.239
1st Qu.:3.300 1st Qu.:3.430 1st Qu.:7.047 1st Qu.:6.573
Median :4.233 Median :4.318 Median :7.426 Median :7.025
Mean :4.467 Mean :4.517 Mean :7.498 Mean :7.137
3rd Qu.:5.573 3rd Qu.:5.550 3rd Qu.:7.941 3rd Qu.:7.650
Max. :7.249 Max. :7.163 Max. :8.533 Max. :8.456
V9
Min. :5.139
1st Qu.:5.663
Median :6.215
Mean :6.363
3rd Qu.:7.016
Max. :8.092
> # The correction for units with unknown fragment lengths
> # equals the median correction factor of all other units
> print(summary(rho[hasUnknownFL,]))
V1 V2 V3 V4 V5
Min. :7.44 Min. :7.983 Min. :6.963 Min. :8.299 Min. :4.233
1st Qu.:7.44 1st Qu.:7.983 1st Qu.:6.963 1st Qu.:8.299 1st Qu.:4.233
Median :7.44 Median :7.983 Median :6.963 Median :8.299 Median :4.233
Mean :7.44 Mean :7.983 Mean :6.963 Mean :8.299 Mean :4.233
3rd Qu.:7.44 3rd Qu.:7.983 3rd Qu.:6.963 3rd Qu.:8.299 3rd Qu.:4.233
Max. :7.44 Max. :7.983 Max. :6.963 Max. :8.299 Max. :4.233
V6 V7 V8 V9
Min. :4.318 Min. :7.426 Min. :7.025 Min. :6.215
1st Qu.:4.318 1st Qu.:7.426 1st Qu.:7.025 1st Qu.:6.215
Median :4.318 Median :7.426 Median :7.025 Median :6.215
Mean :4.318 Mean :7.426 Mean :7.025 Mean :6.215
3rd Qu.:4.318 3rd Qu.:7.426 3rd Qu.:7.025 3rd Qu.:6.215
Max. :4.318 Max. :7.426 Max. :7.025 Max. :6.215
>
> # Plot raw data
> layout(matrix(1:9, ncol=3))
> xlim <- c(0,max(fl, na.rm=TRUE))
> ylim <- c(0,max(y, na.rm=TRUE))
> xlab <- "Fragment length"
> ylab <- expression(log2(theta))
> for (kk in 1:I) {
+ plot(fl, y[,kk], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
+ ok <- (is.finite(fl) & is.finite(y[,kk]))
+ lines(lowess(fl[ok], y[ok,kk]), col="red", lwd=2)
+ }
>
> # Plot normalized data
> layout(matrix(1:9, ncol=3))
> ylim <- c(-1,1)*max(y, na.rm=TRUE)/2
> for (kk in 1:I) {
+ plot(fl, yN[,kk], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
+ ok <- (is.finite(fl) & is.finite(y[,kk]))
+ lines(lowess(fl[ok], yN[ok,kk]), col="blue", lwd=2)
+ }
>
> proc.time()
user system elapsed
0.892 0.061 0.940
aroma.light.Rcheck/tests/normalizeFragmentLength-ex2.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Example 2: Two-enzyme fragment-length normalization of 6 arrays
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> set.seed(0xbeef)
>
> # Number samples
> I <- 5
>
> # Number of loci
> J <- 3000
>
> # Fragment lengths (two enzymes)
> fl <- matrix(0, nrow=J, ncol=2)
> fl[,1] <- seq(from=100, to=1000, length.out=J)
> fl[,2] <- seq(from=1000, to=100, length.out=J)
>
> # Let 1/2 of the units be on both enzymes
> fl[seq(from=1, to=J, by=4),1] <- NA_real_
> fl[seq(from=2, to=J, by=4),2] <- NA_real_
>
> # Let some have unknown fragment lengths
> hasUnknownFL <- seq(from=1, to=J, by=15)
> fl[hasUnknownFL,] <- NA_real_
>
> # Sty/Nsp mixing proportions:
> rho <- rep(1, I)
> rho[1] <- 1/3; # Less Sty in 1st sample
> rho[3] <- 3/2; # More Sty in 3rd sample
>
>
> # Simulate data
> z <- array(0, dim=c(J,2,I))
> maxLog2Theta <- 12
> for (ii in 1:I) {
+ # Common effect for both enzymes
+ mu <- function(fl) {
+ k <- runif(n=1, min=3, max=5)
+ mu <- rep(maxLog2Theta, length(fl))
+ ok <- is.finite(fl)
+ mu[ok] <- mu[ok] - fl[ok]^{1/k}
+ mu
+ }
+
+ # Calculate the effect for each data point
+ for (ee in 1:2) {
+ z[,ee,ii] <- mu(fl[,ee])
+ }
+
+ # Update the Sty/Nsp mixing proportions
+ ee <- 2
+ z[,ee,ii] <- rho[ii]*z[,ee,ii]
+
+ # Add random errors
+ for (ee in 1:2) {
+ eps <- rnorm(J, mean=0, sd=1/sqrt(2))
+ z[,ee,ii] <- z[,ee,ii] + eps
+ }
+ }
>
>
> hasFl <- is.finite(fl)
>
> unitSets <- list(
+ nsp = which( hasFl[,1] & !hasFl[,2]),
+ sty = which(!hasFl[,1] & hasFl[,2]),
+ both = which( hasFl[,1] & hasFl[,2]),
+ none = which(!hasFl[,1] & !hasFl[,2])
+ )
>
> # The observed data is a mix of two enzymes
> theta <- matrix(NA_real_, nrow=J, ncol=I)
>
> # Single-enzyme units
> for (ee in 1:2) {
+ uu <- unitSets[[ee]]
+ theta[uu,] <- 2^z[uu,ee,]
+ }
>
> # Both-enzyme units (sum on intensity scale)
> uu <- unitSets$both
> theta[uu,] <- (2^z[uu,1,]+2^z[uu,2,])/2
>
> # Missing units (sample from the others)
> uu <- unitSets$none
> theta[uu,] <- apply(theta, MARGIN=2, sample, size=length(uu))
>
> # Calculate target array
> thetaT <- rowMeans(theta, na.rm=TRUE)
> targetFcns <- list()
> for (ee in 1:2) {
+ uu <- unitSets[[ee]]
+ fit <- lowess(fl[uu,ee], log2(thetaT[uu]))
+ class(fit) <- "lowess"
+ targetFcns[[ee]] <- function(fl, ...) {
+ predict(fit, newdata=fl)
+ }
+ }
>
>
> # Fit model only to a subset of the data
> subsetToFit <- setdiff(1:J, seq(from=1, to=J, by=10))
>
> # Normalize data (to a target baseline)
> thetaN <- matrix(NA_real_, nrow=J, ncol=I)
> fits <- vector("list", I)
> for (ii in 1:I) {
+ lthetaNi <- normalizeFragmentLength(log2(theta[,ii]), targetFcns=targetFcns,
+ fragmentLengths=fl, onMissing="median",
+ subsetToFit=subsetToFit, .returnFit=TRUE)
+ fits[[ii]] <- attr(lthetaNi, "modelFit")
+ thetaN[,ii] <- 2^lthetaNi
+ }
>
>
> # Plot raw data
> xlim <- c(0, max(fl, na.rm=TRUE))
> ylim <- c(0, max(log2(theta), na.rm=TRUE))
> Mlim <- c(-1,1)*4
> xlab <- "Fragment length"
> ylab <- expression(log2(theta))
> Mlab <- expression(M==log[2](theta/theta[R]))
>
> layout(matrix(1:(3*I), ncol=I, byrow=TRUE))
> for (ii in 1:I) {
+ plot(NA, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, main="raw")
+
+ # Single-enzyme units
+ for (ee in 1:2) {
+ # The raw data
+ uu <- unitSets[[ee]]
+ points(fl[uu,ee], log2(theta[uu,ii]), col=ee+1)
+ }
+
+ # Both-enzyme units (use fragment-length for enzyme #1)
+ uu <- unitSets$both
+ points(fl[uu,1], log2(theta[uu,ii]), col=3+1)
+
+ for (ee in 1:2) {
+ # The true effects
+ uu <- unitSets[[ee]]
+ lines(lowess(fl[uu,ee], log2(theta[uu,ii])), col="black", lwd=4, lty=3)
+
+ # The estimated effects
+ fit <- fits[[ii]][[ee]]$fit
+ lines(fit, col="orange", lwd=3)
+
+ muT <- targetFcns[[ee]](fl[uu,ee])
+ lines(fl[uu,ee], muT, col="cyan", lwd=1)
+ }
+ }
>
> # Calculate log-ratios
> thetaR <- rowMeans(thetaN, na.rm=TRUE)
> M <- log2(thetaN/thetaR)
>
> # Plot normalized data
> for (ii in 1:I) {
+ plot(NA, xlim=xlim, ylim=Mlim, xlab=xlab, ylab=Mlab, main="normalized")
+ # Single-enzyme units
+ for (ee in 1:2) {
+ # The normalized data
+ uu <- unitSets[[ee]]
+ points(fl[uu,ee], M[uu,ii], col=ee+1)
+ }
+ # Both-enzyme units (use fragment-length for enzyme #1)
+ uu <- unitSets$both
+ points(fl[uu,1], M[uu,ii], col=3+1)
+ }
>
> ylim <- c(0,1.5)
> for (ii in 1:I) {
+ data <- list()
+ for (ee in 1:2) {
+ # The normalized data
+ uu <- unitSets[[ee]]
+ data[[ee]] <- M[uu,ii]
+ }
+ uu <- unitSets$both
+ if (length(uu) > 0)
+ data[[3]] <- M[uu,ii]
+
+ uu <- unitSets$none
+ if (length(uu) > 0)
+ data[[4]] <- M[uu,ii]
+
+ cols <- seq_along(data)+1
+ plotDensity(data, col=cols, xlim=Mlim, xlab=Mlab, main="normalized")
+
+ abline(v=0, lty=2)
+ }
>
>
> proc.time()
user system elapsed
0.852 0.086 0.927
aroma.light.Rcheck/tests/normalizeQuantileRank.list.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Simulate ten samples of different lengths
> N <- 10000
> X <- list()
> for (kk in 1:8) {
+ rfcn <- list(rnorm, rgamma)[[sample(2, size=1)]]
+ size <- runif(1, min=0.3, max=1)
+ a <- rgamma(1, shape=20, rate=10)
+ b <- rgamma(1, shape=10, rate=10)
+ values <- rfcn(size*N, a, b)
+
+ # "Censor" values
+ values[values < 0 | values > 8] <- NA_real_
+
+ X[[kk]] <- values
+ }
>
> # Add 20% missing values
> X <- lapply(X, FUN=function(x) {
+ x[sample(length(x), size=0.20*length(x))] <- NA_real_
+ x
+ })
>
> # Normalize quantiles
> Xn <- normalizeQuantile(X)
>
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The normalized distributions")
>
> proc.time()
user system elapsed
0.375 0.052 0.416
aroma.light.Rcheck/tests/normalizeQuantileRank.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Simulate three samples with on average 20% missing values
> N <- 10000
> X <- cbind(rnorm(N, mean=3, sd=1),
+ rnorm(N, mean=4, sd=2),
+ rgamma(N, shape=2, rate=1))
> X[sample(3*N, size=0.20*3*N)] <- NA_real_
>
> # Normalize quantiles
> Xn <- normalizeQuantile(X)
>
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, Xn, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The three original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The three normalized distributions")
>
> proc.time()
user system elapsed
0.302 0.043 0.334
aroma.light.Rcheck/tests/normalizeQuantileSpline.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Simulate three samples with on average 20% missing values
> N <- 10000
> X <- cbind(rnorm(N, mean=3, sd=1),
+ rnorm(N, mean=4, sd=2),
+ rgamma(N, shape=2, rate=1))
> X[sample(3*N, size=0.20*3*N)] <- NA_real_
>
> # Plot the data
> layout(matrix(c(1,0,2:5), ncol=2, byrow=TRUE))
> xlim <- range(X, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The three original distributions")
>
> Xn <- normalizeQuantile(X)
> plotDensity(Xn, lwd=2, xlim=xlim, main="The three normalized distributions")
> plotXYCurve(X, Xn, xlim=xlim, main="The three normalized distributions")
>
> Xn2 <- normalizeQuantileSpline(X, xTarget=Xn[,1], spar=0.99)
> plotDensity(Xn2, lwd=2, xlim=xlim, main="The three normalized distributions")
> plotXYCurve(X, Xn2, xlim=xlim, main="The three normalized distributions")
>
> proc.time()
user system elapsed
0.711 0.067 0.766
aroma.light.Rcheck/tests/normalizeTumorBoost,flavors.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
> library("R.utils")
Loading required package: R.oo
Loading required package: R.methodsS3
R.methodsS3 v1.8.2 (2022-06-13 22:00:14 UTC) successfully loaded. See ?R.methodsS3 for help.
R.oo v1.25.0 (2022-06-12 02:20:02 UTC) successfully loaded. See ?R.oo for help.
Attaching package: 'R.oo'
The following object is masked from 'package:R.methodsS3':
throw
The following objects are masked from 'package:methods':
getClasses, getMethods
The following objects are masked from 'package:base':
attach, detach, load, save
R.utils v2.12.2 (2022-11-11 22:00:03 UTC) successfully loaded. See ?R.utils for help.
Attaching package: 'R.utils'
The following object is masked from 'package:utils':
timestamp
The following objects are masked from 'package:base':
cat, commandArgs, getOption, isOpen, nullfile, parse, warnings
>
> # Load data
> pathname <- system.file("data-ex/TumorBoost,fracB,exampleData.Rbin", package="aroma.light")
> data <- loadObject(pathname)
>
> # Drop loci with missing values
> data <- na.omit(data)
>
> attachLocally(data)
> pos <- position/1e6
>
> # Call naive genotypes
> muN <- callNaiveGenotypes(betaN)
>
> # Genotype classes
> isAA <- (muN == 0)
> isAB <- (muN == 1/2)
> isBB <- (muN == 1)
>
> # Sanity checks
> stopifnot(all(muN[isAA] == 0))
> stopifnot(all(muN[isAB] == 1/2))
> stopifnot(all(muN[isBB] == 1))
>
> # TumorBoost normalization with different flavors
> betaTNs <- list()
> for (flavor in c("v1", "v2", "v3", "v4")) {
+ betaTN <- normalizeTumorBoost(betaT=betaT, betaN=betaN, preserveScale=FALSE, flavor=flavor)
+
+ # Assert that no non-finite values are introduced
+ stopifnot(all(is.finite(betaTN)))
+
+ # Assert that nothing is flipped
+ stopifnot(all(betaTN[isAA] < 1/2))
+ stopifnot(all(betaTN[isBB] > 1/2))
+
+ betaTNs[[flavor]] <- betaTN
+ }
>
> # Plot
> layout(matrix(1:4, ncol=1))
> par(mar=c(2.5,4,0.5,1)+0.1)
> ylim <- c(-0.05, 1.05)
> col <- rep("#999999", length(muN))
> col[muN == 1/2] <- "#000000"
> for (flavor in names(betaTNs)) {
+ betaTN <- betaTNs[[flavor]]
+ ylab <- sprintf("betaTN[%s]", flavor)
+ plot(pos, betaTN, col=col, ylim=ylim, ylab=ylab)
+ }
>
> proc.time()
user system elapsed
0.620 0.048 0.656
aroma.light.Rcheck/tests/normalizeTumorBoost.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
> library("R.utils")
Loading required package: R.oo
Loading required package: R.methodsS3
R.methodsS3 v1.8.2 (2022-06-13 22:00:14 UTC) successfully loaded. See ?R.methodsS3 for help.
R.oo v1.25.0 (2022-06-12 02:20:02 UTC) successfully loaded. See ?R.oo for help.
Attaching package: 'R.oo'
The following object is masked from 'package:R.methodsS3':
throw
The following objects are masked from 'package:methods':
getClasses, getMethods
The following objects are masked from 'package:base':
attach, detach, load, save
R.utils v2.12.2 (2022-11-11 22:00:03 UTC) successfully loaded. See ?R.utils for help.
Attaching package: 'R.utils'
The following object is masked from 'package:utils':
timestamp
The following objects are masked from 'package:base':
cat, commandArgs, getOption, isOpen, nullfile, parse, warnings
>
> # Load data
> pathname <- system.file("data-ex/TumorBoost,fracB,exampleData.Rbin", package="aroma.light")
> data <- loadObject(pathname)
> attachLocally(data)
> pos <- position/1e6
> muN <- genotypeN
>
> layout(matrix(1:4, ncol=1))
> par(mar=c(2.5,4,0.5,1)+0.1)
> ylim <- c(-0.05, 1.05)
> col <- rep("#999999", length(muN))
> col[muN == 1/2] <- "#000000"
>
> # Allele B fractions for the normal sample
> plot(pos, betaN, col=col, ylim=ylim)
>
> # Allele B fractions for the tumor sample
> plot(pos, betaT, col=col, ylim=ylim)
>
> # TumorBoost w/ naive genotype calls
> betaTN <- normalizeTumorBoost(betaT=betaT, betaN=betaN, preserveScale=FALSE)
> plot(pos, betaTN, col=col, ylim=ylim)
>
> # TumorBoost w/ external multi-sample genotype calls
> betaTNx <- normalizeTumorBoost(betaT=betaT, betaN=betaN, muN=muN, preserveScale=FALSE)
> plot(pos, betaTNx, col=col, ylim=ylim)
>
> proc.time()
user system elapsed
0.549 0.072 0.609
aroma.light.Rcheck/tests/robustSmoothSpline.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> data(cars)
> attach(cars)
> plot(speed, dist, main = "data(cars) & robust smoothing splines")
>
> # Fit a smoothing spline using L_2 norm
> cars.spl <- smooth.spline(speed, dist)
> lines(cars.spl, col = "blue")
>
> # Fit a smoothing spline using L_1 norm
> cars.rspl <- robustSmoothSpline(speed, dist)
> lines(cars.rspl, col = "red")
>
> # Fit a smoothing spline using L_2 norm with 10 degrees of freedom
> lines(smooth.spline(speed, dist, df=10), lty=2, col = "blue")
>
> # Fit a smoothing spline using L_1 norm with 10 degrees of freedom
> lines(robustSmoothSpline(speed, dist, df=10), lty=2, col = "red")
>
> # Fit a smoothing spline using Tukey's biweight norm
> cars.rspl <- robustSmoothSpline(speed, dist, method = "symmetric")
> lines(cars.rspl, col = "purple")
>
> legend(5,120, c(
+ paste("smooth.spline [C.V.] => df =",round(cars.spl$df,1)),
+ paste("robustSmoothSpline L1 [C.V.] => df =",round(cars.rspl$df,1)),
+ paste("robustSmoothSpline symmetric [C.V.] => df =",round(cars.rspl$df,1)),
+ "standard with s( * , df = 10)", "robust with s( * , df = 10)"
+ ),
+ col = c("blue","red","purple","blue","red"), lty = c(1,1,1,2,2),
+ bg='bisque')
>
> proc.time()
user system elapsed
0.340 0.049 0.374
aroma.light.Rcheck/tests/rowAverages.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> X <- matrix(1:30, nrow=5L, ncol=6L)
> mu <- rowMeans(X)
> sd <- apply(X, MARGIN=1L, FUN=sd)
>
> y <- rowAverages(X)
> stopifnot(all(y == mu))
> stopifnot(all(attr(y,"deviance") == sd))
> stopifnot(all(attr(y,"df") == ncol(X)))
>
> proc.time()
user system elapsed
0.247 0.031 0.267
aroma.light.Rcheck/tests/sampleCorrelations.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # Simulate 20000 genes with 10 observations each
> X <- matrix(rnorm(n=20000), ncol=10)
>
> # Calculate the correlation for 5000 random gene pairs
> cor <- sampleCorrelations(X, npairs=5000)
> print(summary(cor))
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.869660 -0.240740 -0.005040 -0.002532 0.235558 0.905074
>
>
> proc.time()
user system elapsed
0.432 0.048 0.468
aroma.light.Rcheck/tests/sampleTuples.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> pairs <- sampleTuples(1:10, size=5, length=2)
> print(pairs)
[,1] [,2]
[1,] 3 10
[2,] 1 5
[3,] 8 3
[4,] 6 5
[5,] 5 2
>
> triples <- sampleTuples(1:10, size=5, length=3)
> print(triples)
[,1] [,2] [,3]
[1,] 9 3 10
[2,] 7 1 8
[3,] 7 5 9
[4,] 5 9 2
[5,] 10 3 2
>
> # Allow tuples with repeated elements
> quadruples <- sampleTuples(1:3, size=5, length=4, replace=TRUE)
> print(quadruples)
[,1] [,2] [,3] [,4]
[1,] 2 3 2 1
[2,] 2 3 3 3
[3,] 2 2 1 3
[4,] 2 1 1 3
[5,] 2 2 3 1
>
> proc.time()
user system elapsed
0.225 0.042 0.253
aroma.light.Rcheck/tests/wpca.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> for (zzz in 0) {
+
+ # This example requires plot3d() in R.basic [http://www.braju.com/R/]
+ if (!require(pkgName <- "R.basic", character.only=TRUE)) break
+
+ # -------------------------------------------------------------
+ # A first example
+ # -------------------------------------------------------------
+ # Simulate data from the model y <- a + bx + eps(bx)
+ x <- rexp(1000)
+ a <- c(2,15,3)
+ b <- c(2,3,15)
+ bx <- outer(b,x)
+ eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(x), mean=0, sd=0.1*x))
+ y <- a + bx + eps
+ y <- t(y)
+
+ # Add some outliers by permuting the dimensions for 1/3 of the observations
+ idx <- sample(1:nrow(y), size=1/3*nrow(y))
+ y[idx,] <- y[idx,c(2,3,1)]
+
+ # Down-weight the outliers W times to demonstrate how weights are used
+ W <- 10
+
+ # Plot the data with fitted lines at four different view points
+ N <- 4
+ theta <- seq(0,180,length.out=N)
+ phi <- rep(30, length.out=N)
+
+ # Use a different color for each set of weights
+ col <- topo.colors(W)
+
+ opar <- par(mar=c(1,1,1,1)+0.1)
+ layout(matrix(1:N, nrow=2, byrow=TRUE))
+ for (kk in seq(theta)) {
+ # Plot the data
+ plot3d(y, theta=theta[kk], phi=phi[kk])
+
+ # First, same weights for all observations
+ w <- rep(1, length=nrow(y))
+
+ for (ww in 1:W) {
+ # Fit a line using IWPCA through data
+ fit <- wpca(y, w=w, swapDirections=TRUE)
+
+ # Get the first principal component
+ ymid <- fit$xMean
+ d0 <- apply(y, MARGIN=2, FUN=min) - ymid
+ d1 <- apply(y, MARGIN=2, FUN=max) - ymid
+ b <- fit$vt[1,]
+ y0 <- -b * max(abs(d0))
+ y1 <- b * max(abs(d1))
+ yline <- matrix(c(y0,y1), nrow=length(b), ncol=2)
+ yline <- yline + ymid
+
+ points3d(t(ymid), col=col)
+ lines3d(t(yline), col=col)
+
+ # Down-weight outliers only, because here we know which they are.
+ w[idx] <- w[idx]/2
+ }
+
+ # Highlight the last one
+ lines3d(t(yline), col="red", lwd=3)
+ }
+
+ par(opar)
+
+ } # for (zzz in 0)
Loading required package: R.basic
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, :
there is no package called 'R.basic'
> rm(zzz)
>
> proc.time()
user system elapsed
0.300 0.041 0.324
aroma.light.Rcheck/tests/wpca2.matrix.Rout
R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 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.30.0 (2023-10-15) successfully loaded. See ?aroma.light for help.
>
> # -------------------------------------------------------------
> # A second example
> # -------------------------------------------------------------
> # Data
> x <- c(1,2,3,4,5)
> y <- c(2,4,3,3,6)
>
> opar <- par(bty="L")
> opalette <- palette(c("blue", "red", "black"))
> xlim <- ylim <- c(0,6)
>
> # Plot the data and the center mass
> plot(x,y, pch=16, cex=1.5, xlim=xlim, ylim=ylim)
> points(mean(x), mean(y), cex=2, lwd=2, col="blue")
>
>
> # Linear regression y ~ x
> fit <- lm(y ~ x)
> abline(fit, lty=1, col=1)
>
> # Linear regression y ~ x through without intercept
> fit <- lm(y ~ x - 1)
> abline(fit, lty=2, col=1)
>
>
> # Linear regression x ~ y
> fit <- lm(x ~ y)
> c <- coefficients(fit)
> b <- 1/c[2]
> a <- -b*c[1]
> abline(a=a, b=b, lty=1, col=2)
>
> # Linear regression x ~ y through without intercept
> fit <- lm(x ~ y - 1)
> b <- 1/coefficients(fit)
> abline(a=0, b=b, lty=2, col=2)
>
>
> # Orthogonal linear "regression"
> fit <- wpca(cbind(x,y))
>
> b <- fit$vt[1,2]/fit$vt[1,1]
> a <- fit$xMean[2]-b*fit$xMean[1]
> abline(a=a, b=b, lwd=2, col=3)
>
> # Orthogonal linear "regression" without intercept
> fit <- wpca(cbind(x,y), center=FALSE)
> b <- fit$vt[1,2]/fit$vt[1,1]
> a <- fit$xMean[2]-b*fit$xMean[1]
> abline(a=a, b=b, lty=2, lwd=2, col=3)
>
> legend(xlim[1],ylim[2], legend=c("lm(y~x)", "lm(y~x-1)", "lm(x~y)",
+ "lm(x~y-1)", "pca", "pca w/o intercept"), lty=rep(1:2,3),
+ lwd=rep(c(1,1,2),each=2), col=rep(1:3,each=2))
>
> palette(opalette)
> par(opar)
>
> proc.time()
user system elapsed
0.274 0.049 0.312
aroma.light.Rcheck/aroma.light-Ex.timings
| name | user | system | elapsed | |
| backtransformAffine | 0.003 | 0.000 | 0.003 | |
| backtransformPrincipalCurve | 0.391 | 0.016 | 0.408 | |
| calibrateMultiscan | 0.000 | 0.001 | 0.000 | |
| callNaiveGenotypes | 0.240 | 0.007 | 0.248 | |
| distanceBetweenLines | 0.071 | 0.008 | 0.079 | |
| findPeaksAndValleys | 0.025 | 0.007 | 0.032 | |
| fitPrincipalCurve | 0.521 | 0.008 | 0.528 | |
| fitXYCurve | 0.186 | 0.000 | 0.186 | |
| iwpca | 0.048 | 0.004 | 0.051 | |
| likelihood.smooth.spline | 0.108 | 0.000 | 0.107 | |
| medianPolish | 0.003 | 0.003 | 0.005 | |
| normalizeAffine | 5.653 | 0.029 | 5.684 | |
| normalizeCurveFit | 5.890 | 0.016 | 5.907 | |
| normalizeDifferencesToAverage | 0.266 | 0.008 | 0.274 | |
| normalizeFragmentLength | 1.493 | 0.028 | 1.521 | |
| normalizeQuantileRank | 0.703 | 0.004 | 0.707 | |
| normalizeQuantileRank.matrix | 0.042 | 0.000 | 0.042 | |
| normalizeQuantileSpline | 0.623 | 0.004 | 0.627 | |
| normalizeTumorBoost | 0.325 | 0.004 | 0.329 | |
| robustSmoothSpline | 0.324 | 0.010 | 0.334 | |
| sampleCorrelations | 0.194 | 0.001 | 0.195 | |
| sampleTuples | 0.000 | 0.001 | 0.001 | |
| wpca | 0.047 | 0.005 | 0.052 | |