search

gavinsimpson / gratia

ggplot-based graphics and useful functions for GAMs fitted using the mgcv package
https://gavinsimpson.github.io/gratia/
Other
206 stars 28 forks source link

problem in code reproduction #46

Closed hifrank42 closed 1 year ago

hifrank42 commented 5 years ago

When I followed the code in "Supplementary materials for: Modelling palaeoecological time series using generalized additive models" the result of

sw.cint <- confint(mod, parm = "Year", newdata = newYear, type = "confidence")
sw.sint <- confint(mod, parm = "Year", newdata = newYear, type = "simultaneous")
head(sw.sint)

are different from the guide.

  smooth  by_variable  Year   est    se  crit   lower upper
  <chr>   <fct>       <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl>
1 s(Year) NA          1846. 0.478 0.163  2.92 0.00223 0.954
2 s(Year) NA          1847. 0.481 0.152  2.92 0.0357  0.926
3 s(Year) NA          1848. 0.483 0.143  2.92 0.0665  0.900
4 s(Year) NA          1849. 0.486 0.134  2.92 0.0940  0.878
5 s(Year) NA          1850. 0.489 0.127  2.92 0.118   0.859
6 s(Year) NA          1850. 0.491 0.121  2.92 0.138   0.845

Thus, the d15N values are lower than original values, despite with the same trend. cannot find the bug after a few trying, please help.

gavinsimpson commented 5 years ago

These values are correct; the est column is the value of the smooth, which is centred about 0. To recover the values in the paper you need to add on the model constant / intercept term.

I need to revise the code in the supplement but this change was necessary to handle models with multiple smooths.

hifrank42 commented 5 years ago

Thanks for answering, in addition, when using " intervals(mod$lme, which = "var-cov")$corStruct " to get confidence interval, sometimes creating error "Non-positive definite approximate variance-covariance Consider 'which = "fixed"'", changing values of k can solve, but I don't know why.

gavinsimpson commented 5 years ago

As mentioned in the paper, a model with a wiggly trend and a CAR(1) process is often unidentifiable; both components are effectively the same thing if one or the other is not constrained in some way. The error you are seeing occurs when the fitted model is unstable, likely for this identifiability reason. More generally, this error indicates the model is likely over-fitted; you need only one of the wiggly trend or the CAR(1).

gavinsimpson commented 1 year ago

I have updated the supplementary materials to work with the latest version of gratia (on GitHub). This also "fixes" the issue with the CAR(1) model, wherein I can now reproduce the near-unidentiable fit with default values. See the supplementary materials at the paper GitHub repo: https://github.com/gavinsimpson/frontiers-palaeo-additive-modelling