[R] weights in glmrob
Dear List members, I am runnning a logistic regression that includes sample weights (expressed by a variable called WEIGHT) using glmob (from robustbase). My model looks like this: B12-glmrob(SERE~COMP+COM,family=binomial,weights=WEIGHT) When I include the weights, my standard error columns, as well as the p-values, show only NAs. If I exclude weights=WEIGHTS, everything works. I would like to use the weights. Does anyone know what could be the problem? Thanks in advance, Celso [[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] How to get R to ignore certain values when analyzing a column in a data table ?
On Sun, 9 Jul 2006, Daryl Manning wrote: Apologies if this is in (one of the many) manuals somewhere... Trying to switch to R from other stats programs. Basically, I have a large data table I've dumped from a DB, some of the values which are nulls '-' which I've converted to zeros. I've read it in using read.table fine. I want R to ignore the zero values when graphing or doing various other calculations. Is there a way to do this ? I did try to use NA but kept getting errors that x must be numeric. Did you use numeric NA? That is the way to so this, so you will need to give a reproducible example (as the posting guide asks). -- 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] KS Test Warning Message
Dear Justin Ties means that you have identical values in Year5.lm$residuals. Please remark that you can have a large R^2, but your residuals are not normally distributed. A large R^2 shows a strong linear relationship, but that does not say anything about the error distribution (see example below). So to answer your question. Yes it can take away validity of your model if the residuals are not normally distributed, especially tests and confidence intervals for your parameters are based on the normal assumption. I'd recommend to verify model assumptions by graphical tools, such as qqplot, Tukey-Anscombe Plot, ... Try: plot(Year5.lm) The power of KS-Test is quite small and graphical tools will give you a hint about your true error distribution instead of giving you only a p-value that tells you that the errors are not normal. set.seed(3) x - 1:100 ## t-distributed errors y - x + rt(100,2) ## Strong linear relationship plot(x,y) ## High R^2 due to strong linear relationship summary(reg - lm(y~x)) ## The residuals are not normal distributed qqnorm(resid(reg)) ## Small power of KS-Test. Violation of model assumption is not detected ks.test(resid(reg), pnorm) Best regards, Christoph Buser -- Christoph Buser [EMAIL PROTECTED] Seminar fuer Statistik, LEO C13 ETH Zurich 8092 Zurich SWITZERLAND phone: x-41-44-632-4673 fax: 632-1228 http://stat.ethz.ch/~buser/ -- justin rapp writes: All, Happy World Cup and Wimbledon. This morning finds me with the first of my many daily questions. I am running a ks.test on residuals obtained from a regression model. I use this code: ks.test(Year5.lm$residuals,pnorm) and obtain this output One-sample Kolmogorov-Smirnov test data: Year5.lm$residuals D = 0.7196, p-value 2.2e-16 alternative hypothesis: two.sided Warning message: cannot compute correct p-values with ties in: ks.test(Year5.lm$residuals, pnorm) I am wondering if anybody can tell me what this error message means. Also, could anybody clarify how I could have a regression model with a high Rsquared, rouglhy .67, but with nonnormal residuals? Does this take away from the validity of my model? jdr __ R-help@stat.math.ethz.ch 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] A possible too old question on significant test of correlation matrix
On Mon, 2006-07-10 at 13:27 +0800, Guo Wei-Wei wrote: Dear all, I'm working on a data.frame named en.data, which has n cases and m columns. I generate the correlation matrix of en.data by cor(en.data) I find that there is no p-value on each correlation in the correlation matrix. I searched in the R-help mail list and found some related posts, but I didn't find direct way to solve the problem. Someone said to use cor.test() or t.test(). The problem is that cor.test() and t.test() can only apply on two vectors, not on a data.frame or a matrix. My solution is for (i in 1:(ncol(en.data) -1)) { cor.test(en.data[,i], en.data[, i+1]) } I think it is a stupid way. Is there a direct way to do so? After all, it is a basic function to generate significant level of a correlation in a correlation matrix. Thank you in advance! Wei-Wei Hi, Bill Venables posted a solution to this on the R-Help list in Jan 2000. I made a minor modification to add a class to the result and wrote a print method (which could probably do with some tidying but it works). E.g.: # paste in the functions below, then data(iris) corProb(iris[,1:4]) ## prints Correlations are shown below the diagonal P-values are shown above the diagonal Sepal.Length Sepal.Width Petal.Length Petal.Width Sepal.Length 1. 0.1519 0. 0. Sepal.Width -0.1176 1. 0. 0. Petal.Length 0.8718 -0.4284 1. 0. Petal.Width 0.8179 -0.3661 0.9629 1. Is this what you want? HTH G # correlation function # based on post by Bill Venables on R-Help # Date: Tue, 04 Jan 2000 15:05:39 +1000 # https://stat.ethz.ch/pipermail/r-help/2000-January/009758.html # modified by G L Simpson, September 2003 # version 0.2: added print.cor.prob # added class statement to cor.prob # version 0.1: original function of Bill Venables corProb - function(X, dfr = nrow(X) - 2) { R - cor(X) above - row(R) col(R) r2 - R[above]^2 Fstat - r2 * dfr / (1 - r2) R[above] - 1 - pf(Fstat, 1, dfr) class(R) - corProb R } print.corProb - function(x, digits = getOption(digits), quote = FALSE, na.print = , justify = none, ...) { xx - format(unclass(round(x, digits = 4)), digits = digits, justify = justify) if (any(ina - is.na(x))) xx[ina] - na.print cat(\nCorrelations are shown below the diagonal\n) cat(P-values are shown above the diagonal\n\n) print(xx, quote = quote, ...) invisible(x) } -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% *Note new Address and Fax and Telephone numbers from 10th April 2006* %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Gavin Simpson [t] +44 (0)20 7679 0522 ECRC [f] +44 (0)20 7679 0565 UCL Department of Geography Pearson Building [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street London, UK[w] http://www.ucl.ac.uk/~ucfagls/cv/ WC1E 6BT [w] http://www.ucl.ac.uk/~ucfagls/ %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] A possible too old question on significant test of correlationmatrix
you can use function rcor.test() from package 'ltm', e.g., help(rcor.test, package = ltm) ### library(ltm) dat - data.frame(matrix(rnorm(1000), 100, 10)) rcor.test(dat) rcor.test(dat, method = kendall) rcor.test(dat, method = spearman) I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm - Original Message - From: Guo Wei-Wei [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Monday, July 10, 2006 7:27 AM Subject: [R] A possible too old question on significant test of correlationmatrix Dear all, I'm working on a data.frame named en.data, which has n cases and m columns. I generate the correlation matrix of en.data by cor(en.data) I find that there is no p-value on each correlation in the correlation matrix. I searched in the R-help mail list and found some related posts, but I didn't find direct way to solve the problem. Someone said to use cor.test() or t.test(). The problem is that cor.test() and t.test() can only apply on two vectors, not on a data.frame or a matrix. My solution is for (i in 1:(ncol(en.data) -1)) { cor.test(en.data[,i], en.data[, i+1]) } I think it is a stupid way. Is there a direct way to do so? After all, it is a basic function to generate significant level of a correlation in a correlation matrix. Thank you in advance! Wei-Wei __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.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] about overdispersed poisson model
Dear Chi Kai, Three years ago there was a similar thread. At that time David Firth offered a solution for the quasipoission problem with negative observations, see: http://finzi.psych.upenn.edu/R/Rhelp02a/archive/16143.html I remember that his code gave you slightly different answers than the example in England and Verrall's paper. Kind Regards Markus Gesmann FPMA Lloyd's Market Analysis Lloyd's * One Lime Street * London * EC3M 7HA Telephone +44 (0)20 7327 6472 Facsimile +44 (0)20 7327 5718 http://www.lloyds.com -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of g0354502 Sent: 10 July 2006 05:26 To: r-help@stat.math.ethz.ch Subject: [R] about overdispersed poisson model Dear R users I have been looking for functions that can deal with overdispersed poisson models. According to actuarial literature (England Verall, Stochastic Claims Reserving in General Insurance , Institute of Actiuaries 2002) this can be handled through the use of quasi likelihoods instead of normal likelihoods. However, we see them frequently in this type of data, and we would like to be able to fit the model anyway. If it is possible, would you please show me how to find the corresponding package and utilize them? Best Regards, Chi Kai ** The information in this E-Mail and in any attachments is CON...{{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] A possible too old question on significant test of correlation matrix
Hi, Gavin, your program is excellent. Thank you very much! And I have two further questions. 1. Since it is very possible that the data contains missing value and the program will failed against missing values, I have to delete all the cases contained NA. Can it be done pairwisely? 2. Can the program show t values instead of p values? Best regards, Wei-Wei 2006/7/10, Gavin Simpson [EMAIL PROTECTED]: On Mon, 2006-07-10 at 13:27 +0800, Guo Wei-Wei wrote: Dear all, I'm working on a data.frame named en.data, which has n cases and m columns. I generate the correlation matrix of en.data by cor(en.data) I find that there is no p-value on each correlation in the correlation matrix. I searched in the R-help mail list and found some related posts, but I didn't find direct way to solve the problem. Someone said to use cor.test() or t.test(). The problem is that cor.test() and t.test() can only apply on two vectors, not on a data.frame or a matrix. My solution is for (i in 1:(ncol(en.data) -1)) { cor.test(en.data[,i], en.data[, i+1]) } I think it is a stupid way. Is there a direct way to do so? After all, it is a basic function to generate significant level of a correlation in a correlation matrix. Thank you in advance! Wei-Wei Hi, Bill Venables posted a solution to this on the R-Help list in Jan 2000. I made a minor modification to add a class to the result and wrote a print method (which could probably do with some tidying but it works). E.g.: # paste in the functions below, then data(iris) corProb(iris[,1:4]) ## prints Correlations are shown below the diagonal P-values are shown above the diagonal Sepal.Length Sepal.Width Petal.Length Petal.Width Sepal.Length 1. 0.1519 0. 0. Sepal.Width -0.1176 1. 0. 0. Petal.Length 0.8718 -0.4284 1. 0. Petal.Width 0.8179 -0.3661 0.9629 1. Is this what you want? HTH G # correlation function # based on post by Bill Venables on R-Help # Date: Tue, 04 Jan 2000 15:05:39 +1000 # https://stat.ethz.ch/pipermail/r-help/2000-January/009758.html # modified by G L Simpson, September 2003 # version 0.2: added print.cor.prob # added class statement to cor.prob # version 0.1: original function of Bill Venables corProb - function(X, dfr = nrow(X) - 2) { R - cor(X) above - row(R) col(R) r2 - R[above]^2 Fstat - r2 * dfr / (1 - r2) R[above] - 1 - pf(Fstat, 1, dfr) class(R) - corProb R } print.corProb - function(x, digits = getOption(digits), quote = FALSE, na.print = , justify = none, ...) { xx - format(unclass(round(x, digits = 4)), digits = digits, justify = justify) if (any(ina - is.na(x))) xx[ina] - na.print cat(\nCorrelations are shown below the diagonal\n) cat(P-values are shown above the diagonal\n\n) print(xx, quote = quote, ...) invisible(x) } -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% *Note new Address and Fax and Telephone numbers from 10th April 2006* %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Gavin Simpson [t] +44 (0)20 7679 0522 ECRC [f] +44 (0)20 7679 0565 UCL Department of Geography Pearson Building [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street London, UK[w] http://www.ucl.ac.uk/~ucfagls/cv/ WC1E 6BT [w] http://www.ucl.ac.uk/~ucfagls/ %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R vs. other software? (was: Software for Epidmiologica
On 10-Jul-06 Spencer Graves wrote: [...] Excellent overview of many string reasons for R, Spencer! Many thanks for this well-thought-out advocacy piece. [...] At one time, the market dominance of both SAS and IBM was sustained by the perception that 'nobody ever got fired' for using them. Those days are history now for SAS as well as IBM, [...] Does that imply that, by now, some people *have* been fired for using SAS (implicitly, perhaps, because they should have been using R)? :) Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 10-Jul-06 Time: 09:22:43 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Sjava on Windows?
Hello everyone I am very much a newcomer to R and I need to call R processing routines from java and pass back the results to java. It seems sjava is the way to do this, but there seems to be a lot of complexity surrounding configuring R for this purpose. I need to do this under windows XP. Has anyone else had any experience with this and would they be able to help a beginner to get this to work? A colleague passed me this email address and suggested I contact it, apologies if this is not the usual etiquette. Kind Regards Martin Austwick. __ R-help@stat.math.ethz.ch 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] A possible too old question on significant test of correlationmatrix
Thank you, Dimitris. The function rcor.test() is very nice. It can pass arguments to cor() and solve my first problem of pairwised case deletion. Best regards, Wei-Wei 2006/7/10, Dimitris Rizopoulos [EMAIL PROTECTED]: you can use function rcor.test() from package 'ltm', e.g., help(rcor.test, package = ltm) ### library(ltm) dat - data.frame(matrix(rnorm(1000), 100, 10)) rcor.test(dat) rcor.test(dat, method = kendall) rcor.test(dat, method = spearman) I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm - Original Message - From: Guo Wei-Wei [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Monday, July 10, 2006 7:27 AM Subject: [R] A possible too old question on significant test of correlationmatrix Dear all, I'm working on a data.frame named en.data, which has n cases and m columns. I generate the correlation matrix of en.data by cor(en.data) I find that there is no p-value on each correlation in the correlation matrix. I searched in the R-help mail list and found some related posts, but I didn't find direct way to solve the problem. Someone said to use cor.test() or t.test(). The problem is that cor.test() and t.test() can only apply on two vectors, not on a data.frame or a matrix. My solution is for (i in 1:(ncol(en.data) -1)) { cor.test(en.data[,i], en.data[, i+1]) } I think it is a stupid way. Is there a direct way to do so? After all, it is a basic function to generate significant level of a correlation in a correlation matrix. Thank you in advance! Wei-Wei __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.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] Choice of repository and outdated vs unusable... was (Hunting for snow...)
Brian Lunergan wrote: Uwe Ligges wrote: Where did you find them? At least not in their current versions on the official repository on CRAN master for R-2.3.x, I hope. FracSim_0.2.zip http://cran.stat.ucla.edu/bin/windows/contrib/2.3/ RDCOMClient_0.91-0.zip http://www.omegahat.org/R/bin/windows/contrib/R-2.3.0/ segmented_0.1-4.zip http://cran.r-project.org/bin/windows/contrib/2.1/ VGAM_0.6-9.zip http://www.stat.auckland.ac.nz/~yee/bin/windows/contrib/2.3/ None of them is on CRAN master: - RDCOMClient is on Omegahat - segmented is CRAN but outdated for R-2.1.x - VGAM is in a private repository - FracSim: I cannot access http://cran.stat.ucla.edu/bin/windows/contrib/2.3/ right now, but I guess it does not sync properly. At least, the package should not be there. Uwe Ligges Uwe: So VGAM and RDCOMClient are on private repository sites. Is that so important? If a package is not on the CRAN master site are you supposed to avoid using it? Brian, no, but in the mail before you tols us that the packages are not on the check list, and I gave the reason why ... The FracSim home site indicated that 0.2 was the current version and a google search yielded the site I indicated, among others, as a repository of an apparent windows edition. That it's not on the master is more an issue between the maintainers of the master and the maintainers of the package, not the average end user IMHO. I found a The reason is that it does not pass the checks, and that is only an issue of the maintainer of the package, not the CRAN maintainer. That you found the package is the issue of the maintainer of a CRAN mirror that is out of sync. It is an issue of the average end user to install a package from source if he does not get a binary. source for what is apparently the current edition. What's wrong with that? Using the source is perfect. I must have misunderstood that you have compiled from source. But then, be careful and run R CMD check in order to see whether the package works under your platform. Best, Uwe Ligges With regard to segmented, outdated does not mean unusable unless later editions of R or packages that call it as a dependency use features particular to later editions. If there is a newer windows edition of segmented then please point out a source for it. Regards... __ R-help@stat.math.ethz.ch 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] Weighted histograms
Does anyone have any code for drawing weighted histograms (a la Manet/Mondrian) in R? Thanks, 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] par(mfrow,mai) and multiple plot problem
Hi I'm having difficulty with a multiple plot. What I want is 12 plots, all of which are the same size and shape, differing only in colour. Each one is a square, so there is an asp=1 in the plot call. I'm working in an Sweave environment so I am free to choose the height and width of the plot. I want the columns to be labelled at the head (here t=1,2,3) and the rows to be labelled at the left hand side (here an animal), I don't want axes. I want the plots to occupy the majority of the area of the output postscript. I want the vertical distance between the plots to be the same as the horizontal distance between the plots (not critical but would be nice). Can anyone help me with these requirements? My best attempt is below, but it is no good because the ylab argument that I'm trying to use for the row labelling gets clipped. Also, I can't get my horizontal spacing to be equal to my vertical spacing (although the code below isn't _too_ bad). f - function(...){ jj - expand.grid(1:10,1:10) colnames(jj) - c(,) par(mai=rep(0,4)+0.2) plot(jj,pch=16,asp=1, axes=FALSE, ...) } postscript(file=~/f.ps,width=5,height=7) par(mfrow=c(4,3)) f(main=t=1,ylab=fish) f(main=t=2) f(main=t=3) f(ylab=dog) f() f() f(ylab=slug) f() f() f(ylab=pig) f() f() dev.off() -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Sjava on Windows?
(R)SJava is part of the Omegahat project, not the R project. In theory at least, Omegahat has its own mailing lists. In any case, the R posting guide asks you to contact the maintainer of a contributed package before this list, so please do so now. There are easier-to-setup alternatives to using R from Java: see http://www.rosuda.org/software/ for example (some of which software is on CRAN). (From your description Rserve may suffice: an alternative would be to use DCOM.) On Mon, 10 Jul 2006, Dr Martin Austwick wrote: Hello everyone I am very much a newcomer to R and I need to call R processing routines from java and pass back the results to java. It seems sjava is the way to do this, but there seems to be a lot of complexity surrounding configuring R for this purpose. I need to do this under windows XP. Has anyone else had any experience with this and would they be able to help a beginner to get this to work? A colleague passed me this email address and suggested I contact it, apologies if this is not the usual etiquette. Kind Regards Martin Austwick. __ R-help@stat.math.ethz.ch 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] par(mfrow,mai) and multiple plot problem
[apologies for possible multiple posting] Hi Joris great suggestion! I never thought to use layout() in this way. If I have postscript(file=~/f.ps,width=5,height=8) nf - layout(matrix( c(00,01,00,02,00,03, 04,05,00,06,00,07, 00,00,00,00,00,00, 08,09,00,10,00,11, 00,00,00,00,00,00, 12,13,00,14,00,15, 00,00,00,00,00,00, 16,17,00,18,00,19),8,6,byrow=TRUE), c(1,4, 1,4, 1,4), c(1,4,1,4,1,4,1,4), TRUE) layout.show(nf) dev.off() This is exactly what I want. I hadn't appreciated that the empty rows and columns would be clear. How best to put text in boxes 1,2,3 (that used to be main() and vertical text in boxes 4,8,12,16 (that used to be ylab)? best wishes Robin On 10 Jul 2006, at 10:27, Joris De Wolf wrote: Check ?layout It gives you more flexibility than par(mfrow=c(a,b)) For instance, you can define your margins between the graphs as cells in the layout without filling them by a plot. Joris Robin Hankin wrote: Hi I'm having difficulty with a multiple plot. -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] pvclust missing values problem
Hello all, I posted a question to this list last week and received no response. I am unsure if this means no-one knows the answer or if I posed the question badly. I'm going to assume I posed the question badly and try again. I am new to R so it is quite likely it's a very naive question, however if there is something blindingly obvious that I am missing or if there is another resource I should consult that I haven't seen would someone be kind enough to point it out because it isn't obvious to me. Although my data is from biological experiments I think my problem is with R rather than the nature of the data, but I may be wrong. I am attempting to use the pvclust package to do some hierarchical clustering on some CGH data I have downloaded from the Progenetix database (http://www.progenetix.de/~pgscripts/progenetix/Aboutprogenetix.html). The data is in tab delimited format, each column is a single sample each row is a chromosome band some example dummy data is shown below. band sample1 sample2 sample3 sample4 1p36_33 1 00 1 1p36_32 -1 0-1 0 1p36_31 0 11 1 1p36_220-1 -1 -1 etc where 0 = no change, 1 = gain, -1 = loss I have read this file into R using: ProgenetixCRC.all.noXY - read.table(/home/marraydb/Progenetix/Data/CRCall_noXY.txt, header=TRUE, sep=\t, row.names=band) based on the pvclust documentation I came up with this: ProgenetixCRC.all.pvclust - pvclust(ProgenetixCRC.all, method.dist=cor, method.hclust=average,use.cor=pairwise.complete.obs,nboot=1000) this results in an error Error in hclust(distance, method = method.hclust) : NA/NaN/Inf in foreign function call (arg 11) Digging through the mailing list archives I've discovered this means that my dataset has missing values. This is very confusing because I have checked and there are no missing values. Running is.na() over the data matrix results in all false values which I take to mean none of the values are NA. I tried various options for the use.cor argument all with the same result. Since I originally posted this question I tried changing method.dist to euclidean, in this form the function executes without any errors. This is not to say the results actually mean anything of course. I am at a loss as to how to proceed any input from someone more experienced would be gratefully appreciated. If there is some reason why I should not be doing this analysis this way in the first place then I'd appreciate having that pointed out also. I've tried not to put excess information in here but if more is needed then let me know what and I'll post it. I suspect the problem is me, however if it really is the case that no-one knows how to answer this then could anyone suggest another mailing list where I might get a better response. Would bioconductor be a better option for example? Apologies for any offence caused by posting the same question but it's difficult for me to proceed until I get some kind of response, even if it is that this list is not the right place for this question. Thanks for your patience, Richard Dr Richard Birnie Scientific Officer Section of Pathology and Tumour Biology Welcome Brenner Building, LIMM St James University Hospital Beckett St, Leeds, LS9 7TF Tel:0113 3438624 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
Re: [R] A possible too old question on significant test of correlation matrix
On Mon, 2006-07-10 at 16:22 +0800, Guo Wei-Wei wrote: Hi, Gavin, your program is excellent. Thank you very much! And I have two further questions. 1. Since it is very possible that the data contains missing value and the program will failed against missing values, I have to delete all the cases contained NA. Can it be done pairwisely? Yes, with a modification to accept and pass on argument use, e.g.: data(iris) ## copy data iris2 - iris ## simulate some missing values in Sepal.Length iris2[sample(1:nrow(iris2), 5), 1] - NA ## corProb matrix with missing values temp - corProb(iris2[,1:4], use = pairwise.complete.obs) See ?cor for the options you can specify for use. You'll need to paste in the functions below for this to work. 2. Can the program show t values instead of p values? Yes - this is R! The function Bill Venables wrote uses F-values, so I looked at what cor.test was doing and modified the function to compute either t or F values and to return them or their p-values. temp - corProb(iris[,1:4]) temp - corProb(iris[,1:4], type = t) temp - corProb(iris[,1:4], type = t, pval = FALSE) For t-values you can do the different test as in ?cor.test : ## with different tests temp - corProb(iris[,1:4], type = t, alternative = less) temp - corProb(iris[,1:4], type = t, alternative = greater) Hopefully that does what you wanted. G ## New functions. corProb - function(X, dfr = nrow(X) - 2, use = c(all.obs, complete.obs, pairwise.complete.obs), alternative = c(two.sided, less, greater), type = c(F, t), pval = TRUE) { USE - match.arg(use) ALTERNATIVE - match.arg(alternative) R - cor(X, use = USE) above - row(R) col(R) r2 - R[above]^2 TYPE - match.arg(type) if(TYPE == t) { Tstat - sqrt(dfr) * R[above]/sqrt(1 - r2) if(pval) { p - pt(Tstat, dfr) R[above] - switch(ALTERNATIVE, less = p, greater = 1 - p, two.sided = 2 * min(p, 1 - p)) } else R[above] - Tstat } else { Fstat - r2 * dfr / (1 - r2) if(pval) R[above] - 1 - pf(Fstat, 1, dfr) else R[above] - Fstat } class(R) - corProb attr(R, type) - TYPE attr(R, pval) - pval attr(R, hypoth) - ALTERNATIVE R } print.corProb - function(x, digits = getOption(digits), quote = FALSE, na.print = , justify = none, ...) { xx - format(unclass(round(x, digits = 4)), digits = digits, justify = justify) if (any(ina - is.na(x))) xx[ina] - na.print cat(\nCorrelations are shown below the diagonal\n) if(attr(x, pval)) cat(paste(P-values of the , attr(x, type), -statistics are shown above the diagonal\n\n, sep = )) else cat(paste(attr(x, type), -values are shown above the diagonal\n\n, sep = )) if(attr(x, type) == t) hypoth - switch(attr(x, hypoth), less = less than 0, greater = greater than 0, two.sided = not equal to 0) cat(paste(alternative hypothesis: true correlation is, hypoth, \n\n)) print.default(xx, quote = quote, ...) invisible(x) } -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Gavin Simpson [t] +44 (0)20 7679 0522 ECRC ENSIS, UCL Geography, [f] +44 (0)20 7679 0565 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/cv/ London, UK. WC1E 6BT. [w] http://www.ucl.ac.uk/~ucfagls/ %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] tree-ring response function in R?
Hi, I am looking forward if some of you could help me with this. I have been looking R-functions to conduct tree-ring analysis in R. I know that exists a COFECHA version in R developed by C. Bigler. But I don't know if is there another collection of functions to do dendroclimatological analysis (such as response functions). Has someone developed such rutines in R ?, thanks in advance, Alvaro __ R-help@stat.math.ethz.ch 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] pvclust missing values problem
On Mon, 2006-07-10 at 11:56 +0100, Richard Birnie wrote: Hello all, Hi Richard, Sorry, I know nothing about pvclust and have never used it, but here are a couple of general suggestions/observations. You are asked to contact the package maintainer *not* R-Help for questions relating to contributed packages. I dare say only the maintainer (or some kind soul with too much time on their hands to debug the package) can solve your problem, if you can't work out how to do it with R's debug tools. I have cc'd the maintainer on this reply. When you get an error, it is useful to supply the output from traceback() to see where the error actually occurred. By the way, the error NA/NaN/Inf in foreign function call (arg 11) doesn't necessarily mean you have missing values in your data set. Note the NaN and Inf in that error message. It could just be that one of the calculations resulted in NaN's or Inf's which hclust detected or caught and issued the error. Without traceback(), this is pure speculation. R version info, platform (OS) and version of pvclust are all useful extra bits of info that you could have provided us or the maintainer as well. Hopefully you can sort your problem out with the maintainer. Best, G I posted a question to this list last week and received no response. I am unsure if this means no-one knows the answer or if I posed the question badly. I'm going to assume I posed the question badly and try again. I am new to R so it is quite likely it's a very naive question, however if there is something blindingly obvious that I am missing or if there is another resource I should consult that I haven't seen would someone be kind enough to point it out because it isn't obvious to me. Although my data is from biological experiments I think my problem is with R rather than the nature of the data, but I may be wrong. I am attempting to use the pvclust package to do some hierarchical clustering on some CGH data I have downloaded from the Progenetix database (http://www.progenetix.de/~pgscripts/progenetix/Aboutprogenetix.html). The data is in tab delimited format, each column is a single sample each row is a chromosome band some example dummy data is shown below. band sample1 sample2 sample3 sample4 1p36_33 1 00 1 1p36_32 -1 0-1 0 1p36_31 0 11 1 1p36_220-1 -1 -1 etc where 0 = no change, 1 = gain, -1 = loss I have read this file into R using: ProgenetixCRC.all.noXY - read.table(/home/marraydb/Progenetix/Data/CRCall_noXY.txt, header=TRUE, sep=\t, row.names=band) based on the pvclust documentation I came up with this: ProgenetixCRC.all.pvclust - pvclust(ProgenetixCRC.all, method.dist=cor, method.hclust=average,use.cor=pairwise.complete.obs,nboot=1000) this results in an error Error in hclust(distance, method = method.hclust) : NA/NaN/Inf in foreign function call (arg 11) Digging through the mailing list archives I've discovered this means that my dataset has missing values. This is very confusing because I have checked and there are no missing values. Running is.na() over the data matrix results in all false values which I take to mean none of the values are NA. I tried various options for the use.cor argument all with the same result. Since I originally posted this question I tried changing method.dist to euclidean, in this form the function executes without any errors. This is not to say the results actually mean anything of course. I am at a loss as to how to proceed any input from someone more experienced would be gratefully appreciated. If there is some reason why I should not be doing this analysis this way in the first place then I'd appreciate having that pointed out also. I've tried not to put excess information in here but if more is needed then let me know what and I'll post it. I suspect the problem is me, however if it really is the case that no-one knows how to answer this then could anyone suggest another mailing list where I might get a better response. Would bioconductor be a better option for example? Apologies for any offence caused by posting the same question but it's difficult for me to proceed until I get some kind of response, even if it is that this list is not the right place for this question. Thanks for your patience, Richard -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Gavin Simpson [t] +44 (0)20 7679 0522 ECRC ENSIS, UCL Geography, [f] +44 (0)20 7679 0565 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/cv/ London, UK. WC1E 6BT. [w] http://www.ucl.ac.uk/~ucfagls/ %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ R-help@stat.math.ethz.ch mailing list
[R] How to include NA's of a factor in table?
Dear All, Is there a better way to include NA's of a factor in the output of table() than using as.character()? Admittedly, I do not understand the help page for table concerning the exclude argument applied to factors. I tried in different ways, but could not get NA to be included in the table, if not using as.character() (see example). Greetings, Heinz ## example fcv - factor(c('a', NA, 'c')) table(fcv)# shows a, c table(fcv, exclude='a') # shows c table(fcv, exclude=)# shows a, c table(fcv, exclude=NULL) # shows a, c table(as.character(fcv), exclude=NULL) # shows a, c, NA platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status Patched major 2 minor 3.1 year 2006 month 07 day01 svn rev38471 language R version.string Version 2.3.1 Patched (2006-07-01 r38471) __ R-help@stat.math.ethz.ch 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 include NA's of a factor in table?
fcv - factor(c('a', NA, 'c'), exclude=NULL) table(fcv, exclude=NULL) --- Jacques VESLOT CNRS UMR 8090 I.B.L (2ème étage) 1 rue du Professeur Calmette B.P. 245 59019 Lille Cedex Tel : 33 (0)3.20.87.10.44 Fax : 33 (0)3.20.87.10.31 http://www-good.ibl.fr --- Heinz Tuechler a écrit : Dear All, Is there a better way to include NA's of a factor in the output of table() than using as.character()? Admittedly, I do not understand the help page for table concerning the exclude argument applied to factors. I tried in different ways, but could not get NA to be included in the table, if not using as.character() (see example). Greetings, Heinz ## example fcv - factor(c('a', NA, 'c')) table(fcv)# shows a, c table(fcv, exclude='a') # shows c table(fcv, exclude=)# shows a, c table(fcv, exclude=NULL) # shows a, c table(as.character(fcv), exclude=NULL) # shows a, c, NA platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status Patched major 2 minor 3.1 year 2006 month 07 day01 svn rev38471 language R version.string Version 2.3.1 Patched (2006-07-01 r38471) __ R-help@stat.math.ethz.ch 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] How to include NA's of a factor in table?
If we specify exclude=NULL in factor then it need not also be specified in table: fcv - factor(c('a', NA, 'c'), exclude=NULL) table(fcv) On 7/10/06, Jacques VESLOT [EMAIL PROTECTED] wrote: fcv - factor(c('a', NA, 'c'), exclude=NULL) table(fcv, exclude=NULL) --- Jacques VESLOT CNRS UMR 8090 I.B.L (2ème étage) 1 rue du Professeur Calmette B.P. 245 59019 Lille Cedex Tel : 33 (0)3.20.87.10.44 Fax : 33 (0)3.20.87.10.31 http://www-good.ibl.fr --- Heinz Tuechler a écrit : Dear All, Is there a better way to include NA's of a factor in the output of table() than using as.character()? Admittedly, I do not understand the help page for table concerning the exclude argument applied to factors. I tried in different ways, but could not get NA to be included in the table, if not using as.character() (see example). Greetings, Heinz ## example fcv - factor(c('a', NA, 'c')) table(fcv)# shows a, c table(fcv, exclude='a') # shows c table(fcv, exclude=)# shows a, c table(fcv, exclude=NULL) # shows a, c table(as.character(fcv), exclude=NULL) # shows a, c, NA platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status Patched major 2 minor 3.1 year 2006 month 07 day01 svn rev38471 language R version.string Version 2.3.1 Patched (2006-07-01 r38471) __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] pvclust missing values problem
Thanks for the responses, Now I know why I got no response I'll know better in the future. Apologies all round for any inconvenience caused that was never my intention, just inexperience. Richard Dr Richard Birnie Scientific Officer Section of Pathology and Tumour Biology Welcome Brenner Building, LIMM St James University Hospital Beckett St, Leeds, LS9 7TF Tel:0113 3438624 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
Re: [R] distance in kmeans algorithm?
Hello. Is it possible to choose the distance in the kmeans algorithm? I have m vectors of n components and I want to cluster them using kmeans algorithm but I want to use the Mahalanobis distance or another distance. How can I do it in R? If I use kmeans, I have no option to choose the distance. Thanks in advance, Arnau. You can use Kmeans from the amap package with several distance measures. # example for L1 and L2: x - matrix(c(0,0,0,1.5,1,-1), ncol=2, byrow=TRUE) require(amap) Kmeans(x, x[2:3,], method=manhattan) Kmeans(x, x[2:3,], method=euclidean) Cheers, Timo -- Timo Becker Phonetics Austrian Academy of Sciences Acoustics Research Institute __ R-help@stat.math.ethz.ch 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] A possible too old question on significant test of correlation matrix
On Mon, 2006-07-10 at 12:48 +0100, Gavin Simpson wrote: On Mon, 2006-07-10 at 16:22 +0800, Guo Wei-Wei wrote: Hi, Gavin, your program is excellent. Thank you very much! And I have two further questions. 1. Since it is very possible that the data contains missing value and the program will failed against missing values, I have to delete all the cases contained NA. Can it be done pairwisely? Yes, with a modification to accept and pass on argument use, e.g.: data(iris) ## copy data iris2 - iris ## simulate some missing values in Sepal.Length iris2[sample(1:nrow(iris2), 5), 1] - NA ## corProb matrix with missing values temp - corProb(iris2[,1:4], use = pairwise.complete.obs) See ?cor for the options you can specify for use. You'll need to paste in the functions below for this to work. 2. Can the program show t values instead of p values? Yes - this is R! The function Bill Venables wrote uses F-values, so I looked at what cor.test was doing and modified the function to compute either t or F values and to return them or their p-values. Oops, there was a simple error in the print method. Fixed below: print.corProb - function(x, digits = getOption(digits), quote = FALSE, na.print = , justify = none, ...) { xx - format(unclass(round(x, digits = 4)), digits = digits, justify = justify) if (any(ina - is.na(x))) xx[ina] - na.print cat(\nCorrelations are shown below the diagonal\n) if(attr(x, pval)) cat(paste(P-values of the , attr(x, type), -statistics are shown above the diagonal\n\n, sep = )) else cat(paste(attr(x, type), -values are shown above the diagonal\n\n, sep = )) if(attr(x, type) == t) { hypoth - switch(attr(x, hypoth), less = less than 0, greater = greater than 0, two.sided = not equal to 0) cat(paste(alternative hypothesis: true correlation is, hypoth, \n\n)) } print.default(xx, quote = quote, ...) invisible(x) } -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Gavin Simpson [t] +44 (0)20 7679 0522 ECRC ENSIS, UCL Geography, [f] +44 (0)20 7679 0565 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/cv/ London, UK. WC1E 6BT. [w] http://www.ucl.ac.uk/~ucfagls/ %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] How can I obtain the values of BIAS and STD. ERROR from a bootstrap.
R 2.3.1 Windows XP Question: How can I obtain the values of BIAS and STD. ERROR from a bootstrap. Background: I am running a bootstrap: result2-boot(1:400,regSEvssample,R=5000) and obtain the following results: Bootstrap Statistics : originalbiasstd. error t1* 1.876602 -0.0001368616 0.1630380 I would like to get the values of ORIGINAL and BIAS. I can get the value of ORIGINAL: result2$t0 [1] 1.876602 But I can't find a way to get the values of BIAS or STD. ERROR. Can someone suggest how I can obtain the values for these two statistics? Thanks, John John Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics Baltimore VA Medical Center GRECC, University of Maryland School of Medicine Claude D. Pepper OAIC, University of Maryland Clinical Nutrition Research Unit, and Baltimore VA Center Stroke of Excellence University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 410-605-7119 [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] How can I obtain the values of BIAS and STD. ERROR from a bootstrap.
Look at the code of boot:::print.boot , which does the calculations. On Mon, 10 Jul 2006, John Sorkin wrote: R 2.3.1 Windows XP Question: How can I obtain the values of BIAS and STD. ERROR from a bootstrap. Background: I am running a bootstrap: result2-boot(1:400,regSEvssample,R=5000) and obtain the following results: Bootstrap Statistics : originalbiasstd. error t1* 1.876602 -0.0001368616 0.1630380 I would like to get the values of ORIGINAL and BIAS. I can get the value of ORIGINAL: result2$t0 [1] 1.876602 But I can't find a way to get the values of BIAS or STD. ERROR. Can someone suggest how I can obtain the values for these two statistics? Thanks, John John Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics Baltimore VA Medical Center GRECC, University of Maryland School of Medicine Claude D. Pepper OAIC, University of Maryland Clinical Nutrition Research Unit, and Baltimore VA Center Stroke of Excellence University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 410-605-7119 [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 -- 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 include NA's of a factor in table?
Thank you Jacques and Gabor! Your solution does work and it led me to try also: table(factor(fcv, exclude=NULL)) # shows a, c, NA Greetings, Heinz At 08:34 10.07.2006 -0400, Gabor Grothendieck wrote: If we specify exclude=NULL in factor then it need not also be specified in table: fcv - factor(c('a', NA, 'c'), exclude=NULL) table(fcv) On 7/10/06, Jacques VESLOT [EMAIL PROTECTED] wrote: fcv - factor(c('a', NA, 'c'), exclude=NULL) table(fcv, exclude=NULL) --- Jacques VESLOT CNRS UMR 8090 I.B.L (2ème étage) 1 rue du Professeur Calmette B.P. 245 59019 Lille Cedex Tel : 33 (0)3.20.87.10.44 Fax : 33 (0)3.20.87.10.31 http://www-good.ibl.fr --- Heinz Tuechler a écrit : Dear All, Is there a better way to include NA's of a factor in the output of table() than using as.character()? Admittedly, I do not understand the help page for table concerning the exclude argument applied to factors. I tried in different ways, but could not get NA to be included in the table, if not using as.character() (see example). Greetings, Heinz ## example fcv - factor(c('a', NA, 'c')) table(fcv)# shows a, c table(fcv, exclude='a') # shows c table(fcv, exclude=)# shows a, c table(fcv, exclude=NULL) # shows a, c table(as.character(fcv), exclude=NULL) # shows a, c, NA platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status Patched major 2 minor 3.1 year 2006 month 07 day01 svn rev38471 language R version.string Version 2.3.1 Patched (2006-07-01 r38471) __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Query:chi-squre test
Hi, I have calculated chi-square goodness of fit test,Sample coming from Poisson distribution. please copy this script in R run the script The R script is as follows ## start # No_of_Frauds- c(4,1,6,9,9,10,2,4,8,2,3,0,1,2,3,1,3,4,5,4,4,4,9,5,4,3,11,8,12,3,10,0,7) lambda- mean(No_of_Frauds) # Chi-Squared Goodness of Fit Test # Ho: The data follow a specified distribution Vs H1: Not Ho # observed frequencies variable.cnts - table(No_of_Frauds) variable.cnts variable.cnts.prs - dpois(as.numeric(names(variable.cnts)), lambda) variable.cnts.prs variable.cnts - c(variable.cnts, 0) variable.cnts variable.cnts.prs - c(variable.cnts.prs, 1-sum(variable.cnts.prs)) variable.cnts.prs tst - chisq.test(variable.cnts, p=variable.cnts.prs) Tst # end The result of R is as follows Warning message: Chi-squared approximation may be incorrect in: chisq.test(variable.cnts, p = variable.cnts.prs) tst Chi-squared test for given probabilities data: variable.cnts X-squared = 40.5614, df = 13, p-value = 0.0001122 But I have done calculations in Excel. I am getting different answer. Observed = 2,3,3,5,7,2,1,1,2,3,2,1,1,0 Expected=0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.78 1782103,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.23 4400679,0.095299266,0.035764993 Estimated Parameter =4.878788 Chi square stat = 0.000113 My excel answer tally with the book which I have refer for excel. Please tell me the correct calculation in R. And how to interprit the results in R. Thanks. Regards. Priti. __ R-help@stat.math.ethz.ch 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] retrieving object name
I have a data frame with named columns and I would like to know if it is possible to retrieve a column name once selected: print(colnames(df)) # assumes to print col1 col2 print.name(df$col1) # would like to print col1 print.name(df$col2) # would like to print col2 So what the print.name function should do? My aim is not to print the column name but to select some settings from the column name withing the function (i.e. print.name), while this function is applied to several columns of the list/data.frame. Actually, I solved the problem by providing an extra parameter like: print.name(df$col1, col1) but since I may have many of these columns/parameters combination, this is rather error prone and it would be much better if I could detect which columns of the data frame I am dealing with. Thanks. ld. __ R-help@stat.math.ethz.ch 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] Query: chi square test
Hi, I have calculated chi-square goodness of fit test,Sample coming from Poisson distribution. please copy this script in R run the script The R script is as follows ## start # No_of_Frauds- c(4,1,6,9,9,10,2,4,8,2,3,0,1,2,3,1,3,4,5,4,4,4,9,5,4,3,11,8,12,3,10,0,7) lambda- mean(No_of_Frauds) # Chi-Squared Goodness of Fit Test # Ho: The data follow a specified distribution Vs H1: Not Ho # observed frequencies variable.cnts - table(No_of_Frauds) variable.cnts variable.cnts.prs - dpois(as.numeric(names(variable.cnts)), lambda) variable.cnts.prs variable.cnts - c(variable.cnts, 0) variable.cnts variable.cnts.prs - c(variable.cnts.prs, 1-sum(variable.cnts.prs)) variable.cnts.prs tst - chisq.test(variable.cnts, p=variable.cnts.prs) Tst # end The result of R is as follows Warning message: Chi-squared approximation may be incorrect in: chisq.test(variable.cnts, p = variable.cnts.prs) tst Chi-squared test for given probabilities data: variable.cnts X-squared = 40.5614, df = 13, p-value = 0.0001122 But I have done calculations in Excel. I am getting different answer. Observed = 2,3,3,5,7,2,1,1,2,3,2,1,1,0 Expected=0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.78 1782103,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.23 4400679,0.095299266,0.035764993 Estimated Parameter =4.878788 Chi square stat = 0.000113 My excel answer tally with the book which I have refer for excel. Please tell me the correct calculation in R. And how to interprit the results in R. Because actually data should fit for Poisson dist. Thanks. Regards. Priti. __ R-help@stat.math.ethz.ch 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] retrieving object name
On Mon, 10 Jul 2006, Laurent Deniau wrote: I have a data frame with named columns and I would like to know if it is possible to retrieve a column name once selected: Not really. df$col1 is a new object which does not know where it came from. If you wanted to do this before selection, then print.name - function(df, col) struture(df[[col]], from=col) would do the extraction and set an attribute to be consulted later. print(colnames(df)) # assumes to print col1 col2 print.name(df$col1) # would like to print col1 print.name(df$col2) # would like to print col2 So what the print.name function should do? My aim is not to print the column name but to select some settings from the column name withing the function (i.e. print.name), while this function is applied to several columns of the list/data.frame. Actually, I solved the problem by providing an extra parameter like: print.name(df$col1, col1) but since I may have many of these columns/parameters combination, this is rather error prone and it would be much better if I could detect which columns of the data frame I am dealing with. Thanks. ld. -- 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] parametric proportional hazard regression
On Fri, 7 Jul 2006, Valentin Dimitrov wrote: I do not need a accelerated failure model, but a proportional hazard model with a f0= weibull, exponential, loglogistic or lognormal baseline distribution. The hazard function is lambda(t)=exp(Xi*beta)*lambda0(t), where lambda0 is the baseline hazard lambda0(t)=f0(t)/(1-F0(t)) where f0 and F0 are the baseline density and cumulative distribution functions. This is a proportional hazard model since the ratio lambda(t|Xi)/lambda(t|Xj)=exp(Xi*beta)/exp(Xj*beta) does not depend on t. For a weibull (including exponential) model you can do this with survreg. For the other models you would have to maximize the likelihood directly. This will involve writing the likelihood directly in terms of the hazard and cumulative hazard, since a proportional hazards model that is gaussian at X=0 is not gaussian at any other X. -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
Re: [R] survfit, unused argument(s) (error ...)
On Sat, 8 Jul 2006, S?ren Merser wrote: Hi It seems that survfit() doesn't accept the argumnet 'error' as below survfit(fit, error='greenwood') Error in survfit.coxph(fit, error = greenwood) : unused argument(s) (error ...) Isn't is allowed to do that for a coxph object? Looking at the code it seems that you need vartype=greenwood, so the documentation is wrong. -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
Re: [R] Query:chi-squre test
priti desai [EMAIL PROTECTED] writes: Hi, I have calculated chi-square goodness of fit test,Sample coming from Poisson distribution. please copy this script in R run the script The R script is as follows ## start # No_of_Frauds- c(4,1,6,9,9,10,2,4,8,2,3,0,1,2,3,1,3,4,5,4,4,4,9,5,4,3,11,8,12,3,10,0,7) lambda- mean(No_of_Frauds) # Chi-Squared Goodness of Fit Test # Ho: The data follow a specified distribution Vs H1: Not Ho # observed frequencies variable.cnts - table(No_of_Frauds) variable.cnts variable.cnts.prs - dpois(as.numeric(names(variable.cnts)), lambda) variable.cnts.prs variable.cnts - c(variable.cnts, 0) variable.cnts variable.cnts.prs - c(variable.cnts.prs, 1-sum(variable.cnts.prs)) variable.cnts.prs tst - chisq.test(variable.cnts, p=variable.cnts.prs) Tst # end The result of R is as follows Warning message: Chi-squared approximation may be incorrect in: chisq.test(variable.cnts, p = variable.cnts.prs) tst Chi-squared test for given probabilities data: variable.cnts X-squared = 40.5614, df = 13, p-value = 0.0001122 But I have done calculations in Excel. I am getting different answer. Observed = 2,3,3,5,7,2,1,1,2,3,2,1,1,0 Expected=0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.78 1782103,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.23 4400679,0.095299266,0.035764993 Estimated Parameter =4.878788 Chi square stat = 0.000113 My excel answer tally with the book which I have refer for excel. Please tell me the correct calculation in R. And how to interprit the results in R. As far as I can see, the Chi square stat in Excel is essentially the p-value in R. The slight difference appears to arise from Excel using the point probability rather than the tail ditto in the last cell: O - c(2,3,3,5,7,2,1,1,2,3,2,1,1,0) E - c(0.251005528,1.224602726,2.987288468,4.85811559,5.925428863, + 5.781782103,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456, + 0.234400679,0.095299266,0.035764993) (O-E)^2/E [1] 1.218691e+01 2.573925e+00 5.409021e-05 4.143826e-03 1.948725e-01 [6] 2.473610e+00 2.914053e+00 1.581883e+00 1.465377e-06 3.391598e+00 [11] 4.097178e+00 2.500600e+00 8.588560e+00 3.576499e-02 sum((O-E)^2/E) [1] 40.54315 pchisq(sum((O-E)^2/E), 13,low=F) [1] 0.0001129818 E [1] 0.25100553 1.22460273 2.98728847 4.85811559 5.92542886 5.78178210 [7] 4.70134807 3.27669714 1.99828879 1.08324746 0.52849346 0.23440068 [13] 0.09529927 0.03576499 sum(E) [1] 32.98176 Please don't assume that something is correct, just because it is Excel output and some book mindlessly copied it... The calculations are both wrong, because they ignore the fact that lambda has been estimated from the data, and also because they deal with very small expected cell counts. For a better test, you likely need to simulate the distribution of the chi-square, or, as I'd be inclined to do, go directly for the pretty obvious overdispersion: var(X) [1] 11.17235 var(X)/mean(X) # expected is ca. 1 in the Poisson distrib. [1] 2.289984 r - replicate(10,{X - rpois(33, 4.87879); var(X)/mean(X)}) sum(r 2.289984) [1] 5 -- 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] retrieving object name
Use drop = FALSE. For example using builtin data frame BOD we can display the Time column with its heading: BOD[, Time, drop = FALSE] On 7/10/06, Laurent Deniau [EMAIL PROTECTED] wrote: I have a data frame with named columns and I would like to know if it is possible to retrieve a column name once selected: print(colnames(df)) # assumes to print col1 col2 print.name(df$col1) # would like to print col1 print.name(df$col2) # would like to print col2 So what the print.name function should do? My aim is not to print the column name but to select some settings from the column name withing the function (i.e. print.name), while this function is applied to several columns of the list/data.frame. Actually, I solved the problem by providing an extra parameter like: print.name(df$col1, col1) but since I may have many of these columns/parameters combination, this is rather error prone and it would be much better if I could detect which columns of the data frame I am dealing with. Thanks. ld. __ R-help@stat.math.ethz.ch 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 in my code
On Mon, 10 Jul 2006, Gabor Grothendieck wrote: The problem can be reduced to this: x - 1 x[1] - 2 # error The following are ok: x - 1 x[1] - 3 x - 1 x - 4 x - 1 x - 5 Does anyone know why? Is this a bug in - ? No, it's a feature. The fact that x-5 works is arguably a bug (though probably not worth fixing). x[1] - 2 is equivalent (per section 3.4.4 of the language guide) to `*tmp*` - get(x, envir=parent.env(), inherits=TRUE) `*tmp*`[1] - 2 x - `*tmp*` and the get() fails when you try to do this from the command line. Since the point of superassignment is to assign in a lexical parent environment it makes no sense to do it directly at the command line. -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
Re: [R] package:Matrix handling of data with identical indices
roger == roger koenker [EMAIL PROTECTED] on Sun, 9 Jul 2006 12:31:16 -0500 writes: roger On Jul 9, 2006, at 11:06 AM, Douglas Bates wrote: Your matrix Mc should be flagged as invalid. Martin and I should discuss whether we want to add such a test to the validity method. It is not difficult to add the test but there will be a penalty in that it will slow down all operations on such matrices and I'm not sure if we want to pay that price to catch a rather infrequently occuring problem. roger Elaborating the validity procedure to flag such roger instances seems to be well worth the speed penalty in roger my view. Of course, anticipating every such misstep roger imposes a heavy burden on developers and constitutes roger the real cost of more elaborate validity checking. As I found, we already *have* a validate_dgCMatrix in C code, and adding an improved test for the validity of the 'p' slot, solves ``all problems'' mentioned above --- without any performance penalty. Hence., in the upcoming next version of 'Matrix' (0.95-12), John will get a proper error message immediately from calling new(...) with the wrong 'p' (or 'Dim'). Martin roger [My 2cents based on experience with SparseM.] roger url: www.econ.uiuc.edu/~roger Roger Koenker email roger [EMAIL PROTECTED] Department of Economics vox: roger 217-333-4558 University of Illinois fax: 217-244-6678 roger Champaign, IL 61820 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] retrieving object name
Prof Brian Ripley wrote: On Mon, 10 Jul 2006, Laurent Deniau wrote: I have a data frame with named columns and I would like to know if it is possible to retrieve a column name once selected: Not really. df$col1 is a new object which does not know where it came from. If you wanted to do this before selection, then print.name - function(df, col) struture(df[[col]], from=col) would do the extraction and set an attribute to be consulted later. This is an appropriate solution since I build the list myself (by splitting the data frame following some factors) but I do not call the function myself (e.g. print.name). Here is the final steps of the function which read and set the data.frame/list: magnet.read.rt - function (filename) { # ... # convert to list by magnet states MG - split(MG,MG$magnet_state) # set state attribute for automatic settings for(col in names(MG)) MG[[col]] - structure(MG[[col]], state=col) invisible(MG) } and an example of its use: .wc.scl - function (MG) { c(CC=0.85, CM=1.0)[attr(MG,state)] } wc.offset - function (CM,CR) { scl - .wc.scl(CM) #... } MB - magnet.read.rt(magnet.dat) CC.INJ - wc.offset(MB$CC,MB$INJ) Thanks. a+, ld. __ R-help@stat.math.ethz.ch 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] 10^x instead 10EX on plot axes. How?
Hi, I'm drawing a very simple plot with both axes logarithmic (default base 10). Example: vec=c(1,10,100,1000,1,10,100,1000) plot(vec,vec,log=xy) The axes on the plot now show the technical notation like 1E+3 but I would prefer to have it the notation 10 ^3 i.e. with the exponent here 3 superscript (raised). Any help very much appreciated! Best Regards Tom -- Feel free – 10 GB Mailbox, 100 FreeSMS/Monat ... __ R-help@stat.math.ethz.ch 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] 10^x instead 10EX on plot axes. How?
On 7/10/2006 12:43 PM, [EMAIL PROTECTED] wrote: Hi, I'm drawing a very simple plot with both axes logarithmic (default base 10). Example: vec=c(1,10,100,1000,1,10,100,1000) plot(vec,vec,log=xy) The axes on the plot now show the technical notation like 1E+3 but I would prefer to have it the notation 10 ^3 i.e. with the exponent here 3 superscript (raised). Any help very much appreciated! You can use the plotmath functions for axis labels. For example, vec=c(1,10,100,1000,1,10,100,1000) plot(vec,vec,log=xy, axes=F) axis(1, at=10^c(0,2,4,6), labels=expression(1, 10^2, 10^4, 10^6)) axis(2, at=10^c(0,2,4,6), labels=expression(1, 10^2, 10^4, 10^6)) box() Duncan Murdoch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] 10^x instead 10EX on plot axes. How?
You can always draw the axes by hand. ?par with axes =FALSE ?axis ?plotmath for mathematical notation in R (for exponents) -- Bert Gunter Genentech Non-Clinical Statistics South San Francisco, CA The business of the statistician is to catalyze the scientific learning process. - George E. P. Box -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: Monday, July 10, 2006 9:43 AM To: r-help@stat.math.ethz.ch Subject: [R] 10^x instead 10EX on plot axes. How? Hi, I'm drawing a very simple plot with both axes logarithmic (default base 10). Example: vec=c(1,10,100,1000,1,10,100,1000) plot(vec,vec,log=xy) The axes on the plot now show the technical notation like 1E+3 but I would prefer to have it the notation 10 ^3 i.e. with the exponent here 3 superscript (raised). Any help very much appreciated! Best Regards Tom -- Feel free - 10 GB Mailbox, 100 FreeSMS/Monat ... __ R-help@stat.math.ethz.ch 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] Setting the colors of lines in a trellis plot...
With some help from those with expertise on this list, I managed to produce a plot using trellis that looked like I wanted it to look. Now, I need to take the same plot and make the lines on it color, but I want to specify the color for the lines myself. I've managed to make the key use the colors I want. I've managed to make the symbols of the actual plot use the colors I want. But I have been unable to find the correct incantation to make the lines of the actual plot use the colors I want. Here's the relevant section of code: mycolors - c(black, darkgreen, red) mylines - Rows(superpose.line, 1:numlines); mylines$col - mycolors mysymbols - Rows(superpose.symbol, 1:numlines); mysymbols$pch - c(15:18)[1:numlines] mysymbols$col - mycolors print(xyplot( panel = panel.superpose, log10(states) ~ size, groups=category, data=data, type='b', lwd = 2, par.settings = list(superpose.symbol=mysymbols), ylim=c(y_min, y_max), scales = list(tck=c(1, 0), axs=r, x=list(tick.number=(xmax - xmin + 1), at=xmin:xmax, labels=xmin:xmax, cex=1.75), y=list(axs=r, rot=c(90, 0), labels=y_labels, at=y_at, cex=1.75 ) ), key = list ( text = list(levels(data$category)), lines = list(type=b, lty=mylines$lty, pch=mysymbols$pch, cex=rep(1.25, numlines), col=mylines$col), x = .98, y = .25, corner=c(1,0) ), xlab = list(label=System Size, cex=2), ylab = list(label=States, cex=2), )) Can anyone help me out with this? Thanks! Jamie __ R-help@stat.math.ethz.ch 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] Setting the colors of lines in a trellis plot...
Use col= to specify colors, e.g. library(lattice) x - 1:12 xyplot(x ~ x, group = gl(3,4), col = 1:3, type = l, auto.key = TRUE) If this is not sufficiently close to your problem 1. cut your example down to a *minimal* size and 2. provide it as *self contained* and 3. *reproducible* code. On 7/10/06, Jamieson Cobleigh [EMAIL PROTECTED] wrote: With some help from those with expertise on this list, I managed to produce a plot using trellis that looked like I wanted it to look. Now, I need to take the same plot and make the lines on it color, but I want to specify the color for the lines myself. I've managed to make the key use the colors I want. I've managed to make the symbols of the actual plot use the colors I want. But I have been unable to find the correct incantation to make the lines of the actual plot use the colors I want. Here's the relevant section of code: mycolors - c(black, darkgreen, red) mylines - Rows(superpose.line, 1:numlines); mylines$col - mycolors mysymbols - Rows(superpose.symbol, 1:numlines); mysymbols$pch - c(15:18)[1:numlines] mysymbols$col - mycolors print(xyplot( panel = panel.superpose, log10(states) ~ size, groups=category, data=data, type='b', lwd = 2, par.settings = list(superpose.symbol=mysymbols), ylim=c(y_min, y_max), scales = list(tck=c(1, 0), axs=r, x=list(tick.number=(xmax - xmin + 1), at=xmin:xmax, labels=xmin:xmax, cex=1.75), y=list(axs=r, rot=c(90, 0), labels=y_labels, at=y_at, cex=1.75 ) ), key = list ( text = list(levels(data$category)), lines = list(type=b, lty=mylines$lty, pch=mysymbols$pch, cex=rep(1.25, numlines), col=mylines$col), x = .98, y = .25, corner=c(1,0) ), xlab = list(label=System Size, cex=2), ylab = list(label=States, cex=2), )) Can anyone help me out with this? Thanks! Jamie __ R-help@stat.math.ethz.ch 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] Availability of quadplot3d package (UseR!2006 Four Dimensional Barycentric Plots in 3D)
Carlos and Uwe: I apologize to all, but my server has crashed and I need to order a new computer! I'm at home now, but when I can get to my office I'll try to send each of you a copy. I will be submitting my packages to CRAN soon. Uwe, perhaps the quad3d stuff should just be incorporated into klaR? Geoff Geoffrey Matthews [EMAIL PROTECTED] -Original Message- From: Uwe Ligges [mailto:[EMAIL PROTECTED] Sent: Sat 7/8/2006 8:01 AM To: Carlos Ortega Cc: R-help@stat.math.ethz.ch; Geoffrey Matthews Subject: Re: [R] Availability of quadplot3d package (UseR!2006 Four Dimensional Barycentric Plots in 3D) Carlos Ortega wrote: Dear List, I have been unable fo find the package quadplot3d referred in the Abstract/Presentation Four Dimensional Barycentric Plot in 3D presented in the UserR!2006. Does anyone know if it is available ? And if so, if it is ported to Windows ? I think we should ask the author of that presentation, Geoffrey Matthews (CCing). I'm also interested in adding a corresponding plot feature into package klaR, which currently only covers static quadplots. Hence it would be nice to see the package quadplot3d on CRAN, or at least I'd like to see its functionality contributed to some other CRAN package such as misc3d. Uwe Ligges Thanks in anticipation, Carlos Ortega [[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
Re: [R] Setting the colors of lines in a trellis plot...
The col command worked. Thanks! Jamie On 7/10/06, Gabor Grothendieck [EMAIL PROTECTED] wrote: Use col= to specify colors, e.g. library(lattice) x - 1:12 xyplot(x ~ x, group = gl(3,4), col = 1:3, type = l, auto.key = TRUE) If this is not sufficiently close to your problem 1. cut your example down to a *minimal* size and 2. provide it as *self contained* and 3. *reproducible* code. On 7/10/06, Jamieson Cobleigh [EMAIL PROTECTED] wrote: With some help from those with expertise on this list, I managed to produce a plot using trellis that looked like I wanted it to look. Now, I need to take the same plot and make the lines on it color, but I want to specify the color for the lines myself. I've managed to make the key use the colors I want. I've managed to make the symbols of the actual plot use the colors I want. But I have been unable to find the correct incantation to make the lines of the actual plot use the colors I want. Here's the relevant section of code: mycolors - c(black, darkgreen, red) mylines - Rows(superpose.line, 1:numlines); mylines$col - mycolors mysymbols - Rows(superpose.symbol, 1:numlines); mysymbols$pch - c(15:18)[1:numlines] mysymbols$col - mycolors print(xyplot( panel = panel.superpose, log10(states) ~ size, groups=category, data=data, type='b', lwd = 2, par.settings = list(superpose.symbol=mysymbols), ylim=c(y_min, y_max), scales = list(tck=c(1, 0), axs=r, x=list(tick.number=(xmax - xmin + 1), at=xmin:xmax, labels=xmin:xmax, cex=1.75), y=list(axs=r, rot=c(90, 0), labels=y_labels, at=y_at, cex=1.75 ) ), key = list ( text = list(levels(data$category)), lines = list(type=b, lty=mylines$lty, pch=mysymbols$pch, cex=rep(1.25, numlines), col=mylines$col), x = .98, y = .25, corner=c(1,0) ), xlab = list(label=System Size, cex=2), ylab = list(label=States, cex=2), )) Can anyone help me out with this? Thanks! Jamie __ R-help@stat.math.ethz.ch 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] another tcl/tk query
Try looking at the tkwait.window function, it may be what you need. Create your tkwindow with all its components, then call tkwait.window on the toplevel window and the calling function will wait until that window goes away before continuing execution. Hope this helps, -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare [EMAIL PROTECTED] (801) 408-8111 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Charles Annis, P.E. Sent: Saturday, July 08, 2006 11:11 AM To: r-help@stat.math.ethz.ch Subject: [R] another tcl/tk query Greetings: I wish to use a tcl/tk widget to ask for user-selected parameter values. My widget works - it asks for and returns to my workspace the stuff I need. Here is a snippet of my code: ### OnOK - function() { LOG.X - as.logical(as.character(tclvalue(log.X.buttonValue))) LOG.Y - as.logical(as.character(tclvalue(log.Y.buttonValue))) natural.units.â.decision - as.double(as.character(tclvalue(â.decision))) natural.units.left.censor - as.double(as.character(tclvalue(left.censor))) natural.units.right.censor - as.double(as.character(tclvalue(right.censor))) tkdestroy(t2) } ### My problem is this: I would like to use the new input in the same routine that created, used, and destroyed the widget. I can't seem to do that. The routine executes with what it has. I must wait for the calling routine to end before I can use the new info, which is correctly place in the workspace, in subsequent R routines. Is there a way I can use the updated values in the same routine that created the widget? Thanks for your advice - and patience. Charles Annis, P.E. PS - I did read Prof. Ripley's post of Wed 8/31/2005 Re: [R] tcl/tk return problem but was unable to benefit from it. [EMAIL PROTECTED] phone: 561-352-9699 eFax: 614-455-3265 http://www.StatisticalEngineering.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 __ R-help@stat.math.ethz.ch 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] Setting the colors of lines in a trellis plot...
On 7/10/06, Jamieson Cobleigh [EMAIL PROTECTED] wrote: With some help from those with expertise on this list, I managed to produce a plot using trellis that looked like I wanted it to look. Now, I need to take the same plot and make the lines on it color, but I want to specify the color for the lines myself. I've managed to make the key use the colors I want. I've managed to make the symbols of the actual plot use the colors I want. But I have been unable to find the correct incantation to make the lines of the actual plot use the colors I want. Here's the relevant section of code: mycolors - c(black, darkgreen, red) mylines - Rows(superpose.line, 1:numlines); mylines$col - mycolors mysymbols - Rows(superpose.symbol, 1:numlines); mysymbols$pch - c(15:18)[1:numlines] mysymbols$col - mycolors print(xyplot( panel = panel.superpose, log10(states) ~ size, groups=category, data=data, type='b', lwd = 2, par.settings = list(superpose.symbol=mysymbols), You have forgotten to add 'superpose.line=mylines'. ylim=c(y_min, y_max), scales = list(tck=c(1, 0), axs=r, x=list(tick.number=(xmax - xmin + 1), at=xmin:xmax, labels=xmin:xmax, cex=1.75), y=list(axs=r, rot=c(90, 0), labels=y_labels, at=y_at, cex=1.75 ) ), key = list ( text = list(levels(data$category)), lines = list(type=b, lty=mylines$lty, pch=mysymbols$pch, cex=rep(1.25, numlines), col=mylines$col), You might consider using the 'auto.key' argument instead of 'key' (YMMV). -Deepayan x = .98, y = .25, corner=c(1,0) ), xlab = list(label=System Size, cex=2), ylab = list(label=States, cex=2), )) Can anyone help me out with this? Thanks! Jamie __ R-help@stat.math.ethz.ch 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] Adding Lines to Plot
There are a couple of different options that you can try (each generalize a bit differently to other situations). If you just want lines at the tick marks to go the entire width (and be generated automatically to match the labels), then look at ?par and look at 'tck', you can use this in the barplot call, or in a call to axis after the fact. If you want the bars on top then just issue another call to barplot with add=T. See abline and the h= option for drawing lines across the entire width. See par and the 'usr' entry for how to find the user coordinates of the edges of the plot (for adding lines by hand). See the cnvrt.coords function in the TeachingDemos package for more general conversions between coordinate systems (you can use this to draw a line from half way across to the edge, or other aritrary locations, probably overkill for your problem). The return value for barplot gives the center positions of the bars (the heights match the data you input to barplot), see the first example in the help for barplot. Hope this helps, -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare [EMAIL PROTECTED] (801) 408-8111 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of justin rapp Sent: Saturday, July 08, 2006 1:05 PM To: r-help@stat.math.ethz.ch Subject: [R] Adding Lines to Plot This seems like a question that I should be able to answer on my own but after looking at the documentation I cannot seem to find the correct method. How do I add lines to a bar plot that extend from the vertical axis? For example, my vertical axis is numbered in increments of 10 and I would like these to go across the whole graph. Also, is there a way to have R label the value of each bar so that I know the value of each factor that is being plotted? Thanks in advance. jdr __ R-help@stat.math.ethz.ch 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] Counting observations split by a factor when there are NAs in the data
I am a very novice R user, a social scientist (linguist) who is trying to learn to use R after being very familiar with SPSS. Please be kind! My concern: I cannot figure out a way to get an accurate count of observations of one column of data split by a factor when there are NAs in the data. I know how to use commands like tapply and summaryBy to obtain other summary statistics I am interested in, such as the following: tapply(RLWTEST, list(STATUS), mean, na.rm=T) summaryBy(RLWTEST~STATUS, data=lh.forgotten, FUN=c(mean, sd, min, max), na.rm=T) However, with tapply I know I cannot use length to get a count where there are NAs. summaryBy appears to work the same way. I do know how to get a count of the entire column using sum: sum(!is.na(lh.forgotten$RLWTEST)) However, this does not give me a count split up by my factor (STATUS). I have looked through Daalgard (2002) and Verzani (2005), and have searched the help files, but with no luck. Thank you in advance for your help. I love R and am interested in making it more accessible to social scientist types like me. I know it can do everything SPSS can and more, but sometimes the very simplest things seem to be a lot harder in R. Jenifer Dr. Jenifer Larson-Hall Assistant Professor of Linguistics University of North Texas (940)369-8950 __ R-help@stat.math.ethz.ch 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] Counting observations split by a factor when there are NA s in the data
Wouldn't something like table(status) give you want you want? E.g.: R status - factor(c(A, B, A, NA, A, B)) R table(status) status A B 3 2 Andy From: Jenifer Larson-Hall I am a very novice R user, a social scientist (linguist) who is trying to learn to use R after being very familiar with SPSS. Please be kind! My concern: I cannot figure out a way to get an accurate count of observations of one column of data split by a factor when there are NAs in the data. I know how to use commands like tapply and summaryBy to obtain other summary statistics I am interested in, such as the following: tapply(RLWTEST, list(STATUS), mean, na.rm=T) summaryBy(RLWTEST~STATUS, data=lh.forgotten, FUN=c(mean, sd, min, max), na.rm=T) However, with tapply I know I cannot use length to get a count where there are NAs. summaryBy appears to work the same way. I do know how to get a count of the entire column using sum: sum(!is.na(lh.forgotten$RLWTEST)) However, this does not give me a count split up by my factor (STATUS). I have looked through Daalgard (2002) and Verzani (2005), and have searched the help files, but with no luck. Thank you in advance for your help. I love R and am interested in making it more accessible to social scientist types like me. I know it can do everything SPSS can and more, but sometimes the very simplest things seem to be a lot harder in R. Jenifer Dr. Jenifer Larson-Hall Assistant Professor of Linguistics University of North Texas (940)369-8950 __ R-help@stat.math.ethz.ch 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] Start Model for POLYCLASS
Dear Dr. Graves, Thanks for your information and kindness! I really appreciate it. According to my findings so far, POLYCLASS doesnot seem to allow for variables being forced into the model. Even the commercialized MARS doesnot have this function either. However, one ad hoc approach that they suggested is to fit linear (or logistic) model and form rediduals first and then run MARS with the residuals. Best regards, -XG Su Xiaogang Su, Assistant Professor Department of Statistics and Actuarial Science University of Central Florida Orlando, FL 32816 (407) 823-2940 [O] [EMAIL PROTECTED] http://pegasus.cc.ucf.edu/~xsu/ Spencer Graves [EMAIL PROTECTED] 7/7/2006 5:49:05 AM I have not seen a reply, so I will offer a few suggestions, even though I've never used the 'polspline' package. I scanned several of the help pages and looked in ~\library\polspline' where R is installed on my hard drive and found no further documentation, and I found nothing new from RSiteSearch(polyclass). Googling for 'polyclass' led me to http://bear.fhcrc.org/~clk/soft.html;, the home page for Charles Kooperberg, the author and maintainer. If it were my problem, I might try writing him directly at the email address given in help(package=polyclass); I'm including him as a cc on this reply. I spent time looking at this, because 'polyclass', 'polymars' and 'polspline' seem potentially related related to one of my secondary interests. Unfortunately, I couldn't find sufficient documentation to allow me to proceed with the time I felt I could afford to invest in this right now. If it were my problem, before I wrote another email about this, I'd first list the function 'polyclass', make a local copy, then work through an example line by line using 'debug(polyclass)'. This 'debug' facility is remarkably powerful and easy to use. I've solved many problems like this using 'debug' in this way. If that failed to provide the necessary enlightenment, I'd submit another post including a self-contained example based on a modification of the 'iris' example featured in the 'polyclass' help page. Your example is NOT self contained, so I would NOT use that. Using the 'iris' example would make it much easier to explain clearly what you want. It also makes it much easier for someone like me to experiment with alternatives and describe what I did in terms of a tested example. Hope this helps. p.s. I suggest you also review the posting guide! www.R-project.org/posting-guide.html. Xiaogang Su wrote: Dear all, I have a question on how to set up the starting model in POLYCLASS and make sure the terms in the starting model retained in the final POLYCLASS model. In the function POLYMARS, this can be done using the STARTMODEL option. See below for example, I started with model y= b0 + b1*X1 + b2*X2 + b3*X4 + b4*X5 + b5*X2*X5 + e m00 - matrix(c( 1, NA, 0, NA, 1, 2, NA, 0, NA, 1, 4, NA, 0, NA, 1, 5, NA, 0, NA, 1, 2, NA, 5, NA, 1),nrow = 5, ncol=5, byrow=TRUE); m2 - polymars(response=PID2$y, predictors=PID2[,1:7], startmodel=m00) summary(m2) But I could not figure out how this works for POLYCLASS. There is an option FIT in POLYCLASS, which needs to be a POLYCLASS object though. Any suggestion or information is greatly appreciated. Sincerely, Xiaogang Su Xiaogang Su, Assistant Professor Department of Statistics and Actuarial Science University of Central Florida Orlando, FL 32816 (407) 823-2940 [O] [EMAIL PROTECTED] http://pegasus.cc.ucf.edu/~xsu/ __ R-help@stat.math.ethz.ch 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] Counting observations split by a factor when there are NAs in the data
Jenifer Larson-Hall [EMAIL PROTECTED] writes: I am a very novice R user, a social scientist (linguist) who is trying to learn to use R after being very familiar with SPSS. Please be kind! My concern: I cannot figure out a way to get an accurate count of observations of one column of data split by a factor when there are NAs in the data. I know how to use commands like tapply and summaryBy to obtain other summary statistics I am interested in, such as the following: tapply(RLWTEST, list(STATUS), mean, na.rm=T) summaryBy(RLWTEST~STATUS, data=lh.forgotten, FUN=c(mean, sd, min, max), na.rm=T) However, with tapply I know I cannot use length to get a count where there are NAs. summaryBy appears to work the same way. I do know how to get a count of the entire column using sum: sum(!is.na(lh.forgotten$RLWTEST)) However, this does not give me a count split up by my factor (STATUS). I have looked through Daalgard (2002) and Verzani (2005), and have Ahem! searched the help files, but with no luck. How about with(lh.forgotten, tapply(!is.na(RLWTEST), STATUS, sum) ) or maybe just table(STATUS[!is.na(RLWTEST)]) Thank you in advance for your help. I love R and am interested in making it more accessible to social scientist types like me. I know it can do everything SPSS can and more, but sometimes the very simplest things seem to be a lot harder in R. -- 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
[R] Source code for R graphics devices
Hi Folks, I'm trying to locate the source code for a (typical) R graphics device, in order to study how it's done. The underlying reason is that I'm thinking of trying to create a graphics device for 'pic' (the diagram drawing component of [g]troff). I thought the xfig device would be a good place to look, since the format of an xfig file is similar in nature (though very different in detail) to 'pic' code. When, in R, I type xfig to see the R code, I get the line .Internal(XFig(file, old$paper, old$family, old$bg, old$fg, old$width, old$height, old$horizontal, old$pointsize, old$onefile, old$pagecentre)) so I tried to locate the code for the Internal function Xfig. I can't seem to find it! So where should I be looking? Also -- a more general question on this topic -- presumably any R graphics device is driven by a stream of raw graphics data in some presumably device-independent format, which it then translates. Where can I find this, and information on its structure? With thanks, and best wishes, Ted. PS Total newbie on this front! E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 10-Jul-06 Time: 20:52:46 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] print color
Dear R Users, Is it possible to make R print the largest item in each row of a matrix X with red font? Example: 1247 8431 ... Therefore 7 and 8 should be in red color. I would appreciate any suggestion Robert McFadden [[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] Start Model for POLYCLASS
That's correct. The option is not there. You can force it in a starting model - but it can still be deleted, At some stage Insightful implemented a commercial set of polynomial spline functions in S-Plus that does have that option. They got it as far as a beta-library [quoting from my own homepage http://bear.fhcrc.org/~clk]: Doug Clarkson at Insightful Corp. has implemented commercial versions of hare, heft, logspline and a variety of related methods, with lots of additional features and options. See http://www.insightful.com/downloads/libraries/default.asp. Look for S+BEST. There are remaining bugs in this code Doug Clarkson is no longer at Insightful, and the code seems to be pretty much abbandoned. Charles Quoting Xiaogang Su [EMAIL PROTECTED]: Dear Dr. Graves, Thanks for your information and kindness! I really appreciate it. According to my findings so far, POLYCLASS doesnot seem to allow for variables being forced into the model. Even the commercialized MARS doesnot have this function either. However, one ad hoc approach that they suggested is to fit linear (or logistic) model and form rediduals first and then run MARS with the residuals. Best regards, -XG Su Xiaogang Su, Assistant Professor Department of Statistics and Actuarial Science University of Central Florida Orlando, FL 32816 (407) 823-2940 [O] [EMAIL PROTECTED] http://pegasus.cc.ucf.edu/~xsu/ Spencer Graves [EMAIL PROTECTED] 7/7/2006 5:49:05 AM I have not seen a reply, so I will offer a few suggestions, even though I've never used the 'polspline' package. I scanned several of the help pages and looked in ~\library\polspline' where R is installed on my hard drive and found no further documentation, and I found nothing new from RSiteSearch(polyclass). Googling for 'polyclass' led me to http://bear.fhcrc.org/~clk/soft.html;, the home page for Charles Kooperberg, the author and maintainer. If it were my problem, I might try writing him directly at the email address given in help(package=polyclass); I'm including him as a cc on this reply. I spent time looking at this, because 'polyclass', 'polymars' and 'polspline' seem potentially related related to one of my secondary interests. Unfortunately, I couldn't find sufficient documentation to allow me to proceed with the time I felt I could afford to invest in this right now. If it were my problem, before I wrote another email about this, I'd first list the function 'polyclass', make a local copy, then work through an example line by line using 'debug(polyclass)'. This 'debug' facility is remarkably powerful and easy to use. I've solved many problems like this using 'debug' in this way. If that failed to provide the necessary enlightenment, I'd submit another post including a self-contained example based on a modification of the 'iris' example featured in the 'polyclass' help page. Your example is NOT self contained, so I would NOT use that. Using the 'iris' example would make it much easier to explain clearly what you want. It also makes it much easier for someone like me to experiment with alternatives and describe what I did in terms of a tested example. Hope this helps. p.s. I suggest you also review the posting guide! www.R-project.org/posting-guide.html. Xiaogang Su wrote: Dear all, I have a question on how to set up the starting model in POLYCLASS and make sure the terms in the starting model retained in the final POLYCLASS model. In the function POLYMARS, this can be done using the STARTMODEL option. See below for example, I started with model y= b0 + b1*X1 + b2*X2 + b3*X4 + b4*X5 + b5*X2*X5 + e m00 - matrix(c( 1, NA, 0, NA, 1, 2, NA, 0, NA, 1, 4, NA, 0, NA, 1, 5, NA, 0, NA, 1, 2, NA, 5, NA, 1),nrow = 5, ncol=5, byrow=TRUE); m2 - polymars(response=PID2$y, predictors=PID2[,1:7], startmodel=m00) summary(m2) But I could not figure out how this works for POLYCLASS. There is an option FIT in POLYCLASS, which needs to be a POLYCLASS object though. Any suggestion or information is greatly appreciated. Sincerely, Xiaogang Su Xiaogang Su, Assistant Professor Department of Statistics and Actuarial Science University of Central Florida Orlando, FL 32816 (407) 823-2940 [O] [EMAIL PROTECTED] http://pegasus.cc.ucf.edu/~xsu/ __ R-help@stat.math.ethz.ch 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
[R] R newbie
Hello, I am new to R and still feeling my way thru it. I am trying to plot the values from this file below on the X-axis of a plot. I have attached the graph to the email...the one i am trying to recreate. Exonstart end 5'UTR 2254006022540121 1 2254012222540140 2 2254030322540493 3 2254155222541565 4 2254237322542519 5 2254426522544432 3'UTR 2254443322544856 I would like to create small rectangles on the x-axis as colored boxes from start position to end position of each exon...with the label showing 1, 2 etc under each respective exon. The 5' and 3' UTR would be colored different from the exons 1-5. Any suggestions and ideas would greatly appreciated. Thank you Kiran - This email is intended only for the use of the individual or entity to which it is addressed and may contain information that is privileged and confidential. If the reader of this email message is not the intended recipient, you are hereby notified that any dissemination, distribution, or copying of this communication is prohibited. If you have received this email in error, please notify the sender and destroy/delete all copies of the transmittal. Thank you. - __ R-help@stat.math.ethz.ch 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] Source code for R graphics devices
On 7/10/2006 3:52 PM, (Ted Harding) wrote: Hi Folks, I'm trying to locate the source code for a (typical) R graphics device, in order to study how it's done. The underlying reason is that I'm thinking of trying to create a graphics device for 'pic' (the diagram drawing component of [g]troff). I thought the xfig device would be a good place to look, since the format of an xfig file is similar in nature (though very different in detail) to 'pic' code. When, in R, I type xfig to see the R code, I get the line .Internal(XFig(file, old$paper, old$family, old$bg, old$fg, old$width, old$height, old$horizontal, old$pointsize, old$onefile, old$pagecentre)) so I tried to locate the code for the Internal function Xfig. I can't seem to find it! So where should I be looking? The devices are in the grDevices package. Look in the src directory there for the C code. XFig is in the devPS.c file. Duncan Murdoch Also -- a more general question on this topic -- presumably any R graphics device is driven by a stream of raw graphics data in some presumably device-independent format, which it then translates. Where can I find this, and information on its structure? With thanks, and best wishes, Ted. PS Total newbie on this front! Then it might be worth mentioning that this code isn't distributed with a binary version of R, you need a source tarball, and need to look in src/library/grDevices/src. Duncan Murdoch E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 10-Jul-06 Time: 20:52:46 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch 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 with package sem
Dear Gang, I'm not clear about the form of the model that you're trying to fit, nor what the form of the likelihood would be for such as model, but the short answer is that sem() requires that each covariance among observed variables be for the same observations and hence for the same number (N) of observations. BTW, I didn't respond to your query to r-help because I was out of town and unable to read email. Regards, John On Mon, 10 Jul 2006 11:36:06 -0400 Gang Chen [EMAIL PROTECTED] wrote: Hi Dr. Fox, I posted the following message to the R e-mail list, but didn't get any response. I was wondering whether you could shed some light on my problem. Thank you very much. Gang === From: [EMAIL PROTECTED] Subject: Re: [R] Package sem Date: July 9, 2006 5:49:41 PM EDT To: [EMAIL PROTECTED] Hi, I am trying to run some path analysis with Dr. Fox's sem package. The number N in sem(ram, S, N) is supposed to be the total number of observations, right? However, in my situation the effective number of degrees of freedom for each observed variable is estimated by some auto-regression process, thus each variable bears a different DF. Is there a way I could run such a path analysis with a vector of DF's in sem? Thanks, Gang John Fox Department of Sociology McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox/ __ R-help@stat.math.ethz.ch 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] Source code for R graphics devices
On 10-Jul-06 Ted Harding wrote: Hi Folks, I'm trying to locate the source code for a (typical) R graphics device, in order to study how it's done. [...] When, in R, I type xfig to see the R code, I get the line .Internal(XFig(file, old$paper, old$family, old$bg, old$fg, old$width, old$height, old$horizontal, old$pointsize, old$onefile, old$pagecentre)) so I tried to locate the code for the Internal function Xfig. I can't seem to find it! Many thanks to Duncan and Chuck for help and advice. Which also brought home to me that my problem was that I can't read! Grepping for Xfig didn't work. Grepping for XFig would have worked! Thanks, and best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 10-Jul-06 Time: 21:49:15 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] print color
One option is library(R2HTML) ?HTML.cormat The thing you're after is traffic highlighting (via CSS or HTML tags). If HTML.cormat() doesn't do exactly what you want, modify the source code. (By the way, I haven't used R2HTML so far so maybe there's a more appropriate function.) -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Robert Mcfadden Sent: Monday, July 10, 2006 4:00 PM To: R-help@stat.math.ethz.ch Subject: [R] print color Dear R Users, Is it possible to make R print the largest item in each row of a matrix X with red font? Example: 1247 8431 ... Therefore 7 and 8 should be in red color. I would appreciate any suggestion Robert McFadden [[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
Re: [R] 10^x instead 10EX on plot axes. How?
[EMAIL PROTECTED] wrote: Hi, I'm drawing a very simple plot with both axes logarithmic (default base 10). Example: vec=c(1,10,100,1000,1,10,100,1000) plot(vec,vec,log=xy) The axes on the plot now show the technical notation like 1E+3 but I would prefer to have it the notation 10 ^3 i.e. with the exponent here 3 superscript (raised). Any help very much appreciated! Hi Tom, Have a look at axis.mult in the plotrix package. Jim __ R-help@stat.math.ethz.ch 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] move axis label text
dear R wizards: sorry, I am stumped. what is the parameter to just move the location of the x-label and y-lable (not the labels on the ticks of the x-axis and y-axis)? probably obvious, but not to me right now... regards, /iaw __ R-help@stat.math.ethz.ch 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 newbie
Does this do what you want: x - Exonstart end 5'UTR 2254006022540121 1 2254012222540140 2 2254030322540493 3 2254155222541565 4 2254237322542519 5 2254426522544432 3'UTR 2254443322544856 y - read.table(textConnection(x), quote='', header=TRUE) # create the plot area plot(range(y$start, y$end), c(1,7), type='n', xaxt='n', yaxt='n', ylab='', xlab='') # create the colored rectangles rect(y$start, 1, y$end, 2, col=c('red', rep('green', nrow(y)-2), 'red')) # put labels at midpoints of rectangles axis(1, at=(y$start + y$end)/2, labels=as.character(y$Exon)) On 7/10/06, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: Hello, I am new to R and still feeling my way thru it. I am trying to plot the values from this file below on the X-axis of a plot. I have attached the graph to the email...the one i am trying to recreate. Exonstart end 5'UTR 2254006022540121 1 2254012222540140 2 2254030322540493 3 2254155222541565 4 2254237322542519 5 2254426522544432 3'UTR 2254443322544856 I would like to create small rectangles on the x-axis as colored boxes from start position to end position of each exon...with the label showing 1, 2 etc under each respective exon. The 5' and 3' UTR would be colored different from the exons 1-5. Any suggestions and ideas would greatly appreciated. Thank you Kiran - This email is intended only for the use of the individual or entity to which it is addressed and may contain information that is privileged and confidential. If the reader of this email message is not the intended recipient, you are hereby notified that any dissemination, distribution, or copying of this communication is prohibited. If you have received this email in error, please notify the sender and destroy/delete all copies of the transmittal. Thank you. - __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Jim Holtman Cincinnati, OH +1 513 646 9390 (Cell) +1 513 247 0281 (Home) What is the problem you are trying to solve? [[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] file.choose but for folder/directory?
Hello there, This question is relative to WindowsXP, using R 2.2.1: I am looking for a function that allows a user to interactively choose a directory so I can use list.files to process all the files in that directory. I've looked at getwd, but this is not interactive. The functions file.choose and choose.files are the right idea, but these only permit selection of a file within a folder, not the folder itself. Thanks for any and all suggestions, =Randy= R. Zelick email: [EMAIL PROTECTED] Department of Biology voice: 503-725-3086 Portland State University fax: 503-725-3888 mailing: P.O. Box 751 Portland, OR 97207 shipping: 1719 SW 10th Ave, Room 246 Portland, OR 97201 __ R-help@stat.math.ethz.ch 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] test regression against given slope for reduced major axis regression (RMA)
Hi, for testing if the slope of experimental data differs from a given slope I'm using the function test_regression_against_slope (see below). I am now confronted with the problem that I have data which requires a modelII regression (also called reduced major axes regression (RMA) or geometric mean regression). For this I use the function modelII (see below). What would be a good way of adapting test_regression_against_slope for use with RMA regression? The question I am trying to answer is: Does the slope acquired from experimental data differ significantly from theoretical predictions? Any feedback on this is highly appreciated. And if you need more info do not hesitate to ask. Kind Regards Patrick *test_regression_against_slope* --8schnipp-8--- test_regression_against_slope - function(x,y,slope_2) { ### TEST_REGRESSION_AGAINST_SLOPE tests a measured regression against a ### given regression. ### ### INPUT: ### ### x and y: raw data ### slope_2: the slope you would like to test against (ie 1/3) ### ### OUTPUT: ### ### pvalue: the P-value... ### upperlimit95 and lowerlimit95 give the 95 percent confidence ### intervall (two-tailed). ### ### see Sokal and Rohlf, p. 465/471 n - length(x) mydf - n-2 ## least square fit: x2 - (x-mean(x))^2 y2 - (y-mean(y))^2 ## regression (pedestrian solution): xy - (x-mean(x))*(y-mean(y)) slope1 - sum(xy)/sum(x2) intercept_a - mean(y) - slope1 * mean(x) ## model data y_hat: y_hat - intercept_a + slope1 * x ## least squares of model data: y_hat2 - (y - y_hat)^2 s2yx - sum(y_hat2) / (n-2) sb - sqrt(s2yx/sum(x2)) ts - (slope1 - slope_2) / sb pvalue - 2*(pt(abs(ts), df, lower.tail=FALSE)) ## 0.95 for one-tailed 0.975 for two-tailed t-distribution with ## alpha-5%: tval - qt(.975, df=mydf) ts2 - tval*sb lowerlimit95 - slope1 - ts2 upperlimit95 - slope1 + ts2 list(pvalue = pvalue, lowerlimit95 = lowerlimit95, upperlimit95 = upperlimit95) } --8schnapp-8--- *modelII* --8schnipp-8--- modelII - function(XjArray,YjArray){ ### ### ### Purpose: ### ### Calculates MODEL II Regression paramaters. Also called reduced ### major axis regression (Prentice 1987) or geometric mean ### regression (Webb et al. 1981). ### ### Input: ### ### Two one dimensional arrays XjArray and YjArray containing the X ### and Y vectors. ### ### XjArray = [0 0.9 1.8 2.6 3.3 4.4 5.2 6.1 6.5 7.4] ### YjArray = [5.9 5.4 4.4 4.6 3.5 3.7 2.8 2.8 2.4 1.5] ### ### Output: ### ### A list with the following: ### ### sumXjYj As Double ### sumXj As Double, sumYj As Double ### sumXjSquared As Double, sumYjSquared As Double ### n As Long ### varXj, varYj ### output(7) ### ### sumXjYj - 0 sumXj - 0 sumYj - 0 n - 0 n - length(XjArray) sumXjSquared - 0 sumYjSquared - 0 covariancexy - 0 for(i in 1:n){ sumXjYj - sumXjYj + XjArray[i] * YjArray[i] sumXj - sumXj + XjArray[i] sumYj - sumYj + YjArray[i] sumXjSquared - sumXjSquared + XjArray[i]^2 sumYjSquared - sumYjSquared + YjArray[i]^2 } ## Mean of X and Y vectors meanyj - sumYj / n meanxj - sumXj / n ## Create covariance for(i in 1:n){ covariancexy - covariancexy + ((XjArray[i] - meanxj) * (YjArray[i] - meanyj)) } covariancexy - covariancexy / n ## get variance of X and Y (SD) varXj - (n * sumXjSquared-sumXj^2)/(n*(n - 1)) varYj - (n * sumYjSquared-sumYj^2)/(n*(n - 1)) sdxij - (sumXjSquared)-(sumXj^2/n) sdxik - (sumYjSquared)-(sumYj^2/n) ## make beta'sgn function to return sign with magnitude of 1 betacoeff - sign(covariancexy) * ((varYj^0.5) / (varXj^0.5)) ## 'make intercept Intercept - meanyj - meanxj * betacoeff ## Make R the pearson produce moment correlation coefficient if (varYj==0 | varXj==0){ corrCoeff - 0 }else{ corrCoeff - (sumXjYj - ((sumXj * sumYj) / n)) / ((sdxij * sdxik)^0.5) } ## Make sample variances of betacoefficient and intercept variancebeta - (varYj / varXj) * ((1 - (corrCoeff ^ 2)) / n) varianceintercept - (varYj / n) * (1 - corrCoeff) * (2 + ((meanxj ^ 2) * ((1 + corrCoeff) / varXj))) sdbeta - variancebeta^0.5 sdintercept - varianceintercept^0.5 list(betacoeff=betacoeff, # - Steigung Intercept=Intercept, sdbeta=sdbeta, # standard deviation sdintercept=sdintercept, meanxj=meanxj, meanyj=meanyj, corrCoeff=corrCoeff) # - pearson correlation koeffizient ; } --8schnapp-8--- -- Snoopy (on being house-trained with a rolled-up newspaper): It does tend however to give one a rather distorted view of the press! __
Re: [R] file.choose but for folder/directory?
choose.dir it is documented on the same page as choose.files __ R-help@stat.math.ethz.ch 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] file.choose but for folder/directory?
On 7/10/2006 8:51 PM, Randy Zelick wrote: Hello there, This question is relative to WindowsXP, using R 2.2.1: Time to upgrade. choose.dir (mentioned by Rich) was introduced in the next release. Duncan Murdoch I am looking for a function that allows a user to interactively choose a directory so I can use list.files to process all the files in that directory. I've looked at getwd, but this is not interactive. The functions file.choose and choose.files are the right idea, but these only permit selection of a file within a folder, not the folder itself. Thanks for any and all suggestions, =Randy= R. Zelick email: [EMAIL PROTECTED] Department of Biology voice: 503-725-3086 Portland State University fax: 503-725-3888 mailing: P.O. Box 751 Portland, OR 97207 shipping: 1719 SW 10th Ave, Room 246 Portland, OR 97201 __ R-help@stat.math.ethz.ch 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] --no-save and --save toggle from inside R? + BATCH stderr
dear R wizards: is it possible to instruct R to save or no-save from inside R? or does this have to be given at invokation on the command-line?The same question applies to --no-restore-data, although this presumably would have to be decided in a .First() function or something like it. on a similar note, I would love a CMD BATCH invokation to output just one line to stderr at the end of its run which tells me whether the batch job ended prematurely because of an error. I know this can be undesirable if run in another program that refuses output (e.g., a web browser, where this can cause havoc if issued before the content-type:), so it would require that some switch like --silent would suppress this one line summary/error. This is not urgent---I can use the R return code to write a wrapper for this facility, but this might be nice default behavior. My guess is that most of us won't mind to see a flag in writing if something went badly wrong. sincerely, /iaw __ R-help@stat.math.ethz.ch 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] file.choose but for folder/directory -- THANKS!
Thanks to all who responded. I upgraded to 2.3.1 straight away, and choose.dir does exactly what I need. Cheers, =Randy= R. Zelick email: [EMAIL PROTECTED] Department of Biology voice: 503-725-3086 Portland State University fax: 503-725-3888 mailing: P.O. Box 751 Portland, OR 97207 shipping: 1719 SW 10th Ave, Room 246 Portland, OR 97201 __ R-help@stat.math.ethz.ch 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] --no-save and --save toggle from inside R? + BATCH stderr
ivo welch wrote: dear R wizards: is it possible to instruct R to save or no-save from inside R? You can specify the action when you call q() or quit(), but you can't change the default action. or does this have to be given at invokation on the command-line?The same question applies to --no-restore-data, although this presumably would have to be decided in a .First() function or something like it. This one can't currently be changed. Duncan Murdoch on a similar note, I would love a CMD BATCH invokation to output just one line to stderr at the end of its run which tells me whether the batch job ended prematurely because of an error. I know this can be undesirable if run in another program that refuses output (e.g., a web browser, where this can cause havoc if issued before the content-type:), so it would require that some switch like --silent would suppress this one line summary/error. This is not urgent---I can use the R return code to write a wrapper for this facility, but this might be nice default behavior. My guess is that most of us won't mind to see a flag in writing if something went badly wrong. sincerely, /iaw __ R-help@stat.math.ethz.ch 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] problem with install gstat package
Hi there, I tried to install gstat package to R in Linux. I follow the instruction to use command R CMD INSTALL -1 lib pkgs. But R shows syntax error in R CMD. Can someone help me with the installation of gstat package? Thank you! 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