##############################################################################
##############################################################################
###
### Running command:
###
### C:\Users\biocbuild\bbs-3.11-bioc\R\bin\R.exe CMD check --force-multiarch --install=check:sigaR.install-out.txt --library=C:\Users\biocbuild\bbs-3.11-bioc\R\library --no-vignettes --timings sigaR_1.36.0.tar.gz
###
##############################################################################
##############################################################################
* using log directory 'C:/Users/biocbuild/bbs-3.11-bioc/meat/sigaR.Rcheck'
* using R version 4.0.3 (2020-10-10)
* using platform: x86_64-w64-mingw32 (64-bit)
* using session charset: ISO8859-1
* using option '--no-vignettes'
* checking for file 'sigaR/DESCRIPTION' ... OK
* checking extension type ... Package
* this is package 'sigaR' version '1.36.0'
* 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 ... OK
* checking for portable file names ... OK
* checking whether package 'sigaR' can be installed ... OK
* checking installed package size ... OK
* checking package directory ... OK
* checking 'build' 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
* loading checks for arch 'i386'
** 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
* loading checks for arch 'x64'
** 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 dependencies in R code ... NOTE
There are ::: calls to the package's namespace in its code. A package
almost never needs to use ::: for its own objects:
'.alphabivariate' '.alphaest' '.pretest'
* 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 ... OK
* 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 contents of 'data' directory ... OK
* checking data for non-ASCII characters ... OK
* checking data for ASCII and uncompressed saves ... OK
* checking files in 'vignettes' ... OK
* checking examples ...
** running examples for arch 'i386' ... ERROR
Running examples in 'sigaR-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: RCMestimation
> ### Title: Fitting of the random coefficients model.
> ### Aliases: RCMestimation
>
> ### ** Examples
>
> # load data
> data(pollackCN16)
> data(pollackGE16)
>
> # select features belonging to a region
> ids <- getSegFeatures(20, pollackCN16)
perform input checks...
>
> # extract segmented log2 ratios of the region
> X <- t(segmented(pollackCN16)[ids[1], , drop=FALSE])
>
> # extract segmented log2 ratios of the region
> Y <- exprs(pollackGE16)[ids,]
>
> # center the expression data (row-wise)
> Y <- t(Y - apply(Y, 1, mean))
>
> # specify the linear constraint matrix
> R <- matrix(1, nrow=1)
>
> # fit the random coefficients model to the random data
> RCMresults <- RCMestimation(Y, X, R)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
sigaR
--- call from context ---
RCMestimation(Y, X, R)
--- call from argument ---
if (as.character(class(Y)) != "matrix") {
stop("Input (Y) is of wrong class.")
}
--- R stacktrace ---
where 1: RCMestimation(Y, X, R)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (Y, X, R, hypothesis = "H2", shrinkType = "none", estType = "normal",
corType = "unif", maxNoIt = 100, minSuccDist = 0.005, verbose = FALSE)
{
RCMmlH2 <- function(Y, X, R, maxNoIt, minSuccDist, shrinkType = "none",
estType = "normal", corType = "unif") {
QPest.RCM <- function(cov.mat, ed.Om2, X, Y.as.vec, R,
R.dual, weights, np) {
mat.calc.1 <- function(ev.cov, ed.Om2, X, weights) {
QPmat1 <- matrix(0, ncol = dim(X)[2], nrow = dim(X)[2])
for (i in 1:length(ed.Om2$values)) {
QPmat1 <- QPmat1 + (1/(ev.cov[1] + ed.Om2$values[i])) *
(t((matrix(1, nrow = length(ev.cov[-1]),
ncol = 1) %x% X)) %*% (matrix(ev.cov[-1] %x%
ed.Om2$vectors[, i], ncol = 1) * matrix(weights,
ncol = 1)) %*% (t((matrix(ev.cov[-1] %x%
ed.Om2$vectors[, i], ncol = 1) * matrix(weights,
ncol = 1))) %*% ((matrix(1, nrow = length(ev.cov[-1]),
ncol = 1) %x% X))))
}
return(QPmat1)
}
mat.calc.2 <- function(ev.cov, ed.Om2, X, Y.as.vec,
weights) {
QPmat2 <- matrix(0, ncol = 1, nrow = dim(X)[2])
for (i in 1:length(ed.Om2$values)) {
QPmat2 <- QPmat2 + (1/(ev.cov[1] + ed.Om2$values[i])) *
(t((matrix(1, nrow = length(ev.cov[-1]),
ncol = 1) %x% X)) %*% (matrix(ev.cov[-1] %x%
ed.Om2$vectors[, i], ncol = 1) * matrix(weights,
ncol = 1))) %*% (t(matrix(ev.cov[-1] %x%
ed.Om2$vectors[, i], ncol = 1) * matrix(weights,
ncol = 1)) %*% Y.as.vec)
}
return(QPmat2)
}
ed.cov <- eigen(cov.mat, symmetric = TRUE)
ed.cov <- rbind(ed.cov$values, ed.cov$vectors)
if (np == 1) {
slh.mat1 <- matrix(mean(apply(ed.cov, 2, mat.calc.1,
ed.Om2, X, weights)), ncol = dim(X)[2])
slh.mat2 <- matrix(mean(apply(ed.cov, 2, mat.calc.2,
ed.Om2, X, Y.as.vec, weights)), ncol = 1)
}
if (np > 1) {
slh.mat1 <- matrix(apply(t(apply(ed.cov, 2, mat.calc.1,
ed.Om2, X, weights)), 2, mean), ncol = dim(X)[2])
slh.mat2 <- matrix(apply(t(apply(ed.cov, 2, mat.calc.2,
ed.Om2, X, Y.as.vec, weights)), 2, mean), ncol = 1)
}
Dmat <- R %*% solve(slh.mat1) %*% t(R)
dvec <- R %*% solve(slh.mat1) %*% slh.mat2
qp.sol <- solve.QP(Dmat, -dvec, R.dual, bvec)
lambda <- matrix(qp.sol$solution, ncol = 1)
betas.conc <- solve(slh.mat1) %*% slh.mat2 + solve(slh.mat1) %*%
t(R) %*% lambda
return(as.numeric(betas.conc))
}
cov.mat.par.ests <- cov.par.estimation(projectY2Xcomp(Y,
X), dim(Y)[1], dim(Y)[2], estType, corType)
betas <- rep(0, dim(X)[2])
sigma2s <- cov.mat.par.ests$sigma2s
shrink.par <- shrink.par.calculation(sigma2s, dim(Y)[2],
shrinkType)
sigma2s <- (1 - shrink.par) * sigma2s + shrink.par *
mean(sigma2s)
tau2s <- rep(0, dim(X)[2])
rho <- cov.mat.par.ests$rho
if (rho < 0) {
warning("The estimate of rho is negative. To ensure the positive definiteness of the covariance matrix, it is set to zero in the remainder.")
rho <- 0
}
X.circ <- matrix(rep(1, dim(Y)[1]), ncol = 1) %x% X
Y.as.vec <- matrix(as.numeric(t(Y)), ncol = 1)
R.dual <- diag(rep(1, dim(R)[1]))
bvec <- rep(0, dim(R)[1])
for (i1 in 1:maxNoIt) {
if (verbose) {
cat(paste("ML estimation: iteration ", i1, " of ",
maxNoIt, sep = ""), "\n")
}
cov.mat <- .covMatConstruction(sigma2s, rho, corType)
tau.mat <- matrix(0, nrow = length(as.numeric(tau2s)),
ncol = length(as.numeric(tau2s)))
diag(tau.mat) <- as.numeric(tau2s)
Omega.mat.2 <- X %*% tau.mat %*% t(X)
ed.Om2 <- eigen(Omega.mat.2, symmetric = TRUE)
betas.old <- betas
weights <- matrix(1, ncol = 1, nrow = length(Y.as.vec))
betas <- QPest.RCM(cov.mat, ed.Om2, X, Y.as.vec,
R, R.dual, weights, dim(X)[2])
Y.res <- matrix(Y.as.vec - X.circ %*% matrix(betas,
ncol = 1), nrow = dim(Y)[1], byrow = TRUE)
Y.res <- Y.res - matrix(apply(Y.res, 1, mean), ncol = dim(Y.res)[2],
nrow = dim(Y.res)[1], byrow = FALSE)
tau.estimator <- function(gen.res, Xmat) {
solve(t(Xmat) %*% Xmat) %*% t(Xmat) %*% matrix(gen.res,
ncol = 1) %*% t(matrix(gen.res, ncol = 1)) %*%
Xmat %*% solve(t(Xmat) %*% Xmat)
}
tau2s.old <- tau2s
if (dim(X)[2] == 1) {
tau2s <- apply(Y.res, 1, tau.estimator, Xmat = X)
tau2s <- tau2s - sigma2s * solve(t(X) %*% X)
tau2s <- mean(tau2s)
if (tau2s < 0) {
tau2s <- 0
}
}
else {
tau2s <- apply(Y.res, 1, tau.estimator, Xmat = X)
tau2s.per.gene <- function(st, Xmat, np) {
return(diag(matrix(st[-1], ncol = np) - st[1] *
solve(t(Xmat) %*% Xmat)))
}
tau2s <- apply(rbind(sigma2s, tau2s), 2, tau2s.per.gene,
Xmat = X, np = dim(X)[2])
tau2s <- apply(tau2s, 1, mean)
tau2s[(tau2s < 0)] <- 0
}
if (estType == "robust") {
weights <- matrix(as.numeric(t(obs.weights(Y.res,
dim(Y)[1], dim(Y)[2], shrink.par))), ncol = 1)
rm(Y.res)
cov.mat <- .covMatConstruction(sigma2s, rho,
corType)
tau.mat <- matrix(0, nrow = length(as.numeric(tau2s)),
ncol = length(as.numeric(tau2s)))
diag(tau.mat) <- as.numeric(tau2s)
Omega.mat.2 <- X %*% tau.mat %*% t(X)
ed.Om2 <- eigen(Omega.mat.2, symmetric = TRUE)
betas <- QPest.RCM(cov.mat, ed.Om2, X, Y.as.vec,
R, R.dual, weights, dim(X)[2])
Y.res <- matrix(Y.as.vec - X.circ %*% matrix(betas,
ncol = 1), nrow = dim(Y)[1], byrow = TRUE)
Y.res <- Y.res - matrix(apply(Y.res, 1, mean),
ncol = dim(Y.res)[2], nrow = dim(Y.res)[1],
byrow = FALSE)
tau.estimatorR <- function(gen.res, Xmat, ns) {
solve(t(Xmat) %*% Xmat) %*% t(Xmat) %*% matrix(sqrt(gen.res[(ns +
1):(2 * ns)]) * gen.res[1:ns], ncol = 1) %*%
t(matrix(sqrt(gen.res[(ns + 1):(2 * ns)]) *
gen.res[1:ns], ncol = 1)) %*% Xmat %*%
solve(t(Xmat) %*% Xmat)
}
weights <- matrix(weights, nrow = dim(Y)[1],
byrow = TRUE)
betasPerGene <- apply(Y, 1, lm.uc, Xmat = X,
estType = "robust", mad.times = 4)
betasResiduals <- betasPerGene - betas
if (dim(X)[2] == 1) {
weights <- obs.weights.betas(betasResiduals,
mad.times = 5)
tau2s <- sum(weights * (betasResiduals)^2)/length(sigma2s) -
sum(sigma2s * solve(t(X) %*% X))/length(sigma2s)
tau2s[(tau2s < 0)] <- 0
}
else {
tau2s <- apply(cbind(Y.res, weights), 1, tau.estimatorR,
Xmat = X, dim(Y)[2])
tau2s.per.gene <- function(st, Xmat, np) {
return(diag(matrix(st[-1], ncol = np) - st[1] *
solve(t(Xmat) %*% Xmat)))
}
tau2s <- apply(rbind(sigma2s, tau2s), 2, tau2s.per.gene,
Xmat = X, np = dim(X)[2])
tau2s <- apply(tau2s, 1, mean)
tau2s[(tau2s < 0)] <- 0
}
}
if (max(abs(c(as.numeric(betas), sqrt(tau2s)) - c(as.numeric(betas.old),
sqrt(tau2s.old)))) < minSuccDist) {
break
}
}
estRes <- new("rcmFit", betas = as.numeric(round(betas,
digits = 5)), tau2s = round(tau2s, digits = 10),
sigma2s = as.numeric(round(sigma2s, digits = 5)),
rho = round(rho, digits = 5), shrinkage = round(shrink.par,
digits = 5), av.sigma2s = round(mean(sigma2s),
digits = 5), loglik = 1, corType = corType, X = X)
estRes@loglik <- .RCMloss(estRes, Y)
return(estRes)
}
RCMmlH0orH1 <- function(Y, X, R, shrinkType = "none", estType,
corType, hypothesis) {
cov.mat.par.ests <- cov.par.estimation(projectY2Xcomp(Y,
X), dim(Y)[1], dim(Y)[2], estType, corType)
betas <- rep(0, dim(X)[2])
sigma2s <- cov.mat.par.ests$sigma2s
shrink.par <- shrink.par.calculation(sigma2s, dim(Y)[2],
shrinkType)
sigma2s <- (1 - shrink.par) * sigma2s + shrink.par *
mean(sigma2s)
rho <- cov.mat.par.ests$rho
X.circ <- matrix(rep(1, dim(Y)[1]), ncol = 1) %x% X
Y.as.vec <- matrix(as.numeric(t(Y)), ncol = 1)
R.star.dual <- diag(rep(1, dim(R)[1]))
bvec <- rep(0, dim(R)[1])
if (hypothesis == "H0") {
cov.mat.inv <- solve(.covMatConstruction(sigma2s,
rho, corType))
slh.mat1 <- (matrix(1, ncol = dim(Y)[1], nrow = 1) %*%
cov.mat.inv %*% matrix(1, nrow = dim(Y)[1], ncol = 1)) %x%
(t(X) %*% X)
slh.mat2 <- ((matrix(1, ncol = dim(Y)[1], nrow = 1) %*%
cov.mat.inv) %x% t(X)) %*% Y.as.vec
beta.unconstr <- solve(slh.mat1) %*% slh.mat2
beta.eqconstr <- beta.unconstr - solve(slh.mat1) %*%
(t(R) %*% solve(R %*% solve(slh.mat1) %*% t(R)) %*%
R %*% beta.unconstr)
betas <- matrix(beta.eqconstr, ncol = dim(X)[2],
byrow = TRUE)
}
if (hypothesis == "H1") {
cov.mat.inv <- solve(.covMatConstruction(sigma2s,
rho, corType))
slh.mat <- solve((matrix(1, ncol = dim(Y)[1], nrow = 1) %*%
cov.mat.inv %*% matrix(1, nrow = dim(Y)[1], ncol = 1)) %x%
(t(X) %*% X))
Dmat <- R %*% slh.mat %*% t(R)
dvec <- -R %*% slh.mat %*% ((matrix(1, ncol = dim(Y)[1],
nrow = 1) %*% cov.mat.inv) %x% t(X)) %*% Y.as.vec
qp.sol <- solve.QP(Dmat, dvec, R.star.dual, bvec)
lambda <- matrix(qp.sol$solution, ncol = 1)
betas <- slh.mat %*% ((matrix(1, ncol = dim(Y)[1],
nrow = 1) %*% cov.mat.inv) %x% t(X)) %*% Y.as.vec +
slh.mat %*% t(R) %*% lambda
betas <- matrix(betas, ncol = dim(X)[2], byrow = TRUE)
}
if (estType == "robust") {
res.mat.star <- matrix(Y.as.vec - X.circ %*% as.numeric(t(betas)),
nrow = dim(Y)[1], byrow = TRUE)
weights <- matrix(as.numeric(t(obs.weights(res.mat.star,
dim(Y)[1], dim(Y)[2], shrink.par))), ncol = 1)
X.circ <- matrix(sqrt(weights), ncol = dim(X)[2],
nrow = dim(Y)[2] * dim(Y)[1]) * X.circ
ed.cov <- eigen(cov.mat.inv, symmetric = TRUE)
ed.cov <- rbind(ed.cov$values, ed.cov$vectors)
WY <- sqrt(weights) * Y.as.vec
if (dim(X)[2] == 1) {
slh.mat1 <- matrix(sum(apply(ed.cov, 2, slh.mat.weighted.one,
WX.circ = X.circ, ns = dim(Y)[2])), ncol = dim(X)[2])
if (hypothesis == "H1") {
slh.mat1 <- solve(slh.mat1)
}
slh.mat2 <- matrix(sum(apply(ed.cov, 2, slh.mat.weighted.two,
WX.circ = X.circ, WY = WY, ns = dim(Y)[2])),
nrow = dim(X)[2])
}
else {
slh.mat1 <- matrix(rowSums(apply(ed.cov, 2, slh.mat.weighted.one,
WX.circ = X.circ, ns = dim(Y)[2])), ncol = dim(X)[2])
if (hypothesis == "H1") {
slh.mat1 <- solve(slh.mat1)
}
slh.mat2 <- matrix(rowSums(apply(ed.cov, 2, slh.mat.weighted.two,
WX.circ = X.circ, WY = WY, ns = dim(Y)[2])),
nrow = dim(X)[2])
}
rm(ed.cov)
if (hypothesis == "H0") {
beta.unconstr <- solve(slh.mat1) %*% slh.mat2
beta.eqconstr <- beta.unconstr - solve(slh.mat1) %*%
(t(R) %*% solve(R %*% solve(slh.mat1) %*% t(R)) %*%
R %*% beta.unconstr)
betas <- matrix(beta.eqconstr, ncol = dim(X)[2],
byrow = TRUE)
}
if (hypothesis == "H1") {
Dmat <- R %*% slh.mat1 %*% t(R)
dvec <- -R %*% slh.mat1 %*% slh.mat2
qp.sol <- solve.QP(Dmat, dvec, R.star.dual, bvec)
lambda <- matrix(qp.sol$solution, ncol = 1)
betas <- slh.mat1 %*% slh.mat2 + slh.mat1 %*%
t(R) %*% lambda
betas <- matrix(betas, ncol = dim(X)[2], byrow = TRUE)
}
}
estRes <- new("rcmFit", betas = as.numeric(round(betas,
digits = 5)), tau2s = rep(0, length(betas)), sigma2s = as.numeric(round(sigma2s,
digits = 5)), rho = round(rho, digits = 5), shrinkage = round(shrink.par,
digits = 5), av.sigma2s = round(mean(sigma2s), digits = 5),
loglik = 1, corType = corType, X = X)
estRes@loglik <- .RCMloss(estRes, Y)
return(estRes)
}
projectY2Xcomp <- function(Y, X) {
Uorth <- Null(svd(X)$u)
Porth <- Uorth %*% solve(t(Uorth) %*% Uorth) %*% t(Uorth)
return(t(Porth %*% t(Y)))
}
slh.mat.weighted.one <- function(ev.cov, WX.circ, ns) {
junk <- 0
for (k in 1:ns) {
ev.unit <- matrix(0, ncol = 1, nrow = ns)
ev.unit[k, 1] <- 1
junk <- junk + ev.cov[1] * t(WX.circ) %*% (ev.cov[-1] %x%
ev.unit) %*% t(t(WX.circ) %*% (ev.cov[-1] %x%
ev.unit))
}
return(junk)
}
slh.mat.weighted.two <- function(ev.cov, WX.circ, WY, ns) {
junk <- 0
for (k in 1:ns) {
ev.unit <- matrix(0, ncol = 1, nrow = ns)
ev.unit[k, 1] <- 1
junk <- junk + ev.cov[1] * t(WX.circ) %*% (ev.cov[-1] %x%
ev.unit) %*% (t(ev.cov[-1] %x% ev.unit) %*% WY)
}
return(junk)
}
shrink.par.calculation <- function(sigma2s, ns, shrinkType) {
if (shrinkType == "none") {
shrink.par <- 0
}
if (shrinkType == "full") {
shrink.par <- 1
}
if (shrinkType == "opt") {
teller <- (2 * length(sigma2s) - 3) * sum(sigma2s^2)/((ns -
1) * length(sigma2s))
noemer <- sum((sigma2s - mean(sigma2s))^2) + (2 *
length(sigma2s) - 4) * sum(sigma2s^2)/((ns -
1) * length(sigma2s))
shrink.par <- max(0, min(1, teller/noemer))
}
return(shrink.par)
}
cov.par.estimation <- function(res.mat, ng, ns, estType = "robust",
corType = "unif") {
if (corType == "unif") {
if (estType == "normal") {
res.mat <- res.mat - matrix(apply(res.mat, 1,
mean), ncol = dim(res.mat)[2], nrow = dim(res.mat)[1],
byrow = FALSE)
res.inprods.mat <- res.mat %*% t(res.mat)
rho <- (sum(res.inprods.mat[upper.tri(res.inprods.mat)])/(ng *
(ng - 1)/2)/(sum(diag(res.inprods.mat))/ng))
sigma2s <- diag(res.inprods.mat)/(ns)
}
if (estType == "robust") {
sigma2s <- apply(res.mat, 1, mad)^2
cov.sum <- function(j, ng, X) {
sum(apply(X[(j + 1):ng, , drop = FALSE], 1,
function(x, y) {
(mad(x + y))^2 - (mad(x - y))^2
}, X[j, ]))
}
cov <- sum(sapply(c(1:(ng - 1)), cov.sum, ng = ng,
res.mat))
rho <- (2 * cov/(4 * ng * (ng - 1)))/mean(sigma2s)
}
}
if (corType == "ar1") {
if (estType == "normal") {
res.mat <- res.mat - matrix(apply(res.mat, 1,
mean), ncol = dim(res.mat)[2], nrow = dim(res.mat)[1],
byrow = FALSE)
res.inprods.mat <- res.mat %*% t(res.mat)
sigma2s <- diag(res.inprods.mat)/(ns)
noemer <- (sum(diag(res.inprods.mat))/ng)
res.inprods.mat <- res.inprods.mat[-ng, -1]
rho <- (sum(diag(res.inprods.mat))/(ng - 1)/noemer)
}
if (estType == "robust") {
sigma2s <- apply(res.mat, 1, mad)^2
cov.sum <- function(j, ng, X) {
sum(apply(X[j + 1, , drop = FALSE], 1, function(x,
y) {
(mad(x + y))^2 - (mad(x - y))^2
}, X[j, ]))
}
cov <- sum(sapply(c(1:(ng - 1)), cov.sum, ng = ng,
res.mat))
rho <- (cov/(4 * (ng - 1)))/mean(sigma2s)
}
}
cov.pars <- list()
cov.pars$rho <- rho
cov.pars$sigma2s <- sigma2s
return(cov.pars)
}
obs.weights.betas <- function(res.mat, mad.times = 5) {
res.mat <- res.mat - median(res.mat)
mads <- mad(res.mat)
res.mat <- res.mat * 1/(mad.times * mads)
res.mat[abs(res.mat) > 1] <- 1
weights <- (1 - res.mat^2)^2
return(weights)
}
obs.weights <- function(res.mat, ng, ns, shrink.par = 0,
mad.times = 5) {
res.mat <- res.mat - apply(res.mat, 1, median)
mads <- apply(res.mat, 1, mad)
mads <- (1 - shrink.par) * mads + shrink.par * mean(mads)
res.mat <- t(apply(cbind(res.mat, mads), 1, function(x,
mad.times) {
spread <- x[length(x)]
x <- x[-length(x)]
x[abs(x) > mad.times * spread] <- mad.times * spread
return(x)
}, mad.times))
res.mat <- res.mat * matrix(1/(mad.times * mads), ncol = ns,
nrow = ng, byrow = FALSE)
weights <- (1 - res.mat^2)^2
return(weights)
}
lm.uc <- function(Y, Xmat, estType, mad.times = 3) {
sigma <- 1
cov.mat.inv <- 1/sigma * diag(rep(1, length(Y)))
betas <- solve(t(Xmat) %*% cov.mat.inv %*% Xmat) %*%
t(Xmat) %*% cov.mat.inv %*% Y
if (estType == "robust") {
residuals <- Y - X %*% betas
res.mad <- mad(residuals)
residuals[residuals > mad.times * res.mad] <- mad.times *
res.mad
weights <- (1 - (residuals/(mad.times * mad(residuals)))^2)^2
betas <- solve(t(Xmat) %*% diag(sqrt(as.numeric(weights))) %*%
cov.mat.inv %*% diag(sqrt(as.numeric(weights))) %*%
Xmat) %*% t(Xmat) %*% cov.mat.inv %*% Y
}
return(betas)
}
if (as.character(class(Y)) != "matrix") {
stop("Input (Y) is of wrong class.")
}
if (sum(is.na(Y)) != 0) {
stop("Y contains missings.")
}
if (as.character(class(X)) != "matrix") {
stop("Input (X) is of wrong class.")
}
if (sum(is.na(X)) != 0) {
stop("X contains missings.")
}
if (as.character(class(R)) != "matrix") {
stop("Input (R) is of wrong class.")
}
if (sum(is.na(R)) != 0) {
stop("R contains missings.")
}
if (!(as.character(class(maxNoIt)) == "numeric" | as.character(class(maxNoIt)) ==
"integer")) {
stop("Input (maxNoIt) is of wrong class.")
}
if (as.character(class(minSuccDist)) != "numeric") {
stop("Input (minSuccDist) is of wrong class.")
}
if (as.character(class(corType)) != "character") {
stop("Input (corType) is of wrong class.")
}
if (as.character(class(estType)) != "character") {
stop("Input (estType) is of wrong class.")
}
if (as.character(class(shrinkType)) != "character") {
stop("Input (shrinkType) is of wrong class.")
}
if (as.character(class(hypothesis)) != "character") {
stop("Input (hypothesis) is of wrong class.")
}
if (dim(Y)[1] < 2) {
stop("This is a multivariate procedure, provide a multivariate data set.")
}
if (!(corType %in% c("unif", "ar1"))) {
stop("corType parameter ill-specified.")
}
if (!(estType %in% c("normal", "robust"))) {
stop("estType parameter ill-specified.")
}
if (!(shrinkType %in% c("none", "opt", "full"))) {
stop("shrinkType parameter ill-specified.")
}
if (!(hypothesis %in% c("H0", "H1", "H2"))) {
stop("hypothesis parameter ill-specified.")
}
if (dim(Y)[1] != dim(X)[1]) {
stop("Dimension mismatch between expression and design matrix.")
}
if (dim(X)[2] != dim(R)[2]) {
stop("Dimension mismatch between constraint and design matrix.")
}
if (minSuccDist <= 0) {
stop("Stopping criterion incorrect.")
}
if (maxNoIt < 1) {
stop("Maximum number of iterations smaller than 1.")
}
Y <- t(Y)
if (hypothesis == "H0") {
return(RCMmlH0orH1(Y, X, R, shrinkType = shrinkType,
estType = estType, corType = corType, hypothesis = hypothesis))
}
if (hypothesis == "H1") {
return(RCMmlH0orH1(Y, X, R, shrinkType = shrinkType,
estType = estType, corType = corType, hypothesis = hypothesis))
}
if (hypothesis == "H2") {
return(RCMmlH2(Y, X, R, maxNoIt, minSuccDist, shrinkType = shrinkType,
estType = estType, corType = corType))
}
}
<bytecode: 0x0c2f61b0>
<environment: namespace:sigaR>
--- function search by body ---
Function RCMestimation in namespace sigaR has this body.
----------- END OF FAILURE REPORT --------------
Fatal error: the condition has length > 1
** running examples for arch 'x64' ... ERROR
Running examples in 'sigaR-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: RCMestimation
> ### Title: Fitting of the random coefficients model.
> ### Aliases: RCMestimation
>
> ### ** Examples
>
> # load data
> data(pollackCN16)
> data(pollackGE16)
>
> # select features belonging to a region
> ids <- getSegFeatures(20, pollackCN16)
perform input checks...
>
> # extract segmented log2 ratios of the region
> X <- t(segmented(pollackCN16)[ids[1], , drop=FALSE])
>
> # extract segmented log2 ratios of the region
> Y <- exprs(pollackGE16)[ids,]
>
> # center the expression data (row-wise)
> Y <- t(Y - apply(Y, 1, mean))
>
> # specify the linear constraint matrix
> R <- matrix(1, nrow=1)
>
> # fit the random coefficients model to the random data
> RCMresults <- RCMestimation(Y, X, R)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
sigaR
--- call from context ---
RCMestimation(Y, X, R)
--- call from argument ---
if (as.character(class(Y)) != "matrix") {
stop("Input (Y) is of wrong class.")
}
--- R stacktrace ---
where 1: RCMestimation(Y, X, R)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (Y, X, R, hypothesis = "H2", shrinkType = "none", estType = "normal",
corType = "unif", maxNoIt = 100, minSuccDist = 0.005, verbose = FALSE)
{
RCMmlH2 <- function(Y, X, R, maxNoIt, minSuccDist, shrinkType = "none",
estType = "normal", corType = "unif") {
QPest.RCM <- function(cov.mat, ed.Om2, X, Y.as.vec, R,
R.dual, weights, np) {
mat.calc.1 <- function(ev.cov, ed.Om2, X, weights) {
QPmat1 <- matrix(0, ncol = dim(X)[2], nrow = dim(X)[2])
for (i in 1:length(ed.Om2$values)) {
QPmat1 <- QPmat1 + (1/(ev.cov[1] + ed.Om2$values[i])) *
(t((matrix(1, nrow = length(ev.cov[-1]),
ncol = 1) %x% X)) %*% (matrix(ev.cov[-1] %x%
ed.Om2$vectors[, i], ncol = 1) * matrix(weights,
ncol = 1)) %*% (t((matrix(ev.cov[-1] %x%
ed.Om2$vectors[, i], ncol = 1) * matrix(weights,
ncol = 1))) %*% ((matrix(1, nrow = length(ev.cov[-1]),
ncol = 1) %x% X))))
}
return(QPmat1)
}
mat.calc.2 <- function(ev.cov, ed.Om2, X, Y.as.vec,
weights) {
QPmat2 <- matrix(0, ncol = 1, nrow = dim(X)[2])
for (i in 1:length(ed.Om2$values)) {
QPmat2 <- QPmat2 + (1/(ev.cov[1] + ed.Om2$values[i])) *
(t((matrix(1, nrow = length(ev.cov[-1]),
ncol = 1) %x% X)) %*% (matrix(ev.cov[-1] %x%
ed.Om2$vectors[, i], ncol = 1) * matrix(weights,
ncol = 1))) %*% (t(matrix(ev.cov[-1] %x%
ed.Om2$vectors[, i], ncol = 1) * matrix(weights,
ncol = 1)) %*% Y.as.vec)
}
return(QPmat2)
}
ed.cov <- eigen(cov.mat, symmetric = TRUE)
ed.cov <- rbind(ed.cov$values, ed.cov$vectors)
if (np == 1) {
slh.mat1 <- matrix(mean(apply(ed.cov, 2, mat.calc.1,
ed.Om2, X, weights)), ncol = dim(X)[2])
slh.mat2 <- matrix(mean(apply(ed.cov, 2, mat.calc.2,
ed.Om2, X, Y.as.vec, weights)), ncol = 1)
}
if (np > 1) {
slh.mat1 <- matrix(apply(t(apply(ed.cov, 2, mat.calc.1,
ed.Om2, X, weights)), 2, mean), ncol = dim(X)[2])
slh.mat2 <- matrix(apply(t(apply(ed.cov, 2, mat.calc.2,
ed.Om2, X, Y.as.vec, weights)), 2, mean), ncol = 1)
}
Dmat <- R %*% solve(slh.mat1) %*% t(R)
dvec <- R %*% solve(slh.mat1) %*% slh.mat2
qp.sol <- solve.QP(Dmat, -dvec, R.dual, bvec)
lambda <- matrix(qp.sol$solution, ncol = 1)
betas.conc <- solve(slh.mat1) %*% slh.mat2 + solve(slh.mat1) %*%
t(R) %*% lambda
return(as.numeric(betas.conc))
}
cov.mat.par.ests <- cov.par.estimation(projectY2Xcomp(Y,
X), dim(Y)[1], dim(Y)[2], estType, corType)
betas <- rep(0, dim(X)[2])
sigma2s <- cov.mat.par.ests$sigma2s
shrink.par <- shrink.par.calculation(sigma2s, dim(Y)[2],
shrinkType)
sigma2s <- (1 - shrink.par) * sigma2s + shrink.par *
mean(sigma2s)
tau2s <- rep(0, dim(X)[2])
rho <- cov.mat.par.ests$rho
if (rho < 0) {
warning("The estimate of rho is negative. To ensure the positive definiteness of the covariance matrix, it is set to zero in the remainder.")
rho <- 0
}
X.circ <- matrix(rep(1, dim(Y)[1]), ncol = 1) %x% X
Y.as.vec <- matrix(as.numeric(t(Y)), ncol = 1)
R.dual <- diag(rep(1, dim(R)[1]))
bvec <- rep(0, dim(R)[1])
for (i1 in 1:maxNoIt) {
if (verbose) {
cat(paste("ML estimation: iteration ", i1, " of ",
maxNoIt, sep = ""), "\n")
}
cov.mat <- .covMatConstruction(sigma2s, rho, corType)
tau.mat <- matrix(0, nrow = length(as.numeric(tau2s)),
ncol = length(as.numeric(tau2s)))
diag(tau.mat) <- as.numeric(tau2s)
Omega.mat.2 <- X %*% tau.mat %*% t(X)
ed.Om2 <- eigen(Omega.mat.2, symmetric = TRUE)
betas.old <- betas
weights <- matrix(1, ncol = 1, nrow = length(Y.as.vec))
betas <- QPest.RCM(cov.mat, ed.Om2, X, Y.as.vec,
R, R.dual, weights, dim(X)[2])
Y.res <- matrix(Y.as.vec - X.circ %*% matrix(betas,
ncol = 1), nrow = dim(Y)[1], byrow = TRUE)
Y.res <- Y.res - matrix(apply(Y.res, 1, mean), ncol = dim(Y.res)[2],
nrow = dim(Y.res)[1], byrow = FALSE)
tau.estimator <- function(gen.res, Xmat) {
solve(t(Xmat) %*% Xmat) %*% t(Xmat) %*% matrix(gen.res,
ncol = 1) %*% t(matrix(gen.res, ncol = 1)) %*%
Xmat %*% solve(t(Xmat) %*% Xmat)
}
tau2s.old <- tau2s
if (dim(X)[2] == 1) {
tau2s <- apply(Y.res, 1, tau.estimator, Xmat = X)
tau2s <- tau2s - sigma2s * solve(t(X) %*% X)
tau2s <- mean(tau2s)
if (tau2s < 0) {
tau2s <- 0
}
}
else {
tau2s <- apply(Y.res, 1, tau.estimator, Xmat = X)
tau2s.per.gene <- function(st, Xmat, np) {
return(diag(matrix(st[-1], ncol = np) - st[1] *
solve(t(Xmat) %*% Xmat)))
}
tau2s <- apply(rbind(sigma2s, tau2s), 2, tau2s.per.gene,
Xmat = X, np = dim(X)[2])
tau2s <- apply(tau2s, 1, mean)
tau2s[(tau2s < 0)] <- 0
}
if (estType == "robust") {
weights <- matrix(as.numeric(t(obs.weights(Y.res,
dim(Y)[1], dim(Y)[2], shrink.par))), ncol = 1)
rm(Y.res)
cov.mat <- .covMatConstruction(sigma2s, rho,
corType)
tau.mat <- matrix(0, nrow = length(as.numeric(tau2s)),
ncol = length(as.numeric(tau2s)))
diag(tau.mat) <- as.numeric(tau2s)
Omega.mat.2 <- X %*% tau.mat %*% t(X)
ed.Om2 <- eigen(Omega.mat.2, symmetric = TRUE)
betas <- QPest.RCM(cov.mat, ed.Om2, X, Y.as.vec,
R, R.dual, weights, dim(X)[2])
Y.res <- matrix(Y.as.vec - X.circ %*% matrix(betas,
ncol = 1), nrow = dim(Y)[1], byrow = TRUE)
Y.res <- Y.res - matrix(apply(Y.res, 1, mean),
ncol = dim(Y.res)[2], nrow = dim(Y.res)[1],
byrow = FALSE)
tau.estimatorR <- function(gen.res, Xmat, ns) {
solve(t(Xmat) %*% Xmat) %*% t(Xmat) %*% matrix(sqrt(gen.res[(ns +
1):(2 * ns)]) * gen.res[1:ns], ncol = 1) %*%
t(matrix(sqrt(gen.res[(ns + 1):(2 * ns)]) *
gen.res[1:ns], ncol = 1)) %*% Xmat %*%
solve(t(Xmat) %*% Xmat)
}
weights <- matrix(weights, nrow = dim(Y)[1],
byrow = TRUE)
betasPerGene <- apply(Y, 1, lm.uc, Xmat = X,
estType = "robust", mad.times = 4)
betasResiduals <- betasPerGene - betas
if (dim(X)[2] == 1) {
weights <- obs.weights.betas(betasResiduals,
mad.times = 5)
tau2s <- sum(weights * (betasResiduals)^2)/length(sigma2s) -
sum(sigma2s * solve(t(X) %*% X))/length(sigma2s)
tau2s[(tau2s < 0)] <- 0
}
else {
tau2s <- apply(cbind(Y.res, weights), 1, tau.estimatorR,
Xmat = X, dim(Y)[2])
tau2s.per.gene <- function(st, Xmat, np) {
return(diag(matrix(st[-1], ncol = np) - st[1] *
solve(t(Xmat) %*% Xmat)))
}
tau2s <- apply(rbind(sigma2s, tau2s), 2, tau2s.per.gene,
Xmat = X, np = dim(X)[2])
tau2s <- apply(tau2s, 1, mean)
tau2s[(tau2s < 0)] <- 0
}
}
if (max(abs(c(as.numeric(betas), sqrt(tau2s)) - c(as.numeric(betas.old),
sqrt(tau2s.old)))) < minSuccDist) {
break
}
}
estRes <- new("rcmFit", betas = as.numeric(round(betas,
digits = 5)), tau2s = round(tau2s, digits = 10),
sigma2s = as.numeric(round(sigma2s, digits = 5)),
rho = round(rho, digits = 5), shrinkage = round(shrink.par,
digits = 5), av.sigma2s = round(mean(sigma2s),
digits = 5), loglik = 1, corType = corType, X = X)
estRes@loglik <- .RCMloss(estRes, Y)
return(estRes)
}
RCMmlH0orH1 <- function(Y, X, R, shrinkType = "none", estType,
corType, hypothesis) {
cov.mat.par.ests <- cov.par.estimation(projectY2Xcomp(Y,
X), dim(Y)[1], dim(Y)[2], estType, corType)
betas <- rep(0, dim(X)[2])
sigma2s <- cov.mat.par.ests$sigma2s
shrink.par <- shrink.par.calculation(sigma2s, dim(Y)[2],
shrinkType)
sigma2s <- (1 - shrink.par) * sigma2s + shrink.par *
mean(sigma2s)
rho <- cov.mat.par.ests$rho
X.circ <- matrix(rep(1, dim(Y)[1]), ncol = 1) %x% X
Y.as.vec <- matrix(as.numeric(t(Y)), ncol = 1)
R.star.dual <- diag(rep(1, dim(R)[1]))
bvec <- rep(0, dim(R)[1])
if (hypothesis == "H0") {
cov.mat.inv <- solve(.covMatConstruction(sigma2s,
rho, corType))
slh.mat1 <- (matrix(1, ncol = dim(Y)[1], nrow = 1) %*%
cov.mat.inv %*% matrix(1, nrow = dim(Y)[1], ncol = 1)) %x%
(t(X) %*% X)
slh.mat2 <- ((matrix(1, ncol = dim(Y)[1], nrow = 1) %*%
cov.mat.inv) %x% t(X)) %*% Y.as.vec
beta.unconstr <- solve(slh.mat1) %*% slh.mat2
beta.eqconstr <- beta.unconstr - solve(slh.mat1) %*%
(t(R) %*% solve(R %*% solve(slh.mat1) %*% t(R)) %*%
R %*% beta.unconstr)
betas <- matrix(beta.eqconstr, ncol = dim(X)[2],
byrow = TRUE)
}
if (hypothesis == "H1") {
cov.mat.inv <- solve(.covMatConstruction(sigma2s,
rho, corType))
slh.mat <- solve((matrix(1, ncol = dim(Y)[1], nrow = 1) %*%
cov.mat.inv %*% matrix(1, nrow = dim(Y)[1], ncol = 1)) %x%
(t(X) %*% X))
Dmat <- R %*% slh.mat %*% t(R)
dvec <- -R %*% slh.mat %*% ((matrix(1, ncol = dim(Y)[1],
nrow = 1) %*% cov.mat.inv) %x% t(X)) %*% Y.as.vec
qp.sol <- solve.QP(Dmat, dvec, R.star.dual, bvec)
lambda <- matrix(qp.sol$solution, ncol = 1)
betas <- slh.mat %*% ((matrix(1, ncol = dim(Y)[1],
nrow = 1) %*% cov.mat.inv) %x% t(X)) %*% Y.as.vec +
slh.mat %*% t(R) %*% lambda
betas <- matrix(betas, ncol = dim(X)[2], byrow = TRUE)
}
if (estType == "robust") {
res.mat.star <- matrix(Y.as.vec - X.circ %*% as.numeric(t(betas)),
nrow = dim(Y)[1], byrow = TRUE)
weights <- matrix(as.numeric(t(obs.weights(res.mat.star,
dim(Y)[1], dim(Y)[2], shrink.par))), ncol = 1)
X.circ <- matrix(sqrt(weights), ncol = dim(X)[2],
nrow = dim(Y)[2] * dim(Y)[1]) * X.circ
ed.cov <- eigen(cov.mat.inv, symmetric = TRUE)
ed.cov <- rbind(ed.cov$values, ed.cov$vectors)
WY <- sqrt(weights) * Y.as.vec
if (dim(X)[2] == 1) {
slh.mat1 <- matrix(sum(apply(ed.cov, 2, slh.mat.weighted.one,
WX.circ = X.circ, ns = dim(Y)[2])), ncol = dim(X)[2])
if (hypothesis == "H1") {
slh.mat1 <- solve(slh.mat1)
}
slh.mat2 <- matrix(sum(apply(ed.cov, 2, slh.mat.weighted.two,
WX.circ = X.circ, WY = WY, ns = dim(Y)[2])),
nrow = dim(X)[2])
}
else {
slh.mat1 <- matrix(rowSums(apply(ed.cov, 2, slh.mat.weighted.one,
WX.circ = X.circ, ns = dim(Y)[2])), ncol = dim(X)[2])
if (hypothesis == "H1") {
slh.mat1 <- solve(slh.mat1)
}
slh.mat2 <- matrix(rowSums(apply(ed.cov, 2, slh.mat.weighted.two,
WX.circ = X.circ, WY = WY, ns = dim(Y)[2])),
nrow = dim(X)[2])
}
rm(ed.cov)
if (hypothesis == "H0") {
beta.unconstr <- solve(slh.mat1) %*% slh.mat2
beta.eqconstr <- beta.unconstr - solve(slh.mat1) %*%
(t(R) %*% solve(R %*% solve(slh.mat1) %*% t(R)) %*%
R %*% beta.unconstr)
betas <- matrix(beta.eqconstr, ncol = dim(X)[2],
byrow = TRUE)
}
if (hypothesis == "H1") {
Dmat <- R %*% slh.mat1 %*% t(R)
dvec <- -R %*% slh.mat1 %*% slh.mat2
qp.sol <- solve.QP(Dmat, dvec, R.star.dual, bvec)
lambda <- matrix(qp.sol$solution, ncol = 1)
betas <- slh.mat1 %*% slh.mat2 + slh.mat1 %*%
t(R) %*% lambda
betas <- matrix(betas, ncol = dim(X)[2], byrow = TRUE)
}
}
estRes <- new("rcmFit", betas = as.numeric(round(betas,
digits = 5)), tau2s = rep(0, length(betas)), sigma2s = as.numeric(round(sigma2s,
digits = 5)), rho = round(rho, digits = 5), shrinkage = round(shrink.par,
digits = 5), av.sigma2s = round(mean(sigma2s), digits = 5),
loglik = 1, corType = corType, X = X)
estRes@loglik <- .RCMloss(estRes, Y)
return(estRes)
}
projectY2Xcomp <- function(Y, X) {
Uorth <- Null(svd(X)$u)
Porth <- Uorth %*% solve(t(Uorth) %*% Uorth) %*% t(Uorth)
return(t(Porth %*% t(Y)))
}
slh.mat.weighted.one <- function(ev.cov, WX.circ, ns) {
junk <- 0
for (k in 1:ns) {
ev.unit <- matrix(0, ncol = 1, nrow = ns)
ev.unit[k, 1] <- 1
junk <- junk + ev.cov[1] * t(WX.circ) %*% (ev.cov[-1] %x%
ev.unit) %*% t(t(WX.circ) %*% (ev.cov[-1] %x%
ev.unit))
}
return(junk)
}
slh.mat.weighted.two <- function(ev.cov, WX.circ, WY, ns) {
junk <- 0
for (k in 1:ns) {
ev.unit <- matrix(0, ncol = 1, nrow = ns)
ev.unit[k, 1] <- 1
junk <- junk + ev.cov[1] * t(WX.circ) %*% (ev.cov[-1] %x%
ev.unit) %*% (t(ev.cov[-1] %x% ev.unit) %*% WY)
}
return(junk)
}
shrink.par.calculation <- function(sigma2s, ns, shrinkType) {
if (shrinkType == "none") {
shrink.par <- 0
}
if (shrinkType == "full") {
shrink.par <- 1
}
if (shrinkType == "opt") {
teller <- (2 * length(sigma2s) - 3) * sum(sigma2s^2)/((ns -
1) * length(sigma2s))
noemer <- sum((sigma2s - mean(sigma2s))^2) + (2 *
length(sigma2s) - 4) * sum(sigma2s^2)/((ns -
1) * length(sigma2s))
shrink.par <- max(0, min(1, teller/noemer))
}
return(shrink.par)
}
cov.par.estimation <- function(res.mat, ng, ns, estType = "robust",
corType = "unif") {
if (corType == "unif") {
if (estType == "normal") {
res.mat <- res.mat - matrix(apply(res.mat, 1,
mean), ncol = dim(res.mat)[2], nrow = dim(res.mat)[1],
byrow = FALSE)
res.inprods.mat <- res.mat %*% t(res.mat)
rho <- (sum(res.inprods.mat[upper.tri(res.inprods.mat)])/(ng *
(ng - 1)/2)/(sum(diag(res.inprods.mat))/ng))
sigma2s <- diag(res.inprods.mat)/(ns)
}
if (estType == "robust") {
sigma2s <- apply(res.mat, 1, mad)^2
cov.sum <- function(j, ng, X) {
sum(apply(X[(j + 1):ng, , drop = FALSE], 1,
function(x, y) {
(mad(x + y))^2 - (mad(x - y))^2
}, X[j, ]))
}
cov <- sum(sapply(c(1:(ng - 1)), cov.sum, ng = ng,
res.mat))
rho <- (2 * cov/(4 * ng * (ng - 1)))/mean(sigma2s)
}
}
if (corType == "ar1") {
if (estType == "normal") {
res.mat <- res.mat - matrix(apply(res.mat, 1,
mean), ncol = dim(res.mat)[2], nrow = dim(res.mat)[1],
byrow = FALSE)
res.inprods.mat <- res.mat %*% t(res.mat)
sigma2s <- diag(res.inprods.mat)/(ns)
noemer <- (sum(diag(res.inprods.mat))/ng)
res.inprods.mat <- res.inprods.mat[-ng, -1]
rho <- (sum(diag(res.inprods.mat))/(ng - 1)/noemer)
}
if (estType == "robust") {
sigma2s <- apply(res.mat, 1, mad)^2
cov.sum <- function(j, ng, X) {
sum(apply(X[j + 1, , drop = FALSE], 1, function(x,
y) {
(mad(x + y))^2 - (mad(x - y))^2
}, X[j, ]))
}
cov <- sum(sapply(c(1:(ng - 1)), cov.sum, ng = ng,
res.mat))
rho <- (cov/(4 * (ng - 1)))/mean(sigma2s)
}
}
cov.pars <- list()
cov.pars$rho <- rho
cov.pars$sigma2s <- sigma2s
return(cov.pars)
}
obs.weights.betas <- function(res.mat, mad.times = 5) {
res.mat <- res.mat - median(res.mat)
mads <- mad(res.mat)
res.mat <- res.mat * 1/(mad.times * mads)
res.mat[abs(res.mat) > 1] <- 1
weights <- (1 - res.mat^2)^2
return(weights)
}
obs.weights <- function(res.mat, ng, ns, shrink.par = 0,
mad.times = 5) {
res.mat <- res.mat - apply(res.mat, 1, median)
mads <- apply(res.mat, 1, mad)
mads <- (1 - shrink.par) * mads + shrink.par * mean(mads)
res.mat <- t(apply(cbind(res.mat, mads), 1, function(x,
mad.times) {
spread <- x[length(x)]
x <- x[-length(x)]
x[abs(x) > mad.times * spread] <- mad.times * spread
return(x)
}, mad.times))
res.mat <- res.mat * matrix(1/(mad.times * mads), ncol = ns,
nrow = ng, byrow = FALSE)
weights <- (1 - res.mat^2)^2
return(weights)
}
lm.uc <- function(Y, Xmat, estType, mad.times = 3) {
sigma <- 1
cov.mat.inv <- 1/sigma * diag(rep(1, length(Y)))
betas <- solve(t(Xmat) %*% cov.mat.inv %*% Xmat) %*%
t(Xmat) %*% cov.mat.inv %*% Y
if (estType == "robust") {
residuals <- Y - X %*% betas
res.mad <- mad(residuals)
residuals[residuals > mad.times * res.mad] <- mad.times *
res.mad
weights <- (1 - (residuals/(mad.times * mad(residuals)))^2)^2
betas <- solve(t(Xmat) %*% diag(sqrt(as.numeric(weights))) %*%
cov.mat.inv %*% diag(sqrt(as.numeric(weights))) %*%
Xmat) %*% t(Xmat) %*% cov.mat.inv %*% Y
}
return(betas)
}
if (as.character(class(Y)) != "matrix") {
stop("Input (Y) is of wrong class.")
}
if (sum(is.na(Y)) != 0) {
stop("Y contains missings.")
}
if (as.character(class(X)) != "matrix") {
stop("Input (X) is of wrong class.")
}
if (sum(is.na(X)) != 0) {
stop("X contains missings.")
}
if (as.character(class(R)) != "matrix") {
stop("Input (R) is of wrong class.")
}
if (sum(is.na(R)) != 0) {
stop("R contains missings.")
}
if (!(as.character(class(maxNoIt)) == "numeric" | as.character(class(maxNoIt)) ==
"integer")) {
stop("Input (maxNoIt) is of wrong class.")
}
if (as.character(class(minSuccDist)) != "numeric") {
stop("Input (minSuccDist) is of wrong class.")
}
if (as.character(class(corType)) != "character") {
stop("Input (corType) is of wrong class.")
}
if (as.character(class(estType)) != "character") {
stop("Input (estType) is of wrong class.")
}
if (as.character(class(shrinkType)) != "character") {
stop("Input (shrinkType) is of wrong class.")
}
if (as.character(class(hypothesis)) != "character") {
stop("Input (hypothesis) is of wrong class.")
}
if (dim(Y)[1] < 2) {
stop("This is a multivariate procedure, provide a multivariate data set.")
}
if (!(corType %in% c("unif", "ar1"))) {
stop("corType parameter ill-specified.")
}
if (!(estType %in% c("normal", "robust"))) {
stop("estType parameter ill-specified.")
}
if (!(shrinkType %in% c("none", "opt", "full"))) {
stop("shrinkType parameter ill-specified.")
}
if (!(hypothesis %in% c("H0", "H1", "H2"))) {
stop("hypothesis parameter ill-specified.")
}
if (dim(Y)[1] != dim(X)[1]) {
stop("Dimension mismatch between expression and design matrix.")
}
if (dim(X)[2] != dim(R)[2]) {
stop("Dimension mismatch between constraint and design matrix.")
}
if (minSuccDist <= 0) {
stop("Stopping criterion incorrect.")
}
if (maxNoIt < 1) {
stop("Maximum number of iterations smaller than 1.")
}
Y <- t(Y)
if (hypothesis == "H0") {
return(RCMmlH0orH1(Y, X, R, shrinkType = shrinkType,
estType = estType, corType = corType, hypothesis = hypothesis))
}
if (hypothesis == "H1") {
return(RCMmlH0orH1(Y, X, R, shrinkType = shrinkType,
estType = estType, corType = corType, hypothesis = hypothesis))
}
if (hypothesis == "H2") {
return(RCMmlH2(Y, X, R, maxNoIt, minSuccDist, shrinkType = shrinkType,
estType = estType, corType = corType))
}
}
<bytecode: 0x0000000012b633b0>
<environment: namespace:sigaR>
--- function search by body ---
Function RCMestimation in namespace sigaR has this body.
----------- END OF FAILURE REPORT --------------
Fatal error: the condition has length > 1
* checking for unstated dependencies in vignettes ... OK
* checking package vignettes in 'inst/doc' ... OK
* checking running R code from vignettes ... SKIPPED
* checking re-building of vignette outputs ... SKIPPED
* checking PDF version of manual ... OK
* DONE
Status: 2 ERRORs, 1 NOTE
See
'C:/Users/biocbuild/bbs-3.11-bioc/meat/sigaR.Rcheck/00check.log'
for details.