[R] expression() - Superscript in y-axis, keeping line break in string
I've tried several different ways to accomplish this, but as yet to no avail. My y-axis for a plot has a rather long label, and thus I have been using /n to break it into two lines. However, to make it technically correct for publication, I also need to use superscript in the label. For example: par(oma=c(0,0,2,0),mar=c(5,6,0.25,2),lheight=1) plot(1:10, ylab=14C-glyphosate line1\n line2) will provide the text in two lines as I would like it. However, I am trying to keep those same line breaks when using expression() to get my superscript number. This will not work, as it aligns the 14C section with the bottom line of the expression making little sense to the reader. par(oma=c(0,0,2,0),mar=c(5,6,0.25,2),lheight=1) plot(1:10, ylab=expression( ^14*C*-glyphosate line1\n line2)) Is there a way to align the 14C portion of the expression with the top line of the string rather than the bottom line? Any suggestions are greatly appreciated. Andrew -- Andrew Kniss Assistant Research Scientist University of Wyoming Department of Plant Sciences [EMAIL PROTECTED] Office: (307) 766-3949 Fax:(307) 766-5549 __ R-help@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 and provide commented, minimal, self-contained, reproducible code.
RE: [R] power.anova.test for interaction effects
Using the data set posted at http://uwstudentweb.uwyo.edu/A/AKNISS/sbxherb.txt I obtained this output from R: data(sbxherb) attach(sbxherb) sbxherb$rep-factor(c(I,II,III)) sbxherb$var-factor(sbxherb$var) sbxherb$trt-factor(sbxherb$trt, labels=c(Check,Early,Late,Micro)) detach(sbxherb) attach(sbxherb) aov.sbx-aov(yield~rep + trt * var + Error(rep/(trt+var))) summary(aov.sbx) Error: rep Df Sum Sq Mean Sq rep 2 122.158 61.079 Error: rep:trt Df Sum Sq Mean Sq F value Pr(F) trt3 201.131 67.044 1.3096 0.355 Residuals 6 307.160 51.193 Error: rep:var Df Sum Sq Mean Sq F valuePr(F) var 36 4613.9 128.2 8.4392 1.425e-14 *** Residuals 72 1093.515.2 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Error: Within Df Sum Sq Mean Sq F value Pr(F) trt:var 108 956.628.86 0.9731 0.5574 Residuals 216 1966.049.10 I would like to use this output to help design next years study, which will be carried out with the same design. However, I would like to know what the sample size should be in order to have an 80% chance of finding a 10% difference in the interaction term trt:var at alpha=0.05. Is there a way to coerce power.anova.test() to make this calculation for me? I have tried using this function to give me the sample size needed to find differences given the within and between groups variance from the previous study, but with no success. If I have it calculate the power from this data set, it tells me I have a power of 1, which I doubt is true given the high p-value and wide confidence intervals around the estimates. power.anova.test(groups=148, + sig.level=0.05, + power=0.8, + within.var=9.10, + between.var=8.86) Error in uniroot(function(n) eval(p.body) - power, c(2, 1e+05)) : f() values at end points not of opposite sign power.anova.test(groups=148, + sig.level=0.05, + n=3, + within.var=9.10, + between.var=8.86) Balanced one-way analysis of variance power calculation groups = 148 n = 3 between.var = 8.86 within.var = 9.1 sig.level = 0.05 power = 1 NOTE: n is number in each group Any help is greatly appreciated. I apologize if I seem to be nagging on this issue, but I thought this post was much better at describing what I am after. Thanks in advance. Andrew Kniss University of Wyoming [EMAIL PROTECTED] -Original Message- From: Andrew Criswell [mailto:[EMAIL PROTECTED] Sent: Tuesday, February 22, 2005 9:30 PM To: [EMAIL PROTECTED] Subject: Re: [R] power.anova.test for interaction effects Dear Andrew, If you could supply us with your data, more help would possibly be forthcoming. Best wishes, Andrew Andrew Criswell, PhD Graduate School, Bangkok University __ R-help@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] RE: R-help Digest, Vol 24, Issue 22
I used SAS to analyze the data initially, since the data set was made up of several files when I received it, and I'm still not very good at manipulating data in R. I have posted the data set from one location at the following address: http://uwstudentweb.uwyo.edu/A/AKNISS/sxherb.txt var=cultivar trt=herbicide treatment yield=response variable of interest All plot# from 101 to 104 are rep 1, 201-204 rep 2, and 301 to 304 rep 3. It was the only file that was in an easy format for R to read at the moment, and was probably the most reliable trial of the two locations. I would like to use power.anova.test() with this data set to plan next years study (to get a sample size for each herb*var combination), but I'm not quite sure how that is done for an interaction effect. Do I just use the MS for herb*var as the between group variance and the MSE as the within group variance? Or do I need to somehow include other variance parameters in the model? The model for this location (split-block design): yield = rep + herb + var + herb*var ## all are fixed effects rep*herb = error term for herb rep*var = error term for cultivar residual = error term for herb*var I hope this attempt at my question was a little more clear. I appreciate any help that is offered. Andrew Kniss Assistant Research Scientist University of Wyoming Department of Plant Sciences 1000 E. Univesity Ave Laramie, WY 82071 USA [EMAIL PROTECTED] -Original Message- From: John Maindonald [mailto:[EMAIL PROTECTED] Sent: Tuesday, February 22, 2005 3:37 PM To: r-help@stat.math.ethz.ch Cc: [EMAIL PROTECTED] Subject: Re: R-help Digest, Vol 24, Issue 22 You need to give the model formula that gave your output. There are two sources of variation (at least), within and between locations; though it looks as though your analysis may have tried to account for this (but if so, the terms are not laid out in a way that makes for ready interpretation. The design is such (two locations) that you do not have much of a check that effects are consistent over locations. You need to check whether results really are similar for all cultivars and for all herbicides, so that it is legitimate to pool as happens in the overall analysis. If a herbicide:cultivar combination has little effect the variability may be large, while if it has a dramatic effect (kills everything!), there may be no variability to speak of. John Maindonald. On 22 Feb 2005, at 10:06 PM, [EMAIL PROTECTED] wrote: To: 'Bob Wheeler' [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Subject: RE: [R] power.anova.test for interaction effects Reply-To: [EMAIL PROTECTED] It's a rather complex model. A 37*4 factorial (37 cultivars[var]; 4 herbicide treatments[trt]) with three replications[rep] was carried out at two locations[loc], with different randomizations within each rep at each location. Source DF Error Term MS Loc 1 Trt*rep(loc)12314 Rep(loc) 4 Trt*rep(loc)1230.5 Trt 3 Trt*rep(loc)64.72 Trt*loc 3 Trt*rep(loc)33.42 Trt*rep(loc) 12 Residual76.78 Var 36 Var*trt*loc 93.91 Var*trt 108 Var*trt*loc 12.06 Var*trt*loc 144 Residual43.09 Residual575 NA 21.23 -Original Message- From: Bob Wheeler [mailto:[EMAIL PROTECTED] Sent: Monday, February 21, 2005 4:33 PM To: [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Subject: Re: [R] power.anova.test for interaction effects Your F value is so low as to make me suspect your model. Where did the 144 denominator degrees of freedom come from? John Maindonald email: [EMAIL PROTECTED] phone : +61 2 (6125)3473fax : +61 2(6125)5549 Centre for Bioinformation Science, Room 1194, John Dedman Mathematical Sciences Building (Building 27) Australian National University, Canberra ACT 0200. __ R-help@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] power.anova.test for interaction effects
This question will probably get me in trouble on theoretical grounds, but I will pose it anyway. The situation: I recently ran a field study looking for differences in sugarbeet cultivar tolerance to a specific herbicide. The study was set up so that 37 cultivars were treated with 4 different applications of the herbicide (37*4 factorial). In doing so, we found that the interaction effect was highly insignificant (ndf=108, ddf=144, F=0.28, p=1.). Now my problem is this... the study takes up an enormous amount of time, energy, and money (as you could guess with 37 cultivars in a field study). We need to determine weather it is worth the effort to repeat the study this summer (practically, it is not, but our funding source would like a more concrete demonstration). I decided to try using power.anova.test just as a starting point to see what our power was. My question is: is this valid to do on an interaction term? If I use power.anova.test with on the interaction term, this is what I get: ~ power.anova.test(groups=(37*4), n=3, between.var=12.06, ~+ within.var=21.23, sig.level=0.05) ~ ~ Balanced one-way analysis of variance power calculation ~ ~ groups = 148 ~ n = 3 ~between.var = 12.06 ~ within.var = 21.23 ~ sig.level = 0.05 ~ power = 1 ~ ~ NOTE: n is number in each group This would imply that given the variability we observed with 3 replications, we almost certainly would have found differences if they existed. But given what I have read on power analysis, a high p-value and wide confidence intervals nearly always suggest inadequate sample size. (Our 90% confidence intervals differed from the estimates by as much as 28%, when a 10% difference would be significant from a practical perspective.) So is this a valid approach? Or does the power.anova.test fall apart if using an interaction effect? Thank you in advance for any help or references you are willing to point me to. Best regards, Andrew Kniss Assistant Research Scientist University of Wyoming Department of Plant Sciences 1000 E. University Ave. Laramie, WY 82071 USA [EMAIL PROTECTED] __ R-help@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
RE: [R] power.anova.test for interaction effects
It's a rather complex model. A 37*4 factorial (37 cultivars[var]; 4 herbicide treatments[trt]) with three replications[rep] was carried out at two locations[loc], with different randomizations within each rep at each location. Source DF Error Term MS Loc 1 Trt*rep(loc)12314 Rep(loc) 4 Trt*rep(loc)1230.5 Trt 3 Trt*rep(loc)64.72 Trt*loc 3 Trt*rep(loc)33.42 Trt*rep(loc) 12 Residual76.78 Var 36 Var*trt*loc 93.91 Var*trt 108 Var*trt*loc 12.06 Var*trt*loc 144 Residual43.09 Residual575 NA 21.23 -Original Message- From: Bob Wheeler [mailto:[EMAIL PROTECTED] Sent: Monday, February 21, 2005 4:33 PM To: [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Subject: Re: [R] power.anova.test for interaction effects Your F value is so low as to make me suspect your model. Where did the 144 denominator degrees of freedom come from? Andrew Kniss wrote: This question will probably get me in trouble on theoretical grounds, but I will pose it anyway. The situation: I recently ran a field study looking for differences in sugarbeet cultivar tolerance to a specific herbicide. The study was set up so that 37 cultivars were treated with 4 different applications of the herbicide (37*4 factorial). In doing so, we found that the interaction effect was highly insignificant (ndf=108, ddf=144, F=0.28, p=1.). Now my problem is this... the study takes up an enormous amount of time, energy, and money (as you could guess with 37 cultivars in a field study). We need to determine weather it is worth the effort to repeat the study this summer (practically, it is not, but our funding source would like a more concrete demonstration). I decided to try using power.anova.test just as a starting point to see what our power was. My question is: is this valid to do on an interaction term? If I use power.anova.test with on the interaction term, this is what I get: ~ power.anova.test(groups=(37*4), n=3, between.var=12.06, ~+ within.var=21.23, sig.level=0.05) ~ ~ Balanced one-way analysis of variance power calculation ~ ~ groups = 148 ~ n = 3 ~between.var = 12.06 ~ within.var = 21.23 ~ sig.level = 0.05 ~ power = 1 ~ ~ NOTE: n is number in each group This would imply that given the variability we observed with 3 replications, we almost certainly would have found differences if they existed. But given what I have read on power analysis, a high p-value and wide confidence intervals nearly always suggest inadequate sample size. (Our 90% confidence intervals differed from the estimates by as much as 28%, when a 10% difference would be significant from a practical perspective.) So is this a valid approach? Or does the power.anova.test fall apart if using an interaction effect? Thank you in advance for any help or references you are willing to point me to. Best regards, Andrew Kniss Assistant Research Scientist University of Wyoming Department of Plant Sciences 1000 E. University Ave. Laramie, WY 82071 USA [EMAIL PROTECTED] __ R-help@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 -- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches. __ R-help@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] Storing loop output from a function
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))}
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
[EMAIL PROTECTED]
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)
__
[EMAIL PROTECTED] 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] tsdiag() titles
I am using the ts package to fit ARIMA models, and the tsdiag() function to plot diagnostics. In doing so I'm generating an awful lot of diagnostic plots of different models and different data sets all within the same R session. So my question is, is there an option in tsdiag() similar to main=Title that I can use? This would be quite helpful when I print out the plots, so I can tell which plot goes with a particular data set and model. I can't seem to find any examples where this has been done, and no options (other than gof.lag) are listed in the R manual. library(ts) data(tbills) #Treasury Bills attach(tbills) ts.tbills-ts(tbills) diff.tbills-diff(ts.tbills) #Differenced Series arima.diff.tbills.100-arima(ts.tbills, order=c(1,0,0)) win.metafile(HW_ARIMA/tbill1.wmf) tsdiag(arima.diff.tbills.100, main=Treasury Bills) #main= does not work, is there a way to name the plot? dev.off() Thanks for any help or ideas. Andrew Kniss Assistant Research Scientist University of Wyoming Dept. 3354 1000 E. University Ave. Laramie, WY 82071 (307) 766-3949 [EMAIL PROTECTED] __ [EMAIL PROTECTED] 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] X11 Time Series Decomposition
I apologize if a similar question has been asked previously, but I have been unable to find anything relevant in searches of previous posts. Is there an available package that will perform a census X-11 (or X-12) time series decomposition? Thanks, Andrew Kniss [EMAIL PROTECTED] Associate Research Scientist University of Wyoming [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
