Comment trouver un bon ajustement pour un modèle semi-sinusoïdal en R?

37

Je veux supposer que la température de surface de la mer Baltique est la même année après année, puis la décrire avec un modèle fonction / linéaire. L'idée que j'avais était de simplement entrer l'année sous forme de nombre décimal (ou num_months / 12) et de déterminer la température à ce moment-là. En le lançant dans la fonction lm () dans R, il ne reconnaît pas les données sinusoïdales, il produit donc une ligne droite. J'ai donc placé la fonction sin () dans un crochet I () et essayé quelques valeurs pour l'ajuster manuellement, ce qui se rapproche de ce que je veux. Mais la mer se réchauffe plus rapidement en été et se refroidit plus lentement à l’automne ... Le modèle a donc tort la première année, puis est devenu plus correct au bout de quelques années, puis à l’avenir, il deviendra de plus en plus et plus faux encore.

Comment puis-je demander à R d'estimer le modèle pour moi, de sorte que je n'ai pas à deviner les chiffres moi-même? La clé ici est que je veux qu’il produise les mêmes valeurs année après année, et ne soit pas juste correct pour un an. Si j'en savais plus sur les mathématiques, je pourrais peut-être l'assimiler à quelque chose comme un poisson ou un gaussien au lieu de sin (), mais je ne sais pas comment le faire non plus. Toute aide pour se rapprocher d'une bonne réponse serait grandement appréciée.

Voici les données que j'utilise et le code pour afficher les résultats jusqu'à présent:

# 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))
GaRyu
la source

Réponses:

44

Cela peut être fait avec une régression linéaire -

Vous avez juste besoin d’un terme et d’un terme à chaque fréquence.cossincos

La raison pour laquelle vous pouvez utiliser un terme et dans une régression linéaire pour gérer la saisonnalité avec toute amplitude et phase est due à l' identité trigonométrique suivante :cossincos

Une onde sinusoïdale «générale» d'amplitude et de phase , , peut être écrite comme la combinaison linéaire où et sont tels que et . Voyons que les deux sont équivalents:φ Un péché ( x + φ ) un péché x + b cos x a b A = AφAsin(x+φ)asinx+bcosxab sinφ=bA=a2+b2sinφ=ba2+b2

asin(x)+bcos(x)=a2+b2(aa2+b2sin(x)+ba2+b2cos(x))=A[sin(x)cos(φ)+cos(x)sin(φ)]=Asin(x+φ).

Voici le modèle de base:

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

[couper]

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)

ajustement du péché

Éditer: Remarque importante - le terme fonctionne car la période de la fonction a été configurée de manière à ce qu'une période = 1 unité de . Si la période est différente de 1, disons que la période est , alors vous avez besoin de place.2πttω(2π/ω)t

Voici le modèle avec le deuxième harmonique:

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

[couper]

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

... et ainsi de suite, avec 6*pi*ToYetc. S'il y avait un peu de bruit dans les données, je m'arrêterais probablement avec ce deuxième modèle.

Avec suffisamment de termes, vous pouvez ajuster exactement des séquences périodiques asymétriques et même irrégulières, mais les ajustements qui en résultent risquent de "bouger". Voici une fonction asymétrique (c'est une dent de scie - en dents de scieajoutée à une version mise à l'échelle de votre fonction périodique), avec une troisième harmonique (rouge) et une quatrième harmonique (verte). La coupe verte est en moyenne un peu plus serrée mais "ondulée" (même si la crise passe par tous les points, elle peut être très ondulée entre les points).

sin fit 3 et 4

La périodicité signifie ici qu'il n'y a que 12 df disponibles pour un modèle saisonnier dans les données. Avec l'interception dans le modèle, vous ne disposez que de suffisamment de degrés de liberté pour 11 paramètres saisonniers supplémentaires. Puisque vous ajoutez deux termes à chaque harmonique, le dernier harmonique que vous pouvez adapter ne vous permettra que l’ un d’eux pour le dernier terme, le sixième harmonique (et il faut que ce soit un ; le terme zéro, tandis que le cos alterne entre 1 et -1).cossin

Si vous souhaitez des ajustements plus lisses que ce que cette approche produit sur des séries non lisses, vous souhaiterez peut-être examiner l' ajustement périodique des splines .

Une autre approche consiste à utiliser des variables nominales saisonnières, mais l'approche sin / cos est souvent préférable s'il s'agit d'une fonction périodique lisse.

Ce type d’approche de la saisonnalité peut également s’adapter à des situations dans lesquelles la saisonnalité évolue, comme l’utilisation de la saisonnalité trigonométrique ou fictive avec des modèles d’espace à états.


