[R-sig-ME] Analyzing hierarchical MRI data

[Stephen Weigand]

Hello,

I'd like to explore ways to analyze data from brain MRIs using lmer.

My data are typically of the following form:

* There are a small number of groups (e.g., disease 1, disease 2,
control) with from 10 to 50 human subjects in each group.

* Each subject has had an MRI of the head. From the MRI I get one
value at each of 30 different regions of the brain (e.g., a value for
the amygdala, a value for the hippocampus, etc.). This value is
typically an estimate of the region size. So I get 30 values for each
subject. These values are typically conditionally Gaussian. My data
frame may look like this

group subject region y
A Subj1 Amygdala 2.5
A Subj1 Hippocampus 2.8
A Subj2 Amygdala 3.2
A Subj2 Hippocampus 4.8
...
B Subj3 Amygdala 1.7
B Subj3 Hippocampus 4.9
B Subj4 Amygdala 2.2
B Subj4 Hippocampus 3.5
...

The question is, In which regions do the groups differ and by how much?

I can treat group and region as fixed effects by arguing that they are
the only groups and regions of interest and fit a random intercept
model with a group by region interaction of the form:

lmer(y ~ group*region + (1 | subject))

But I would like to take advantage of pooling/penalization/shrinkage
to get more reliable estimates of the differences between groups at
each region.  So I think I want to turn region (and maybe group?) into
random effects.

I want to try

lmer(y ~ (1 | group) + (1 | region) + (1 | subject))

but is that ignoring the nested structure of the data? I would greatly
appreciate suggestions.

Thank you,

Stephen

PS I'm OK with ignoring spatial correlation for now although I expect
that regions that are close to one another in the brain are going to
be more correlated than regions at opposite sides of the brain.

PPS I would guess that to the experts these questions seem all the
same but to the uninitiated, every problem seems like a unique case!!

-- 
Rochester, Minn. USA