Thank you for your response, Dennis!

I have gone through some literature and examples and think I have chosen correct model structure(s). My data is structurally similar to much used data sets on student scores in different schools, i.e. a multilevel or hierarchical structure. E.g. in Faraway's Extending Linear Model in R (page 174-181), there is an example:

mmodc <- lmer(math ~craven*social + schraven*social + (1|school) + (1|school:class))

where math is a test score, craven is a centered Raven score and social is father's social class, while schraven is the average Raven score of the school (=compositional effect).

In my data set, stations correspond to schools, individual fish to students, math score to weight, age to social class, and density to schraven. Ignoring the class in school random effect, this model is similar to m1.

In the book example, the factor 'social' and both interactions have df of 926, and not < 50 as I would have expected based on the number of schools (=50). Hence, I realize that it is my understanding of df's in relation to mixed effects models that are incorrect. I assume I have to dig into the formulas...

Regards,

Eli


From: Dennis Murphy [mailto:djmuser@gmail.com]
Sent: 10. mars 2010 16:01
To: Kvingedal, Eli
Subject: Re: [R-sig-ME] mixed effects models and pseudo replication

Hi:

See inline...
On Wed, Mar 10, 2010 at 6:07 AM, Kvingedal, Eli <Eli.Kvingedal@nina.no<mailto:Eli.Kvingedal@nina.no>> wrote:
Hi,

I am analysing effects of local population density on fish performance (e.g. weight). My dataset is based on fish sampled from different sites (17 stations) and in addition to measures on individual performance, I have information on age (0 and 1). On site level, I have information on fish densities for both age groups. I am interesting in estimating the effects of fish density on performance and particularly interested in determining possible differences between age groups in the density response.

Traditionally, these kind of data are analysed based on mean values (ancovas). However, based on mixed effects model, the among individual variance will be included in the analysis and not just averaged out. I started by using lmer (lme4 package), but realizing that the variance is increasing with density, I switched to lme (nlme package) and applied variance structures.

My starting model is thus:

m1 <- lme(weight ~ age*density0 + age*density1, random = ~1|station, weights=....)

with station and age as factors.

Now, my issue is pseudo-replication. The summary table shows that the factors age and age*density have very high degrees of freedom (~700) and accordingly low p-values. It seems to me like age and the interactions between age and density are analysed as if the samples were independent, and if so, it means pseudo-replication, doesn't it?

>From what I can tell, the experimental units are sites and fish are the
observational units. A similar situation arises when trees are experimental
units but measures are taken on the apples (or pick your choice of abundant fruit)
from a tree. This is a design feature, not a model feature. Is it reasonable to assume
that fish are nested within site? If so, is there a unique ID for the individual fish or
were measurements taken in bulk? These considerations, to the extent that they
exist in your data, should somehow be represented in the model.

If I set up an alternative random structure allowing for random variance between age classes within station:
m2 <- lme(weight ~ age*density0 + age*density1, random = ~1|station/age, weights=....)

the summary table is more like I think it should be: 14 df for all fixed effects parameters and interactions, and the p-values seem more realistic.

Age is an attribute of an individual fish, right? If so, it should be treated as a fixed
factor. I don't see how age is nested within station, but I'm not familiar with your
data. Perhaps fishID (if available) should be nested within station rather than age?


When comparing m1 and m2 (REML estimation), however, m2 do not provide better fit, and based on literature (e.g. Zuur et al. 2009), then I should use m1.

You start by representing the correct structure of your data in the model (treatments,
blocks, crossing/nesting, etc.) and then adjust from there (e.g., testing different
covariance models or adding/deleting terms). If one model that appears scientifically
unreasonable is no different in fit quality from one that is, which one would you choose?


Testing the significance of the interaction terms by model comparisons (which is what I do to find the optimal model), the significance levels of the likelihood ratio test for specific interaction terms are equivalent whether I use station or station/age as random factors. Which is sort of comforting.

Pseudo-replication (aka 'duplication' in some other fields) generally produces false
significance of effects due to the inflated power arising from additional 'error' df.
In linear models, it uses the wrong error term (sampling error rather than
experimental error). The software doesn't know the correct error term - it deduces
it from the model it is fed - so it's your responsibility to provide the appropriate
model. It seems reasonable to me that the fixed effects tests on density be made
with respect to site as the error term, since density is dependent on site.
I presume that age is a property of an individual fish; if so, then fish need to be
considered as a sub-unit nested within site, assuming that the individual fish
can be uniquely identified in the data set. Age might be considered a 'split-plot'
type of treatment with fish as its appropriate unit.

So, my question is, do I really control for pseudo-replication in the estimation of all fixed effects and interactions when using m1? If so, why these high dfs in the summary table??

If at all possible, yes. The high df in the first model appear to derive from a
failure to distinguish between site as the unit of analysis for density and fish as
the unit of analysis for age. It might be worth spending a little time deciphering
the structure of your data with respect to the effects you want to compare and
then define a model that reflects that structure.

HTH,
Dennis

I would really appreciate if someone could enlighten me!

Regards,

Eli


________________________________________________________________

Eli Kvingedal
PhD Student

Norwegian Institute for Nature Research - NINA
Postal address: NO-7485 Trondheim, NORWAY
Delivery/Visiting address: Tungasletta 2, NO-7047 Trondheim, NORWAY
Phone: +47 73 80 14 00 * Fax: +47 73 80 14 01 * www.nina.no<http://www.nina.no>

_______________________________________________
R-sig-mixed-models@r-project.org<mailto:R-sig-mixed-models@r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


	[[alternative HTML version deleted]]


