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.