40
$\begingroup$

I want to assume that the sea surface temperature of the Baltic Sea is the same year after year, and then describe that with a function / linear model. The idea I had was to just input year as a decimal number (or num_months/12) and get out what the temperature should be about that time. Throwing it into lm() function in R, it doesn't recognize sinusoidal data so it just produces a straight line. So I put the sin() function within a I() bracket and tried a few values to manually fit the function, and that gets close to what I want. But the sea is warming up faster in the summer and then cooling off slower in the fall... So the model is wrong the first year, then gets more correct after a couple of years, and then in the future I guess it becomes more and more wrong again.

How can I get R to estimate the model for me, so I don't have to guess numbers myself? The key here is that I want it to produce the same values year after year, not just be correct for one year. If I knew more about math, maybe I could guesstimate it as something like a Poisson or Gaussian instead of sin(), but I don't know how to do that either. Any help to get closer to a good answer would be greatly appreciated.

Here is the data I use, and the code to show results so far:

# SST from Bradtke et al 2010
ToY <- c(1/12,2/12,3/12,4/12,5/12,6/12,7/12,8/12,9/12,10/12,11/12,12/12,13/12,14/12,15/12,16/12,17/12,18/12,19/12,20/12,21/12,22/12,23/12,24/12,25/12,26/12,27/12,28/12,29/12,30/12,31/12,32/12,33/12,34/12,35/12,36/12,37/12,38/12,39/12,40/12,41/12,42/12,43/12,44/12,45/12,46/12,47/12,48/12)
Degrees <- c(3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5)
SST <- data.frame(ToY, Degrees)
SSTlm <- lm(SST$Degrees ~ I(sin(pi*2.07*SST$ToY)))
summary(SSTlm)
plot(SST,xlim=c(0,4),ylim=c(0,17))
par(new=T)
plot(data.frame(ToY=SST$ToY,Degrees=8.4418-6.9431*sin(2.07*pi*SST$ToY)),type="l",xlim=c(0,4),ylim=c(0,17))
$\endgroup$
48
$\begingroup$

It can be done with linear regression -

You just need both a $\sin$ and a $\cos$ term at each frequency.

The reason why you can use a $\sin$ and $\cos$ term in a linear regression to handle seasonality with any amplitude and phase is because of the following trigonometric identity:

A 'general' sine wave with amplitude $A$ and phase $\varphi$, $A \sin (x + \varphi)$, can be written as the linear combination $a\sin x+b\cos x$ where $a$ and $b$ are such that $A=\sqrt{a^2+b^2}$ and $\sin\varphi = \frac{b}{\sqrt{a^2+b^2}}$. Let's see that the two are equivalent:

\begin{eqnarray} a \sin(x) + b \cos(x) &=& \sqrt{a^2+b^2} \left(\frac{a}{\sqrt{a^2+b^2}} \sin(x) + \frac{b}{\sqrt{a^2+b^2}} \cos(x)\right)\\ &=& A\left[\sin(x)\cos(\varphi) + \cos(x)\sin(\varphi)\right]\\ &=& A\sin(x+\varphi)\,\text{.} \end{eqnarray}

Here's the 'basic' model:

 SSTlm <- lm(Degrees ~ sin(2*pi*ToY)+cos(2*pi*ToY),data=SST)
 summary(SSTlm)

[snip]

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)              8.292      0.135   61.41   <2e-16 *** 
sin(2 * pi * ToY)       -5.916      0.191  -30.98   <2e-16 ***  
cos(2 * pi * ToY)       -4.046      0.191  -21.19   <2e-16 *** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.9355 on 45 degrees of freedom
Multiple R-squared: 0.969,      Adjusted R-squared: 0.9677 
F-statistic: 704.3 on 2 and 45 DF,  p-value: < 2.2e-16 

 plot(Degrees~ToY,ylim=c(1.5,16.5),data=SST)
 lines(SST$ToY,SSTlm$fitted,col=2)

sin fit

Edit: Important note - the $2\pi\,t$ term works because the period of the function has been set up so that one period = 1 unit of $t$. If the period is different from 1, say the period is $\omega$, then you need $(2\pi/\omega)\, t$ instead.

