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 = '?')
y is a variable (EEG power) measured at 10Hz which I expect to be roughly exponentially distributed, 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
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, :
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)
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:
descdist from the
fitdistrplus package suggests that while both
ynew have high kurtosis and skewness, this is even more extreme for
> 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?