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… et nous allons découvrir quelques fonctions utilies dans ce contexte pour les séries chronologiques.

Comme nous l’avons déjà vu, le graphe de base (obtenu par plot()) est une ligne brisée qui relie les différents points de la séries, avec le temps sur l’axe des abscisses, et l’axe des ordonnées représentant les observation de la variable quantitative mesurée au fil du temps. Le choix d’utiliser une ligne brisée qui relie les points au lieu d’une nuage de points dans le cas présent est délibéré : les segments de droites qui relient les points matérialisent la relation entre eux, c’est-à-dire, l’existence d’une interdépendance entre ces observations.

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 temporelle 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 groupe de 10 stations, les imprimer et faire un graphique de la “moyenne glissante” à l’aide des instructions suivantes :

library(pastecs)
marbio <- ts(read("marbio", package = "pastecs"))
statsl <- stat.slide(1:68, marbio[, "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.4231
plot(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[, "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+03
plot(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é :

statsl2$stat[c("mean", "median", "min", "max"), ]
#          [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
statsl2$y
#  [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
statsl2$x
#  [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
statsl2$xcut
# [1]  0 17 25 30 41 46 70
À vous de jouer !

Effectuez maintenant les exercices du tutoriel C04La_stat_slide (Statistiques glissantes).

BioDataScience3::run("C04La_stat_slide")

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ême 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 construre une série de données mensuelles 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 :

tser <- ts(sin((1:100) / 6 * pi) + rnorm(100, sd = 0.5), start = c(1998, 4), 
  frequency = 12)
tser
#              Jan         Feb         Mar         Apr         May         Jun
# 1998                                      0.50866418  0.49137435  1.06735226
# 1999 -1.00133606 -0.60315118  0.50978783  0.77068655  0.96721919  0.45190207
# 2000 -0.89586557 -0.64357732 -0.10595280 -0.15452341  0.74620082  0.54052216
# 2001 -1.87925525  0.29233418 -0.16451998 -0.24854969  1.76666170  0.84885954
# 2002 -1.55561175 -1.62245955  0.90775059  0.33646151  0.56603842  1.07829134
# 2003 -1.01395187 -0.25563117  0.28302433  1.74126810  1.14084862  0.94946647
# 2004 -0.74889821 -0.38303223 -0.87389794  0.71325345  1.15839849  0.81901088
# 2005 -0.32828743 -0.36370109  0.71687847 -0.51522619  0.91122437  0.92440137
# 2006 -1.71666568  0.37342409 -0.06235499  0.26600770  1.29986498  1.11581354
#              Jul         Aug         Sep         Oct         Nov         Dec
# 1998  1.01681560  0.72096206 -0.35896182 -0.37786153 -1.50990396 -0.53682306
# 1999  0.86305699  0.75284390  0.20850766 -0.95769063 -2.18768239 -1.40006630
# 2000  1.50177838  1.23268167 -0.28785286 -0.09577899 -1.22082755 -0.48826730
# 2001  2.07783483  0.96922023 -0.73671157  0.04406268 -0.52502562 -0.14592694
# 2002  0.32143148  0.70017313  0.35741249 -0.54061332 -1.22432460 -0.68802281
# 2003  0.89650846  1.93609038 -0.08725773  0.03834280 -0.45657037 -1.04534435
# 2004  0.82754482  0.56605294 -0.51691120 -0.94398869 -0.64701855 -1.05509701
# 2005  1.63354553  0.61192625  0.10040011  0.25238306 -1.04620868 -0.47318402
# 2006  1.05181960

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

À vous de jouer !

Effectuez maintenant les exercices du tutoriel C04Lb_ts_create (Création d’objet ts).

BioDataScience3::run("C04Lb_ts_create")

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 :

tsp(tser)
# [1] 1998.25 2006.50   12.00
start(tser)
# [1] 1998    4
end(tser)
# [1] 2006    7
frequency(tser)
# [1] 12

Nous 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 :

tser.cycle <- cycle(tser)
tser.cycle
#      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
boxplot(split(tser, tser.cycle), names = month.abb, col = "cornsilk")

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 de un mois en une série de pas de temps de un trimestre, en calculant des moyennes trimestrielles, nous ferons (puisqu’il y a quatre trimestres dans une année, nous indiquons 4 ici) :

aggregate(tser, 4, mean)
#              Qtr1         Qtr2         Qtr3         Qtr4
# 1998               0.689130266  0.459605283 -0.808196183
# 1999 -0.364899804  0.729935935  0.608136184 -1.515146439
# 2000 -0.548465230  0.377399857  0.815535731 -0.601624614
# 2001 -0.583813682  0.788990517  0.770114498 -0.208963292
# 2002 -0.756773570  0.660263757  0.459672367 -0.817653577
# 2003 -0.328852905  1.277194393  0.915113705 -0.487857306
# 2004 -0.668609459  0.896887605  0.292228853 -0.882034752
# 2005  0.008296654  0.440133181  0.781957299 -0.422336547
# 2006 -0.468532197  0.893895408

Notez 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 :

window(tser, start = c(1999, 1), end = c(2001, 12))
#              Jan         Feb         Mar         Apr         May         Jun
# 1999 -1.00133606 -0.60315118  0.50978783  0.77068655  0.96721919  0.45190207
# 2000 -0.89586557 -0.64357732 -0.10595280 -0.15452341  0.74620082  0.54052216
# 2001 -1.87925525  0.29233418 -0.16451998 -0.24854969  1.76666170  0.84885954
#              Jul         Aug         Sep         Oct         Nov         Dec
# 1999  0.86305699  0.75284390  0.20850766 -0.95769063 -2.18768239 -1.40006630
# 2000  1.50177838  1.23268167 -0.28785286 -0.09577899 -1.22082755 -0.48826730
# 2001  2.07783483  0.96922023 -0.73671157  0.04406268 -0.52502562 -0.14592694

La 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 fonction 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.

À vous de jouer !

Effectuez maintenant les exercices du tutoriel C04Lc_ts_manip (Manipulation des séries chronologiques).

BioDataScience3::run("C04Lc_ts_manip")
Exercez-vous !

Manipulez et décrivez les différentes séries d’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 un). Comment feriez-vous pour transformer cette série avec une unité de temps de la seconde ? Ou de la milliseconde ?