R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("pastecs-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('pastecs') Loading required package: boot > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "AutoD2" > > ### * AutoD2 > > flush(stderr()); flush(stdout()) > > ### Name: AutoD2 > ### Title: AutoD2, CrossD2 or CenterD2 analysis of a multiple time-series > ### Aliases: AutoD2 CrossD2 CenterD2 > ### Keywords: ts multivariate htest > > ### ** Examples > > data(marphy) > marphy.ts <- as.ts(as.matrix(marphy[, 1:3])) > AutoD2(marphy.ts) $lag [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 $D2 [1] 4.22573 10.58976 23.47267 43.38305 153.44283 433.54182 [7] 1064.54132 1802.43083 1940.63240 2111.62737 2298.56609 2166.24216 [13] 2443.59445 2708.61304 3121.54365 3267.95116 2370.99649 1772.94559 [19] 1366.08938 1239.05153 1045.82545 163.43575 $call AutoD2(series = marphy.ts) $data [1] "marphy.ts" $type [1] "AutoD2" $units.text [1] "" attr(,"class") [1] "D2" > marphy.ts2 <- as.ts(as.matrix(marphy[, c(1, 4, 3)])) > CrossD2(marphy.ts, marphy.ts2) $lag [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 $D2 [1] 4.323043 10.155231 45.764444 124.524831 395.341890 [6] 1008.216533 2314.110604 3780.247450 4391.117099 5289.482066 [11] 6527.234759 7341.803004 9887.267753 13544.128938 18768.866938 [16] 24708.252337 23961.156382 23500.686505 23796.894859 25932.921344 [21] 28394.452862 22836.606581 $call CrossD2(series = marphy.ts, series2 = marphy.ts2) $data [1] "marphy.ts" $data2 [1] "marphy.ts2" $type [1] "CrossD2" $units.text [1] "" attr(,"class") [1] "D2" > # This is not identical to: > CrossD2(marphy.ts2, marphy.ts) $lag [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 $D2 [1] 30.76921 62.83820 112.50639 215.13122 475.90799 979.41968 [7] 1907.94064 2927.17901 3434.67167 4209.00250 4918.96322 5091.38349 [13] 5842.69358 6713.53204 7942.49390 8759.94122 7834.21656 7034.27185 [19] 6140.98599 5689.15330 5334.47184 3841.11059 $call CrossD2(series = marphy.ts2, series2 = marphy.ts) $data [1] "marphy.ts2" $data2 [1] "marphy.ts" $type [1] "CrossD2" $units.text [1] "" attr(,"class") [1] "D2" > marphy.d2 <- CenterD2(marphy.ts, window=16) > lines(c(17, 17), c(-1, 15), col=4, lty=2) > lines(c(25, 25), c(-1, 15), col=4, lty=2) > lines(c(30, 30), c(-1, 15), col=4, lty=2) > lines(c(41, 41), c(-1, 15), col=4, lty=2) > lines(c(46, 46), c(-1, 15), col=4, lty=2) > text(c(8.5, 21, 27.5, 35, 43.5, 57), 11, labels=c("Peripheral Zone", "D1", + "C", "Front", "D2", "Central Zone")) # Labels > time(marphy.ts)[marphy.d2$D2 > marphy.d2$chisq] [1] 17 26 32 42 43 46 60 > > > > cleanEx(); ..nameEx <- "GetUnitText" > > ### * GetUnitText > > flush(stderr()); flush(stdout()) > > ### Name: GetUnitText > ### Title: Format a nice time units for labels in graphs > ### Aliases: GetUnitText > ### Keywords: ts > > ### ** Examples > > timeser <- ts(1:24, frequency=12) # 12 observations per year > if (exists("is.R") && is.function(is.R) && is.R()) { # We are in R + attr(timeser, "units") <- "years" # time in years for 'ts' object + } else { # We are in Splus + attr(attr(timeser, "tspar"), "units") <- "years" # Idem for Splus 'rts' object + } > GetUnitText(timeser) # formats unit (1/12 year=months) [1] "months" > > > > cleanEx(); ..nameEx <- "abund" > > ### * abund > > flush(stderr()); flush(stdout()) > > ### Name: abund > ### Title: Sort variables by abundance > ### Aliases: abund extract.abund identify.abund lines.abund plot.abund > ### print.abund print.summary.abund summary.abund > ### Keywords: multivariate > > ### ** Examples > > data(bnr) > bnr.abd <- abund(bnr) > summary(bnr.abd) Sorting of descriptors according to abundance for: bnr Coefficient f: 0.2 163 variables sorted Number of individuals (% of most abundant in log): S8 S2 S3 S4 S6 S13 S1 72.273641 90.069317 88.739428 78.019800 76.302796 68.235798 100.000000 S5 S10 S14 S9 S15 S21 S17 77.095883 70.713663 67.746617 71.405852 65.019607 56.765809 61.013010 S22 S39 S12 S26 S25 S11 S38 56.217342 42.731127 68.629878 53.892117 55.223940 68.848694 43.220902 S19 S20 S29 S37 S41 S27 S16 60.065505 56.854955 50.145792 43.416565 42.561076 51.790497 62.798247 S45 S23 S43 S33 S36 S24 S47 41.106427 55.899523 42.254638 47.328915 44.929107 55.322688 39.646573 S49 S52 S31 S32 S46 S67 S58 38.667117 36.354344 49.292179 48.648509 39.888968 27.616484 32.237324 S50 S71 S54 S35 S61 S18 S34 36.889898 25.804639 33.768615 45.427018 31.764202 60.458819 47.328915 S68 S7 S69 S64 S72 S62 S48 27.056287 74.463292 26.760995 28.868131 25.459071 30.908713 39.396718 S70 S59 S78 S51 S42 S80 S53 26.454557 32.082659 23.004660 36.803057 42.431159 20.176537 35.976618 S73 S88 S100 S84 S77 S55 S56 25.459071 17.840663 14.548434 19.470036 23.468766 32.538086 32.538086 S30 S112 S135 S136 S137 S138 S139 49.694585 11.256205 0.000000 0.000000 0.000000 0.000000 0.000000 S140 S141 S142 S143 S144 S145 S146 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 S147 S148 S149 S150 S151 S152 S153 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 S154 S155 S156 S157 S158 S159 S160 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 S161 S162 S163 S92 S74 S114 S90 0.000000 0.000000 0.000000 15.800081 24.324255 8.920332 16.884308 S91 S107 S119 S120 S121 S125 S127 16.884308 13.068049 5.628103 5.628103 5.628103 5.628103 5.628103 S128 S131 S132 S133 S134 S93 S96 5.628103 5.628103 5.628103 5.628103 5.628103 14.548434 14.548434 S98 S85 S83 S113 S115 S117 S102 14.548434 18.696152 19.470036 8.920332 8.920332 8.920332 13.068049 S103 S108 S122 S123 S124 S126 S129 13.068049 13.068049 5.628103 5.628103 5.628103 5.628103 5.628103 S130 S89 S94 S95 S97 S101 S110 5.628103 17.840663 14.548434 14.548434 14.548434 14.548434 11.256205 S111 S81 S116 S118 S104 S105 S106 11.256205 20.176537 8.920332 8.920332 13.068049 13.068049 13.068049 S63 S99 S87 S109 S60 S75 S82 29.096868 14.548434 18.696152 11.256205 31.764202 24.324255 19.470036 S76 S79 S86 S57 S65 S66 S44 23.468766 20.826454 18.696152 32.538086 28.632763 28.390368 41.700812 S28 S40 50.964058 42.688947 Percent of non-zero values: S8 S2 S3 S4 S6 S13 S1 95.1456311 95.1456311 93.2038835 89.3203883 86.4077670 82.5242718 89.3203883 S5 S10 S14 S9 S15 S21 S17 83.4951456 81.5533981 79.6116505 78.6407767 76.6990291 70.8737864 70.8737864 S22 S39 S12 S26 S25 S11 S38 68.9320388 64.0776699 67.9611650 64.0776699 64.0776699 66.0194175 59.2233010 S19 S20 S29 S37 S41 S27 S16 63.1067961 58.2524272 54.3689320 51.4563107 49.5145631 50.4854369 51.4563107 S45 S23 S43 S33 S36 S24 S47 43.6893204 44.6601942 39.8058252 38.8349515 35.9223301 37.8640777 33.0097087 S49 S52 S31 S32 S46 S67 S58 30.0970874 28.1553398 31.0679612 30.0970874 23.3009709 19.4174757 20.3883495 S50 S71 S54 S35 S61 S18 S34 21.3592233 18.4466019 20.3883495 23.3009709 19.4174757 26.2135922 22.3300971 S68 S7 S69 S64 S72 S62 S48 16.5048544 28.1553398 15.5339806 15.5339806 14.5631068 15.5339806 16.5048544 S70 S59 S78 S51 S42 S80 S53 12.6213592 12.6213592 9.7087379 12.6213592 13.5922330 7.7669903 11.6504854 S73 S88 S100 S84 S77 S55 S56 8.7378641 6.7961165 5.8252427 6.7961165 7.7669903 9.7087379 9.7087379 S30 S112 S135 S136 S137 S138 S139 13.5922330 3.8834951 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 S140 S141 S142 S143 S144 S145 S146 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 S147 S148 S149 S150 S151 S152 S153 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 S154 S155 S156 S157 S158 S159 S160 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 S161 S162 S163 S92 S74 S114 S90 0.9708738 0.9708738 0.9708738 4.8543689 6.7961165 2.9126214 4.8543689 S91 S107 S119 S120 S121 S125 S127 4.8543689 3.8834951 1.9417476 1.9417476 1.9417476 1.9417476 1.9417476 S128 S131 S132 S133 S134 S93 S96 1.9417476 1.9417476 1.9417476 1.9417476 1.9417476 3.8834951 3.8834951 S98 S85 S83 S113 S115 S117 S102 3.8834951 4.8543689 4.8543689 1.9417476 1.9417476 1.9417476 2.9126214 S103 S108 S122 S123 S124 S126 S129 2.9126214 2.9126214 0.9708738 0.9708738 0.9708738 0.9708738 0.9708738 S130 S89 S94 S95 S97 S101 S110 0.9708738 3.8834951 2.9126214 2.9126214 2.9126214 2.9126214 1.9417476 S111 S81 S116 S118 S104 S105 S106 1.9417476 3.8834951 0.9708738 0.9708738 1.9417476 1.9417476 1.9417476 S63 S99 S87 S109 S60 S75 S82 5.8252427 1.9417476 2.9126214 0.9708738 5.8252427 3.8834951 1.9417476 S76 S79 S86 S57 S65 S66 S44 2.9126214 1.9417476 0.9708738 3.8834951 1.9417476 0.9708738 3.8834951 S28 S40 4.8543689 0.9708738 > plot(bnr.abd, dpos=c(105, 100)) > bnr.abd$n <- 26 > # To identify a point on the graph, use: bnr.abd$n <- identify(bnr.abd) > lines(bnr.abd) > bnr2 <- extract(bnr.abd) Warning in UseMethod("extract", e, n, ...) : arguments after the first two are ignored > names(bnr2) [1] "S8" "S2" "S3" "S4" "S6" "S13" "S1" "S5" "S10" "S14" "S9" "S15" [13] "S21" "S17" "S22" "S39" "S12" "S26" "S25" "S11" "S38" "S19" "S20" "S29" [25] "S37" "S41" > > > > cleanEx(); ..nameEx <- "buysbal" > > ### * buysbal > > flush(stderr()); flush(stdout()) > > ### Name: buysbal > ### Title: Buys-Ballot table > ### Aliases: buysbal > ### Keywords: ts > > ### ** Examples > > data(releve) > buysbal(releve$Day, releve$Melosul, frequency=4, units="days", + datemin="21/03/1989", dateformat="d/m/Y") 1 2 3 4 1989 10500.000 14000.000 4666.667 23530 1990 5406.667 2550.000 5950.000 41920 1991 8800.000 21166.667 5668.333 9900 1992 7566.667 6816.667 4000.000 364 > buysbal(releve$Day, releve$Melosul, frequency=4, units="days", + datemin="21/03/1989", dateformat="d/m/Y", count=TRUE) 1 2 3 4 1989 1 1 3 2 1990 3 2 6 5 1991 3 3 6 6 1992 6 6 6 2 > > > > cleanEx(); ..nameEx <- "daystoyears" > > ### * daystoyears > > flush(stderr()); flush(stdout()) > > ### Name: daystoyears > ### Title: Convert time units from "days" to "years" or back > ### Aliases: daystoyears yearstodays > ### Keywords: ts > > ### ** Examples > > # A vector with a "days" time-scale (25 values every 30 days) > A <- (1:25)*30 > # Convert it to a "years" time-scale, using 23/05/2001 (d/m/Y) as first value > B <- daystoyears(A, datemin="23/05/2001", dateformat="d/m/Y") > B [1] 2001.389 2001.472 2001.554 2001.636 2001.718 2001.800 2001.882 2001.964 [9] 2002.047 2002.129 2002.211 2002.293 2002.375 2002.457 2002.539 2002.621 [17] 2002.704 2002.786 2002.868 2002.950 2003.032 2003.114 2003.196 2003.279 [25] 2003.361 > # Convert it back to "days" time-scale > yearstodays(B, xmin=30) [1] 30 60 90 120 150 180 210 240 270 300 330 360 390 420 450 480 510 540 570 [20] 600 630 660 690 720 750 > > # Here is an example showing the shift introduced, and its consequence: > C <- daystoyears(unclass(as.Date(c("1970-1-1","1971-1-1","1972-1-1","1973-1-1")))) > C [1] 1970.000 1970.999 1971.999 1973.001 > > > > cleanEx(); ..nameEx <- "decaverage" > > ### * decaverage > > flush(stderr()); flush(stdout()) > > ### Name: decaverage > ### Title: Time series decomposition using a moving average > ### Aliases: decaverage > ### Keywords: ts smooth > > ### ** Examples > > data(marbio) > ClausoB.ts <- ts(log(marbio$ClausocalanusB + 1)) > ClausoB.dec <- decaverage(ClausoB.ts, order=2, times=10, sides=2, ends="fill") > plot(ClausoB.dec, col=c(1, 3, 2), xlab="stations") > # A stacked graph is more representative in this case > plot(ClausoB.dec, col=c(1, 3), xlab="stations", stack=FALSE, resid=FALSE, + lpos=c(53, 4.3)) > > > > cleanEx(); ..nameEx <- "deccensus" > > ### * deccensus > > flush(stderr()); flush(stdout()) > > ### Name: deccensus > ### Title: Time decomposition using the CENSUS II method > ### Aliases: deccensus > ### Keywords: ts smooth > > ### ** Examples > > data(releve) > # Get regulated time series with a 'years' time-scale > rel.regy <- regul(releve$Day, releve[3:8], xmin=6, n=87, units="daystoyears", frequency=24, tol=2.2, methods="linear", datemin="21/03/1989", dateformat="d/m/Y") A 'tol' of 2.2 in 'days' is 0.00602327173169062 in 'years' 'tol' was adjusted to 0.00595238095238083 > rel.ts <- tseries(rel.regy) > # We must have complete cycles to allow using deccensus() > start(rel.ts) [1] 1989 6 > end(rel.ts) [1] 1992 20 > rel.ts2 <- window(rel.ts, end=c(1992,5)) > rel.dec2 <- deccensus(rel.ts2[, "Melosul"], trend=TRUE) > plot(rel.dec2, col=c(1,4,3,2)) > > > > cleanEx(); ..nameEx <- "decdiff" > > ### * decdiff > > flush(stderr()); flush(stdout()) > > ### Name: decdiff > ### Title: Time series decomposition using differences (trend elimination) > ### Aliases: decdiff > ### Keywords: ts smooth > > ### ** Examples > > data(marbio) > ClausoB.ts <- ts(log(marbio$ClausocalanusB + 1)) > ClausoB.dec <- decdiff(ClausoB.ts, lag=1, order=2, ends="fill") > plot(ClausoB.dec, col=c(1, 4, 2), xlab="stations") > > > > cleanEx(); ..nameEx <- "decevf" > > ### * decevf > > flush(stderr()); flush(stdout()) > > ### Name: decevf > ### Title: Time series decomposition using eigenvector filtering (EVF) > ### Aliases: decevf > ### Keywords: ts smooth > > ### ** Examples > > data(releve) > melo.regy <- regul(releve$Day, releve$Melosul, xmin=9, n=87, + units="daystoyears", frequency=24, tol=2.2, methods="linear", + datemin="21/03/1989", dateformat="d/m/Y") A 'tol' of 2.2 in 'days' is 0.00602327173169062 in 'years' 'tol' was adjusted to 0.00595238095238083 > melo.ts <- tseries(melo.regy) > acf(melo.ts) > # Autocorrelation is not significant after 0.16 > # That corresponds to a lag of 0.16*24=4 (frequency=24) > melo.evf <- decevf(melo.ts, lag=4, axes=1) > plot(melo.evf, col=c(1, 4, 2)) > # A superposed graph is better in the present case > plot(melo.evf, col=c(1, 4), xlab="stations", stack=FALSE, resid=FALSE, + lpos=c(0, 60000)) > > > > cleanEx(); ..nameEx <- "decloess" > > ### * decloess > > flush(stderr()); flush(stdout()) > > ### Name: decloess > ### Title: Time series decomposition by the LOESS method > ### Aliases: decloess > ### Keywords: ts smooth > > ### ** Examples > > data(releve) > melo.regy <- regul(releve$Day, releve$Melosul, xmin=9, n=87, + units="daystoyears", frequency=24, tol=2.2, methods="linear", + datemin="21/03/1989", dateformat="d/m/Y") A 'tol' of 2.2 in 'days' is 0.00602327173169062 in 'years' 'tol' was adjusted to 0.00595238095238083 > melo.ts <- tseries(melo.regy) > melo.dec <- decloess(melo.ts, s.window="periodic") > plot(melo.dec, col=1:3) > > > > cleanEx(); ..nameEx <- "decmedian" > > ### * decmedian > > flush(stderr()); flush(stdout()) > > ### Name: decmedian > ### Title: Time series decomposition using a running median > ### Aliases: decmedian > ### Keywords: ts smooth > > ### ** Examples > > data(marbio) > ClausoB.ts <- ts(log(marbio$ClausocalanusB + 1)) > ClausoB.dec <- decmedian(ClausoB.ts, order=2, times=10, ends="fill") > plot(ClausoB.dec, col=c(1, 4, 2), xlab="stations") > # This is a transect across a frontal zone: > plot(ClausoB.dec, col=c(0, 2), xlab="stations", stack=FALSE, resid=FALSE) > lines(c(17, 17), c(0, 10), col=4, lty=2) > lines(c(25, 25), c(0, 10), col=4, lty=2) > lines(c(30, 30), c(0, 10), col=4, lty=2) > lines(c(41, 41), c(0, 10), col=4, lty=2) > lines(c(46, 46), c(0, 10), col=4, lty=2) > text(c(8.5, 21, 27.5, 35, 43.5, 57), 8.7, labels=c("Peripheral Zone", "D1", + "C", "Front", "D2", "Central Zone")) > > > > cleanEx(); ..nameEx <- "decreg" > > ### * decreg > > flush(stderr()); flush(stdout()) > > ### Name: decreg > ### Title: Time series decomposition using a regression model > ### Aliases: decreg > ### Keywords: ts smooth > > ### ** Examples > > data(marphy) > density <- ts(marphy[, "Density"]) > plot(density) > Time <- time(density) > > # Linear model to represent trend > density.lin <- lm(density ~ Time) > summary(density.lin) Call: lm(formula = density ~ Time) Residuals: Min 1Q Median 3Q Max -0.052226 -0.018014 0.001945 0.017672 0.058895 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.884e+01 6.585e-03 4379.27 <2e-16 *** Time 3.605e-03 1.659e-04 21.73 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.02685 on 66 degrees of freedom Multiple R-Squared: 0.8774, Adjusted R-squared: 0.8755 F-statistic: 472.3 on 1 and 66 DF, p-value: < 2.2e-16 > xreg <- predict(density.lin) > lines(xreg, col=3) > density.dec <- decreg(density, xreg) > plot(density.dec, col=c(1, 3, 2), xlab="stations") > > # Order 2 polynomial to represent trend > density.poly <- lm(density ~ Time + I(Time^2)) > summary(density.poly) Call: lm(formula = density ~ Time + I(Time^2)) Residuals: Min 1Q Median 3Q Max -0.067079 -0.010665 0.001498 0.014777 0.045341 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.880e+01 8.374e-03 3439.404 < 2e-16 *** Time 6.593e-03 5.600e-04 11.773 < 2e-16 *** I(Time^2) -4.330e-05 7.866e-06 -5.505 6.73e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.02234 on 65 degrees of freedom Multiple R-Squared: 0.9164, Adjusted R-squared: 0.9138 F-statistic: 356.2 on 2 and 65 DF, p-value: < 2.2e-16 > xreg2 <- predict(density.poly) > plot(density) > lines(xreg2, col=3) > density.dec2 <- decreg(density, xreg2) > plot(density.dec2, col=c(1, 3, 2), xlab="stations") > > # Fit a sinusoidal model on seasonal (artificial) data > tser <- ts(sin((1:100)/12*pi)+rnorm(100, sd=0.3), start=c(1998, 4), + frequency=24) > Time <- time(tser) > tser.sin <- lm(tser ~ I(cos(2*pi*Time)) + I(sin(2*pi*Time))) > summary(tser.sin) Call: lm(formula = tser ~ I(cos(2 * pi * Time)) + I(sin(2 * pi * Time))) Residuals: Min 1Q Median 3Q Max -0.73366 -0.16820 -0.01390 0.16019 0.66442 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.03423 0.02694 1.271 0.207 I(cos(2 * pi * Time)) -0.47737 0.03852 -12.392 <2e-16 *** I(sin(2 * pi * Time)) 0.81252 0.03766 21.576 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.269 on 97 degrees of freedom Multiple R-Squared: 0.8623, Adjusted R-squared: 0.8594 F-statistic: 303.6 on 2 and 97 DF, p-value: < 2.2e-16 > tser.reg <- predict(tser.sin) > tser.dec <- decreg(tser, tser.reg) > plot(tser.dec, col=c(1, 4), xlab="stations", stack=FALSE, resid=FALSE, + lpos=c(0, 4)) > plot(tser.dec, col=c(1, 4, 2), xlab="stations") > > # One can also use nonlinear models (see 'nls') > # or autoregressive models (see 'ar' and others in 'ts' library) > > > > cleanEx(); ..nameEx <- "disjoin" > > ### * disjoin > > flush(stderr()); flush(stdout()) > > ### Name: disjoin > ### Title: Complete disjoined coded data (binary coding) > ### Aliases: disjoin > ### Keywords: manip > > ### ** Examples > > # Artificial data with 1/5 of zeros > Z <- c(abs(rnorm(8000)), rep(0, 2000)) > # Let the program chose cuts > table(cut(Z, breaks=5)) (-0.00381,0.76] (0.76,1.52] (1.52,2.29] (2.29,3.05] (3.05,3.81] 6361 2560 886 176 17 > # Create one class for zeros, and 4 classes for the other observations > Z2 <- Z[Z != 0] > cuts <- c(-1e-10, 1e-10, quantile(Z2, 1:5/5, na.rm=TRUE)) > cuts 20% 40% 60% -0.0000000001 0.0000000001 0.2542817121 0.5275370301 0.8636987590 80% 100% 1.3156956655 3.8102766807 > table(cut(Z, breaks=cuts)) (-1e-10,1e-10] (1e-10,0.254] (0.254,0.528] (0.528,0.864] (0.864,1.32] 2000 1600 1600 1600 1600 (1.32,3.81] 1600 > # Binary coding of these data > disjoin(cut(Z, breaks=cuts))[1:10, ] (-1e-10,1e-10] (1e-10,0.254] (0.254,0.528] (0.528,0.864] (0.864,1.32] [1,] 0 0 0 1 0 [2,] 0 1 0 0 0 [3,] 0 0 0 1 0 [4,] 0 0 0 0 0 [5,] 0 0 1 0 0 [6,] 0 0 0 1 0 [7,] 0 0 1 0 0 [8,] 0 0 0 1 0 [9,] 0 0 0 1 0 [10,] 0 0 1 0 0 (1.32,3.81] [1,] 0 [2,] 0 [3,] 0 [4,] 1 [5,] 0 [6,] 0 [7,] 0 [8,] 0 [9,] 0 [10,] 0 > > > > cleanEx(); ..nameEx <- "disto" > > ### * disto > > flush(stderr()); flush(stdout()) > > ### Name: disto > ### Title: Compute and plot a distogram > ### Aliases: disto > ### Keywords: multivariate ts > > ### ** Examples > > data(bnr) > disto(bnr) distance distogram 1 1 0.3581454 2 2 0.4797966 3 3 0.5173366 4 4 0.5501900 5 5 0.4942032 6 6 0.4922207 7 7 0.5422460 8 8 0.5774747 9 9 0.6148662 10 10 0.6674661 11 11 0.6634647 12 12 0.6457200 13 13 0.6567708 14 14 0.6706163 15 15 0.6743301 16 16 0.6943444 17 17 0.7143640 18 18 0.7279900 19 19 0.7162906 20 20 0.6569722 21 21 0.6708284 22 22 0.6974970 23 23 0.6976051 24 24 0.7008844 25 25 0.6854326 26 26 0.6526971 > > > > cleanEx(); ..nameEx <- "escouf" > > ### * escouf > > flush(stderr()); flush(stdout()) > > ### Name: escouf > ### Title: Choose variables using the Escoufier's equivalent vectors method > ### Aliases: escouf extract.escouf identify.escouf lines.escouf plot.escouf > ### print.escouf print.summary.escouf summary.escouf > ### Keywords: multivariate > > ### ** Examples > > data(marbio) > marbio.esc <- escouf(marbio) Variable 5 incorporated, RV = 0.5938735 Variable 9 incorporated, RV = 0.7607149 Variable 3 incorporated, RV = 0.8280316 Variable 17 incorporated, RV = 0.8620636 Variable 11 incorporated, RV = 0.8981089 Variable 13 incorporated, RV = 0.9039953 Variable 10 incorporated, RV = 0.9109498 Variable 7 incorporated, RV = 0.917168 Variable 23 incorporated, RV = 0.9254146 Variable 1 incorporated, RV = 0.932456 Variable 19 incorporated, RV = 0.9399177 Variable 4 incorporated, RV = 0.9469293 Variable 21 incorporated, RV = 0.952844 Variable 20 incorporated, RV = 0.9584662 Variable 6 incorporated, RV = 0.9621404 Variable 14 incorporated, RV = 0.9679197 Variable 22 incorporated, RV = 0.9727763 Variable 24 incorporated, RV = 0.9773901 Variable 15 incorporated, RV = 0.981146 Variable 18 incorporated, RV = 0.9850263 Variable 8 incorporated, RV = 0.9882042 Variable 16 incorporated, RV = 0.9921904 Variable 2 incorporated, RV = 0.995268 Variable 12 incorporated, RV = 1 > summary(marbio.esc) Escoufier's method of equivalent vectors for: marbio Calculation level: 1 24 variable(s) calculated on a total of 24 RV: Copepodits3 ClausocalanusB Copepodits1 0.5938735 0.7607149 0.8280316 EchinodermsLarvae AdultsOfCentropages Nauplii 0.8620636 0.8981089 0.9039953 ClausocalanusC Copepodits5 Siphonophores 0.9109498 0.9171680 0.9254146 Acartia GasteropodsLarvae Copepodits2 0.9324560 0.9399177 0.9469293 Ostracods EggsOfCrustaceans Copepodits4 0.9528440 0.9584662 0.9621404 Oithona Pteropods BellsOfCalycophores 0.9679197 0.9727763 0.9773901 Acanthaires DecapodsLarvae ClausocalanusA 0.9811460 0.9850263 0.9882042 Cladocerans AdultsOfCalanus JuvenilesOfCentropages 0.9921904 0.9952681 1.0000000 > plot(marbio.esc) > # The x-axis has short labels. For more info., enter: > marbio.esc$vr Copepodits3 ClausocalanusB Copepodits1 5 9 3 EchinodermsLarvae AdultsOfCentropages Nauplii 17 11 13 ClausocalanusC Copepodits5 Siphonophores 10 7 23 Acartia GasteropodsLarvae Copepodits2 1 19 4 Ostracods EggsOfCrustaceans Copepodits4 21 20 6 Oithona Pteropods BellsOfCalycophores 14 22 24 Acanthaires DecapodsLarvae ClausocalanusA 15 18 8 Cladocerans AdultsOfCalanus JuvenilesOfCentropages 16 2 12 > # Define a level at which to extract most significant variables > marbio.esc$level <- 0.90 > # Show it on the graph > lines(marbio.esc) > # This can also be done interactively on the plot using: > # marbio.esc$level <- identify(marbio.esc) > # Finally, extract most significant variables > marbio2 <- extract(marbio.esc) Warning in UseMethod("extract", e, n, ...) : arguments after the first two are ignored > names(marbio2) [1] "Copepodits3" "ClausocalanusB" "Copepodits1" [4] "EchinodermsLarvae" "AdultsOfCentropages" > > > > cleanEx(); ..nameEx <- "first" > > ### * first > > flush(stderr()); flush(stdout()) > > ### Name: first > ### Title: Get the first element of a vector > ### Aliases: first > ### Keywords: manip > > ### ** Examples > > a <- c(NA, 1, 2, NA, 3, 4, NA) > first(a) [1] NA > first(a, na.rm=TRUE) [1] 1 > > > > cleanEx(); ..nameEx <- "is.tseries" > > ### * is.tseries > > flush(stderr()); flush(stdout()) > > ### Name: is.tseries > ### Title: Is this object a time series? > ### Aliases: is.tseries > ### Keywords: ts > > ### ** Examples > > tser <- ts(sin((1:100)/6*pi)+rnorm(100, sd=0.5), start=c(1998, 4), frequency=12) > is.tseries(tser) # TRUE [1] TRUE > not.a.ts <- c(1,2,3) > is.tseries(not.a.ts) # FALSE [1] FALSE > > > > cleanEx(); ..nameEx <- "last" > > ### * last > > flush(stderr()); flush(stdout()) > > ### Name: last > ### Title: Get the last element of a vector > ### Aliases: last > ### Keywords: manip > > ### ** Examples > > a <- c(NA, 1, 2, NA, 3, 4, NA) > last(a) [1] NA > last(a, na.rm=TRUE) [1] 4 > > > > cleanEx(); ..nameEx <- "local.trend" > > ### * local.trend > > flush(stderr()); flush(stdout()) > > ### Name: local.trend > ### Title: Calculate local trends using cumsum > ### Aliases: local.trend identify.local.trend > ### Keywords: ts > > ### ** Examples > > data(bnr) > # Calculate and plot cumsum for the 8th series > bnr8.lt <- local.trend(bnr[,8]) > # To identify local trends, use: > # identify(bnr8.lt) > # and click points between which you want to compute local linear trends... > > > > cleanEx(); ..nameEx <- "match.tol" > > ### * match.tol > > flush(stderr()); flush(stdout()) > > ### Name: match.tol > ### Title: Determine matching observation with a tolerance in time-scale > ### Aliases: match.tol > ### Keywords: ts manip > > ### ** Examples > > X <- 1:5 > Table <- c(1, 3.1, 3.8, 4.4, 5.1, 6) > match.tol(X, Table) [1] 1 NA NA NA NA > match.tol(X, Table, tol=0.2) [1] 1 NA 2 NA 5 > match.tol(X, Table, tol=0.55) 'tol' was adjusted to 0.5 [1] 1 NA 2 3 5 > > > > cleanEx(); ..nameEx <- "pennington" > > ### * pennington > > flush(stderr()); flush(stdout()) > > ### Name: pennington > ### Title: Calculate Pennington statistics > ### Aliases: pennington > ### Keywords: univar > > ### ** Examples > > data(marbio) > pennington(marbio[, "Copepodits2"]) mean var mean.var 38.60456 1302.88215 18.47618 > pennington(marbio[, "Copepodits2"], calc="mean", na.rm=TRUE) [1] 38.60456 > > > > cleanEx(); ..nameEx <- "pgleissberg" > > ### * pgleissberg > > flush(stderr()); flush(stdout()) > > ### Name: pgleissberg > ### Title: Gleissberg distribution probability > ### Aliases: pgleissberg > ### Keywords: distribution > > ### ** Examples > > # Until n=50, the exact probability is returned > pgleissberg(20, 10, lower.tail=TRUE, two.tailed=FALSE) [1] 0.2011741 > # For higher n values, it is approximated by a normal distribution > pgleissberg(60, 33, lower.tail=TRUE, two.tailed=FALSE) [1] 0.03904556 > > > > cleanEx(); ..nameEx <- "regarea" > > ### * regarea > > flush(stderr()); flush(stdout()) > > ### Name: regarea > ### Title: Regulate a series using the area method > ### Aliases: regarea > ### Keywords: ts manip chron smooth > > ### ** Examples > > data(releve) > # The 'Melosul' series is regulated with a 25-days window > reg <- regarea(releve$Day, releve$Melosul, window=25) > # Then with a 50-days window > reg2 <- regarea(releve$Day, releve$Melosul, window=50) > # The initial and both regulated series are shown on the graph for comparison > plot(releve$Day, releve$Melosul, type="l") > lines(reg$x, reg$y, col=2) > lines(reg2$x, reg2$y, col=4) > > > > cleanEx(); ..nameEx <- "regconst" > > ### * regconst > > flush(stderr()); flush(stdout()) > > ### Name: regconst > ### Title: Regulate a series using the constant value method > ### Aliases: regconst > ### Keywords: ts manip chron smooth > > ### ** Examples > > data(releve) > reg <- regconst(releve$Day, releve$Melosul) > plot(releve$Day, releve$Melosul, type="l") > lines(reg$x, reg$y, col=2) > > > > cleanEx(); ..nameEx <- "reglin" > > ### * reglin > > flush(stderr()); flush(stdout()) > > ### Name: reglin > ### Title: Regulation of a series using a linear interpolation > ### Aliases: reglin > ### Keywords: ts manip chron smooth > > ### ** Examples > > data(releve) > reg <- reglin(releve$Day, releve$Melosul) > plot(releve$Day, releve$Melosul, type="l") > lines(reg$x, reg$y, col=2) > > > > cleanEx(); ..nameEx <- "regspline" > > ### * regspline > > flush(stderr()); flush(stdout()) > > ### Name: regspline > ### Title: Regulation of a time series using splines > ### Aliases: regspline > ### Keywords: ts manip chron smooth > > ### ** Examples > > data(releve) > reg <- regspline(releve$Day, releve$Melosul) > plot(releve$Day, releve$Melosul, type="l") > lines(reg$x, reg$y, col=2) > > > > cleanEx(); ..nameEx <- "regul" > > ### * regul > > flush(stderr()); flush(stdout()) > > ### Name: regul > ### Title: Regulation of one or several time series using various methods > ### Aliases: regul extract.regul hist.regul identify.regul lines.regul > ### plot.regul print.regul print.specs.regul print.summary.regul > ### specs.regul summary.regul > ### Keywords: ts manip chron smooth > > ### ** Examples > > data(releve) > # The series in this data frame are very irregularly sampled in time: > releve$Day [1] 1 51 108 163 176 206 248 315 339 356 389 449 480 493 501 [16] 508 522 554 568 597 613 624 639 676 697 723 751 786 814 842 [31] 863 877 891 906 922 940 954 971 983 999 1010 1027 1038 1054 1069 [46] 1081 1094 1110 1129 1143 1156 1173 1186 1207 1226 1235 1249 1271 1290 1314 [61] 1325 > length(releve$Day) [1] 61 > intervals <- releve$Day[2:61]-releve$Day[1:60] > intervals [1] 50 57 55 13 30 42 67 24 17 33 60 31 13 8 7 14 32 14 29 16 11 15 37 21 26 [26] 28 35 28 28 21 14 14 15 16 18 14 17 12 16 11 17 11 16 15 12 13 16 19 14 13 [51] 17 13 21 19 9 14 22 19 24 11 > range(intervals) [1] 7 67 > mean(intervals) [1] 22.06667 > # The series must be regulated to be converted in a 'rts' or 'ts object > rel.reg <- regul(releve$Day, releve[3:8], xmin=9, n=63, deltat=21, + tol=1.05, methods=c("s","c","l","a","s","a"), window=21) > rel.reg Regulation of, by "method" : Astegla Chae Dity Gymn Melosul Navi "s" "c" "l" "a" "s" "a" Arguments for "methods" : tol.type tol rule f periodic window split "both" "1.05" "1" "0" "FALSE" "21" "100" 44 interpolated values on 63 ( 0 NAs padded at ends ) Time scale : start deltat frequency 9.00000000 21.00000000 0.04761905 Time units : days call : regul(x = releve$Day, y = releve[3:8], xmin = 9, n = 63, deltat = 21, tol = 1.05, methods = c("s", "c", "l", "a", "s", "a"), window = 21) > plot(rel.reg, 5) > specs(rel.reg) Warning in UseMethod("specs", x, ...) : arguments after the first two are ignored xmin n frequency deltat 9 63 0.0476190476190476 21 units dateformat tol tol.type days m/d/Y 1.05 both methods1 methods2 methods3 methods4 s c l a methods5 methods6 rule f s a 1 0 periodic window split FALSE 21 100 > # Now we can extract one or several regular time series > melo.ts <- extract(rel.reg, series="Melosul") Warning in UseMethod("extract", e, n, ...) : arguments after the first two are ignored > is.tseries(melo.ts) [1] TRUE > > # One can convert time-scale from "days" to "years" during regulation > # This is most useful for analyzing seasonal cycles in a second step > melo.regy <- regul(releve$Day, releve$Melosul, xmin=6, n=87, + units="daystoyears", frequency=24, tol=2.2, methods="linear", + datemin="21/03/1989", dateformat="d/m/Y") A 'tol' of 2.2 in 'days' is 0.00602327173169062 in 'years' 'tol' was adjusted to 0.00595238095238083 > melo.regy Regulation of, by "method" : Series "linear" Arguments for "methods" : tol.type tol rule f "both" "2.2" "1" "0" periodic window split "FALSE" "15.3953488372093" "100" 63 interpolated values on 87 ( 0 NAs padded at ends ) Time scale : start deltat frequency 2.109000e+03 4.166667e-02 2.400000e+01 Time units : years call : regul(x = releve$Day, y = releve$Melosul, xmin = 6, n = 87, units = "daystoyears", frequency = 24, datemin = "21/03/1989", dateformat = "d/m/Y", tol = 2.2, methods = "linear") > plot(melo.regy, main="Regulation of Melosul") > # In this case, we have only one series in 'melo.regy' > # We can use also 'tseries()' instead of 'extract()' > melo.tsy <- tseries(melo.regy) > is.tseries(melo.tsy) [1] TRUE > > > > cleanEx(); ..nameEx <- "regul.adj" > > ### * regul.adj > > flush(stderr()); flush(stdout()) > > ### Name: regul.adj > ### Title: Adjust regulation parameters > ### Aliases: regul.adj > ### Keywords: ts chron > > ### ** Examples > > # This example follows the example for regul.screen() > # where we determined that xmin=9, deltat=21, n=63, with tol=1.05 > # is a good choice to regulate the irregular time series in 'releve' > data(releve) > regul.adj(releve$Day, xmin=9, deltat=21) $params xmin n deltat tol 9 63 21 21 $match [1] 59 $exact.match [1] 10 $match.counts 0 1 2 3 4 5 6 7 8 9 10 12 13 14 15 16 19 20 Inf 10 19 21 23 24 29 32 41 46 48 49 51 53 55 56 57 58 59 63 > # The histogram indicates that it is not useful to increase tol > # more than 1.05, because few observations will be added > # except if we increase it to 5-7, but this value could be > # considered to be too large in comparison with deltat=22 > # On the other hand, with tol <= 1, the number of matching > # observations will be almost divided by two! > > > > cleanEx(); ..nameEx <- "regul.screen" > > ### * regul.screen > > flush(stderr()); flush(stdout()) > > ### Name: regul.screen > ### Title: Test various regulation parameters > ### Aliases: regul.screen > ### Keywords: ts chron > > ### ** Examples > > data(releve) > # This series is very irregular, and it is difficult > # to choose the best regular time-scale > releve$Day [1] 1 51 108 163 176 206 248 315 339 356 389 449 480 493 501 [16] 508 522 554 568 597 613 624 639 676 697 723 751 786 814 842 [31] 863 877 891 906 922 940 954 971 983 999 1010 1027 1038 1054 1069 [46] 1081 1094 1110 1129 1143 1156 1173 1186 1207 1226 1235 1249 1271 1290 1314 [61] 1325 > length(releve$Day) [1] 61 > intervals <- releve$Day[2:61]-releve$Day[1:60] > intervals [1] 50 57 55 13 30 42 67 24 17 33 60 31 13 8 7 14 32 14 29 16 11 15 37 21 26 [26] 28 35 28 28 21 14 14 15 16 18 14 17 12 16 11 17 11 16 15 12 13 16 19 14 13 [51] 17 13 21 19 9 14 22 19 24 11 > range(intervals) [1] 7 67 > mean(intervals) [1] 22.06667 > # A combination of xmin=1, deltat=22 and n=61 seems correct > # But is it the best one? > regul.screen(releve$Day, xmin=0:11, deltat=16:27, tol=1.05) $tol d=16 d=17 d=18 d=19 d=20 d=21 d=22 d=23 1.066667 1.062500 1.058824 1.055556 1.052632 1.050000 1.047619 1.045455 d=24 d=25 d=26 d=27 1.043478 1.041667 1.040000 1.038462 $n d=16 d=17 d=18 d=19 d=20 d=21 d=22 d=23 d=24 d=25 d=26 d=27 x=0 83 78 74 70 67 64 61 58 56 54 51 50 x=1 83 78 74 70 67 64 61 58 56 53 51 50 x=2 83 78 74 70 67 64 61 58 56 53 51 50 x=3 83 78 74 70 67 63 61 58 56 53 51 49 x=4 83 78 74 70 67 63 61 58 56 53 51 49 x=5 83 78 74 70 67 63 61 58 56 53 51 49 x=6 83 78 74 70 66 63 60 58 55 53 51 49 x=7 83 78 74 70 66 63 60 58 55 53 51 49 x=8 83 78 74 70 66 63 60 58 55 53 51 49 x=9 83 78 74 70 66 63 60 58 55 53 51 49 x=10 83 78 74 70 66 63 60 58 55 53 51 49 x=11 83 78 74 70 66 63 60 58 55 53 51 49 $nbr.match d=16 d=17 d=18 d=19 d=20 d=21 d=22 d=23 d=24 d=25 d=26 d=27 x=0 9 12 13 6 8 5 5 5 10 12 10 12 x=1 10 15 12 6 8 9 5 9 11 9 7 12 x=2 13 9 11 10 11 11 8 11 14 9 8 13 x=3 13 10 10 7 9 12 10 11 13 6 8 7 x=4 14 9 8 10 7 9 11 7 11 5 7 9 x=5 10 11 6 6 7 6 8 6 8 10 5 5 x=6 12 13 4 7 7 1 7 8 5 7 6 5 x=7 9 11 5 8 11 4 9 6 6 9 7 2 x=8 10 11 9 11 11 14 12 6 4 4 7 5 x=9 13 11 9 13 12 19 9 6 4 6 4 3 x=10 14 10 13 13 13 17 7 10 4 5 5 7 x=11 14 9 14 12 8 9 3 9 5 5 7 6 $nbr.exact.match d=16 d=17 d=18 d=19 d=20 d=21 d=22 d=23 d=24 d=25 d=26 d=27 x=0 3 7 4 1 2 1 2 2 2 1 2 3 x=1 3 3 6 4 3 3 2 1 5 7 3 7 x=2 4 5 2 1 3 5 1 6 4 2 2 2 x=3 6 1 3 5 5 3 5 4 5 0 3 4 x=4 3 4 5 1 1 5 4 1 4 4 3 2 x=5 5 4 0 4 1 1 2 2 2 1 1 3 x=6 2 3 1 1 5 0 2 3 2 5 1 0 x=7 5 6 3 2 2 0 4 3 2 1 4 2 x=8 2 2 1 5 4 4 3 0 2 3 2 0 x=9 3 3 5 4 5 10 5 3 0 0 1 3 x=10 8 6 3 4 3 5 1 3 2 3 1 0 x=11 3 1 5 5 5 2 1 4 2 2 3 4 > # Now we can tell that xmin=9, deltat=21, n=63, with tol=1.05 > # is a much better choice! > > > > cleanEx(); ..nameEx <- "stat.desc" > > ### * stat.desc > > flush(stderr()); flush(stdout()) > > ### Name: stat.desc > ### Title: Descriptive statistics on a data frame or time series > ### Aliases: stat.desc > ### Keywords: ts > > ### ** Examples > > data(marbio) > stat.desc(marbio[,13:16], basic=TRUE, desc=TRUE, norm=TRUE, p=0.95) Nauplii Oithona Acanthaires Cladocerans nbr.val 6.800000e+01 6.800000e+01 6.800000e+01 6.800000e+01 nbr.null 0.000000e+00 0.000000e+00 8.000000e+00 4.500000e+01 nbr.na 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 min 3.000000e+00 3.900000e+01 0.000000e+00 0.000000e+00 max 6.360000e+02 7.792000e+03 5.200000e+01 2.000000e+02 range 6.330000e+02 7.753000e+03 5.200000e+01 2.000000e+02 sum 9.018000e+03 1.051590e+05 9.870000e+02 3.940000e+02 median 1.040000e+02 1.228000e+03 1.200000e+01 0.000000e+00 mean 1.326176e+02 1.546456e+03 1.451471e+01 5.794118e+00 SE.mean 1.254198e+01 1.554392e+02 1.458903e+00 3.035968e+00 CI.mean.0.95 2.503389e+01 3.102580e+02 2.911983e+00 6.059818e+00 var 1.069648e+04 1.642972e+06 1.447311e+02 6.267629e+02 std.dev 1.034238e+02 1.281785e+03 1.203042e+01 2.503523e+01 coef.var 7.798644e-01 8.288530e-01 8.288439e-01 4.320802e+00 skewness 2.037948e+00 2.579061e+00 9.098177e-01 6.945251e+00 skew.2SE 3.504462e+00 4.434961e+00 1.564525e+00 1.194307e+01 kurtosis 6.433595e+00 8.866041e+00 5.389480e-01 5.039083e+01 kurt.2SE 5.604129e+00 7.722966e+00 4.694629e-01 4.389407e+01 normtest.W 8.239310e-01 7.538450e-01 9.182656e-01 2.253480e-01 normtest.p 1.507169e-07 2.457429e-09 2.736212e-04 3.745690e-17 > > > > cleanEx(); ..nameEx <- "stat.pen" > > ### * stat.pen > > flush(stderr()); flush(stdout()) > > ### Name: stat.pen > ### Title: Pennington statistics on a data frame or time series > ### Aliases: stat.pen > ### Keywords: ts > > ### ** Examples > > data(marbio) > stat.pen(marbio[,c(4, 14:16)], basic=TRUE, desc=TRUE) Copepodits2 Oithona Acanthaires Cladocerans nbr.val 68.000000 68.000 68.000000 68.000000 nbr.null 5.000000 0.000 8.000000 45.000000 percnull 7.352941 0.000 11.764706 66.176471 nbr.na 0.000000 0.000 0.000000 0.000000 median 24.500000 1228.000 12.000000 0.000000 mean 39.632353 1546.456 14.514706 5.794118 var 1619.877744 1642971.655 144.731124 626.762950 std.dev 40.247705 1281.785 12.030425 25.035234 pos.median 28.000000 1228.000 13.000000 4.000000 pos.mean 42.777778 1546.456 16.450000 17.130435 pos.var 1613.788530 1642971.655 131.980508 1705.754941 pos.std.dev 40.171987 1281.785 11.488277 41.300786 geo.mean 31.504949 1161.198 12.030807 6.482851 pen.mean 38.604561 1624.345 16.011929 4.183929 pen.var 1302.882154 2463485.626 389.856084 144.116540 pen.std.dev 36.095459 1569.549 19.744774 12.004855 pen.mean.var 18.476177 34230.002 5.325181 1.900941 > > > > cleanEx(); ..nameEx <- "stat.slide" > > ### * stat.slide > > flush(stderr()); flush(stdout()) > > ### Name: stat.slide > ### Title: Sliding statistics > ### Aliases: stat.slide lines.stat.slide plot.stat.slide print.stat.slide > ### Keywords: ts > > ### ** Examples > > data(marbio) > # Sliding statistics with fixed-length blocs > 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.000 50.000 60.0000 xmax 10.000000 20.0000 30.0000 40.0000 50.000 60.000 70.0000 nbr.val 9.000000 10.0000 10.0000 10.0000 10.000 10.000 9.0000 nbr.null 2.000000 0.0000 0.0000 0.0000 0.000 0.000 0.0000 nbr.na 0.000000 0.0000 0.0000 0.0000 0.000 0.000 0.0000 min 0.000000 160.0000 158.0000 96.0000 136.000 56.000 52.0000 max 12.000000 832.0000 1980.0000 1204.0000 2484.000 824.000 656.0000 median 4.000000 344.0000 732.0000 752.0000 260.000 340.000 120.0000 mean 4.777778 350.8000 810.1000 682.0000 494.400 364.800 219.1111 std.dev 4.790036 196.9708 619.8578 350.0616 711.051 234.915 202.4231 > plot(statsl, stat="mean", leg=TRUE, lpos=c(55, 2500), xlab="Station", + ylab="ClausocalanusA") > > # More information on the series, with predefined blocs > 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.300444e-01 8.162687e-01 9.560890e-01 normtest.p 7.467714e-04 5.943045e-02 1.092400e-01 7.221453e-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.00000000 4.600000e+01 xmax 46.00000000 7.000000e+01 nbr.val 5.00000000 2.300000e+01 nbr.null 0.00000000 0.000000e+00 nbr.na 0.00000000 0.000000e+00 min 136.00000000 5.200000e+01 max 264.00000000 2.484000e+03 range 128.00000000 2.432000e+03 sum 1088.00000000 9.236000e+03 median 256.00000000 3.080000e+02 mean 217.60000000 4.015652e+02 SE.mean 27.29395537 1.048611e+02 CI.mean.0.95 75.78016880 2.174685e+02 var 3724.80000000 2.529043e+05 std.dev 61.03113959 5.028960e+02 coef.var 0.28047399 1.252339e+00 skewness -0.35978254 3.033043e+00 skew.2SE -0.19706101 3.150646e+00 kurtosis -2.10864917 9.865540e+00 kurt.2SE -0.52716229 5.277023e+00 normtest.W 0.78640766 6.024972e-01 normtest.p 0.06250056 9.610843e-07 pos.median 256.00000000 3.080000e+02 pos.mean 217.60000000 4.015652e+02 geo.mean 209.92260248 2.527195e+02 pen.mean 218.05538851 3.905421e+02 pen.var 4518.02570358 1.913652e+05 pen.std.dev 67.21626071 4.374531e+02 pen.mean.var 903.39387927 7.745030e+03 > plot(statsl2, stat="median", xlab="Stations", ylab="Counts", + main="Clausocalanus A") # Median > lines(statsl2, stat="min") # Minimum > lines(statsl2, stat="max") # Maximum > lines(c(17, 17), c(-50, 2600), col=4, lty=2) # Cuts > 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.5, 21, 27.5, 35, 43.5, 57), 2300, labels=c("Peripheral Zone", "D1", + "C", "Front", "D2", "Central Zone")) # Labels > legend(0, 1900, c("series", "median", "range"), col=1:3, lty=1) > # Get cuts back from the object > statsl2$xcut [1] 0 17 25 30 41 46 70 > > > > cleanEx(); ..nameEx <- "trend.test" > > ### * trend.test > > flush(stderr()); flush(stdout()) > > ### Name: trend.test > ### Title: Test if an increasing or decreasing trend exists in a time > ### series > ### Aliases: trend.test > ### Keywords: ts htest > > ### ** Examples > > data(marbio) > trend.test(marbio[, 8]) Warning in cor.test.default(x[, i], Time, alternative = "two.sided", method = "spearman") : p-values may be incorrect due to ties $"Series 1" Spearman's rank correlation rho data: Series 1 and time(Series 1) S = 43853, p-value = 0.1836 alternative hypothesis: true rho is not equal to 0 sample estimates: rho 0.1630113 > # Run a bootstrap test on the same series > marbio8.trend.test <- trend.test(marbio[, 8], R=99) > # R=999 is a better value... but it is very slow! > marbio8.trend.test BLOCK BOOTSTRAP FOR TIME SERIES Fixed Block Length of 1 Call: tsboot(tseries = x, statistic = test.trend, R = R, l = 1, sim = "fixed") Bootstrap Statistics : original bias std. error t1* 0.1630113 -0.1664767 0.1123223 > plot(marbio8.trend.test) > marbio8.trend.test$p.value [1] 0.0909091 > > > > cleanEx(); ..nameEx <- "tsd" > > ### * tsd > > flush(stderr()); flush(stdout()) > > ### Name: tsd > ### Title: Decomposition of one or several regular time series using > ### various methods > ### Aliases: tsd extract.tsd plot.tsd print.specs.tsd print.summary.tsd > ### print.tsd specs.tsd summary.tsd > ### Keywords: ts smooth loess nonparametric > > ### ** Examples > > data(releve) > # Regulate the series and extract them as a time series object > rel.regy <- regul(releve$Day, releve[3:8], xmin=6, n=87, units="daystoyears", + frequency=24, tol=2.2, methods="linear", datemin="21/03/1989", + dateformat="d/m/Y") A 'tol' of 2.2 in 'days' is 0.00602327173169062 in 'years' 'tol' was adjusted to 0.00595238095238083 > rel.ts <- tseries(rel.regy) > > # Decompose all series in the set with the "loess" method > rel.dec <- tsd(rel.ts, method="loess", s.window=13, trend=FALSE) > rel.dec Call: tsd(x = rel.ts, method = "loess", s.window = 13, trend = FALSE) Series: [1] "Astegla" "Chae" "Dity" "Gymn" "Melosul" "Navi" Components for Astegla [1] "deseasoned" "seasonal" Components for Chae [1] "deseasoned" "seasonal" Components for Dity [1] "deseasoned" "seasonal" Components for Gymn [1] "deseasoned" "seasonal" Components for Melosul [1] "deseasoned" "seasonal" Components for Navi [1] "deseasoned" "seasonal" > plot(rel.dec, series=5, col=1:3) # An plot series 5 > > # Extract "deseasoned" components > rel.des <- extract(rel.dec, series=3:6, components="deseasoned") Warning in UseMethod("extract", e, n, ...) : arguments after the first two are ignored > rel.des[1:10,] Dity.deseasoned Gymn.deseasoned Melosul.deseasoned Navi.deseasoned [1,] -1856.450218 728.08507 11232.346 1263.204 [2,] -326.203029 925.87962 4857.551 2193.775 [3,] -5.397257 958.98367 12218.422 2141.072 [4,] 305.799403 725.75255 19485.965 2042.397 [5,] 510.738352 389.66511 19241.367 2524.486 [6,] 673.469410 46.52693 16485.644 2800.179 [7,] 666.811431 -205.96157 17601.801 2505.898 [8,] 668.430789 928.90368 16997.344 2273.421 [9,] 694.149206 -1949.28870 17367.641 1813.558 [10,] 696.564421 -1890.15039 16520.938 1603.852 > > # Further decompose these components with a moving average > rel.des.dec <- tsd(rel.des, method="average", order=2, times=10) > plot(rel.des.dec, series=3, col=c(2, 4, 6)) > # In this case, a superposed graph is more appropriate: > plot(rel.des.dec, series=3, col=c(2,4), stack=FALSE, resid=FALSE, + labels=c("without season cycle", "trend"), lpos=c(0, 55000)) > # Extract residuals from the latter decomposition > rel.res2 <- extract(rel.des.dec, components="residuals") Warning in UseMethod("extract", e, n, ...) : arguments after the first two are ignored > > > > cleanEx(); ..nameEx <- "tseries" > > ### * tseries > > flush(stderr()); flush(stdout()) > > ### Name: tseries > ### Title: Convert a 'regul' or a 'tsd' object into a time series > ### Aliases: tseries > ### Keywords: ts manip > > ### ** Examples > > data(releve) > rel.regy <- regul(releve$Day, releve[3:8], xmin=6, n=87, units="daystoyears", + frequency=24, tol=2.2, methods="linear", datemin="21/03/1989", + dateformat="d/m/Y") A 'tol' of 2.2 in 'days' is 0.00602327173169062 in 'years' 'tol' was adjusted to 0.00595238095238083 > # This object is not a time series > is.tseries(rel.regy) # FALSE [1] FALSE > # Extract all time series contained in the 'regul' object > rel.ts <- tseries(rel.regy) > # Now this is a time series > is.tseries(rel.ts) # TRUE [1] TRUE > > > > cleanEx(); ..nameEx <- "turnogram" > > ### * turnogram > > flush(stderr()); flush(stdout()) > > ### Name: turnogram > ### Title: Calculate and plot a turnogram for a regular time series > ### Aliases: turnogram extract.turnogram identify.turnogram plot.turnogram > ### print.summary.turnogram print.turnogram summary.turnogram > ### Keywords: ts htest > > ### ** Examples > > data(bnr) > # Let's transform series 4 into a time series (supposing it is regular) > bnr4 <- as.ts(bnr[, 4]) > plot(bnr4, type="l", main="bnr4: raw data", xlab="Time") > # A simple turnogram is calculated > bnr4.turno <- turnogram(bnr4) > summary(bnr4.turno) Simple turnogram for: bnr4 options : mean / two-tailed probability call : turnogram(series = bnr4) max. info. : 7.903453 at interval 3 (P = 0.004176608: 15 turning points for 35 observations) extract level: 3 () interval n turns info 1 1 103 61 2.8849788 2 2 52 32 0.6097183 3 3 35 15 7.9034527 4 4 26 11 5.8490835 5 5 21 11 1.2432696 6 6 18 11 -0.3733604 7 7 15 8 0.4169020 8 8 13 7 0.0000000 9 9 12 7 -0.4792596 10 10 11 7 -1.1157070 11 11 10 5 0.0000000 12 12 9 5 -0.5790350 13 13 8 5 -1.3917540 14 14 8 2 3.7388752 15 15 7 3 0.0000000 16 16 7 4 -1.7427020 17 17 7 2 1.7427020 18 18 6 2 0.7776076 19 19 6 2 0.7776076 20 20 6 3 -0.7776076 > # A complete turnogram confirms that "level=3" is a good value: > turnogram(bnr4, complete=TRUE) Complete turnogram for: bnr4 options : mean / two-tailed probability intervals : 1 .. 20 / step = 1 nbr of obs. : 103 .. 6 max. info. : 6.64499 at interval 3 (P = 0.009992147: 15.66667 turning points for 35 observations) extract level: 3 () > # Data with maximum info. are extracted (thus taking 1 every 3 observations) > bnr4.interv3 <- extract(bnr4.turno) Warning in UseMethod("extract", e, n, ...) : arguments after the first two are ignored > plot(bnr4, type="l", lty=2, xlab="Time") > lines(bnr4.interv3, col=2) > title("Original bnr4 (dotted) versus max. info. curve (plain)") > # Choose another level (for instance, 6) and extract the corresponding series > bnr4.turno$level <- 6 > bnr4.interv6 <- extract(bnr4.turno) Warning in UseMethod("extract", e, n, ...) : arguments after the first two are ignored > # plot both extracted series on top of the original one > plot(bnr4, type="l", lty=2, xlab="Time") > lines(bnr4.interv3, col=2) > lines(bnr4.interv6, col=3) > legend(70, 580, c("original", "interval=3", "interval=6"), col=1:3, lty=c(2, 1, 1)) > # It is hard to tell on the graph which series contains more information > # The turnogram shows us that it is the "interval=3" one! > > > > cleanEx(); ..nameEx <- "turnpoints" > > ### * turnpoints > > flush(stderr()); flush(stdout()) > > ### Name: turnpoints > ### Title: Analyze turning points (peaks or pits) > ### Aliases: turnpoints extract.turnpoints lines.turnpoints plot.turnpoints > ### print.summary.turnpoints print.turnpoints summary.turnpoints > ### Keywords: ts > > ### ** Examples > > data(marbio) > plot(marbio[, "Nauplii"], type="l") > # Calculate turning points for this series > Nauplii.tp <- turnpoints(marbio[, "Nauplii"]) > summary(Nauplii.tp) Turning points for: marbio[, "Nauplii"] nbr observations : 68 nbr ex-aequos : 2 nbr turning points: 45 (first point is a peak) E(p) = 44 Var(p) = 11.76667 (theoretical) point type proba info 1 3 peak 0.250000000 2.0000000 2 4 pit 0.666666667 0.5849625 3 5 peak 0.666666667 0.5849625 4 6 pit 0.250000000 2.0000000 5 8 peak 0.027777778 5.1699250 6 11 pit 0.100000000 3.3219281 7 12 peak 0.666666667 0.5849625 8 14 pit 0.250000000 2.0000000 9 16 peak 0.027777778 5.1699250 10 19 pit 0.007936508 6.9772799 11 22 peak 0.027777778 5.1699250 12 24 pit 0.250000000 2.0000000 13 25 peak 0.666666667 0.5849625 14 26 pit 0.666666667 0.5849625 15 27 peak 0.666666667 0.5849625 16 28 pit 0.666666667 0.5849625 17 29 peak 0.666666667 0.5849625 18 30 pit 0.666666667 0.5849625 19 31 peak 0.666666667 0.5849625 20 32 pit 0.666666667 0.5849625 21 33 peak 0.666666667 0.5849625 22 34 pit 0.666666667 0.5849625 23 35 peak 0.666666667 0.5849625 24 36 pit 0.666666667 0.5849625 25 37 peak 0.666666667 0.5849625 26 38 pit 0.666666667 0.5849625 27 39 peak 0.013888889 6.1699250 28 43 pit 0.027777778 5.1699250 29 44 peak 0.666666667 0.5849625 30 45 pit 0.013888889 6.1699250 31 49 peak 0.007936508 6.9772799 32 51 pit 0.250000000 2.0000000 33 52 peak 0.666666667 0.5849625 34 53 pit 0.666666667 0.5849625 35 54 peak 0.666666667 0.5849625 36 55 pit 0.666666667 0.5849625 37 56 peak 0.250000000 2.0000000 38 58 pit 0.250000000 2.0000000 39 59 peak 0.666666667 0.5849625 40 60 pit 0.666666667 0.5849625 41 62 peak 0.666666667 0.5849625 42 63 pit 0.666666667 0.5849625 43 64 peak 0.666666667 0.5849625 44 65 pit 0.250000000 2.0000000 45 67 peak 0.250000000 2.0000000 > plot(Nauplii.tp) > # Add envelope and median line to original data > plot(marbio[, "Nauplii"], type="l") > lines(Nauplii.tp) > # Note that lines() applies to the graph of original dataset!!! > title("Raw data, envelope maxi., mini. and median line") > > > > cleanEx(); ..nameEx <- "vario" > > ### * vario > > flush(stderr()); flush(stdout()) > > ### Name: vario > ### Title: Compute and plot a semi-variogram > ### Aliases: vario > ### Keywords: ts > > ### ** Examples > > data(bnr) > vario(bnr[, 4]) distance semivario 1 1 8933.284 2 2 9450.129 3 3 10609.430 4 4 11919.288 5 5 13466.546 6 6 15853.098 7 7 15666.052 8 8 15818.568 9 9 15811.266 10 10 15756.645 11 11 15066.837 12 12 17747.786 13 13 16389.750 14 14 18679.669 15 15 17961.426 16 16 18236.943 17 17 19139.413 18 18 18463.482 19 19 18376.732 20 20 18767.747 21 21 17480.317 22 22 17754.858 23 23 18407.075 24 24 18675.184 25 25 19656.635 26 26 17710.909 27 27 19933.375 28 28 21158.673 29 29 19820.446 30 30 20525.432 31 31 21885.569 32 32 20792.761 33 33 21654.493 34 34 22013.464 > > > > ### *