I’m having a hard time figuring out what type of model I should use to analyse data with a rather complexe design.

Experimental goal and protocol: The goal is to study the effect of beehive color and inner temperature on bee activity. The experimental layout is: The number of bees at the entrance to 48 beehives was sampled over several months. Beehives were of different colors (4 colors) and the outside temperature of the beehives was controlled (3 temperatures). Because it was not possible to set this up in a random manner, the beehives were 1/ first grouped into three lots (each lot had 16 beehives), with three different outside temperatures (TempOut: 10, 15, 20) ; and 2/ in each lot, the 16 beehives were grouped into 4 Groups (4 beehives in each group), painted in different colors (Color: factor noted 1 to 4). I would call this a nested plot design. In addition, there was a wasp infestation close by, during month 6 of the experiment, causing the number of bees to drop. We consider here that this variable has three modalities (Wasp: "Before" infestation: 5 months of data, "During" infestation: one month of data and "After" infestation: 6 months of data). Observations (number of bees at entrance, inner temperature of beehive noted "TempIn") was done every week during 12 months.

Here’s an illustation, which shows the experimental design.

Experimental setup

We would like to analyse the effect of (1) beehive color "Color" (categorical, 4 levels), (2) beehive inner temperature "TempIn" (continuous), and (3) wasp infestation "Wasp" (categorical, 3 levels) on bee activity "BeeActivity" (continuous).

Questions : (1) What type of statistical model should I test?

There is clearly a temporal autocorrelation in the data, as beehives were sampled multiple times. Also, the set-up is clearly not randomized.

I’m thinking of three solutions:

i. applying a mixed model effect, with a temporal autocorrelation, with three fixed effects (ColorTempInWasp), one random effect (Group) and specifying that data from each beehive is temporally correlated.

The problem I see is that the Groups (12 of them, with 4 beehives in each group) received different treatments (Color*TempIn). This is different from random effects such as « subject » (different patients in a test for example) which all receive the same treatments.

ii. applying a mixed model effect, with a temporal autocorrelation, with three fixed effects (ColorTempInWasp), one random effect (TempOut) and specifying that data from each beehive is temporally correlated. The problem is that TempIn and TempOut will be well correlated.

iii. The other thing I was thinking of was to average bee activity for the 4 beehives in each group, so that I would have 12 average measures of bee activity per sampling date. I could then test a linear model of the type: gls(BeeActivity ∼ ColorTempInWasp, na.action = na.omit, data = Beehive, correlation = corAR1(form =∼ DoY)) The assumption here, I think, would be that these average measures of bee activity randomly received the different treatments. However, I don’t think this works because the hives were still grouped into « lots » with constant temperatures within the « lots ». These kind of look like blocs in a randomized block design, but the blocs are not repeated.

(2) How would the model(s) be correctly written in R?

For the mixed model, I was thinking of something like:

M5 <- lme(BeeActivity ∼ Color*TempIn*Wasp, 
        random= ∼1 | Group, method = "REML", 
        correlation = corAR1(form =∼ DoY), data = Beehive)

Can I set both a random structure and an autocorrelation structure in the same model?

Help much appreciated! Thanks.

Example dataset attached here (partially re-constructed)

  • 1
    $\begingroup$ Don't average. M5 is the right direction but I don't see why you use Group as the grouping factor. I believe you should group by the hives. The auto-correlation term needs the same grouping structure as the random effect: random= ∼1 | HiveID, correlation = corCAR1(form = ~ DoY | HiveID). You should also test random slopes. $\endgroup$ – Roland Jan 4 '18 at 8:06
  • $\begingroup$ Thanks. I agree that there should be a term saying that measurements in each beehive are correlated temporally, so I changed the autocorrelation formulation (HiveID instead of DoY). Or do you actually mean I should test a spatial correlation (all measures on same DoY are correlated) as well as a temporal correlation (all measured from the same beehive are correlated)? I understand that the same grouping structure is required, but I’m not sure how to do that (depending on the formulation here below). $\endgroup$ – Christine H. Jan 4 '18 at 19:53
  • $\begingroup$ Concerning the use of Group as a random variable, the idea was that the beehives are not just randomly distributed, they were first grouped in similar TempOut, then with similar Color conditions. Beehives within each Group are pseudo-replications in the experimental design, and each of the measurements at each of the four pseudo-replicates of each group are not independent samples. Maybe a better formulation would be: M5 <- lme(BeeActivity ∼ 1+Color*TempIn*Wasp, random= ∼1 | BlocID / SubBlocID, method = "REML", correlation = corCAR1(form =∼ HiveID), data = Beehive) $\endgroup$ – Christine H. Jan 4 '18 at 19:53
  • $\begingroup$ BlocID : grouping into « lots » (categorical, 3 levels : temp1, temp2, temp3) ; SubBlocID: grouping into Groups within each « lot » (categorical, 4 levels : col1,col2,col3,col4). In this case: (i) is the formulation of the random intercept effect correct, considering I have two random variables? « ∼1 | BlocID / SubBlocID » $\endgroup$ – Christine H. Jan 4 '18 at 19:54
  • $\begingroup$ (ii) how do I then formulate that I what to test random slopes ? I should test all slopes between fixed and random effects. Should I put something like « ∼1+ ColorTempInWasp | BlocID / SubBlocID » (iii) Doing so, I get many fixed and random variables, and testing all the interactions seems tedious. Using a top-down strategy, I should first find the optimal structure of the random component of the « beyond optimal models » (all fixed variables with their interactions). How would I deconvulate « ∼1+ ColorTempInWasp | BlocID / SubBlocID » and would this be the right thing to do ? Thks agn $\endgroup$ – Christine H. Jan 4 '18 at 19:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.