[BioC] need help with mixed effects model

[Mark W Kimpel]

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