[BioC] need help with mixed effects model

[Mark W Kimpel]

Having not gotten a reply and having found a list dedicated to just 
these types of questions on mixed models, I am going to repost on
R-SIG-ME.

Mark
Mark W. Kimpel MD  ** Neuroinformatics ** Dept. of Psychiatry
Indiana University School of Medicine

15032 Hunter Court, Westfield, IN  46074

(317) 490-5129 Work, & Mobile & VoiceMail
(317) 204-4202 Home (no voice mail please)

mwkimpel<at>gmail<dot>com

******************************************************************


Mark W Kimpel wrote:
> I would like to look at gene-gene correlations within a multi-factorial,
> mixed effects experiment. Here are the factors, with levels:
> 
> Animals: 12 (6 in each strain, see below) note the animals are numbered
> 1:12, but I have tried 1:6 for each strain and it did not address my
> problem. I still wonder, however, how they should be numbered or if it
> makes a difference.
> 
> Tissues: 3
> Animal Strain: 2
> 
> I thus have 6*3*2 = 36 samples, from a total of 12 animals.
> 
> Tissue and strain are fixed effects and, for other purposes, of some
> interest, the animals are random effects and not of interest. All I want
> to look at now is the correlation of the continuous variable of
> expression between 2 genes. In biological terms, I want to look for
> evidence of co-regulation.
> 
> It looks to me like I should use lmer/lmer2 for this, but I have not set
> up a model with this function before. I tried to follow an example in
> "The R Book" by Crawley, pp. 648-50, without success (see my code chunk
> below).
> 
> mod <- lmer(gene.mat[2,] ~  + gene.mat[1,] + (1|Strain/Tissue/Rat))
> 
> I strongly suspect my error is in the model, but in case it might help,
> below is the actual output. Note my model terms above have been
> decrypted. I did check to make sure that the phenoData columns are
> correct, and they are.
> 
> Thanks, Mark
> 
> output:
> Error in pData(exprSet)$Treatment:factor : NA/NaN argument
> In addition: Warning messages:
> 1: In pData(exprSet)$Treatment:factor :
>    numerical expression has 36 elements: only the first used
> 2: In inherits(x, "factor") : NAs introduced by coercion
> 
> Enter a frame number, or 0 to exit
> 
> 1: lmer(gene.mat[2, ] ~ +gene.mat[1, ] + (1 |
> factor(pData(exprSet)$Treatment)
> 2: lmerFactorList(formula, mf, fltype)
> 3: lapply(bars, function(x) eval(substitute(as.factor(fac)[, drop =
> TRUE], lis
> 4: FUN(X[[1]], ...)
> 5: eval(substitute(as.factor(fac)[, drop = TRUE], list(fac = x[[3]])), mf)
> 6: eval(expr, envir, enclos)
> 7:
> as.factor(factor(pData(exprSet)$Rat):(factor(pData(exprSet)$Tissue):(pData(
> 8: is.factor(x)
> 9: inherits(x, "factor")
> 
>  > sessionInfo()
> R version 2.7.0 Under development (unstable) (2008-02-17 r44506)
> x86_64-unknown-linux-gnu
> 
> locale:
> LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C
> 
> attached base packages:
> [1] tools     stats     graphics  grDevices datasets  utils     methods
> [8] base
> 
> other attached packages:
> [1] lme4_0.99875-9       Matrix_0.999375-4    lattice_0.17-6
> [4] affy_1.17.5          preprocessCore_1.1.5 affyio_1.7.12
> [7] Biobase_1.17.13
> 
> loaded via a namespace (and not attached):
> [1] grid_2.7.0  tcltk_2.7.0
>  >
> 
> --
> 
> Mark W. Kimpel MD  ** Neuroinformatics ** Dept. of Psychiatry
> Indiana University School of Medicine
> 
> 15032 Hunter Court, Westfield, IN  46074
> 
> (317) 490-5129 Work, & Mobile & VoiceMail
> (317) 204-4202 Home (no voice mail please)
> 
> mwkimpel<at>gmail<dot>com
> 
> ******************************************************************
>