do not compare bio and batch

[zhou13]

I am not familiar with R but I think even if nlevels are equal, you still cannot compare bio and batch.

    if(nlevels(bio)==nlevels(batch)) {
      if(all(bio==batch)) {
        stop("Biological conditions and batches are confounded. They cannot both be included in the model, please set at least one of 'adjust_bio' and 'adjust_batch' to 'no.'")
      }
    }

[mbcole]

Hey Yichao - you are correct: alone, this code would only stop the user from using the same factors to describe both biology and batch. However, prior to execution of this line, the levels of batch have been re-ordered according to bio factor. If bio and batch are completely confounded, levels for each factor should be in the same order.

[drisso]

I think that turning this into

all(as.numeric(bio)==as.numeric(batch))

should work.

I can try it out later.

[drisso]

This is harder than I thought.

Consider the following to scenarios.

Scenario 1. Two identical factors, but with different levels.

bio <- gl(2, 5)
batch <- as.factor(rep(c("one", "two"), each=5))

> table(bio, batch)
   batch
bio one two
  1   5   0
  2   0   5

In this case, the test in the if statement will throw an error:

> all(bio==batch)
Error in Ops.factor(bio, batch) : level sets of factors are different

My proposed solution would work

> all(as.numeric(bio)==as.numeric(batch))
[1] TRUE

Scenario 2. Two confounded factors with the same levels in a different order.

bio <- gl(2, 5)
batch <- factor(rep(c(1, 2), each=5), levels=2:1)
> table(batch, bio)
     bio
batch 1 2
    2 0 5
    1 5 0

In this case the current code will work but my proposed solution fails.

> all(batch==bio)
[1] TRUE
> all(as.numeric(batch)==as.numeric(bio))
[1] FALSE

Can we think of other scenarios?

[drisso]

Stepping back, what we want to test is really the collinearity of the design matrix.
So I think that we should test directly on the design matrix.

This should work

dm <- model.matrix(~bio+batch)
dm_rank <- qr(dm)$rank
if(NCOL(dm)>dm_rank) {
stop("Biological conditions and batches are confounded. They cannot both be included in the model, please set at least one of 'adjust_bio' and 'adjust_batch' to 'no.'")
}

[drisso]

It turns out that we were already testing for nested bio and batch, if they are nested and have the same number of factors, then they are confounded.