title: How to specify the non-linear interaction of two factor variables in generalised additive models [R]
body: <p>I have a time-series data set that contains an outcome variable which is continuous and two factor predictors (one with 6 levels and one with 2 levels).</p>
  <p>I would like to model the non-linear interaction of the two factor variables on the continuous variable.</p>
  <p>This is the model I have so far:</p>
  <pre><code>library(mgcv)
  
  model &lt;- bam(
      outcome ~
          factor_1 + factor_2 +
          s(time, k = 9) +
          s(time, by = factor_1, k = 9) +
          s(time, by = factor_2, k = 9),
      data = df
  )
  
  summary(model)
  Family: gaussian 
  Link function: identity 
  
  Formula:
  outcome ~ factor_1 + factor_2 + s(time, k = 9) + s(time, by = factor_1, 
      k = 9) + s(time, by = factor_2, k = 9)
  
  Parametric coefficients:
              Estimate Std. Error t value Pr(&gt;|t|)    
  (Intercept)  2612.72      23.03 113.465   &lt;2e-16 ***
  factor_1b      33.19      27.00   1.229     0.22    
  factor_2z    -488.52      27.00 -18.093   &lt;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(time)           2.564  3.184  6.408 0.000274 ***
  s(time):factor_1b 1.000  1.001  0.295 0.587839    
  s(time):factor_2z 2.246  2.792 34.281  &lt; 2e-16 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  R-sq.(adj) =  0.679   Deviance explained = 69.1%
  fREML = 1359.6  Scale est. = 37580     n = 207
  </code></pre>
  <p>Now I would like to add a non-linear interaction of <code>factor_1</code> and <code>factor_2</code> with <code>time</code> for the effect on <code>outcome</code>, so that the smoothers in every combination could differ (for example: <code>factor_2</code> has a stronger non-linear effect for some levels of <code>factor_1</code>). Something like <code>s(time, factor_1, factor_2)</code> or <code>s(time, factor_1, by = factor_2)</code> does not work.</p>