[R-sig-ME] Question on mixed model design to analyze microarray data

[Tom Wenseleers]

Dear Ben,
Many thanks for your message - after some experimentation I did indeed confirm that the correct specification in lmer was the one you mentioned, i.e. the correct code to get almost identical output as in SAS/glimmix was (as I thought it might be useful to share this for others I attach the datafile, and also the PDF with the original SAS analysis)

microarray=read.csv("microarray_dataset1612_R.csv",colClasses=c(rep("factor",8),"numeric"))
library(lme4)
library(lmerTest)
#options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly")) # not sure if this would be necessary
fit=lmer(log_intensity~(1|array/gene)+(1|array:pin)+(1|array:dip)+dye+trt+gene+dye:gene+trt:gene+pin, data=microarray) 
summary(fit)
anova(fit,ddf="Kenward-Roger")

(the degrees of freedom are slightly different though as in the SAS output, not sure why that is, but the p vals are practically identical)

The only remaining question I now have is how I could obtain the least square means, equivalent to the SAS code
lsmeans trt*gene / slice=gene
slicediff=gene
slicedifftype=control('0' '2');
(cf PDF in attachment)

Would anyone happen to know perhaps how I could obtain these using package lsmeans? (I am new to that package, so apologies if this would happen to be trivial)

Cheers,
Tom


-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Ben Bolker
Sent: 29 June 2013 17:58
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Question on mixed model design to analyze microarray data

Tom Wenseleers <Tom.Wenseleers at ...> writes:

> 
> Dear all,
> I was just trying to analyse some microarray gene expression data 
> using mixed models in R.
> Based on the SAS code given here
> http://support.sas.com/documentation/cdl/en/
> statug/63033/HTML/default/viewer.htm#statug_hpmixed_sect035.htm

 (URL broken, sorry)

>    proc hpmixed data=microarray;
>       class marray dye trt gene pin dip;
>       model log2i = dye trt gene dye*gene trt*gene pin;
>       random marray marray*gene dip(marray) pin*marray;
>       test trt;
>    run;
> 
> I wanted to try to run the equivalent model in R using lmer and then 
> test significance using lmerTest and/or lsmeans. Am I correct that the 
> correct syntax would be 
> lmer(log2i~dye+trt+gene+dye:gene+trt:gene+pin+(1|marray)+
> (gene|marray)+(1|marray/dip)+(pin|marray), data=microarray) (gene has 
> about 10 000 levels and marray 10, would this run in lmer?)

  I think you probably want a random effect of 

(1|marray) + (1|marray:gene) + (1|marray:pin)+(1|marray:dip)

or equivalently

(1|marray/gene) + (1|marray:pin) + (1|marray:dip)

 you almost certainly _don't_ want (gene|marray), which will try to fit correlations among gene effects by microarray (i.e a 10K by 10K correlation matrix), which will explode your computer.  At a quick read I'm not quite sure what distinction the SAS authors are making between "X by Y" and "X within Y": I would code "X by Y" as
(1|X) + (1|Y) + (1|X:Y)  [ (1|X*Y) might work, but I haven't tried it ] and "X within Y" as (1|Y/X) or (1|Y) + (1|Y:X)
> 
> Also, would it be possible in any way to allow variances across 
> treatment (trt) groups to be unequal? Or is that only possible in 
> nlme/lme?
> (but that wouldn't allow for such complex random effect designs, 
> right?)

  Not at present in lmer.  You might be able to code these random effects in nlme (see the relevant pages of Pinheiro & Bates 2000, referenced in http://glmm.wikidot.com/faq under "Crossed random effects", but lme4 is likely to be a lot faster.

_______________________________________________
R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-------------- next part --------------
A non-text attachment was scrubbed...
Name: littell et al 2006_microarray example.pdf
Type: application/pdf
Size: 148460 bytes
Desc: littell et al 2006_microarray example.pdf
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20130630/e9c2cba1/attachment-0001.pdf>