Bien que l’approche du modèle linéaire décrite ici soit simple à utiliser, l’un des avantages de la méthode de régression non linéaire de @ COOLSerdash est qu’elle peut traiter un plus grand nombre de situations - vous n’avez pas à changer beaucoup avant de vous retrouver dans une situation où la régression ne convient plus, mais des moindres carrés non linéaires peuvent toujours être utilisés (une période inconnue en serait un exemple).

Glen_b -Reinstate Monica
la source
Impressionnant! Merci, je devrais vraiment essayer d’en apprendre davantage sur les méthodes de traitement des fréquences. Je ne comprends pas très bien pourquoi le cos est nécessaire, mais connaître le principe le rend facile à mettre en œuvre.
GaRyu
@COOLSerdash - en fait, j'aimerais que vous n'ayez pas supprimé votre réponse (en fait, je l'ai votée); il a l'avantage de fonctionner dans un éventail de situations beaucoup plus large; modifiez quelques points à propos du problème et vous risquez de perdre la linéarité - et alors mon approche est inutile, mais la vôtre fonctionne toujours. Je pense qu'il y a beaucoup à dire pour pouvoir le faire de cette façon.
Glen_b -Reinstate Monica
@Glen_b Ah désolé, je pensais que votre message rendait le mien superflu, car je n'utilisais pas la méthode standard pour traiter le problème. Je l'ai supprimé.
COOLSerdash
@GaRyu vois mon édition, vers le haut de ma réponse, où je donne un aperçu de la raison pour laquelle l'ajout de fait l'affaire. cos
Glen_b -Reinstate Monica
1
Ce n'était pas moi ... Vous dites décalage de phase comme si cela nommait ce qui se passait, et c'est mathématique. Mais pour vous, le point le plus important est probablement que le 31 décembre / 1er janvier est une origine arbitraire pour la période de l’année étant donné que la réponse en température est décalée par rapport aux variations de la réception du rayonnement. Le décalage de phase est donc ici un nom pour quelque chose de climatologique, le réglage de la température minimale et maximale par rapport à votre système d’enregistrement. (C'est un détail mineur, mais je préfère quantifier la période de l'année sur 12 mois: 1/24, 3/24, ..., 23/24.)
Nick Cox
10

La température que vous indiquez dans votre question se répète exactement chaque année. Je soupçonne que ce ne sont pas vraiment des températures mesurées sur quatre ans. Dans votre exemple, vous n’auriez pas besoin d’un modèle, car les températures se répètent exactement. Mais sinon, vous pourriez utiliser la nlsfonction pour ajuster une courbe sinusoïdale:

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")

Ajustement NLS

Mais l'ajustement n'est pas très bon, surtout au début. Il semble que vos données ne puissent pas être modélisées de manière adéquate par une simple courbe sinusoïdale. Peut-être qu'une fonction trigonométrique plus complexe fera l'affaire?

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

La courbe rouge correspond mieux aux données. Avec la nlsfonction, vous pouvez insérer le modèle que vous jugez approprié.

Ou peut-être pourriez-vous utiliser le forecastpaquet. Dans l'exemple ci-dessous, j'ai supposé que la série chronologique avait commencé en janvier 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

Comme les données sont déterministes, aucune bande de confiance n’est montrée.

COOLSerdash
la source
4
Il n'y a aucune raison pour les moindres carrés non linéaires ici, pas que cela ne fonctionnera pas raisonnablement bien. Calculez à l'avance le péché (2 * pi * ToY), cos (2 * pi * ToY) et indiquez lm()-le comme tout autre prédicteur. En d'autres termes, il lm()n'est pas nécessaire de voir la trigonométrie du tout. Cependant, vous aurez peut-être besoin d'un autre modèle pour bien capturer une asymétrie marquée. Je ne suis pas un utilisateur régulier de R, mais j'ai souvent utilisé cette approche ailleurs (voir stata-journal.com/sjpdf.html?articlenum=st0116 ).
Nick Cox
@ NickCox Merci Nick, c'est un conseil très utile. Je mettrai à jour ma réponse dans un instant.
COOLSerdash
Glen était plus rapide :)
COOLSerdash
1
@COOLserdash Je n'y ai même pas vu le commentaire de Nick Cox; il est venu pendant que je produisais ma réponse. (Cette approche est assez évidente si vous avez vu une série de Fourier.)
Glen_b -Reinstate Monica
2
Comme @Glen_b l’implique, c’est une approche standard, mais qui n’est pas universellement connue.
Nick Cox