0
$\begingroup$

I am currently writing my master thesis about the effect of an insecticide (clothianidin) on the microflora of bumblebees. I received the bumblebees from an experiment with a nested study design. 16 fields were paired according to land use of the surroundings etc. In each field 2 boxes were placed, which contained 2 hives (colonies) each.

I was trying to determine if the treatment affects the prevalence of certain organisms including ABPV, N.bombi and Snodgrasella as you can see in this data frame:

    structure(list(treatment = structure(c(2L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L), .Label = c("Clothianidin", 
"Control"), class = "factor"), pair = structure(c(1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L), .Label = c("P01", 
"P02", "P03", "P04", "P05", "P10", "P11", "P12"), class = "factor"), 
    field = structure(c(6L, 6L, 6L, 6L, 12L, 12L, 12L, 12L, 1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 10L, 10L, 10L, 10L, 13L, 13L, 
    13L, 13L, 7L, 7L, 7L, 7L, 16L, 16L, 16L, 16L, 8L, 8L, 8L, 
    8L, 9L, 9L, 9L, 9L, 3L, 3L, 3L, 3L, 11L, 11L, 11L, 11L, 4L, 
    4L, 4L, 4L, 5L, 5L, 5L, 5L, 14L, 14L, 14L, 14L, 15L, 15L, 
    15L, 15L), .Label = c("VR02", "VR03", "VR04", "VR05", "VR06", 
    "VR07", "VR09", "VR12", "VR13", "VR14", "VR16", "VR17", "VR18", 
    "VR20", "VR21", "VR23"), class = "factor"), box.nested = c(12, 
    11, 12, 11, 23, 24, 23, 24, 1, 1, 2, 2, 4, 3, 3, 4, 20, 20, 
    19, 19, 25, 26, 25, 26, 14, 14, 13, 13, 31, 31, 32, 32, 16, 
    15, 15, 16, 18, 17, 17, 18, 6, 5, 5, 6, 21, 22, 22, 21, 8, 
    8, 7, 7, 9, 9, 10, 10, 28, 27, 28, 27, 29, 30, 30, 29), hive.nested = c(24L, 
    21L, 23L, 22L, 46L, 48L, 45L, 47L, 2L, 1L, 4L, 3L, 8L, 5L, 
    6L, 7L, 40L, 39L, 38L, 37L, 49L, 52L, 50L, 51L, 27L, 28L, 
    26L, 25L, 62L, 61L, 64L, 63L, 31L, 29L, 30L, 32L, 36L, 34L, 
    33L, 35L, 12L, 9L, 10L, 11L, 41L, 43L, 44L, 42L, 15L, 16L, 
    13L, 14L, 17L, 18L, 20L, 19L, 55L, 54L, 56L, 53L, 58L, 59L, 
    60L, 57L), ABPV.detected = structure(c(0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
    1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0), .Dim = c(64L, 1L), .Dimnames = list(NULL, "ABPV.detected")), 
    N.bombi.detected = structure(c(0, 0, 0, 0, 0, 0, 0, 0, 1, 
    1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), .Dim = c(64L, 
    1L), .Dimnames = list(NULL, "N.bombi.detected")), Snodgrasella.detected = structure(c(0, 
    1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 
    1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1), .Dim = c(64L, 1L), .Dimnames = list(NULL, 
        "Snodgrasella.detected"))), .Names = c("treatment", "pair", 
"field", "box.nested", "hive.nested", "ABPV.detected", "N.bombi.detected", 
"Snodgrasella.detected"), class = "data.frame", row.names = c(NA, 
-64L))

I was trying estimate the treatment effect with a model that included pair, field, box and hive as (nested) random effects:

library(lme4)
ABPV.prev <- glmer(ABPV.detected ~ treatment 
                   + (1|pair/field/box.nested/hive.nested)
                   ,data=data.f, 
                   family=binomial)
summary(ABPV.prev)