Here's the model with the second harmonic:

 SSTlm2 <- lm(Degrees ~ sin(2*pi*ToY)+cos(2*pi*ToY)
                        +sin(4*pi*ToY)+cos(4*pi*ToY),data=SST)
 summary(SSTlm2)

[snip]

Coefficients:
                  Estimate Std. Error  t value Pr(>|t|)    
(Intercept)        8.29167    0.02637  314.450  < 2e-16 ***  
sin(2 * pi * ToY) -5.91562    0.03729 -158.634  < 2e-16 ***  
cos(2 * pi * ToY) -4.04632    0.03729 -108.506  < 2e-16 ***  
sin(4 * pi * ToY)  1.21244    0.03729   32.513  < 2e-16 ***  
cos(4 * pi * ToY)  0.33333    0.03729    8.939 2.32e-11 ***  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1827 on 43 degrees of freedom
Multiple R-squared: 0.9989,     Adjusted R-squared: 0.9988 
F-statistic:  9519 on 4 and 43 DF,  p-value: < 2.2e-16 

 plot(Degrees~ToY,ylab="Degrees",xlab="ToY",ylim=c(1.5,16.5),data=SST)
 lines(SSTlm2$fitted~ToY,col=2,data=SST)

sin fit 2

... and so forth, with 6*pi*ToY etc. If there was a tiny bit of noise in the data I'd probably stop with this second model though.

With enough terms, you can exactly fit asymmetric and even jagged periodic sequences, but the resulting fits may 'wiggle'. Here's an asymmetric function (it's a sawtooth -sawtooth) added to a scaled version of your periodic function), with third (red) and fourth (green) harmonics. The green fit is on average a little closer but "wiggly" (even when the fit goes through every point, the fit may be very wiggly between points).

sin fit 3&4

The periodicity here means there's only 12 d.f. available for a seasonal model in the data. With the intercept in the model, you only have enough degrees of freedom for 11 additional seasonal parameters. Since you are adding two terms with each harmonic, the last harmonic you can fit will only allow you one of them for the last term, the sixth harmonic (and that one has to be a $\cos$; the $\sin$ term will be all-zero, while the cos alternates between 1 and -1).

If you want fits that are smoother than this approach produces on non-smooth series, you may want to look into periodic spline fits.

Yet another approach is to use seasonal dummies, but the sin/cos approach is often better if it's a smooth periodic function.

This kind of approach to seasonality can also adapt to situations where seasonality is changing, such as using trigonometric or dummy seasonality with state-space models.


While the linear model approach discussed here is simple to use, one advantage of @COOLSerdash's nonlinear regression approach is that it can deal with a much wider range of situations - you don't have to change much before you're in a situation where linear regression is no longer suitable but nonlinear least-squares may still be used (having an unknown period would be one such case).

$\endgroup$
11
  • $\begingroup$ Awesome! Thank you, I really should try to learn more about methods to deal with frequencies. I don't quite understand why the cos part is needed, but knowing the principle makes it easy to implement. $\endgroup$ – GaRyu May 31 '13 at 9:17
  • $\begingroup$ @COOLSerdash - actually, I wish you hadn't deleted your answer (indeed I upvoted it); it has the advantage of working in a much wider array of circumstances; tweak a few things about the problem and you can lose linearity - and then my approach is useless, but yours still works. I think there's a lot to be said for being able to do it that way. $\endgroup$ – Glen_b May 31 '13 at 9:24
  • $\begingroup$ @Glen_b Ah sorry, I thought that your post made mine redundant because I didn't use the standard way of dealing with the problem. I undeleted it. $\endgroup$ – COOLSerdash May 31 '13 at 9:29
  • $\begingroup$ @GaRyu see my edit, near the top of my answer, where I give an outline of why adding in the $\cos$ does the trick. $\endgroup$ – Glen_b May 31 '13 at 9:50
  • 1
    $\begingroup$ That wasn't me.... You say phase offset as if that named what is going on, and it does mathematically. But for you the key point is more likely to be that Dec 31/Jan 1 is an arbitrary origin for time of year given lags in temperature response to variations in radiation receipt. So phase offset is a name here for something climatological too, the timing of minimum and maximum temperature relative to your recording system. (It's a minor detail but I prefer quantifying time of year for 12 months as 1/24, 3/24, ..., 23/24.) $\endgroup$ – Nick Cox May 31 '13 at 12:15
