powerSim will not run and repeatedly gives me the warning: boundary (singular) fit: see help('isSingular')

Question

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"))
View(gdata)
glimpse(gdata)

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)
glimpse(gdata)

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")
summary(mod1)

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:

powerSim(mod1)

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.

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!