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
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.
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"?
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.
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.
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)  "gamm" "list"
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:  “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 :)
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.