My apologies, this is a long question: the TL;DR version is "my GAMs don't converge, is this a distribution issue and if so can I fix it?"

Edit: fixed the distribution information. ynew is imported in 10*log10 (dB) format, and so to put it on the same scale as the original y, the reverse operation has to be performed (ynew <- 10^(ynew/10)). I did check the new data (visually) and they looked fine, but it turns out that in the process of merging the old with the new data I accidentally had the reverse operation performed again, i.e., double exponentiation.

Thanks a lot for pointing out that the range was off - I was so mystified by the fact that it didn't work that for three days I couldn't see what was right in front of me...

For a research project I am trying to model responses to three types of stimuli using GAMs. I expect these responses to be long and therefore to overlap (in time) to some degree. I am therefore fitting models of the following general form using mgcv::bam in R:

mod <- bam(y <- s(x1) + s(x2) + s(x3), family = ...(link = '?')

Where y is a variable (EEG power) measured at 10Hz which I expect to be roughly exponentially distributed, and x1, x2, and x3 are (constructed) predictors which are set to a small negative value for many of the samples. The data look something like this:

    # A tibble: 6,184 x 1
   forcv$time      $x1      $x2     $x3    $y
        <dbl>    <dbl>    <dbl>   <dbl> <dbl>
 1        0   -0.00588 -0.00588 -0.0333  3.37
 2        0.1 -0.00588 -0.00588 -0.0333  2.12
 3        0.2 -0.00588 -0.00588 -0.0333  1.47
 4        0.3 -0.00588 -0.00588 -0.0333 13.3 
 5        0.4 -0.00588 -0.00588 -0.0333 10.0 
 6        0.5 -0.00588 -0.00588 -0.0333  3.38
 7        0.6 -0.00588 -0.00588 -0.0333 10.3 
 8        0.7 -0.00588 -0.00588 -0.0333 16.8 
 9        0.8 -0.00588 -0.00588 -0.0333  3.26
10        0.9 -0.00588 -0.00588 -0.0333  1.33
# ... with 6,174 more rows

At intervals, the predictors will increase linearly from 0 to a maximum of 1, like this: plot of x1 vs time

These non-zero intervals overlap between the three predictor variables, to account for the expected overlap between the responses.

I have used this approach with one method of measuring the response, using family = Gamma(link = 'log'), or family = gaussian(link = 'log'). This yields smooths for the predictor variables that look extremely plausible (although the model is clearly far from complete, e.g., observed vs fitted plots show virtually no correlation).

Now, I have implemented a new method of calculating y, i.e., using FFT rather than the Hilbert transform, and suddenly my models do not converge anymore. I get errors like

Error in bgam.fit(G, mf, chunk.size, gp, scale, gamma, method = method, : non-finite deviance

sometimes with the additional warning

In addition: Warning message: In log(ifelse(y == 0, 1, y/mu)) : NaNs produced


Possible divergence detected in fast.REML

In some but not all cases I can get the model to run by using

mod <- bam(ynew <- s(x1) + s(x2) + s(x3), family = ...(link = '?'), discrete = TRUE)

although the results often look completely implausible. By playing around with the data I've discovered that I can also get around these warnings by switching to:

mod <- bam(log(ynew) <- s(x1) + s(x2) + s(x3), family = gaussian(link = 'identity')

which produces roughly plausible results.

Density plots of the distributions of (log) y and ynew do not make it immediately apparent to me that there are large differences between the two distributions, although ynew does seem to have considerably heavier tails:

density of log(y) and log(ynew)

Using descdist from the fitdistrplus package suggests that while both y and ynew have high kurtosis and skewness, this is even more extreme for ynew:

> descdist(y)
summary statistics
min:  0.01059139   max:  35.1851 
median:  0.3602276 
mean:  0.6296528 
estimated sd:  1.00021 
estimated skewness:  12.29973 
estimated kurtosis:  332.6891

> descdist(ynew)
summary statistics
min:  0.0001107598   max:  29.6122 
median:  0.5384826 
mean:  1 
estimated sd:  1.478953 
estimated skewness:  5.454704 
estimated kurtosis:  58.4672 

In summary, I have three (related) questions:

  • Why do these problems occur?
  • Is modeling the log-transformed data appropriate, or does that just mask problems with the data?
  • Can I fix this by some form of scaling?

put on hold as off-topic by mkt, mdewey, Michael Chernick, BruceET, jbowman 2 days ago

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – mkt, Michael Chernick, jbowman
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 2
    $\begingroup$ Why does the new method yield such high values, like the maximum value of 7.251297e+12? Are these high values even plausible or are they a sign that something is off with the way they are computed? If the new method is supposed to measure the same thing as the old one, how come the new method produces these humangous values? How many values of ynew exceed the highest value of 35.1851 for y? Before fitting any model to ynew data, you should do this type of common sense check to determine what is going on. $\endgroup$ – Isabella Ghement Aug 13 at 15:02
  • 1
    $\begingroup$ Also, from a practical point of view, having data that varies over 12 orders of magnitude to such large numbers is asking for trouble. Perhaps rescale, by dividing by a million or some such number. $\endgroup$ – Gavin Simpson Aug 13 at 17:11