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)
anova(glmm1VarB,glmm1)
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.
Questions:
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.
lmerTest
to get the p-value, as opposed to relying on LRT. $\endgroup$