We were unable to load Disqus. If you are a moderator please see our troubleshooting guide.

Join the discussion…

  • in this conversation
      Media preview placeholder
      Log in with
      or sign up with Disqus or pick a name

      Disqus is a discussion network

      • Disqus never moderates or censors. The rules on this community are its own.
      • Your email is safe with us. It's only used for moderation and optional notifications.
      • Don't be a jerk or do anything illegal. Everything is easier that way.

      Read full terms and conditions

      By signing up, you agree to the Disqus Basic Rules, Terms of Service, and Privacy Policy.
      By posting, you agree to the Disqus Basic Rules, Terms of Service, and Privacy Policy.
      • Avatar

        Interesting. Good to see paleos using R! Thanks for the post. I'm going to need to find the time to play around with this a bit...

        • Avatar

          Interesting post and nice usage of GAMM. What was your basis for choosing the smoothing parameter, k? And why does it change from 20 (m1) to 30 (m2, m3)?

            • Thanks David. That is a typo, but it doesn't affect the results IIRC. If you use the default k, then the smoothness selected in m1 is close to the maximum allowed by k. In such circumstances, the advice is to fit the model allowing for a larger basis dimension which the model will penalise back to an "optimal" value for the smoothness. I realised whilst writing the post that setting k = 20 made no difference to the fit obtained using k = 30 so went with the smaller starting value, but forgot to update the code for m2 and m3. In general I use the defaults in gam() / gamm() unless I expect more variation in the trend/level than the defaults allow. To some extent this is dependent upon the data you have available, also.

              • This post is absolutely outstanding. Thanks for teaching me everything I needed to know about additive models in just one post. I used the R code extensively to fit an identical model over health expenditure data, with interesting results.

                  • Avatar

                    Thanks for the informative post! Couple of observations:

                    1) The downloads from github didn't work for me on Windows at first, but removing the option (...,method = "wget") made it work just fine.

                    2) Didn't know about "month.abb". Very handy!

                    I'm getting more acquainted with GAMMs, and your example was really helpful. Cheers!

                      • Avatar

                        Hi Gavin
                        I have downloaded your code but cannot get the Deriv function to work. I copy and paste all your code from the main post but get an error at this step:
                        > m2.d <- Deriv(m2, n = 200)
                        Error in terms.default(mod) : no terms component nor attribute

                        Do you know what the problem might be? I am using R version 2.14.0 for Windows.

                        Another question, you write:
                        "note that I only intend to use the annual means in this posting — dealing with monthly data needs a few extra steps to model the seasonal variation in the data"

                        As I understand it calculating monthly anomalies removes the avarage seasonal variation in the data. However if the seasonal cycle changes over time there might be what Tamino* calls a "residual seasonal cycle". Is this what you mean by "seasonal variation in the data"?

                        *) See this:


                          • The code works for me on my laptop with R 2.14.0 and mgcv 1.7-6. That is an old version of mgcv and I needed to install that over the one supplied with R 2.14.0 as there was a bug in te() smooths I was using for another project. In looking at the ChangeLog for mgcv I note that there were issues with terms etc in returned objects, so it may be that you need to update mgcv via update.packages(ask = FALSE) and retry the code.

                            As for what I meant be seasonal variation, the post only looks at the annual data so there isn't an issue. If one were to model the monthly data series, then the model would need a seasonal smoother at least. What Tanimo is referring to is that if you just fit a seasonal smoother or some other unchanging seasonal term in the model, if the seasonal pattern of variation in the model changes with the trend/through time, your model won't capture all the cyclic variation due to season and the residuals will more than likely contain seasonal temporal autocorrelation. I would prefer to model the variation in the data rather than remove signals from the data. How you might deal with this will depend on the data; you could fit a trends for each month for quarter, instead of a single trend (see the by argument to s()), but that won't necessarily model a change in the seasonal pattern of temperatures throughout the year. An alternative is to fit a multivariate smoother for trend and seasonal components, which is where the te() smooths I mentioned above come in. I'm swamped with work at the moment, but I hope to get a chance to write a blog post on that methodology over the Christmas vacation as I'm currently writing a methods paper using these ideas.

                              • Avatar

                                Updating the mgcv package to the newest version (ver. 1.7-12) did not solve the problem, However I found a simple solution. I commented this line out in the derivFun.R file:
                                if(isTRUE(all.equal(class(mod), "list")))
                                then the next line (mod class(m2)
                                [1] "gamm" "list"

                                  • Avatar

                                    Sorry, some words dissappeared in my comment, Is should read:
                                    Updating the mgcv package to the newest version (ver. 1.7-12) did not solve the problem. However I found a simple solution. I commented this line out in the derivFun.R file:
                                    if(isTRUE(all.equal(class(mod), “list”)))
                                    then the next line is always executed.
                                    The output from class(m2) is:
                                    [1] “gamm” “list”

                              • May I ask another question: Has this approach been covered somewhere in the literature? I wonder where I can get a full-fledged explanation on the functions used, in order to apply it to health expenditure data. I have started looking at the GLAMM books cited in the R documentation, but any other pointer would be of great help. Thanks :)

                                  • Avatar

                                    Hi Gavin
                                    Thanks for your reply at Realclimate - I posted there as SRJ mentioning this post. I have also mentioned a while ago at Skeptical Science.
                                    I would appreciate it if you would show how to fit monthly data using a cyclic smoother as you mention. I tried extending your model from this post to monthly data but couldn't get it to work.
                                    By the way, updating this analysis with the data for 2011 makes the fitted model significant only to 2003, since 2011 was cooler than 2010.