[BioC] covariate information

W. Evan Johnson wej at bu.edu
Fri Oct 26 15:04:12 CEST 2012 

Chol-hee, 

So I think I now better understand your issue. Why do you have individual ID in the model as a covariate? Do you have repeated measures in your dataset? If not, you should NOT be including individual ID in as a covariate. If it is a repeated measures experiment, then is seems that individual is nested within your treatments/covariates (do you only have one individual ID per treatment combination?). Also, note that you should be treating age as a numerical covariate. 

I would strongly advise to you that you go meet with a statistician with experience in linear models and discuss with her/him what covariates should be included in your model (and why!). It seems that your question is more than just a ComBat question but rather a design and linear modeling question. Any good statistician should be able to help you with this.

Thanks!

Evan



On Oct 26, 2012, at 12:46 AM, Cholhee Jung wrote:

> Dear Evan,
> 
> Thank you for your reply.
> 
> So, my understanding is that you suggest my first attempt as the more
> appropriate way to run ComBat than the second attempt.
> However, removing covariate one after another until I don't get the
> singularity error spared just one covariate.
> 
> While the only remaining covariate is the most important one (as it
> is the individual ID), I wonder if it's really OK to lose all other
> covariate information, such as sample preparation methods, age,
> disease etc, for the technical artefact removal.
> 
> I'm just worried about losing actual signals coming up between these
> covariates by ComBat.
> 
> Regards,
> Chol-hee
> 
> 
> On Fri, Oct 26, 2012 at 12:19 AM, W. Evan Johnson <wej at bu.edu> wrote:
>> Chol-hee,
>> 
>> Notice the simple example:
>> 
>>> x=as.factor(c(1,0,1,0));y=as.factor(c(1,2,1,2));z=rnorm(4)
>> 
>> Notice that x and y are the same covariate. Now:
>> 
>>> design=model.matrix(z~x+y)
>>> design
>>  (Intercept) x1 y2
>> 1           1  1  0
>> 2           1  0  1
>> 3           1  1  0
>> 4           1  0  1
>> attr(,"assign")
>> [1] 0 1 2
>> attr(,"contrasts")
>> attr(,"contrasts")$x
>> [1] "contr.treatment"
>> 
>> attr(,"contrasts")$y
>> [1] "contr.treatment"
>> 
>>> solve(t(design)%*%design)
>> Error in solve.default(t(design) %*% design) :
>>  Lapack routine dgesv: system is exactly singular
>> 
>> You get the singularity error because your covariates are exactly the same
>> (or not linearly independent). Now if you concatenate the variables like you
>> did:
>> 
>>> xy=as.factor(paste(x,y,sep='.'))
>>> xy
>> [1] 1.1 0.2 1.1 0.2
>> Levels: 0.2 1.1
>> 
>> Which is clearly different than the original x+y. Now:
>> 
>>> design=model.matrix(z~xy)
>>> design
>>  (Intercept) xy1.1
>> 1           1     1
>> 2           1     0
>> 3           1     1
>> 4           1     0
>> attr(,"assign")
>> [1] 0 1
>> attr(,"contrasts")
>> attr(,"contrasts")$xy
>> [1] "contr.treatment"
>> 
>>> solve(t(design)%*%design)
>>            (Intercept) xy1.1
>> (Intercept)         0.5  -0.5
>> xy1.1              -0.5   1.0
>> 
>> Which now works. However, note that I wouldn't use the latter. You need to
>> find out which of your six covariates are linearly dependent with each other
>> and remove one or more so they are NOT linearly dependent. This will be
>> different from your second attempt but will be equivalent to what you were
>> trying to accomplish in your first attempt
>> 
>> Let me know if this doesn't work!
>> 
>> Thanks!
>> 
>> Evan
>> 
>> --
>> W. Evan Johnson
>> Assistant Professor
>> Division of Computational Biomedicine
>> Boston University School of Medicine
>> 72 East Concord St., E-645
>> Boston, MA 02118
>> Phone: (617) 638-2541
>> 
>> 
>> 
>> On Oct 25, 2012, at 6:00 AM, bioconductor-request at r-project.org wrote:
>> 
>> ------------------------------
>> 
>> Message: 8
>> Date: Wed, 24 Oct 2012 13:40:33 -0400
>> From: Achilleas Pitsillides <anp4r at virginia.edu>
>> To: Cholhee Jung <jung.cholhee at gmail.com>, Bioconductor mailing list
>> <bioconductor at stat.math.ethz.ch>
>> Subject: Re: [BioC] Fwd: covariate information
>> Message-ID:
>> <CABDy-6=dP5-Wu+Zqo4-UpMBPX51zeCRPQCwdSV5yGYHg2Z821A at mail.gmail.com>
>> Content-Type: text/plain
>> 
>> Dear Chol-hee,
>> 
>> The short answer is that the two model matrices are different and they have
>> different dimensions; you can verify this by using the dim(mod_mat) to see
>> the dimension of the model matrix.
>> 
>> Here is my understanding: If you have a factor  f1 with 3 levels and a
>> factor f2 with 2 levels (where all the possible level combinations exist),
>> then f1:f2 is all the combinations of f1 and f2 (an equivalent  factor
>> with 6 levels) and  the  model matrix ~f1:f2 would have six columns (i.e.
>> fit a model with 6 coefficients).
>> However, the model matrix ~f1+f2 will have 4 columns ( i.e. fit 4
>> coefficients: constant, two for f1 and one for f2).
>> The model ~f1*f2 will fit 6 coefficients and have a model matrix with the
>> same column space as the model matrix for ~f1:f2.
>> 
>> 
>> I hope this helps,
>> 
>> cheers,
>> Achilleas
>> 
>> 
>> On Tue, Oct 23, 2012 at 10:10 PM, Cholhee Jung
>> <jung.cholhee at gmail.com>wrote:
>> 
>> 
>> 
>> Dear users,
>> 
>> 
>> Below is the question I posted originally in the 'ComBat user forum' of
>> 
>> Google group.
>> 
>> But, as I was suggested to forward my question to Bioconductor mailing
>> 
>> list, I'm doing it now.
>> 
>> 
>> Please find my question, below.
>> 
>> 
>> 
>> Regards,
>> 
>> Chol-hee
>> 
>> 
>> On Tuesday, October 23, 2012 10:13:26 AM UTC+11, Cholhee Jung wrote:
>> 
>> 
>> 
>> 
>> Dear users,
>> 
>> 
>> I was trying ComBat on ~1,000 samples.
>> 
>> Samples are spread over 12 batches and each batch contains 4 technical
>> 
>> replicates that are identical across all batches.
>> 
>> The number of covariates is 5, and I was using the ComBat implemented in
>> 
>> the 'sva' package.
>> 
>> 
>> I tried ComBat with two model matrix built from the same covariate
>> 
>> information.
>> 
>> 
>> First model matrix was constructed as below:
>> 
>> mod_mat = model.matrix(~as.factor(cov1) + as.factor(cov2) +
>> 
>> as.factor(cov3) + as.factor(cov4) + as.factor(cov5), data=pheno_data )
>> 
>> 
>> Second one was built as below:
>> 
>> mod_mat = model.matrix(~as.factor(paste(pheno_data$cov1,
>> 
>> pheno_data$cov2, pheno_data$cov3, pheno_data$cov4, pheno_data$cov5,
>> 
>> sep=":")))
>> 
>> 
>> Basically, covariates were concatenated into one string for the the
>> 
>> second model matrix.
>> 
>> 
>> ComBat with the first model matrix raised the 'singular' error like
>> 
>> below:
>> 
>> 
>> Error in solve.default(t(design) %*% design) :
>> 
>>  Lapack routine dgesv: system is exactly singular
>> 
>> 
>> But, ComBat run without error with the second model matrix.
>> 
>> 
>> 
>> Now I wonder if the two different model matrices are same?
>> 
>> 
>> Regards,
>> 
>> Chol-hee
>> 
>> 
>> 
>> _______________________________________________
>> 
>> Bioconductor mailing list
>> 
>> Bioconductor at r-project.org
>> 
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> 
>> Search the archives:
>> 
>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>> 
>> 
>> 
>> [[alternative HTML version deleted]]
>> 
>> 
>

More information about the Bioconductor mailing list