Manipulación de datos de serie temporal en R
F. Tusell
Esta versión generada el: 2021-11-02
Lectura
Veremos a continuación cómo leer series en R y fecharlas. Por el momento, trabajaremos sólo con series mensuales; en R, sin embargo, podemos trabajar con series de cualquier frecuencia de observación —horaria, cada minuto, cada segundo, cada trimestre— e incluso con series que no son regularmente observadas (procesos puntuales).
setwd("/home/etptupaf/ft/docencia/MasterModMat/MaterialClase/Ejemplos")
system(command="head -n 12 AccidentesMortales.csv", intern=TRUE)## [1] "Transportes"
## [2] " Otras informaciones relacionadas con el transporte terrestre"
## [3] ""
## [4] "Accidentes de circulaci\xf3n y v\xedctimas"
## [5] "Units:unidades"
## [6] ""
## [7] ""
## [8] ";Accidentes con v\xedctimas mortales. Total;Accidentes en zona urbana con v\xedctimas mortales;Accidentes en carretera con v\xedctimas mortales;"
## [9] "1995M01;\"321.0\";\"58.0\";\"263.0\";"
## [10] "1995M02;\"258.0\";\"46.0\";\"212.0\";"
## [11] "1995M03;\"354.0\";\"73.0\";\"281.0\";"
## [12] "1995M04;\"370.0\";\"66.0\";\"304.0\";"acc <- read.csv(file="AccidentesMortales.csv", skip=8, header=FALSE, sep=";",
encoding="latin1")
tail(acc, n=12)## V1
## 210 2012M06
## 211 2012M07
## 212 2012M08
## 213 2012M09
## 214 2012M10
## 215 2012M11
## 216 2012M12
## 217 Notes:
## 218 1.- Fuente de información: Dirección General de Tráfico
## 219 Source:Boletín Mensual de Estadística. INE
## 220 Copyright INE 2014
## 221 Paseo de la Castellana, 183 - 28071 - Madrid - España Teléfono: (+34) 91 583 91 00 - Contacta:www.ine.es/infoine
## V2 V3 V4 V5
## 210 136 35 101 NA
## 211 133 20 113 NA
## 212 146 23 123 NA
## 213 156 33 123 NA
## 214 116 29 87 NA
## 215 106 28 78 NA
## 216 120 33 87 NA
## 217 NA NA NA NA
## 218 NA NA NA NA
## 219 NA NA NA NA
## 220 NA NA NA NA
## 221 NA NA NA NAacc <- read.csv(file="AccidentesMortales.csv",skip=8,nrows=216,header=FALSE,sep=";")El resultado es una dataframe:
class(acc)## [1] "data.frame"Fechado
Formato ts
Nos interesará nombrar las columnas y convertir a un formato de serie temporal. El “standard” de R es ts, aunque es bastante limitado:
acc <- acc[,2:4]
colnames(acc) <- c("MuertosTot","MuertosUrb","MuertosCtra")
acc.ts <- ts(acc,start=c(1995,1),frequency=12)Tenemos ahora los datos fechados, lo que facilita su legibilidad y la interpetación de los gráficos.
acc.ts## MuertosTot MuertosUrb MuertosCtra
## Jan 1995 321 58 263
## Feb 1995 258 46 212
## Mar 1995 354 73 281
## Apr 1995 370 66 304
## May 1995 320 54 266
## Jun 1995 307 41 266
## Jul 1995 421 65 356
## Aug 1995 392 52 340
## Sep 1995 396 71 325
## Oct 1995 323 68 255
## Nov 1995 323 64 259
## Dec 1995 349 68 281
## Jan 1996 269 59 210
## Feb 1996 284 59 225
## Mar 1996 299 65 234
## Apr 1996 308 50 258
## May 1996 328 68 260
## Jun 1996 367 67 300
## Jul 1996 384 63 321
## Aug 1996 374 62 312
## Sep 1996 338 55 283
## Oct 1996 343 64 279
## Nov 1996 316 69 247
## Dec 1996 352 50 302
## Jan 1997 302 61 241
## Feb 1997 263 53 210
## Mar 1997 337 63 274
## Apr 1997 306 55 251
## May 1997 309 44 265
## Jun 1997 313 54 259
## Jul 1997 403 64 339
## Aug 1997 456 72 384
## Sep 1997 348 53 295
## Oct 1997 329 52 277
## Nov 1997 327 56 271
## Dec 1997 369 81 288
## Jan 1998 310 63 247
## Feb 1998 300 60 240
## Mar 1998 334 57 277
## Apr 1998 268 42 226
## May 1998 359 75 284
## Jun 1998 383 68 315
## Jul 1998 394 60 334
## Aug 1998 444 67 377
## Sep 1998 379 65 314
## Oct 1998 382 65 317
## Nov 1998 369 65 304
## Dec 1998 397 65 332
## Jan 1999 318 55 263
## Feb 1999 296 50 246
## Mar 1999 313 56 257
## Apr 1999 321 53 268
## May 1999 340 50 290
## Jun 1999 312 54 258
## Jul 1999 411 54 357
## Aug 1999 422 67 355
## Sep 1999 376 44 332
## Oct 1999 395 64 331
## Nov 1999 378 71 307
## Dec 1999 357 66 291
## Jan 2000 352 65 287
## Feb 2000 301 53 248
## Mar 2000 340 78 262
## Apr 2000 348 72 276
## May 2000 341 62 279
## Jun 2000 348 55 293
## Jul 2000 438 61 377
## Aug 2000 406 56 350
## Sep 2000 375 56 319
## Oct 2000 387 75 312
## Nov 2000 351 61 290
## Dec 2000 385 54 331
## Jan 2001 312 51 261
## Feb 2001 316 69 247
## Mar 2001 321 57 264
## Apr 2001 323 66 257
## May 2001 337 62 275
## Jun 2001 385 63 322
## Jul 2001 387 60 327
## Aug 2001 422 58 364
## Sep 2001 347 48 299
## Oct 2001 325 61 264
## Nov 2001 329 65 264
## Dec 2001 366 58 308
## Jan 2002 356 64 292
## Feb 2002 292 50 242
## Mar 2002 320 51 269
## Apr 2002 274 42 232
## May 2002 301 43 258
## Jun 2002 356 54 302
## Jul 2002 362 69 293
## Aug 2002 428 55 373
## Sep 2002 345 62 283
## Oct 2002 327 56 271
## Nov 2002 331 60 271
## Dec 2002 339 48 291
## Jan 2003 324 58 266
## Feb 2003 293 57 236
## Mar 2003 301 53 248
## Apr 2003 294 32 262
## May 2003 296 38 258
## Jun 2003 375 68 307
## Jul 2003 393 70 323
## Aug 2003 437 58 379
## Sep 2003 342 63 279
## Oct 2003 340 56 284
## Nov 2003 339 72 267
## Dec 2003 350 44 306
## Jan 2004 290 48 242
## Feb 2004 249 54 195
## Mar 2004 281 57 224
## Apr 2004 296 61 235
## May 2004 309 51 258
## Jun 2004 327 56 271
## Jul 2004 347 59 288
## Aug 2004 356 55 301
## Sep 2004 280 45 235
## Oct 2004 357 63 294
## Nov 2004 263 49 214
## Dec 2004 288 53 235
## Jan 2005 268 51 217
## Feb 2005 259 52 207
## Mar 2005 289 54 235
## Apr 2005 262 49 213
## May 2005 273 42 231
## Jun 2005 307 65 242
## Jul 2005 336 50 286
## Aug 2005 306 45 261
## Sep 2005 275 34 241
## Oct 2005 280 38 242
## Nov 2005 246 38 208
## Dec 2005 276 46 230
## Jan 2006 286 48 238
## Feb 2006 218 38 180
## Mar 2006 265 44 221
## Apr 2006 268 55 213
## May 2006 278 40 238
## Jun 2006 277 44 233
## Jul 2006 289 50 239
## Aug 2006 252 36 216
## Sep 2006 262 41 221
## Oct 2006 240 41 199
## Nov 2006 244 36 208
## Dec 2006 240 45 195
## Jan 2007 225 39 186
## Feb 2007 210 46 164
## Mar 2007 271 57 214
## Apr 2007 241 49 192
## May 2007 250 44 206
## Jun 2007 241 40 201
## Jul 2007 288 50 238
## Aug 2007 273 44 229
## Sep 2007 270 47 223
## Oct 2007 261 49 212
## Nov 2007 190 36 154
## Dec 2007 237 41 196
## Jan 2008 204 51 153
## Feb 2008 179 35 144
## Mar 2008 189 30 159
## Apr 2008 191 44 147
## May 2008 200 33 167
## Jun 2008 192 40 152
## Jul 2008 239 36 203
## Aug 2008 227 41 186
## Sep 2008 178 43 135
## Oct 2008 210 38 172
## Nov 2008 190 35 155
## Dec 2008 186 31 155
## Jan 2009 168 35 133
## Feb 2009 169 41 128
## Mar 2009 179 34 145
## Apr 2009 157 34 123
## May 2009 186 32 154
## Jun 2009 188 33 155
## Jul 2009 196 39 157
## Aug 2009 209 36 173
## Sep 2009 163 32 131
## Oct 2009 184 45 139
## Nov 2009 151 33 118
## Dec 2009 167 27 140
## Jan 2010 159 35 124
## Feb 2010 120 34 86
## Mar 2010 136 34 102
## Apr 2010 123 26 97
## May 2010 162 21 141
## Jun 2010 159 32 127
## Jul 2010 196 41 155
## Aug 2010 204 37 167
## Sep 2010 176 38 138
## Oct 2010 194 38 156
## Nov 2010 163 33 130
## Dec 2010 161 37 124
## Jan 2011 129 28 101
## Feb 2011 117 26 91
## Mar 2011 118 21 97
## Apr 2011 122 33 89
## May 2011 149 24 125
## Jun 2011 127 27 100
## Jul 2011 183 29 154
## Aug 2011 176 43 133
## Sep 2011 150 30 120
## Oct 2011 139 32 107
## Nov 2011 139 25 114
## Dec 2011 134 22 112
## Jan 2012 124 30 94
## Feb 2012 116 28 88
## Mar 2012 134 39 95
## Apr 2012 108 17 91
## May 2012 128 31 97
## Jun 2012 136 35 101
## Jul 2012 133 20 113
## Aug 2012 146 23 123
## Sep 2012 156 33 123
## Oct 2012 116 29 87
## Nov 2012 106 28 78
## Dec 2012 120 33 87Observemos que mandatos como head y tail no funcionan del mejor modo sobre un objeto de tipo ts:
head(acc.ts)## MuertosTot MuertosUrb MuertosCtra
## [1,] 321 58 263
## [2,] 258 46 212
## [3,] 354 73 281
## [4,] 370 66 304
## [5,] 320 54 266
## [6,] 307 41 266Si queremos un subconjunto y que el output conserve información sobre las fechas, debemos recurrir al mandato window:
window(acc.ts, start=c(1995,1), end=c(1995,8))## MuertosTot MuertosUrb MuertosCtra
## Jan 1995 321 58 263
## Feb 1995 258 46 212
## Mar 1995 354 73 281
## Apr 1995 370 66 304
## May 1995 320 54 266
## Jun 1995 307 41 266
## Jul 1995 421 65 356
## Aug 1995 392 52 340El formato ts es el más antiguo en R (deriva de su homólogo en S). Sólo considera datos regulares, no admite años bisiestos ni datos horarios, etc. Su utilidad en general vendrá limitada al empleo de funciones de R que requieren este formato (notablemente, spectrum, si estamos interesados en análisis espectral). Para la mayor parte de nuestro trabajo preferiremos los formatos zoo y su extensión xts que se describen a continuación.
Formato zoo y xts
Podemos leer datos fechados e interpretar la información sobre fechas directamente, si el formato en que están escritas es más o menos fácil de reconocer.
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericacc.zoo <- read.zoo(file="AccidentesMortales.csv", skip=8, nrow=216, header=FALSE, sep=";",
index=1, format="%YM%m", FUN=as.yearmon)
head(acc.zoo)## V2 V3 V4 V5
## ene 1995 321 58 263 NA
## feb 1995 258 46 212 NA
## mar 1995 354 73 281 NA
## abr 1995 370 66 304 NA
## may 1995 320 54 266 NA
## jun 1995 307 41 266 NASe han fechado bien los datos, pero se ha leído una columna de más (por el ; al fin de cada línea).
Habitualmente preferiremos leer en una operación y fechar en otra. Podríamos emplear read.csv en la lectura, como hemos hecho para el objeto acc, y generar las fechas por ejemplo así:
fecha <- seq.Date(from=as.Date("15-01-1995",format="%d-%m-%Y"),
to=as.Date("15-12-2012","%d-%m-%Y"),
by="month")
head(fecha)## [1] "1995-01-15" "1995-02-15" "1995-03-15" "1995-04-15" "1995-05-15"
## [6] "1995-06-15"(Las especificaciones de formato para pasar de una cadena de caracteres a una fecha pueden consultarse mediante help(strptime).)
Al dato de cada mes le atribuiremos la abscisa temporal que corresponde al día 15. Esto es innecesariamente complejo: sólo tenemos año y mes. Podemos emplear un tipo de fecha que carece de días,
fecha <- as.yearmon(fecha)
head(fecha)## [1] "ene 1995" "feb 1995" "mar 1995" "abr 1995" "may 1995" "jun 1995"y emplearlo para fechar el objeto acc.
acc.zoo <- zoo(x=acc,order.by=fecha)
head(acc.zoo)## MuertosTot MuertosUrb MuertosCtra
## ene 1995 321 58 263
## feb 1995 258 46 212
## mar 1995 354 73 281
## abr 1995 370 66 304
## may 1995 320 54 266
## jun 1995 307 41 266plot(acc.zoo)
Tenemos también cierto control sobre la representación:
plot(acc.zoo, plot.type="single", col=1:3)
Leamos ahora series de consumo de gasoil y gasolina (variables delegadas de la de verdadero interés, que sería “tráfico por carretera”"). Lo hacemos directamente mediante read.zoo:
# gasolina <- read.zoo(file="gasolina.txt",skip=2, header=FALSE) # Error
# gasolina <- read.zoo(file="gasolina.txt",skip=2, header=FALSE,nrow=445) # Error
gasolina <- read.zoo(file="gasolina.txt",skip=2, header=FALSE,nrow=445,
FUN=function(x) { as.yearmon(x,format="%Y:%m") } )
# gasoleo <- read.zoo(file="gasoleo.txt",skip=2, header=FALSE,nrow=445,
# FUN=function(x) { as.yearmon(x,format="%Y:%m") } )
# gasoleo <- read.zoo(file="gasoleo.txt",skip=2, header=FALSE,nrow=445,
# FUN=function(x) { as.yearmon(x,format="%Y:%m") })
gasoleo <- read.zoo(file="gasoleo.txt",skip=2, header=FALSE,nrow=445,
FUN=function(x) { as.yearmon(x,format="%Y:%m") },
fill=TRUE)
plot(gasolina)
plot(gasoleo)
carburante <- cbind(gasolina,gasoleo)
plot(carburante)
Las ventajas del fechado se ponen de manifiesto al reunir en una sóla serie temporal multivariante series univariantes definidas sobre diferentes periodos: el alineamiento es automático y libre de preocupaciones para el usuario:
datos <- cbind(carburante, acc.zoo)
plot(datos)
El formato xts se basa en el excelente diseño de zoo, añade alguna funcionalidad y altera otra. El cambio es transparente, mediante funciones como as.xts o as.zoo:
library(xts)
acc.xts <- as.xts(acc.zoo)
head(acc.xts)## MuertosTot MuertosUrb MuertosCtra
## ene 1995 321 58 263
## feb 1995 258 46 212
## mar 1995 354 73 281
## abr 1995 370 66 304
## may 1995 320 54 266
## jun 1995 307 41 266La función plot de xts representa las series de forma diferente, añadiendo una retícula.
plot(acc.xts)
Frecuentemente puede ser preciso hacer alguna manipulación sobre los datos, especialmente cuando se leen ficheros Excel (o LibreOffice) en que las fechas pueden no reconocerse bien. Lo que sigue es un ejemplo con datos epidemiológicos diarios.
library(readxl)
covid <- read_excel("situacion-epidemiologica.xlsx", skip = 1)## New names:
## * `` -> ...26class(covid)## [1] "tbl_df" "tbl" "data.frame"covid <- as.data.frame(covid)
summary(covid)## Data / Fecha PCR guztira / PCR totales
## Length:252 Min. : 8
## Class :character 1st Qu.: 85098
## Mode :character Median : 287871
## Mean : 380943
## 3rd Qu.: 604697
## Max. :1225477
##
## Test azkarrak guztira / Test rápidos totales
## Min. : 26
## 1st Qu.: 60421
## Median :116544
## Mean : 94979
## 3rd Qu.:131648
## Max. :140946
## NA's :136
## Pertsona bakarrak guztira / Personas únicas totales
## Min. : 63428
## 1st Qu.:148526
## Median :209267
## Mean :224272
## 3rd Qu.:256027
## Max. :657320
## NA's :155
## PCR-dun pertsona bakarrak / Personas únicas con PCR
## Min. :124591
## 1st Qu.:162973
## Median :197287
## Mean :230724
## 3rd Qu.:231194
## Max. :626122
## NA's :178
## Milio bat biztanleko PCR-dun kopurua / Personas con PCR por millón de habitantes
## Min. : 56942
## 1st Qu.: 74484
## Median : 90167
## Mean :105318
## 3rd Qu.:105664
## Max. :284697
## NA's :178
## PCR kasu positiboak / Casos positivos PCR
## Min. :12455
## 1st Qu.:13572
## Median :15943
## Mean :26181
## 3rd Qu.:38884
## Max. :70006
## NA's :61
## Serologia bidez baieztatutako kasuak / Casos positivos detectados por serología
## Min. :5216
## 1st Qu.:6474
## Median :7038
## Mean :6962
## 3rd Qu.:7635
## Max. :7872
## NA's :175
## Kasu positiboak guztira / Casos positivos totales
## Min. : 841
## 1st Qu.:18896
## Median :20320
## Mean :22733
## 3rd Qu.:21227
## Max. :69507
## NA's :155
## Araban PCR kasu positiboak / Casos positivos PCR Álava
## Min. : 2962
## 1st Qu.: 3346
## Median : 3949
## Mean : 5343
## 3rd Qu.: 7490
## Max. :10819
## NA's :81
## Bizkaian PCR kasu positiboak / Casos positivos PCR Bizkaia
## Min. : 7531
## 1st Qu.: 7907
## Median :10903
## Mean :15389
## 3rd Qu.:23266
## Max. :35714
## NA's :81
## Gipuzkoan PCR kasu positiboak / Casos positivos PCR Gipuzkoa
## Min. : 2263
## 1st Qu.: 2314
## Median : 3540
## Mean : 6272
## 3rd Qu.: 9612
## Max. :21627
## NA's :81
## Beste PCR kasu positiboak / Otros casos positivos PCR
## Min. : 156.0
## 1st Qu.: 200.0
## Median : 589.0
## Mean : 735.5
## 3rd Qu.:1181.5
## Max. :1846.0
## NA's :81
## Sendatutakoak edo ospitaleko altak / Recuperados o altas Hospitalarias
## Min. : 0.0
## 1st Qu.: 177.8
## Median : 5620.5
## Mean : 7368.6
## 3rd Qu.:14700.0
## Max. :17895.0
## NA's :150
## Ez berreskuratuak / No recuperados Hildakoak / Fallecidos
## Min. : 0.0 Min. : 0
## 1st Qu.: 617.2 1st Qu.: 477
## Median :1810.5 Median :1442
## Mean :2299.3 Mean :1104
## 3rd Qu.:4115.0 3rd Qu.:1604
## Max. :5206.0 Max. :2153
## NA's :150 NA's :95
## PCR positibodun ospitaleratze berriak plantan / Nuevos ingresos planta con PCR positivo
## Min. : 0.00
## 1st Qu.: 3.00
## Median :14.00
## Mean :20.67
## 3rd Qu.:36.00
## Max. :81.00
## NA's :81
## ZIUn ospitaleratuak / Hospitalizados en UCI R0 Euskadin / R0 en Euskadi
## Min. : 2.00 Length:252
## 1st Qu.: 8.25 Class :character
## Median : 32.50 Mode :character
## Mean : 36.33
## 3rd Qu.: 55.00
## Max. :115.00
## NA's :62
## Negatibizatuak / Negativizados Ospitaleko altak / Altas hospitalarias
## Min. :11902 Min. :4184
## 1st Qu.:13575 1st Qu.:4284
## Median :13912 Median :4339
## Mean :14408 Mean :4742
## 3rd Qu.:14312 3rd Qu.:4468
## Max. :19278 Max. :8068
## NA's :198 NA's :198
## ZIUN ospitaleratze berriak azken 14 egunetan / Nuevos ingresos en UCI en los últimos 14 días
## Min. :1.000
## 1st Qu.:2.500
## Median :3.000
## Mean :3.484
## 3rd Qu.:4.000
## Max. :9.000
## NA's :221
## Positiboen tasa 100.000 biztanleko Araban/Tasa de positivos por 100.000 habitantes en Álava
## Length:252
## Class :character
## Mode :character
##
##
##
##
## Positiboen tasa 100.000 biztanleko Bizkaian/Tasa de positivos por 100.000 habitantes en Bizkaia
## Length:252
## Class :character
## Mode :character
##
##
##
##
## Positiboen tasa 100.000 biztanleko Gizpuzkoan/Tasa de positivos por 100.000 habitantes en Gipuzkoa
## Length:252
## Class :character
## Mode :character
##
##
##
##
## ...26
## Mode:logical
## NA's:252
##
##
##
##
## Las fechas se han leido como cadenas de caracteres. Podemos convertirlas a fechas diariaqs (de tipo Date) así:
covid[,1] <- as.Date(covid[,1], format="%d/%m/%Y")
class(covid[,1])## [1] "Date"Algunas otras variables, que deben ser numéricas, han sido leidas como caracteres, en algunos casos porque emplean coma en lugar de punto decimal. Las convertimos así:
for (i in c(19,23:25)) {
covid[,i] <- gsub(",", ".", covid[,i])
}
for (i in 2:ncol(covid)) {
covid[,i] <- as.numeric(covid[,i])
}Finalmente seleccionamos las variables que nos interesan y les onemos nombres más cortos:
ncols <- c("Fecha","PCRtot","Rapidos","PCRpos","TotPos",
"Recup", "NoRecup", "Muertes","Ingresos","UCI",
"R0","Altas","UCI14")
nloc <- c(1:3,7,9,14,15,16,17:19,21,22)
covid <- covid[,nloc]
colnames(covid) <- ncols
head(covid)## Fecha PCRtot Rapidos PCRpos TotPos Recup NoRecup Muertes Ingresos UCI R0
## 1 2020-02-24 8 NA NA NA 0 0 0 NA NA NA
## 2 2020-02-25 15 NA NA NA 0 0 0 NA NA NA
## 3 2020-02-26 19 NA NA NA 0 0 0 NA NA NA
## 4 2020-02-27 30 NA NA NA 0 0 0 NA NA NA
## 5 2020-02-28 60 NA NA NA 0 0 0 NA NA NA
## 6 2020-02-29 79 NA NA NA 0 0 0 NA NA NA
## Altas UCI14
## 1 NA NA
## 2 NA NA
## 3 NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NAAhora podemos crear una serie multivariante diaria así,
covid <- zoo(x=covid[,-1], order.by=covid[,1])y dibujar cualquiera de sus componentes o todas ellas:
plot(covid[,"UCI"])
Manipulaciones simples
Hay manipulaciones simples muy frecuentes que podemos necesitar antes de comenzar el tratamiento de una serie temporal. Muchas de ellas las podremos hacer tanto en formato zoo como en formato xts, otras pueden requerir un formato específico. Interesa saber que xts es un objeto que “hereda” de zoo. Podemos sustituir un objeto xts encualquier lugar en que uno zoo sea requerido.
Por ejemplo, podemos querer obtener los totales trimestrales de muertos en los datos sobre accidentes de tráfico. El resultado es una serie de periodicidad trimestral, fechada el último mes de cada trimestre:
acc.TotTrim <- apply.quarterly(acc.xts$MuertosTot, FUN=sum)
class(acc.TotTrim)## [1] "xts" "zoo"acc.TotTrim## MuertosTot
## mar 1995 933
## jun 1995 997
## sep 1995 1209
## dic 1995 995
## mar 1996 852
## jun 1996 1003
## sep 1996 1096
## dic 1996 1011
## mar 1997 902
## jun 1997 928
## sep 1997 1207
## dic 1997 1025
## mar 1998 944
## jun 1998 1010
## sep 1998 1217
## dic 1998 1148
## mar 1999 927
## jun 1999 973
## sep 1999 1209
## dic 1999 1130
## mar 2000 993
## jun 2000 1037
## sep 2000 1219
## dic 2000 1123
## mar 2001 949
## jun 2001 1045
## sep 2001 1156
## dic 2001 1020
## mar 2002 968
## jun 2002 931
## sep 2002 1135
## dic 2002 997
## mar 2003 918
## jun 2003 965
## sep 2003 1172
## dic 2003 1029
## mar 2004 820
## jun 2004 932
## sep 2004 983
## dic 2004 908
## mar 2005 816
## jun 2005 842
## sep 2005 917
## dic 2005 802
## mar 2006 769
## jun 2006 823
## sep 2006 803
## dic 2006 724
## mar 2007 706
## jun 2007 732
## sep 2007 831
## dic 2007 688
## mar 2008 572
## jun 2008 583
## sep 2008 644
## dic 2008 586
## mar 2009 516
## jun 2009 531
## sep 2009 568
## dic 2009 502
## mar 2010 415
## jun 2010 444
## sep 2010 576
## dic 2010 518
## mar 2011 364
## jun 2011 398
## sep 2011 509
## dic 2011 412
## mar 2012 374
## jun 2012 372
## sep 2012 435
## dic 2012 342acc.TotTrim <- zoo(x=coredata(acc.TotTrim), order.by=as.yearqtr(index(acc.TotTrim)))
plot(acc.TotTrim)
Del mismo modo podríamos calcular el promedio mensual (FUN=mean) o cualquier otra cosa, pues podemos incluso definir la función que pasamos como argumento. Observemos que aggregate es una función que trabaja sobre objetos zoo, pero como xts es una especialización de zoo, podríamos invocarla también sobre objetos de clase xts. Una forma alternativa de hacer lo anterior sería:
acc.TotTrim.alt <- aggregate(acc.xts$MuertosTot, as.yearqtr, sum)
plot(acc.TotTrim.alt)
Un último tipo de manipulación que examinaremos se parece al precedente: son agregaciones en periodos móviles. Por ejemplo, la suma de los tres últimos meses, o la media móvil trimestral centrada. La función rollapply permite calcular fácilmente una media móvil como la requerida:
acc.media3M <- rollapply(acc.zoo$MuertosTot, width=3, FUN=mean, align="center")
head(acc.media3M)## feb 1995 mar 1995 abr 1995 may 1995 jun 1995 jul 1995
## 311.0000 327.3333 348.0000 332.3333 349.3333 373.3333acc.Tot3M <- rollapply(acc.zoo$MuertosTot, width=3, FUN=sum, align="right")
head(acc.Tot3M)## mar 1995 abr 1995 may 1995 jun 1995 jul 1995 ago 1995
## 933 982 1044 997 1048 1120Al igual que antes, la posibilidad de definir la función que deseamos aplicar ensancha mucho el ámbito de utilización. Por ejemplo, sin en el caso anterior en lugar de la media aritmética hubiéramos querido por algún motivo la media armónica, hubiéramos podido hacer:
acc.3Msumlog <- rollapply(log(acc.zoo$MuertosTot), width=3, FUN=sum, align="center")
acc.Tot3M.armon <- exp( (1/3)*acc.3Msumlog )
head(acc.Tot3M.armon)## feb 1995 mar 1995 abr 1995 may 1995 jun 1995 jul 1995
## 308.3491 323.3019 347.3641 331.2557 345.8253 370.0288Algunas manipulaciones simples como la toma de diferencias pueden hacerse directamente, tanto en formato zoo como xts:
dPCRtot <- diff(covid[,"PCRtot"])
plot(dPCRtot, main="Pruebas PCR diarias en la CAPV")
Otros paquetes especializados
Los formatos anteriores, con ser flexibles y potentes, de ningún modo agotan las facilidades disponibles en R. Algunos otros paquetes que pueden ser de interés para utilizaciones especiales incluyen los siguientes.
lubridate
Es un excelente paquete que incluye muchas funciones de utilidad para hacer la vida más fácil con datos temporales. Añade a los tipos de fechado más habituales en R (Date, POSIX, etc.) la posibilidad de usar resoluciones mucho mas pequeñas: series cuya abscisa temporal se cuenta en microsegundos, por ejemplo. Define también duraciones. Conviene leer su viñeta.
timeDate
Pertenece al conjunto de paquetes Rmetrics y es de utilidad para manejo de series financieras. (Sobre el empleo de R con series financieras puede verse Working with Financial Time Series Data in R.) Facilita la vida cuando hay que tener presentes festividades en distintos países, difedrentes zonas horarias o fiestas móviles.
Consideremos por ejemplo la serie siguiente de consumo de carne de cordero en un gran supermercado en San Sebastián/Donostia:
cordero <- dget(file="/home/etptupaf/ft/docencia/MasterModMat/MaterialClase/Datos/Cordero.dge")
class(cordero)## [1] "ts"plot(cordero)
Hay que saber que el formato ts especifica las abscisas como número fraccionarios; o sea, que el segundo mes del año 2002, por ejemplo, se representa en la abscisa 2002 + 2/12 = 2002.167, y así los demás.
Observemos que en cada año aparecen tres picos, el más notorio claramente asociado al consumo navideño, otro en mitad de agosto asociado a festividades locales y un tercero que podemos sospechar asociado al consumo en la temporada de Pascua.
Podemos recurrir al paquete timeDate para darnos las fechas de Pascua de los años concernidos:
library(timeDate)
Pascua <- Easter(2001:2008)
Pascua## GMT
## [1] [2001-04-15] [2002-03-31] [2003-04-20] [2004-04-11] [2005-03-27]
## [6] [2006-04-16] [2007-04-08] [2008-03-23]class(Pascua)## [1] "timeDate"
## attr(,"package")
## [1] "timeDate"Para poder representar esas abscisas temporales en el plot anterior, necesitamos convertirlas a años decimales. El paquete lubridate nos proporciona una función para ello: decimal_date. No admite como input objetos de tipo timeDate, hemos de convertir a Date.
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionmarcas <- decimal_date(as.Date(Pascua))Ahora podemos señalar las fectividades de Pascua:
plot(cordero)
abline(v=marcas, col="red")
Nuestra sospecha queda corroborada.
anytime
El propósito de anytime es convertir casi cualquier cosa que tenga apariencia de representación textual de una fecha a un formato de fecha utilizado por R; y viceversa. Tiene una viñeta bastante informativa.
Algunos ejemplos:
library(anytime)
anydate(20160101 + 0:2)## [1] "2016-01-01" "2016-01-02" "2016-01-03"a <- anydate(20160101 + 0:2)
a## [1] "2016-01-01" "2016-01-02" "2016-01-03"class(a)## [1] "Date"a <- anytime("20160203 10:23")
a## [1] "2016-02-03 10:23:00 CET"class(a)## [1] "POSIXct" "POSIXt"Cambios de zona horaria
a## [1] "2016-02-03 10:23:00 CET"attr(a, "tzone") <- "UTC"
a## [1] "2016-02-03 09:23:00 UTC"b <- a <- anytime("2018100100") + (1:(30*24))*3600
attr(b, "tzone") <- "UTC"
head(a)## [1] "2018-10-01 01:00:00 CEST" "2018-10-01 02:00:00 CEST"
## [3] "2018-10-01 03:00:00 CEST" "2018-10-01 04:00:00 CEST"
## [5] "2018-10-01 05:00:00 CEST" "2018-10-01 06:00:00 CEST"class(a)## [1] "POSIXct" "POSIXt"Observemos que el cambio de hora se toma en cuenta automáticamente:
ab <- data.frame(CET=a, UTC=b)
ab[640:660,]## CET UTC
## 640 2018-10-27 16:00:00 2018-10-27 14:00:00
## 641 2018-10-27 17:00:00 2018-10-27 15:00:00
## 642 2018-10-27 18:00:00 2018-10-27 16:00:00
## 643 2018-10-27 19:00:00 2018-10-27 17:00:00
## 644 2018-10-27 20:00:00 2018-10-27 18:00:00
## 645 2018-10-27 21:00:00 2018-10-27 19:00:00
## 646 2018-10-27 22:00:00 2018-10-27 20:00:00
## 647 2018-10-27 23:00:00 2018-10-27 21:00:00
## 648 2018-10-28 00:00:00 2018-10-27 22:00:00
## 649 2018-10-28 01:00:00 2018-10-27 23:00:00
## 650 2018-10-28 02:00:00 2018-10-28 00:00:00
## 651 2018-10-28 02:00:00 2018-10-28 01:00:00
## 652 2018-10-28 03:00:00 2018-10-28 02:00:00
## 653 2018-10-28 04:00:00 2018-10-28 03:00:00
## 654 2018-10-28 05:00:00 2018-10-28 04:00:00
## 655 2018-10-28 06:00:00 2018-10-28 05:00:00
## 656 2018-10-28 07:00:00 2018-10-28 06:00:00
## 657 2018-10-28 08:00:00 2018-10-28 07:00:00
## 658 2018-10-28 09:00:00 2018-10-28 08:00:00
## 659 2018-10-28 10:00:00 2018-10-28 09:00:00
## 660 2018-10-28 11:00:00 2018-10-28 10:00:00