Although I have used R to explain my question, I think and hope this is of general interest and is not off-topic here.

I want to test whether the continuous variable “varA” is different between the two levels of the categorical variable “varB”. I know that “varA” is affected also by “varC” and “varD” that are categorical variable with, respectively, 5 and 12 levels each one. Data have been collected over n years and I assume that there is not a systematic (fixed) effect of the variable “year” on “varA”, however my expectation is that “year” may add some random noise to “varA”. Therefore, I decide to use a GLMM approach (lmer function from package lme4 in R). My syntaxis in R is this:

glmm1<- lmer(varA~varB+varC+varD+(1|Year),data=dati,REML = F)

where “dati” is the dataframe containing all the data.

Then I obtain the p value by using a Likelihood Ratio Test:

glmm1VarB<- lmer(varA~ varC+varD+(1|Year),data=dati,REML = F)

I get a highly significant effect of varB on varA and, by using the fitted values of the model, I calculate that the decrease in varA when varB= “1”, is about 30% with respect to the value of varA when varB= “0”.

However, my data are highly unbalanced at such an extent that I have about 300 observations for varB= “0” and just 20 for varB= “1”. Then, I decide to see what I would get if I sampled, by random, the same number of observations (20) for varB= “0” and repeated the entire (above) analysis 1000 times. In R:

dati0<- dati[dati$varB==”0”,]
dati1<- dati[dati$varB==”1”,]
p<- diffPerc<- vector()
for (i in 1:1000){
  datiboot<- rbind(dati0[sample(nrow(dati0),nrow(dati1)),],dati1)
  glmm1boot<- lmer(varA~ varB+varC+varD+(1|Year),data=dati,REML = F)
  glmm1bootVarB<- lmer(varA~ varC+varD+(1|Year),data=dati,REML = F)
  p[i]<- (anova(glmm1boot,glmm1boot))[2,8]
  diffPerc[i]<-((effect("varB",glmm1boot)$fit[1] effect("varB",glmm1boot)$fit[2])/effect("varB",glmm1boot)$fit[1])*100

Then, by having a look at the 0.05, 0.5 and 0.95 quantiles of diffPerc I found that varA decreases between varB= “0” and varB= “1” of about 23% (between 15% and 30%). Aso, I found that 8% of times I did not get a significant (at alpha= 0.05) result.


1) I know that glmm is robust to unbalanced data, however, what is the effect of unbalanced data on the statistical power of a glmm like this? I have googled quite a bit on this question but I have not made a clear idea on that. In general I have found that it should lower your statistical power but sometimes it is mentioned that increases also the type I error.

2) if you were in my shoes, which analysis would you be more confident with (the first approach or the bootstrap)?

3) would it be more appropriate to assess the significance of the bootstrap analysis by using the Fisher's method that combines the p-values by the sum of logs? in that case (using the function "sumlog" from package "metap") I obtain that the p is still very significant in a way similar to the original analysis without bootstrapping.

  • 1
    $\begingroup$ (1) You don't need to care about power after you already found a significant result. (2) I definitely prefer the first analysis where you use all available data at once. (3) This wouldn't make sense because the result is going to depend on the number of bootstrap iterations with combined p->0 when n_iter->infinity which does not make sense. $\endgroup$ – amoeba says Reinstate Monica Jan 18 '18 at 12:25
  • 1
    $\begingroup$ PS. You might consider modeling varC and varD as random, perhaps even including random slopes as well: (1|varC) or (varB|varC) and similarly for varD. You might also consider using lmerTest to get the p-value, as opposed to relying on LRT. $\endgroup$ – amoeba says Reinstate Monica Jan 18 '18 at 12:27
  • $\begingroup$ It is not clear how you want to define variable B? two levels or it could be eithe 0 or 1. Also, it seems that you have several issues in the single question. Try to be brief and specific. $\endgroup$ – Subhash C. Davar Jan 18 '18 at 13:40
  • 1
    $\begingroup$ Thank you. I did not realize that using the Fisher's method does not make any sense here because the number of p-values is completely "artificial" and it strives with the logic behind the method. To see it I have tried this: > sumlog(rnorm(10,0.3,0.1)) chisq = 26.66289 with df = 20 p = 0.1450174 > sumlog(rnorm(100,0.3,0.1)) chisq = 257.7205 with df = 200 p = 0.003680841 $\endgroup$ – user3844454 Jan 18 '18 at 13:42
  • 1
    $\begingroup$ @amoeba has really answered all your questions. Just for the record the behaviour of some methods for combining $p$-values as $N\rightarrow\infty$ has received some discussion in this Q&A stats.stackexchange.com/questions/243003/… $\endgroup$ – mdewey Jan 18 '18 at 13:59

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.