I am struggling with understanding how to run a simulated power analysis for GAMMs. Though my planned studies will include more variables, the core issue I'm trying to work out is to simulate the following:

  • A binary outcome y. However, since I don't know how to correlate y here, I have opted to just use a Gaussian response as it is easier to model.
  • Two predictors, modeled with s(x, by = group) to estimate this for two groups.
  • The variables should be moderately correlated, with x2 having a negative association with y.
  • An interaction will be modeled as ti(x1,x2,by=group).
  • Subject and item random effects (with no specific prior information).
  • For Group 1, x1 and x2 should have a larger effect on y.


The best I could come up with was something like this for a basic GAM:

#### Load Libraries #### 

#### Perform Simulation PA ####
B <- 10 # number of simulations (10 for now so its faster)
p1 <- NULL # store p values
p2 <- NULL
p3 <- NULL 
p4 <- NULL
p5 <- NULL 
p6 <- NULL

for (i in 1:B) {
  #### Create RE Grid ####
  sub <- factor(1:50)
  item <- factor(1:20)
  obs <- length(sub)*length(item)
  grid <- expand.grid(sub=sub,
  #### Set Up Simulation Parameters ####
  n <- 50 # sample size
  alpha <- 0.05 # significance level
  power <- 0.8 # desired power
  #### Simulate Data ####
  dat <- rnorm_multi(n = n, 
                     mu = c(20, 20, 20),
                     sd = c(10, 10, 10),
                     r = c(0.30, 0, -.30), 
                     varnames = c("x1", "x2", "y"),
                     empirical = FALSE) 
  #### Make Grouping Factor ####
  group <- factor(
    labels = c("Group 1", "Group 2")
  ) # not sure how to show differences in groups
  #### Merge Data ####
  df <- data.frame(dat,group,grid) %>% 
  #### Build Model ####
  model_sim <- bam(y 
                   ~ s(x1, by = group)
                   + s(x2, by = group) 
                   + ti(x1,x2, by = group)
                   + s(sub, bs = "re")
                   + s(item, bs = "re"),
                   data = df) # fit GAM model to simulated data
  #### Store P Values ####
  p1[i] <- as.vector(summary(model_sim)$s.table[1,"p-value"])
  p2[i] <- as.vector(summary(model_sim)$s.table[2,"p-value"])
  p3[i] <- as.vector(summary(model_sim)$s.table[3,"p-value"])
  p4[i] <- as.vector(summary(model_sim)$s.table[4,"p-value"])
  p5[i] <- as.vector(summary(model_sim)$s.table[5,"p-value"])
  p6[i] <- as.vector(summary(model_sim)$s.table[6,"p-value"])
  p.data <- data.frame(p1,p2,p3,p4,p5,p6)

#### Inspect P Values for Variables 1:6 ####


Running this as-is takes forever, even with the very small number of boots here $(B = 10)$ and often doesn't converge, with this error popping up:

Warning message:
In bgam.fitd(G, mf, gp, scale, nobs.extra = 0, rho = rho, coef = coef,  :
  algorithm did not converge

I'm assuming it has to do with the data I have simulated. I also have no way of getting the groups to do anything other than arbitrarily split, while tinkering sometimes leads to complete separation in the GAMM (which I obviously don't want). As mentioned in the Question section, one group should have a predictably larger effect on y for each x, so modeling this would be helpful, Finally, achieving all of this somehow with binary data seems hard, so I have no way of figuring that part out.

So to summarize, I need to solve the following:

  • What is causing the slow/weird estimaton?
  • How do I get the groups to have a meaningful distinction from each other?
  • Is there a way to make the outcome data binary? This is not my greatest priority, but it would be a nice bonus.


(continuation of first edit)

Because I don't want to make this question verbose, I've only included my recent changes. To follow the advise in the comments (hopefully I've done that), I tried estimating $y$ in the way prescribed in the data simulation functions of mgcv and gratia. I have kept it simple and made one "positive" function that varies slightly by group and one "negative" function that varies slightly by group.

#### Load Libraries #### 

#### N and Scale ####
n <- 100
scale <- 1

#### Create X and F ####
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
f1.x1 <- exp(3 * x1) 
f2.x1 <- exp(2 * x1)
f1.x2 <- exp(3 * x2) * -1
f2.x2 <- exp(2 * x2) * -1 
e <- rnorm(n, 0, scale)

#### Create Factors ####
fac <- as.factor(sample(1:2,n,replace=TRUE))
fac.1 <- as.numeric(fac==1)
fac.2 <- as.numeric(fac==2)

#### Estimate Y ####
y <- x1 + x2 + f1.x1*fac.1 + f2.x1*fac.2 + f1.x2*fac.1 + f2.x2*fac.2 + e 

#### Merge Data ####
df <- data.frame(y=y, # y values
                 x1=x1, # x values
                 fac=fac) # factors

#### Convert and Inspect ####
data <- df %>% 

#### Fit Model ####
fit <- gam(y
           ~ fac 
           + s(x1, by = fac)
           + s(x2, by = fac)
           + ti(x1, x2, by = fac),
           method = "REML",
           data = data)

Which gives me this:

enter image description here

This gets me close to what I want for now with the main effects and their grouping (though the error seems to always be higher for the second factor). The major parts I am still missing:

  1. How do I code the interaction between the two main effects? I am aware of how you do this in a vanilla regression, but not GAMs. I've tried assembling functions to do this like ti.function <- function(x1,x2){x1*x2} into y, but this doesn't seem to do anything different. I also tried correlating the variables beforehand using this method, but this in turn changes the direction of the plotted function. To summarize, I need increases in $x_1$ to be associated with increases in $y$, increases in $x_2$ to be associated with decreases in $y$, and and a linear interaction between $x_1$ and $x_2$ When $x_1$ is high and $x_2$ is low, this should result in higher $y$.
  2. I still don't know how to add in subject/item random variance in a meaningful way. The faux package seems promising, but so far my attempts to add REs have failed.
  3. I'm having issues adjusting the error variance with the scale object in this simulation. The standard error seems to always be relatively low but I expect my actual data to have more, so I would like to have better ways to tweak this.

Another Edit...(Dec 2023)

TLDR: I have solved the perfect prediction issues and have included crossed random effects like I originally wanted, but the relationships between the predictors and outcome are still not ideal.

I have still not solved my problem after repeated attempts, so I am setting up another bounty. I realize part of my issue previously was not providing enough information about the variables. Unfortunately, there is a serious lack of a priori information for my case, so I can only fill in what I know from some work in my area. However, I can offer a lot more information based off some previous knowledge of the research and what previous answerers have suggested adding. The variables should have the following features:

  • A binary dependent variable called RD for "reading success" with values $0$ for failure to read and $1$ for success.
  • A group variable called LG for "language group", where the L1 group has a higher rate of reading success than the L2 group. Its unclear what this means from a logistic sense. In one analysis the L1 group scored a mean of $120/150$ items correct and the L2 group merely got $30/150$, but the L2 group came from a poor disadvantaged background. Another study only found a $15$ point difference in scores on average (with the L2 group again scoring lower), but the L2 group in this case was from a much higher SES group. The actual reading test is a scale which ranges from a total score of $[0,150]$.
  • A discrete subject-level variable called OA which has a positive linear effect on RD for the L2 group and a more weak curvilinear effect for the L1 group. There does not appear to be strong differences between the groups in at least one study (L1 OA = $12/16$, L2 OA = $11/16$) but I imagine those differences also vary wildly in reality and I have a strong feeling that the differences are more extreme than what is reported in just one study. As one can see, the scores range between $[0, 16]$. However, another study uses a similar measure range between $[0,20]$ with reported scores of $11/20$ for the L2 group but no reported scores of L1 groups with this measure.
  • A discrete item-level variable called GC which has a negative linear effect on reading for the L2 group and a more weak curvilinear effect for the L1 group. The mean of the GC variable should be about $115$ and has no ceiling value, so it spans the range $[0, \infty)$. A previous analysis showed that GC was linear for reading and curvilinear for writing for the L1 group, so to be careful I am just going to assume the relationship is generally curvilinear for the L1 group, but that is of course just an assumption.
  • A subject x item crossed random effects design, with a starting estimate of $50$ subjects per group and $150$ reading items. However, because this is a power analysis, the items will be essentially fixed and I am only after the number of subjects required for adequate power.
  • Something else to add is that the estimates are likely to be fairly noisy. In the previous simulations in my question and the generous answer by Gavin, the standard error is very narrow, but I anticipate that this will not be the case given how measures in my field are not strongly correlated with the response.

The only thing I have been able to do to modify the code so far is something like this, which at least includes the random effects and logistic regression I need (thus solving the perfect prediction issues and added in random effects like I needed before), but the simulation still doesn't model the data like I would hope. Note for simplicity I just included continuous versions of the predictors because I ran into some more issues with Poisson-distributed variables:

#### Load Libraries ####

#### Set Seed ####

#### Define REs ####
n_subjects <- 100
n_items <- 150

#### Create Subject Data ####
subject_data <- data.frame(
  subject = rep(1:n_subjects, each = n_items),
  LG = rep(c("L1", "L2"), each = n_items * 50),
  OA = rnorm(n_subjects * n_items, mean = 11, sd = 1)

#### Create Item Data ####
item_data <- data.frame(
  item = rep(1:n_items, times = n_subjects),
  GC = rnorm(n_subjects * n_items, mean = 100, sd = 10)

#### Combine Data ####
data <- cbind(subject_data, item_data)

#### Define Y Variable ####
data$RD <- with(data, rbinom(n_subjects * n_items, 1, plogis(
  (LG == "L1") * (0.05 * OA - 0.0001 * GC^2) +
    (LG == "L2") * (0.03 * OA - 0.0001 * GC) +
    (LG == "L2") * (0.0002 * OA * GC)

data.fix <- data %>% 
  mutate(subject = factor(subject),
         item = factor(item))

#### Fit GAMM ####
fit <- bam(
  ~ s(OA, by = LG)
  + s(GC, by = LG)
  + ti(OA, GC, by =LG)
  + s(subject, bs = "re")
  + s(item, bs = "re"),
  family = binomial

#### Plot Effects ####

enter image description here

I'm not looking for an exact answer to my problem, so the solution doesn't need to fit all of these criteria exactly. I just need a workable logistic GAMM with crossed random effects, and if it achieves most of what I am trying to get here then I will be happy.

  • 2
    $\begingroup$ You're model is missing the group means for a start; you must include a parametric effect group if using a factor by term. If you're doing a power analysis, you probably should use te() in place of the three separate smooths as the te() form is a simpler model $\endgroup$ Commented Jun 16, 2023 at 7:28
  • 1
    $\begingroup$ I really don't follow what you are trying to do; I appreciate that faux makes some of this easy, but it is obfuscating what you are actually (trying to) doing or want to do. The mgcv manual and the tests in gratia (see the setup.R file) contain lots of examples of reference models and simulating from them. You need to specify shapes for the functions, plus write out the linear predictor Lp, and then simulate data using a suitable RNG function. If you can spell out those things, I can see about producing an answer $\endgroup$ Commented Jun 17, 2023 at 12:11
  • 1
    $\begingroup$ I have attempted to do what you say with the formulation shown in gratia and mgcv. That seems to give me the main effects and grouping effects I want, but I'm still struggling with coding in the interactions and random effects. I tried multiplying the x1 and x2 variables in the y object but that didn't seem to fix anything. See my edit for more deets. $\endgroup$ Commented Jun 17, 2023 at 16:44
  • 2
    $\begingroup$ That's a bit rough @FrankHarrell. The OP isn't approximating a binary with a Gaussian. They don't know how to simulate binary data so they're figuring out the first step (simulating data from a GAM with known functions) for Gaussian Y (to which they fit a Gaussian GAM), and simulating a binary $Y$ is an extra ask on top of how to simulate a smooth interaction. In a learning situation, it can often be cognitively simpler to use inefficient code that more closely follows simple math steps than highly optimized code. Once they understand how to do the first bit, they can optimize. $\endgroup$ Commented Jun 18, 2023 at 11:40
  • 2
    $\begingroup$ @FrankHarrell I apologize if my code looks poor. I'm self-taught in R and try to first understand every piece of my code before trying to make it more concise. I'm also not amazing at math (believe me I try to learn as much as I can on the side while I'm in grad school), so if I look clueless it's because this is one of the few places I can get feedback to get better. Feel free to provide a solution if you feel you can use this as a teaching moment. I would appreciate whatever wisdom you can share. $\endgroup$ Commented Jun 18, 2023 at 11:53

1 Answer 1


This is a work in progress...

Using your example:

#### Load packages #### 

#### N and Scale ####
n <- 100
scale <- 2

#### Create X and F ####
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
f1 <- function(x, b, c = 1) exp(b * x) * c
f2 <- function(x, b, c = 1) exp(b * x) * c
f1_x1 <- f1(x1, b = 3) 
f2_x1 <- f2(x1, b = 2)
f1_x2 <- f1(x2, b = 3, c = -1)
f2_x2 <- f2(x2, b = 2, c = -1)
e <- rnorm(n, 0, scale)

#### Create Factors ####
fac <- as.factor(sample(1:2,n,replace=TRUE))
fac1 <- as.numeric(fac==1)
fac2 <- as.numeric(fac==2)

#### Estimate Y ####
y <- x1 + x2 + f1_x1*fac1 + f2_x1*fac2 + f1_x2*fac1 + f2_x2*fac2 + 
     ((f1_x1*f1_x2) * fac1) + ((f2_x1*f2_x2) * fac2) + # interactions

#### Merge Data ####
df <- tibble(y = y, x1 = x1, x2 = x2, fac = fac)

Absent any information on the nature of the interactions, I simply created the interactions $f_{\text{fac}}(x_1,x_2)$ in the same way you'd create an interaction in a linear model, you multiply the values of the two predictors together, hence:

((f1_x1*f1_x2) * fac1) + ((f2_x1*f2_x2) * fac2)
#### Fit Model ####
fit <- gam(y ~ fac + s(x1, by = fac) + s(x2, by = fac) + ti(x1, x2, by = fac),
           method = "REML",
           data = df)

Now we get a significant pair of pure interaction tensor products:

> summary(fit)                                                                

Family: gaussian 
Link function: identity 

y ~ fac + s(x1, by = fac) + s(x2, by = fac) + ti(x1, x2, by = fac)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -39.3796     0.5046  -78.05   <2e-16 ***
fac2         30.0963     0.5935   50.71   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                  edf Ref.df       F p-value    
s(x1):fac1      5.406  6.299  939.12  <2e-16 ***
s(x1):fac2      2.231  2.761   33.18  <2e-16 ***
s(x2):fac1      5.973  6.805 1260.83  <2e-16 ***
s(x2):fac2      3.451  4.194  138.91  <2e-16 ***
ti(x1,x2):fac1 11.600 13.004  358.29  <2e-16 ***
ti(x1,x2):fac2  4.154  5.728   23.10  <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.997   Deviance explained = 99.8%
-REML = 245.19  Scale est. = 4.1068    n = 100

And we can inspect the fitted smooths with


It can be hard to see if the estimated functions do what you want from that decomposed fit, so fitting with te() should help:

#### Fit Model ####
m <- gam(y ~ fac + te(x1, x2, by = fac),
         method = "REML", data = df)

enter image description here

But you can also use vis.gam() from {mgcv} to visualise the estimated surfaces:

layout(matrix(c(1,2), ncol = 2))
vis.gam(m, view = c("x1", "x2"), theta = 180, phi = 20,
        cond = list(fac = 1), main = "fac == 1")
vis.gam(m, view = c("x1", "x2"), theta = 180, phi = 20,
        cond = list(fac = 2), main = "fac == 2")

which produces:

enter image description here

So the estimated functions fit the requirements, I believe.

To check all of this, you want to use regular data over the range of the covariates so that you can actually plot the true surfaces you are trying to estimate:

g1 <- g2 <- seq(0, 1, by = 0.01)
gg <- expand.grid(g1, g2)
f1f2 <- function(x1, x2, b = c(3,2), c = c(1,-1)) {


For example, f1 for x1:

g1 <- seq(0, 1, by = 0.01)
ggplot(data = data.frame(x = g1), aes(x)) + 
    geom_function(fun = f1, args = list(b = 3, c = 1))

And this is the linear interaction of the two non-linear functions:

ff <- function(x, b) {
  f1_x1 <- f1(x[,1], b = b)
  f1_x2 <- f1(x[,2], b = b, c = -1)
  f1_x1 * f1_x2

gg <- expand.grid(x1 = g1, x2 = g1)
dt <- cbind(gg, y = ff(gg, b = 3))

ggplot(dt, aes(x = x1, y = x2)) + 
    geom_raster(aes(fill = y)) + 
    geom_contour(aes(z = y), colour = "black")

enter image description here

Or, for the full effect specified in the linear predictor

ff_full <- function(x, b) {
  f1_x1 <- f1(x[,1], b = b)
  f1_x2 <- f1(x[,2], b = b, c = -1)
  x[,1] + x[,2] + f1_x1 + f1_x1 + (f1_x1 * f1_x2)

gg <- expand.grid(x1 = g1, x2 = g1)
dt <- cbind(gg, y = ff_full(gg, b = 3))

ggplot(dt, aes(x = x1, y = x2)) + 
    geom_raster(aes(fill = y)) + 
    geom_contour(aes(z = y), colour = "black")

enter image description here

Being able to plot the truth (as I've done at the end (change the value of b in the creation of gg to get the second group)) is important and you should really start there so that you can get the shapes of the functions you anticipate correct. Then it is pretty simple to then simulate data from those functions and fit a model to them.

Simulating binary data

This becomes easier when you consider the form of a generalized linear model. There we would have

$$y_i \sim \text{Bernoulli}(\mu_i)$$


$$\text{logit}(\mu_i) = \eta_i$$

where $\mu_i$ is the expected value on the response scale (a probability that $y_i = 1$, give the value of $\eta_i$, the linear predictor, and $\text{logit}()$ is the canonical link function for the binomial GLM (of which the Bernoulli distribution is a special case).

There's no error term $\varepsilon$, your e in a GLM; the data are a random draw from the specified distribution with expectations $\mu_i$.

Assuming your linear predictor is as you wrote for the Gaussian model, then we have

$$\eta_i = x_{1i} + x_{2i} + f_1(x_{1i}) + f_2(x_{2i}) +( f_1(x_{1i}) \times f_2(x_{2i}))$$

Now, we want $\mu_i$, but $\eta_i$ gives us $\text{logit}(\mu_i)$. We need the inverse of the link function such that we have:

$$\mu_i = \text{logit}^{-1}(\eta_i)$$

where $\text{logit}^{-1}$ is the inverse of the logit link function. Then we just simulate from rbinom() with p being the value of $\mu_i$ we want and size = 1).

There's some rescaling perhaps needed on the linear predictor to get it with a suitable range, but it's not much more difficult than that.

In R code, we might encapsulate the simulations in a function:

`sim_binary` <- function(x, scale = 2, ...) {
    ilink <- inv_link(binomial())
    x <- (x - 5) * scale
    p <- ilink(x)
    tibble(y = rbinom(n = p, size = 1, prob = p), f = x)

The rescaling I mentioned is that subtraction of 5 from the values of the linear predictor, which are the x in that function. And notice that p is the value of the linear predictor on the response scale (after application of the inverse of the link function), but because of how rbinom() and co work, we can specify p as the argument n, so that we get the right number of simulated values, each with their own probability that $y = 1$ as given by the model (the inverse of the linear predictor).

Putting this together, reusing things from earlier, we have

# the known linear predictor
lp <- x1 + x2 + f1_x1*fac1 + f2_x1*fac2 + f1_x2*fac1 + f2_x2*fac2 + 
     ((f1_x1*f1_x2) * fac1) + ((f2_x1*f2_x2) * fac2)

# simulate binary data
yb <- sim_binary(lp, scale = scale)

which doesn't work:

> head(yb)                                                                    
# A tibble: 6 × 2
      y      f
  <int>  <dbl>
1     0  -44.3
2     0  -49.7
3     0  -39.6
4     0 -509. 
5     0 -125. 
6     0  -67.9
> any(y == 1)                                                                 

In the output f is the value of $\eta$ on the link scale. These are tiny values, most of the action in the logit (== logistic) function happens between -4 and 4:

enter image description here

(Image from Wikipedia: https://en.wikipedia.org/wiki/Logistic_function)

So we need to recentre the values of our linear predictor around say 0. The -5 in gratia::sim_binary() is doing this for the function you can simulate from in gratia.

`sim_binary2` <- function(x, scale = 5, shift = 5, ...) {
    ilink <- inv_link(binomial())
    x <- (x - shift) * scale
    p <- ilink(x)
    tibble(y = rbinom(n = p, size = 1, prob = p), f = x)

# the known linear predictor
lp <- x1 + x2 + f1_x1*fac1 + f2_x1*fac2 + f1_x2*fac1 + f2_x2*fac2 + 
     ((f1_x1*f1_x2) * fac1) + ((f2_x1*f2_x2) * fac2)

# simulate binary data
yb <- sim_binary2(lp, scale = scale, shift = -15)

After a bit of trial and error I used a value of -15 for the shift.

Now we can fit our model:

m_bin <- gam(y ~ fte(x1, x2, by = fac),
         method = "REML", data = yb, family = binomial())

but we now also get a warning:

> m_bin <- gam(y ~ te(x1, x2, by = fac), 
             method = "REML", data = yb, family = binomial())                   
Warning message:
In newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,  :
  Fitting terminated with step failure - check results carefully
> summary(m_bin)                                                              

Family: binomial 
Link function: logit 

y ~ fac + te(x1, x2, by = fac)

Parametric coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    -151.6   623644.4       0        1
fac2            982.6  5116057.1       0        1

Approximate significance of smooth terms:
                 edf Ref.df Chi.sq p-value
te(x1,x2):fac1 3.981   4.05      0       1
te(x1,x2):fac2 3.000   3.00      0       1

R-sq.(adj) =      1   Deviance explained =  100%
-REML = -123.35  Scale est. = 1         n = 100

wherein the problem is immediately clear: we can perfectly predict the data, the massive standard errors on the parametric terms are an indicator of (quasi) complete separation. I haven't had time yet to look into why, but I suspect it is the values of the functions you are using the two groups. I think they are just too different on the link scale so it is trivial to guess the 0s and 1s once you know the fac.

Something to think about; you don't actually include any group means in your linear predictor (simulated data). The only differences between the groups comes in the form of the values 3 and 2 used in the exponential functions.

  • 2
    $\begingroup$ Thank you for your input. Since you put in a lot of effort, I have awarded the bounty. I will try to work on the other unanswered parts of this question and hopefully update it when I have a solution to some of the separation issues you noted. $\endgroup$ Commented Jun 24, 2023 at 3:20
  • 2
    $\begingroup$ @ShawnHemelstrand Sorry, ran out of time to return to the random effect bit before vacation beckons. I'll circle back, but you can look again at that setup.R file in gratia's tests/testthat folder on GitHub as it has an example for random effects based on Simon's example. The general idea is break the n data into g groups and provide say a group mean plus random noise to those observations in the gth group. $\endgroup$ Commented Jun 24, 2023 at 7:27
  • 2
    $\begingroup$ No worries. Enjoy your vacation. I'll check that out in the meantime. $\endgroup$ Commented Jun 24, 2023 at 7:48

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.