The models of ABPV and N. bombi failed to converge, because I have so many zeros. ABPV was only found in one pair and N. bombi was only found in 2 pairs.

Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.100788 (tol = 0.001, component 2)

The p-values of the models did indicate significant treatment effects, but I guess it's pair or field effects that cause the variation...

In another forum I read that it can be tested whether the failure to converge represents a real problem using this function..

relgrad <- with(ABPV.prev@optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad))

... and it does as the p-value is rather large p=0.1 (for ABPV).

I removed some of the random effects and it works and it works e.g. when I remove both hive.nested and box.nested from the model, though increasing the AIC (why?).

I also tried to include pair as a fixed effect:

library(lme4)
ABPV.prev <- glmer(ABPV.detected ~ treatment + pair 
                   + (1|field/box.nested/hive.nested)
                   ,data=data.f, 
                   family=binomial)

but it produced an error:

Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

My questions are

  1. How do I judge which random effects to include and which not?
  2. How do I test if it is a pair (or field effect) not a treatment effect that causes the variation?

N.B. In the Snodgrasella model it makes absolutely no difference if I exclude, any or all of the random effects. Why is that?

Hope you can help me, thanks!

$\endgroup$
1
  • 2
    $\begingroup$ This should be asked on the CrossValidated.com website. (It's not really a coding question but rather a statistical methods question.) $\endgroup$
    – DWin
    Aug 29 '15 at 18:42
2
$\begingroup$

There are too many questions here. It would be a really good idea if you were able to get some local statistical help, as you're operating at the edge of what you can do with generalized linear mixed models. Frank Harrell's rule of thumb is that if you have p "effective" samples (which in your case would be the number of detections for a particular species) you shouldn't be trying to fit a model with more than p/10 parameters (optimistically; the denominator should perhaps even be 20).

colSums(data.f[,6:8])
##   ABPV.detected      N.bombi.detected Snodgrasella.detected 
##               5                     5                    50 

p/10 is less than 1 for the first two cases (ABPV and N. bombi). That means all you can really do with your smaller samples is to collapse your data to the simplest possible case --- a straight-up treatment comparison --- and explain to your audience that it's just not feasible to try to control for non-independence:

ABPV.tab <- with(data.f,table(ABPV.detected,treatment))
##              treatment
## ABPV.detected Clothianidin Control
##             0           31      28
##             1            1       4
chisq.test(ABPV.tab,simulate.p.value=TRUE)
## data:  ABPV.tab
## X-squared = 1.9525, df = NA, p-value = 0.3603

You'd get a similar answer from prop.test(x=c(1,4),n=c(32,32)).

In general, people are pretty ready to let you simplify your model when you give them a "negative" (p>0.05) answer, as they figure you're not trying to get away with something. What you should not do is accept the null hypothesis (make strong claims that the treatment does not have an effect): prop.test tells you

95 percent confidence interval:
-0.25447656  0.06697656

i.e. the 95% CI on the difference of the probability of detection is between -25% and +6% (you might be anywhere from 25% less likely to 6% more likely to detect ABPV in the Clothianidin treatment).

There are a number of reasons you shouldn't believe the results from the GLMMs in these cases -- even if the model converged correctly, inference from GLMMs is most reliable when the data sets are "large" in some sense.

N.B. In the Snodgrasella model it makes absolutely no difference if I exclude, any or all of the random effects. Why is that?

Because most (not all) of the random effects are exactly or nearly zero.

Snodgrasella.prev <- glmer(Snodgrasella.detected ~ treatment 
               + (1|pair/field/box.nested/hive.nested)
               ,data=data.f, 
               family=binomial)
VarCorr(Snodgrasella.prev)   
##  Groups                                Name        Std.Dev.  
##  hive.nested:(box.nested:(field:pair)) (Intercept) 0.00047247
##  box.nested:(field:pair)               (Intercept) 0.00000000
##  field:pair                            (Intercept) 1.08244578
##  pair                                  (Intercept) 0.00000000

This would suggest that excluding field:pair would make a difference, and excluding hive.nested:... would make a very small difference.

summary(Snodgrasella.prev)

Suggests a weak effect (decrease of 1.84 log-odds units from Clothianidin to Control, p = 0.057). This is backed up by a slightly more accurate likelihood ratio test:

Snodgrasella.null <- update(Snodgrasella.prev,.~.-treatment)
anova(Snodgrasella.prev,Snodgrasella.null)

The weird thing is that the overall effect seems to be in the opposite direction

(stab <- with(data.f,table(Snodgrasella.detected,treatment)))
## Snodgrasella.detected Clothianidin Control
##                     0            3      11
##                     1           29      21
chisq.test(stab,simulate=TRUE)

I would again suggest collapsing the data to the most aggregated level you can (see Murtaugh "Simplicity and complexity in ecological data analysis"):

library("plyr")
d2 <- ddply(data.f,c("field","pair","treatment"),
          summarise,
          tot=sum(Snodgrasella.detected),
          n=length(Snodgrasella.detected),
          prop=tot/n)
summary(g1 <- glm(prop~treatment,weights=n,data=d2))
anova(g1,update(g1,.~1),test="Chisq")

The more reliable (LRT) p-value = 0.059, is similar to the results above ...

Some further comments, pasted here in case this gets migrated to CrossValidated:

  • if you want individual confidence intervals for each treatment separately then binom.test() will give you better results ("Clopper-Wilson" CIs);
  • however, CIs on the individual groups are answering a different question from testing the differences in proportions between groups, which is what I'm trying to address above;
  • the examples above seem to show that the results for ABPV and N.bombi are not significant at $\alpha<0.05$, but very low-powered; the results for Snodgrasella are marginally significant
  • if you're worried about the chi-squared warnings from the proportion tests, Fisher's test is OK -- some good statisticians have issues with it but in this case I think the practical differences will be small. You could check out the version of chisq_test in the coin package for another alternative ...
$\endgroup$
6
  • $\begingroup$ Thanks for the excellent answer. I was trying to find the confidence intervals of each of the proportions like this: library(plyr) `ABPV<-ddply(data.f,.(treatment),summarise, prop=sum(ABPV.detected)/length(ABPV.detected), low=prop.test(sum(ABPV.detected),length(ABPV.detected))$conf.int[1], upper=prop.test(sum(ABPV.detected),length(ABPV.detected))$conf.int[2])`` $\endgroup$
    – bee guy
    Aug 31 '15 at 11:37
  • $\begingroup$ It does not look significant to me when I plot that. Why is that? Should I be worried about the warnings that the Chi^2 approximation of prop.test may be incorrect? I read fishers test can be used for small sample sizes. What do you think about that? $\endgroup$
    – bee guy
    Aug 31 '15 at 11:52
  • $\begingroup$ (1) if you want confidence intervals for each treatment then binom.test() will give you better results; (2) that's different from testing the differences in proportions between groups; (3) I don't think I said it was significant above (the p-value from chisq.test() with simulations was 0.3603 ...); (4) Fisher's test is OK -- it has some definitional issues <andrewgelman.com/2009/05/15/i_hate_the_so-c> -- but I think all the different tests are going to tell you in this case that at least for ABPV/N. bombi, you just don't have enough information to say much ... $\endgroup$
    – Ben Bolker
    Aug 31 '15 at 15:41
  • $\begingroup$ Thanks a lot for the answers. I couldnt comment before as a new user. I have plotted now both the proportions and the differences, both with CIs. Of course you were right that you didn’t say ABPV was significant, I looked at a different organism. One more question: I used a lmer to determine the effect of the treatment on the weight of the bees. I have ten bees per colony so 640 “effective” samples. Can I use the lmer in that case? $\endgroup$
    – bee guy
    Sep 3 '15 at 7:38
  • $\begingroup$ Hello again, I still have a question about the LRT. Does it matter if the normality assumption (of anovas) is violated in g1? $\endgroup$
    – bee guy
    Feb 4 '17 at 18:07

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.