4.3 Manipulation et description
Tout comme pour les statistiques plus classiques, nous commençons toujours notre exploration de séries chronologiques par une partie purement descriptive, soit numérique, soit graphique. De plus, nous savons que nous ne devons pas espérer obtenir d’emblée des données qui se présentent de manière parfaite. Un petit peu de manipulation de données est parfois nécessaire… Nous allons découvrir quelques fonctions utiles dans ce contexte pour les séries chronologiques.
Comme nous l’avons déjà vu, le graphe de base (obtenu par plot() ou en utilisant geom_line()) est une ligne brisée qui relie les différents points de la série, avec le temps sur l’axe des abscisses, et l’axe des ordonnées représentant les observations de la variable quantitative mesurée au fil du temps.
4.3.1 Statistiques glissantes
Les descripteurs statistiques tels que la moyenne, la médiane, l’écart type, etc. restent utiles pour les séries temporelles. Par contre, la valeur de ces descripteurs est susceptible de varier au cours du temps, et cette variation donne des indications utiles sur les mécanismes biologiques sous-jacents (par exemple, des fluctuations, des changements brutaux de régimes, etc.). Il nous faut donc calculer ces descripteurs pour des intervalles de temps précis le long de l’axe temporel de notre série. Les statistiques glissantes (sliding statistics en anglais) correspondent précisément à l’analyse de blocs successifs de données suivant un axe spatio-temporel. Dans {pastecs}, elles se calculent à l’aide de la fonction stat.slide().
Exemple
Nous reprenons ici l’exemple du manuel de {pastecs}. On peut calculer des statistiques glissantes pour la série ClausocalanusA de marbio par groupes de 10 stations, les imprimer et faire un graphique de la “moyenne glissante” à l’aide des instructions suivantes :
marbio_ts <- ts(read("marbio", package = "pastecs"))
statsl <- stat.slide(1:68, marbio_ts[, "ClausocalanusA"],
xmin = 0, n = 7, deltat = 10)
statsl# [0,10[ [10,20[ [20,30[ [30,40[ [40,50[ [50,60[ [60,70[
# xmin 0.000000 10.0000 20.0000 30.0000 40.0000 50.000 60.0000
# xmax 10.000000 20.0000 30.0000 40.0000 50.0000 60.000 70.0000
# nbr.val 9.000000 10.0000 10.0000 10.0000 10.0000 10.000 9.0000
# nbr.null 2.000000 0.0000 0.0000 0.0000 0.0000 0.000 0.0000
# nbr.na 0.000000 0.0000 0.0000 0.0000 0.0000 0.000 0.0000
# min 0.000000 160.0000 158.0000 96.0000 136.0000 56.000 52.0000
# max 12.000000 832.0000 1980.0000 1204.0000 2484.0000 824.000 656.0000
# median 4.000000 344.0000 732.0000 752.0000 260.0000 340.000 120.0000
# mean 4.777778 350.8000 810.1000 682.0000 494.4000 364.800 219.1111
# std.dev 4.790036 196.9708 619.8578 350.0616 711.0511 234.915 202.4231plot(statsl, stat = "mean", leg = TRUE, lpos = c(55, 2500),
xlab = "Station", ylab = "ClausocalanusA")
Toujours pour la même série, on peut calculer toutes les statistiques sur des intervalles irréguliers, puis représenter l’étendue (minimum, maximum) et la médiane pour chaque intervalle comme suit (ces intervalles correspondent aux différentes masses d’eaux croisées le long du transect, voir légende du graphique) :
statsl2 <- stat.slide(1:68, marbio_ts[, "ClausocalanusA"],
xcut = c(0, 17, 25, 30, 41, 46, 70),
basic = TRUE, desc = TRUE, norm = TRUE, pen = TRUE, p = 0.95)
statsl2# [0,17[ [17,25[ [25,30[ [30,41[
# xmin 0.000000e+00 1.700000e+01 2.500000e+01 3.000000e+01
# xmax 1.700000e+01 2.500000e+01 3.000000e+01 4.100000e+01
# nbr.val 1.600000e+01 8.000000e+00 5.000000e+00 1.100000e+01
# nbr.null 2.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
# nbr.na 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
# min 0.000000e+00 1.580000e+02 1.040000e+03 9.600000e+01
# max 8.320000e+02 4.240000e+02 1.980000e+03 1.204000e+03
# range 8.320000e+02 2.660000e+02 9.400000e+02 1.108000e+03
# sum 2.785000e+03 2.151000e+03 6.716000e+03 7.060000e+03
# median 1.200000e+01 2.120000e+02 1.204000e+03 7.440000e+02
# mean 1.740625e+02 2.688750e+02 1.343200e+03 6.418182e+02
# SE.mean 6.066894e+01 3.704796e+01 1.668120e+02 1.078927e+02
# CI.mean.0.95 1.293128e+02 8.760450e+01 4.631443e+02 2.403999e+02
# var 5.889153e+04 1.098041e+04 1.391312e+05 1.280492e+05
# std.dev 2.426758e+02 1.047875e+02 3.730029e+02 3.578396e+02
# coef.var 1.394188e+00 3.897255e-01 2.776973e-01 5.575404e-01
# skewness 1.248551e+00 4.181831e-01 8.469063e-01 -4.141160e-02
# skew.2SE 1.106268e+00 2.780098e-01 4.638697e-01 -3.133978e-02
# kurtosis 7.020336e-01 -1.852389e+00 -1.193537e+00 -1.424298e+00
# kurt.2SE 3.218053e-01 -6.254351e-01 -2.983842e-01 -5.566203e-01
# normtest.W 7.557989e-01 8.300443e-01 8.162687e-01 9.560890e-01
# normtest.p 7.467713e-04 5.943034e-02 1.092399e-01 7.221457e-01
# pos.median 8.600000e+01 2.120000e+02 1.204000e+03 7.440000e+02
# pos.mean 1.989286e+02 2.688750e+02 1.343200e+03 6.418182e+02
# geo.mean 3.935062e+01 2.522408e+02 1.307640e+03 5.181136e+02
# pen.mean 3.351203e+02 2.683489e+02 1.340628e+03 6.787440e+02
# pen.var 1.907110e+06 1.039770e+04 1.119418e+05 3.111941e+05
# pen.std.dev 1.380981e+03 1.019691e+02 3.345771e+02 5.578477e+02
# pen.mean.var 5.558396e+04 1.298434e+03 2.238628e+04 2.768567e+04
# [41,46[ [46,70[
# xmin 41.0000000 4.600000e+01
# xmax 46.0000000 7.000000e+01
# nbr.val 5.0000000 2.300000e+01
# nbr.null 0.0000000 0.000000e+00
# nbr.na 0.0000000 0.000000e+00
# min 136.0000000 5.200000e+01
# max 264.0000000 2.484000e+03
# range 128.0000000 2.432000e+03
# sum 1088.0000000 9.236000e+03
# median 256.0000000 3.080000e+02
# mean 217.6000000 4.015652e+02
# SE.mean 27.2939554 1.048611e+02
# CI.mean.0.95 75.7801688 2.174685e+02
# var 3724.8000000 2.529043e+05
# std.dev 61.0311396 5.028960e+02
# coef.var 0.2804740 1.252339e+00
# skewness -0.3597825 3.033043e+00
# skew.2SE -0.1970610 3.150646e+00
# kurtosis -2.1086492 9.865540e+00
# kurt.2SE -0.5271623 5.277023e+00
# normtest.W 0.7864078 6.024972e-01
# normtest.p 0.0625007 9.610838e-07
# pos.median 256.0000000 3.080000e+02
# pos.mean 217.6000000 4.015652e+02
# geo.mean 209.9226025 2.527195e+02
# pen.mean 218.0553885 3.905421e+02
# pen.var 4518.0257036 1.913652e+05
# pen.std.dev 67.2162607 4.374531e+02
# pen.mean.var 903.3938793 7.745030e+03plot(statsl2, stat = "median", xlab = "Stations", ylab = "Effectifs",
main = "Clausocalanus A") # Médiane
lines(statsl2, stat = "min") # Minimum
lines(statsl2, stat = "max") # Maximum
lines(c(17, 17), c(-50,2600), col = 4, lty = 2) # Séparations des masses d'eaux
lines(c(25, 25), c(-50,2600), col = 4, lty = 2)
lines(c(30, 30), c(-50,2600), col = 4, lty = 2)
lines(c(41, 41), c(-50,2600), col = 4, lty = 2)
lines(c(46, 46), c(-50,2600), col = 4, lty = 2)
text(c(8.4, 21, 27.5, 35, 43.5, 57.2), 2300,
labels = c("Zone périphérique", "D1", "C", "Front", "D2", "Zone centrale")) # Labels
legend(0, 1900, c("série", "médian", "étendue"), col = 1:3, lty = 1)
Enfin, on peut extraire différentes statistiques, les valeurs de la série initiale (y), les valeurs de temps de la série (x), et le vecteur temps de coupure des périodes (xcut) à partir de l’objet stat.slide renvoyé :
# [0,17[ [17,25[ [25,30[ [30,41[ [41,46[ [46,70[
# mean 174.0625 268.875 1343.2 641.8182 217.6 401.5652
# median 12.0000 212.000 1204.0 744.0000 256.0 308.0000
# min 0.0000 158.000 1040.0 96.0000 136.0 52.0000
# max 832.0000 424.000 1980.0 1204.0000 264.0 2484.0000# [1] 0 4 2 8 4 0 1 12 12 196 160 406 448 328 372
# [16] 832 360 200 206 218 158 193 424 392 1980 1348 1144 1204 1040 1204
# [31] 1080 424 744 96 888 608 776 240 760 240 256 264 136 264 168
# [46] 516 136 480 2484 824 376 352 160 640 280 488 144 328 56 112
# [61] 408 308 68 52 84 164 656 120# [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
# [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
# [51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68# [1] 0 17 25 30 41 46 70À vous de jouer !
Effectuez maintenant les exercices du tutoriel C04La_ts_intro (Introduction aux séries chronologiques).
BioDataScience3::run("C04La_ts_intro")4.3.2 Manipulations de ts
Les jeux de données bricolés artificiellement sont parfois bien utiles pour explorer de nouvelles techniques. En effet, puisque nous les construisons nous-mêmes, nous savons dès le départ de quoi ils sont faits. Et donc, nous pouvons vérifier que le résultat obtenu est bien conforme. Par exemple, nous pouvons construire une série de données mensuelle commençant en avril 1998 (start = c(1998, 4)), contenant 100 valeurs et comportant une composante sinusoïdale d’une période annuelle ainsi qu’une composante aléatoire ayant une distribution Normale de moyenne nulle et d’écart type 0.5 :
# Jan Feb Mar Apr May Jun
# 1998 1.48630078 0.56621520 1.57091092
# 1999 -0.09873557 -1.34364712 0.70308821 1.40131844 1.64247788 1.00379275
# 2000 -0.80396834 -0.33091303 -0.26818113 1.11392471 0.62067360 0.36806062
# 2001 -0.17172956 -0.68706093 0.11277892 0.60971508 1.28549977 1.32428601
# 2002 -1.35724550 -1.14569981 -0.17574209 0.24383190 0.48141119 0.39119521
# 2003 -0.91051530 -0.84930152 -0.40608321 0.24615059 0.30719730 1.61492378
# 2004 -0.75535281 -0.10220204 -0.09589435 0.43469167 -0.36483508 0.25936417
# 2005 -0.64354205 -0.56978719 -0.06054512 0.37333726 0.56956888 1.32893013
# 2006 -1.21564155 -0.31461409 0.76251652 0.24327239 0.72441114 1.22925712
# Jul Aug Sep Oct Nov Dec
# 1998 0.93872568 0.32163815 0.38611825 -1.27957234 -0.69607075 -0.51056820
# 1999 0.07782303 0.07160203 -0.16533890 -1.02251410 -1.34738011 -1.02040394
# 2000 -0.57223290 0.36600443 -0.88374263 -1.54797313 0.26505529 -1.04725899
# 2001 0.64012054 0.86225793 0.23649716 -1.01348947 -0.79993086 -0.65016015
# 2002 0.72213869 0.65964592 -0.20911195 -0.79500964 -0.47702753 -0.76954370
# 2003 0.88187397 0.41665300 0.37422613 -0.77932827 -0.38461946 -0.45941514
# 2004 0.61963082 0.83463118 0.72905351 0.03245230 -0.62866707 -0.55709260
# 2005 0.55346512 0.70731697 0.27609633 0.22317115 0.08664955 -1.11204005
# 2006 0.42253963Créons maintenant une série multiple en prenant notre série de départ, ainsi qu’une version décalée de cinq mois à l’aide de la fonction stats::lag() et représentons cela graphiquement :
mtser <- ts.intersect(tser, stats::lag(tser, 5))
plot(mtser, xlab = "Temps", main = "Série tser et série décalée de 5 mois")
4.3.2.1 Manipulation du temps
Les fonctions tsp(), start(), end() et frequency() permettent de récupérer les différents paramètres de temps de la série :
# [1] 1998.25 2006.50 12.00# [1] 1998 4# [1] 2006 7# [1] 12Nous avons déjà vu la fonction time() qui permet de reconstituer le vecteur temps pour la série. La fonction cycle() indique l’appartenance de chaque donnée de la série à un cycle. Elle permet, par exemple, de séparer les données par mois si la base de temps est l’année à l’aide de la fonction split(). Ensuite, il est possible de traiter ou de représenter séparément les statistiques mois par mois pour en faire, par exemple, des boites de dispersions parallèles :
# Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
# 1998 4 5 6 7 8 9 10 11 12
# 1999 1 2 3 4 5 6 7 8 9 10 11 12
# 2000 1 2 3 4 5 6 7 8 9 10 11 12
# 2001 1 2 3 4 5 6 7 8 9 10 11 12
# 2002 1 2 3 4 5 6 7 8 9 10 11 12
# 2003 1 2 3 4 5 6 7 8 9 10 11 12
# 2004 1 2 3 4 5 6 7 8 9 10 11 12
# 2005 1 2 3 4 5 6 7 8 9 10 11 12
# 2006 1 2 3 4 5 6 7
À vous de jouer !
Une autre manipulation courante de l’axe temporel consiste à agréger les données en des pas de temps plus longs. aggregate() permet de réduire le nombre d’observations en diminuant la fréquence. Par exemple pour transformer la série tser qui a un pas de temps d’un mois en une série de pas de temps d’un trimestre, en calculant des moyennes trimestrielles, nous ferons (puisqu’il y a quatre trimestres dans une année, nous indiquons 4 ici) :
# Qtr1 Qtr2 Qtr3 Qtr4
# 1998 1.207808964 0.548827357 -0.828737098
# 1999 -0.246431494 1.349196355 -0.005304612 -1.130099385
# 2000 -0.467687502 0.700886310 -0.363323700 -0.776725610
# 2001 -0.248670523 1.073166956 0.579625210 -0.821193492
# 2002 -0.892895802 0.372146100 0.390890887 -0.680526956
# 2003 -0.721966677 0.722757223 0.557584367 -0.541120956
# 2004 -0.317816400 0.109740249 0.727771839 -0.384435788
# 2005 -0.424624787 0.757278758 0.512292803 -0.267406451
# 2006 -0.255913038 0.732313549Notez encore une fois que R assume que l’unité de temps est l’année et il crée des intitulés particuliers pour des séries de fréquence égale à 4 (Qtr1 -> Qtr4) ou à 12 (intitulé des mois en abrégé).
Si nous ne souhaitons pas utiliser la série sur toute sa longueur, nous pouvons employer window() pour extraire une sous-série contenue dans une fenêtre de temps spécifique sans modifier la fréquence des observations. Par exemple, pour extraire les années 1999 à 2001 complètes de tser, on utilisera :
# Jan Feb Mar Apr May Jun
# 1999 -0.09873557 -1.34364712 0.70308821 1.40131844 1.64247788 1.00379275
# 2000 -0.80396834 -0.33091303 -0.26818113 1.11392471 0.62067360 0.36806062
# 2001 -0.17172956 -0.68706093 0.11277892 0.60971508 1.28549977 1.32428601
# Jul Aug Sep Oct Nov Dec
# 1999 0.07782303 0.07160203 -0.16533890 -1.02251410 -1.34738011 -1.02040394
# 2000 -0.57223290 0.36600443 -0.88374263 -1.54797313 0.26505529 -1.04725899
# 2001 0.64012054 0.86225793 0.23649716 -1.01348947 -0.79993086 -0.65016015La fonction stats::lag() permet de décaler la série en arrière dans le temps (ou en avant si on fournit une valeur de décalage négative) ; nous l’avons déjà utilisée plus haut. La fonction diff() calcule la différence entre une série et elle-même décalée dans le temps de k observations. Les fonctions ts.intersect() et ts.union() permettent de réunir deux ou plusieurs séries au sein d’une même matrice multivariée, avec une échelle de temps commune. Les fréquences respectives des séries à assembler doivent être identiques. Là où ts.intersect() retient uniquement l’intervalle de temps commun à toutes les séries, ts.union() conservera l’ensemble de l’intervalle temporel, toutes séries confondues.
Exercez-vous !
Manipulez et décrivez les différentes séries en exemples. Notez par exemple que EEG est constituée de mesures effectuées à un rythme de 256 par seconde, mais que l’axe du temps reprend simplement les observations les unes après les autres (ou si vous préférez, une “unité” de temps de 1/256 sec et une fréquence de 1). Comment feriez-vous pour transformer cette série avec une unité de temps de la seconde ? Ou de la milliseconde ?