I am designing an experiment with juvenile fishes. There are 12 replicated tanks. In each of the first 6 tanks there is a even mix of male and female fish (100 of each sex in each tank); 3 of these tanks will be treated with chemical A, and 3 treated with chemical B. The remaining 6 tanks will be set up identically but using a different genetic family of fish. The measured response is male and female survival in each tank/treatment.

When constructing a mixed model, I am assuming that this is a simple nested design. Would it be correct to include treatment as a model factor, along with tank, fish sex nested in tank, and genetic ID nested within fish sex? How best to deal with sex and genetic ID is what I am most uncertain about. Any suggestions for coding this in R?

  • 1
    $\begingroup$ What is your main hypothesis? Are you trying to test the interaction between genetic family and treatment options? Also, from your description, it seems like gender is 50/50 split for all 12 tanks. Is this true? $\endgroup$
    – Mark White
    Jul 7 '17 at 13:12
  • $\begingroup$ Hi there - the prediction is that chemical treatment differentially affects survival of male and female fish. But I am also exploring whether different genetic lines respond the same way. I should point out that survival is measured as time alive - a continuous variable, not a binomial dead or alive. And yes, sex is split 50/50 in all 12 tanks. Thanks. $\endgroup$
    – GTC
    Jul 7 '17 at 13:20
  • $\begingroup$ Also - I will code in R. $\endgroup$
    – GTC
    Jul 7 '17 at 13:21

Structure of data

Your data are nested, and here is how I would code it, given your research question:

Level 2: You have 12 different tanks, each of which should have a unique ID number coded as factor (I will code this as tankID). At the tank level, you have two independent variables: genetic family (I will code it as genfam, which is a dichotomous factor) and treatment (I will code it as treat, which is a dichotomous factor).

Level 1: You have $n$ fish, each of which should have a unique ID number coded as factor (I will code this as fishID). You have one independent variable at this level: sex (which I will code as sex, a dichotomous factor). Your dependent variable is at this level (as is the case with multilevel models generally): survival (which I will code as survival, which I assume is either int, num, or dbl).

Your data should look something like:

fishID  tankID  genfam  treat  sex  survival
1       1       A       A      M    12
2       2       A       B      F    8
...     ...     ...     ...    ...  ...

...and so on.

Model code

I will skip the mathematical/Greek notation, unless you are into that. Since you said you are coding in R, I'll go right to how to code it. I would use the lme4 and lmerTest packages to run your model. I will code your data as data. The saturated model would be:

lmer(survival ~ sex * treat * genfam + (1 + sex | tankID), data=data)

What this is saying is that survival is predicted by a three-way interaction between your independent variables of interest (R will automatically fill in all main effects and two-way interactions for you). The stuff in parentheses are the random parts of the model. 1 specifies that there is a different intercept for each tank—called a "random intercept." The sex here means that the effect of sex is allowed to differ across different tanks as well—called a "random slope." The pipe | means that these effects are nested within tank, called tankID. I should also note that genfam and treat are not included here, because they are at Level 2, so the effects of genfam and treat cannot vary by Level 2 clusters.

lme4 does not like testing hypotheses with p-values and t-tests, especially since the degrees of freedom for multilevel models aren't straightforward, so the package lmerTest will automatically add degrees of freedom approximations and p-values. Another option is running likelihood ratio tests with the anova function and using nested models, but with so many predictors in your model, this would likely be too much of a hassle (and I find that the likelihood ratio and lmerTest rarely disagree).

  • $\begingroup$ +1. Given that everything is balanced here, one should be able to set up a classical ANOVA, but it's of course much easier with lmer. $\endgroup$
    – amoeba
    Jul 7 '17 at 13:48
  • $\begingroup$ Yes, it could be done with ANOVA, but I am partial to a regression approach, especially in R, which does not have the clearest ANOVA functions (and I do not particularly like the ez package). $\endgroup$
    – Mark White
    Jul 7 '17 at 13:54
  • $\begingroup$ Thank you very much indeed - I appreciate your input. I will code this and run it. I will let you know how it turns out. Cheers... $\endgroup$
    – GTC
    Jul 7 '17 at 13:58
  • $\begingroup$ @GTC Great! I'd be glad to look at your R output or if you have any other questions. $\endgroup$
    – Mark White
    Jul 7 '17 at 14:08
  • $\begingroup$ Cheers mate. Will do. GTC $\endgroup$
    – GTC
    Jul 7 '17 at 14:17

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.