Backwash throws errors with more than one covariate of interest

[scottgigante]

I am having trouble running backwash with length(cov_of_interest) > 1. With 102 samples of gene expression data in 18 groups of affymetrix arrays, I get the following error:

ematrix <- load("ematrix.RData")
groups <- as.factor(sub("_[_0-9\\.]*$", "", rownames(ematrix)))
X <- model.matrix(~groups)
mout <- backwash(ematrix, X, cov_of_interest=2:3, scale_var=TRUE)

Number of confounders not provided so being estimated with package sva.
longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length
Error in eigen(crossprod(alpha_tilde, alpha_tilde), symmetric = TRUE) : infinite or missing values in 'x'
4.
stop("infinite or missing values in 'x'")
3.
eigen(crossprod(alpha_tilde, alpha_tilde), symmetric = TRUE)
2.
backwash_second_step(betahat_ols = betahat_ols_star, S_diag = S_diag_star, alpha_tilde = alpha_tilde_star, tau2_seq = tau2_seq, lambda_seq = lambda_seq, pi_init_type = pi_init_type, scale_var = scale_var, sprop = sprop, var_inflate_pen = var_inflate_pen)
1.
backwash(ematrix, X, cov_of_interest = 2:3, scale_var = TRUE)

Expression matrix is here: ematrix.zip

I tried simulating some simpler data to recreate the issue and I got another issue.

Y <- matrix(rnbinom(n=100*20000, size=200, mu=1000), ncol=20000)
groups <- as.factor(rep(1:10, times=10))
X <- model.matrix(~groups)
mout <- backwash(Y, X, cov_of_interest=2:3, scale_var=TRUE)

Number of confounders not provided so being estimated with package sva.

Error: alpha_tilde is not a matrix
4.
stop(assertError(attr(res, "msg")))
3.
assertthat::assert_that(is.matrix(alpha_tilde))
2.
backwash_second_step(betahat_ols = betahat_ols_star, S_diag = S_diag_star, alpha_tilde = alpha_tilde_star, tau2_seq = tau2_seq, lambda_seq = lambda_seq, pi_init_type = pi_init_type, scale_var = scale_var, sprop = sprop, var_inflate_pen = var_inflate_pen)
1.
backwash(Y, X, cov_of_interest = 2:3, scale_var = TRUE)

Any ideas how I can use mouthwash or backwash to estimate effects of more than one covariate?

Thanks!

[dcgerard]

Thanks for trying out vicar!

Unfortunately, vicar doesn't support having more than one covariate of interest yet. The documentation did say you could put in a vector, but I was inheriting that documentation from another function vruv4. I've modified vicar so the documentation is correct and so that the error message is more informative.

It's not ideal, and I have no idea about how it will work, but you could try running mouthwash or backwash one time for each covariate of interest, while still controlling for the other covariates.

The second issue occurs because the package estimated that there were no confounders in the simulated data. Running mouthwash without any confounders is not supported so I made the error message in this situation explicit. If k is estimated to be 0 in real data, then I would recommend trying out https://github.com/stephens999/ashr

[scottgigante]

Oh, that makes a lot of sense. Thanks for clarifying, and separately thank you very much for putting this package together! I'll have a go at looking at the covariates one by one and see what happens.