11
$\begingroup$

The temperature you provide in your question repeats exactly every year. I suspect this aren't really measured temperatures over four years. In your example, you wouldn't need a model, because the temperatures just repeat exactly. But otherwise you could use the nls function to fit a sine curve:

ToY <- c(1/12,2/12,3/12,4/12,5/12,6/12,7/12,8/12,9/12,10/12,11/12,12/12,13/12,14/12,15/12,16/12,17/12,18/12,19/12,20/12,21/12,22/12,23/12,24/12,25/12,26/12,27/12,28/12,29/12,30/12,31/12,32/12,33/12,34/12,35/12,36/12,37/12,38/12,39/12,40/12,41/12,42/12,43/12,44/12,45/12,46/12,47/12,48/12)
Degrees <- c(3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5)
SST <- data.frame(ToY, Degrees)

par(cex=1.5, bg="white")
plot(Degrees~ToY,xlim=c(0,4),ylim=c(0,17), pch=16, las=1)

nls.mod <-nls(Degrees ~ a + b*sin(2*pi*c*ToY), start=list(a = 1, b = 1, c=1))

co <- coef(nls.mod) 
f <- function(x, a, b, c) {a + b*sin(2*pi*c*x) }

curve(f(x, a=co["a"], b=co["b"], c=co["c"]), add=TRUE ,lwd=2, col="steelblue")

NLS fit

But the fit isn't very good, especially at the beginning. It seems that your data cannot be adequately modelled by a simple sine curve. Maybe a more complex trigonometrical function will do the trick?

nls.mod2 <-nls(Degrees ~ a + b*sin(2*pi*c*ToY)+d*cos(2*pi*e*ToY), start=list(a = 1, b = 1, c=1, d=1, e=1))

co2 <- coef(nls.mod2) 
f <- function(x, a, b, c, d, e) {a + b*sin(2*pi*c*x)+d*cos(2*pi*e*x) }

curve(f(x, a=co2["a"], b=co2["b"], c=co2["c"], d=co2["d"], e=co2["e"]), add=TRUE ,lwd=2, col="red")

NLS fit 2

The red curve fits the data better. With the nls function, you can put in the model that you think is appropriate.

Or maybe you could make use the the forecast package. In the example below, I have assumed that the time series started in January 2010:

library(forecast)

Degrees.ts <- ts(Degrees, start=c(2010,1), frequency=12)

Degree.trend <- auto.arima(Degrees.ts)

degrees.forecast <- forecast(Degree.trend, h=12, level=c(80,95), fan=F)

plot(degrees.forecast, las=1, main="", xlab="Time", ylab="Degrees")

ARIMA

Because the data is deterministic, no confidence bands are shown.

$\endgroup$
6
  • 5
    $\begingroup$ There is no reason for nonlinear least-squares here, not that it won't work reasonably well. Calculate sin(2 * pi * ToY), cos(2 * pi * ToY) in advance and feed them to lm() just like any other predictors. In other words, lm() need not see any trigonometry at all. However, you may need another model to capture marked asymmetry well. I am not a regular R user but I've often used this approach elsewhere (see stata-journal.com/sjpdf.html?articlenum=st0116). $\endgroup$ – Nick Cox May 31 '13 at 7:31
  • $\begingroup$ @NickCox Thanks Nick, that's very helpful advice. I will update my answer in a bit. $\endgroup$ – COOLSerdash May 31 '13 at 7:35
  • $\begingroup$ Glen was faster :) $\endgroup$ – COOLSerdash May 31 '13 at 7:48
  • 1
    $\begingroup$ @COOLserdash I didn't even see Nick Cox's comment there; it came while I was generating my answer. (This approach is pretty obvious if you've seen any Fourier series.) $\endgroup$ – Glen_b May 31 '13 at 7:58
  • 3
    $\begingroup$ As @Glen_b implies, this is a standard approach, just not universally known. $\endgroup$ – Nick Cox May 31 '13 at 8:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.