Re: [R] When calling external C-function repeatedly I get differentresults; can't figure out why..
Not an expert in programming either, but to me it seems like you've forgotten to initialize the variable tr. It just picks up garbage from allocated memory previously initialized by other processes. Med venlig hilsen Frede Aakmann Tøgersen -Oprindelig meddelelse- Fra: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] På vegne af Søren Højsgaard Sendt: 9. marts 2006 02:15 Til: r-help@stat.math.ethz.ch Emne: [R] When calling external C-function repeatedly I get differentresults; can't figure out why.. Dear all, I need to calculate tr(xyxy) fast for matrices x and y. Inspired by the R-source code, I've created the following functions (I am *new* to writing external C-functions, so feel free to laugh at my code - or, perhaps, suggest changes): #include Rinternals.h #include R_ext/Applic.h /* for dgemm */ static void matprod(double *x, int nrx, int ncx, double *y, int nry, int ncy, double *z) { char *transa = N, *transb = N; double one = 1.0, zero = 0.0; F77_CALL(dgemm)(transa, transb, nrx, ncy, ncx, one, x, nrx, y, nry, zero, z, nrx); } SEXP trProd2(SEXP x, SEXP y) { int nrx, ncx, nry, ncy, mode, i; SEXP xdims, ydims, ans, ans2, tr; xdims = getAttrib(x, R_DimSymbol); ydims = getAttrib(y, R_DimSymbol); mode = REALSXP; nrx = INTEGER(xdims)[0]; ncx = INTEGER(xdims)[1]; nry = INTEGER(ydims)[0]; ncy = INTEGER(ydims)[1]; PROTECT(ans = allocMatrix(mode, nrx, ncy)); PROTECT(ans2 = allocMatrix(mode, nrx, ncy)); PROTECT(tr = allocVector(mode, 1)); matprod(REAL(x), nrx, ncx, REAL(y), nry, ncy, REAL(ans)); matprod(REAL(ans), nrx, ncy, REAL(ans), nrx, ncy, REAL(ans2)); for (i=0; i nrx; i++){ REAL(tr)[0] = REAL(tr)[0] + REAL(ans2)[i*(nrx+1)]; } UNPROTECT(3); return(tr); } In R I do: trProd2 - function(x, y) { .Call(trProd2, x, y) } x - y - matrix(1:4,ncol=2) storage.mode(x) - storage.mode(y) - double for (i in 1:10) print(trProd2(x, y)) [1] 835 [1] 833 [1] 833 [1] 862 [1] 834 [1] 835 [1] 834 [1] 836 [1] 862 [1] 833 The correct answer is 833. Can anyone give me a hint to what is wrong? I am on a windows xp platform. Thanks in advance 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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] problem installing RNetCDF
You forgot to tell us you are running MacOS X. Looks like a problem with your compiler installation (and are you using a pre-compiled R from CRAN and an incompatible compiler?). The installation procedure is complicated because the system dependencies need to be checked. It works well under Unix and not well under Windows. Some MacOS X users seem to think MacOS X is 'Unix-alike', but it vies with AIX as being the most extreme outlier in that set. (For a start, it is based on BSD, not the Single Unix specification, and no other (recent) Unix-alike has bundles.) On Wed, 8 Mar 2006, Zepu Zhang wrote: Hello all, I set 'UDUNITS_PATH' and 'NETCDF_PATH' successfully to my custom places and then % R CMD INSTALL RNetCDF_1.1-3.tar.gz and got this: ... checking for executable suffix... checking for object suffix... o checking whether we are using the GNU C compiler... yes checking whether gcc accepts -g... yes checking for main in -lnetcdf... yes checking for main in -ludunits... yes checking for /Users/zpzhang/research/library/include/netcdf.h... yes checking for /Users/zpzhang/research/library/include/udunits.h... yes configure: creating ./config.status config.status: creating R/load.R config.status: creating src/Makevars ** libs gcc -no-cpp-precomp -I/sw/Library/Frameworks/R.framework/Resources/ include -I/Users/zpzhang/research/library/include -I/Users/zpzhang/research/ library/include -I/sw/include -fno-common -g -O2 -c RNetCDF.c -o RNetCDF.o gcc: unrecognized option '-no-cpp-precomp' gcc -bundle -flat_namespace -undefined suppress -lcc_dynamic -L/sw/lib -o RNetCDF.so RNetCDF.o -L/Users/zpzhang/research/library/lib -lnetcdf -L/ Users/zpzhang/research/library/lib -ludunits -lcc_dynamic -F/sw/Library/ Frameworks -framework R gcc: couldn't run 'undle-gcc-4.0.2': No such file or directory make: *** [RNetCDF.so] Error 1 ERROR: compilation failed for package 'RNetCDF' Does anyone have any insight to this 'undle-gcc-4.02' problem? By the way, the package contains one C source file only but the installation seems quite complicated. I understand it uses 'standard' procedures--- configure, etc.---but wouldn't a simple Makefile be much easier to cope with? I would be pretty straightfoward what variables and paths to change, etc. By the way again, any comments on ncdf vs RNetCDF? I just couldn't like the inconsistent naming in ncdf (put.var.ncdf, att.put.ncdf, dim.get.ncdf, get.var.ncdf, ...), otherwise its very recent update would have attracted me to it. Anyway these are not big packages and should be all fine. Thanks for the help! __ 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 -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] memory limits in Windows
Is your OS set to allow 3GB users address space for processes? Since you are running a 32-bit executable under `Windows XP-64' I don't know if that is possible or how to do it, so please ask whoever advised you to buy that OS. You might want to send some sample code to the mgcv maintainer to see if he has any suggestions for how to do this more frugally. On Wed, 8 Mar 2006, [EMAIL PROTECTED] wrote: Hello, I apologize, since I know that variations of this question come up often, but I have searched both the archives and the rw-FAQ, and haven't had any luck solving my problem. I am trying to run some generalized additive mixed models (using the mgcv library) with some large datasets, and get error messages regarding memory allocation. An example dataset has around 45,000 points contributed by about 2200 individuals (so about 20 observations per individual). Errors within individuals are modeled as AR(1). Currently, I can run the models on a random subset of about 25% of the data in a reasonable amount of time, so I think that memory is the only major issue here. I have used the --max-mem-size command line option, and set it to the maximum allowable, which is 3Gb (any larger and I get a message that it is too large and is ignored when I open R). I also run the R command memory.limit(4095), again the maximum allowable without receiving a don't be silly message. I am running this on a brand new computer with Windows XP-64 OS and 4GB RAM (it was bought primarily to be able to handle these models!) Do I have any options to increase memory? Any advice would be hugely appreciated. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to plot the xaxis label at 45 degree angle?
Hello, Check under par() the option srt and las. Although the way to rotate the labels is via manipulating the axis() of your plot separately. Regards, Carlos Ortega. On 3/8/06, Lisa Wang [EMAIL PROTECTED] wrote: Hello there, I would like to plot a graph with the x axis's label displayed at a 45 angle to the x axis instead of horizontal to it as the label is very long. What should I do? Thank you for your help in advance Lisa Wang Princess Margaret Hospital Toronto, Ca __ 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 [[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
Re: [R] Debugging a program written in the R language
On Wed, 8 Mar 2006, Liaw, Andy wrote: [...] Now, Roger Peng has kindly made available a note on debugging R code: http://www.biostat.jhsph.edu/~rpeng/docs/R-debug-tools.pdf There is also a new chapter in the `Writing R Extensions' manual in the R-devel version of R (to be R 2.3.0). So now is a very good time for people to offer further information/suggestions for that chapter. Also, for compiled code linked to R, Duncan Murdoch's notes will be helpful: http://www.stats.uwo.ca/faculty/murdoch/software/debuggingR/ That's Windows-specific: the new chapter is mainly Unix-based. [...] -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] coding problems
Please look at the help page for arima.sim(). You first argument is being matched to 'rand.gen', and is not a function. I suspect you meant to call arima(). On Tue, 7 Mar 2006, (s) Richard Nuttall wrote: hi, I am trying to fit an ARIMA model to some time series data, I have used differencing to make the data stationary. dailyibm is the data I am using, could someone please help out in identifying the coding as I can't seem to identify the problem many thanks fit.ma - arima.sim(dailyibm - mean(dailyibm), model=list(order=c(0,1,0)),n = ) Error in inherits(x, data.frame) : couldn't find function rand.gen __ 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 -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Making an S3 object act like a data.frame
On Tue, 7 Mar 2006, hadley wickham wrote: I tend to have to use trial and error myself. Here is another possibility. That's got the subsetting solved, so here's the next challenge lm(x ~ y, z) Error in as.data.frame.default(data) : cannot coerce class myobj into a data.frame as.data.frame.myobj - function(x) x[[1]] lm(x ~ y, z) Error in eval(expr, envir, enclos) : numeric 'envir' arg not of length one I'm guessing this is pretty much impossible to get around, because there is no way to tell eval how to deal with myobj type objects, and lm only dispatches based on the type of the first argument. Did you write an as.data.frame method? From ?model.frame data: 'data.frame', list, 'environment' or object coercible to 'data.frame' containing the variables in 'formula'. Neither a matrix nor an array will be accepted. so I believe that an as.data.frame method is all that is required. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] 'less' for R?
you can also use str(data.frame) which allows you to access the structure of your object. Florence. On 3/8/06, Federico Calboli [EMAIL PROTECTED] wrote: Hi All, is there an equivalent of the Unix command 'less' (or 'more'), so I can look at what's inside a data.frame or a matrix without having it printed out on console? I am using R on Debian Linux and Mac OS 10.4.5 Cheers, F -- Federico C. F. Calboli Department of Epidemiology and Public Health Imperial College, St Mary's Campus Norfolk Place, London W2 1PG Tel +44 (0)20 7594 1602 Fax (+44) 020 7594 3193 f.calboli [.a.t] imperial.ac.uk f.calboli [.a.t] gmail.com __ 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 [[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
[R] Multivariate Autoregressive Model calibration and residual testing
Hi, I am using the mAr package to calibrate an Multivariate model (size 3, order 12). I am trying to do the two following things: 1. I would like to calibrate the model using not a single time series, but several of them: each time series should be seen as one independent realisation of the mAr process; for instance this happens when you have a time series with lacking data (holes) : if I have two time series and a hiatus in the middle, I would like my obective function in the calibration to be the sum of the two classical objective function, one for each continuous time series. Is there a function in R that could do this type of estimation? Or is there a nice way to do that by modifying the objective function within the mAr.est function (I tried but the code is too opaque..)? 2. Second, I would like to use the multivariate Ljung-Box test for the residuals of this multivariate AR, but I couldn't find it - (I found just the univariate case) - does it exist somewhere? Thanks a lot for your help!! Stephane *** Bear Stearns is not responsible for any recommendation, solicitation, offer or agreement or any information about any transaction, customer account or account activity contained in this communication. *** __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] how to make some characters in the xlab of a plot to be superscripted or subscript?
Dear R-users, I am trying to make my plot nicer. Does anyboby know how to make some characters in the xlab of a plot to be superscripted or subscript? ? for example m3/s in the x-axis label, how to make 3 supersripted? Best,Jing __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Linking Rblas
Dear all, when making a DLL via Rcmd SHLIB is there a way to link against a library such as Rblas (I am on a Windows platform) on a case to case basis? I played a bit around with some self-written C-code which should call the function 'dasum' defined in Blas.h. I encountered the following problem: u:\codeplayRcmd SHLIB mysumming.c Rcmd SHLIB mysumming.c gcc --shared -s -o mysumming.dll mysumming.def mysumming.a -LU:/R/R-2.2.0/src/gnuwin32 -lg2c -lR mysumming.a(mysumming.o)(.text+0x1b):mysumming.c: undefined reference to `dasum_' collect2: ld returned 1 exit status make: *** [mysumming.dll] Error 1 My (very limited) knowledge on C tells me that the function 'dasum' is actually defined (due to #include R_ext/BLAS.h in my C-file) but it needs to be linked against the correct library during the linking procedure. I guess the correct library should be Rblas.dll since I could do it by hand calling gcc like this: gcc --shared -s -o mysumming.dll mysumming.def mysumming.a -LU:/R/R-2.2.0/src/gnuwin32 -lg2c -lR -lRblas and it worked. Thank you very much for your help in advance. If you need more details such as the C-Code, the R-Code calling C, ..: just let me know. Thanks, Roland version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.0 year 2005 month10 day 06 svn rev 35749 language R u:\codeplaygcc --version gcc --version gcc (GCC) 3.4.2 (mingw-special) Copyright (C) 2004 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + This mail has been sent through the MPI for Demographic Rese...{{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
[R] R + calculations on clusters
Dear community, if I install R locally as some user on a server without root rights, will R be cabable of performing calculations on a cluster spreading the job on several cpu's and nodes? Or will R just simulate that it is doing so while just using one cpu in reality? regards Benjamin __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Linking Rblas [under Windows]
Yes, lots of packages link to Rblas. The help page for SHLIB tells you Please consult section 'Creating shared objects' in the manual 'Writing R Extensions' for how to customize it (for example to add 'cpp' flags and to add libraries to the link step) and for details of some of its quirks. so it should not have been at all hard to find the information. A portable way would be to have a file Makevars containing PKG_LIBS=$(BLAS_LIBS) Note that your precise example will not work in R-2.3.0 when it is released, whereas the portable form will. (In fact it only works in R-2.2.0 if you have previously make libRblas.a.) On Thu, 9 Mar 2006, Rau, Roland wrote: Dear all, when making a DLL via Rcmd SHLIB is there a way to link against a library such as Rblas (I am on a Windows platform) on a case to case basis? I played a bit around with some self-written C-code which should call the function 'dasum' defined in Blas.h. I encountered the following problem: u:\codeplayRcmd SHLIB mysumming.c Rcmd SHLIB mysumming.c gcc --shared -s -o mysumming.dll mysumming.def mysumming.a -LU:/R/R-2.2.0/src/gnuwin32 -lg2c -lR mysumming.a(mysumming.o)(.text+0x1b):mysumming.c: undefined reference to `dasum_' collect2: ld returned 1 exit status make: *** [mysumming.dll] Error 1 My (very limited) knowledge on C tells me that the function 'dasum' is actually defined (due to #include R_ext/BLAS.h in my C-file) but it needs to be linked against the correct library during the linking procedure. I guess the correct library should be Rblas.dll since I could do it by hand calling gcc like this: gcc --shared -s -o mysumming.dll mysumming.def mysumming.a -LU:/R/R-2.2.0/src/gnuwin32 -lg2c -lR -lRblas and it worked. Thank you very much for your help in advance. If you need more details such as the C-Code, the R-Code calling C, ..: just let me know. Thanks, Roland version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.0 year 2005 month10 day 06 svn rev 35749 language R u:\codeplaygcc --version gcc --version gcc (GCC) 3.4.2 (mingw-special) Copyright (C) 2004 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + This mail has been sent through the MPI for Demographic Rese...{{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 -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] POSIX time zone codes
On Tue, 7 Mar 2006, Jason Horn wrote: Thank you again Gabor, that did the trick. Any thoughts on where I go go for a reference for these time codes? Where did you get CDT6CST from? Or is this just one of those things that is common knowledge in UNIX circles. To the R developers: I recommend a sentence be added to the manual for as.POSIXxx such as vales for tz are system dependent, examples for common systems are. You mean repeating what is already there: tz: A timezone specification to be used for the conversion, _if one is required_. System-specific, but '' is the current timezone, and 'GMT' is UTC (Coordinated Universal Time, in French). ? What are `common systems'? If you have access to such a set of systems and can provide a comprehensive set of (tested) examples please provide a patch. The POSIX specification is at http://www.opengroup.org/onlinepubs/009695399/basedefs/xbd_chap08.html under TZ. You are expecting R to give you information about your OS, and that's an unreasonable expectation. Thanks! On Mar 7, 2006, at 9:00 AM, Gabor Grothendieck wrote: Only and GMT are really guaranteed to work on all systems since the time zones are system dependent but try: CDT6CST and see if that works on your system. On 3/7/06, Jason Horn [EMAIL PROTECTED] wrote: Whoops, [EDIT] as.POSIX(x, tz=UTC) ... works, gives UTC times as.POSIX(x, tz=EST) ... works, gives EST times as.POSIX(x, tz=CST) ... does NOT work, gives UTC times [/EDIT] On Mar 7, 2006, at 8:05 AM, Jason Horn wrote: as.POSIX(x, tz=UTC) ... works, gives UTC times as.POSIX(x, tz=UTC) ... works, gives EST times as.POSIX(x, tz=CST) ... does NOT work, gives UTC times -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R + calculations on clusters
On 3/9/06 5:58 AM, Benjamin Otto [EMAIL PROTECTED] wrote: Dear community, if I install R locally as some user on a server without root rights, will R be cabable of performing calculations on a cluster spreading the job on several cpu's and nodes? Or will R just simulate that it is doing so while just using one cpu in reality? Benjamin, R will not spread the job to different nodes on a cluster without you writing parallelized code (and this is not a function of where or how R is installed). Rmpi, RPVM, and the snow packages are useful for doing this. Sean __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Linking Rblas [under Windows]
Dear Prof. Ripley, thank you very much. It works now absolutely fine using a Makevars file with the contents you suggested. Thanks again, Roland -Original Message- From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] Sent: Thursday, March 09, 2006 12:09 PM To: Rau, Roland Cc: r-help@stat.math.ethz.ch Subject: Re: [R] Linking Rblas [under Windows] Yes, lots of packages link to Rblas. The help page for SHLIB tells you Please consult section 'Creating shared objects' in the manual 'Writing R Extensions' for how to customize it (for example to add 'cpp' flags and to add libraries to the link step) and for details of some of its quirks. so it should not have been at all hard to find the information. A portable way would be to have a file Makevars containing PKG_LIBS=$(BLAS_LIBS) Note that your precise example will not work in R-2.3.0 when it is released, whereas the portable form will. (In fact it only works in R-2.2.0 if you have previously make libRblas.a.) On Thu, 9 Mar 2006, Rau, Roland wrote: Dear all, when making a DLL via Rcmd SHLIB is there a way to link against a library such as Rblas (I am on a Windows platform) on a case to case basis? I played a bit around with some self-written C-code which should call the function 'dasum' defined in Blas.h. I encountered the following problem: u:\codeplayRcmd SHLIB mysumming.c Rcmd SHLIB mysumming.c gcc --shared -s -o mysumming.dll mysumming.def mysumming.a -LU:/R/R-2.2.0/src/gnuwin32 -lg2c -lR mysumming.a(mysumming.o)(.text+0x1b):mysumming.c: undefined reference to `dasum_' collect2: ld returned 1 exit status make: *** [mysumming.dll] Error 1 My (very limited) knowledge on C tells me that the function 'dasum' is actually defined (due to #include R_ext/BLAS.h in my C-file) but it needs to be linked against the correct library during the linking procedure. I guess the correct library should be Rblas.dll since I could do it by hand calling gcc like this: gcc --shared -s -o mysumming.dll mysumming.def mysumming.a -LU:/R/R-2.2.0/src/gnuwin32 -lg2c -lR -lRblas and it worked. Thank you very much for your help in advance. If you need more details such as the C-Code, the R-Code calling C, ..: just let me know. Thanks, Roland version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.0 year 2005 month10 day 06 svn rev 35749 language R u:\codeplaygcc --version gcc --version gcc (GCC) 3.4.2 (mingw-special) Copyright (C) 2004 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + This mail has been sent through the MPI for Demographic Rese...{{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 -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 + This mail has been sent through the MPI for Demographic Rese...{{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
Re: [R] how to make some characters in the xlab of a plot to be superscripted or subscript?
Jing Yang wrote: Dear R-users, I am trying to make my plot nicer. Does anyboby know how to make some characters in the xlab of a plot to be superscripted or subscript? ? for example m3/s in the x-axis label, how to make 3 supersripted? Best,Jing See ?plotmath. plot(1:10, xlab = expression(m^3/5)) or plot(1:10, xlab = expression(over(m^3, 5))) 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
Re: [R] how to make some characters in the xlab of a plot to be superscripted or subscript?
Hi see ?plotmath # not tested expression(paste(m^3, /s)) but I do not know how / is printed. HTH Petr On 9 Mar 2006 at 11:48, Jing Yang wrote: Date sent: Thu, 9 Mar 2006 11:48:45 +0100 From: Jing Yang [EMAIL PROTECTED] To: r-help r-help@stat.math.ethz.ch Subject:[R] how to make some characters in the xlab of a plot to be superscripted or subscript? Send reply to: [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] mailto:[EMAIL PROTECTED] Dear R-users, I am trying to make my plot nicer. Does anyboby know how to make some characters in the xlab of a plot to be superscripted or subscript? ? for example m3/s in the x-axis label, how to make 3 supersripted? Best,Jing Petr Pikal [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Error (Anova)
Hello I realized Anova repeated mesur And I have result following: summary(a) Error: indiv Df Sum Sq Mean Sq F value Pr(F) Alt 2 28.377 14.188 63.5236 2e-16 *** Site 9 75.238 8.360 37.4279 2e-16 *** Espèce1 1.002 1.002 4.4867 0.03514 * Competition 1 0.228 0.228 1.0202 0.31344 Residuals 252 56.286 0.223 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error: indiv:time Df Sum Sq Mean Sq F valuePr(F) time 3 68.328 22.776 762.3477 2.2e-16 *** Alt:time 6 1.513 0.252 8.4387 6.998e-09 *** time:Site 27 6.872 0.255 8.5192 2.2e-16 *** time:Espèce3 7.403 2.468 82.5979 2.2e-16 *** time:Competition 3 0.121 0.040 1.35360.2559 Residuals756 22.586 0.030 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 But I do not know ( Error)? - [[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
Re: [R] power and sample size for a GLM with Poisson response variable
The problem with the simulation is with the second group where there is a high probability of obtaining all zeroes for the sample and this is causing problems with the Wald statistics you are using to check for a difference. Here is an example. Browse[1] summary(res) Call: glm(formula = y[i, ] ~ group, family = poisson()) Deviance Residuals: Min 1Q Median 3Q Max -1.290e-01 -1.290e-01 -1.231e-05 -1.231e-05 2.756e+00 Coefficients: Estimate Std. Error z value Pr(|z|) (Intercept) -4.7884 0. -14.365 2e-16 *** group2 -18.5142 1235.1829 -0.0150.988 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 110.881 on 4260 degrees of freedom Residual deviance: 86.192 on 4259 degrees of freedom AIC: 108.19 Number of Fisher Scoring iterations: 21 On Wed, 8 Mar 2006 Wassell, James T., Ph.D. wrote: Subject: Re: [R] power and sample size for a GLM with Poisson response variable To: r-help@stat.math.ethz.ch Cc: Craig A. Faulhaber [EMAIL PROTECTED] Message-ID: [EMAIL PROTECTED] Content-Type: text/plain; charset=US-ASCII Craig, Thanks for your follow-up note on using the asypow package. My problem was not only constructing the constraints vector but, for my particular situation (Poisson regression, two groups, sample sizes of (1081,3180), I get very different results using asypow package compared to my other (home grown) approaches. library(asypow) pois.mean-c(0.0065,0.0003) info.pois - info.poisson.kgroup(pois.mean, group.size=c(1081,3180)) constraints - matrix(c(2,1,2), ncol = 3, byrow=T) poisson.object - asypow.noncent(pois.mean, info.pois,constraints) power.pois - asypow.power(poisson.object, c(1081,3180), 0.05) print(power.pois) [1] 0.2438309 asy.pwr() # the function is shown below. $power [1] 0.96 sim.pwr() # the function is shown below. 4261000 Poisson random variates simulated $power [1] 0.567 I tend to think the true power is greater than 0.567 but less than 0.96 (not as small as 0.2438309). Maybe I am still not using the asypow functions correctly? Suggestions/comments welcome. -Original Message- From: Craig A. Faulhaber [mailto:[EMAIL PROTECTED] Sent: Saturday, March 04, 2006 12:04 PM To: Wassell, James T., Ph.D. Subject: Re: [R] power and sample size for a GLM with Poisson response variable Hi James, Thanks again for your help. With the assistance of a statistician colleauge of mine, I figured out the constraints vector. It is a 3-column vector describing the null hypothesis. For example, let's say you have 3 populations and your null hypothesis is that there is no difference between the 3. The first row of the matrix would be 3 1 2 indicating that you have 3 populations and that populations 1 and 2 are equal. The second row would be 3 2 3 indicating that you have 3 populations and that populations 2 and 3 are equal. If you had only 2 populations, there would be only one row (2 1 2, indicating 2 populations with 1 and 2 equal). I hope this helps. Craig Wassell, James T., Ph.D. wrote: Craig, I found the package asypow difficult to use and it did not yield results in the ballpark of other approaches. (could not figure out the constraints vector). I wrote some simple functions, one asy.pwr uses the non-central chi-square distn. asy.pwr-function(counts=c(7,1),py=c(1081,3180)) { # a two group Poisson regression power computation # py is person-years or person-time however measured group-gl(2,1) ncp-summary(glm(counts~group+offset(log(py)),family=poisson)) $null.de viance q.tile-qchisq(.95,1) # actually just the X2 critical value of 3.841459 list(power=round(1-pchisq(q.tile,df=1,ncp),2))} The second function, sim.pwr, estimates power using simulated Poisson random variates. The most time consuming step is the call to rpois. (Maybe someone knows a more efficient way to accomplish this?). The for loop is rather quick in comparison. I hope you may find this helpful, or if you have solved your problem some other way, please pass along your approach.Note, that for this problem, very small values of lambda, the two approaches give much different power estimates (96% vs. 55% or so). My problem may be better addressed as binomial logistic regression, maybe then the simulation and the asymptotic estimates my agree better. sim.pwr-function(means=c(0.0065,0.0003),ptime=c (1081,3180),nsim=1000) { # a two group poisson regression power computation # based simulating lots of Poisson r.v.'s # input rates followed by a vector of the corresponding person times # the most time consuming part is the r.v. generation. # power is determined by counting the how often p-values are = 0.05 group-as.factor(rep(c(1,2),ptime)) rej-vector(length=nsim)
[R] how to use the randomForest and rpart function?
Michael - I recall reading something Breiman wrote that said essentially don't skimp on the number of trees - they are cheap to build and it makes for a better model. Also, look at your error rates (using plot), and make sure you run enough trees so that the error settles down. You'll likely be building 1000 or so trees. Tim Hi Andy, Does the randomForest have a Cross Validation built-in to decide what is the best number of trees or I have to find the best number manually by myself? Thanks a lot! Michael. On 3/7/06, Liaw, Andy [EMAIL PROTECTED] wrote: Yes, I do know. That's why I pointed you to the reference linked from the help page. BTW, there's also an R News article describing the initial version of the package. Have you perused that? Andy __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Identifying or searching for labels in a hclust/dendrogram/heatmap
Hi Sorry if this is in the help :-S I've looked at example(dendrogram) and though it gives some indication of what I want, it doesn't do all. OK, so here is what I want to do: draw a tree, and then have an action, on user-click, to either draw a sub tree or a plot of the data. I also want users to be able to search for a particular label and have it highlighted on the tree, say in red, where all the other labels are black. Now, the only tree I can use the user-click with is an hclust object, with the identify.hclust() function. As far as I know, neither the dendrogram objects or the output of heatmap is cut-able in this way. So I have the sub-tree and plot drawing set up on user-click using an hclust object and identify.hclust - good :) Therefore I am working with hclust objects (which is a shame as the dendrogram and heatmap objects look prettier) but I can cope. How do I then go on and highlight a single label, or group of labels, when I have plotted the hclust object? Can I highlight labels on dendrograms and heatmaps 2? I can just imagine that after a user has drawn the tree, they will want to know where their genes of interest actually lie within that tree. I'm happy to share the completed code when I have sorted all of this out Many thanks Mick __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Numerical Mathematics Consortium
Dear all, I found a link to the Numerical Mathematics Consortium http://www.nmconsortium.org/index.aspx Their rationale can be summarized (as far as I understood) by the paragraphs which I copied from their homepage and pasted below. Maybe this could be of interest also for the development of R? Best, Roland And here the copy+paste from their homepage: As modern numerical programming languages have replaced FORTRAN as the language of choice, and as new interactive design and development tools were introduced over the years, the standardized representation of common functions began to disperse again. The problem is as simple as looking at the Fast Fourier Transform (FFT) function: Mathcad t := 0..127 x := sin(2.p..01.t) w :=fft(x) MATLAB t = 0:127 x = sin(2*pi*.01*t) w = fft(x) Size of w: Mathcad = 65 MATLAB = 128 As in this case, even if the syntax is identical, the calculations yield different results! By definition, the FFT function of real data produces a result where the second half of the data is the complex conjugate of the first half. Mathcad uses this and only returns the first half of the data (i.e. N/2 + 1). MATLAB returns all of the data (N). Both are valid interpretations. The Numerical Mathematics Consortium is proposing a semantic standard for core numerical functions so as to help solve issues such as this. + This mail has been sent through the MPI for Demographic Rese...{{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
Re: [R] Identifying or searching for labels in a hclust/dendrogram/heatmap
On 3/9/06 7:47 AM, michael watson (IAH-C) [EMAIL PROTECTED] wrote: Hi Sorry if this is in the help :-S I've looked at example(dendrogram) and though it gives some indication of what I want, it doesn't do all. OK, so here is what I want to do: draw a tree, and then have an action, on user-click, to either draw a sub tree or a plot of the data. I also want users to be able to search for a particular label and have it highlighted on the tree, say in red, where all the other labels are black. Now, the only tree I can use the user-click with is an hclust object, with the identify.hclust() function. As far as I know, neither the dendrogram objects or the output of heatmap is cut-able in this way. So I have the sub-tree and plot drawing set up on user-click using an hclust object and identify.hclust - good :) Therefore I am working with hclust objects (which is a shame as the dendrogram and heatmap objects look prettier) but I can cope. How do I then go on and highlight a single label, or group of labels, when I have plotted the hclust object? Can I highlight labels on dendrograms and heatmaps 2? I can just imagine that after a user has drawn the tree, they will want to know where their genes of interest actually lie within that tree. I'm happy to share the completed code when I have sorted all of this out Mick, I personally would make use of the ctc package (bioconductor) to output data to Treeview (and/or cluster). R may be able to do this, but interaction with graphics still leaves a bit to be desired in R (at least as compared to a dedicated, specialized Java App), particularly once you get beyond about 50 genes. Treeview is just much better set up for this kinda stuff. You could probably automate the process given, for example, a limma MAList or fit object, or something like that. That said, I'd love to see what other people come up with to accomplish what you are asking. Sean __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Making an S3 object act like a data.frame
I'm guessing this is pretty much impossible to get around, because there is no way to tell eval how to deal with myobj type objects, and lm only dispatches based on the type of the first argument. Did you write an as.data.frame method? From ?model.frame My as.data.frame is : as.data.frame.ggobiDataset - function(x, ...) { as.data.frame(.GGobiCall(getData, x)) } (which in the context of my problem retrieves the data corresponding to x, an external pointer to a ggobi dataset) I still get the error reported above. Hadley data: 'data.frame', list, 'environment' or object coercible to 'data.frame' containing the variables in 'formula'. Neither a matrix nor an array will be accepted. so I believe that an as.data.frame method is all that is required. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] substituting values
dear list, i have a matrix with missing values like 1234x x2222 222x2 the x stands for the missing value. i have to substitute it to NA substitute or replace didn´t work out, as they are for vectors only however matrix[,i] also didn´t work. can anybody tell me, who i can change these values easily? thanks stefan [[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
Re: [R] substituting values
X - matrix(c( 1, 2, 3, 4, x, x, 2, 2, 2, 2, 2, 2, 2, x, 2), ncol=5, byrow=TRUE) X.new - matrix(as.numeric(X), ncol=ncol(X)) or X.new - apply(X, 2, function(x){replace(x, x == x, NA)}) Stefan Semmeling wrote: dear list, i have a matrix with missing values like 1234x x2222 222x2 the x stands for the missing value. i have to substitute it to NA substitute or replace didn´t work out, as they are for vectors only however matrix[,i] also didn´t work. can anybody tell me, who i can change these values easily? thanks stefan [[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 -- Chuck Cleland, Ph.D. NDRI, Inc. 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 452-1424 (M, W, F) fax: (917) 438-0894 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] substituting values
One way: A - matrix(round(rnorm(20)), ncol=4) A[2,1] - x A[which(A == x, arr.ind=TRUE)] - NA A - matrix(as.numeric(A), ncol=dim(A)[2]) Best, Matthias -Ursprüngliche Nachricht- Von: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Im Auftrag von Stefan Semmeling Gesendet: Donnerstag, 09. März 2006 14:32 An: r-help@stat.math.ethz.ch Betreff: [R] substituting values Vertraulichkeit: Persönlich dear list, i have a matrix with missing values like 1234x x2222 222x2 the x stands for the missing value. i have to substitute it to NA substitute or replace didn´t work out, as they are for vectors only however matrix[,i] also didn´t work. can anybody tell me, who i can change these values easily? thanks stefan [[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
Re: [R] POSIX time zone codes
On Wed, 8 Mar 2006, Don MacQueen wrote: I would have suggested US/Central, by analogy with US/Pacific which I use on my unix-like Mac OS X 10.3.9 system. I don't know what might be common knowledge in UNIX circles, but here's a bit of information from my system. It has a directory named /usr/share/zoneinfo. zoneinfo[166]% pwd /usr/share/zoneinfo zoneinfo[167]% ls | grep DT CST6CDT EST5EDT MST7MDT PST8PDT There are files named GMT, UTC, GB, and a number of others. There are numerous subdirectories, including for example, Asia, Australia, Europe, US. zoneinfo[168]% cd US US[169]% ls ./ AleutianEast-IndianaIndiana-Starke Pacific ../ Arizona Eastern MichiganPacific-New Alaska Central Hawaii MountainSamoa This suggests to me that valid values for the tz argument are relative paths to files /usr/share/zoneinfo, but I don't know that for a fact. That is documented as _one_ of the allowed forms on MacOS X, and 'man tzset' will show you others (I looked at online version for your OS). Strictly, that form should be preceded by a colon, that is ':US/Central'. But the information is only in 'man tzset' on some systems. According to http://www.twinsun.com/tz/tz-link.htm this is the case on most Unix-alikes, although the base path differs considerably. However, :US/Central is not valid on Windows, and in fact as far as I know Windows only allows you to set (via TZ) timezones with no DST or the US rules for DST. The Control Panel knows some other sets of rules. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] substituting values
On Thu, Mar 09, 2006 at 08:44:27AM -0500, Chuck Cleland wrote: ChuckX - matrix(c( 1, 2, 3, 4, x, Chuck x, 2, 2, 2, 2, Chuck 2, 2, 2, x, 2), Chuck ncol=5, byrow=TRUE) Why not: X[X==x] - NA mode(X) - numeric ## to get a numeric matrix Btw if you got X with read.table you can specify the missing string. ?read.table Ste Chuck ChuckX.new - matrix(as.numeric(X), ncol=ncol(X)) Chuck Chuckor Chuck ChuckX.new - apply(X, 2, function(x){replace(x, x == x, NA)}) Chuck ChuckStefan Semmeling wrote: Chuck dear list, Chuck Chuck i have a matrix with missing values like Chuck Chuck 1234x Chuck x2222 Chuck 222x2 Chuck Chuck the x stands for the missing value. Chuck i have to substitute it to NA Chuck Chuck substitute or replace didn?t work out, as they are for vectors only Chuck Chuck however matrix[,i] also didn?t work. Chuck Chuck can anybody tell me, who i can change these values easily? Chuck Chuck thanks stefan Chuck Chuck[[alternative HTML version deleted]] Chuck Chuck Chuck Chuck Chuck Chuck __ Chuck R-help@stat.math.ethz.ch mailing list Chuck https://stat.ethz.ch/mailman/listinfo/r-help Chuck PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Chuck Chuck-- ChuckChuck Cleland, Ph.D. ChuckNDRI, Inc. Chuck71 West 23rd Street, 8th floor ChuckNew York, NY 10010 Chucktel: (212) 845-4495 (Tu, Th) Chucktel: (732) 452-1424 (M, W, F) Chuckfax: (917) 438-0894 Chuck Chuck__ ChuckR-help@stat.math.ethz.ch mailing list Chuckhttps://stat.ethz.ch/mailman/listinfo/r-help ChuckPLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] lsa and Rstem?
Dear r-helpers, I can't get lsa to run because: library(lsa) Loading required package: Rstem Error in library(pkg, character.only = TRUE, logical = TRUE, lib.loc = lib.loc) : 'Rstem' is not a valid package -- installed 2.0.0? In addition: Warning message: cannot create HTML package index in: make.packages.html() install.packages('Rstem') Warning in install.packages(Rstem) : argument 'lib' is missing: using /Users/mk/Library/R/library Warning in download.packages(pkgs, destdir = tmpd, available = available, : no package 'Rstem' at the repositories Since Rstem is on http://www.omegahat.org/, I then followed the advice at http://www.omegahat.org/cranRepository.html: options(CRAN = c(getOption(CRAN), http://www.omegahat.org/R;)) install.packages('Rstem') Warning in install.packages(Rstem) : argument 'lib' is missing: using /Users/mk/Library/R/library Warning in download.packages(pkgs, destdir = tmpd, available = available, : no package 'Rstem' at the repositories Running platform powerpc-apple-darwin7.9.0 arch powerpc os darwin7.9.0 system powerpc, darwin7.9.0 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Debugging a program written in the R language
Thomas, This is what I use for writing R programs. I use it with Linux but adapting it for Windows shouldn't be that much of a problem for us MIT guys. In ~/.Rprofile (so it gets loaded every time I start R) I have myedit - function(object) { system(if [ ! -d $HOME/stat-misc/Rsrc ]; then mkdir $HOME/stat-misc/Rsrc; fi) system(paste(vim $HOME/stat-misc/Rsrc/,object,.R,sep=)) source(paste(~/stat-misc/Rsrc/,object,.R,sep=)) } To write a program called 'functionName' I do, within R, myedit('functionName'). This opens up an editor, vim in my case, for writing the function and saves the result in the ~/stat-misc/Rsrc directory and in the current R session. I can now repeat myedit('functionName') until I finally get it right. The function is also accessible with myedit() from any other R session at any time. Good luck, Robert Burrows, PhD New England Biometrics [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lsa and Rstem?
Michael-- I had emailed Duncan Lang about this a while back (the maintainer of Rstem). He gave me the following suggestion (which worked for me with R v2.2.1 and Windows XP): install.packages(Rstem, repos = http://www.omegahat.org/R;) cheers, Dave Dear r-helpers, I can't get lsa to run because: library(lsa) Loading required package: Rstem Error in library(pkg, character.only = TRUE, logical = TRUE, lib.loc = lib.loc) : 'Rstem' is not a valid package -- installed 2.0.0? In addition: Warning message: cannot create HTML package index in: make.packages.html() install.packages('Rstem') Warning in install.packages(Rstem) : argument 'lib' is missing: using /Users/mk/Library/R/library Warning in download.packages(pkgs, destdir = tmpd, available = available, : no package 'Rstem' at the repositories Since Rstem is on http://www.omegahat.org/, I then followed the advice at http://www.omegahat.org/cranRepository.html: options(CRAN = c(getOption(CRAN), http://www.omegahat.org/R;)) install.packages('Rstem') Warning in install.packages(Rstem) : argument 'lib' is missing: using /Users/mk/Library/R/library Warning in download.packages(pkgs, destdir = tmpd, available = available, : no package 'Rstem' at the repositories Running platform powerpc-apple-darwin7.9.0 arch powerpc os darwin7.9.0 system powerpc, darwin7.9.0 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R -- Dave Atkins, PhD [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Making an S3 object act like a data.frame
Can you create a small self contained reproducible example that does not work? The reproducible example I provided earlier on this thread worked fine. One idea is to check what the class is of the output of your .GGobiCall. If it were of class ggobiDataset then it would in turn be calling as.data.frame.ggobiDataset again. I mention this since the problem in your reproducible example was precisely of that sort. On 3/9/06, hadley wickham [EMAIL PROTECTED] wrote: I'm guessing this is pretty much impossible to get around, because there is no way to tell eval how to deal with myobj type objects, and lm only dispatches based on the type of the first argument. Did you write an as.data.frame method? From ?model.frame My as.data.frame is : as.data.frame.ggobiDataset - function(x, ...) { as.data.frame(.GGobiCall(getData, x)) } (which in the context of my problem retrieves the data corresponding to x, an external pointer to a ggobi dataset) I still get the error reported above. Hadley data: 'data.frame', list, 'environment' or object coercible to 'data.frame' containing the variables in 'formula'. Neither a matrix nor an array will be accepted. so I believe that an as.data.frame method is all that is required. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] help about TYPE III SS
Hi, I am very confused about those note that I found in this web page about TYPE III SS in ANOVA. In my mind Type III is not useful when there are interaction. Can someone help me ? web site : http://www.xlstat.com/fr/support/tutorials/ano2.htm __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] help about TYPE III SS
I suppose you're referring to the sentence (or whatever it corresponds to in the French version you cited): This means that the order in which the variables are selected will not have any effect on the values in the Type III SS. The Type III SS is generally the best method to use to interpret results when an interaction is part of the model. I would suggest that you take it up with XLSTAT support. Andy From: justin bem Hi, I am very confused about those note that I found in this web page about TYPE III SS in ANOVA. In my mind Type III is not useful when there are interaction. Can someone help me ? web site : http://www.xlstat.com/fr/support/tutorials/ano2.htm __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Identifying or searching for labels in a hclust/dendrogram/heatmap
Hi Sean Thanks for the help, but I really wanted to do this in R :) Any suggestions? Mick From: Sean Davis [mailto:[EMAIL PROTECTED] Sent: Thu 09/03/2006 1:01 PM To: michael watson (IAH-C); r-help Subject: Re: [R] Identifying or searching for labels in a hclust/dendrogram/heatmap On 3/9/06 7:47 AM, michael watson (IAH-C) [EMAIL PROTECTED] wrote: Hi Sorry if this is in the help :-S I've looked at example(dendrogram) and though it gives some indication of what I want, it doesn't do all. OK, so here is what I want to do: draw a tree, and then have an action, on user-click, to either draw a sub tree or a plot of the data. I also want users to be able to search for a particular label and have it highlighted on the tree, say in red, where all the other labels are black. Now, the only tree I can use the user-click with is an hclust object, with the identify.hclust() function. As far as I know, neither the dendrogram objects or the output of heatmap is cut-able in this way. So I have the sub-tree and plot drawing set up on user-click using an hclust object and identify.hclust - good :) Therefore I am working with hclust objects (which is a shame as the dendrogram and heatmap objects look prettier) but I can cope. How do I then go on and highlight a single label, or group of labels, when I have plotted the hclust object? Can I highlight labels on dendrograms and heatmaps 2? I can just imagine that after a user has drawn the tree, they will want to know where their genes of interest actually lie within that tree. I'm happy to share the completed code when I have sorted all of this out Mick, I personally would make use of the ctc package (bioconductor) to output data to Treeview (and/or cluster). R may be able to do this, but interaction with graphics still leaves a bit to be desired in R (at least as compared to a dedicated, specialized Java App), particularly once you get beyond about 50 genes. Treeview is just much better set up for this kinda stuff. You could probably automate the process given, for example, a limma MAList or fit object, or something like that. That said, I'd love to see what other people come up with to accomplish what you are asking. Sean __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Debugging a program written in the R language
Prof Brian Ripley [EMAIL PROTECTED] wrote in message news:[EMAIL PROTECTED] On Wed, 8 Mar 2006, Liaw, Andy wrote: So now is a very good time for people to offer further information/suggestions for that chapter. This may be useful: Mark Bravington, Debugging Without (Too Many) Tears, R News, Vol 3, No 3, Dec 2003, pp. 29-32. R-News: http://cran.r-project.org/doc/Rnews/ Exact link: http://cran.r-project.org/doc/Rnews/Rnews_2003-3.pdf efg __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] HCLUST subroutine question -- FORTRAN DO loops
Shown below is most of the FORTRAN subroutine named HCLUST. My question concerns the DO loop labeled as '10'. What happened to its CONTINUE statement? I will assume that after FLAG(I)=.TRUE. is executed that control returns to DO 10 I=1,N. Am I correct? Dave C Initializations C DO 10 I=1,N CWe do not initialize MEMBR in order to be able to restart the Calgorithm from a cut. CMEMBR(I)=1. 10 FLAG(I)=.TRUE. NCL=N C C Carry out an agglomeration - first create list of NNs C Note NN and DISNN are the nearest neighbour and its distance C TO THE RIGHT of I. C DO 30 I=1,N-1 DMIN=INF DO 20 J=I+1,N IND=IOFFST(N,I,J) IF (DISS(IND).GE.DMIN) GOTO 20 DMIN=DISS(IND) JM=J 20 CONTINUE NN(I)=JM DISNN(I)=DMIN 30 CONTINUE C 400 CONTINUE C Next, determine least diss. using list of NNs DMIN=INF DO 600 I=1,N-1 IF (.NOT.FLAG(I)) GOTO 600 IF (DISNN(I).GE.DMIN) GOTO 600 DMIN=DISNN(I) IM=I JM=NN(I) 600CONTINUE NCL=NCL-1 C C This allows an agglomeration to be carried out. C I2=MIN0(IM,JM) J2=MAX0(IM,JM) IA(N-NCL)=I2 IB(N-NCL)=J2 CRIT(N-NCL)=DMIN FLAG(J2)=.FALSE. C C Update dissimilarities from new cluster. C DMIN=INF DO 50 K=1,N IF (.NOT.FLAG(K)) GOTO 50 IF (K.EQ.I2) GOTO 50 IF (I2.LT.K) THEN IND1=IOFFST(N,I2,K) ELSE IND1=IOFFST(N,K,I2) ENDIF IF (J2.LT.K) THEN IND2=IOFFST(N,J2,K) ELSE IND2=IOFFST(N,K,J2) ENDIF IND3=IOFFST(N,I2,J2) D12=DISS(IND3) C C WARD'S MINIMUM VARIANCE METHOD - IOPT=1. C IF (IOPT.EQ.1) THEN DISS(IND1)=(MEMBR(I2)+MEMBR(K))*DISS(IND1)+ X (MEMBR(J2)+MEMBR(K))*DISS(IND2) - MEMBR(K)*D12 DISS(IND1)=DISS(IND1) / (MEMBR(I2)+MEMBR(J2)+MEMBR(K)) ENDIF C C SINGLE LINK METHOD - IOPT=2. C IF (IOPT.EQ.2) THEN DISS(IND1)=MIN(DISS(IND1),DISS(IND2)) ENDIF C C COMPLETE LINK METHOD - IOPT=3. C IF (IOPT.EQ.3) THEN DISS(IND1)=MAX(DISS(IND1),DISS(IND2)) ENDIF C C AVERAGE LINK (OR GROUP AVERAGE) METHOD - IOPT=4. C IF (IOPT.EQ.4) THEN DISS(IND1)=(MEMBR(I2)*DISS(IND1)+MEMBR(J2)*DISS(IND2))/ X (MEMBR(I2)+MEMBR(J2)) ENDIF C C MCQUITTY'S METHOD - IOPT=5. C IF (IOPT.EQ.5) THEN DISS(IND1)=0.5*DISS(IND1)+0.5*DISS(IND2) ENDIF C C MEDIAN (GOWER'S) METHOD - IOPT=6. C IF (IOPT.EQ.6) THEN DISS(IND1)=0.5*DISS(IND1)+0.5*DISS(IND2)-0.25*D12 ENDIF C C CENTROID METHOD - IOPT=7. C IF (IOPT.EQ.7) THEN DISS(IND1)=(MEMBR(I2)*DISS(IND1)+MEMBR(J2)*DISS(IND2)- X MEMBR(I2)*MEMBR(J2)*D12/(MEMBR(I2)+MEMBR(J2)))/ X (MEMBR(I2)+MEMBR(J2)) ENDIF C 50 CONTINUE MEMBR(I2)=MEMBR(I2)+MEMBR(J2) C C Update list of NNs C DO 900 I=1,N-1 IF (.NOT.FLAG(I)) GOTO 900 C(Redetermine NN of I:) DMIN=INF DO 870 J=I+1,N IF (.NOT.FLAG(J)) GOTO 870 IND=IOFFST(N,I,J) IF (DISS(IND).GE.DMIN) GOTO 870 DMIN=DISS(IND) JJ=J 870 CONTINUE NN(I)=JJ DISNN(I)=DMIN 900CONTINUE C C Repeat previous steps until N-1 agglomerations carried out. C IF (NCL.GT.1) GOTO 400 C C RETURN END C of HCLUST() __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] HCLUST subroutine question -- FORTRAN DO loops
David Emmith [EMAIL PROTECTED] writes: Shown below is most of the FORTRAN subroutine named HCLUST. My question concerns the DO loop labeled as '10'. What happened to its CONTINUE statement? I will assume that after FLAG(I)=.TRUE. is executed that control returns to DO 10 I=1,N. Am I correct? Yes. CONTINUE is purely used for legibility in such contexts; you can put the label on the last statement of the block instead (as Google might have told you soon enough...) Dave C Initializations C DO 10 I=1,N CWe do not initialize MEMBR in order to be able to restart the Calgorithm from a cut. CMEMBR(I)=1. 10 FLAG(I)=.TRUE. -- 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
Re: [R] Making an S3 object act like a data.frame
Can you create a small self contained reproducible example that does not work? The reproducible example I provided earlier on this thread worked fine. I wish I could, and I'm very grateful for the help, but because the data is an external pointer it's not easy to make a self-contained example. str.default(x) Classes ggobiDataset and `data.frame': 32 obs. of 1 variable: Classes 'ggobiDataset', 'data.frame' externalptr - attr(*, ggobi)=Class 'ggobi' externalptr getAnywhere(str.default)$objs[[1]](x) Classes 'ggobiDataset', 'data.frame' externalptr - attr(*, ggobi)=Class 'ggobi' externalptr One idea is to check what the class is of the output of your .GGobiCall. If it were of class ggobiDataset then it would in turn be calling as.data.frame.ggobiDataset again. I mention this since the problem in your reproducible example was precisely of that sort. class(x) [1] ggobiDataset data.frame class(as.data.frame(x)) [1] data.frame So I don't think that's the problem. Thanks again for your help. 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
[R] Call user written R functions in C
I have written some R functions and want to call them in C. Is there any way to do it? I have read through the Writing R Extensionsdocument. It seems that the only R functions that we can call in C is wrapped in Rmath.h or by embedding R.dll. Thanks. [[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
[R] count pixels of same color in pixmap object?
Dear all, I try to figure out how to use R to count the number of pixels of the same color in some gray-level picture. I managed to read it in either tiff or jpeg format, but the returned pixmap object keeps its information out of (my) reach. Is there an easy way to tabulate the different color/graylevel pixels and their numbers? Or should I use a completely different (free) software? Thanks for any hint, Christian. -- *** http://cognition.ups-tlse.fr/vas-y.php?id=chj [EMAIL PROTECTED] Christian Jost (PhD, MdC) Centre de Recherches sur la Cognition Animale Universite Paul Sabatier, Bat IV R3 118 route de Narbonne 31062 Toulouse cedex 9, France Tel: +33 5 61 55 64 37 Fax: +33 5 61 55 61 54 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] count pixels of same color in pixmap object?
On Thu, 9 Mar 2006, Christian Jost wrote: Dear all, I try to figure out how to use R to count the number of pixels of the same color in some gray-level picture. I managed to read it in either tiff or jpeg format, but the returned pixmap object keeps its information out of (my) reach. Is there an easy way to tabulate the different color/graylevel pixels and their numbers? Or should I use a completely different (free) software? The objects are new-style class objects, so they are documented internally - query by getSlots() or from ?pixmap-class: library(pixmap) x - read.pnm(system.file(pictures/logo.ppm, package=pixmap)[1]) z - as(x, pixmapGrey) class(z) getSlots(class(z)) # lets you know what the slots are greydata - slot(z, grey) # extract the slot object summary(c(greydata)) # c() to flatten a matrix length(unique(c(greydata))) length((c(greydata))) # rather a lot of different values, let's round: rgreydata - round(greydata*10) length(unique(c(rgreydata))) # looks good table(c(rgreydata)) I assume that you converted from tiff or jpeg to pnm outside pixmap. Thanks for any hint, Christian. -- 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
[R] Error with write.table
I'm having this error when I try to output out a gene expression data frame Error in file(file, ifelse(append, a, w)) : unable to open connection In addition: Warning message: cannot open file 'SAMPLES.txt' when issuing the command: write.table(SAMPLES, file = SAMPLES.txt) write.table(SAMPLES) will work fine, but when I specify a location it gives me the above error [[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
Re: [R] Error with write.table
mark salsburg [EMAIL PROTECTED] writes: I'm having this error when I try to output out a gene expression data frame Error in file(file, ifelse(append, a, w)) : unable to open connection In addition: Warning message: cannot open file 'SAMPLES.txt' when issuing the command: write.table(SAMPLES, file = SAMPLES.txt) write.table(SAMPLES) will work fine, but when I specify a location it gives me the above error The first guess would be that you do not have write permission to the working directory. If so, either obtain permissions or place the file in a directory where you can write: Use setwd() or give the full path to the file. -- 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
Re: [R] Making an S3 object act like a data.frame
I think the problem is that your ggobiDataset objects are also data.frame objects. They must NOT be. For example, note how the following example fails once we add data.frame to the class vector of x: The reason is that ordinary inheritance is not used by model.frame; rather, it uses is.data.frame to make the determination so by having data.frame in your class vector it does not try to convert it. # constructor myobj - function(...) +structure(list(value = data.frame(...)), class = myobj) $.myobj - function(obj, x) .subset2(obj, 1)[[x]] [[.myobj - function(obj, ...) .subset2(obj, 1)[[...]] [.myobj - function(obj, ...) .subset2(obj, 1)[...] as.data.frame.myobj - function(x) .subset2(x, 1) # test x - myobj(x = 1:3, y = 4:6) class(x) [1] myobj lm(y ~ x, x) Call: lm(formula = y ~ x, data = x) Coefficients: (Intercept)x 31 # now change class of x class(x) - c(myobj, data.frame) lm(y ~ x, x) Error in eval(expr, envir, enclos) : object y not found On 3/9/06, hadley wickham [EMAIL PROTECTED] wrote: Can you create a small self contained reproducible example that does not work? The reproducible example I provided earlier on this thread worked fine. I wish I could, and I'm very grateful for the help, but because the data is an external pointer it's not easy to make a self-contained example. str.default(x) Classes ggobiDataset and `data.frame': 32 obs. of 1 variable: Classes 'ggobiDataset', 'data.frame' externalptr - attr(*, ggobi)=Class 'ggobi' externalptr getAnywhere(str.default)$objs[[1]](x) Classes 'ggobiDataset', 'data.frame' externalptr - attr(*, ggobi)=Class 'ggobi' externalptr One idea is to check what the class is of the output of your .GGobiCall. If it were of class ggobiDataset then it would in turn be calling as.data.frame.ggobiDataset again. I mention this since the problem in your reproducible example was precisely of that sort. class(x) [1] ggobiDataset data.frame class(as.data.frame(x)) [1] data.frame So I don't think that's the problem. Thanks again for your help. 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
[R] Suppress legend in plotting groupedData
Dear All, I would be grateful if you can tell how can I suppress the legend (automatically created) when I plot a groupedData. suppose that I have data farme df_0 which contains the following columns: ID represents the clusters, t represents time of observation, Y: the variable of interest and X a binary covariate. Then: df_1 - groupedData(Y ~ t | ID) plot(df_1, outer ~X) This gives me the plot (Y ~ t) for each cluster with a legend containing all the clusters. How can I drop the legend? Thanks you very much, Bernard, - [[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
[R] tcltk loading in R-2.2.1 from src
Hi, Having trouble loading tcltk in R 2.2.1 built from source. ./configure, make, make check, and make install run ok. library(tcltk) Error in firstlib(which.lib.loc, package) : Tcl/Tk support is not available on this system Error in library(tcltk) : .First.lib failed for 'tcltk' even though it is listed in library() output. I have the same problem even if i compile with options: ./configure --with-tcltk --with-tcl-config=/usr/lib/tclConfig.sh --with-tk-config=/usr/lib/tkConfig.sh Is there a dep for R 2.2.1 on a specific version of tcl? Any hints on this issue appreciated. running on linux (2.4.21-4) version: rpm -qa | grep tcl tcl-devel-8.3.5-92 tcl-8.3.5-92 Specifically, the package pbatR loads this library during installation. Thanks, Patrice -- Patrice Seyed Linux System Administrator - LinGA RHCE, SCSA Boston University Medical Campus __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] RES: survival
On Wed, 8 Mar 2006, Paulo Brando wrote: summary(model.fit) # just one species from one treatment shown below Call: survfit(formula = Surv(time, censo) ~ treatment + species, data = wsuv) treatment=0, species=1 time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 15440 3860.975 0.001260.9730.977 2 15054 3360.953 0.001700.9500.957 3 14668 3020.934 0.002000.9300.938 4 14282 2960.914 0.002260.9100.919 5 13896 2810.896 0.002470.8910.901 6 13510 2640.878 0.002640.8730.883 7 13124 2510.861 0.002800.8560.867 8 12738 2320.846 0.002930.8400.852 9 12352 2160.831 0.003050.8250.837 10 11966 2060.817 0.003150.8110.823 11 11580 1900.803 0.003250.7970.810 12 11194 1790.790 0.003330.7840.797 13 10808 1670.778 0.003410.7720.785 14 10422 1670.766 0.003490.7590.773 15 10036 1450.755 0.003560.7480.762 16 9650 1420.744 0.003630.7370.751 17 9264 1350.733 0.003690.7260.740 18 8878 1220.723 0.003750.7150.730 19 8492 990.714 0.003800.7070.722 20 8106 840.707 0.003850.6990.714 21 7720 680.701 0.003890.6930.708 22 7334 660.694 0.003930.6870.702 23 6948 510.689 0.003970.6810.697 24 6562 400.685 0.004000.6770.693 25 6176 380.681 0.004030.6730.689 26 5790 370.676 0.004070.6690.684 27 5404 330.672 0.004110.6640.680 28 5018 310.668 0.004150.6600.676 29 4632 260.664 0.004190.6560.673 30 4246 220.661 0.004230.6530.669 31 3860 150.658 0.004270.6500.667 32 3474 140.656 0.004310.6470.664 33 3088 140.653 0.004360.6440.661 34 2702 130.650 0.004430.6410.658 35 2316 120.646 0.004510.6380.655 36 1930 110.643 0.004620.6340.652 37 1544 120.638 0.004800.6280.647 38 1158 100.632 0.005070.6220.642 39772 90.625 0.005570.6140.636 40386 80.612 0.007090.5980.626 I don't get why with 8 leaves remaining (out of 384), the survival is about 0.6??? It looks as though the majority of your leaves are censored, especially at later time points. At each of your 40 time points about 1-2% of the leaves under observation die, so the survival curve should end up somewhere between 0.98^40 =0.45 and 0.99^40=0.67, and it does. -thomas __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] tcltk loading in R-2.2.1 from src
Patrice - I had a very similar problem using TCL/TK 8.3. Below is the email I sent to my computing group at work about how I fixed it. Note that since my TCL/TK header (.h) files were in an odd location, the first step probably isn't relevant for you. But I bet the second step is. -- I believe I have found the solution to this problem. There were 2 steps I took to get a build of R in my home directory that properly uses Tcl/Tk 8.3 under Linux. The first was setting an environment variable to let R know where the tcl.h and tk.h files reside. This environment variable was TCLTK_CPPFLAGS and is set to -I/s/include . This can bet set in the config.site in R's build directory also. (NOTE to R-help: The above is site specific to my location, /s is for software on the network.) Second, there is some problem with the way R interacts with the tkConfig.sh file in /s/lib (or wherever your .sh file is located). It comes from the following line in tkConfig.sh. TK_XINCLUDES='# no special path needed' The # character, indicating a comment, is somehow misinterpreted by the configure script which breaks the Tcl/Tk functionality. I was able to get around this by simply removing the comment and leaving it as TK_XINCLUDES='' That got me a working version of the newest R using Tcl/Tk. I'm not sure if Tcl/Tk version 8.4 would still put that comment in there, I believe it's changed though, so you'd only have to follow the first step to get R compiled with Tcl/Tk support. I found reference to this problem on the R mailing list, it appears to only affect certain installations of Tcl/Tk. - HTH, Erik Iverson Patrice Seyed wrote: Hi, Having trouble loading tcltk in R 2.2.1 built from source. ./configure, make, make check, and make install run ok. library(tcltk) Error in firstlib(which.lib.loc, package) : Tcl/Tk support is not available on this system Error in library(tcltk) : .First.lib failed for 'tcltk' even though it is listed in library() output. I have the same problem even if i compile with options: ./configure --with-tcltk --with-tcl-config=/usr/lib/tclConfig.sh --with-tk-config=/usr/lib/tkConfig.sh Is there a dep for R 2.2.1 on a specific version of tcl? Any hints on this issue appreciated. running on linux (2.4.21-4) version: rpm -qa | grep tcl tcl-devel-8.3.5-92 tcl-8.3.5-92 Specifically, the package pbatR loads this library during installation. Thanks, Patrice __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Need help on Coxph
Hi all, I have a dataset which includes 84 rows and 4313 columns. Starting from the 2nd row, each row represents a patient. The 1st column is for arrayID 2nd column is for time 3rd column is for cancer 4th column is for patientID Starting for the 5th columns, each column's ID is like CTC_232B23, RP11_181G12, RP11_62M23... I extract first 5 rows and 7 columns and get the data in below: arrayID time cancer patientID CTC_232B23 RP11_181G12 RP11_62M23 X1747_4224 37.633 0 30635 -0.02665-0.02665 -0.038025 X1750_4214 89.300 0 22158 -0.02665-0.02665 -0.038025 X1751_4208 53.333 0 31669 -0.02665-0.02665 -0.038025 X1754_4194 39.467 0 32849 -0.02665-0.02665 -0.038025 X1775_4497 84.900 0 33563 -0.02665-0.02665 -0.038025 Finally I would like to build the cox model for each column (starting from the 5th column). The code would be like: mod.allison - coxph( Surv(time, cancer) ~ CTC_232B23 , data=Rossi2) mod.allison - coxph( Surv(time, cancer) ~ RP11_181G12 , data=Rossi2) mod.allison - coxph( Surv(time, cancer) ~ RP11_62M23 , data=Rossi2) I have no problem getting the result. However, since there are 4309 columns (counted begining at the 5th column), I decide to make a 4309*1 matrix (called BACs) to store the name of each column. Then run a loop like: for (i in 1:4309){ mod.allison - coxph( Surv(time, cancer) ~ BAC[i,1], data=Rossi2) } mod.allison However, I got the error. Would you like to give me some suggestions about this? Many thanks in advance! [[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
[R] R2.2.1-patched build failed with PGI 6.1 on x86-64
Hi All, While attempting to build R-patched using the Portland Group compiler suite on our dual Opteron 250 box running Scyld 29cz5 (based on RH, kernel 2.4.29), I get: make[5]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library/tools/src' make[5]: Entering directory `/home/andy/Rbuild/R221-pgi/src/library/tools/src' pgcc -I../../../../include -I/include -I/include/CC -fpic -O2 -Kieee -c /home/andy/Rbuild/R-patched/src/library/tools/src/text.c -o text.o pgcc -I../../../../include -I/include -I/include/CC -fpic -O2 -Kieee -c /home/andy/Rbuild/R-patched/src/library/tools/src/init.c -o init.o pgcc -I../../../../include -I/include -I/include/CC -fpic -O2 -Kieee -c /home/andy/Rbuild/R-patched/src/library/tools/src/Rmd5.c -o Rmd5.o pgcc -I../../../../include -I/include -I/include/CC -fpic -O2 -Kieee -c /home/andy/Rbuild/R-patched/src/library/tools/src/md5.c -o md5.o pgcc -shared -L/libso -L/usr/lib64 -o tools.so text.o init.o Rmd5.o md5.o mkdir -p -- ../../../../library/tools/libs make[5]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library/tools/src' make[4]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library/tools/src' **ERROR: in routine alloca() there is a stack overflow: thread 0, max 8180KB, used 0KB, request 16B make[3]: *** [all] Error 1 make[3]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library/tools' make[2]: *** [R] Error 1 make[2]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library' make[1]: *** [R] Error 1 make[1]: Leaving directory `/home/andy/Rbuild/R221-pgi/src' make: *** [R] Error 1 I configured it as follows: ~/Rbuild/R-patched/configure \ PG_HOME=/usr/pgi/linux86-64/6.1 \ CC=pgcc \ CFLAGS=-O2 -Kieee \ CPPFLAGS=-I$PG_HOME/include -I$PG_HOME/include/CC \ CPICFLAGS=-fpic \ F77=pgf77 \ FFLAGS=-Kieee\ FPICFLAGS=-fpic \ CXX=pgCC \ CXXFLAGS=-Kieee\ CXXPICFLAGS=-fpic \ SHLIB_CXXLDFLAGS=-shared \ SHLIB_LDFLAGS=-shared \ LDFLAGS=-L$PG_HOME/libso -L/usr/lib64 \ --with-blas=-L/usr/pgi/linux86-64/6.1/libso -lacml (Note that pgcc no longer recognize -mieee-fp, so I changed to -Kieee instead, and had to manually remove -mieee-fp in the resulting Makeconf and config.status.) I'd very much appreciate any pointer as to how to proceed... Best, Andy Andy Liaw, PhD Biometrics ResearchPO Box 2000 RY33-300 Merck Research LabsRahway, NJ 07065 andy_liaw(a)merck.com 732-594-0820 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Suppress legend in plotting groupedData
On 3/9/06, Marc Bernard [EMAIL PROTECTED] wrote: Dear All, I would be grateful if you can tell how can I suppress the legend (automatically created) when I plot a groupedData. suppose that I have data farme df_0 which contains the following columns: ID represents the clusters, t represents time of observation, Y: the variable of interest and X a binary covariate. Then: df_1 - groupedData(Y ~ t | ID) plot(df_1, outer ~X) This gives me the plot (Y ~ t) for each cluster with a legend containing all the clusters. How can I drop the legend? key=FALSE. See ?plot.nfnGroupedData or ?plot.nffGroupedData (whichever is appropriate according to the class(df_1)) for details. 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
Re: [R] R2.2.1-patched build failed with PGI 6.1 on x86-64
For some reason I had to see my own post to realize my stupidity... Defining PG_HOME _before_ running the configure solved the problem. Sorry for wasting the bandwidth. Andy From: Liaw, Andy Hi All, While attempting to build R-patched using the Portland Group compiler suite on our dual Opteron 250 box running Scyld 29cz5 (based on RH, kernel 2.4.29), I get: make[5]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library/tools/src' make[5]: Entering directory `/home/andy/Rbuild/R221-pgi/src/library/tools/src' pgcc -I../../../../include -I/include -I/include/CC -fpic -O2 -Kieee -c /home/andy/Rbuild/R-patched/src/library/tools/src/text.c -o text.o pgcc -I../../../../include -I/include -I/include/CC -fpic -O2 -Kieee -c /home/andy/Rbuild/R-patched/src/library/tools/src/init.c -o init.o pgcc -I../../../../include -I/include -I/include/CC -fpic -O2 -Kieee -c /home/andy/Rbuild/R-patched/src/library/tools/src/Rmd5.c -o Rmd5.o pgcc -I../../../../include -I/include -I/include/CC -fpic -O2 -Kieee -c /home/andy/Rbuild/R-patched/src/library/tools/src/md5.c -o md5.o pgcc -shared -L/libso -L/usr/lib64 -o tools.so text.o init.o Rmd5.o md5.o mkdir -p -- ../../../../library/tools/libs make[5]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library/tools/src' make[4]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library/tools/src' **ERROR: in routine alloca() there is a stack overflow: thread 0, max 8180KB, used 0KB, request 16B make[3]: *** [all] Error 1 make[3]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library/tools' make[2]: *** [R] Error 1 make[2]: Leaving directory `/home/andy/Rbuild/R221-pgi/src/library' make[1]: *** [R] Error 1 make[1]: Leaving directory `/home/andy/Rbuild/R221-pgi/src' make: *** [R] Error 1 I configured it as follows: ~/Rbuild/R-patched/configure \ PG_HOME=/usr/pgi/linux86-64/6.1 \ CC=pgcc \ CFLAGS=-O2 -Kieee \ CPPFLAGS=-I$PG_HOME/include -I$PG_HOME/include/CC \ CPICFLAGS=-fpic \ F77=pgf77 \ FFLAGS=-Kieee\ FPICFLAGS=-fpic \ CXX=pgCC \ CXXFLAGS=-Kieee\ CXXPICFLAGS=-fpic \ SHLIB_CXXLDFLAGS=-shared \ SHLIB_LDFLAGS=-shared \ LDFLAGS=-L$PG_HOME/libso -L/usr/lib64 \ --with-blas=-L/usr/pgi/linux86-64/6.1/libso -lacml (Note that pgcc no longer recognize -mieee-fp, so I changed to -Kieee instead, and had to manually remove -mieee-fp in the resulting Makeconf and config.status.) I'd very much appreciate any pointer as to how to proceed... Best, Andy Andy Liaw, PhD Biometrics ResearchPO Box 2000 RY33-300 Merck Research LabsRahway, NJ 07065 andy_liaw(a)merck.com 732-594-0820 __ 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 -- Notice: This e-mail message, together with any attachments, contains information of Merck Co., Inc. (One Merck Drive, Whitehouse Station, New Jersey, USA 08889), and/or its affiliates (which may be known outside the United States as Merck Frosst, Merck Sharp Dohme or MSD and in Japan, as Banyu) that may be confidential, proprietary copyrighted and/or legally privileged. It is intended solely for the use of the individual or entity named on this message. If you are not the intended recipient, and have received this message in error, please notify us immediately by reply e-mail and then delete it from your system. -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Multiple logistic regression
In addition to the multinom(nnet) function mentioned below there is some literature on how one can divide such polytomous problems into an set of dichotomous classifications and then aggregate the results, e.g.: 1) one-vs-all 2) pairwise comparisons (aka [double] round-robin) (Führnkranz) 3) nested dichotomies 3) ensembles of nested dichotomies (aka ENDs) (Frank Kramer) The article by Eibe Frank Stefan Kramer, Ensembles of nested dichotomies for multi-class problems http://wwwkramer.in.tum.de/kramer/frankkramer_icml04.pdf firstly gives an concise overview of the various above strategies and compares their performance, arguing for the use of the method they have themselves devised, i.e. ENDs, and secondly provides references for articles describing the other methods in detail (e.g. Führnkranz). The strategies mentioned above have the advantage that they do not have a default class, in contrast to the multinom function. Another question is whether any of these strategies have been implemented in a publicly avaiblable library? At least my recent cursory searches in the R-help archives and with help.search(...) have not produced any tangible results. I've managed to concoct a set of R-functions which crudely implement the strategies 1) one-vs-all and 2) pairwise comparisons, which I attach below. They are probably too much geared to my own research question and cut a few too many corners to be used more generally without substantial modification, and they could most probably be implemented in a more elegant manner, but they might nevertheless be of some inspiration. Having hacked these solutions on my own it would be all too typical that some of the above multilevel classification strategies have in fact already been in implemented in an available library. So, is anyone on this list aware of such functions/libraries? Regards, -Antti Arppe == Antti Arppe - Master of Science (Engineering) Researcher doctoral student (Linguistics) E-mail: [EMAIL PROTECTED] WWW: http://www.ling.helsinki.fi/~aarppe 13. Multiple logistic regression (Stephanie Delalieux) Date: Wed, 8 Mar 2006 14:15:58 +0100 From: Stephanie Delalieux [EMAIL PROTECTED] Subject: [R] Multiple logistic regression To: r-help@stat.math.ethz.ch Is there a function in R that classifies data in more than 2 groups using logistic regression/classification? I want to compare the c-indices of earlier research (lrm, binary response variables) with new c-indices obtained from 'multiple' (more response variables) logistic regression. Message: 23 Date: Wed, 8 Mar 2006 22:26:24 +0800 From: ronggui [EMAIL PROTECTED] Subject: Re: [R] Multiple logistic regression To: Stephanie Delalieux [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Do you mean multinomial logistic model? If it is,the multinom function in nnet package and multinomial function in VGAM(http://www.stat.auckland.ac.nz/~yee) package can do it. 8- 1) dat: data (with the first column containing the multiclass variable which is being predicted) 2) fn: predictor variables as a string, e.g. fn - A + B + C. In this implementation, the predictor variables are assumed to be logical (and thus binary); therefore, the GLM model family=binomial, and should be changed if the data is of another sort. 3) lex: list with multiple classes being predicted, e.g. lex - c(a, b, c, d) 4) freq: a Nx1 vector mapping frequency order of predicted classes to their actual order in (3) lex, needed for the double-round method for determining ties (- alternative with the highest frequency selected) 5) teach.test.ratio: a list of length(2) indicating the proportions of the data to be used for teaching the models and testing, e.g. c(1,1) - 50% teach vs. 50% testing, c(2,1) - 66.6% vs. 33.3% 6) iter: number of iteration rounds in evaluating the accuracy of classication performance 7) classifier: either 'double.round.robin' or 'one.vs.all' repeated.tests - function(dat,fn,lex,freq,teach.test.ratio=c(1,1),iter=1,hold.out=FALSE,classifier=double.round.robin, ...) { n.tot = nrow(dat); if(length(teach.test.ratio)==2) n.teach=round(teach.test.ratio[1]*n.tot/sum(teach.test.ratio)); n.test = n.tot - n.teach; nlex - length(lex); success - matrix(c(n.teach,round(n.teach*100/n.tot,2),n.test,round(n.test*100/n.tot,2),0,0,0),iter,7,byrow=TRUE); colnames(success) - c(Teach,%,Test,%,Success,%,tau (Kendall)); test.lx - matrix(0,iter,nlex); colnames(test.lx) - lex; success.lx - guess.lx - test.lx; for(i in 1:iter) { selected - sample(seq(1:n.tot),n.teach,replace=hold.out); teach - dat[selected,]; test - dat[-selected,]; result - switch(classifier, double.round.robin = double.round.robin(teach,test,fn,lex,freq), one.vs.all = one.vs.all(teach,test,fn,lex)); for(j in 1:n.test) {
Re: [R] tcltk loading in R-2.2.1 from src
Patrice Seyed [EMAIL PROTECTED] writes: Hi, Having trouble loading tcltk in R 2.2.1 built from source. ./configure, make, make check, and make install run ok. library(tcltk) Error in firstlib(which.lib.loc, package) : Tcl/Tk support is not available on this system Error in library(tcltk) : .First.lib failed for 'tcltk' even though it is listed in library() output. Well, the package is there, as a stub to tell you that it doesn't work I have the same problem even if i compile with options: ./configure --with-tcltk --with-tcl-config=/usr/lib/tclConfig.sh --with-tk-config=/usr/lib/tkConfig.sh In both cases, do check the configure output for clues. Is there a dep for R 2.2.1 on a specific version of tcl? Any hints on this issue appreciated. Anything later than 8.1 or so should work. Some old versions had bugs in tclConfig.sh/tkConfig.sh which didn't quite give out enough info to tell configure where to find the include files and libraries. I'm not quite sure when that was, but we're up to 8.4.9 now and tcl/tk has a notoriously slow release cycle, so 8.3.5 sound quite ancient running on linux (2.4.21-4) version: Of which distribution? rpm -qa | grep tcl tcl-devel-8.3.5-92 tcl-8.3.5-92 and rpm -qa | grep tk says what? -- 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
Re: [R] install.packages saying the car package is not in therepositories
Thank you, this did the trick. In which case $ apt-get install r-cran-car is your friend. Is there anyway to upgrade my R install through CRAN? Or is there some Debian repository that has an upgraded version of R? Thanks again. Jeremy __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] shift / rota
How to do a shift/rotate os a list? if a = c(1,2,3) what is the best way to make a equal c(3,1,2) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] install.packages saying the car package is not in therepositories
On 9 March 2006 at 20:47, Jeremy Morris wrote: | Thank you, this did the trick. | | In which case | | $ apt-get install r-cran-car | | is your friend. | | Is there anyway to upgrade my R install through CRAN? Or is there | some Debian repository that has an upgraded version of R? Yes. Please re-read the entire message I sent you yesterday as it already answered this question. (Hint: The answer is in the R FAQ too.) Hope this helps, Dirk -- Hell, there are no rules here - we're trying to accomplish something. -- Thomas A. Edison __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] shift / rota
On Thu, 9 Mar 2006, Omar Lakkis wrote: How to do a shift/rotate os a list? if a = c(1,2,3) what is the best way to make a equal c(3,1,2) a - c(a[length(a)],a[-length(a)]) or n - length(a) a - c(a[n],a[-n]) David Scott _ David Scott Department of Statistics, Tamaki Campus The University of Auckland, PB 92019 AucklandNEW ZEALAND Phone: +64 9 373 7599 ext 86830 Fax: +64 9 373 7000 Email: [EMAIL PROTECTED] Graduate Officer, Department of Statistics __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] shift / rota
a = c(1,2,3) a [1] 1 2 3 rev(a) [1] 3 2 1 PS: a in your example is not a list; i.e class(a) From: Omar Lakkis [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Subject: [R] shift / rota Date: Thu, 9 Mar 2006 15:51:51 -0500 How to do a shift/rotate os a list? if a = c(1,2,3) what is the best way to make a equal c(3,1,2) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Optimal platform for R
I'm looking to buy a new desktop which will primarily be used for analyses of large datasets (100s of MB). I've seen postings from several years back re the 'optimal' platform for running R, but nothing more recently. Specifically, I want to know: 1) if I run R under Windows, does having a dual-processor machine help speed things up? And 2) is it still true that R performs about as well under Windows as Linux? Thanks, Greg Gregory Wellenius, ScD Cardiology Research Fellow Beth Israel Deaconess Medical Center 330 Brookline Avenue, Deaconess 306 Boston, MA 02215 617-632-7680 (phone) 617-632-7698 (fax) [[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
Re: [R] Optimal platform for R
On 3/9/2006 4:47 PM, [EMAIL PROTECTED] wrote: I'm looking to buy a new desktop which will primarily be used for analyses of large datasets (100s of MB). I've seen postings from several years back re the 'optimal' platform for running R, but nothing more recently. Specifically, I want to know: 1) if I run R under Windows, does having a dual-processor machine help speed things up? And 2) is it still true that R performs about as well under Windows as Linux? For a big dataset, you're better off with a 64 bit version of R. There isn't one for Windows yet, and won't be for quite a while, since the build tools (gcc, etc.) don't exist. So you're probably better off in Linux. Duncan Murdoch Thanks, Greg Gregory Wellenius, ScD Cardiology Research Fellow Beth Israel Deaconess Medical Center 330 Brookline Avenue, Deaconess 306 Boston, MA 02215 617-632-7680 (phone) 617-632-7698 (fax) [[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 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] newbie question: grouping rows
Hi all, I have a very simple question that I can't seem to find the answer to. How do I extract rows that meet a certain criteria from a data frame and group them into a new data frame? For example, if I want to make a new data frame that only includes rows of data for which the p values (given by one of the columns in the data frame) are less than a certain value, how do I do this? It seems that there should be a simple function that does this. I looked into getGroups from the nmle package, but am not sure how to construct the form argument correctly or even if it's the appropriate way to tackle this. Thanks in advance of your answer, Matt __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Trellis - setting xlim or ylim by data range in whole column or row
Dear List-mates, I have been trying to set up a 4x8 trellis plot (that is, 4 columns, 8 rows), and would like to have the axis limits set based on the data range of rows (for ylim) and columns (for xlim). I've been using the call: foo-xyplot(y~x|Epoch+Subject, type=c(l,r), par.strip.text=list(cex=0.5), ...) and updating to see different effects of scale adjustments, etc. foo-update(foo, scale=list(relation=sliced)) #etc. I can have each panel with its own limits using relation=free, limits wide enough to accommodate the whole data range with same or the limits of the widest panel for all panels with split. I have not, however, figured out a way to have the limits for x y grouped by their column or row. Documentation points to 3 alternate ways of trying to set stuff in panels. (1) a prepanel function, (2) using scales, or (3) focusing in on each panel individually changing some setting. I've played around with accessing different panels individually with foo[column, row], and using a list to determine which get displayed (instead of using skip because I can't get skip to work). Would I be able to set the xlim values in a similar way? foo[1,]$ylim-newvalues to set a whole columns ylims (e.g by data range of y in conditioning variable 2 (subject in my case)) and foo[,1]$xlim-newvalues to get a whole rows xlims by the data range of x in each level of conditioning variable 1 (Epoch in my formula))? If so, what attribute should I access, and if not what would you recommend? I've been reading posts, working examples in Venables Ripley 4th Ed., and experimenting with different things for the last 4 days. I'm still not used to the lattice terminology, so I could have easily miss interpreted what something was meant for (example- prepanel makes no sense to me yet). Conversely, I got a lot farther, a lot faster, using Lattice than I did using plot or ts.plot. In addition to a much shorter list of attributes that don't make sense to me yet than otherwise, I have been really tickled with the Lattice package. Thanks in advance for your feedback. KeithC. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] newbie question: grouping rows
Matt, Have a look at subset specially the examples at the end. I use it a lot. Hope it helps, Augusto Augusto Sanabria. MSc, PhD. Mathematical Modeller Risk Research Group Geospatial Earth Monitoring Division Geoscience Australia (www.ga.gov.au) Cnr. Jerrabomberra Av. Hindmarsh Dr. Symonston ACT 2609 Ph. (02) 6249-9155 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Matthew Scholz Sent: Friday, 10 March 2006 9:18 AM To: r-help@stat.math.ethz.ch Subject: [R] newbie question: grouping rows Hi all, I have a very simple question that I can't seem to find the answer to. How do I extract rows that meet a certain criteria from a data frame and group them into a new data frame? For example, if I want to make a new data frame that only includes rows of data for which the p values (given by one of the columns in the data frame) are less than a certain value, how do I do this? It seems that there should be a simple function that does this. I looked into getGroups from the nmle package, but am not sure how to construct the form argument correctly or even if it's the appropriate way to tackle this. Thanks in advance of your answer, Matt __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] newbie question: grouping rows
You can try: new.dataframe - my.dataframe[my.dataframe$p.value 0.05, ] This will select all columns. Alternatively, you can specify the columns that you want after the ,. -Christos -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Matthew Scholz Sent: Thursday, March 09, 2006 5:18 PM To: r-help@stat.math.ethz.ch Subject: [R] newbie question: grouping rows Hi all, I have a very simple question that I can't seem to find the answer to. How do I extract rows that meet a certain criteria from a data frame and group them into a new data frame? For example, if I want to make a new data frame that only includes rows of data for which the p values (given by one of the columns in the data frame) are less than a certain value, how do I do this? It seems that there should be a simple function that does this. I looked into getGroups from the nmle package, but am not sure how to construct the form argument correctly or even if it's the appropriate way to tackle this. Thanks in advance of your answer, Matt __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] newbie question: grouping rows
?subset Marc -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Matthew Scholz Sent: Thursday, March 09, 2006 2:18 PM To: r-help@stat.math.ethz.ch Subject: [R] newbie question: grouping rows Hi all, I have a very simple question that I can't seem to find the answer to. How do I extract rows that meet a certain criteria from a data frame and group them into a new data frame? For example, if I want to make a new data frame that only includes rows of data for which the p values (given by one of the columns in the data frame) are less than a certain value, how do I do this? It seems that there should be a simple function that does this. I looked into getGroups from the nmle package, but am not sure how to construct the form argument correctly or even if it's the appropriate way to tackle this. Thanks in advance of your answer, Matt __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] import from LISREL output of parameter estimates
I am using R and LISREL for simulation studies. R generates the data that is analyzed with LISREL. In LISREL I use PV in the LISREL output statement to request estimated variances. LISREL writes these in a file that looks like this: 1 0 0 0.100331D+01 0.144845D+01 0.141009D+01 0.214423D+01 0.214129D+01 0.194464D+01 0.191531D+01 0.198328D+01 0.100683D+00-0.236392D-01 0.200655D+01 0.100142D+01 0.572501D+01 0.131652D+01 0.112175D+01 0.186140D+01 0.212321D+01 0.207455D+01 0.201915D+01-0.224004D+01 0.966613D+00 0.459977D+01 0.133921D+01 0.613532D+01 0.201852D+01 0.198881D+01 0.203360D+01 0.200934D+01-0.115235D+01 0.975292D+00 -0.756997D-01 0.341922D+01 0.334352D+01 0.360463D-01-0.112819D+00 I need a function that reads these numbers back into R. I tried around with read.fwf(pv.lpr, widths=rep(13, 6), skip=1) read.fortran() scan() without success. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] import from LISREL output of parameter estimates
The ouput from LISREL is a matrix,so the read.matrix(tseries) will do the job in this situation. ?read.matrix read.matrix(tseries) R Documentation Read Matrix Data Description Reads a matrix data file. Usage read.matrix(file, header = FALSE, sep = , skip = 0) 2006/3/10, Felix Flory [EMAIL PROTECTED]: I am using R and LISREL for simulation studies. R generates the data that is analyzed with LISREL. In LISREL I use PV in the LISREL output statement to request estimated variances. LISREL writes these in a file that looks like this: 1 0 0 0.100331D+01 0.144845D+01 0.141009D+01 0.214423D+01 0.214129D+01 0.194464D+01 0.191531D+01 0.198328D+01 0.100683D+00-0.236392D-01 0.200655D+01 0.100142D+01 0.572501D+01 0.131652D+01 0.112175D+01 0.186140D+01 0.212321D+01 0.207455D+01 0.201915D+01-0.224004D+01 0.966613D+00 0.459977D+01 0.133921D+01 0.613532D+01 0.201852D+01 0.198881D+01 0.203360D+01 0.200934D+01-0.115235D+01 0.975292D+00 -0.756997D-01 0.341922D+01 0.334352D+01 0.360463D-01-0.112819D+00 I need a function that reads these numbers back into R. I tried around with read.fwf(pv.lpr, widths=rep(13, 6), skip=1) read.fortran() scan() without success. __ 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 -- 黄荣贵 Deparment of Sociology Fudan University __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] 2nd R console?
hello all, i'm forwarding this question for a colleague. Is it possible to open a 2nd R Console? regards, mark+ -- mark garey ucsf department of epidemiology and biostatistics division of biostatistics 185 berry street, suite 5700 san francisco, ca. 94107-1739 415.514.8147 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] 2nd R console?
On Thu, 9 Mar 2006, mark garey wrote: hello all, i'm forwarding this question for a colleague. Is it possible to open a 2nd R Console? The answer is probably No, but since it isn't clear what your colleague means by a 2nd R Console or what OS this is, it's hard to be sure. -thomas Thomas Lumley Assoc. Professor, Biostatistics [EMAIL PROTECTED] University of Washington, Seattle __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] coloring leaves in a hclust or dendrogram plot
Greetings, I have perused the r-help mailing list archives for an answer to this question, without avail. I would like to color the leaves of a dendrogram plot based on a cutoff in one of the variables involved in the initial clustering. My input data is in the form of: B K Alameda 0.2475770 0.7524230 Alpine0.4546784 0.5453216 Amador0.6278610 0.3721390 essentially rows labeled by county name, with two variables: percent voted for B and percent voted for K. While it is obvious that this is somewhat of a contrived example, I intend to use this as a learning device. Here is the code used to create and plot the dendrogram: hc - hclust(dist(y), ave) dend - as.dendrogram(hc) plot(dend, main=CA 2004 Election Results by County) An example of the output can be found here: http://casoilresource.lawr.ucdavis.edu/drupal/node/206?size=_original I have experimented with the edgePar and nodePar parameters for the plot.dendrogram() method, but have not been able to make sense of the output. The basis for setting the colors of the leaves in the dendrogram is a simple majority calculation: reds - y[y$B 0.5, ] blues - y[y$K 0.5, ] Such that leaves in the tree will be colored based on the membership in either of the two above groups. Is there a resource documenting how this might be accomplished? Any thoughts or ideas would be greatly appreciated. Cheers, Dylan -- Dylan Beaudette Soils and Biogeochemistry Graduate Group University of California at Davis 530.754.7341 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] maximum likelihood estimate
I suggest you go to www.r-project.org - CRAN - (select a local mirror) - Packages. Among the offerings there, you will find the following: distr Object orientated implementation of distributions distrEx Extensions of package distr distrSimSimulation classes based on package distr distrTEst Estimation and Testing classes based on package distr I don't know that any of these will solve your problem, but if you try them and still have problems, I suggest you then PLEASE do read the posting guide! www.R-project.org/posting-guide.html, and submit another question based on what you've learned in the interim. Hope this helps. spencer graves [EMAIL PROTECTED] wrote: Hi! Recently I try to find the method maximum likelihood for gamma,weibull,Pearson type III,Kappa Distribution, mixed exponential distribution, skew distribution. I have tried function ms() for gamma two parameters and weibull two parameters.It works but not for Pearson type III. I have problem to find the likelihood function for mixed exponential distribution and kappa distribution. So can anyone tell me the programming for those functions? especially in finding their maximum likelihood,because I have trouble to write down the likelihood function for mixed exponential distribution and the other two.If possible ,how to get the hessian matrix? Currently I am using SPlus 6.2. Thanks. suhaila __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Date and Times a la Dalgaard
Does anyone know of a resource for learning the basics of how to manage and manipulate dates and times in R? I have been reading Introductory Statistics with R by Peter Dalgaard which is fantastic. But alas, I could find no reference to date and time. I have looked at the reference manual but it is particularly unapproachable. So rather than dense technical talk I would rather see a dataset, several arguments and commentary of what is happening. The alternative to is to make Excel or Microsoft Access handle all the date and time work (I have had about 16 years experience with how to calculate and handle date and time in those programs) and then use RODBC to read it. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] interrupted time series analysis using ARIMA models
I'm familiar with Box and Tiao (1975) intervention analysis; I studied time series under Box and Tiao. I don't know how to do that in R, but there must be a way. Have you looked at the 'dse' bundle? That comes with vignettes that make it relatively easy to learn (or at least to learn the capabilities covered in the vignettes). The models you want may not be identified by the names with which you are familiar, but I believe they are probably available. If you try that and still have questions, I suggest you consult the posting guide (www.R-project.org/posting-guide.html) for help in crafting another question that may attract quicker and more useful replies. I also highly recommend the zoo package. It won't help you solve the problem you mentioned, but it might help you keep time stamps with your data. It, too, has a vignette to help people learn the capabilities. hope this helps. spencer graves Berta wrote: Dear R-users, Thanks Spencer for your suggestion, i think we are near but still that is not what i am looking for. I think I was not clear using that notation for the impact: (yt= d * yt-1 + w * It ), this yt is not my original series, it is only the impact, the series would be modeled as Yt=yt +Nt, with yt the impact written above and Nt the ARIMA part of the model. Hence, Yt is the series (your lh), and yt the impact. With your suggestion IntReg - cbind(It=(1:48)20, It.w=((1:48)20)*(1:48), It.lh=((1:48)20)*c(0, lh[-48]) ) arima(lh, order = c(1,0,0), xreg=IntReg) I would have for the original series Yt=lh(t) lh(20)=0 + Nt. lh(21)=w + beta1*21 + beta2*lh(20) + Nt lh(22)=w + beta1*22 + beta2*lh(21) + Nt etc. What I am trying to model is a gradual permanent impact, which would lead to: lh(t)= impact(t) + Nt lh(t)= w*It + d*yt-1 + Nt lh(20)= 0+ Nt lh(21)= w + Nt lh(22)= d*w + w + Nt lh(22)= (d^2)*w + d*w + w + Nt ... lh(n)=(d^n)*w +(d^(n-1))*w ++(d^2)*w + d*w + w + Nt, which asymptoticaly would be = w/(1-d) + Nt. In that way, I can model the impact not only as an abrupt permanent impact (like a step) but also as a gradual permanent impact (which grows gradually, as a linear trend or as a parabolic grow trend, or whatever) with just two parameters. In SAS they are called denominator factors for transfer functions for an input series. I also would like to modelize an abrupt temporary impact (a high pick in the moment of the impact decreasing gradually after it), but hopefully that will be easy after knowing the first. Any suggestion for implementing this would be very very well received!! Thank a lot in advance, Berta. answer of Spencer . Does the following illustrate the kind of interevention model you want IntReg - cbind(It=(1:48)20, It.w=((1:48)20)*(1:48), It.lh=((1:48)20)*c(0, lh[-48]) ) arima(lh, order = c(1,0,0), xreg=IntReg) hope this helps. spencer graves Berta wrote: Hi R-users, I am using arima to fit a time series. Now I would like to include an intervention component It (0 before intervention, 1 after) using different types of impacts, that is, not only trying the simple abrupt permanent impact (yt = w It ) with the xreg option but also trying with a gradual permanent impact (yt= d * yt-1 + w * It ), following the filosophy of Box and Tiao (1975). Intervention analysis with applications to economic and environmental problems. JASA 70: 70-92. Does anybody know where could I find how to incorporate them using the arima comand (or other), or a statistical package which can incorporate it? Thanks, Berta. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Date and Times a la Dalgaard
To manipulate date/time I use packages zoo and survival. Hope it helps, Augusto Augusto Sanabria. MSc, PhD. Mathematical Modeller Risk Research Group Geospatial Earth Monitoring Division Geoscience Australia (www.ga.gov.au) Cnr. Jerrabomberra Av. Hindmarsh Dr. Symonston ACT 2609 Ph. (02) 6249-9155 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Farrel Buchinsky Sent: Friday, 10 March 2006 2:39 PM To: r-help@stat.math.ethz.ch Subject: [R] Date and Times a la Dalgaard Does anyone know of a resource for learning the basics of how to manage and manipulate dates and times in R? I have been reading Introductory Statistics with R by Peter Dalgaard which is fantastic. But alas, I could find no reference to date and time. I have looked at the reference manual but it is particularly unapproachable. So rather than dense technical talk I would rather see a dataset, several arguments and commentary of what is happening. The alternative to is to make Excel or Microsoft Access handle all the date and time work (I have had about 16 years experience with how to calculate and handle date and time in those programs) and then use RODBC to read it. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Date and Times a la Dalgaard
On Thu, 2006-03-09 at 22:38 -0500, Farrel Buchinsky wrote: Does anyone know of a resource for learning the basics of how to manage and manipulate dates and times in R? I have been reading Introductory Statistics with R by Peter Dalgaard which is fantastic. But alas, I could find no reference to date and time. I have looked at the reference manual but it is particularly unapproachable. So rather than dense technical talk I would rather see a dataset, several arguments and commentary of what is happening. The alternative to is to make Excel or Microsoft Access handle all the date and time work (I have had about 16 years experience with how to calculate and handle date and time in those programs) and then use RODBC to read it. You might want to look at the following articles from R News: author = {Brian D. Ripley and Kurt Hornik}, title= {Date-Time Classes}, journal = {R News}, year = 2001, volume = 1, number = 2, pages= {8--11}, month= {June}, url = {http://CRAN.R-project.org/doc/Rnews/} author = {Gabor Grothendieck and Thomas Petzoldt}, title= {R Help Desk: Date and Time Classes in {R}}, journal = {R News}, year = 2004, volume = 4, number = 1, pages= {29--32}, month= {June}, url = {http://CRAN.R-project.org/doc/Rnews/} 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
Re: [R] 2nd R console?
On Windows, I often have more than one R console (instance of Rgui.exe) open at the same time. Unless by mistake, each instance is open on a unique working directory. MHP mark garey wrote on 3/9/2006 8:05 PM: hello all, i'm forwarding this question for a colleague. Is it possible to open a 2nd R Console? regards, mark+ -- mark garey ucsf department of epidemiology and biostatistics division of biostatistics 185 berry street, suite 5700 san francisco, ca. 94107-1739 415.514.8147 __ 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 -- Michael H. Prager, Ph.D. Population Dynamics Team NOAA Center for Coastal Habitat and Fisheries Research NMFS Southeast Fisheries Science Center Beaufort, North Carolina 28516 USA http://shrimp.ccfhrb.noaa.gov/~mprager/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Optimal platform for R
[EMAIL PROTECTED] wrote on 3/9/2006 4:47 PM: I'm looking to buy a new desktop which will primarily be used for analyses of large datasets (100s of MB). I've seen postings from several years back re the 'optimal' platform for running R, but nothing more recently. Specifically, I want to know: 1) if I run R under Windows, does having a dual-processor machine help speed things up? And 2) is it still true that R performs about as well under Windows as Linux? Duncan Murdoch has already answered your questions about operating systems. I would like to add that if there is significant I/O in what you are doing, fast SCSI disks are a worthwhile investment. The speed increase over ATA or SATA disks has been, in my experience, quite noticeable. -- Michael H. Prager, Ph.D. Population Dynamics Team NOAA Center for Coastal Habitat and Fisheries Research NMFS Southeast Fisheries Science Center Beaufort, North Carolina 28516 USA http://shrimp.ccfhrb.noaa.gov/~mprager/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Trellis - setting xlim or ylim by data range in whole column or row
On 3/9/06, Keith Chamberlain [EMAIL PROTECTED] wrote: Dear List-mates, I have been trying to set up a 4x8 trellis plot (that is, 4 columns, 8 rows), and would like to have the axis limits set based on the data range of rows (for ylim) and columns (for xlim). I've been using the call: foo-xyplot(y~x|Epoch+Subject, type=c(l,r), par.strip.text=list(cex=0.5), ...) and updating to see different effects of scale adjustments, etc. foo-update(foo, scale=list(relation=sliced)) #etc. I can have each panel with its own limits using relation=free, limits wide enough to accommodate the whole data range with same or the limits of the widest panel for all panels with split. I have not, however, figured out a way to have the limits for x y grouped by their column or row. Documentation points to 3 alternate ways of trying to set stuff in panels. (1) a prepanel function, (2) using scales, or (3) focusing in on each panel individually changing some setting. I've played around with accessing different panels individually with foo[column, row], and using a list to determine which get displayed (instead of using skip because I can't get skip to work). Would I be able to set the xlim values in a similar way? foo[1,]$ylim-newvalues to set a whole columns ylims (e.g by data range of y in conditioning variable 2 (subject in my case)) and foo[,1]$xlim-newvalues to get a whole rows xlims by the data range of x in each level of conditioning variable 1 (Epoch in my formula))? If so, what attribute should I access, and if not what would you recommend? What you want is not impossible, but not simple either. There's a good reason for that---this behavior makes sense only in a very special situation: when you have exactly two conditioning variables (a and b, say) and your layout has exactly nlevels(a) columns and nlevels(b) rows. To fix ideas, consider the following example: library(lattice) d - expand.grid(a = gl(4, 1, 20), b = gl(8, 1, 40)) d$x - with(d, rnorm(nrow(d), mean = as.numeric(a) + as.numeric(b))) d$y - with(d, rnorm(nrow(d), mean = as.numeric(a) + as.numeric(b))) foo is as in your example, foo - xyplot(y ~ x | a + b, d) At this point, foo has no layout defined (foo$layout is NULL). The layout will be determined when it is plotted. Since there are exactly two conditioning variables, the default here is layout = dim(foo) = c(4, 8). So plotting foo is the same as update(foo, layout = dim(foo)) However, if you instead do update(foo, layout = c(0, prod(dim(foo (which is essentially what happens when there is one conditioning variable) you will get a different layout, probably 6x6 (unless you have resized the device). In this case, choosing limits by row or column no longer makes sense. You might say, why not do this whatever the layout, whether it makes sense or not. The problem is, the axis limits for each panel needs to be determined _before_ the layout, because the layout may depend on the aspect ratio and the aspect ratio may depend on the axis limits (e.g. when aspect = xy). Workarounds: It is possible to explicitly specify a limit for each panel as a list. For example, you could have (update doesn't work well for this): xyplot(y ~ x | a + b, d, scales = list(x = list(relation = free, limits = rep(list(c(0, 11), c(1, 12), c(2, 13), c(4, 14)), 8))), par.strip.text = list(cex = 0.5)) If you are happy with this, that's great, but you will most likely want to omit all but one set of labels and use the freed up space. To control the labelling, you can specify the tick mark locations as a list just like the limits. An undocumented trick to simplify this is to have TRUE for the defaults, and NULL to omit them, so you could do: xyplot(y ~ x | a + b, d, scales = list(x = list(relation = free, limits = rep(list(c(0, 11), c(1, 12), c(2, 13), c(4, 14)), 8), at = rep(list(TRUE, NULL), c(4, 28, par.strip.text = list(cex = 0.5)) Finally, to make use of the space: xyplot(y ~ x | a + b, d, scales = list(x = list(relation = free, limits = rep(list(c(0, 11), c(1, 12), c(2, 13), c(4, 14)), 8), at = rep(list(TRUE, NULL), c(4, 28, par.settings = list(layout.heights = list(axis.panel = rep(c(1, 0), c(1, 7, par.strip.text = list(cex = 0.5)) I've been reading posts, working examples in Venables Ripley 4th Ed., and experimenting with different things for the last 4 days. I'm still not used to the lattice terminology, so I could have easily miss interpreted what something was meant for (example- prepanel makes no sense to me yet). It's a function that determines the limits for a panel given the data in it (these limits are then combined according to the value of
Re: [R] lsa and Rstem?
Hi Dave, On Mar 9, 2006, at 9:25 AM, Dave Atkins wrote: install.packages(Rstem, repos = http://www.omegahat.org/R;) This works. install.packages(lsa) Also works. But: library(lsa) Loading required package: Rstem Error in dyn.load(x, as.logical(local), as.logical(now)) : unable to load shared library '/Library/Frameworks/R.framework/ Resources/library/Rstem/libs/Rstem.so': dlopen(/Library/Frameworks/R.framework/Resources/library/Rstem/ libs/Rstem.so, 6): Library not loaded: libR.dylib Referenced from: /Library/Frameworks/R.framework/Resources/library/ Rstem/libs/Rstem.so Reason: image not found In addition: Warning message: package 'Rstem' was built under R version 2.3.0 Error: package 'Rstem' could not be loaded Running platform powerpc-apple-darwin7.9.0 arch powerpc os darwin7.9.0 system powerpc, darwin7.9.0 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R Does anyone have further advice? _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Trellis - setting xlim or ylim by data range in whole column or row
Dear Deepayan, My deepest thanks! Your example code worked perfectly. I understood the caveats you detailed about figuring out the axis limits before resolving the layout, which depend on a bunch of other stuff. I'm glad this is possible. Yes, this particular layout does only make sense in this special case (e.g. 2 conditioning variables exactly n and m-levels). Thank you again for your feedback. I am excited to go put your examples to work. Rgds, KeithC. -Original Message- From: Deepayan Sarkar [mailto:[EMAIL PROTECTED] Sent: Thursday, March 09, 2006 9:55 PM To: Keith Chamberlain Cc: r-help@stat.math.ethz.ch Subject: Re: Trellis - setting xlim or ylim by data range in whole column or row On 3/9/06, Keith Chamberlain [EMAIL PROTECTED] wrote: Dear List-mates, I have been trying to set up a 4x8 trellis plot (that is, 4 columns, 8 rows), and would like to have the axis limits set based on the data range of rows (for ylim) and columns (for xlim). I've been using the call: foo-xyplot(y~x|Epoch+Subject, type=c(l,r), par.strip.text=list(cex=0.5), ...) and updating to see different effects of scale adjustments, etc. foo-update(foo, scale=list(relation=sliced)) #etc. I can have each panel with its own limits using relation=free, limits wide enough to accommodate the whole data range with same or the limits of the widest panel for all panels with split. I have not, however, figured out a way to have the limits for x y grouped by their column or row. Documentation points to 3 alternate ways of trying to set stuff in panels. (1) a prepanel function, (2) using scales, or (3) focusing in on each panel individually changing some setting. I've played around with accessing different panels individually with foo[column, row], and using a list to determine which get displayed (instead of using skip because I can't get skip to work). Would I be able to set the xlim values in a similar way? foo[1,]$ylim-newvalues to set a whole columns ylims (e.g by data range of y in conditioning variable 2 (subject in my case)) and foo[,1]$xlim-newvalues to get a whole rows xlims by the data range of x in each level of conditioning variable 1 (Epoch in my formula))? If so, what attribute should I access, and if not what would you recommend? What you want is not impossible, but not simple either. There's a good reason for that---this behavior makes sense only in a very special situation: when you have exactly two conditioning variables (a and b, say) and your layout has exactly nlevels(a) columns and nlevels(b) rows. To fix ideas, consider the following example: library(lattice) d - expand.grid(a = gl(4, 1, 20), b = gl(8, 1, 40)) d$x - with(d, rnorm(nrow(d), mean = as.numeric(a) + as.numeric(b))) d$y - with(d, rnorm(nrow(d), mean = as.numeric(a) + as.numeric(b))) foo is as in your example, foo - xyplot(y ~ x | a + b, d) At this point, foo has no layout defined (foo$layout is NULL). The layout will be determined when it is plotted. Since there are exactly two conditioning variables, the default here is layout = dim(foo) = c(4, 8). So plotting foo is the same as update(foo, layout = dim(foo)) However, if you instead do update(foo, layout = c(0, prod(dim(foo (which is essentially what happens when there is one conditioning variable) you will get a different layout, probably 6x6 (unless you have resized the device). In this case, choosing limits by row or column no longer makes sense. You might say, why not do this whatever the layout, whether it makes sense or not. The problem is, the axis limits for each panel needs to be determined _before_ the layout, because the layout may depend on the aspect ratio and the aspect ratio may depend on the axis limits (e.g. when aspect = xy). Workarounds: It is possible to explicitly specify a limit for each panel as a list. For example, you could have (update doesn't work well for this): xyplot(y ~ x | a + b, d, scales = list(x = list(relation = free, limits = rep(list(c(0, 11), c(1, 12), c(2, 13), c(4, 14)), 8))), par.strip.text = list(cex = 0.5)) If you are happy with this, that's great, but you will most likely want to omit all but one set of labels and use the freed up space. To control the labelling, you can specify the tick mark locations as a list just like the limits. An undocumented trick to simplify this is to have TRUE for the defaults, and NULL to omit them, so you could do: xyplot(y ~ x | a + b, d, scales = list(x = list(relation = free, limits = rep(list(c(0, 11), c(1, 12), c(2, 13), c(4, 14)), 8), at = rep(list(TRUE, NULL), c(4, 28, par.strip.text = list(cex = 0.5)) Finally, to make use of the space: xyplot(y ~ x | a + b, d, scales = list(x = list(relation =
Re: [R] lsa and Rstem?
On Thu, 9 Mar 2006, Michael Kubovy wrote: Does anyone have further advice? Build Rstem from the sources on your own machine. install.packages(Rstem, repos = http://www.omegahat.org/R;, type = source) should work if you have your compilers etc set up to build packages from sources. Hi Dave, On Mar 9, 2006, at 9:25 AM, Dave Atkins wrote: install.packages(Rstem, repos = http://www.omegahat.org/R;) This works. install.packages(lsa) Also works. But: library(lsa) Loading required package: Rstem Error in dyn.load(x, as.logical(local), as.logical(now)) : unable to load shared library '/Library/Frameworks/R.framework/ Resources/library/Rstem/libs/Rstem.so': dlopen(/Library/Frameworks/R.framework/Resources/library/Rstem/ libs/Rstem.so, 6): Library not loaded: libR.dylib Referenced from: /Library/Frameworks/R.framework/Resources/library/ Rstem/libs/Rstem.so Reason: image not found In addition: Warning message: package 'Rstem' was built under R version 2.3.0 Error: package 'Rstem' could not be loaded Running platform powerpc-apple-darwin7.9.0 arch powerpc os darwin7.9.0 system powerpc, darwin7.9.0 status major2 minor2.1 year 2005 month12 day 20 svn rev 36812 language R Does anyone have further advice? _ Professor Michael Kubovy University of Virginia Department of Psychology USPS: P.O.Box 400400Charlottesville, VA 22904-4400 Parcels:Room 102Gilmer Hall McCormick RoadCharlottesville, VA 22903 Office:B011+1-434-982-4729 Lab:B019+1-434-982-4751 Fax:+1-434-982-4766 WWW:http://www.people.virginia.edu/~mk9y/ __ 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 -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] 2nd R console?
On Thu, 9 Mar 2006, Thomas Lumley wrote: On Thu, 9 Mar 2006, mark garey wrote: hello all, i'm forwarding this question for a colleague. Is it possible to open a 2nd R Console? The answer is probably No, but since it isn't clear what your colleague means by a 2nd R Console or what OS this is, it's hard to be sure. I would say the answer is almost certainly 'yes'. In all the cases I know of if you launch 2 (or more) R consoles you will get separate R processes with separate workspaces. The only problem I can envisage is that they might share a working directory and so the last one shut down could overwrite the workspace image and history saved from other consoles. (But that is only a problem if you save your work that way, and I rarely do.) -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] what's wrong with my cov?
Hi all, Why cov(y, y) only gives one value, and cov(t(y), t(y)) gives 3x3 NA matrix? Here my y is listed below and it is a 3x1 matrix. I am expecting that if I have a random vector y=[y1 y2 y3]', here ' denotes a transposition so that y is a column vector, where y1, y2, y3 are independent random variables... then cov(y, y) should be E[ y * y' ] - E[y] * E[y] ', and it should generate an 3x3 identity matrix. In fact, as you can see, my y below is generated by y=[f(1)+e1, f(2)+e2, f(3)+e3]', =[ 1^2+e1, 2^2+e2, 3^2+e3]', =[ 1 + e1, 4 + e1, 9 + e1]. where f(x)=x^2, e1, e2, e3 are three independent normal(0, 1) sample values... so I expect the cov(y, y) will be a 3x3 identity matrix. What's wrong here? cov(y, y) [,1] [1,] 13.78024 cov(t(y), t(y)) [,1] [,2] [,3] [1,] NA NA NA [2,] NA NA NA [3,] NA NA NA y [,1] [1,] 1.146766 [2,] 5.412608 [3,] 8.542067 [[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
Re: [R] Optimal platform for R
On Thu, 9 Mar 2006, [EMAIL PROTECTED] wrote: I'm looking to buy a new desktop which will primarily be used for analyses of large datasets (100s of MB). I've seen postings from several years back re the 'optimal' platform for running R, but nothing more recently. It is a subject which comes up every few months. Many of the developers are running dual (or dual-core) Opterons/Athlon 64s under Linux these days. Specifically, I want to know: 1) if I run R under Windows, does having a dual-processor machine help speed things up? And 2) is it still true that R performs about as well under Windows as Linux? Duncan Murdoch has already mentioned the 64-bit advantage if you need large datasets, but there is also a speed penalty if you do not. Your description seems on the margins (depends how many 100s and what the format is and what you want to do). One advantage of AMD64 Linux is that I can run either 32- or 64-bit versions of R and choose to have speed or space for any given task. A dual processor will be of little help in running R faster. R's interpreter is single-threaded, and although you can get some advantage in using multi-threaded BLAS libraries in large matrix computations these are not readily available for R under Windows, and the advantage is often small under Linux. Running two or more instances of R will take advantage of dual processers, and I have been running dual CPU machines for a decade. As for Windows vs Linux, R runs on the same hardware at about the same speed when comparing the standard Windows build with a shared library version on Linux (standard for e.g. the RH RPMs), but the standard Linux build is 10-20% faster. For one set of comparisons see http://sekhon.berkeley.edu/macosx/ -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] count pixels of same color in pixmap object?
Another solution is to use external tools. It really depends on what you have to do: if it is to count pixels of one, or a few gray levels, it could be fine to do it in R,... but if you want to count *all* pixels, this could be more efficient using a C program, especially if you work with 16bit images! You can find something in ImageMagick or in Netpbm. Here is a solution using pgmhist from Netpbm (http://netpbm.sourceforge.net/). This implies that you first converted your image in PGM format (the Netpbm library provides all tools required for that conversion). Here is the R code, assuming pgmhist or pgmhist.exe (Windows) is installed and accessible from the path: # Get statistics about pixel in graylevel pgm images pixelStats - function(imagePgm, pgmhist = if (.Platform$OS.type == windows) pgmhist.exe else pgmhist) { # Get pixel count for each gray level pix - system(paste(pgmhist, imagePgm), intern = TRUE, invisible = TRUE) if (pix[1] != value\tcount\tb%\tw%) stop(Error when running 'pgmhist' on the image , imagePgm, \n\n, pix) pix - pix[-(1:2)] getPixCount - function(str) as.numeric(strsplit(str, \t)[[1]][1:2]) pix2 - t(sapply(pix, getPixCount, USE.NAMES = FALSE)) colnames(pix2) - c(pixel.value, count) return(pix2) } pixelStats(myimage.pgm) Depending on the treatment you have to do, it is sometimes better to delegate it to specialized programs (ImageMagick, Netpbm, ImageJ, etc.) instead of using pixmap. The best is to try both and to determine which solution is faster for your particular application. I am working with very large images (hundreds of megabytes) for the ZooImage application (http://www.sciviews.org/zooimage), and I delegate 100% of the work done on these images to external programs, because R is not designed to handle them efficiently. Best, Philippe Grosjean Roger Bivand wrote: On Thu, 9 Mar 2006, Christian Jost wrote: Dear all, I try to figure out how to use R to count the number of pixels of the same color in some gray-level picture. I managed to read it in either tiff or jpeg format, but the returned pixmap object keeps its information out of (my) reach. Is there an easy way to tabulate the different color/graylevel pixels and their numbers? Or should I use a completely different (free) software? The objects are new-style class objects, so they are documented internally - query by getSlots() or from ?pixmap-class: library(pixmap) x - read.pnm(system.file(pictures/logo.ppm, package=pixmap)[1]) z - as(x, pixmapGrey) class(z) getSlots(class(z)) # lets you know what the slots are greydata - slot(z, grey) # extract the slot object summary(c(greydata)) # c() to flatten a matrix length(unique(c(greydata))) length((c(greydata))) # rather a lot of different values, let's round: rgreydata - round(greydata*10) length(unique(c(rgreydata))) # looks good table(c(rgreydata)) I assume that you converted from tiff or jpeg to pnm outside pixmap. Thanks for any hint, Christian. __ 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