# What kind of statistical method do I need for comparing element concentrations?

I have assayed the concentrations of various metal elements in different anatomic regions of two strains of mice. Now, for each element, in each region, I want to find whether there is any difference between the two strains.

The concentrations of the different metals within an animal are correlated, it is said that doing pairwise t tests will almost certainly result in nonsense P values.

What kind of statistical method do I need here? Preferably, please suggest in the context of the R statistical language.

Here is what I did, perhaps incorrectly (using simulated data as an example, in reality I have 7 mice in each group, 6 anatomic regions for each mouse, 13 elements for each region):

# create the data frame
> elemconc = data.frame(expand.grid(id=1:3, geno=c('exp', 'wt'), region=c('brain', 'spine'), elem=c('fe', 'cu', 'zn')), conc=rnorm(36, 10))
> elemconc
id geno region elem      conc
1   1  exp  brain   fe  8.497498
2   2  exp  brain   fe  9.280944
3   3  exp  brain   fe  9.726271
4   1   wt  brain   fe 11.556397
5   2   wt  brain   fe 10.992550
6   3   wt  brain   fe  9.711200
7   1  exp  spine   fe 11.168603
8   2  exp  spine   fe  9.331127
9   3  exp  spine   fe 11.048226
10  1   wt  spine   fe  8.480867
11  2   wt  spine   fe  8.887062
12  3   wt  spine   fe  8.329797
13  1  exp  brain   cu 10.242652
14  2  exp  brain   cu  9.865984
15  3  exp  brain   cu  9.163728
16  1   wt  brain   cu 11.099385
17  2   wt  brain   cu  9.364261
18  3   wt  brain   cu  9.718322
19  1  exp  spine   cu 10.720157
20  2  exp  spine   cu 11.505430
21  3  exp  spine   cu  9.499359
22  1   wt  spine   cu  9.855950
23  2   wt  spine   cu 10.120489
24  3   wt  spine   cu  9.526252
25  1  exp  brain   zn  9.736196
26  2  exp  brain   zn 11.938710
27  3  exp  brain   zn  9.668625
28  1   wt  brain   zn  9.961574
29  2   wt  brain   zn 10.461621
30  3   wt  brain   zn  9.873667
31  1  exp  spine   zn  9.708067
32  2  exp  spine   zn 10.109309
33  3  exp  spine   zn 10.973387
34  1   wt  spine   zn  8.406536
35  2   wt  spine   zn  7.797746
36  3   wt  spine   zn 11.127984

# use tapply to aggregate
> tapply(elemconc$conc, elemconc[c('elem', 'region')], function(x) x) region elem brain spine fe Numeric,6 Numeric,6 cu Numeric,6 Numeric,6 zn Numeric,6 Numeric,6 # check whether the order of data has been preserved after aggregation > x['fe', 'brain'] [[1]] [1] 8.497498 9.280944 9.726271 11.556397 10.992550 9.711200 # create an external factor for strain grouping > tmpgeno = rep(c('exp', 'wt'), each=3) > tmpgeno [1] "exp" "exp" "exp" "wt" "wt" "wt" # do the t test using the grouping factor > x = tapply(elemconc$conc, elemconc[c('elem', 'region')], function(x) t.test(x~tmpgeno) )
> x
region
elem brain  spine
fe List,9 List,9
cu List,9 List,9
zn List,9 List,9
• Roughly how many elements, and how many mice? – onestop Dec 4 '11 at 16:46
• Please see the edit. – qed Dec 4 '11 at 17:09
• I assume ANOVA would help, but I'm unsure of how to build the model, something like conc~geno+region+elem, random=id (which means relate conc to geno, region and elem, and allow for the random factor id, I think !). Check out personality-project.org/r/r.anova.html and cran.r-project.org/doc/contrib/Faraway-PRA.pdf – PaulHurleyuk Dec 4 '11 at 21:14
• I've never done this myself but I think its possible to run a pair-wise t test controlled by some covariate, in this case mouse speciment. – Dan Dec 6 '11 at 22:57

Yes, this is an analysis of variance problem. Check out the aov function. You might want something like this ("geno" is the thing you want to test, right?):

> model1 <- aov(conc ~ region*elem*geno+Error(id), data=elemconc)
> model.null <- aov(conc ~ region*elem+Error(id), data=elemconc)
> summary(model1)

Error: id
Df Sum Sq Mean Sq F value Pr(>F)
Residuals  1  0.312   0.312

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
region            1   0.04   0.039    0.04   0.84
elem              2   1.70   0.848    0.87   0.43
geno              1   0.02   0.020    0.02   0.89
region:elem       2   3.23   1.614    1.66   0.21
region:geno       1   0.86   0.862    0.89   0.36
elem:geno         2   3.56   1.779    1.83   0.18
region:elem:geno  2   2.32   1.158    1.19   0.32
Residuals        23  22.37   0.973
> summary(model.null)

Error: id
Df Sum Sq Mean Sq F value Pr(>F)
Residuals  1  0.312   0.312

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
region       1   0.04   0.039    0.04   0.84
elem         2   1.70   0.848    0.84   0.44
region:elem  2   3.23   1.614    1.61   0.22
Residuals   29  29.13   1.005
> pf(((29.13-22.37)/6) / (29.13/29),6,29, lower.tail=F)
[1] 0.3744
>