Re: [R] ANOVA non-sphericity test and corrections (eg, Greenhouse-Geisser)
Darren, Further to Peter Dalgaard's help; Take a look at the example in --- library(car) ?Anova # note upper case 'A' The example in the Anova help page following the --- ## a multivariate linear model for repeated-measures data ## See ?OBrienKaiser for a description of the data set used in this example gives a workings for -- ## Greenhouse-Geisser and Huynh-Feldt Corrections ## for Departure from Compound Symmetry You might find that example helpful. John __ 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] extracting the mode of a vector
Beno?t L?t? wrote: Hello, I have an elementary question (for which I couldn't find the answer on the web or the help): how can I extract the mode (modal score) of a vector? Assuming that your vector contains only integers: v - sample(1:5, size=20, replace=T) v [1] 1 1 1 1 2 3 5 1 1 5 2 4 1 3 1 1 5 4 1 5 vt - table(v) as.numeric(names(vt[vt == max(vt)])) [1] 1 Cheers, Gad # or more succinctly, names(vt[which.max(vt)]) [1] 1 John __ 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] ordering boxplots according to median
Use reorder # boxplot with increasing order of medians s2-with(InsectSprays,reorder(spray,count,median)) with(InsectSprays,boxplot(count~s2)) # boxplot with decreasing order of medians s2-with(InsectSprays,reorder(spray,-count,median)) with(InsectSprays,boxplot(count~s2)) John Talloen, Willem wrote-- Dear R-users, Does anyone knows how I can order my serie of boxplots from lowest to highest median (which is much better for visualization purposes). thanks in advance, willem -- __ 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] How to find statistics like that.
Adai, I recently came across the following definition of a statistic which may be relevent to the discussion. John - Berans (2003) provocative definition of statistics as the study of algorithms for data analysis elevates computational considerations to the forefront of the field. It is apparent that the evolutionary success of statistical methods is to a significant degree determined by considerations of computational convenience. As a result,design and dissemination of statistical software has become an integral part of statistical research. from this it follows that a 'Statistic' is A mathematical function or algorithm for data analysis Duncan Murdoch wrote On 11/9/2005 10:01 PM, Adaikalavan Ramasamy wrote: I think an alternative is to use a p-value from F distribution. Even tough it is not a statistics, it is much easier to explain and popular than 1/F. Better yet to report the confidence intervals. Just curious about your usage: why do you say a p-value is not a statistic? Duncan Murdoch Adaikalavan Ramasamy replied - If my usage is wrong please correct me. Thank you. Here are my reason : 1. p-value is a (cumulative) probability and always ranges from 0 to 1. A test statistic depending on its definition can wider range of possible values. 2. A test statistics is one that is calculated from the data without the need of assuming a null distribution. Whereas to calculate p-values, you need to assume a null distribution or estimate it empirically using permutation techniques. 3. The directionality of a test statistics may be ignored. For example a t-statistics of -5 and 5 are equally interesting in a two-sided testing. But the smaller the p-value, more evidence against the null hypothesis. Regards, Adai __ 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] Extracting Variance Components
Mike, use --- VarCorr(lme.object) or for a user friendly output use varcomp from the 'ape' package-- require(ape) varcomp(lme.object) varcomp also allows scaling of components to unity (*100 gives %) and also allows for cumulative sum of components. Note. varcomp doesn't work for lmer objects. HTH, John -- Michel Friesenhahn wrote- Dear List, Is there a way to extract variance components from lmeObjects or summary.lme objects without using intervals()? For my purposes I don't need the confidence intervals which I'm obtaining using parametric bootstrap. Thanks, Mike __ 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] lme and lmer syntax
Ronaldo, According to Douglas Bates's paper in 'R' News, It would seem that the correct model for nested split plot random effects with lmer , in your example ,with x2 nested within x1, would be -- lmer(y~x1 + x2 +(1|x1)+(1|x1:x2)) Try it with your model any see how it compares with your aov and lme models, John m Seg 24 Out 2005 18:08, Doran, Harold escreveu: Ronaldo See the article on lmer pasted below for syntax. It is the only current source documenting the code. In lmer(), the nesting structure for the ranmdom effects is handled in a slightly different way. If your observations are nested as you note, then you can use lmer(y~x1 + x2 +(1|x1) + (1|x2), data) @Article{Rnews:Bates:2005, author = {Douglas Bates}, title = {Fitting Linear Mixed Models in {R}}, journal = {R News}, year = 2005, volume = 5, number = 1, pages = {27--30}, month = {May}, url= {http://CRAN.R-project.org/doc/Rnews/}, } Hi, I try this with a splitsplitplot example. I make the correct model with aov, lme do compare with lmer. But I cant make a correct model in lmer. Look that the aov and lme results are similars, but very different from lmer. In aov and lme is used the correct DF for each variable, in lmer it use a same DF for all? Denom=54. What is my mistake? Thanks Ronaldo __ 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] Subsetting a list
Dennis Try TEST[-3] [[1]] [1] A1 A2 [[2]] [1] B1 B2 for removing more than one element from the list (say 2 3) -- TEST[-c(2,3)] [[1]] [1] A1 A2 HTH John Dennis Fisher wrote--- Colleagues, I have created a list in the following manner: TEST- list(c(A1, A2), c(B1, B2), c(C1, C2)) I now want to delete one element from the list, e.g., the third. The command TEST[[3]] yields (as expected): [1] C1 C2 The command TEST[[-3]] yields: Error: attempt to select more than one element How can I accomplish delete one or more elements from this list? I am running R2.2.0 on a Linux platform. Dennis __ 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] Sorting a data frame by one of the variables
Leaf, using your example data as 'dat' below -- dat-read.table(clipboard,header=T) dat XYZ 1 22.0 24.0 4.3 2 2.3 3.4 5.3 3 57.2 23.4 34.0 #to order the data frame by say X (for column 1)-- dat1-dat[order(dat[,1]),] dat1 XYZ 2 2.3 3.4 5.3 1 22.0 24.0 4.3 3 57.2 23.4 34.0 By way of interest if you wanted to order EVERY column in ascending order then you could do a loop --- # to order all cols of dat by rows (ascending) dat2-dat for (i in 1:3) dat2[,i]-dat[order(dat[,i]),i] dat2 XYZ 1 2.3 3.4 4.3 2 22.0 23.4 5.3 3 57.2 24.0 34.0 I hope that helps, John = Leaf Sun wrote--- Dear all, I have a date frame like this: X Y Z 22 24 4.3 2.3 3.4 5.3 . 57.223.434 What my purpose is: to sort the data frame by either X, Y or Z. sample output is (sorted by X) : X Y Z 2.3 3.4 5.3 . .. 22 24 4.3 ... 57.2 23.4 34 I have no idea how to use sort, order or rank functions. Please help me out. Thanks! __ 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] Doubt about nested aov output
Ronaldo , It looks as though you have specified you model incorrectly. In the Rats example ,the Treatment is the only fixed effect,Rat and Liver are random effects In aov testing for sig of 'Means' of Random Effects is pointless and that is why 'p' values are not given.Further more the interaction between a Random Effect and a Fixed Effect is also a Random Effect. The 'aov' with error structure terms output reflects this by only giving 'p' values to Fixed Effects and their interactions. Note That the Fixed Effects of Block, Variety and their interaction Block:Variety are given p values while the Field'Random Effects have not. model - aov(Glycogen~Treatment+Error(Rat/Liver)) summary(model) Error: Rat Df Sum Sq Mean Sq F value Pr(F) #Rat is random effect Residuals 1 413.44 413.44 Error: Rat:Liver #Rat:Liver is Random effect Df Sum Sq Mean Sq F value Pr(F) Residuals 4 164.444 41.111 Error: Within Df Sum Sq Mean Sq F valuePr(F) Treatment 2 1557.56 778.78 18.251 8.437e-06 *** #Fixed effect Residuals 28 1194.78 42.67 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 I hope that this is of help. John __ 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] FW: Re: Doubt about nested aov output
Ronaldo, Further to my previous posting on your Glycogen nested aov model. Having read Douglas Bates' response and Reflected on his lmer analysis output of your aov nested model example as given.The Glycogen treatment has to be a Fixed Effect.If a 'treatment' isn't a Fixed Effect what is ? If Douglas Bates' lmer model is modified to treat Glycogen Treatment as a purely Fixed Effect,with Rat and the interaction Rat:Liver as random effects then-- model.lmer-lmer(Glycogen~Treatment+(1|Rat)+(1|Rat:Liver)) summary(model.lmer) Linear mixed-effects model fit by REML Formula: Glycogen ~ Treatment + (1 | Rat) + (1 | Rat:Liver) AIC BIClogLik MLdeviance REMLdeviance 239.095 248.5961 -113.5475 238.5439 227.095 Random effects: GroupsNameVariance Std.Dev. Rat:Liver (Intercept) 2.1238e-08 0.00014573 Rat (Intercept) 2.0609e+01 4.53976242 Residual 4.2476e+01 6.51733769 # of obs: 36, groups: Rat:Liver, 6; Rat, 2 Fixed effects: Estimate Std. Error DF t value Pr(|t|) (Intercept) 140.5000 3.7208 33 37.7607 2.2e-16 *** Treatment2 10.5000 2.6607 33 3.9463 0.0003917 *** Treatment3 -5. 2.6607 33 -2.0045 0.0532798 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) Trtmn2 Treatment2 -0.358 Treatment3 -0.358 0.500 anova(model.lmer) Analysis of Variance Table Df Sum Sq Mean Sq Denom F valuePr(F) Treatment 2 1557.56 778.78 33.00 18.335 4.419e-06 *** which agrees with the aov model below. model - aov(Glycogen~Treatment+Error(Rat/Liver)) summary(model) John -Original Message- From: John Wilkinson (pipex) [mailto:[EMAIL PROTECTED] Sent: 07 September 2005 12:04 PM To: Ronaldo Reis-Jr. Cc: r-help Subject: Re: [R] Doubt about nested aov output Ronaldo , It looks as though you have specified your model incorrectly. In the Rats example ,the Treatment is the only fixed effect,Rat and Liver are random effects In aov testing for sig of 'Means' of Random Effects is pointless and that is why 'p' values are not given.Further more the interaction between a Random Effect and a Fixed Effect is also a Random Effect. The 'aov' with error structure terms output reflects this by only giving 'p' values to Fixed Effects and their interactions model - aov(Glycogen~Treatment+Error(Rat/Liver)) summary(model) Error: Rat Df Sum Sq Mean Sq F value Pr(F) #Rat is random effect Residuals 1 413.44 413.44 Error: Rat:Liver #Rat:Liver is Random effect Df Sum Sq Mean Sq F value Pr(F) Residuals 4 164.444 41.111 Error: Within Df Sum Sq Mean Sq F valuePr(F) Treatment 2 1557.56 778.78 18.251 8.437e-06 *** #Fixed effect Residuals 28 1194.78 42.67 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 I hope that this is of help. John __ 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] stepAIC invalid scope argument
Adai, The following works.Perhaps you should define your 'upper' and 'lower' in the list as aov's, as you have done with your lo,hi and mid. John stepAIC( mid, scope=list(upper = mid , lower = lo) ) Start: AIC= -594.66 y ~ x2 + x3 Df Sum of Sq RSS AIC - x21 0.11 548.56 -596.45 - x31 0.95 549.40 -594.93 none 548.45 -594.66 Step: AIC= -596.45 y ~ x3 Adaikalavan Ramasamy wrote --- I am trying to replicate the first example from stepAIC from the MASS package with my own dataset but am running into error. If someone can point where I have gone wrong, I would appreciate it very much. Here is an example : set.seed(1) df - data.frame( x1=rnorm(1000), x2=rnorm(1000), x3=rnorm(1000) ) df$y - 0.5*df$x1 + rnorm(1000, mean=8, sd=0.5) # pairs(df); head(df) lo - aov( y ~ 1, data=df ) hi - aov( y ~ .^2, data=df ) mid - aov( y ~ x2 + x3, data=df ) Running any of the following commands stepAIC( mid, scope=list(upper = ~x1 + x2 + x3 , lower = ~1) ) stepAIC( mid, scope=list(upper = hi , lower = lo) ) addterm( mid, ~ x1 + x2 + x3 ) addterm( lo, hi ) gives the same error message : Error in eval(expr, envir, enclos) : invalid second argumen __ 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] Question about 'text'
Dan, Another tweak ! If you want the 'legend' to look pretty you can resize it by adding,say, 'cex=0.6' into the legend code; try--- legend(topleft, #inset=-1, legend = do.call(expression, L), bg='white', ncol = 2, pch=c('','','',':',':',':'), x.intersp = 0.4, title=Yay! Thank You!, cex=0.6 ) John __ 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] plot question
a.d. I refer you to ?title and its given examples. try this -- plot(rnorm(10),rnorm(10),xlab= ,ylab= ) title(xlab=year, ylab=expression(paste('M x'*10^{3},)),font=2) note that 'title()' will alos accept a list for x and y labs, for additional parameters,e.g., 'col' and 'cex' John a.d wrote--- dear list: in the following plot: plot(rnorm(10),rnorm(10),xlab=year,ylab=expression (paste('M x'*10^{3},)),font.lab=2) font.lab=2, but xlab and ylab are different. I want both labels in the same way. help? a.d. __ 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] plot question
Gabor, I thought that I had worked around the 'expression' format problem, but if the x-y labels are to be bold, then using,say, cex.lab=1.25in the title(), appears to simulate 'bold' font very well, both for the ylab maths expression and xlab text. Your solution is the rigorous one! John for example -- plot(rnorm(10),rnorm(10),xlab= ,ylab= ) title(xlab=year, ylab=expression(paste(M x*10^{3})),font.lab=1,col.lab=4, cex.lab=1.25) Gabor Grothendieck wrote --- Trying characters and expressions variously it seems that font.lab applies to character strings but not to expressions so if you want to use an expression just use bold (or whatever) explicitly on the expression. One gotcha is that bold will not work as one might have expected on numbers so they must be represented as character strings -- which is why we have used 3 rather than 3 below. plot(rnorm(10),rnorm(10),xlab=quote(bold(year)),ylab=quote(bold(Mx10^3)) ) On 7/2/05, alex diaz [EMAIL PROTECTED] wrote: dear list: in the following plot: plot(rnorm(10),rnorm(10),xlab=year,ylab=expression (paste('M x'*10^{3},)),font.lab=2) font.lab=2, but xlab and ylab are different. I want both labels in the same way. help? a.d. __ 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 Components of variance
John, In addition to 'VarCorr(nlme) and VarCorr(Matrix)', you could also try 'varcomp' function in the 'ape' package.This requires an 'lme' class of file, from the 'nlme package as input.It additionally gives the options to scale the components by normalizing them (both cumulatively and in total). I hope this helps, John John Sorkin Wrote- Could someone identify a function that I might use to perform a components of variance analysis? In addition to the variance attributable to each factor, I would also like to obtain the SE of the variances. Thank you, John John Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics Baltimore VA Medical Center GRECC and University of Maryland School of Medicine Claude Pepper OAIC University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 410-605-7119 -- NOTE NEW EMAIL ADDRESS: [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] extracting components of a list
Dimitris, wouldn't this be more precise --- sapply(jj,function(x) which(x$b[1]==4)) [[1]] [1] 1 [[2]] numeric(0) [[3]] [1] 1 John Dimitris wrote --- maybe something like this: jj - list(list(a = 1, b = 4:7), list(a = 5, b = 3:6), list(a = 10, b = 4:5)) ### jj[sapply(jj, function(x) x$b[1] == 4)] I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/16/336899 Fax: +32/16/337015 Web: http://www.med.kuleuven.ac.be/biostat/ http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm - Original Message - From: Robin Hankin [EMAIL PROTECTED] To: r-help R-help@stat.math.ethz.ch Sent: Monday, June 13, 2005 4:23 PM Subject: [R] extracting components of a list Hi how do I extract those components of a list that satisfy a certain requirement? If jj - list(list(a=1,b=4:7),list(a=5,b=3:6),list(a=10,b=4:5)) I want just the components of jj that have b[1] ==4 which in this case would be the first and third of jj, vizlist (jj[[1]],jj[[3]]). How to do this efficiently? My only idea was to loop through jj, and set unwanted components to NULL, but FAQ 7.1 warns against this. -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ 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] identify label format problem
Hi , I am attempting to format the 'identify' function labels. I can format the colour but the 'cex' parameter appears not to work for me. example-- x-1:5 y-1:5 plot(x,y) identify(x,y,cex=0.5,col=2) [1] 3 The label is coloured red but the 'cex=0.5' does not reduce the label size. Why not,, am I doing it wrong ? version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor1.0 year 2005 month04 day 18 language R Thanks for any help, John __ 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] two line plot title with expression problem
Dear R-users, I am having a problem formatting a two line plot title with the first line text and the second line an expression. I cant get the second line expression to line up with the first by starting at the LHS of the line. A simple example illustrates the situation. x-seq(-5,5,length=100) f-function(x) x^2-x-1 For a one line title the following works fine plot(x,f(x),type=l) title(expression(paste(Golden Section ,fn:(phi^2-phi-1)),sep=)) but when attempting a two line title as follows-- plot(x,f(x),type=l) title(expression(paste(Golden Section\n ,fn:(phi^2-phi-1)),sep=)) the second line containing the expression for phi is indented, to the right, starting immediately after the position of the end letter of the first line. If the second line were to be text, this would not be a problem. How can I format the title so that the expression starts immediately under the first line? Thanks for help. John. [John Wilkinson : [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
[R] Re:how to compute a mode of a distribution
On 08/24/04 13:50, Paolo Tommasini wrote: Hi my name is Paolo Tommasini does anyone know how to compute a mode ( most frequent element ) for a distribution ? try median(density(my.vector)) John John Wilknson [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
Re: [R] Re:how to compute a mode of a distribution
you are quite right ,I forgot to insert $x, I should have written median(density(my.vector)$x) that gives a median result but it is still not the proper result. density(my.vector)$y) gives the probabilites y of each x value The correst answer would be the value of x corresponding to y=max(density(my.vector)$y) see summary(density(my.vector) of course a plot(density(my.vector)) would reveal the mode easily by inspection. John - Original Message - From: Kjetil Brinchmann Halvorsen [EMAIL PROTECTED] To: John Wilkinson [EMAIL PROTECTED] Cc: Paolo Tommasini [EMAIL PROTECTED]; [EMAIL PROTECTED] Sent: Thursday, August 26, 2004 1:49 PM Subject: Re: [R] Re:how to compute a mode of a distribution John Wilkinson wrote: On 08/24/04 13:50, Paolo Tommasini wrote: Hi my name is Paolo Tommasini does anyone know how to compute a mode ( most frequent element ) for a distribution ? try median(density(my.vector)) That cannot be right. Try: test - rnorm(1000) d - density(test) median(d) Error in median(d) : need numeric data str(d) List of 7 $ x: num [1:512] -3.64 -3.63 -3.62 -3.60 -3.59 ... $ y: num [1:512] 2.16e-05 2.60e-05 3.11e-05 3.71e-05 4.43e-05 ... $ bw : num 0.221 $ n: int 1000 $ call : language density(x = test) $ data.name: chr test $ has.na : logi FALSE - attr(*, class)= chr density i - which.max(d$y) d$x[i] [1] 0.05625893 Kjetil halvorsen John John Wilknson [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 __ [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] Unloading packages
In order to update the packages from Cran, it is first necessary to unload all packages from the R Console. Can anyone, please, inform me, how to unload all loaded packages with one command? Thanks for any help, John Wilkinson __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help