[R] Storing loop output from a function
Liaw, Andy
andy_liaw at merck.com
Sat Nov 27 17:24:35 CET 2004
> From: Uwe Ligges
>
> Andrew Kniss wrote:
>
> > I am attempting to write an R function to aid in time
> series diagnostics.
> > The tsdiag() works well, but I would prefer to get a plot
> with ACF, PACF,
> > and Ljung-Box p-values of residuals. In my attempt to get
> to that point, I
> > am writing a function to calculate Ljung-Box p-values for
> residuals at
> > various lag distances.
> >
> > ljung<-function(fit)
> > for(i in c(1:24,36,48))
> > {box<-(Box.test(fit$residuals, lag=i, type="Ljung-Box"))
> > print(c(i, box$p.value))}
> >
>
> You need to return() rather than print() the object.
Yes, but print() is supposed to return (invisibly) its argument(s)...
Andy
> Uwe Ligges
>
>
> > This is one of my first R function writing attempts, so my
> apologies if
> > there is an obvious mistake. The above function produces
> the desired effect
> > in printing the lags and p-values to be plotted [where fit
> is the result of
> > arima()]; however I cannot seem to get the output stored in
> a data.frame for
> > subsequent plotting. I have tried storing the output using
> various methods
> > including data.frame, write.table, cat, capture.output, all
> with no success.
> >
> >
> > e.g:
> > ljung.out<-capture.output(print(c(i, box$p.value)))
> > #I saw a suggestion similar to this one in a previous post
> to the list...
> >
> > Any hints on how I can get the output stored so that I may
> plot it later? I
> > am using v 1.9.1, RGui for Windows.
> > Thanks,
> >
> > Andrew Kniss
> > Assistant Research Scientist
> > University of Wyoming
> > Dept. 3354
> > 1000 E. University Ave.
> > Laramie, WY 82071
> > (307) 766-3949
> > akniss at uwyo.edu
> >
> >
> > Below is an arima fit that I have used in testing the function.
> >
> >
> >>dump("coal.fit", file=stdout())
> >
> > coal.fit <-
> > structure(list(coef = structure(c(0.69539198614687, 484.903589344614
> > ), .Names = c("ar1", "intercept")), sigma2 = 3109.45476265185,
> > var.coef = structure(c(0.0104024111201598, 0.302125085158484,
> > 0.302125085158484, 634.39981868771), .Dim = as.integer(c(2,
> > 2)), .Dimnames = list(c("ar1", "intercept"), c("ar1",
> "intercept"
> > ))), mask = c(TRUE, TRUE), loglik = -266.892361376605, aic =
> > 539.784722753209,
> > arma = as.integer(c(1, 0, 0, 0, 1, 0, 0)), residuals =
> > structure(c(60.4342567589239,
> > -127.383559378086, -14.9885854976147, 123.839062585504,
> > -56.6019914334983,
> > 35.7247594443981, 63.6906479431108, -28.1651273226733,
> > -6.91856808459544,
> > 8.90309567990135, -30.8784722646861, -91.1489687772519,
> > -103.345257968621,
> > -29.2770349660464, -20.9664426335713, -25.3512422872431,
> > 32.6086618928476, -6.98260117899269, -108.850345082021,
> > 4.60267757422564,
> > 38.6146462114696, 42.7187751257762, 79.9491758184327,
> 36.880952815858,
> > 62.0132089128299, -0.848550671576206, -15.6420872534077,
> > 111.955160137055, 13.5021374808082, -126.940710948639,
> 63.7127908071542,
> >
> > 27.4722158876983, -52.0448398629454, -15.4535767911051,
> > -73.4996569296364,
> > 46.7008221699102, 27.5464232088949, -2.40151233395178,
> > -80.5337684309237,
> > -20.8162335807335, -18.2070175530272, -33.9885854976147,
> > -5.94848967770536, 17.8390625855041, 0.109559098069901,
> > 39.5464232088949,
> > 30.2537838322858, 32.9551601370546, 13.4381043864110),
> .Tsp = c(1,
> > 49, 1), class = "ts"), call = stats:::arima(x = x,
> order = order,
> > seasonal = seasonal, include.mean = include.mean),
> series = "x",
> > code = as.integer(0), n.cond = 0, model = structure(list(
> > phi = 0.69539198614687, theta = numeric(0), Delta =
> numeric(0),
> > Z = 1, a = 60.0964106553856, P = structure(0, .Dim
> = as.integer(c(1,
> >
> > 1))), T = structure(0.69539198614687, .Dim =
> as.integer(c(1,
> > 1))), V = structure(1, .Dim = as.integer(c(1, 1))), h = 0,
> > Pn = structure(1, .Dim = as.integer(c(1, 1)))),
> .Names = c("phi",
> > "theta", "Delta", "Z", "a", "P", "T", "V", "h", "Pn")), x =
> > structure(as.integer(c(569,
> > 416, 422, 565, 484, 520, 573, 518, 501, 505, 468, 382, 310,
> > 334, 359, 372, 439, 446, 349, 395, 461, 511, 583, 590, 620,
> > 578, 534, 631, 600, 438, 516, 534, 467, 457, 392, 467, 500,
> > 493, 410, 412, 416, 403, 422, 459, 467, 512, 534, 552, 545
> > )), .Dim = as.integer(c(49, 1)), .Dimnames = list(NULL,
> "Coal"), .Tsp =
> > c(1,
> > 49, 1), class = "ts")), .Names = c("coef", "sigma2",
> "var.coef",
> > "mask", "loglik", "aic", "arma", "residuals", "call", "series",
> > "code", "n.cond", "model", "x"), class = "Arima")
> >
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
>
More information about the R-help
mailing list