Re: [R] searching for data.frame rows / processing of rows
Grouping the data frame by the first two columns, apply colMeans and then rbind the resulting by-structure together: do.call(rbind, by(DF, DF[2:1], colMeans, na.rm = TRUE)) On 10/5/06, Greg Tarpinian [EMAIL PROTECTED] wrote: R 2.3.1, WinXP: I have a puzzling problem that I suspect may be solved using grep or a regular expression but am at a loss how to actually do it... My data.frame looks like Location TimeXY --- --- 1 0 1.6 9.3 1 3 4.2 10.4 1 6 2.7 16.3 2 0 0.5 2.1 2 3 NA 3.6 2 3 5.0 0.06 2 6 3.4 14.0 and so forth. I would like to search for duplicate Time values within a Location and take the numerical average (where possible) of the elements in X and Y. These numerical averages should then be used to create a single row where multiple rows once existed. So, I would like to obtain 2 3 5.0 1.83 for the two Time = 3 rows for Location = 2 and use it to replace these two rows. Ideally, avoiding for(i in 1:blah) loops would be nice because the data.frame has about 10,000 rows that need to be searched and processed. My intent is to do some comparing of SAS to R -- the DATA step processing in SAS is quite fast and using the RETAIN statement along with the LAG( ) function allows this sort of thing to be done rapidly. Thanks in advance, Greg __ 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. __ 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] help with script
On 10/4/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: Also see package caTools. On 10/4/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: See: http://tolstoy.newcastle.edu.au/R/help/04/10/5161.html On 10/4/06, JOHN VOIKLIS [EMAIL PROTECTED] wrote: Hello, I wrote the function, below, in the hope of _quickly_ generating a sliding-window time series of the mean, sd, median, and mad values from a matrix of data. The script below is anything but quick; it has been working away on a 2500 x 2500 matrix with a sliding window of 100 x 100 for five days, with no end in sight. Obviously, I am not the best of programmers. Can anyone offer suggestions as to how I might optimize this script. Thank you, John common.ground-function(inMatrix,window){ cleanMatrix-as.matrix(inMatrix) diag(cleanMatrix)-NA outMatrix-matrix(0,dim(inMatrix)[1]-window,4) colnames(outMatrix)-c(mean,SD, median,mad) for(i in 1:dim(outMatrix)[1]){ for(j in window:dim(cleanMatrix)[2]){ outMatrix[i,1]-mean(cleanMatrix[c(i:j),c(i:j)],na.rm=TRUE) outMatrix[i,2]-sd(c(cleanMatrix[c(i:j),c(i:j)]),na.rm=TRUE) outMatrix[i,3]-median(cleanMatrix[c(i:j),c(i:j)],na.rm=TRUE) outMatrix[i,4]-mad(cleanMatrix[c(i:j),c(i:j)],na.rm=TRUE) } } return(outMatrix) } Also you could look at rollmax, rollmean and rollmedian in the zoo package. rollapply in zoo with FUN = sd could be used for the sd. __ 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.
[R] A statement over multiple lines (i.e. the ... feature in Matlab)
Hello again list, I thought I'd start a new thread, seeing as it's completely different from my previous question. Some functions I have written require many parameters, and so do not fit nicely into an 80 column width display. This is usually avoided, by spreading that particular statement over a few lines. This is something that I do in Matlab with the following: myFunc( parameter1, ... parameter2, ... parameter3, ... parameter4) The ... operator in Matlab allows me to spread a statement over several lines. The ... operator in R seems to be more like the ... operator in C, which allows for a variable argument list. How do I accomplish this task in R? It's not a show stopper, but it would make reading my code much much easier. Cheers, Wee-Jin __ 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.
[R] Variables in RODBC environment
Hello Experts, how can I use variables in the RODBC environment. Example which does not work: Thanks for your help. Thorsten pn - '39R5238'; library(RODBC); odbcobj - odbcConnect(SQUIT21C,uid=muehge,pwd=xxx); sql - select u.unitid, from test where part in ('pn') ; parameter - sqlQuery(odbcobj,sql); odbcClose(odbcobj); __ 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] A statement over multiple lines (i.e. the ... feature in Matlab)
Hi Wee-Jin you can block out bits of R code with if(FALSE){ code not executed } For the line breaking, R deals with incomplete lines by not executing the statement until you finish it. In the function case, it waits for you to close a bracket. If you type: myFunc(a=3, b=5, c=6, d=7 ) myFunc() will only execute when you close the bracket HTH rksh On 5 Oct 2006, at 08:28, Wee-Jin Goh wrote: Hello again list, I thought I'd start a new thread, seeing as it's completely different from my previous question. Some functions I have written require many parameters, and so do not fit nicely into an 80 column width display. This is usually avoided, by spreading that particular statement over a few lines. This is something that I do in Matlab with the following: myFunc( parameter1, ... parameter2, ... parameter3, ... parameter4) The ... operator in Matlab allows me to spread a statement over several lines. The ... operator in R seems to be more like the ... operator in C, which allows for a variable argument list. How do I accomplish this task in R? It's not a show stopper, but it would make reading my code much much easier. Cheers, Wee-Jin __ 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. -- 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] R Graphics: Saving PDF and other formats from Windows Graphic Device for LaTeX
pdf(sgr6100.pdf, horizontal=FALSE, onefile=FALSE, height=3, width=3, pointsize=6) Reducing point-size below 6 does not seem to make any difference to the size of text and symbols. Any suggestions to get smaller font sizes? I have never used the pointsize option. Increasing the height and the width decreases the relative size of the font so if you scale down the figure in LaTeX you would get smaller fonts. I am using WinEdt with MikTeX set-up. as well do I ... Latex complains about the [scale...] part for any scale in \includegraphics[scale=1]{} A scale of one does not make any sense since it is 100% or with other words you could leave that out. If you want to scale the figure to the text width you could write \begin{figure} \includegraphic[width=\textwidth]{sgr6100.pdf } \end{figure} btw. I use the postscript device to produce an eps file and dvipdfm to convert the dvi to pdf postscript(sgr6100.eps, width = 6, height = 6, horizontal = FALSE, onefile =T, paper=special) and then in LaTeX: \begin{figure} \includegraphic[width=\textwidth]{sgr6100} \end{figure} and pauses. Suggestions will be appreciated about what is the best way to scale R graphics for inclusion in LaTeX. see above. Have a look at a good LaTeX guide like lshort http://tobi.oetiker.ch/lshort/lshort.pdf (e.g. p.73/73) Stefan Grosse __ 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] A statement over multiple lines (i.e. the ... feature in Matlab)
On 5 Oct 2006, at 08:39, Robin Hankin wrote: Hi Wee-Jin you can block out bits of R code with if(FALSE){ code not executed } For the line breaking, R deals with incomplete lines by not executing the statement until you finish it. In the function case, it waits for you to close a bracket. If you type: myFunc(a=3, b=5, c=6, d=7 ) myFunc() will only execute when you close the bracket HTH rksh Great stuff. That's exactly what I was looking for. Thanks a bunch :) Wee-Jin __ 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] A statement over multiple lines (i.e. the ... feature in
On 05-Oct-06 Wee-Jin Goh wrote: Hello again list, I thought I'd start a new thread, seeing as it's completely different from my previous question. Some functions I have written require many parameters, and so do not fit nicely into an 80 column width display. This is usually avoided, by spreading that particular statement over a few lines. This is something that I do in Matlab with the following: myFunc( parameter1, ... parameter2, ... parameter3, ... parameter4) The ... operator in Matlab allows me to spread a statement over several lines. The ... operator in R seems to be more like the ... operator in C, which allows for a variable argument list. How do I accomplish this task in R? It's not a show stopper, but it would make reading my code much much easier. Cheers, Wee-Jin Just spread it over several lines, without the ..., which in R you do not need for this purpose: myFunc(parameter1, parameter2, parameter3, parameter4) The point is that R uses its syntactic rules to determine when a component of an expression is complete. Therefore (in this instance) the opening ( must at some stage be matched by a ), so, whether or not you interpolate newlines, so R will continue to expect you to input more items until ) is encountered. Since you can (at any stage of input) have incomplete expressions within incomplete expressions, the same principle applies until every thing is closed up, to the outermost level. Since R *does* use ... for a variable argument list, you can of course use this as well, for that purpose. The important thing to remember, though, is that an expression which is incomplete from the point of view of your intentions may look complete to R. So, if you want to break it in the middle, make sure you do so at a point where it is incomplete from R's point of view. For example, suppose you intend x - 3*4*5*6*7*8 If, in R, you put x - 2*3*4* 5*6*7*8 in your input file, it will work: x - 2*3*4* + 5*6*7*8 x [1] 40320 because, syntactically, the * at the end of the first line is expected to be followed by something, so R waits for that. But if you put it in as x - 2*3*4 *5*6*7*8 it will not, because the first line is already complete when the newline is encountered, and then the second line will generate a syntax error because the initial * is required to be preceded by something. Having come to the end of a complete valid expression, R is now waiting for the start of a new valid expression; and an expression which starts with a * is not valid. Hence: x - 2*3*4 *5*6*7*8 Error: syntax error x [1] 24 Best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 05-Oct-06 Time: 09:18:51 -- XFMail -- __ 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] A statement over multiple lines (i.e. the ... feature in Matlab)
Robin Hankin wrote: For the line breaking, R deals with incomplete lines by not executing the statement until you finish it. Beware, however, that syntactically valid lines do get executed immediately (at least at the prompt). So 1 + 2 - 3 will be interpreted as two commands (returning 3 and -3, respectively), while 1 + 2 - 3 will be interpreted as a single command (returnig 0). -- Bjørn-Helge Mevik __ 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] A statement over multiple lines (i.e. the ... feature in Matlab)
On 5 Oct 2006, at 09:34, Bjørn-Helge Mevik wrote: Robin Hankin wrote: For the line breaking, R deals with incomplete lines by not executing the statement until you finish it. Beware, however, that syntactically valid lines do get executed immediately (at least at the prompt). yes, you're quite right. Your observation leads to my number 1 gotcha: if(12) print(12) else print(12) which (properly) returns an error because the first line is syntactically complete, so the else statement is processed ab initio. So now I always always always always type if(12){ print(12) } else { print(12) } if by some chance I ignore emacs/ESS's indentation telling me I've forgotten a brace. best wishes rksh So 1 + 2 - 3 will be interpreted as two commands (returning 3 and -3, respectively), while 1 + 2 - 3 will be interpreted as a single command (returnig 0). -- Bjørn-Helge Mevik __ 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. -- 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Block comments in R?
Wee-Jin Goh wrote: Hello list, Is there any way to perform a block comment in R? In C++, anything in between a /* and */ is considered a comment, and it allows programmers to comment out chunks of code for testing and debugging. Is there such a feature in R? This has frequently been asked on the list. Try to google for it. You will find answers like if(FALSE){ code block commented } or use a good editor that supports commenting and uncommenting blocks. Uwe Ligges Cheers, Wee-Jin __ 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. __ 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] compiling rgdal package on windows / macos
On Wed, 4 Oct 2006, Dylan Beaudette wrote: Greetings: As I am not a windows user, I cannot try this: is it possible to install rgdal on windows without having to compile it from source ? Andy Jaworski already replied that rgdal Windows binaries are available from CRAN mirrors, thanks to Uwe Ligges. Compilation on MacOS is within my abilities, however each time i try and install the rgdal package it dies complaining that it cannot find gdal-config --- which was recently installed with GRASS. I have updated my PATH environment variable, logged out, but R still cannot find the gdal-config program. For Mac OSX, please consult the detailed instructions on: http://www.r-project.org/Rgeo - Maps (in the left navigation bar). Hint: see the configure.args= argument of install.packages() and use --with-gdal-config= to get the address right. No Mac OSX binaries will be provided from CRAN, there are simply too many Mac OSX varieties out there to be sure that the installed external software harmonises with the rgdal binary build. However, one way of getting there is referenced in Rgeo, using William Kyngesburye's frameworks (rgdal binaries are provided in harmony with PROJ.4 and GDAL). By the way: RSiteSearch(rgdal macosx) first hit takes you to a recent reply on this topic on the R-sig-geo list. any tips on getting the rgdal package up and running on MacOS or Windows would be greatly appreciated. Cheers, -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Block comments in R?
On 5 Oct 2006, at 10:05, Uwe Ligges wrote: Wee-Jin Goh wrote: Hello list, Is there any way to perform a block comment in R? In C++, anything in between a /* and */ is considered a comment, and it allows programmers to comment out chunks of code for testing and debugging. Is there such a feature in R? This has frequently been asked on the list. Try to google for it. You will find answers like if(FALSE){ code block commented } That method doesn't work for me: if(FALSE){ if(12) print(12) else print(12) } returns an error. How would I comment out that block of (incorrect) code? or use a good editor that supports commenting and uncommenting blocks. Uwe Ligges -- 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] Block comments in R?
Robin Hankin wrote: On 5 Oct 2006, at 10:05, Uwe Ligges wrote: Wee-Jin Goh wrote: Hello list, Is there any way to perform a block comment in R? In C++, anything in between a /* and */ is considered a comment, and it allows programmers to comment out chunks of code for testing and debugging. Is there such a feature in R? This has frequently been asked on the list. Try to google for it. You will find answers like if(FALSE){ code block commented } That method doesn't work for me: if(FALSE){ if(12) print(12) else print(12) } Use an editor that comments out a whole block which is what I do all the time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of them. Or you can go ahead and use if(FALSE){ if(12) print(12) else print(12) } or if(FALSE){' if(12) print(12) else print(12) '} Uwe returns an error. How would I comment out that block of (incorrect) code? or use a good editor that supports commenting and uncommenting blocks. Uwe Ligges -- 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] searching for data.frame rows / processing of rows
I tought that aggregate was the way to go, but only for large dataframes it is faster. df - read.table(stdin(),header=TRUE) 0: Location TimeXY 1: 1 0 1.6 9.3 2: 1 3 4.2 10.4 3: 1 6 2.7 16.3 4: 2 0 0.5 2.1 5: 2 3 NA 3.6 6: 2 3 5.0 0.06 7: 2 6 3.4 14.0 8: aggregate(df[,3:4],df[,1:2],FUN=mean,na.rm=TRUE) Location Time X Y 110 1.6 9.30 220 0.5 2.10 313 4.2 10.40 423 5.0 1.83 516 2.7 16.30 626 3.4 14.00 system.time( aggregate(df[,3:4],df[,1:2],FUN=mean,na.rm=TRUE) ) [1] 0.008 0.000 0.008 0.000 0.000 system.time( do.call(rbind, by(df, df[2:1], colMeans, na.rm = TRUE))) [1] 0.005 0.000 0.005 0.000 0.000 df - data.frame(Location=rep(1:50,50), + Time=sample(rep(1:50,each=10),2500,replace=TRUE), + X=runif(2500),Y=runif(2500)) system.time( aggregate(df[,3:4],df[,1:2],FUN=mean,na.rm=TRUE) ) [1] 0.162 0.000 0.163 0.000 0.000 system.time( do.call(rbind, by(df, df[2:1], colMeans, na.rm = TRUE))) [1] 2.179 0.006 2.216 0.000 0.000 Kees Gabor Grothendieck wrote: Grouping the data frame by the first two columns, apply colMeans and then rbind the resulting by-structure together: do.call(rbind, by(DF, DF[2:1], colMeans, na.rm = TRUE)) On 10/5/06, Greg Tarpinian [EMAIL PROTECTED] wrote: R 2.3.1, WinXP: I have a puzzling problem that I suspect may be solved using grep or a regular expression but am at a loss how to actually do it... My data.frame looks like Location TimeXY --- --- 1 0 1.6 9.3 1 3 4.2 10.4 1 6 2.7 16.3 2 0 0.5 2.1 2 3 NA 3.6 2 3 5.0 0.06 2 6 3.4 14.0 and so forth. I would like to search for duplicate Time values within a Location and take the numerical average (where possible) of the elements in X and Y. These numerical averages should then be used to create a single row where multiple rows once existed. So, I would like to obtain 2 3 5.0 1.83 for the two Time = 3 rows for Location = 2 and use it to replace these two rows. Ideally, avoiding for(i in 1:blah) loops would be nice because the data.frame has about 10,000 rows that need to be searched and processed. My intent is to do some comparing of SAS to R -- the DATA step processing in SAS is quite fast and using the RETAIN statement along with the LAG( ) function allows this sort of thing to be done rapidly. Thanks in advance, Greg __ 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. __ 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. __ 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] Variables in RODBC environment
hi, use paste: # a string/number/date pn - '39R5238'; sql - paste(select u.unitid from test where part =, pn) ^^ # an array of strings/numbers/dates pn=c(1,2,3,4) sql - paste(select u.unitid from test where part in (, paste(pn, collapse=','), )) ^^ regards soren obs the use of '=' and in - Original Message - From: Thorsten Muehge [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Thursday, October 05, 2006 9:35 AM Subject: [R] Variables in RODBC environment Hello Experts, how can I use variables in the RODBC environment. Example which does not work: Thanks for your help. Thorsten pn - '39R5238'; library(RODBC); odbcobj - odbcConnect(SQUIT21C,uid=muehge,pwd=xxx); sql - select u.unitid, from test where part in ('pn') ; parameter - sqlQuery(odbcobj,sql); odbcClose(odbcobj); __ 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. __ 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.
[R] VGAM Package ?
Hi! R users I would like to ask you where could we find the VGAM Package. I don't find it in the list of packages. Thak you for your help Lassana KOITA Etudes de Sécurité et d'Exploitation aéroportuaires / Safety Study Statistical analysis Service Technique de l'Aviation Civile (STAC) / Civil Aviation Technical Department Direction Générale de l'Aviation Civile (DGAC) / French Civil Aviation Authority Tel: 01 49 56 80 60 Fax: 01 49 56 82 14 E-mail: [EMAIL PROTECTED] http://www.stac.aviation-civile.gouv.fr/ __ 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] VGAM Package ?
google VGAM On Thu, 5 Oct 2006, KOITA Lassana - STAC/ACE wrote: Hi! R users I would like to ask you where could we find the VGAM Package. I don't find it in the list of packages. Thak you for your help Lassana KOITA Etudes de Sécurité et d'Exploitation aéroportuaires / Safety Study Statistical analysis Service Technique de l'Aviation Civile (STAC) / Civil Aviation Technical Department Direction Générale de l'Aviation Civile (DGAC) / French Civil Aviation Authority Tel: 01 49 56 80 60 Fax: 01 49 56 82 14 E-mail: [EMAIL PROTECTED] http://www.stac.aviation-civile.gouv.fr/ __ 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. _ David Scott Visiting (Until January 07) Department of Probability and Statistics The University of Sheffield The Hicks Building Hounsfield Road Sheffield S3 7RH United Kingdom Phone: +44 114 222 3908 Email: [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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] (not so real) function (more like #include)
Thanks a lot! longfun2 longfun3 work perfect for my very untypical problem. As I have many local variables, usual functions with parameters are very uncomfortable, but the code you gave me is great! Meinhard On Oct 4, 2006, at 5:25 PM, Gabor Grothendieck wrote: longfun could just pass a, b and d to each of the individual functions and each of the individual functions could pass out back as a return value. f1 - f2 - function(a, b, d) a+b+d longfun1 - function() { a - b - d - 1 out - f1(a, b, d) out - f2(a, b, d) + out out } longfun1() # 6 If the problem is that a, b and d actually represent 100 variables, say, and its a nuisance to pass them all explicitly you could do this: f1 - f2 - function() with(parent.frame(), a + b + d) longfun2 - function() { a - b - d - 1 out - f1() out - f2() + out out } longfun2() # 6 or you could do it this way if you would prefer to have longfun control the scoping rather than having it done within each individual function: f1 - f2 - function() a+b+d longfun3 - function() { a - b - d - 1 environment(f1) - environment(f2) - environment() out - f1() f2() + out } longfun3() # 6 On 10/4/06, Meinhard Ploner [EMAIL PROTECTED] wrote: Hello R users, I have a very long function with parts of that looking like a list of jobs. The first part of the function defines a lot of local variables, which are used by the jobs. Now I extracted the job part und putted them into an external file, used by source(, local=T), which is more comfortable to me. Apart from that it gives me the opportunity that more users can work on the project. Is it possible to define a function for that code and passing to it the environment of f() without save(list=ls()) and load() ? Thanks in advance Meinhard Ploner longfun - f() { ## lot of local variables: a - ... b - ... d - ... ... out - ... source(job1.txt, local=TRUE) #changes out, uses a, b, d, ... source(job2.txt, local=TRUE) # changes out, uses a, b, d, ... ... } version _ platform i386-apple-darwin8.6.1 arch i386 os darwin8.6.1 system i386, darwin8.6.1 status major 2 minor 3.1 year 2006 month 06 day01 svn rev38247 language R version.string Version 2.3.1 (2006-06-01) __ 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. __ 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] Variables in RODBC environment
how can I use variables in the RODBC environment. Example which does not work: Thanks for your help. Thorsten pn - '39R5238'; library(RODBC); odbcobj - odbcConnect(SQUIT21C,uid=muehge,pwd=xxx); sql - select u.unitid, from test where part in ('pn') ; parameter - sqlQuery(odbcobj,sql); odbcClose(odbcobj); You can compose your query simply using paste. In your example it would be pn - '39R5238'; sql - paste(select u.unitid, from test where part in (', pn, ')); or, to avoid problems with the newline character: sql - paste(select, u.unitid,, from test, where part in (', pn, ')); In this case you'd have: sql [1] select u.unitid, from test where part in (' 39R5238 ') It's clear that, if the spaces in (' 39R5238 ') are a problem for you , you can use sql - paste(select , u.unitid, , from test , where part in (', pn, '), sep = ); which lends to sql [1] select u.unitid, from test where part in ('39R5238') HTH Sandro Bonfigli __ 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] VGAM Package ?
On Thu, 2006-10-05 at 11:39 +0200, KOITA Lassana - STAC/ACE wrote: Hi! R users I would like to ask you where could we find the VGAM Package. I don't find it in the list of packages. Thak you for your help Lassana KOITA That's because it isn't on CRAN yet. You can find a link to VGAM on the Environmetrics Task View: http://www.stats.bris.ac.uk/R/src/contrib/Views/Environmetrics.html HTH G -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Gavin Simpson [t] +44 (0)20 7679 0522 ECRC ENSIS, UCL Geography, [f] +44 (0)20 7679 0565 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/ UK. WC1E 6BT. [w] http://www.freshwaters.org.uk %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ 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.
[R] format of data for use in clim.pact
Goodmorning everyone! Just a quick question: I want to apply eof analysis and downscaling using the clim.pact package. What form should my data have? Only netCDF? Right now my data are of the form of ascii files for each month with 3 columns corresponding to latitude, longitude and value of my observations. Any idea on how to handle them would be welcome. Thanks a lot Isidora __ 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.
[R] glm with nesting
I just had a manuscript returned with the biggest problem being the analysis. Instead of using principal components in a regression I've been asked to analyze a few variables separately. So that's what I'm doing. I pulled a feather from young birds and we quantified certain aspects of the color of those feathers. Since I often have more than one sample from a nest, I thought I should use a nested design. Here's the code I've been using and I'd appreciate if someone could look it over and see if it was correct. bb.glm1 - glm(rtot ~ box/(julian +purbank), data=bbmale, family=gaussian, na.action=na.omit) where rtot = total reflectance, box = nest box (i.e., birdhouse), julian = day of the year and purbank = the proportion of urban cover in a 1 km buffer around the nest box. I'm not interested in the box effect and I've seperated males and female chicks. I've asked about nestedness before and I was given code that included | to indicate nestedness but this indicates a grouping does it not? I suspect that there is something wrong. In the summary I get Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 2.880e-01 3.224e-03 89.322 2e-16 *** box -3.219e-05 6.792e-05 -0.4740.636 box:julian 7.093e-08 3.971e-07 0.1790.859 box:purbank -1.735e-05 1.502e-04 -0.1150.908 The other question I have is how do I test a null hypothesis - no explanatory variables? [rtot ~ NULL?] Many thanks, Jeff Jeffrey A. Stratford, Ph.D. Postdoctoral Associate 331 Funchess Hall Department of Biological Sciences Auburn University Auburn, AL 36849 334-329-9198 FAX 334-844-9234 http://www.auburn.edu/~stratja __ 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] glm with nesting
Jeffrey Stratford [EMAIL PROTECTED] writes: I just had a manuscript returned with the biggest problem being the analysis. Instead of using principal components in a regression I've been asked to analyze a few variables separately. So that's what I'm doing. I pulled a feather from young birds and we quantified certain aspects of the color of those feathers. Since I often have more than one sample from a nest, I thought I should use a nested design. Notwithstanding comments below, that quote could be aiming for the fortunes package... Here's the code I've been using and I'd appreciate if someone could look it over and see if it was correct. bb.glm1 - glm(rtot ~ box/(julian +purbank), data=bbmale, family=gaussian, na.action=na.omit) where rtot = total reflectance, box = nest box (i.e., birdhouse), julian = day of the year and purbank = the proportion of urban cover in a 1 km buffer around the nest box. I'm not interested in the box effect and I've seperated males and female chicks. I've asked about nestedness before and I was given code that included | to indicate nestedness but this indicates a grouping does it not? I suspect that there is something wrong. In the summary I get Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 2.880e-01 3.224e-03 89.322 2e-16 *** box -3.219e-05 6.792e-05 -0.4740.636 box:julian 7.093e-08 3.971e-07 0.1790.859 box:purbank -1.735e-05 1.502e-04 -0.1150.908 Several things look wrong here. Most importantly, you appear to have single-degree of freedom effects (t tests) of things that appear not to be linear effects: Certainly, you have more than two nest boxes, but also day of year as a linear term looks suspicious to me. Unless there is something I have missed completely, box should be a factor variable, and you might also need trigonometric terms for the julian effect (depending on what sort of time spans we are talking about.) Secondly, notation like box/julian suggests that julian only makes sense within a nest box i.e. 1st of March in one box is completely different from 1st of March in another box (the notation is more commonly used to describe bird number within nests and the like). And with purbank presumably constant for measurements from the same box, the box:purbank term looks strange indeed. If you want to take account of a between-box variation in the effect of covariates, you probably need to add them as variance components, but this requires non-glm software, either lme() or lmer(). However, instructing you on those is outside the scope of this mailing list, and you may need to find a local consultant. The other question I have is how do I test a null hypothesis - no explanatory variables? [rtot ~ NULL?] Many thanks, Jeff Jeffrey A. Stratford, Ph.D. Postdoctoral Associate 331 Funchess Hall Department of Biological Sciences Auburn University Auburn, AL 36849 334-329-9198 FAX 334-844-9234 http://www.auburn.edu/~stratja __ 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. -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ 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.
[R] [Fwd: Re: Block comments in R?]
Ooops! Sorry, I send it only to Uwe Ligges the first time. Best, Philippe Grosjean This is perhaps another solution, more elegant in the way the block comment is written... but it requires to redefine `!` and slows it a little bit because it tests first its arguments before calling .Primitive(!): It takes advantage of `!` being not defined for character arguments: !2 [1] FALSE !some text Error in !some text : invalid argument type So, now, we will define it for character arguments (invisibly returns the argument) `!`- function(x) + if (inherits(x, character) == FALSE) + .Primitive(!)(x) else invisible(x) Now, `!` can be used to construct a block comment: A R script with block comments = 1+1 # This is a line comment !' This is a block comment spread on several lines... ' ls() !!' This is another block comment, possibly of higher importance than the previous one ' search() !!!' For color syntax highlighting and to better detect the end of block comments, one may also decide to use the same code for opening and closing the comments like it is the present case !!!' !' Note that the only constraint is to escape single quotes in block comments (or not to use single quotes) Of course, one could also decide to use double quotes instead of single quotes !' !' Now, it would be nice to have a little patch of .Primitive(!) that simply displays no error message in case the argument of `!`is a character sting. So, the hock would not be required any more !' And of the R script = Best, Philippe Grosjean Uwe Ligges wrote: Robin Hankin wrote: On 5 Oct 2006, at 10:05, Uwe Ligges wrote: Wee-Jin Goh wrote: Hello list, Is there any way to perform a block comment in R? In C++, anything in between a /* and */ is considered a comment, and it allows programmers to comment out chunks of code for testing and debugging. Is there such a feature in R? This has frequently been asked on the list. Try to google for it. You will find answers like if(FALSE){ code block commented } That method doesn't work for me: if(FALSE){ if(12) print(12) else print(12) } Use an editor that comments out a whole block which is what I do all the time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of them. Or you can go ahead and use if(FALSE){ if(12) print(12) else print(12) } or if(FALSE){' if(12) print(12) else print(12) '} Uwe returns an error. How would I comment out that block of (incorrect) code? or use a good editor that supports commenting and uncommenting blocks. Uwe Ligges -- 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 and provide commented, minimal, self-contained, reproducible code. __ 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.
[R] Partition into quantiles
Is there any function that divides a sample into N quantiles? For example, for N = 2, this would be the solution: x - rnorm(100) m - median(x) q - ifelse(x = median, 1, 2) Alberto Monteiro __ 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] Partition into quantiles
Alberto Monteiro [EMAIL PROTECTED] writes: Is there any function that divides a sample into N quantiles? For example, for N = 2, this would be the solution: x - rnorm(100) m - median(x) q - ifelse(x = median, 1, 2) Have a look at N - 2 table(cut(x,quantile(x,seq(0,1,1/N)), include.lowest=TRUE)) [-2.78,0.205] (0.205,2.22] 5050 table(cut(x,c(-Inf,quantile(x,(1:(N-1))/N),Inf))) (-Inf,0.205] (0.205, Inf] 50 50 -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ 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] Partition into quantiles
cut.quantile - function(x, N, use.ppoints=TRUE) { qq - if (use.ppoints) c(0,ppoints(N-1),1) else seq(0, 1, 1/N) breaks - quantile(x, qq) breaks[1] - breaks[1] - 1 breaks[N+1] - breaks[N+1]+1 cut(x, breaks) } tmpT - cut.quantile(rnorm(100), 10, TRUE) table(tmpT) tmpF - cut.quantile(rnorm(100), 10, FALSE) table(tmpF) __ 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] Partition into quantiles
You might also want to look at the function quantcut in the gtools package (part of the gregmisc bundle). On 05/10/06, Alberto Monteiro [EMAIL PROTECTED] wrote: Is there any function that divides a sample into N quantiles? For example, for N = 2, this would be the solution: x - rnorm(100) m - median(x) q - ifelse(x = median, 1, 2) Alberto Monteiro __ 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. -- = David Barron Said Business School University of Oxford Park End Street Oxford OX1 1HP __ 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] [Fwd: Re: Block comments in R?]
Philippe Grosjean wrote: It takes advantage of `!` being not defined for character arguments: *gasp* how much R code is destined to feature on www.thedailywtf.com in the future? whats the chances of block commenting being included in a future version? and more generally, is there a feature request system, or discussion of whats going in new releases? the nearest I can find is the 'ideas' file: http://developer.r-project.org/ideas.txt which has a bit of a 1990's feel to it (okay, it is in the 'Older Material' section). Barry __ 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] glm with nesting
Peter and list, Thanks for the response. A did add box as a factor (box - factor(box)). Julian should be linear - bluebird chicks are bluer as the season progresses from March to August. I did try the following rtot.lme - lmer(rtot ~ sex +(purban|box:chick) + (purban|box), data=bb, na.action=na.omit) # from H. Doran and rtot.lme2 - lme(fixed=rtot ~ sex + purban + sexv:purban, data = bb, random = ~1 |box) # from K. Jones [EMAIL PROTECTED] but these did not work (months ago and I don't remember exactly why) and I have since seperated males and females and added day of the year (julian). But | does indicate grouping not nested, correct? Could someone suggest some coding that might work? Thanks again, Jeff Peter Dalgaard [EMAIL PROTECTED] 10/05/06 7:14 AM Jeffrey Stratford [EMAIL PROTECTED] writes: I just had a manuscript returned with the biggest problem being the analysis. Instead of using principal components in a regression I've been asked to analyze a few variables separately. So that's what I'm doing. I pulled a feather from young birds and we quantified certain aspects of the color of those feathers. Since I often have more than one sample from a nest, I thought I should use a nested design. Notwithstanding comments below, that quote could be aiming for the fortunes package... Here's the code I've been using and I'd appreciate if someone could look it over and see if it was correct. bb.glm1 - glm(rtot ~ box/(julian +purbank), data=bbmale, family=gaussian, na.action=na.omit) where rtot = total reflectance, box = nest box (i.e., birdhouse), julian = day of the year and purbank = the proportion of urban cover in a 1 km buffer around the nest box. I'm not interested in the box effect and I've seperated males and female chicks. I've asked about nestedness before and I was given code that included | to indicate nestedness but this indicates a grouping does it not? I suspect that there is something wrong. In the summary I get Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 2.880e-01 3.224e-03 89.322 2e-16 *** box -3.219e-05 6.792e-05 -0.4740.636 box:julian 7.093e-08 3.971e-07 0.1790.859 box:purbank -1.735e-05 1.502e-04 -0.1150.908 Several things look wrong here. Most importantly, you appear to have single-degree of freedom effects (t tests) of things that appear not to be linear effects: Certainly, you have more than two nest boxes, but also day of year as a linear term looks suspicious to me. Unless there is something I have missed completely, box should be a factor variable, and you might also need trigonometric terms for the julian effect (depending on what sort of time spans we are talking about.) Secondly, notation like box/julian suggests that julian only makes sense within a nest box i.e. 1st of March in one box is completely different from 1st of March in another box (the notation is more commonly used to describe bird number within nests and the like). And with purbank presumably constant for measurements from the same box, the box:purbank term looks strange indeed. If you want to take account of a between-box variation in the effect of covariates, you probably need to add them as variance components, but this requires non-glm software, either lme() or lmer(). However, instructing you on those is outside the scope of this mailing list, and you may need to find a local consultant. The other question I have is how do I test a null hypothesis - no explanatory variables? [rtot ~ NULL?] Many thanks, Jeff Jeffrey A. Stratford, Ph.D. Postdoctoral Associate 331 Funchess Hall Department of Biological Sciences Auburn University Auburn, AL 36849 334-329-9198 FAX 334-844-9234 http://www.auburn.edu/~stratja 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. -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ 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] [Fwd: Re: Block comments in R?]
On 2006-10-5 9:20, Barry Rowlingson wrote: Philippe Grosjean wrote: It takes advantage of `!` being not defined for character arguments: *gasp* how much R code is destined to feature on www.thedailywtf.com in the future? whats the chances of block commenting being included in a future version? and more generally, is there a feature request system, or discussion of whats going in new releases? the nearest I can find is the 'ideas' file: There's a feature request system (one of the categories of bug report is wishlist). I think a specific proposal with some work behind it would probably get action in a case like this. That is, block comments would be nice likely won't get acted on, but block comments in the form ... which is used by other language ..., which has the fantastic properties ... and no downside, and which could be implemented by the following patch to the parser,, likely would. Duncan Murdoch http://developer.r-project.org/ideas.txt which has a bit of a 1990's feel to it (okay, it is in the 'Older Material' section). Barry __ 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. __ 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] glm with nesting
It's not really possible to help without knowing what errors you received and maybe some reproducible code. I think I remember this, though. From what I recall, there was no distinction between box and chick, so you cannot estimate both variance components. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Jeffrey Stratford Sent: Thursday, October 05, 2006 9:27 AM To: [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Subject: Re: [R] glm with nesting Peter and list, Thanks for the response. A did add box as a factor (box - factor(box)). Julian should be linear - bluebird chicks are bluer as the season progresses from March to August. I did try the following rtot.lme - lmer(rtot ~ sex +(purban|box:chick) + (purban|box), data=bb, na.action=na.omit) # from H. Doran and rtot.lme2 - lme(fixed=rtot ~ sex + purban + sexv:purban, data = bb, random = ~1 |box) # from K. Jones [EMAIL PROTECTED] but these did not work (months ago and I don't remember exactly why) and I have since seperated males and females and added day of the year (julian). But | does indicate grouping not nested, correct? Could someone suggest some coding that might work? Thanks again, Jeff Peter Dalgaard [EMAIL PROTECTED] 10/05/06 7:14 AM Jeffrey Stratford [EMAIL PROTECTED] writes: I just had a manuscript returned with the biggest problem being the analysis. Instead of using principal components in a regression I've been asked to analyze a few variables separately. So that's what I'm doing. I pulled a feather from young birds and we quantified certain aspects of the color of those feathers. Since I often have more than one sample from a nest, I thought I should use a nested design. Notwithstanding comments below, that quote could be aiming for the fortunes package... Here's the code I've been using and I'd appreciate if someone could look it over and see if it was correct. bb.glm1 - glm(rtot ~ box/(julian +purbank), data=bbmale, family=gaussian, na.action=na.omit) where rtot = total reflectance, box = nest box (i.e., birdhouse), julian = day of the year and purbank = the proportion of urban cover in a 1 km buffer around the nest box. I'm not interested in the box effect and I've seperated males and female chicks. I've asked about nestedness before and I was given code that included | to indicate nestedness but this indicates a grouping does it not? I suspect that there is something wrong. In the summary I get Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 2.880e-01 3.224e-03 89.322 2e-16 *** box -3.219e-05 6.792e-05 -0.4740.636 box:julian 7.093e-08 3.971e-07 0.1790.859 box:purbank -1.735e-05 1.502e-04 -0.1150.908 Several things look wrong here. Most importantly, you appear to have single-degree of freedom effects (t tests) of things that appear not to be linear effects: Certainly, you have more than two nest boxes, but also day of year as a linear term looks suspicious to me. Unless there is something I have missed completely, box should be a factor variable, and you might also need trigonometric terms for the julian effect (depending on what sort of time spans we are talking about.) Secondly, notation like box/julian suggests that julian only makes sense within a nest box i.e. 1st of March in one box is completely different from 1st of March in another box (the notation is more commonly used to describe bird number within nests and the like). And with purbank presumably constant for measurements from the same box, the box:purbank term looks strange indeed. If you want to take account of a between-box variation in the effect of covariates, you probably need to add them as variance components, but this requires non-glm software, either lme() or lmer(). However, instructing you on those is outside the scope of this mailing list, and you may need to find a local consultant. The other question I have is how do I test a null hypothesis - no explanatory variables? [rtot ~ NULL?] Many thanks, Jeff Jeffrey A. Stratford, Ph.D. Postdoctoral Associate 331 Funchess Hall Department of Biological Sciences Auburn University Auburn, AL 36849 334-329-9198 FAX 334-844-9234 http://www.auburn.edu/~stratja 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. -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED])
[R] Changes in console preferences won't take
Greetings, I've attempted to change the font size directly in Rconsole as well as through the GUI preferences. Neither seems to take even when I save the preference changes. When I make the change through the gui preference, and click apply the change is made, but when I save it, close, and reopen, there is no change. Any thouughts? I'm running R.2.40. David -- === David Kaplan, Ph.D. Professor Department of Educational Psychology University of Wisconsin - Madison Educational Sciences, Room, 1061 1025 W. Johnson Street Madison, WI 53706 email: [EMAIL PROTECTED] homepage: http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm Phone: 608-262-0836 __ 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] Block comments in R?
Uwe Ligges [EMAIL PROTECTED] writes: Use an editor that comments out a whole block which is what I do all the time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of them. This, of course, works. The if(FALSE) approach does not because it requires the comment to be syntactically correct. The multiline string trick is _almost_ there, the problem is that both and ' are fairly common in text and having to escape that is a pain. My wtf feature request is to add a multiline string delimiter ala Python like . 8-P + seth __ 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] Partition into quantiles
David Barron wrote: You might also want to look at the function quantcut in the gtools package (part of the gregmisc bundle). Also, cut2 in Hmisc will do this and will label the intervals compactly. Frank On 05/10/06, Alberto Monteiro [EMAIL PROTECTED] wrote: Is there any function that divides a sample into N quantiles? For example, for N = 2, this would be the solution: x - rnorm(100) m - median(x) q - ifelse(x = median, 1, 2) Alberto Monteiro __ 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. -- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt 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 and provide commented, minimal, self-contained, reproducible code.
[R] research
dear sir iam algerian student iam 30 years old i finished my postgraduate in applied economic and statistics .i prepare my doctorat about effect of kyoto protocol on our economic. i need the document and the articals and the adress of doctor in this domain in your university i shall be looking forward and hoping that give my request forward consideration rachid toumache adress: 11 rue doudou mothtar ben aknoun algeria institut national de la planification et de la statistique inps - [[alternative HTML version deleted]] __ 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] Partition into quantiles
Sorry, folks. Thanks for the help, but none of the suggestions worked quite as I wanted :-/ So I wrote the code. It seems to work: gera_particao - function(x, n) { y - sort(x) icut - as.integer(seq(1, length(x)+1, length = n + 1)) icut - icut[c(-(n+1))] ycut - y[icut] for (i in 1:length(x)) xpart[i] - sum(x[i] = ycut) return(xpart) } First, x is sorted to y. Then, I select the minima y of each of the n segments. Then, an ugly and slow loop, counts for each x how many ycut lie below them. x - runif(12) x [1] 0.2971455266 0.6112766485 0.5571073645 0.5886481798 0.7499023860 [6] 0.1681732289 0.6319822536 0.0005354732 0.8055324992 0.8841625380 [11] 0.0726578285 0.6250309648 gera_particao(x, 2) [1] 1 2 1 1 2 1 2 1 2 2 1 2 gera_particao(x, 3) [1] 1 2 2 2 3 1 3 1 3 3 1 2 gera_particao(x, 4) [1] 2 3 2 2 4 1 3 1 4 4 1 3 gera_particao(x, 12) [1] 4 7 5 6 10 3 9 1 11 12 2 8 Alberto Monteiro __ 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] warning message in nlm
Without 'msc39', I can't say for sure, but I doubt if it's anything to worry about. It looks like 'nlm' tests values for theta and len that produce either NA or Inf in computing 'loglikcs'. Since you got an answer, it does not look to me like it's anything worth worrying about; just make sure you are minimizing (-log(likelihood)), not maximizing it. For more detail, you can run 'nlm' with print.level = 2 rather than the default of 0 -- and make contour plots of loglikcs in an appropriate region of the optimum, as suggested in Venables and Ripley (2002) Modern Applied Statistics with S (Springer), discussing 'expand.grid' and 'contour'. Hope this helps. Spencer Graves singyee ling wrote: Dear R-users, I am trying to find the MLEs for a loglikelihood function (loglikcs39) and tried using both optim and nlm. fredcs39-function(b1,b2,x){return(exp(b1+b2*x))} loglikcs39-function(theta,len){ sum(mcs39[1:len]*fredcs39(theta[1],theta[2],c(8:(7+len))) - pcs39[1:len] * log(fredcs39(theta[1],theta[2],c(8:(7+len) } theta.start-c(0.1,0.1) 1. The output from using optim is as follow -- optcs39-optim(theta.start,loglikcs39,len=120,method=BFGS) optcs39 $par [1] -1.27795226 -0.03626846 $value [1] 7470.551 $counts function gradient 133 23 $convergence [1] 0 $message NULL 2. The output from using nlm is as follow --- outcs39-nlm(loglikcs39,theta.start,len=120) Warning messages: 1: NA/Inf replaced by maximum positive value 2: NA/Inf replaced by maximum positive value 3: NA/Inf replaced by maximum positive value 4: NA/Inf replaced by maximum positive value 5: NA/Inf replaced by maximum positive value 6: NA/Inf replaced by maximum positive value 7: NA/Inf replaced by maximum positive value outcs39 $minimum [1] 7470.551 $estimate [1] -1.27817854 -0.03626027 $gradient [1] -8.933577e-06 -1.460512e-04 $code [1] 1 $iterations [1] 40 As you can see, the values obtained from using both functions are very similar. But, what puzzled is the warning message that i got from using nlm. Could anyone please shed some light on how this warning message come about and whether it is a cause for concern? Many thanks in advance for any advice! singyee [[alternative HTML version deleted]] __ 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. __ 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] Block comments in R?
Seth Falcon [EMAIL PROTECTED] writes: Uwe Ligges [EMAIL PROTECTED] writes: Use an editor that comments out a whole block which is what I do all the time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of them. This, of course, works. The if(FALSE) approach does not because it requires the comment to be syntactically correct. The multiline string trick is _almost_ there, the problem is that both and ' are fairly common in text and having to escape that is a pain. My wtf feature request is to add a multiline string delimiter ala Python like . Still gives you problems with nested comments. *IF* we want to go down that route, we should have directional symbols like dsaldfysdfk (which would likely drive version control users up the wall...) or /# kdlsfk dfkjsdfl #/ (but those character combinations have a marginal risk of being used in existing code). -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ 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.
[R] matrix multiplication
Dear all, I have 2 matrices, one is a 8x3 matrix, called X; the other matrix is a 3x3 indicator matrix with the diagonal element as 0 or 1. when a variable is included in the model, the corresponding diagonal element is 1, otherwise, it is 0. Let A be a set of matrices that contain the possible indicator matrix. e.g., A= [A1, A2, A3], where A1= [1,0,0;0,0,0,0,0,0], A2 =[1,0,0;0,1,0,0,0,0], A3 =[1,0,0;0,0,0,0,1,0] In order to derive the new X matries depending on the indicator matrices, I use for loops to multiply both X and A as following p-3 for ( i in 1:p){ XX- X%*%A[i] } However, it only shows the result when i=p. How can I derive the results which include all possible value of XX[i], i=1,..,p, rather than just i=p thanks for any help Ya-Hsiu [[alternative HTML version deleted]] __ 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] Changes in console preferences won't take
On 10/5/2006 9:41 AM, David Kaplan wrote: Greetings, I've attempted to change the font size directly in Rconsole as well as through the GUI preferences. Neither seems to take even when I save the preference changes. When I make the change through the gui preference, and click apply the change is made, but when I save it, close, and reopen, there is no change. Any thouughts? I'm running R.2.40. I would guess you've got conflicting settings in multiple Rconsole files, or you're saving the change to the wrong place. See the ?Rconsole man page for a description of how R looks for the file (and the example section of it to see what would normally be loaded on your own system). Duncan Murdoch __ 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] glm with nesting
Harold and list, I've changed a few things since the last time so I'm really starting from scratch. I start with bbmale - read.csv(c:\\eabl\\2004\\feathers\\male_feathers2.csv, header=TRUE) box -factor(box) chick - factor(chick) Here's a sample of the data box,chick,julian,cltchsz,mrtot,cuv,cblue,purbank,purban2,purban1,pgrassk,pgrass2,pgrass1,grassdist,grasspatchk 1,2,141,2,21.72290152,0.305723811,0.327178868,0.003813435,0.02684564,0.06896552,0.3282487,0.6845638,0.7586207,0,3.73 4,1,164,4,18.87699007,0.281863299,0.310935559,0.06072162,0.2080537,0.06896552,0.01936052,0,0,323.1099,0.2284615 4,2,164,4,19.64359348,0.294117388,0.316049817,0.06072162,0.2080537,0.06896552,0.01936052,0,0,323.1099,0.2284615 7,1,118,4,13.48699876,0.303649408,0.31765218,0.3807568,0.4362416,0.6896552,0.06864183,0.03355705,0,94.86833,0.468 12,1,180,4,21.42196378,0.289731361,0.317562294,0.09238011,0.1342282,0,0.2430127,0.8322148,1,0,1.199032 12,2,180,4,18.79487905,0.286052077,0.316367349,0.09238011,0.1342282,0,0.2430127,0.8322148,1,0,1.199032 12,3,180,4,12.83127682,0.260197475,0.292636914,0.09238011,0.1342282,0,0.2430127,0.8322148,1,0,1.199032 15,1,138,4,20.07161467,0.287632782,0.318671887,0.07046477,0.03355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818 15,2,138,4,17.61146256,0.305581768,0.315848051,0.07046477,0.03355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818 15,3,138,4,20.36397134,0.271795667,0.30539683,0.07046477,0.03355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818 15,4,138,4,20.81940158,0.269468041,0.304160648,0.07046477,0.03355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818 As you can see I have multiple boxes ( 70). Sometimes I have multiple chicks per box each having their own response to mrtot, cuv, and cblue but the same landscape variables for that box. Chick number is randomly assigned and is not an effect I'm interested in. I'm really not interested in the box effect either. I would like to know if landscape affects the color of chicks (which may be tied into chick health/physiology). We also know that chicks get bluer as the season progresses and that clutch size (cltchsz) has an effect so I'm including that as covariates. Hopefully, this clears things up a bit. I do have the MASS and MEMS (Pineiro's) texts in hand. Many thanks, Jeff __ 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] Plotting text with lattice
Hi, On a related note, if I wanted to add different texts to different panels, should I stick to using trellis.focus() for each text in each panel? I cannot figure out a way to do it using a panel function. Ritwik. On 9/29/06, Deepayan Sarkar [EMAIL PROTECTED] wrote: On 9/29/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: Here are two possibilities. The first use trellis.focus/trellis/unfocus to add text subsequent to drawing the xyplot and the second uses a custom panel: xyplot(x ~ x, data = data.frame(x = 1:10)) trellis.focus(panel, 1, 1) panel.text(x=2, y=4, labels=Text) trellis.unfocus() xyplot(x ~ x, data = data.frame(x = 1:10), panel = function(...) { panel.xyplot(...) panel.text(x=2, y=4, labels=Text) }) Right, on a more general note, this is necessary because it is not clear what ltext(x=2, y=4, labels=Text) should do for a multi-panel plot. On an even more general note, you will keep getting in trouble when using lattice if you (1) try to follow the incremental addition approach of standard graphics or (2) use the par() system in any way. Deepayan On 9/29/06, McGehee, Robert [EMAIL PROTECTED] wrote: Hello, I've decided to take the leap and try my hand at the lattice package, though I am getting stuck at what one might consider a trivial problem, plotting text at a point in a graph. Apologies in advance if (that) I'm missing something extremely basic. Consider in base graphics: plot(1:10) text(2, 4, Text) In the above you will see text centered at the point (2, 4) on the graph. Now I would like to try to do the same thing using the lattice package: xyplot(x ~ x, data = data.frame(x = 1:10)) ltext(x=2, y=4, labels=Text) panel.text(x=2, y=4, labels=Text) grid.text(label=Text, x=2, y=4) grid.text(label=Text, x=unit(2, native), y=unit(4, native)) None of the above four commands puts the Text at the (2, 4) point on the graph. Any help with this would be appreciated! Also, if I have more than one panel and would like to place text at different points on different panels how would I do this? Also, note that I'm hoping to use text to label interesting points in a levelplot, but am using the above xyplot as an example. Thanks, Robert Robert McGehee Quantitative Analyst Geode Capital Management, LLC 53 State Street, 5th Floor | Boston, MA | 02109 Tel: 617/392-8396Fax:617/476-6389 mailto:[EMAIL PROTECTED] This e-mail, and any attachments hereto, are intended for us...{{dropped}} __ 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. __ 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. -- http://www.stat.wisc.edu/~deepayan/ __ 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. -- Ritwik Sinha Graduate Student Epidemiology and Biostatistics Case Western Reserve University [EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha __ 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] Block comments in R?
Still gives you problems with nested comments. *IF* we want to go down that route, we should have directional symbols like dsaldfysdfk What about heredocs? (http://en.wikipedia.org/wiki/Heredoc) Hadley __ 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] Block comments in R?
PD == Peter Dalgaard [EMAIL PROTECTED] on 05 Oct 2006 16:12:50 +0200 writes: PD Seth Falcon [EMAIL PROTECTED] writes: Uwe Ligges [EMAIL PROTECTED] writes: Use an editor that comments out a whole block which is what I do all the time, e.g. use Tinn-R, Emacs or WinEdt, to mention just a few of them. This, of course, works. The if(FALSE) approach does not because it requires the comment to be syntactically correct. The multiline string trick is _almost_ there, the problem is that both and ' are fairly common in text and having to escape that is a pain. My wtf feature request is to add a multiline string delimiter ala Python like . PD Still gives you problems with nested comments. *IF* we want to go down PD that route, we should have directional symbols like PD PD dsaldfysdfk PD (which would likely drive version control users up the wall...) i.e. all of R-core ... PD or PD /# kdlsfk PD dfkjsdfl #/ PD (but those character combinations have a marginal risk of being used PD in existing code). or ## ... ## but I'm not yet fond of the general idea. Martin PD -- PD O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B PD c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K PD (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 PD ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ 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] treatment effect at specific time point within mixed effects model
Hi Spencer, Thanks for your reply. I don't think this answers my question. If I understand correctly, your model simply removes the intercept and thus the intercept in fm1 is the same as the first time factor in fm1a ... but am I confused as to why the other coefficient estimates are now different for the time factor if this is just a re-naming. The coefficient estimates for the interactions are the same for fm1 and fm1a, as expected. But my question relates to the signifcance of drug at a specific time point, e.g., time = 3. The coeffecieint for say factor(time)3:drugP measures the interaction of the effect of drug=P and time=3, which is not testing what I want to test. Based on the info below, I want to compare 3) versus 4). 1) time=1, Drug=I : Intercept 2) time=1, Drug=P : Intercept + DrugP 3) time=3, Drug=I : Intercept + factor(time)3 4) time=3, Drug=P : Intercept + factor(time)3 + DrugP + factor(time)3:drugP I'm surprised this isn't simple or maybe I'm missing something competely. thanks dave -Original Message- From: Spencer Graves [mailto:[EMAIL PROTECTED] Sent: Wednesday, October 04, 2006 7:11 PM To: Afshartous, David Cc: r-help@stat.math.ethz.ch Subject: Re: [R] treatment effect at specific time point within mixed effects model Consider the following modification of your example: fm1a = lme(z ~ (factor(Time)-1)*drug, data = data.grp, random = list(Patient = ~ 1) ) summary(fm1a) snip Value Std.Error DFt-value p-value factor(Time)1 -0.6238472 0.7170161 10 -0.8700602 0.4047 factor(Time)2 -1.0155283 0.7170161 10 -1.4163256 0.1871 factor(Time)30.1446512 0.7170161 10 0.2017405 0.8442 factor(Time)40.7751736 0.7170161 10 1.0811105 0.3050 factor(Time)50.1566588 0.7170161 10 0.2184871 0.8314 factor(Time)60.0616839 0.7170161 10 0.0860286 0.9331 drugP1.2781723 1.0140139 3 1.2605077 0.2966 factor(Time)2:drugP 0.4034690 1.4340322 10 0.2813528 0.7842 factor(Time)3:drugP -0.6754441 1.4340322 10 -0.4710104 0.6477 factor(Time)4:drugP -1.8149720 1.4340322 10 -1.2656424 0.2343 factor(Time)5:drugP -0.6416580 1.4340322 10 -0.4474502 0.6641 factor(Time)6:drugP -2.1396105 1.4340322 10 -1.4920240 0.1666 Does this answer your question? Hope this helps. Spencer Graves Afshartous, David wrote: All, The code below is for a pseudo dataset of repeated measures on patients where there is also a treatment factor called drug. Time is treated as categorical. What code is necessary to test for a treatment effect at a single time point, e.g., time = 3? Does the answer matter if the design is a crossover design, i.e, each patient received drug and placebo? Finally, what would be a good response to someone that suggests to do a simple t-test (paired in crossover case) instead of the test above within a mixed model? thanks! dave z = rnorm(24, mean=0, sd=1) time = rep(1:6, 4) Patient = rep(1:4, each = 6) drug = factor(rep(c(I, P), each = 6, times = 2)) ## P = placebo, I = Ibuprofen dat.new = data.frame(time, drug, z, Patient) data.grp = groupedData(z ~ time | Patient, data = dat.new) fm1 = lme(z ~ factor(time) + drug + factor(time):drug, data = data.grp, random = list(Patient = ~ 1) ) __ 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. __ 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] Block comments in R?
On 10/5/2006 10:57 AM, Martin Maechler wrote: ## ... ## but I'm not yet fond of the general idea. Martin Is that just because you don't need it, or that you see something objectionable in it? Duncan Murdoch __ 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] Plotting text with lattice
This requires R 2.4.0. Its the same as my earlier example except the labels= arg has been changed to labels=which.packet(). xyplot(x ~ x | g, data = data.frame(x = 1:12), panel = function(...) { panel.xyplot(...) panel.text(x=2, y=4, labels=which.packet()) }) On 10/5/06, Ritwik Sinha [EMAIL PROTECTED] wrote: Hi, On a related note, if I wanted to add different texts to different panels, should I stick to using trellis.focus() for each text in each panel? I cannot figure out a way to do it using a panel function. Ritwik. On 9/29/06, Deepayan Sarkar [EMAIL PROTECTED] wrote: On 9/29/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: Here are two possibilities. The first use trellis.focus/trellis/unfocus to add text subsequent to drawing the xyplot and the second uses a custom panel: xyplot(x ~ x, data = data.frame(x = 1:10)) trellis.focus(panel, 1, 1) panel.text(x=2, y=4, labels=Text) trellis.unfocus() xyplot(x ~ x, data = data.frame(x = 1:10), panel = function(...) { panel.xyplot(...) panel.text(x=2, y=4, labels=Text) }) Right, on a more general note, this is necessary because it is not clear what ltext(x=2, y=4, labels=Text) should do for a multi-panel plot. On an even more general note, you will keep getting in trouble when using lattice if you (1) try to follow the incremental addition approach of standard graphics or (2) use the par() system in any way. Deepayan On 9/29/06, McGehee, Robert [EMAIL PROTECTED] wrote: Hello, I've decided to take the leap and try my hand at the lattice package, though I am getting stuck at what one might consider a trivial problem, plotting text at a point in a graph. Apologies in advance if (that) I'm missing something extremely basic. Consider in base graphics: plot(1:10) text(2, 4, Text) In the above you will see text centered at the point (2, 4) on the graph. Now I would like to try to do the same thing using the lattice package: xyplot(x ~ x, data = data.frame(x = 1:10)) ltext(x=2, y=4, labels=Text) panel.text(x=2, y=4, labels=Text) grid.text(label=Text, x=2, y=4) grid.text(label=Text, x=unit(2, native), y=unit(4, native)) None of the above four commands puts the Text at the (2, 4) point on the graph. Any help with this would be appreciated! Also, if I have more than one panel and would like to place text at different points on different panels how would I do this? Also, note that I'm hoping to use text to label interesting points in a levelplot, but am using the above xyplot as an example. Thanks, Robert Robert McGehee Quantitative Analyst Geode Capital Management, LLC 53 State Street, 5th Floor | Boston, MA | 02109 Tel: 617/392-8396Fax:617/476-6389 mailto:[EMAIL PROTECTED] This e-mail, and any attachments hereto, are intended for us...{{dropped}} __ 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. __ 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. -- http://www.stat.wisc.edu/~deepayan/ __ 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. -- Ritwik Sinha Graduate Student Epidemiology and Biostatistics Case Western Reserve University [EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha __ 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. __ 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] Block comments in R?
Seth Falcon wrote: My wtf feature request is to add a multiline string delimiter ala Python like . Anything that makes R more like Python syntax gets a +1 from me. How about getting rid of curly brackets for blocks and using indenting in R? Time to attack gram.y with a sledgehammer... Barry __ 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.
[R] lmer BIC changes between output and anova
list, i am using lmer to fit multilevel models and trying to use anova to compare the models. however, whenever i run the anova, the AIC, BIC and loglik are different from the original model output- as below. can someone help me out with why this is happening? (i'm hoping the output assocaited with the anova is right!). thank you, darren unconditional-lmer(log50 ~ 1 + (1 | Stream:Site) + (1|Stream), data) summary(unconditional) Linear mixed-effects model fit by REML Formula: log50 ~ 1 + (1 | Stream:Site) + (1 | Stream) Data: data AICBIC logLik MLdeviance REMLdeviance -138.8 -132.8 72.42 -150.4 -144.8 nosection-lmer(log50 ~ 1 + meanlogATS + residuallogATS + (1 | Stream:Site) + (1|Stream), data) summary(nosection) Linear mixed-effects model fit by REML Formula: log50 ~ 1 + meanlogATS + residuallogATS + (1 | Stream:Site) + (1 | Stream) Data: data AICBIC logLik MLdeviance REMLdeviance -140.8 -130.7 75.4 -168.9 -150.8 anova(unconditional, nosection) Data: data Models: unconditional: log50 ~ 1 + (1 | Stream:Site) + (1 | Stream) nosection: log50 ~ 1 + meanlogATS + residuallogATS + (1 | Stream:Site) + unconditional: (1 | Stream) Df AIC BIC logLik Chisq Chi Df Pr(Chisq) unconditional 3 -144.370 -138.294 75.185 nosection 5 -158.861 -148.734 84.430 18.491 2 9.657e-05 *** __ 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] glm with nesting
Jeffrey: Please... May I repeat what Peter Dalgaard already said: consult a local statistician. The structure of your study is sufficiently complicated that your stat 101 training is inadequate. Get professional help, which this list is not set up to provide (though it often does, through the good offices and patience of many wise contributors). Bert Gunter Genentech Nonclinical Statistics South San Francisco, CA 94404 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Jeffrey Stratford Sent: Thursday, October 05, 2006 7:46 AM To: [EMAIL PROTECTED]; r-help@stat.math.ethz.ch Subject: Re: [R] glm with nesting Harold and list, I've changed a few things since the last time so I'm really starting from scratch. I start with bbmale - read.csv(c:\\eabl\\2004\\feathers\\male_feathers2.csv, header=TRUE) box -factor(box) chick - factor(chick) Here's a sample of the data box,chick,julian,cltchsz,mrtot,cuv,cblue,purbank,purban2,purba n1,pgrassk,pgrass2,pgrass1,grassdist,grasspatchk 1,2,141,2,21.72290152,0.305723811,0.327178868,0.003813435,0.02 684564,0.06896552,0.3282487,0.6845638,0.7586207,0,3.73 4,1,164,4,18.87699007,0.281863299,0.310935559,0.06072162,0.208 0537,0.06896552,0.01936052,0,0,323.1099,0.2284615 4,2,164,4,19.64359348,0.294117388,0.316049817,0.06072162,0.208 0537,0.06896552,0.01936052,0,0,323.1099,0.2284615 7,1,118,4,13.48699876,0.303649408,0.31765218,0.3807568,0.43624 16,0.6896552,0.06864183,0.03355705,0,94.86833,0.468 12,1,180,4,21.42196378,0.289731361,0.317562294,0.09238011,0.13 42282,0,0.2430127,0.8322148,1,0,1.199032 12,2,180,4,18.79487905,0.286052077,0.316367349,0.09238011,0.13 42282,0,0.2430127,0.8322148,1,0,1.199032 12,3,180,4,12.83127682,0.260197475,0.292636914,0.09238011,0.13 42282,0,0.2430127,0.8322148,1,0,1.199032 15,1,138,4,20.07161467,0.287632782,0.318671887,0.07046477,0.03 355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818 15,2,138,4,17.61146256,0.305581768,0.315848051,0.07046477,0.03 355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818 15,3,138,4,20.36397134,0.271795667,0.30539683,0.07046477,0.033 55705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818 15,4,138,4,20.81940158,0.269468041,0.304160648,0.07046477,0.03 355705,0.03448276,0.2755622,0.6577181,0.8275862,0,1.503818 As you can see I have multiple boxes ( 70). Sometimes I have multiple chicks per box each having their own response to mrtot, cuv, and cblue but the same landscape variables for that box. Chick number is randomly assigned and is not an effect I'm interested in. I'm really not interested in the box effect either. I would like to know if landscape affects the color of chicks (which may be tied into chick health/physiology). We also know that chicks get bluer as the season progresses and that clutch size (cltchsz) has an effect so I'm including that as covariates. Hopefully, this clears things up a bit. I do have the MASS and MEMS (Pineiro's) texts in hand. Many thanks, Jeff __ 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. __ 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] treatment effect at specific time point within mixedeffects model
Hi David: In looking at your original post it is a bit difficult to ascertain exactly what your null hypothesis was. That is, you want to assess whether there is a treatment effect at time 3, but compared to what. I think your second post clears this up. You should refer to pages 224- 225 of Pinhiero and Bates for your answer. This shows how to specify contrasts. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Afshartous, David Sent: Thursday, October 05, 2006 11:08 AM To: Spencer Graves Cc: r-help@stat.math.ethz.ch Subject: Re: [R] treatment effect at specific time point within mixedeffects model Hi Spencer, Thanks for your reply. I don't think this answers my question. If I understand correctly, your model simply removes the intercept and thus the intercept in fm1 is the same as the first time factor in fm1a ... but am I confused as to why the other coefficient estimates are now different for the time factor if this is just a re-naming. The coefficient estimates for the interactions are the same for fm1 and fm1a, as expected. But my question relates to the signifcance of drug at a specific time point, e.g., time = 3. The coeffecieint for say factor(time)3:drugP measures the interaction of the effect of drug=P and time=3, which is not testing what I want to test. Based on the info below, I want to compare 3) versus 4). 1) time=1, Drug=I : Intercept 2) time=1, Drug=P : Intercept + DrugP 3) time=3, Drug=I : Intercept + factor(time)3 4) time=3, Drug=P : Intercept + factor(time)3 + DrugP + factor(time)3:drugP I'm surprised this isn't simple or maybe I'm missing something competely. thanks dave -Original Message- From: Spencer Graves [mailto:[EMAIL PROTECTED] Sent: Wednesday, October 04, 2006 7:11 PM To: Afshartous, David Cc: r-help@stat.math.ethz.ch Subject: Re: [R] treatment effect at specific time point within mixed effects model Consider the following modification of your example: fm1a = lme(z ~ (factor(Time)-1)*drug, data = data.grp, random = list(Patient = ~ 1) ) summary(fm1a) snip Value Std.Error DFt-value p-value factor(Time)1 -0.6238472 0.7170161 10 -0.8700602 0.4047 factor(Time)2 -1.0155283 0.7170161 10 -1.4163256 0.1871 factor(Time)30.1446512 0.7170161 10 0.2017405 0.8442 factor(Time)40.7751736 0.7170161 10 1.0811105 0.3050 factor(Time)50.1566588 0.7170161 10 0.2184871 0.8314 factor(Time)60.0616839 0.7170161 10 0.0860286 0.9331 drugP1.2781723 1.0140139 3 1.2605077 0.2966 factor(Time)2:drugP 0.4034690 1.4340322 10 0.2813528 0.7842 factor(Time)3:drugP -0.6754441 1.4340322 10 -0.4710104 0.6477 factor(Time)4:drugP -1.8149720 1.4340322 10 -1.2656424 0.2343 factor(Time)5:drugP -0.6416580 1.4340322 10 -0.4474502 0.6641 factor(Time)6:drugP -2.1396105 1.4340322 10 -1.4920240 0.1666 Does this answer your question? Hope this helps. Spencer Graves Afshartous, David wrote: All, The code below is for a pseudo dataset of repeated measures on patients where there is also a treatment factor called drug. Time is treated as categorical. What code is necessary to test for a treatment effect at a single time point, e.g., time = 3? Does the answer matter if the design is a crossover design, i.e, each patient received drug and placebo? Finally, what would be a good response to someone that suggests to do a simple t-test (paired in crossover case) instead of the test above within a mixed model? thanks! dave z = rnorm(24, mean=0, sd=1) time = rep(1:6, 4) Patient = rep(1:4, each = 6) drug = factor(rep(c(I, P), each = 6, times = 2)) ## P = placebo, I = Ibuprofen dat.new = data.frame(time, drug, z, Patient) data.grp = groupedData(z ~ time | Patient, data = dat.new) fm1 = lme(z ~ factor(time) + drug + factor(time):drug, data = data.grp, random = list(Patient = ~ 1) ) __ 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. __ 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. __ 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.
[R] unexpected behavior of boxplot(x, notch=TRUE, log=y)
A function I've been using for a while returned a surprising [to me, given the data] error recently: Error in plot.window(xlim, ylim, log, asp, ...) : Logarithmic axis must have positive limits After some digging I realized what was going on: x - c(10460.97, 10808.67, 29499.98, 1, 35818.62, 48535.59, 1, 1, 42512.1, 1627.39, 1, 7571.06, 21479.69, 25, 1, 16143.85, 12736.96, 1, 7603.63, 1, 33155.24, 1, 1, 50, 3361.78, 1, 37781.84, 1, 1, 1, 46492.05, 22334.88, 1, 1) summary(x) boxplot(x,notch=TRUE,log=y) #unexpected boxplot(x) #ok boxplot(x,log=y) #ok boxplot(x,notch=TRUE) #aha I can get around this, but thought that maybe boxplot() should be adjusted to deal with something like this on its own. Thank you, b. platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major 2 minor 4.0 year 2006 month 10 day03 svn rev39566 language R version.string R version 2.4.0 (2006-10-03) __ 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.
[R] Fw: Bivariate Weibull distribution -- Copula
I repost this in order for more responses. Thanks! From: Jenny Stadt Sent: 2006-10-04 16:49:33 To: r-help@stat.math.ethz.ch CC: Subject: [R] Bivariate Weibull distribution -- Copula Hi All, I am struggling in a bivariate Weibull distribution although I searched R-Site-Help and found suggestion with Copula. Seems the maximum likelihood estimate is beyond what I can understand. My case is: given two known marginal distribution (both are Weibull), and the correlation between them. How can I get the bivariate distribution based on these two? I wish if there is anybody has the experience in bivariate distribution could give me some hints. Thanks! Jen. [[alternative HTML version deleted]] __ 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. [[alternative HTML version deleted]] __ 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] Block comments in R?
Duncan == Duncan Murdoch [EMAIL PROTECTED] on Thu, 05 Oct 2006 11:13:03 -0400 writes: Duncan On 10/5/2006 10:57 AM, Martin Maechler wrote: ## ... ## but I'm not yet fond of the general idea. Martin Duncan Is that just because you don't need it, or that you see something Duncan objectionable in it? The combination: It complicates the S language (slightly) which I think is too expensive given that I don't see the need for it (smart editors ...; if(FALSE) { .. }). {But my feelings are not very strong here probably mostly because I haven't thought very long about it.} Martin __ 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.
[R] mixed models: correlation between fixed and random effects??
Hello, I built 4 mixed models using different data sets and standardized variables as predictors. In all the models each of the fixed effects has an associated random effect (same predictor). What I find is that fixed effects with larger (absolute) standardized parameter estimates have also a higher estimate of the related random effect. In other words, the higher the average of the absolute value of the BLUPs for a given standardized parameter, the higher its variance. Is this a common situation or I am missing some additional normalization factor necessary to compare the different random effects? Thanks a lot, Bruno ~~ Bruno L. Giordano, Ph.D. CIRMMT Schulich School of Music, McGill University 555 Sherbrooke Street West Montréal, QC H3A 1E3 Canada http://www.music.mcgill.ca/~bruno/ __ 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] Block comments in R?
There are two places that I find the current way it works to be less than ideal although its not that bad either: 1 .when one wants to define strings that have quotes then one must be careful to use double quoted strings to contain single quotes and single quoted strings to contain double quotes or if one needs both then one is into the complexities of multiple backslashes. Enhancing sQuote and dQuote might help somewhat in this regards to allow them to be used in more situations but having an extended syntax might be the best. 2. Sometime I like to define a file inline so that it is defined like this: Lines - 1 2 3 4 read.table(textConnection(Lines)) I haven't hit the size limit on the character string yet (I usually do this with small config style files) but suspect that I may if I start using this for data. What is the limit anyways? On 10/5/06, Martin Maechler [EMAIL PROTECTED] wrote: Duncan == Duncan Murdoch [EMAIL PROTECTED] on Thu, 05 Oct 2006 11:13:03 -0400 writes: Duncan On 10/5/2006 10:57 AM, Martin Maechler wrote: ## ... ## but I'm not yet fond of the general idea. Martin Duncan Is that just because you don't need it, or that you see something Duncan objectionable in it? The combination: It complicates the S language (slightly) which I think is too expensive given that I don't see the need for it (smart editors ...; if(FALSE) { .. }). {But my feelings are not very strong here probably mostly because I haven't thought very long about it.} Martin __ 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. __ 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] treatment effect at specific time point within mixedeffects model
Hi Harold, Thanks for your response. I'll check out p.224 in PB, thanks. The null hypothesis is that there is no difference between say A=[time=3, drug=I] and B=[time=3, drug=P], or mu_A = mu_B. If the study is a crossover design, i.e., each patient receives drug=I and drug=P, I assume that a simple paried t-test could also be employed at time=3. However, I'd like to test this within a mixed effects model; With respect to 3) and 4) below, it seems somewhat difficult to express this specific hypothesis in terms of the model paramaters. Ways in which this null are violated under the mixed effects models could be: 1) there is no interaction between time and Drug, i.e., there is a drug effect but it is the same at all time points. (the specific interaction in 3) below represents the shift of the effect of drug=P from time=1 to time=3 ... so the lack of significance of the paramater factor(time)3:drugP doesn't capture what I want) 2) there is neither interaction nor drug effect (variable Drug not significant). But both these violations are more general than my null; I think testing fixed effects 3) versus 4) below is what I want, but this also seems strange since possibly the drug effect and drug:time interaction as defined in the model are signicant (with time=1 as the reference baseline). Regardless, I assume I would need to employ coef() and vcov() to obtain the needed info ... but I notice that coef() produces 4 values for the intercept of fm1 below, does anyone know why this occurs? I apologize if my explanation above is confusing, I've tried to make it as clear as possible. thanks again, dave -Original Message- From: Doran, Harold [mailto:[EMAIL PROTECTED] Sent: Thursday, October 05, 2006 11:40 AM To: Afshartous, David; Spencer Graves Cc: r-help@stat.math.ethz.ch Subject: RE: [R] treatment effect at specific time point within mixedeffects model Hi David: In looking at your original post it is a bit difficult to ascertain exactly what your null hypothesis was. That is, you want to assess whether there is a treatment effect at time 3, but compared to what. I think your second post clears this up. You should refer to pages 224- 225 of Pinhiero and Bates for your answer. This shows how to specify contrasts. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Afshartous, David Sent: Thursday, October 05, 2006 11:08 AM To: Spencer Graves Cc: r-help@stat.math.ethz.ch Subject: Re: [R] treatment effect at specific time point within mixedeffects model Hi Spencer, Thanks for your reply. I don't think this answers my question. If I understand correctly, your model simply removes the intercept and thus the intercept in fm1 is the same as the first time factor in fm1a ... but am I confused as to why the other coefficient estimates are now different for the time factor if this is just a re-naming. The coefficient estimates for the interactions are the same for fm1 and fm1a, as expected. But my question relates to the signifcance of drug at a specific time point, e.g., time = 3. The coeffecieint for say factor(time)3:drugP measures the interaction of the effect of drug=P and time=3, which is not testing what I want to test. Based on the info below, I want to compare 3) versus 4). 1) time=1, Drug=I : Intercept 2) time=1, Drug=P : Intercept + DrugP 3) time=3, Drug=I : Intercept + factor(time)3 4) time=3, Drug=P : Intercept + factor(time)3 + DrugP + factor(time)3:drugP I'm surprised this isn't simple or maybe I'm missing something competely. thanks dave -Original Message- From: Spencer Graves [mailto:[EMAIL PROTECTED] Sent: Wednesday, October 04, 2006 7:11 PM To: Afshartous, David Cc: r-help@stat.math.ethz.ch Subject: Re: [R] treatment effect at specific time point within mixed effects model Consider the following modification of your example: fm1a = lme(z ~ (factor(Time)-1)*drug, data = data.grp, random = list(Patient = ~ 1) ) summary(fm1a) snip Value Std.Error DFt-value p-value factor(Time)1 -0.6238472 0.7170161 10 -0.8700602 0.4047 factor(Time)2 -1.0155283 0.7170161 10 -1.4163256 0.1871 factor(Time)30.1446512 0.7170161 10 0.2017405 0.8442 factor(Time)40.7751736 0.7170161 10 1.0811105 0.3050 factor(Time)50.1566588 0.7170161 10 0.2184871 0.8314 factor(Time)60.0616839 0.7170161 10 0.0860286 0.9331 drugP1.2781723 1.0140139 3 1.2605077 0.2966 factor(Time)2:drugP 0.4034690 1.4340322 10 0.2813528 0.7842 factor(Time)3:drugP -0.6754441 1.4340322 10 -0.4710104 0.6477 factor(Time)4:drugP -1.8149720 1.4340322 10 -1.2656424 0.2343 factor(Time)5:drugP -0.6416580 1.4340322 10 -0.4474502 0.6641 factor(Time)6:drugP -2.1396105 1.4340322 10 -1.4920240 0.1666 Does this answer your question? Hope this helps.
Re: [R] (not so real) function (more like #include)
I thought of a few more possibilities. If options 1, 2 and 3 correspond to longfun, longfun2 and longfun3 then: 4. We could arrange it so that the variables are stored outside of longfun and then longfun, f1 and f2 reference them symmetrically. a - b - d - 1 longfun4 - function() { out - f1() # since a is not in longfun4 use - to set it a - a + 1 f2() + out } f1 - f2 - function() a + b + d longfun4() Given the problem as stated I don't think this approach is really the best but if the problem gets extended to one where there is not just one longfun but several longfuns then it will not be feasible to store all the variables in longfun since the other longfuns won't be able to access them and this approach becomes desirable. 5. We could also use the proto package for this purpose placing a, b, d, longfun, f1 and f2 into an object p. If we wanted to be able to override f1, f2, a, b and d in delegated subobjects and be able to refer to f1 and f2 outside of longfun, i.e. outside of the context of object p, then we would have to refer to them as .$f1, .$f2, .$a, .$b and .$d and also pass the target object as arg 1 of f1 and f2 but we have not bothered with that here since it seems not to be a requirement: library(proto) p - proto(a = 1, b = 1, d = 1) p$longfun - function(.) { out - f1() # since a is not in longfun4 use - to set it a - a + 1 f2() + out } p$f1 - p$f2 - function() a + b + d p$longfun() # 7 This is probably overkill here but it might be that if one thinks about the design some more that there are some opportunities to define subobjects of p in which case the ability to do inheritance would come into play. 6. I think the R.oo package would work too although I am not sure its as suitable here since it would require the definition of classes which you don't really need. proto lets you define objects directly even without classes. On 10/5/06, Meinhard Ploner [EMAIL PROTECTED] wrote: Thanks a lot! longfun2 longfun3 work perfect for my very untypical problem. As I have many local variables, usual functions with parameters are very uncomfortable, but the code you gave me is great! Meinhard On Oct 4, 2006, at 5:25 PM, Gabor Grothendieck wrote: longfun could just pass a, b and d to each of the individual functions and each of the individual functions could pass out back as a return value. f1 - f2 - function(a, b, d) a+b+d longfun1 - function() { a - b - d - 1 out - f1(a, b, d) out - f2(a, b, d) + out out } longfun1() # 6 If the problem is that a, b and d actually represent 100 variables, say, and its a nuisance to pass them all explicitly you could do this: f1 - f2 - function() with(parent.frame(), a + b + d) longfun2 - function() { a - b - d - 1 out - f1() out - f2() + out out } longfun2() # 6 or you could do it this way if you would prefer to have longfun control the scoping rather than having it done within each individual function: f1 - f2 - function() a+b+d longfun3 - function() { a - b - d - 1 environment(f1) - environment(f2) - environment() out - f1() f2() + out } longfun3() # 6 On 10/4/06, Meinhard Ploner [EMAIL PROTECTED] wrote: Hello R users, I have a very long function with parts of that looking like a list of jobs. The first part of the function defines a lot of local variables, which are used by the jobs. Now I extracted the job part und putted them into an external file, used by source(, local=T), which is more comfortable to me. Apart from that it gives me the opportunity that more users can work on the project. Is it possible to define a function for that code and passing to it the environment of f() without save(list=ls()) and load() ? Thanks in advance Meinhard Ploner longfun - f() { ## lot of local variables: a - ... b - ... d - ... ... out - ... source(job1.txt, local=TRUE) #changes out, uses a, b, d, ... source(job2.txt, local=TRUE) # changes out, uses a, b, d, ... ... } version _ platform i386-apple-darwin8.6.1 arch i386 os darwin8.6.1 system i386, darwin8.6.1 status major 2 minor 3.1 year 2006 month 06 day01 svn rev38247 language R version.string Version 2.3.1 (2006-06-01) __ 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. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide
[R] any package can visualize original affymetrix data?
Dear All, I'm thinking of visualizing the original binary data from affymetrix. Would you recommend any package? Not necessarily limited to R packages. Thanks! Best, -Cao __ 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.
[R] convert day of week from number to character and include in lm
All, I am trying to include a day of week variable (1-7) in in a regression model. I would like to have the day of week treated as a categorical variable rather than a number the code looks like lm( dep ~ WKDY) I know this is a basic question, but help would be appreciated thanks spencer [[alternative HTML version deleted]] __ 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] treatment effect at specific time point within mixedeffects model
Afshartous, David wrote: Hi Harold, Thanks for your response. I'll check out p.224 in PB, thanks. The null hypothesis is that there is no difference between say A=[time=3, drug=I] and B=[time=3, drug=P], or mu_A = mu_B. If the study is a crossover design, i.e., each patient receives drug=I and drug=P, I assume that a simple paried t-test could also be employed at time=3. However, I'd like to test this within a mixed effects model; With respect to 3) and 4) below, it seems somewhat difficult to express this specific hypothesis in terms of the model paramaters. Ways in which this null are violated under the mixed effects models could be: 1) there is no interaction between time and Drug, i.e., there is a drug effect but it is the same at all time points. (the specific interaction in 3) below represents the shift of the effect of drug=P from time=1 to time=3 ... so the lack of significance of the paramater factor(time)3:drugP doesn't capture what I want) If there is no evidence for an interaction between drug and time, do you still want to ask about the drug effect at time=3, or would you then want to ask about the time-averaged drug effect? 2) there is neither interaction nor drug effect (variable Drug not significant). But both these violations are more general than my null; I think testing fixed effects 3) versus 4) below is what I want, but this also seems strange since possibly the drug effect and drug:time interaction as defined in the model are signicant (with time=1 as the reference baseline). If you fit a model with an intercept, main effects for drug and time, and an interaction, would'nt the coefficient for the drug main effect test the drug effect at a particular time? Perhaps you only need to change the contrasts for time so that time=3 is the reference category? Regardless, I assume I would need to employ coef() and vcov() to obtain the needed info ... but I notice that coef() produces 4 values for the intercept of fm1 below, does anyone know why this occurs? I think Harold was getting at the fact that you could get an estimate of the drug effect at time=3 simply by setting the contrasts for time in the right way. I apologize if my explanation above is confusing, I've tried to make it as clear as possible. thanks again, dave -Original Message- From: Doran, Harold [mailto:[EMAIL PROTECTED] Sent: Thursday, October 05, 2006 11:40 AM To: Afshartous, David; Spencer Graves Cc: r-help@stat.math.ethz.ch Subject: RE: [R] treatment effect at specific time point within mixedeffects model Hi David: In looking at your original post it is a bit difficult to ascertain exactly what your null hypothesis was. That is, you want to assess whether there is a treatment effect at time 3, but compared to what. I think your second post clears this up. You should refer to pages 224- 225 of Pinhiero and Bates for your answer. This shows how to specify contrasts. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Afshartous, David Sent: Thursday, October 05, 2006 11:08 AM To: Spencer Graves Cc: r-help@stat.math.ethz.ch Subject: Re: [R] treatment effect at specific time point within mixedeffects model Hi Spencer, Thanks for your reply. I don't think this answers my question. If I understand correctly, your model simply removes the intercept and thus the intercept in fm1 is the same as the first time factor in fm1a ... but am I confused as to why the other coefficient estimates are now different for the time factor if this is just a re-naming. The coefficient estimates for the interactions are the same for fm1 and fm1a, as expected. But my question relates to the signifcance of drug at a specific time point, e.g., time = 3. The coeffecieint for say factor(time)3:drugP measures the interaction of the effect of drug=P and time=3, which is not testing what I want to test. Based on the info below, I want to compare 3) versus 4). 1) time=1, Drug=I : Intercept 2) time=1, Drug=P : Intercept + DrugP 3) time=3, Drug=I : Intercept + factor(time)3 4) time=3, Drug=P : Intercept + factor(time)3 + DrugP + factor(time)3:drugP I'm surprised this isn't simple or maybe I'm missing something competely. thanks dave -Original Message- From: Spencer Graves [mailto:[EMAIL PROTECTED] Sent: Wednesday, October 04, 2006 7:11 PM To: Afshartous, David Cc: r-help@stat.math.ethz.ch Subject: Re: [R] treatment effect at specific time point within mixed effects model Consider the following modification of your example: fm1a = lme(z ~ (factor(Time)-1)*drug, data = data.grp, random = list(Patient = ~ 1) ) summary(fm1a) snip Value Std.Error DFt-value p-value factor(Time)1 -0.6238472 0.7170161 10 -0.8700602 0.4047 factor(Time)2 -1.0155283 0.7170161 10 -1.4163256 0.1871
Re: [R] Block comments in R?
+1 for Python style comments. your comments here and some more ... and some related info jab -- John Bollinger, CFA, CMT www.BollingerBands.com If you advance far enough, you arrive at the beginning. __ 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.
[R] The W statistic in wilcox.exact
Does anyone know why wilcox.exact gives W-statistic 6 instead of 12 as indicated below. 12 is the rank sum of group 0 of x, which is the linear statistic computed by wilcox_test. y-c(1,2,3,4,5) x-c(1,1,0,0,0) (a) wilcox.exact wilcox.exact(y~x) Exact Wilcoxon rank sum test data: y by x W = 6, p-value = 0.2 alternative hypothesis: true mu is not equal to 0 (b) wilcox_test tt-wilcox_test(y~factor(x),distribution=exact) statistic(tt,linear) 0 12 Jue Wang, Biostatistician Contracted Position for Preclinical Research Biostatistics PrO Unlimited (908) 231-3022 __ 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.
[R] Multivariate AR - prediction
Hi, does anybody know how to predict a multivariate AR within R? If I just estimate a multi AR-object and plug it into predict I get an error from the aperm - just works for univariates. thx alex __ 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] Plotting text with lattice
I seem to have omitted g library(lattice) xyplot(x ~ x | g, data = data.frame(x = 1:12, g = gl(3,4)), panel = function(...) { panel.xyplot(...) panel.text(x=2, y=4, labels=which.packet()) }) On 10/5/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: This requires R 2.4.0. Its the same as my earlier example except the labels= arg has been changed to labels=which.packet(). xyplot(x ~ x | g, data = data.frame(x = 1:12), panel = function(...) { panel.xyplot(...) panel.text(x=2, y=4, labels=which.packet()) }) On 10/5/06, Ritwik Sinha [EMAIL PROTECTED] wrote: Hi, On a related note, if I wanted to add different texts to different panels, should I stick to using trellis.focus() for each text in each panel? I cannot figure out a way to do it using a panel function. Ritwik. On 9/29/06, Deepayan Sarkar [EMAIL PROTECTED] wrote: On 9/29/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: Here are two possibilities. The first use trellis.focus/trellis/unfocus to add text subsequent to drawing the xyplot and the second uses a custom panel: xyplot(x ~ x, data = data.frame(x = 1:10)) trellis.focus(panel, 1, 1) panel.text(x=2, y=4, labels=Text) trellis.unfocus() xyplot(x ~ x, data = data.frame(x = 1:10), panel = function(...) { panel.xyplot(...) panel.text(x=2, y=4, labels=Text) }) Right, on a more general note, this is necessary because it is not clear what ltext(x=2, y=4, labels=Text) should do for a multi-panel plot. On an even more general note, you will keep getting in trouble when using lattice if you (1) try to follow the incremental addition approach of standard graphics or (2) use the par() system in any way. Deepayan On 9/29/06, McGehee, Robert [EMAIL PROTECTED] wrote: Hello, I've decided to take the leap and try my hand at the lattice package, though I am getting stuck at what one might consider a trivial problem, plotting text at a point in a graph. Apologies in advance if (that) I'm missing something extremely basic. Consider in base graphics: plot(1:10) text(2, 4, Text) In the above you will see text centered at the point (2, 4) on the graph. Now I would like to try to do the same thing using the lattice package: xyplot(x ~ x, data = data.frame(x = 1:10)) ltext(x=2, y=4, labels=Text) panel.text(x=2, y=4, labels=Text) grid.text(label=Text, x=2, y=4) grid.text(label=Text, x=unit(2, native), y=unit(4, native)) None of the above four commands puts the Text at the (2, 4) point on the graph. Any help with this would be appreciated! Also, if I have more than one panel and would like to place text at different points on different panels how would I do this? Also, note that I'm hoping to use text to label interesting points in a levelplot, but am using the above xyplot as an example. Thanks, Robert Robert McGehee Quantitative Analyst Geode Capital Management, LLC 53 State Street, 5th Floor | Boston, MA | 02109 Tel: 617/392-8396Fax:617/476-6389 mailto:[EMAIL PROTECTED] This e-mail, and any attachments hereto, are intended for us...{{dropped}} __ 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. __ 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. -- http://www.stat.wisc.edu/~deepayan/ __ 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. -- Ritwik Sinha Graduate Student Epidemiology and Biostatistics Case Western Reserve University [EMAIL PROTECTED] | +12163682366 | http://darwin.cwru.edu/~rsinha __ 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. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read
[R] [R-pkgs] new package proptest
Dear R users, The package proptest (http://www.davidkraus.net/proptest/) has been recently released to CRAN. The package provides functions for testing the proportional hazards assumption in the Cox model for right censored survival data. Two types of tests for identifying nonproportional covariates are implemented: * data-driven Neyman type smooth tests (using both nested and all-subsets variant of the BIC), * tests based on the score process (KS, CM, and AD type, with the LWY simulation). Bug reports, comments or suggestions are welcome. Regards, David. -- David Kraus Institute of Information Theory and Automation Pod Vodarenskou vezi 4 18208 Praha 8 Czechia Charles University in Prague Department of Statistics Sokolovska 83 18675 Praha 8 Czechia kraus #at# karlin.mff.cuni.cz http://www.davidkraus.net/ ___ R-packages mailing list R-packages@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-packages __ 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] R Graphics: Saving PDF and other formats from Windows Graphic Device for LaTeX
Anupam Tyagi [EMAIL PROTECTED] writes: [...] What is the best format to save R graphics for inclusion into a LaTeX documents? When using pdflatex use pdf for graphics as reference format. Using ps2pdf or some such may have some problems when it comes to alpha channels and transparency. [...] I will be grateful if someone can share their experience. Indispensable is having a good editing cycle. I.e. compile the latex, jump to the position where you are editing the file in the preview, double click somewhere in the preview, move to that position in your editor (or at least near that position). As I am not aware of any tools doing this with pdf, you may want to render dvi instead. When rendering DVI, you probably want to generate ps. My 2c, Jens __ 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] littler release 0.0.6
For those of you who are littler fans, you'll be pleased to know that we've already decommissioned version 0.0.6 in light of the new and improved 0.0.7 (actually 0.0.6 was broken because of a missing file). You can get it from the links below. Jeff -- http://biostat.mc.vanderbilt.edu/JeffreyHorner Dirk Eddelbuettel wrote: What ? -- We are pleased to announce version 0.0.6 of littler What's new ? This version includes a bug fix or two as well as a number of small enhancements to the documentation. For OS X and the r/R confusion, our recommended suggestion is to call configure using either the --program-suffix=X or --program-prefix=Y option to have the binary and manual page renamed when 'make install' runs. In this example, r would be renamed to rX and Yr, resptively. Suitable choices of X and Y are left to the user. Where ? --- You can find the source tarball at both http://dirk.eddelbuettel.com/code/littler/ and http://biostat.mc.vanderbilt.edu/LittleR Comments are welcome, as are are suggestions, bug fixes, or patches. - Jeffrey Horner [EMAIL PROTECTED] - Dirk Eddelbuettel [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 and provide commented, minimal, self-contained, reproducible code.
[R] Help on plot multiple plot(hexbin) in the same page
Dear R-users, I try to plot multiple plots of hexbin in the same page. However, par does not work when hexbin is used. Neither do xlim and ylim. Any ideas? Thanks, Sue [[alternative HTML version deleted]] __ 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.
[R] Writing Text into Plots
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 Hello, I have a simple question on the text() method in plots, e.g.: text(3,4,adj=1,cex=1.0, expression(alpha == beta)) I know there exists a lot more like frac(), etc which could be used for expression. But a help(frac) doesn't return any results - where do I have to look for a documentation of possible text commands? More in detail I'm searching for how to set an index like mu_{1} and mu_{2} ? And how to get something like a bar over a character like \bar x ? TIA, Lothar -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.5 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org iD8DBQFFJWX0HRf7N9c+X7sRAtyiAJ4yssSZtt/DGzWVTfVI2qvzoO9FigCfe7Fm NTOUZN7jL/KdthkNEzM8S0M= =tBvl -END PGP SIGNATURE- __ 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] matrix multiplication
Is this what you want ? A1-matrix(c(1,0,0,0,0,0,0,0,0),3) A2-matrix(c(1,0,0,0,1,0,0,0,0),3) A3-matrix(c(1,0,0,0,1,0,0,1,0),3) X - matrix(1:24,8) XX-list() for(i in 1:3){ + XX[[i]]-X%*%A[[i]] + } XX [[1]] [,1] [,2] [,3] [1,]100 [2,]200 [3,]300 [4,]400 [5,]500 [6,]600 [7,]700 [8,]800 [[2]] [,1] [,2] [,3] [1,]190 [2,]2 100 [3,]3 110 [4,]4 120 [5,]5 130 [6,]6 140 [7,]7 150 [8,]8 160 [[3]] [,1] [,2] [,3] [1,]199 [2,]2 10 10 [3,]3 11 11 [4,]4 12 12 [5,]5 13 13 [6,]6 14 14 [7,]7 15 15 [8,]8 16 16 Ya-Hsiu Chuang wrote: Dear all, I have 2 matrices, one is a 8x3 matrix, called X; the other matrix is a 3x3 indicator matrix with the diagonal element as 0 or 1. when a variable is included in the model, the corresponding diagonal element is 1, otherwise, it is 0. Let A be a set of matrices that contain the possible indicator matrix. e.g., A= [A1, A2, A3], where A1= [1,0,0;0,0,0,0,0,0], A2 =[1,0,0;0,1,0,0,0,0], A3 =[1,0,0;0,0,0,0,1,0] In order to derive the new X matries depending on the indicator matrices, I use for loops to multiply both X and A as following p-3 for ( i in 1:p){ XX- X%*%A[i] } However, it only shows the result when i=p. How can I derive the results which include all possible value of XX[i], i=1,..,p, rather than just i=p thanks for any help Ya-Hsiu [[alternative HTML version deleted]] __ 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. -- Ferdinand Alimadhi Programmer / Analyst Harvard University The Institute for Quantitative Social Science (617) 496-0187 [EMAIL PROTECTED] www.iq.harvard.edu [[alternative HTML version deleted]] __ 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] R Graphics: Saving PDF and other formats from Windows Graphic Device for LaTeX
Jens Scheidtmann [EMAIL PROTECTED] wrote: Indispensable is having a good editing cycle. I.e. compile the latex, jump to the position where you are editing the file in the preview, double click somewhere in the preview, move to that position in your editor (or at least near that position). As I am not aware of any tools doing this with pdf, you may want to render dvi instead. Agreed. As a point of interest, the newer Windows versions of vTeX (by MicroPress) will do this (in both directions) with PDF. This is a commercial TeX distribution with value added. They have a free distribution suitable for Unix or Linux, but I doubt that it does that particular trick. -- Mike Prager, NOAA, Beaufort, NC * Opinions expressed are personal and not represented otherwise. * Any use of tradenames does not constitute a NOAA endorsement. __ 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.
[R] How to get the function names
I've defined the function getFunNames - function(FUN){ if (!is.list(FUN)) fun.names - paste(deparse(substitute(FUN)), collapse = ) else fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a))) fun.names } which gives what I want : getFunNames(mean) [1] mean getFunNames(ff) [1] ff getFunNames(c(mean,ff)) [1] mean ff If I call this within a function, things go wrong: 1] FUN foo(ff) [1] FUN foo(c(mean,ff)) Error in substitute(FUN)[-1] : object is not subsettable Obviously there are some things (quite a few things) which I have not understood. Can anyone help? Thanks Søren __ 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] The W statistic in wilcox.exact
Probably because of the offset: U = W - n*(n+1)/2 In your example, W=12 (=3+4+5) as reported by wilcox.test. The offset is 6 (=3*4/2) and therefore U=6. I am not certain as I haven't installed the exactRankTests package, but it seems that wilcox.exact reports U instead of W. -Christos Hatzis -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: Thursday, October 05, 2006 1:35 PM To: r-help@stat.math.ethz.ch Subject: [R] The W statistic in wilcox.exact Does anyone know why wilcox.exact gives W-statistic 6 instead of 12 as indicated below. 12 is the rank sum of group 0 of x, which is the linear statistic computed by wilcox_test. y-c(1,2,3,4,5) x-c(1,1,0,0,0) (a) wilcox.exact wilcox.exact(y~x) Exact Wilcoxon rank sum test data: y by x W = 6, p-value = 0.2 alternative hypothesis: true mu is not equal to 0 (b) wilcox_test tt-wilcox_test(y~factor(x),distribution=exact) statistic(tt,linear) 0 12 Jue Wang, Biostatistician Contracted Position for Preclinical Research Biostatistics PrO Unlimited (908) 231-3022 __ 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. __ 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.
[R] xyplot
Hi, for the data below: time-c(rep(1:10,5)) y-time+rnorm(50,5,2) subject-c(rep('a',10),rep('b',10),rep('c',10),rep('d',10),rep('e',10)) group-c(rep('A',30),rep('B',20)) df-data.frame(subject,group,time,y) I'd like to produce a plot with a single pannel with two loess curves one for each group. the code below does that but on two different pannels, I'd like to have both curves on the same pannel with different colors. xyplot(y~time|group, groups=subject, panel=function(...){ panel.loess(...) } ,data=df) -- Osman O. Al-Radi, MD, MSc, FRCSC Fellow, Cardiovascular Surgery The Hospital for Sick Children University of Toronto, Canada [[alternative HTML version deleted]] __ 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] Writing Text into Plots
On 10/5/2006 4:10 PM, Lothar Botelho-Machado wrote: -BEGIN PGP SIGNED MESSAGE- Hash: SHA1 Hello, I have a simple question on the text() method in plots, e.g.: text(3,4,adj=1,cex=1.0, expression(alpha == beta)) I know there exists a lot more like frac(), etc which could be used for expression. But a help(frac) doesn't return any results - where do I have to look for a documentation of possible text commands? More in detail I'm searching for how to set an index like mu_{1} and mu_{2} ? And how to get something like a bar over a character like \bar x ? See ?plotmath, and demo(plotmath). Duncan Murdoch __ 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] Writing Text into Plots
On Thu, 2006-10-05 at 22:10 +0200, Lothar Botelho-Machado wrote: Hello, I have a simple question on the text() method in plots, e.g.: text(3,4,adj=1,cex=1.0, expression(alpha == beta)) I know there exists a lot more like frac(), etc which could be used for expression. But a help(frac) doesn't return any results - where do I have to look for a documentation of possible text commands? More in detail I'm searching for how to set an index like mu_{1} and mu_{2} ? And how to get something like a bar over a character like \bar x ? See ?plotmath, which BTW, is listed in the See Also section of ?text HTH, Marc Schwartz __ 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.
[R] which command
I obtained error messages when I run these commands in UNIX, but I obtained correct result when I run these command in WINDOWS. Can somebody point out the problem and give the solution. Thanks. dt-read.table(file=Fall.dat) dim(dt) [1] 19415 table(dt$V2) 0 1 2 3 220 989 639 93 Favg-as.matrix(c(1:max(dt$V2))) for(i in 1:max(dt$V2)){Favg[i]-mean(dt[which(dt[,2] == i),5])} Favg [,1] [1,] NA [2,] NA [3,] NA Warning messages: argument is not numeric or logical: returning NA (etc) __ 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.
[R] Matrix input problem
Hi, Included is an R.script that came from a much large date set being read in to R from a .txt file. Why is it that the first line codes with an error, but the second line works fine? Is there some line length limit in R? This happens at random places all through my data. Any help would be appreciated. Bill Wyatt Associate Instructor Ergonomics Graduate Student Department of Kinesiology School of Health, Physical Education, and Recreation Indiana University O:(812)856-5924 [EMAIL PROTECTED] //orginal line b1x-cbind(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-132.5,-132.5,-132.5,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-133.6,-133.6,-133.05,-133.05,-131.35,-131.35,-129.7,-129.7,-126.9,-126.9,-124.7,-124.7,-121.35,-121.35,-118,-118,-114.7,-114.7,-110.8,-110.8,-106.9,-106.9,-102.45,-102.45,-98,-98,-92.95,-92.95,-87.95,-87.95,-82.95,-82.95,-77.95,-77.95,-72.35,-67.3 5,-67.35,-61.8,-61.8,-55.65,-55.65)) //copy of line with end )) moved one data point forward b1x-cbind(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-132.5,-132.5,-132.5,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-133.6,-133.6,-133.05,-133.05,-131.35,-131.35,-129.7,-129.7,-126.9,-126.9,-124.7,-124.7,-121.35,-121.35,-118,-118,-114.7,-114.7,-110.8,-110.8,-106.9,-106.9,-102.45,-102.45,-98,-98,-92.95,-92.95,-87.95,-87.95,-82.95,-82.95,-77.95,-77.95,-72.35,-67.3 5,-67.35,-61.8,-61.8,-55.65)) ,-55.65 __ 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] xyplot
Osman Al-Radi said the following on 10/5/2006 3:43 PM: Hi, for the data below: time-c(rep(1:10,5)) y-time+rnorm(50,5,2) subject-c(rep('a',10),rep('b',10),rep('c',10),rep('d',10),rep('e',10)) group-c(rep('A',30),rep('B',20)) df-data.frame(subject,group,time,y) I'd like to produce a plot with a single pannel with two loess curves one for each group. the code below does that but on two different pannels, I'd like to have both curves on the same pannel with different colors. xyplot(y~time|group, groups=subject, panel=function(...){ panel.loess(...) } ,data=df) Simple: xyplot(y ~ time | group, data = df, groups = subject, panel = function(...) { panel.superpose(...) panel.superpose(panel.groups = panel.loess, ...) }) HTH, --sundar __ 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.
[R] which command
I obtained error messages when I run these commands in UNIX, but I obtained correct result when I run these command in WINDOWS. Can somebody point out the problem and give the solution. Thanks. dt-read.table(file=Fall.dat) dim(dt) [1] 19415 table(dt$V2) 0 1 2 3 220 989 639 93 Favg-as.matrix(c(1:max(dt$V2))) for(i in 1:max(dt$V2)){Favg[i]-mean(dt[which(dt[,2] == i),5])} Favg [,1] [1,] NA [2,] NA [3,] NA Warning messages: argument is not numeric or logical: returning NA (etc) Cheers, @ __ 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] How to get the function names
On 10/5/2006 4:41 PM, Søren Højsgaard wrote: I've defined the function getFunNames - function(FUN){ if (!is.list(FUN)) fun.names - paste(deparse(substitute(FUN)), collapse = ) else fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a))) fun.names } which gives what I want : getFunNames(mean) [1] mean getFunNames(ff) [1] ff getFunNames(c(mean,ff)) [1] mean ff If I call this within a function, things go wrong: 1] FUN foo(ff) [1] FUN foo(c(mean,ff)) Error in substitute(FUN)[-1] : object is not subsettable Obviously there are some things (quite a few things) which I have not understood. Can anyone help? I don't think you'll be able to do what you want. The problem is that R objects don't know their own name. substitute just gives you the unevaluated expression that was passed as an arg; so getFunNames is reporting on the expression in foo() that was used when it was called. You'll need to call substitute directly on the args to get the expressions, rather than passing them to another function to do that. In C code you have a bit more flexibility for non-standard evaluation, but not in R code. Duncan Murdoch __ 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] xyplot
On 10/5/06, Osman Al-Radi [EMAIL PROTECTED] wrote: Hi, for the data below: time-c(rep(1:10,5)) y-time+rnorm(50,5,2) subject-c(rep('a',10),rep('b',10),rep('c',10),rep('d',10),rep('e',10)) group-c(rep('A',30),rep('B',20)) df-data.frame(subject,group,time,y) I'd like to produce a plot with a single pannel with two loess curves one for each group. the code below does that but on two different pannels, I'd like to have both curves on the same pannel with different colors. xyplot(y~time|group, groups=subject, panel=function(...){ panel.loess(...) } ,data=df) Well, the fix for your code is xyplot(y ~ time | group, data = df, groups = subject, panel = panel.superpose, panel.groups = function(...) panel.loess(...) ) Which, incidentally is effectively equivalent to xyplot(y ~ time | group, data = df, groups = subject, type = smooth) However, this is not what you describe as your goal. If you want a single smooth for each group combining all subjects within that group, use something like xyplot(y ~ time, data = df, groups = group. ... The panel part does not change. -Deepayan __ 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] Matrix input problem
On 10/5/2006 4:52 PM, Bill Wyatt wrote: Hi, Included is an R.script that came from a much large date set being read in to R from a .txt file. Why is it that the first line codes with an error, but the second line works fine? I get syntax errors on both, due to the missing comma at the end of the very long line. Is there some line length limit in R? Yes, there's a 1024 character length limit (which is documented in the Intro to R manual in the current release; not sure how long it's been there). Why not put in some line breaks? R source is supposed to be human readable, not just machine readable. Duncan Murdoch This happens at random places all through my data. Any help would be appreciated. Bill Wyatt Associate Instructor Ergonomics Graduate Student Department of Kinesiology School of Health, Physical Education, and Recreation Indiana University O:(812)856-5924 [EMAIL PROTECTED] //orginal line b1x-cbind(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-132.5,-132.5,-132.5,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-133.6,-133.6,-133.05,-133.05,-131.35,-131.35,-129.7,-129.7,-126.9,-126.9,-124.7,-124.7,-121.35,-121.35,-118,-118,-114.7,-114.7,-110.8,-110.8,-106.9,-106.9,-102.45,-102.45,-98,-98,-92.95,-92.95,-87.95,-87.95,-82.95,-82.95,-77.95,-77.95,-72.35,-6 7.3 5,-67.35,-61.8,-61.8,-55.65,-55.65)) //copy of line with end )) moved one data point forward b1x-cbind(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-132.5,-132.5,-132.5,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-134.15,-133.6,-133.6,-133.05,-133.05,-131.35,-131.35,-129.7,-129.7,-126.9,-126.9,-124.7,-124.7,-121.35,-121.35,-118,-118,-114.7,-114.7,-110.8,-110.8,-106.9,-106.9,-102.45,-102.45,-98,-98,-92.95,-92.95,-87.95,-87.95,-82.95,-82.95,-77.95,-77.95,-72.35,-6 7.3 5,-67.35,-61.8,-61.8,-55.65)) ,-55.65 __ 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. __ 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] How to get the function names
On Thu, 2006-10-05 at 22:41 +0200, Søren Højsgaard wrote: I've defined the function getFunNames - function(FUN){ if (!is.list(FUN)) fun.names - paste(deparse(substitute(FUN)), collapse = ) else fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a))) fun.names } Hi, Try this: getFunNames - function(x) + sapply(as.list(sys.call()[[2]][-1]),as.character) getFunNames(c(mean,ff)) [1] mean ff foo - function() getFunNames(c(mean,ff)) foo() [1] mean ff HTH, Jerome -- Jerome Asselin, M.Sc., Agent de recherche, RHCE CHUM -- Centre de recherche 3875 rue St-Urbain, 3e etage // Montreal QC H2W 1V1 Tel.: 514-890-8000 Poste 15914; Fax: 514-412-7106 __ 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] Writing Text into Plots
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 Duncan Murdoch wrote: On 10/5/2006 4:10 PM, Lothar Botelho-Machado wrote: -BEGIN PGP SIGNED MESSAGE- Hash: SHA1 Hello, I have a simple question on the text() method in plots, e.g.: text(3,4,adj=1,cex=1.0, expression(alpha == beta)) I know there exists a lot more like frac(), etc which could be used for expression. But a help(frac) doesn't return any results - where do I have to look for a documentation of possible text commands? More in detail I'm searching for how to set an index like mu_{1} and mu_{2} ? And how to get something like a bar over a character like \bar x ? See ?plotmath, and demo(plotmath). Duncan Murdoch Thanx a lot!! -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.5 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org iD8DBQFFJXsfHRf7N9c+X7sRAkfGAKCUMsskOqH5pDA+ToiQRDAEaxpw9gCghDjd X5wETK17s21+u255hQYcqWU= =1JXU -END PGP SIGNATURE- __ 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] The W statistic in wilcox.exact
Jue, On a second look, it appears that wilcox.test does report the offset-adjusted statistic U, as also mentioned in the help page. wilcox.test returns W=6 (instead of 12 as your example showed, unless wilcox_test is a different function). wilcox.test( 1:5 ~ c(1,1,0,0,0) )$statistic # or wilcox.test( 1:5 ~ factor(c(1,1,0,0,0)) )$statistic W 6 So there does not appear to be a difference between the two methods. Did I miss something? -Christos -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Christos Hatzis Sent: Thursday, October 05, 2006 4:41 PM To: [EMAIL PROTECTED]; r-help@stat.math.ethz.ch Subject: Re: [R] The W statistic in wilcox.exact Probably because of the offset: U = W - n*(n+1)/2 In your example, W=12 (=3+4+5) as reported by wilcox.test. The offset is 6 (=3*4/2) and therefore U=6. I am not certain as I haven't installed the exactRankTests package, but it seems that wilcox.exact reports U instead of W. -Christos Hatzis -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: Thursday, October 05, 2006 1:35 PM To: r-help@stat.math.ethz.ch Subject: [R] The W statistic in wilcox.exact Does anyone know why wilcox.exact gives W-statistic 6 instead of 12 as indicated below. 12 is the rank sum of group 0 of x, which is the linear statistic computed by wilcox_test. y-c(1,2,3,4,5) x-c(1,1,0,0,0) (a) wilcox.exact wilcox.exact(y~x) Exact Wilcoxon rank sum test data: y by x W = 6, p-value = 0.2 alternative hypothesis: true mu is not equal to 0 (b) wilcox_test tt-wilcox_test(y~factor(x),distribution=exact) statistic(tt,linear) 0 12 Jue Wang, Biostatistician Contracted Position for Preclinical Research Biostatistics PrO Unlimited (908) 231-3022 __ 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. __ 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. __ 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] which command
Are you sure Fall.dat is identical on both platforms, and that dt is also? Did you intend to ignore the cases where dt$V2 equals zero? Try, for example, table(which(dt[,2] == 0)) ## also ==1, ==2, ==3 unique( dt[which(dt[,2] == 0),5]## also ==1, ==2, ==3 unique( dt$V5 ) adding na.rm=TRUE to mean() Also, try this: tmp - lapply( split( dt$V5, dt$V2) , mean ) Favg - unlist(tmp) -Don At 4:57 PM -0400 10/5/06, AgusSusanto wrote: I obtained error messages when I run these commands in UNIX, but I obtained correct result when I run these command in WINDOWS. Can somebody point out the problem and give the solution. Thanks. dt-read.table(file=Fall.dat) dim(dt) [1] 19415 table(dt$V2) 0 1 2 3 220 989 639 93 Favg-as.matrix(c(1:max(dt$V2))) for(i in 1:max(dt$V2)){Favg[i]-mean(dt[which(dt[,2] == i),5])} Favg [,1] [1,] NA [2,] NA [3,] NA Warning messages: argument is not numeric or logical: returning NA (etc) __ 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. -- -- Don MacQueen Environmental Protection Department Lawrence Livermore National Laboratory Livermore, CA, USA __ 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.
[R] emacs-ess tab indentation
Hi, I have a 2 tab-identation in my emacs-ess and I would like to make it 8. Can somebody help me with this ? thanks Johan __ 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] which command
On Thu, 2006-10-05 at 16:57 -0400, AgusSusanto wrote: I obtained error messages when I run these commands in UNIX, but I obtained correct result when I run these command in WINDOWS. Can somebody point out the problem and give the solution. Thanks. Not without Fall.dat I suspect. please read the posting guide and act accordingly. (And why post twice?) I suspect some of your data are missing in Fall.dat (or end up missing because you didn't check what had been read in by R) and that this has propagated warnings in say mean, when called within na.rm = TRUE. But I'm guessing as you fail to provide a reproducible example. Why this is different between Unix/windows versions of R, I don't know, but are they running the same versions of R? G dt-read.table(file=Fall.dat) dim(dt) [1] 19415 table(dt$V2) 0 1 2 3 220 989 639 93 Favg-as.matrix(c(1:max(dt$V2))) for(i in 1:max(dt$V2)){Favg[i]-mean(dt[which(dt[,2] == i),5])} Favg [,1] [1,] NA [2,] NA [3,] NA Warning messages: argument is not numeric or logical: returning NA (etc) Cheers, @ __ 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. -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% *Note new Address and Fax and Telephone numbers from 10th April 2006* %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Gavin Simpson [t] +44 (0)20 7679 0522 ECRC [f] +44 (0)20 7679 0565 UCL Department of Geography Pearson Building [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street London, UK[w] http://www.ucl.ac.uk/~ucfagls/cv/ WC1E 6BT [w] http://www.ucl.ac.uk/~ucfagls/ %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ 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.
[R] [OT] testing for synchronicity
Greetings, friends in the R community, this is an OT question about statistics. Given four time series of events, what possibilities do I have to test for synchronicity? e.g. times - data.frame(year= c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10), event.1=c(1, 0, 0, 1, 2, 4, 1, 0, 0, 0), event.2=c(0, 0, 0, 1, 0, 2, 1, 0, 0, 0), event.3=c(1, 0, 0, 0, 1, 2, 4, 1, 0, 0), event.4=c(0, 1, 0, 1, 0, 0, 1, 1, 1, 0)) I have about 100 years of each, and my null hypothesis is that the events are not synchronous. Right now my thinking is to focus on pairwise comparisons: 1) ignore the magnitude and convert the series to binary 2) the sum of the product of the events of any two years is then the number of overlapping occurences. 3) I think that, in the absence of temporal autocorrelation, I could assume that this sum is hypergeometrically distributed. 4) I can test for statistical significance of this number by a moving blocks monte-carlo simulation. I will do this by taking blocks of contiguous years with a random start and reordering them randomly. This conditions on the number of event occurrences, which I would rather do than have it be random, and partially preserves the temporal autocorrelation. If anyone has any thoughts, or pointers, they'd be very welcome. Cheers Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 Email: [EMAIL PROTECTED] http://www.ms.unimelb.edu.au __ 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.
[R] Aggregate Values for All Levels of a Factor
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 Hello, I'm a novice user trying to figure out how to retain NA aggregate values. For example, given a data frame with data for 3 of the 4 possible factor colors(orange is omitted from the data frame), I want to calculate the average height by color, but I'd like to retain the knowledge that orange is a possible factor, its just missing. Here is the example code: data - data.frame(color = factor(c(blue,red,red,green,blue), levels = c(blue,red,green,orange)), height = c(2,8,4,4,5)) aggregate(data$height, list(color = data$color), mean) color x 1 blue 3.5 2 red 6.0 3 green 4.0 Instead I would like to get color x 1 blue 3.5 2red 6.0 3 green 4.0 4 orange NA Is this possible. I've read as much documentation as I can find, but am unable to find the solution. It seems like something people would need to do. So I would assume it must be built in somewhere or do I need to write my own version of aggregate? Thanks in advance, Kaom -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.5 (MingW32) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org iD8DBQFFJYrLaaZgZdCbWv4RApNoAJ9jqKXne3IlQnd+PprS+7Kz1l4oRACfeu5I Nv/xYWVsSGJD5+fdCP+02jk= =b5TI -END PGP SIGNATURE- __ 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.
[R] solaris 64 build?
We have a solaris/sparc machine that has been running an old version of R-devel: Version 2.2.0 Under development (unstable) (2005-06-04 r34577) which was built as m64 from sources. Attempting to upgrade to 2.4.0 the configure step goes ok, but I'm getting early on from make: gcc -m64 -L/opt/sfw/lib/sparcv9 -L/usr/lib/sparcv9 -L/usr/openwin/lib/sparcv9 -L/usr/local/lib -o R.bin Rmain.o CConverters.o CommandLineArgs.o Rdynload.o Renviron.o RNG.o apply.o arithmetic.o apse.o array.o attrib.o base.o bind.o builtin.o character.o coerce.o colors.o complex.o connections.o context.o cov.o cum.o dcf.o datetime.o debug.o deparse.o deriv.o dotcode.o dounzip.o dstruct.o duplicate.o engine.o envir.o errors.o eval.o format.o fourier.o gevents.o gram.o gram-ex.o graphics.o identical.o internet.o iosupport.o lapack.o list.o localecharset.o logic.o main.o mapply.o match.o memory.o model.o names.o objects.o optim.o optimize.o options.o par.o paste.o pcre.o platform.o plot.o plot3d.o plotmath.o print.o printarray.o printvector.o printutils.o qsort.o random.o regex.o registration.o relop.o rlocale.o saveload.o scan.o seq.o serialize.o size.o sort.o source.o split.o sprintf.o startup.o subassign.o subscript.o subset.o summary.o sysutils.o unique.o util.o version.o vfonts.o xxxpr.o mkdtemp.o ../unix/libunix.a ../appl/libappl.a ../nmath/libnmath.a -L../../lib -lRblas -L/usr/local/encap/gf7764-3.4.3+2/lib/gcc/sparc64-sun-solaris2.9/3.4.3 -L/usr/ccs/bin/sparcv9 -L/usr/ccs/bin -L/usr/ccs/lib -L/usr/local/encap/gf7764-3.4.3+2/lib/sparcv9 -L/usr/local/encap/gf7764-3.4.3+2/lib -lg2c -lm -lgcc_s ../extra/zlib/libz.a ../extra/bzip2/libbz2.a ../extra/pcre/libpcre.a ../extra/intl/libintl.a -lreadline -ltermcap -lnsl -lsocket -ldl -lm Undefined first referenced symbol in file __builtin_isnan arithmetic.o ld: fatal: Symbol referencing errors. No output written to R.bin collect2: ld returned 1 exit status I've tried to look at the difference in outcomes in the old R-devel version -- if I touch arithmetic.c there and then type make I get a something almost the same as above except for the following bits that are new to 2.4.0 (this diff is after replacing spaces with linebreaks obviously.) ysidro.econ.uiuc.edu% diff t0 t1 54a55 localecharset.o 81a83 rlocale.o 101a104 mkdtemp.o 104a108,109 -L../../lib -lRblas Has there been some change in the way that Rblas is used, or in isnan? It didn't seem so from a look at arithmetic.c, but this is well beyond me. I hope that someone sees something suspicious, or could point me toward a better diagnostic. Thanks, Roger url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 __ 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] xyplot
If you include a conditioning factor (here it is the variable groups) you will get more than one panel. To get loess curves for the two groups, the easiest way is to use the type argument: xyplot(y ~ time, groups=group,data=df,type=c(p,smooth)) On 05/10/06, Osman Al-Radi [EMAIL PROTECTED] wrote: Hello, Thanks for your answer, however, the resulting graph has two panels with a loess curve per subject. I need a single panel with a loess curve per group.. Osman On 10/5/06, Sundar Dorai-Raj [EMAIL PROTECTED] wrote: Osman Al-Radi said the following on 10/5/2006 3:43 PM: Hi, for the data below: time-c(rep(1:10,5)) y-time+rnorm(50,5,2) subject-c(rep('a',10),rep('b',10),rep('c',10),rep('d',10),rep('e',10)) group-c(rep('A',30),rep('B',20)) df-data.frame(subject,group,time,y) I'd like to produce a plot with a single pannel with two loess curves one for each group. the code below does that but on two different pannels, I'd like to have both curves on the same pannel with different colors. xyplot(y~time|group, groups=subject, panel=function(...){ panel.loess(...) } ,data=df) Simple: xyplot(y ~ time | group, data = df, groups = subject, panel = function(...) { panel.superpose(...) panel.superpose(panel.groups = panel.loess, ...) }) HTH, --sundar -- Osman O. Al-Radi, MD, MSc, FRCSC Fellow, Cardiovascular Surgery The Hospital for Sick Children University of Toronto, Canada [[alternative HTML version deleted]] __ 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. -- = David Barron Said Business School University of Oxford Park End Street Oxford OX1 1HP __ 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] Block comments in R?
Please let's not waste any time trying to figure out how to add block comments to R. In any guise they are highly error prone. Although its precursors (PL/I and Pascal) had block comments and did not have end-of-line comments, the programming language Ada was explicitly designed to have end-of-line comments and not block comments, as block comments were thought too easy to write incorrectly and too easy to misread. As just one example, the version of vim I have does syntax highlighting, BUT Vim doesn't parse the whole file, so if you have code in the middle of a block comment it may be highlighed incorrectly. Commenting code out and providing documentation comments are easily done with a good editor, although R documentation comments really belong in files where help() can find them. In the text editor I use, comment-region (using end-of-line comments) and uncomment-region (ditto) are just two keystrokes each; it's LESS effort than writing block comment delimiters would be. Also, when writing a big comment, the newline key automatically copies white space*(end of line comment mark white space*)? to the beginning of a new line, so again it is less typing to make a block comment using end-of-line comments than /* ... */ or anything like that, and of course fill-paragraph and justify-paragraph (two keystrokes each) can cope with the comment marks (a change from lines with comment marks to ones without, or vice versa, constitutes a paragraph break). __ 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] How to get the function names
Probably the best you can hope for is to cover most cases. This one uses match.call and handles a number of cases and perhaps if you spend more time on it might be able to add some cases where it fails such as the second L below: f - function(x) { if (!is.list(x)) x - list(x) if (is.null(names(x))) names(x) - names(x)[names(x) == ] - NA mc - match.call()[-1][[1]] if (length(mc) 1) mc - mc[-1] ifelse(is.na(names(x)), as.character(mc), names(x)) } f(c(a = mean)) f(list(a = mean, b = sd)) f(c(f = function(x)x*x)) f(list(f = function(x)x*x, function(x)1-x)) L - list(a = mean, b = sd) f(L) L - list(a = mean, function(x)x) f(L) f(mean) f(list(a = mean, sd)) f(list(mean, sd)) f(function(x)x*x) f(list(function(x)x*x, function(y)y-1)) On 10/5/06, Søren Højsgaard [EMAIL PROTECTED] wrote: I've defined the function getFunNames - function(FUN){ if (!is.list(FUN)) fun.names - paste(deparse(substitute(FUN)), collapse = ) else fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a))) fun.names } which gives what I want : getFunNames(mean) [1] mean getFunNames(ff) [1] ff getFunNames(c(mean,ff)) [1] mean ff If I call this within a function, things go wrong: 1] FUN foo(ff) [1] FUN foo(c(mean,ff)) Error in substitute(FUN)[-1] : object is not subsettable Obviously there are some things (quite a few things) which I have not understood. Can anyone help? Thanks Søren __ 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. __ 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] Aggregate Values for All Levels of a Factor
On Thu, 2006-10-05 at 15:44 -0700, Kaom Te wrote: -BEGIN PGP SIGNED MESSAGE- Hash: SHA1 Hello, I'm a novice user trying to figure out how to retain NA aggregate values. For example, given a data frame with data for 3 of the 4 possible factor colors(orange is omitted from the data frame), I want to calculate the average height by color, but I'd like to retain the knowledge that orange is a possible factor, its just missing. Here is the example code: data - data.frame(color = factor(c(blue,red,red,green,blue), levels = c(blue,red,green,orange)), height = c(2,8,4,4,5)) aggregate(data$height, list(color = data$color), mean) color x 1 blue 3.5 2 red 6.0 3 green 4.0 Instead I would like to get color x 1 blue 3.5 2red 6.0 3 green 4.0 4 orange NA Is this possible. I've read as much documentation as I can find, but am unable to find the solution. It seems like something people would need to do. So I would assume it must be built in somewhere or do I need to write my own version of aggregate? Thanks in advance, Kaom If you review the Details section of ?aggregate, you will note: Empty subsets are removed, ... Thus, one approach is: tmp - tapply(data$height, data$color, mean, na.rm = TRUE) tmp bluered green orange 3.56.04.0 NA DF - data.frame(color = names(tmp), mean.height = tmp, row.names = seq(along = tmp)) DF color mean.height 1 blue 3.5 2red 6.0 3 green 4.0 4 orange NA HTH, Marc Schwartz __ 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] How to get the function names
I should have mentioned is that the way it works is that it uses the name of the list component, if any, otherwise it uses the name of the function if its given as a number and otherwise it uses the function itself or possibly the name of the list. On 10/5/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: Probably the best you can hope for is to cover most cases. This one uses match.call and handles a number of cases and perhaps if you spend more time on it might be able to add some cases where it fails such as the second L below: f - function(x) { if (!is.list(x)) x - list(x) if (is.null(names(x))) names(x) - names(x)[names(x) == ] - NA mc - match.call()[-1][[1]] if (length(mc) 1) mc - mc[-1] ifelse(is.na(names(x)), as.character(mc), names(x)) } f(c(a = mean)) f(list(a = mean, b = sd)) f(c(f = function(x)x*x)) f(list(f = function(x)x*x, function(x)1-x)) L - list(a = mean, b = sd) f(L) L - list(a = mean, function(x)x) f(L) f(mean) f(list(a = mean, sd)) f(list(mean, sd)) f(function(x)x*x) f(list(function(x)x*x, function(y)y-1)) On 10/5/06, Søren Højsgaard [EMAIL PROTECTED] wrote: I've defined the function getFunNames - function(FUN){ if (!is.list(FUN)) fun.names - paste(deparse(substitute(FUN)), collapse = ) else fun.names - unlist(lapply(substitute(FUN)[-1], function(a) paste(a))) fun.names } which gives what I want : getFunNames(mean) [1] mean getFunNames(ff) [1] ff getFunNames(c(mean,ff)) [1] mean ff If I call this within a function, things go wrong: 1] FUN foo(ff) [1] FUN foo(c(mean,ff)) Error in substitute(FUN)[-1] : object is not subsettable Obviously there are some things (quite a few things) which I have not understood. Can anyone help? Thanks Søren __ 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. __ 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.