6
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Roughly how many elements, and how many mice? $\endgroup$ – onestop Dec 4 '11 at 16:46
  • $\begingroup$ Please see the edit. $\endgroup$ – qed Dec 4 '11 at 17:09
  • 1
    $\begingroup$ 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 $\endgroup$ – PaulHurleyuk Dec 4 '11 at 21:14
  • $\begingroup$ 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. $\endgroup$ – Dan Dec 6 '11 at 22:57
1
$\begingroup$

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
> 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.