4
$\begingroup$

I would like to fit a mixed effect model to the following dataset, but I am having difficulties figuring out the best way to define the random effects.

For each subjects (ratID, N=10), I measure the same variable cc_marg using six different instruments (the six levels of the factor mPair) for each cycle of observations (cycles $c=1...N_c$, $N_c>>N$, $N_c$ different across subjects). For each cycle, the six measurements are taken simultaneously (therefore these measurements are correlated by the cycle at which they are taken). For each subject, I repeat this experiment three times (in random order across subjects), one for each level of the controlled variable spd_des (factor with three levels). I am interested in studying the effect of spd_des and mPair (and their possible interaction) on the variable cc_marg. I am not interested in the effect of cycle on the output variable.

There are two sources of randomness: ratID and cycles. However, I am confused on how to nest the latter into the former. There are several cycles for each subject, which would make me think that I should simply need to do ~1|ratID/cycle. However, the cycles obtained at a given level of spd_des (for each subject) are unrelated to those obtained at another level (even though they have the same identifiers $c=1...N_c$). Should I then nest cycle within spd_des within ratID, using ~1|ratID/spd_des/cycle? If I do so, however, I am also defining a random effect of spd_des, which I wasn't planning on doing actually. How do you think I should define the random effects in this design? (this is my main question).

If I do not nest cycle, I obtain unreasonably high numbers of denominator DF when I run anova, increasing the probability of false positive results. Here are the results if I do not nest:

> linM3 <- lme(cc_marg ~ mPair*spd_des , random = ~1|ratID, data=dat_trf, na.action=na.omit, method = "ML", control=lCtr )
> anova.lme(linM3,type="marginal")

              numDF denDF   F-value p-value
(Intercept)       1 14540  128.5679  <.0001
mPair             5 14540 2405.9828  <.0001
spd_des           2 14540    5.4406  0.0043
mPair:spd_des    10 14540   42.7502  <.0001

If I nest cycle within ratID, I obtain:

> linM3n <- lme(cc_marg ~ mPair*spd_des , random = ~1|ratID/cycle, data=dat_trf, na.action=na.omit, method = "ML", control=lCtr )
> anova.lme(linM3n,type="marginal")

              numDF denDF   F-value p-value
(Intercept)       1 12843  128.7659  <.0001
mPair             5 12843 2563.1850  <.0001
spd_des           2 12843    5.0572  0.0064
mPair:spd_des    10 12843   43.9206  <.0001

If I nest cycle within spd_des within ratID, I obtain:

> linM3n4 <- lme(cc_marg ~ mPair*spd_des , random = ~1|ratID/spd_des/cycle, data=dat_trf, na.action=na.omit, method = "ML", control=lCtr )
> anova.lme(linM3n4,type="marginal")

              numDF denDF   F-value p-value
(Intercept)       1 11503  120.7824  <.0001
mPair             5 11503 2803.9750  <.0001
spd_des           2    15    0.8420  0.4502
mPair:spd_des    10 11503   35.2944  <.0001

The results between the first and the second model above are not very different, but the third model provides different results in terms of spd_des. It is therefore important choosing the right model. How should I define the random effects considering the experimental design and the research question? Thanks!

[UPDATE]

I have tried to create a variable 'exp_spd' that keeps track of the experimental session. As stated, there is one experimental session for each level of spd_des, but the order in which I performed the experiments is randomized across subjects. The model is as follows:

linM1n <- lme(cc_marg ~ mPair*spd_des , random = ~1|ratID/exp_spd/cycle, data=dat_trf, na.action=na.omit, method = "ML", control=lCtr )
anova.lme(linM1n,type="marginal")
              numDF denDF   F-value p-value
(Intercept)       1 11528  122.3557  <.0001
mPair             5 11528 2802.2565  <.0001
spd_des           2 11528    0.3990   0.671
mPair:spd_des    10 11528   35.1272  <.0001

In terms of significance of the fixed-effect, the results are equivalent to the model linM3n4 above, where I nested spd_des instead of exp_spd. However, the denDF are different. In particular, the denDF of spd_des changes drastically. Shouldn't it be the same, given that each level of exp_spd is associated to only one level of spd_des (and vice-versa)? This issue is very obscure to me, and any help is highly appreciated.

$\endgroup$

1 Answer 1

3
$\begingroup$

It seems to me that the random grouping factors in your study are:

  • rat;
  • experiment;
  • cycle.

So the "random" portion of your model would be listed as random = ~1|ratID/experimentID/cycleID.

The predictor variable spd_des is an experiment-level predictor which can appear in the fixed effects portion of your model.

Usually, people advise that you should have at least 5 levels for a random grouping variable in your model, so that you can meaningfully estimate the standard deviation (or variance) of the random effects associated with that random grouping variable. However, you only have 3 such levels for the experiment variable. If this poses an issue for the estimation of the standard deviation (or variance) of the random intercepts associated with the experiment variable, you will need to respecify your model.

$\endgroup$
10
  • 1
    $\begingroup$ Thanks! I do not really have an experimentID variable. I repeat the measurements three times to assess the three levels of spd_des. So basically, for each "experiment" (group of cycles), the value of spd_des will be different. I could create an experimentID variable, but wouldn't that be equivalen of doing random = ~1|ratID/spd_des/cycle? $\endgroup$
    – Cristiano
    Commented May 17, 2019 at 15:22
  • 1
    $\begingroup$ No, it wouldn't be equivalent, because the experimentID would be taking values like "1", "2", "3", and appear in the random portion of your model to keep track of which experiment you are dealing with, whereas the spd_des could take the same values, though these would have a totally different meaning, and would be listed in the fixed effects portion of your model (it's a predictor whose value you measure for each experiment). There is a big difference between a random grouping factor and a predictor variable you are measuring for each level of that random grouping factor! $\endgroup$ Commented May 17, 2019 at 15:27
  • 1
    $\begingroup$ I would suggest to number the experiments for each rat in the exact order in which you performed them - so that 1 corresponds to the first experiment for that rat, 2 corresponds to the second experiment, etc. It was my understanding that you had 3 experiments per rat? But that each experiment included multiple cycles? That's what your initial post seemed to suggest. $\endgroup$ Commented May 17, 2019 at 15:51
  • 1
    $\begingroup$ ok, sorry I misunderstood. Yes, you are right: there are three experiments for subject (one for each level of spd_des), and within each experiment there are multiple cycles. In truth, there are 6 experiments per rat; three of them are to assess the effect of spd_des, and the other three to assess three levels of another variable inc_des. I was going to study spd_des and inc_des in two different models because I am not evaluating all the possible combinations and I am not interested in their interactions. The experiments are randomized across the level of both variables though. $\endgroup$
    – Cristiano
    Commented May 17, 2019 at 16:31
  • 1
    $\begingroup$ The names you use for the experiments are just labels but these labels should have a clear meaning - whether you call an experiment 4, D or IV for a rat, you should be able to describe what differentiates that experiment from the other experiments you conducted for that rat. The label is something that helps you distinguish between experiments, while also characterizing the nature of each experiment. $\endgroup$ Commented May 17, 2019 at 17:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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