I am attempting to perform a power test to determine the minimum sample size of my next upcoming experiment. First, I made a data set:

gdata = expand.grid((rep(1:10,4)), Lineage=c("A","B","C", 
         "D", "E", "F","G", "H", "I", "J"))
gdata$Treatment <- (Treatment=c("low","high"))
gdata$Genome = expand.grid(rep (1:40,5), Genome=c("top","bottom"))

Then, I changed some items in the dataset to factors. Lineage is already a factor:

gdata$Genome <- factor(gdata$Genome$Genome)
gdata$Treatment <- factor(gdata$Treatment)

I also made some factors as objects just to make the generalized linear mixed model easier to run:

Genome <- gdata$Genome
Lineage <- gdata$Lineage

I then added gamma distributed dependent variable to the data set since I am expecting our data to be gamma distributed: If there is another package that runs general linear mixed models I am also willing to do that. I am assuming my data will be gamma distributed because past data have been, but my data may be normally distributed in the future OR I may be able to perform a transform on my future dependent variable.

gdata$G = rgamma(100, 2, 11)

Here, I want to run a generalized linear mixed model with a gamma distribution (the default link function is fine). The dependent variable is by treatment, genome, an interaction of treatment and genome (all fixed factors). There is a random effect of lineage and lineage is nested within genome. This nesting can be seen in my dataset - Lineages A to E are Top genomes and Lineage F to J are Bottom genomes. Please let me know if my model is incorrect.

mod1 <- glmer (G ~ Treatment*Genome + ((1|Lineage)%in%Genome), 
         data = gdata, family = "Gamma")

I receive the warning: boundary (singular) fit: see help('isSingular'). I assume my random effect size is very small and/or the variance among my lineages are near 0.

I set the effect sizes of my fixed factors from the GLM (based on fixed effect from past work):

fixef(mod1)["Treatmentlow"] <- 0.4476
fixef(mod1)["Genomebottom"] <- 0.4476
fixef(mod1)["Treatmentlow:Genomebottom"] <- 0.4476

I continue to perform a power analysis with the simr package:


I receive an error that I don't have a default fixed effect, which I believe I do in my glmer. So I decided to run the power analysis on my random effect:

powerSim(mod1, test = random(((1|Lineage)%in%Genome)))

But I receive the error message of of my random effect being an unused argument.

This is one of the first times I am making a dataset in R so my problem might be my dataset. There is definitely a problem with the glmer, but I thought that powerSim would still allow me to use my model. Does anyone have any advice?

More information:

If I do specify a fixed effect like: powerSim(mod1, test = fixed(gdata$Genome)) I receive the warning: boundary (singular) fit: see help('isSingular') several times over a few minutes until I stop the code.

  • $\begingroup$ There are two errors in your code. (I assume that fixing those will still not fix the issue as @Sointu explains.) 1) expand.grid(rep (1:40,5), Genome=c("top","bottom")) generates data frame with 400 rows but gdata$G = rgamma(100, 2, 11) generates 100 random variables (R recycles them 4 times). b) ((1|Lineage)%in%Genome) is not valid syntax. It's not clear what you want to achieve with it. $\endgroup$
    – dipetkov
    Commented Aug 28, 2022 at 17:22

1 Answer 1


I think you get the error message from powerSim because you haven't specified the fixed effect you want to test in the powerSim command. So you need to specify:

powerSim(mod1, test=fixed("Treatment*Genome"))

For what it's worth, I've often had issues with singular fit when generating discrete or count data and using it to experiment with different types of glmer models and my problem also seems to be that my generated variables have too little variance per cluster(s), just like you have.

Edit to add: I've found this tutorial helpful in running power analyses: https://humburg.github.io/Power-Analysis/simr_power_analysis.html

Edit 2: I can't comment so as a reply to your reply: I don't know why powerSim doesn't pick up your default fixed effect, but specifying the effect seems to remove the error. I do think you getting singular fit is due to low variance, maybe try increasing the cluster-specific variance in your data? Good luck!

  • $\begingroup$ Thank you for your response. I do not think that specifying my fixed effect is the problem here since the CRAN package says that the powerSim function “By default, the first fixed effect in fit will be tested.” If I do specify a fixed effect like: powerSim(mod1, test = fixed(gdata$Genome)) I receive the warning: boundary (singular) fit: see help('isSingular') several times over a few minutes until I stop the code. If I make a fixed factor an object: Genome <- gdata$Genome. I still receive the same isSingular warning over and over again. I will take a look at that link you posted. $\endgroup$
    – bribina
    Commented Aug 2, 2022 at 19:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.