Re: [R] nonlinear modeling with rational functions?
On Wed, Jun 16, 2004 at 10:28:39AM -0500, Spencer Graves wrote:
Rational functions (ratios of polynomials) often provide good
approximations to many functions. Does anyone know of any literature on
nonlinear modeling with rational functions, sequential estimation,
diagnostics, etc.? I know I can do it with nls and other nonlinear
regression functions, but I'm wondering what literature might exist
discussing how a search for an appropriate rational approximation?
The book
@Book{judd98,
author = {Judd, Kenneth L},
title ={Numerical methods in economics},
publisher ={MIT Press},
year = 1998
}
provides some information about constructing rational approximations
(aka Padé approximations) for functions with known algebraic forms,
but not for the estimation from data.
However, the text has some references, in particular
Cuyt, A and L Wuytack. 1986. Nonlinear Numerical Methods: Theory and
Practice. Amsterdam: North-Holland
which you might find helpful.
Best,
Tamas
--
Tamás K. Papp
E-mail: [EMAIL PROTECTED]
Please try to send only (latin-2) plain text, not HTML or other garbage.
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] PCA and AR on circular data
Could anyone advise the best method of performing PCA and/or AR modelling on circular data? Until now I have been working with vectors (speed*direction), but am now wanting to look purely at the directional aspect. Thanks in advance.. Laura __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Date Calculation
Hi, I've been playing with: joinDate - format(strptime(as.vector(forum[,2]), %d-%b-%y), +%d-%b-%Y) today - format(strptime(as.vector(14-Jun-04), %d-%b-%Y), + %d-%b-%Y) joinDate [1] 04-Feb-2004 13-Feb-2004 26-Feb-2004 27-Feb-2004 27-Feb-2004 [6] 27-Feb-2004 29-Feb-2004 01-Mar-2004 02-Mar-2004 07-Mar-2004 [11] 08-Mar-2004 17-Mar-2004 20-Mar-2004 22-Mar-2004 22-Mar-2004 [16] 23-Mar-2004 23-Mar-2004 24-Mar-2004 01-Apr-2004 01-Apr-2004 [21] 01-Apr-2004 01-Apr-2004 02-Apr-2004 06-Apr-2004 09-Apr-2004 [26] 11-Apr-2004 14-Apr-2004 03-May-2004 04-May-2004 30-May-2004 [31] 01-Jun-2004 10-Jun-2004 14-Jun-2004 17-Jun-2004 17-Jun-2004 today [1] 14-Jun-0004 joinDate - today Error in joinDate - today : non-numeric argument to binary operator But it didn't quite work. What I'd like joinDate - today to return is the number of days to today, since joinDate. I'm sure it has been asked before however a search on r-help didn't found me any relevant information *_*. Cheers, Kevin Ko-Kang Kevin Wang, MSc(Hon) SLC Stats Workshops Co-ordinator The University of Auckland New Zealand __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] subset(..., drop=TRUE) doesn't seem to work.
Hoi Peter,
--On woensdag 16 juni 2004 17:35 +0200 Peter Dalgaard
[EMAIL PROTECTED] wrote:
Anyways, the way out is
d2 - subset(dd,c==1)
ifac - sapply(dd,is.factor)
d2[ifac] - lapply(d2[ifac],factor)
or
d2 - subset(dd,c==1)
d2[] - lapply(d2, function(x) if (is.factor(x)) factor(x) else x)
My toolbox.r now contains:
my.subset - function(x, drop.unused.levels=FALSE, ...) {
subsetted - subset(x, ...)
if (drop.unused.levels) {
subsetted[] - lapply(subsetted,
function(x) if (is.factor(x)) factor(x) else x)
}
subsetted
}
Thnx for all the help!
--
Paul Lemmens
NICI, University of Nijmegen ASCII Ribbon Campaign /\
Montessorilaan 3 (B.01.05)Against HTML Mail \ /
NL-6525 HR Nijmegen X
The Netherlands / \
Phonenumber+31-24-3612648
Fax+31-24-3616066
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] nonlinear modeling with rational functions?
Tamas == Tamas Papp [EMAIL PROTECTED]
on Thu, 17 Jun 2004 08:40:13 +0200 writes:
Tamas On Wed, Jun 16, 2004 at 10:28:39AM -0500, Spencer Graves wrote:
Rational functions (ratios of polynomials) often provide good
approximations to many functions. Does anyone know of any literature on
nonlinear modeling with rational functions, sequential estimation,
diagnostics, etc.? I know I can do it with nls and other nonlinear
regression functions, but I'm wondering what literature might exist
discussing how a search for an appropriate rational approximation?
Tamas The book
Tamas @Book{judd98,
Tamas author = {Judd, Kenneth L},
Tamas title ={Numerical methods in economics},
Tamas publisher ={MIT Press},
Tamas year = 1998
Tamas }
Tamas provides some information about constructing rational approximations
Tamas (aka Padé approximations) for functions with known algebraic forms,
Tamas but not for the estimation from data.
Tamas However, the text has some references, in particular
Tamas Cuyt, A and L Wuytack. 1986. Nonlinear Numerical Methods: Theory and
Tamas Practice. Amsterdam: North-Holland
Tamas which you might find helpful.
But be aware:
In these contexts, one is interested in L_inf() norm
approximations, not in L_2 or even L_1 as we are in statistics.
L_inf aka sup-norm aka minimax or worst case approximation
makes much sense in the context of numerical approximation:
e.g. you might want a Padé approximation of the sin() function with
a maximal (absolute) error of 1e-12 to be implemented for
your pocket calculator.
OTOH, the sup-norm is quite a bad idea for data fitting, since
there, errors are typically either normal or more heavy tailed,
i.e. if using an L_p norm, it should be p = 2.
Hoping this helps even more,
Martin Maechler
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Date Calculation
KKWa == Ko-Kang Kevin Wang [EMAIL PROTECTED] on Thu, 17 Jun 2004 19:57:39 +1200 writes: KKWa Hi, KKWa I've been playing with: joinDate - format(strptime(as.vector(forum[,2]), %d-%b-%y), KKWa +%d-%b-%Y) today - format(strptime(as.vector(14-Jun-04), %d-%b-%Y), KKWa + %d-%b-%Y) joinDate KKWa [1] 04-Feb-2004 13-Feb-2004 26-Feb-2004 27-Feb-2004 KKWa 27-Feb-2004 KKWa [6] 27-Feb-2004 29-Feb-2004 01-Mar-2004 02-Mar-2004 KKWa 07-Mar-2004 KKWa today KKWa [1] 14-Jun-0004 joinDate - today KKWa Error in joinDate - today : non-numeric argument to binary operator KKWa But it didn't quite work. What I'd like joinDate - KKWa today to return is the number of days to today, since KKWa joinDate. I'm sure it has been asked before however a KKWa search on r-help didn't found me any relevant KKWa information *_*. Kevin, [[did you have tough day? usually your Q/A are much better ;-()]] both joinDate and today are results of format() and hence character vectors (with no time information left). You didn't really expect - to work with characters? OTOH, - properly does work with POSIX*t objects, see, e.g., ?difftime Regards, Martin __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Date Calculation
You have formatted the dates, so you are trying to subtract character strings. Convert to Date instead (although you can subtract POSIXlt dates). On Thu, 17 Jun 2004, Ko-Kang Kevin Wang wrote: Hi, I've been playing with: joinDate - format(strptime(as.vector(forum[,2]), %d-%b-%y), +%d-%b-%Y) today - format(strptime(as.vector(14-Jun-04), %d-%b-%Y), + %d-%b-%Y) joinDate [1] 04-Feb-2004 13-Feb-2004 26-Feb-2004 27-Feb-2004 27-Feb-2004 [6] 27-Feb-2004 29-Feb-2004 01-Mar-2004 02-Mar-2004 07-Mar-2004 [11] 08-Mar-2004 17-Mar-2004 20-Mar-2004 22-Mar-2004 22-Mar-2004 [16] 23-Mar-2004 23-Mar-2004 24-Mar-2004 01-Apr-2004 01-Apr-2004 [21] 01-Apr-2004 01-Apr-2004 02-Apr-2004 06-Apr-2004 09-Apr-2004 [26] 11-Apr-2004 14-Apr-2004 03-May-2004 04-May-2004 30-May-2004 [31] 01-Jun-2004 10-Jun-2004 14-Jun-2004 17-Jun-2004 17-Jun-2004 today [1] 14-Jun-0004 joinDate - today Error in joinDate - today : non-numeric argument to binary operator But it didn't quite work. What I'd like joinDate - today to return is the number of days to today, since joinDate. I'm sure it has been asked before however a search on r-help didn't found me any relevant information *_*. Cheers, Kevin Ko-Kang Kevin Wang, MSc(Hon) SLC Stats Workshops Co-ordinator The University of Auckland New Zealand __ [EMAIL PROTECTED] mailing list https://www.stat.math.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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Tutorial for graphics
On Wed, Jun 16, 2004 at 10:12:49PM -0400, Phil wrote: I'm coming to R from Matlab and I'm finding it difficult to find a good introduction to graphics in R. (The best I've found so far is Ch. 4 of Emmanuel Paradis R for Beginners. Still, I have been unable to discover simple things like how to resize the axes on an existing plot I'm afraid it's not possible to resize existing axis, by mouse. But if you want to add the another plot with differnet scale, you could do it this way: plot(...)# you normally plot 1st graph par(new = TRUE) # 2nd graph won't clean the frame plot(..., axes = FALSE, xlab = , ylab = ) # plotting 2ng graph axis(4) how to add (or change) axis labels on an existing plot, etc. You can plot without axes plot(..., axes = FALSE), then add axes, with labels or not, ticks with any length, and so on. Can anyone point me to a suitable tutorial, or even tell me how to perform those tasks? help(plot.default) help(par) Also, Matlab's graphical widget has the ability to zoom (and unzoom) by drawing a rectangle on the graph with the mouse. Is there anything similar in R? Probably there are R packages perfoming such a thing, but I do not use R interactively at all. I change xlim() and/or ylim() in R script, and source() it again :) -- WBR, Timur __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Date Calculation
Hi, -Original Message- From: Martin Maechler [mailto:[EMAIL PROTECTED] Kevin, [[did you have tough day? usually your Q/A are much better ;-()]] Thanks to those who have replied, and yes shame on me.. [I also realised I can just use Sys.Date() to get today's date, instead of typing it in..I really had a tough day *_*] Cheers, Kevin __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] packages data-sets
It's possible to create a package with functions and data, from which the use library(pkg-name) attaches not only the functions, but also the data? I want avoid to use data(dataset, package=name) because this makes a global copy of the data-set ... Anyone could help me? Meinhard __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] packages data-sets
On Thu, 17 Jun 2004, Meinhard Ploner wrote: It's possible to create a package with functions and data, from which the use library(pkg-name) attaches not only the functions, but also the data? I want avoid to use data(dataset, package=name) because this makes a global copy of the data-set ... What do you mean by `global'? You cannot have access to R data without having it in memory and having it visible. You don't need to have it in the user workspace (.GlobalEnv), though. Suggestions: 1) See how MASS does it. Data objects are loaded into the MASS namespace on first use. 2) See ?attach and attach an R save file containing your datasets -- wasteful unless you use them all at once. It is planned that data() will be superseded by a better mechanism in R 2.0.0 (but that plan has slipped a bit already, which is why MASS does it differently). -- 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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Question : simtest result
Dear Mr.Torsten: At Wed, 16 Jun 2004 08:04:08 +0200 (CEST), Torsten Hothorn wrote: Call: simtest.formula(formula = y4 ~ f2, data = dat2, type = Dunnett) Dunnett contrasts for factor f2 Contrast matrix: f21 f22 f23 f24 f25 f22-f21 0 -1 1 0 0 0 f23-f21 0 -1 0 1 0 0 f24-f21 0 -1 0 0 1 0 f25-f21 0 -1 0 0 0 1 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj f25-f215.167 -4.6441.022 0.000 0.000 0.000 f23-f212.875 -2.8131.022 0.008 0.024 0.022 f24-f212.625 -2.5691.022 0.015 0.029 0.028 f22-f212.125 -2.0791.113 0.045 0.045 0.045 - Chihiro, Frank and I used your data to check the program and example with an independent algorithm and implementation (Westfall-Young stepdown resampling procedure). Theory suggests that the results should be similar to (but not necessarily the same as) those obtained with multcomp in this special case. These are the adjusted p-values obtained with the Westfall-Young approach for 100,000 replications: which fit nicely with the ones obtained from multcomp. Hope this helps sorry for the delay, Thank you. At Mon, 07 Jun 2004 12:14:26 +0900, myself-oak wrote: dunnett(dat2,1,2) rho=0.426 group:5 t=4.644 p=0.000 group:3 t=2.813 p=0.028 group:4 t=2.569 p=0.051 --- (B) group:2 t=2.079 p=0.145 (sorted in order of p values.) p values are different although t values are equal. I got the following inequality from the appended chart of a book. hm, without knowing what 2.558 d(5, 35, 0.4263464, 0.05) 2.598 means it is hard to tell what the problem is. Could you please explain it further? The alternative hypothesis is two sided. When significant level is equal to 0.05 , number of groups=5, df of error=35 and rho=0.426, I think that absolute t-value should be between 2.558 and 2.598. So, (B) is easy to understand for me than (A). You said adj p values are ... 0.0437 0.0260 0.0204 0.0001 I used SPSS ver.10 and got the following result. Dunnett t (two sided) | --- | -- | | | -- | | | mean diff. | s.e. | adj p| 95% C.I. | | -- | -- | (I-J) | | | -- | - | | (I) V3 | (J) V3 || | || | | -- | -- | -- | | | -- | - | | 2 | 1 | 2.125 | 1.022| .145 | -.507 | 4.757 | | -- | -- | -- | | | -- | - | | 3 | 1 | 2.875(*) | 1.022| .028 | .243 | 5.507 | | -- | -- | -- | | | -- | - | | 4 | 1 | 2.625 | 1.022| .051 | -7.325E-03 | 5.257 | | -- | -- | -- | | | -- | - | | 5 | 1 | 5.167(*) | 1.113| .000 | 2.301 | 8.032 | | -- | -- | -- | | | -- | - | I might make a mistake in the way to use the simtest()...How should I think? Best regards, -- kuroki GnuPG fingerprint = 90FD FE79 905F 26F9 29C4 096F 8AA2 2C42 5130 1469 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] How to order a vector
Hi all I have a vector like this 2003 2002 2001 2000 1999 1998 1997 1996 106 105 106 106105 106 101 107 How can I get it sorted right(19962003)? Thank you Luis Ridao Cruz Fiskirannsóknarstovan Nóatún 1 P.O. Box 3051 FR-110 Tórshavn Faroe Islands Phone: +298 353900 Phone(direct): +298 353912 Mobile: +298 580800 Fax: +298 353901 E-mail: [EMAIL PROTECTED] Web:www.frs.fo __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to order a vector
[EMAIL PROTECTED] wrote on 17/06/2004 13:10:35: Hi all I have a vector like this 2003 2002 2001 2000 1999 1998 1997 1996 106 105 106 106105 106 101 107 How can I get it sorted right(19962003)? rev(1:4) [1] 4 3 2 1 HTH, Tobias __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to order a vector
Hi! I assume that 2003 2002 are the names of your vector myvector myvector[order(names(myvector))] Eryk *** REPLY SEPARATOR *** On 6/17/2004 at 12:10 PM Luis Rideau Cruz wrote: Hi all I have a vector like this 2003 2002 2001 2000 1999 1998 1997 1996 106 105 106 106105 106 101 107 How can I get it sorted right(19962003)? Thank you Luis Ridao Cruz Fiskirannsóknarstovan Nóatún 1 P.O. Box 3051 FR-110 Tórshavn Faroe Islands Phone: +298 353900 Phone(direct): +298 353912 Mobile: +298 580800 Fax: +298 353901 E-mail: [EMAIL PROTECTED] Web:www.frs.fo __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Dipl. bio-chem. Eryk Witold Wolski@MPI-Moleculare Genetic Ihnestrasse 63-73 14195 Berlin 'v' tel: 0049-30-83875219 / \ mail: [EMAIL PROTECTED]---W-Whttp://www.molgen.mpg.de/~wolski __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to order a vector
See ?order. In particular, if you do: a - order(names(vec)), a will contain the indices of the vector in order, so vec[a] will be in the correct order. Sean On 6/17/04 7:10 AM, Luis Rideau Cruz [EMAIL PROTECTED] wrote: Hi all I have a vector like this 2003 2002 2001 2000 1999 1998 1997 1996 106 105 106 106105 106 101 107 How can I get it sorted right(19962003)? Thank you Luis Ridao Cruz Fiskirannsóknarstovan Nóatún 1 P.O. Box 3051 FR-110 Tórshavn Faroe Islands Phone: +298 353900 Phone(direct): +298 353912 Mobile: +298 580800 Fax: +298 353901 E-mail: [EMAIL PROTECTED] Web:www.frs.fo __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Using predict.lm()
Greetings, Following the example in help(predict.lm): x - rnorm(15) y - x + rnorm(15) new - data.frame(x = seq(-3, 3, 0.5)) predict(lm(y ~ x), new) predicts the response elements corresponding to new$x as can be viewed by: plot(x,y) lines(new$x,predict(lm(y ~ x), new)) I am trying to extend this fitting and prediction over a variety of factors as follows: f-rep(c(FIRST,SECOND),each=15) f-as.factor(f) x-rep(rnorm(15),2) y-x+rnorm(length(x)) old-data.frame(f=f,x=x,y=y) new-data.frame(f=rep(levels(f),each=length(seq(-4,4,0.2))),x=seq(-4,4,0.2)) ...where variable new simply substitutes a differing domain than old. When I try to predict on the frame new using x y, I get a response that corresponds to the length of new: predict(lm(y~x),new) but when I use the same variables from within the frame old, the frame new is ignored: predict(lm(old$y~old$x),new) ...results in a response the length of old$x (presumably predicting over the values of old$x). Furthermore, this behavior also precludes using something more useful, i.e.: predict(lm(old$y~old$f/(1+old$x)-1),new) to return predictions over a number of factors over redefined domains. In my case, I am attempting to do 2nd order polynomial fitting over noisy data collected for a large number of factors (~85). The data were collected for each factor at convenient (and therefore dissimilar) points within a common domain, but I need to compare the responses of each factor at similar points within the common domain. I am obviously missing something here because I continue to be puzzled by the result. I had thought (perhaps erroneously) that lm() would return a model object that would permit prediction. Indeed: lm(old$y~old$f/(1+old$x)-1) ...results in: Call: lm(formula = old$y ~ old$f/(1 + old$x) - 1) Coefficients: old$fFIRSTold$fSECOND old$fFIRST:old$x old$fSECOND:old$x -0.08489 -0.058391.153510.72981 which clearly provides a model fit for each factor, and identifies the factor from which each model coefficient was extracted, so lm() does provide the capability to predict over the factors. It seems however (as nearly as I can tell), that predict simply ignores the frame new altogether, failing even to provide a warning. Is this the intended behavior? Have I missed something very simple or have a fundamental misunderstanding of how this should work? Lastly, I'd appreciate any suggestions that avoid the lengthy and wholly undesirable brute force approach I an now considering. Thanks Best Regards, Steve __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Using predict.lm()
On Thu, 17 Jun 2004, Steven White wrote: Following the example in help(predict.lm): x - rnorm(15) y - x + rnorm(15) new - data.frame(x = seq(-3, 3, 0.5)) predict(lm(y ~ x), new) predicts the response elements corresponding to new$x as can be viewed by: plot(x,y) lines(new$x,predict(lm(y ~ x), new)) Note that the model is fitted to `x' and new contains `x'. You haven't copied that. I am trying to extend this fitting and prediction over a variety of factors as follows: f-rep(c(FIRST,SECOND),each=15) f-as.factor(f) x-rep(rnorm(15),2) y-x+rnorm(length(x)) old-data.frame(f=f,x=x,y=y) new-data.frame(f=rep(levels(f),each=length(seq(-4,4,0.2))),x=seq(-4,4,0.2)) ...where variable new simply substitutes a differing domain than old. When I try to predict on the frame new using x y, I get a response that corresponds to the length of new: predict(lm(y~x),new) but when I use the same variables from within the frame old, That you have not done correctly: see ?lm. the frame new is ignored: No, it is not ignored but it does not contain a variable named `old$x' and your workspace does. newdata is the first place to look for variables, but not the only place. predict(lm(old$y~old$x),new) ...results in a response the length of old$x (presumably predicting over the values of old$x). Furthermore, this behavior also precludes using something more useful, i.e.: predict(lm(old$y~old$f/(1+old$x)-1),new) to return predictions over a number of factors over redefined domains. In my case, I am attempting to do 2nd order polynomial fitting over noisy data collected for a large number of factors (~85). The data were collected for each factor at convenient (and therefore dissimilar) points within a common domain, but I need to compare the responses of each factor at similar points within the common domain. I am obviously missing something here because I continue to be puzzled by the result. I had thought (perhaps erroneously) that lm() would return a model object that would permit prediction. Indeed it does. Indeed: lm(old$y~old$f/(1+old$x)-1) ...results in: Call: lm(formula = old$y ~ old$f/(1 + old$x) - 1) Coefficients: old$fFIRSTold$fSECOND old$fFIRST:old$x old$fSECOND:old$x -0.08489 -0.058391.153510.72981 which clearly provides a model fit for each factor, and identifies the factor from which each model coefficient was extracted, so lm() does provide the capability to predict over the factors. It seems however (as nearly as I can tell), that predict simply ignores the frame new altogether, failing even to provide a warning. Nope. You just haven't set new to match your fit. Is this the intended behavior? Have I missed something very simple or have a fundamental misunderstanding of how this should work? Yes, yes. You should be using lm(y ~ f/(1+x)-1, data=old) etc, although in your example you could omit data=old. That is in all good books on the S language -- 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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] non-linear binning? power-law in R
Why not try to avoid binning (and density plot) at all? An alternative could be a qqplot (as a log-log-plot), e.g. plot(ppoints(length(x4)), x4[order(x4)], log=xy) abline(lm(log(x4[order(x4)])~log(ppoints(length(x4, col=red) If the assumptions of uniform distribution and power transformation y=a*x**b are true, the coefficient of lm estimates the exponent b. Herwig -- Dr. Herwig Meschke Wissenschaftliche Beratung Hagsbucher Weg 27 D-89150 Laichingen phone +49 7333 210 417 / fax +49 7333 210 418 email [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Resolution of plots
On Wed, 16 Jun 2004, Prof Brian Ripley wrote: You will have to tell us more. Exporting how: to what format using what device and what exact command on what operating system? The only device I know of that even knows about dpi is bitmap() and that has no such limit unless imposed by your implementation of ghostscript. There is an issue with PNG. libpng provides png_set_pHYs to set resolution (in pixels/metre) but provides a default if it is not set. We don't set it, and so get the default resolution. The default is good for screen display, and irrelevant if you are including the file in a document, but photo editing software will show the true resolution, and publishers may care. One way to change the resoluton of PNG files is to use SNG (sng.sourceforge.net) to translate the PNG to a text format, edit the pHYs chunk and translate it back. -thomas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] 2D Kolmogorov-Smirnov test: solution
Hi - A little while ago I posted a question about the implementation of a two-dimensional analog of the Kolmogorov-Smirnov test in R[1][2]. As there isn't one, as far as I know, people might be interested in a very fast C++ implementation called MUAC which is available as a function and as a standalone program from http://www.acooke.org/jara/muac/index.html. Apparently the code is interesting from a computer science viewpoint, and an MPI-powered parallel version is available from http://beowulf.lcs.mit.edu/18.337-2002/projects-2002/ianchan/KS2D/ Project%20Page.htm If anyone has any experience of this code, and in particular words of warning, I'd love to hear about them (off-list would be best). Rich [1] Peacock J, Monthly Notices of the Royal Astronomical Society, 1983, vol 202 p615: Two-Dimensional Goodness-of-Fit Testing in Astronomy [2] Fasano G, Franceschini A. Monthly Notices of the Royal Astronomical Society, 1987 vol 225 p 155: A Multidimensional Version of the Kolmogorov-Smirnov Test Rich Grenyer, Ph.D. Biology Department - University of Virginia Gilmer Hall Charlottesville, Virginia VA 22904 United States of America tel: (+1) 434 982 5629 fax: (+1) 434 982 5626 http://faculty.virginia.edu/gittleman/rich __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Re: Clustering in R
Thanks a lot, Michael!
I cc to R-help, where this question really belongs {as the
'Subject' suggests itself...} -- please drop 'bioconductor' from
CC'ing further replies.
michael == michael watson (IAH-C) [EMAIL PROTECTED]
on Thu, 17 Jun 2004 09:16:59 +0100 writes:
michael OK, admittedly it is not incredibly simple, but it
michael is not *that* difficult.
michael If you are familiar with R, it should take you an
michael hour or two; if unfamiliar, perhaps a day or two.
michael The commands you want (and need to read the help on) are:
michael hclust
michael plclust
michael cutree
and I would add identify.hclust() {and rect.hclust()}
a very neat but not known / used enough function
a link to which I have just added to the help(hclust) page.
Look at its examples {not with example() since they are
dontrun} correcting the extraneous . in the last (and
coolest!) example!
michael dendrogram
michael as.dendrogram
michael heatmap
where you use dendrograms produced from hclust objects via
as.dendrogram(hc-obj) or also twins objects produced from
package cluster's agnes() or diana() via
as.dendrogram(as.hclust( twins-obj ) )
help(dendrogram) also mentions
[[ (and shows examples) and cut() for cutting dendrograms and shows
how you can depict dendrograms into its parts.
michael With intelligent use of hclust - cutree - subsetting - hclust
michael (in that order) you will be able to drill down
michael into your dendrogram and create sub-trees - until
michael you get to the level where you can see your gene
michael names.
or also
hclust - as.dendrogram - cut - ..
- [[ -
Note that there also is reorder.dendrogram() for reordering
dendrogram nodes ``sensibly'' --- something that heatmap() does,
but you can play with quite a bit.
Further, note Catherine Hurley's gclus package which
orders/reorders 'hclust' objects directly, but with a more
interesting algorithm.
Note that I'd strongly recommend to use R 1.9.1 beta for these,
since I know which bugs in the dendrogram code I have fixed
since R 1.9.0...
michael An important message to take home here is that if
michael you have 14000 genes and therefore 14000 labels,
michael it's going to be difficult to display your tree in
michael ANY software, including the expensive commercial products.
not showing the labels and using identify.hclust() and the
command line to extract the indices of observations in
clusters (and subclusters) and visualize them in other, non-dendrogram plots,
might well be feasible.
michael Let me know how you get on
michael Thanks
michael Mick
michael -Original Message-
michael From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
michael Sent: 16 June 2004 21:26
michael To: [EMAIL PROTECTED]
michael Subject: [BioC] Clustering in R
Dear list members,
I'm an undergrad and I work in a lab at Brandeis.
I am trying to cluster around 14,000 genes across 6
microarray experiments. Two of these experiments
are replicates. I have decided to use R since it
seems to be the most complete and flexible software
package for normalization and clustering of
microarray data.
The problem is that I am new to clustering and to
R. Just to mention of a few of the problems I'm
having: the dendrogram that is drawn by R from the
agnes object is far too dense to see any of the
gene names; kmeans won't work, returning an error
saying that my data has NAs in it (there weren't
any missing values in the original table though);
I'd like to be able to see a heatmap or a
cumulative plot of expression profiles for genes
that are clustered together or are on the same
branch of the dendrogram.
I know that these questions are probably very
simple, but I can't seem to find the answer to them
online or in the documentation. If anyone can
answer these questions or direct me toward
resources that deal with clustering in R or
BioConductor, a basic tutorial that takes a
practical approach to it, I would really appreciate
it. Any other reading material that isn't too
heavy on statistics that deals with clustering for
that matter, would be very helpful.
Thank you in advance,
Wayne Mak
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Error with arima()
Could someone please give a brief explanation, or pointer to an explanation, of the following error: arima(ts.growth, order = c(1,0,0),include.mean=T) Error in arima(ts.growth, order = c(1, 0, 0), include.mean = T) : non-stationary AR part from CSS and why it does not arise with arima0(ts.growth, order = c(1,0,0)) Many thanks Dr. Daniel P. Bebber Department of Plant Sciences University of Oxford South Parks Road Oxford OX1 3RB [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] disappointed with x-axes in hist and density plots
I've got a few issues with the x-axes in the histogram and density plots. First, often the default x-axis doesn't even extend to the length of my data. R often draws histogram bars (or density lines) farther than the drawn x-axis extends. For example, I might have a histogram bar at -15,000. But I wouldn't know that, because the most negative number on the x-axis is -10,000. The second issue is the use of scientific notation. Yes I can read it, but I don't prefer it. Is there any way for R just to print out 100 and not 1e+6 on these charts? Thanks for your help. Rishi [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] L é vy distribution?
Hi all, Does R support the Lévy distribution? I looked everywhere and couldn't find it. Thanks for any help! Gabriel __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] disappointed with x-axes in hist and density plots
On Thu, 17 Jun 2004, Rishi Ganti wrote: I've got a few issues with the x-axes in the histogram and density plots. First, often the default x-axis doesn't even extend to the length of my data. R often draws histogram bars (or density lines) farther than the drawn x-axis extends. For example, I might have a histogram bar at -15,000. But I wouldn't know that, because the most negative number on the x-axis is -10,000. The second issue is the use of scientific notation. Yes I can read it, but I don't prefer it. Is there any way for R just to print out 100 and not 1e+6 on these charts? Thanks for your help. Rishi You can use the axis() function to draw axes with any set of labels you want. -thomas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] L é vy distribution?
Gabriel Erbano wrote: Hi all, Does R support the Lévy distribution? I looked everywhere and couldn't find it. Thanks for any help! Gabriel I believe Jim Lindsey's rmutil package supports this though I've never used it: http://popgen0146uns50.unimaas.nl/~jlindsey/rcode.html BTW, I found this out by doing an R site search at http://cran.r-project.org/ - Search - R site search and typing in Levy distribution. I got 13 hits the first of which was Lindsey's help page. --sundar __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Error with arima()
It means that the CSS estimates from the model are invalid. That's possible as CSS does not enforce validity. It probably means the model is inappropriate. On Thu, 17 Jun 2004, Dan Bebber wrote: Could someone please give a brief explanation, or pointer to an explanation, of the following error: arima(ts.growth, order = c(1,0,0),include.mean=T) Error in arima(ts.growth, order = c(1, 0, 0), include.mean = T) : non-stationary AR part from CSS and why it does not arise with arima0(ts.growth, order = c(1,0,0)) Not the same algorithm. These algorithms only find a local minimum. You are asking the question the wrong way round: unless you know there is uniquely defined estimator, why would you expect the same results? Especially if the model is a poor fit. -- 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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] R help in Firefox on Windows XP
I had to reinstall my machine, so I installed Firefox 0.9 as browser I am using WinXP and R 1.9.1 beta. Now search in R html help does not work. I checked that the Java VM is working correctlt, Sun's test site says my installation is OK. Firefoxalso tells me that Applet Searchengine loaded Applet Searchengine started it just does not find anything. Does anybody know how to solve this? Erich __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] beta regression in R
Hello,
I'm using optim to program a set of mle regression procedures for non-normal
disturbances. This is for teaching and expository purposes only. I've
successfully programmed the normal, generalized gamma, gamma, weibull,
exponential, and lognormal regression functions. And optim returns
reasonable answers for all of these compared with the identical optimization
problems in STATA and LIMDEP. However, I'm having a problem with beta
regression. The parameterization of beta regression that I am using solves
for the second distributional parameter in terms of the first and xbeta.
Therefore, there is only a single distribution parameter. This
parameterization works fine in STATA, LIMDEP, and SHAZAM. However, the
results returned by optim in R are not reasonable in terms of the value of
the log likelihood and parameter estimates.
Here is my code. Does anyone see a problem with what I'm doing here? Any
advice would be appreciated.
Thanks.
B. Dan Wood
# Define the function to be optimized
llik.beta - function(par,X,YP) {
X - as.matrix(x)
YP - as.vector(y)
xbeta - X %*% par[1:K]
p - par[K1:K1]
sum(
-lgamma(p)
+lgamma(p+(p/xbeta-p))
-lgamma(p/xbeta-p)
+(p-1)*log(YP)
+log(1-YP)*(p/xbeta-p-1)
)
}
llik.beta
# Now use the above function to estimate the model. First, create a set of
reasonable start values.
startvalues - c(coefficients(ols.model),74)
startvalues
# Now call optim
mod.beta - optim(startvalues,llik.beta, method = BFGS, control =
list(trace,
maxit=1000,fnscale = -1), hessian = TRUE)
mod.beta
B. Dan Wood, Professor and University Faculty Fellow
Texas AM University
Department of Political Science
4348 TAMU
College Station, TX 77843-4348
Office: (979) 845-1610
Home: (979) 690-0390
Voice Mail: (979) 492-7599
FAX: (979) 847-8924
http://www-polisci.tamu.edu/bdanwood http://www-polisci.tamu.edu/bdanwood
[[alternative HTML version deleted]]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] L é vy distribution?
Actually, I tried searching the R site, but it just returned 4 hits, and one of them was about Lindsey's package. But since it was not listed in the packages section, I thought it was broken... Well, my bad! Thanks! From: Sundar Dorai-Raj [EMAIL PROTECTED] Organization: PDF Solutions, Inc. Reply-To: [EMAIL PROTECTED] Date: Thu, 17 Jun 2004 10:00:25 -0700 To: Gabriel Erbano [EMAIL PROTECTED] Cc: [EMAIL PROTECTED] Subject: Re: [R] L é vy distribution? Gabriel Erbano wrote: Hi all, Does R support the Lévy distribution? I looked everywhere and couldn't find it. Thanks for any help! Gabriel I believe Jim Lindsey's rmutil package supports this though I've never used it: http://popgen0146uns50.unimaas.nl/~jlindsey/rcode.html BTW, I found this out by doing an R site search at http://cran.r-project.org/ - Search - R site search and typing in Levy distribution. I got 13 hits the first of which was Lindsey's help page. --sundar __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Question on lists and vectors of lists
I have an elementary programming question. Could someone please point
me in the right direction?
I have a function which will run for thousands of companies. At each
invocation, it returns 2 numbers. I plan to do something like:
think_one_firm - function(filename) {
# Do stuff
return(list(x=x,y=y))
}
So for each of the firms in my dataset, I will call
think_one_firm(datafilename) and it will give me back two numbers x
and y. I could say:
l = think_one_firm(blah)
print(l$x); print(l$y);
and all would be fine.
What I want to do is: To tuck away the returned lists into a vector or
a data frame.
At the end, I would like to endup with a data structure allfirms
containing one 'row' for each firm. I guess this could be a data
frame, or a matrix, or a vector of lists (if such a thing exists). In
case it's a data frame, I would say allfirms$x[400] to access the 'x'
value returned for the 400th firm.
How would I do this?
--
Ajay Shah Consultant
[EMAIL PROTECTED] Department of Economic Affairs
http://www.mayin.org/ajayshah Ministry of Finance, New Delhi
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] disappointed with x-axes in hist and density plots
TL == Thomas Lumley [EMAIL PROTECTED]
on Thu, 17 Jun 2004 09:53:33 -0700 (PDT) writes:
TL On Thu, 17 Jun 2004, Rishi Ganti wrote:
I've got a few issues with the x-axes in the histogram
and density plots. First, often the default x-axis
doesn't even extend to the length of my data. R often
draws histogram bars (or density lines) farther than the
drawn x-axis extends. For example, I might have a
histogram bar at -15,000. But I wouldn't know that,
because the most negative number on the x-axis is
-10,000. The second issue is the use of scientific
notation. Yes I can read it, but I don't prefer it. Is
there any way for R just to print out 100 and not
1e+6 on these charts? Thanks for your help. Rishi
TL You can use the axis() function to draw axes with any set of labels you
TL want.
and I've recently written the following to show someone how to
beautify the situation (when 1en labels appear: BTW, under
windows there's an etra space in there which IMO makes
the labels much uglier) :
###- Do a 10^k labeling instead of a ek ---
x - 1e7*(-10:50)
y - dnorm(x, m=10e7, s=20e7)
plot(x,y)
axTexpr - function(side, at = axTicks(side, axp=axp, usr=usr, log=log),
axp = NULL, usr = NULL, log = NULL)
{
## Purpose: Do a 10^k labeling instead of a ek
##this auxiliary should return 'at' and 'label' (expression)
## --
## Arguments: as for axTicks()
## --
## Author: Martin Maechler, Date: 7 May 2004, 18:01
eT - floor(log10(abs(at)))# at == 0 case is dealt with below
mT - at / 10^eT
ss - lapply(seq(along = at),
function(i) if(at[i] == 0) quote(0) else
substitute(A %*% 10^E, list(A=mT[i], E=eT[i])))
do.call(expression, ss)
}
plot(x,y, axes= FALSE, frame=TRUE)
aX - axTicks(1); axis(1, at=aX, label= axTexpr(1, aX))
if(FALSE) # rather the next one
aY - axTicks(2); axis(2, at=aY, label= axTexpr(2, aY))
## or rather (horizontal labels on y-axis):
aY - axTicks(2); axis(2, at=aY, label= axTexpr(2, aY), las=2)
--
I hope this decreases your deception..
Further note that you can always call box() after hist()
which may also improve the picture to your eyes.
Regards,
Martin Maechler
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: Re: [R] disappointed with x-axes in hist and density plots
Thanks, but even with axis() I can't get the x-axis to extend to the sides. Try, e.g., x = rnorm(1000) you should have some values in excess of 3 (or below -3). I want to draw the x-axis from -4 to 4, thus encapsulating all points. axis(1,-4:4) but it won't draw. It TRIES to draw it, but I don't see a -4 or 4 on the plot. - Original Message - From: Thomas Lumley Sent: 6/17/2004 9:53:33 AM To: [EMAIL PROTECTED] Cc: [EMAIL PROTECTED] Subject: Re: [R] disappointed with x-axes in hist and density plots On Thu, 17 Jun 2004, Rishi Ganti wrote: I've got a few issues with the x-axes in the histogram and density plots. First, often the default x-axis doesn't even extend to the length of my data. R often draws histogram bars (or density lines) farther than the drawn x-axis extends. For example, I might have a histogram bar at -15,000. But I wouldn't know that, because the most negative number on the x-axis is -10,000. The second issue is the use of scientific notation. Yes I can read it, but I don't prefer it. Is there any way for R just to print out 100 and not 1e+6 on these charts? Thanks for your help. Rishi You can use the axis() function to draw axes with any set of labels you want. -thomas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R help in Firefox on Windows XP
On Thu, 2004-06-17 at 12:06, Erich Neuwirth wrote: I had to reinstall my machine, so I installed Firefox 0.9 as browser I am using WinXP and R 1.9.1 beta. Now search in R html help does not work. I checked that the Java VM is working correctlt, Sun's test site says my installation is OK. Firefoxalso tells me that Applet Searchengine loaded Applet Searchengine started it just does not find anything. Does anybody know how to solve this? Erich Erich, Do you also have JavaScript enabled in the Firefox Tools - Options settings? Both Java and JavaScript need to be enabled for the help.start() search engine to function properly. I reviewed the release notes at http://www.mozilla.org/products/firefox/releases/ and did not see anything relating to Java there, as had been the case with prior releases. The message that you are getting on the status line suggests that the R search applet is being found and properly enabled, which is typically the primary source of problems. Check the above and let us know. Marc Schwartz __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Question on lists and vectors of lists
Try this :
silly.fn - function(x){ return( c(x^2, x^3) ) }
output - matrix( nr=100, nc=2 )
for(i in 1:100){
my.x - rnorm(1)
output[i, ] - silly.fn( my.x )
rownames(output)[i] - paste(Dataset, i, sep=)
}
colnames(output) - c(squared, cubed)
Notice that silly.fn returns a vector which is set to be the ith row of
output.
It is also possible to speed this with the apply family but this depends
on your input format. For example if your input was a list where the
first element corresponds to 1st company information in dataframe format
sapply( my.list.of.matrices, function(single.mat) some.fn(single.mat) )
Another option is to use output[i, ] - unlist(think_one_firm(blah))
where unlist hopefully will convert the list to vector.
On Thu, 2004-06-17 at 18:26, Ajay Shah wrote:
I have an elementary programming question. Could someone please point
me in the right direction?
I have a function which will run for thousands of companies. At each
invocation, it returns 2 numbers. I plan to do something like:
think_one_firm - function(filename) {
# Do stuff
return(list(x=x,y=y))
}
So for each of the firms in my dataset, I will call
think_one_firm(datafilename) and it will give me back two numbers x
and y. I could say:
l = think_one_firm(blah)
print(l$x); print(l$y);
and all would be fine.
What I want to do is: To tuck away the returned lists into a vector or
a data frame.
At the end, I would like to endup with a data structure allfirms
containing one 'row' for each firm. I guess this could be a data
frame, or a matrix, or a vector of lists (if such a thing exists). In
case it's a data frame, I would say allfirms$x[400] to access the 'x'
value returned for the 400th firm.
How would I do this?
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Question on lists and vectors of lists
On Thu, 17 Jun 2004 22:56:29 +0530 Ajay Shah wrote:
I have an elementary programming question. Could someone please point
me in the right direction?
I have a function which will run for thousands of companies. At each
invocation, it returns 2 numbers. I plan to do something like:
think_one_firm - function(filename) {
# Do stuff
return(list(x=x,y=y))
}
So for each of the firms in my dataset, I will call
think_one_firm(datafilename) and it will give me back two numbers x
and y. I could say:
l = think_one_firm(blah)
print(l$x); print(l$y);
and all would be fine.
What I want to do is: To tuck away the returned lists into a vector or
a data frame.
At the end, I would like to endup with a data structure allfirms
containing one 'row' for each firm. I guess this could be a data
frame, or a matrix, or a vector of lists (if such a thing exists). In
case it's a data frame, I would say allfirms$x[400] to access the 'x'
value returned for the 400th firm.
How would I do this?
If I understand correctly what you are trying to do, you could do the
following:
Suppose you have some firms:
firms - c(blah, foo, bar)
and a function which computes something based on these firms
thinkOneFirm - function(firm) list(x = rnorm(1), y = rnorm(1))
(which for this example just returns random numbers, but let's imagine
it's something more sensible :))
And then
rval1 - sapply(firms, thinkOneFirm)
gives you a matrix with the numeric results. To get the desired
data.frame you need to transpose and coerce
rval2 - as.data.frame(t(rval1))
and then you can get value x for firm foo by
rval2$x[foo]
or
rval2[foo, x]
hth,
Z
--
Ajay Shah Consultant
[EMAIL PROTECTED] Department of Economic Affairs
http://www.mayin.org/ajayshah Ministry of Finance, New Delhi
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R help in Firefox on Windows XP
Erich Neuwirth wrote: I had to reinstall my machine, so I installed Firefox 0.9 as browser I am using WinXP and R 1.9.1 beta. Now search in R html help does not work. A workaround (which is slow, provisional, and as yet untested on your configuration) is to use a different help engine, namely source(http://www.statslab.cam.ac.uk/~djw1005/Stats/Interests/search.R;) helpHTML() Damon. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] problem with restore and some .RData
Hi, I have problem with the restore function in some .RData using R 1.9.0 Look the error: [EMAIL PROTECTED] RAnalise]$ R ... Error: object 'family' not found whilst loading namespace 'MASS' Fatal error: unable to restore saved data in .RData But if I load this .RData with the load() function all objects are recovered. [EMAIL PROTECTED] RAnalise]$ R --no-restore-data ... ls() character(0) load(.RData) [1] .Traceback RespY [3] Trat1Trat2 ... [63] varied.371.84210116.9varied.371.842101169.5 [65] varied.371.8421011695.12 All recovered objects work fine. I can use q() and save workspace image to quit session. But to load this .RData I need to use R with --no-restore-data option and load the .RData manually. I try to find the problem, but I dont find this. Anybody have any idea? Thanks Ronaldo ps. how to initiate R in xemacs with the --no-restore-data option? M-x R --no-restore-data dont work -- Hlade's Law: If you have a difficult task, give it to a lazy person -- they will find an easier way to do it. -- | // | \\ [***] | ( õ õ ) [Ronaldo Reis Júnior] | V [UFV/DBA-Entomologia] |/ \ [36571-000 Viçosa - MG ] | /(.''`.)\ [Fone: 31-3899-2532 ] | /(: :' :)\ [EMAIL PROTECTED]] |/ (`. `'` ) \[ICQ#: 5692561 | LinuxUser#: 205366 ] |( `- ) [***] | _/ \_Powered by GNU/Debian Woody/Sarge __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: Re: [R] disappointed with x-axes in hist and density plots
You haven't read ?axis, I guess. Try using the at= argument. Whatever you do not like about the default output, chances are you can customize it to your heart's content... Andy From: Rishi Ganti Thanks, but even with axis() I can't get the x-axis to extend to the sides. Try, e.g., x = rnorm(1000) you should have some values in excess of 3 (or below -3). I want to draw the x-axis from -4 to 4, thus encapsulating all points. axis(1,-4:4) but it won't draw. It TRIES to draw it, but I don't see a -4 or 4 on the plot. - Original Message - From: Thomas Lumley Sent: 6/17/2004 9:53:33 AM To: [EMAIL PROTECTED] Cc: [EMAIL PROTECTED] Subject: Re: [R] disappointed with x-axes in hist and density plots On Thu, 17 Jun 2004, Rishi Ganti wrote: I've got a few issues with the x-axes in the histogram and density plots. First, often the default x-axis doesn't even extend to the length of my data. R often draws histogram bars (or density lines) farther than the drawn x-axis extends. For example, I might have a histogram bar at -15,000. But I wouldn't know that, because the most negative number on the x-axis is -10,000. The second issue is the use of scientific notation. Yes I can read it, but I don't prefer it. Is there any way for R just to print out 100 and not 1e+6 on these charts? Thanks for your help. Rishi You can use the axis() function to draw axes with any set of labels you want. -thomas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] disappointed with x-axes in hist and density plots
On Thu, 17 Jun 2004 10:28:57 -0700 Rishi Ganti wrote: Thanks, but even with axis() I can't get the x-axis to extend to the sides. Try, e.g., x = rnorm(1000) you should have some values in excess of 3 (or below -3). I want to draw the x-axis from -4 to 4, thus encapsulating all points. axis(1,-4:4) but it won't draw. It TRIES to draw it, but I don't see a -4 or 4 on the plot. Well you need to make enough space before! When you have got R x - rnorm(1000) R y - rnorm(1000) you need to make sure that the desired range is covered by the plot: R plot(x, y, axes = FALSE, xlim = c(-4, 4)) then you add the x-axis R axis(1, at = -4:4) and y-axis and a box. R axis(2) R box() Consult the manuals and ?plot, ?par for more information. hth, Z - Original Message - From: Thomas Lumley Sent: 6/17/2004 9:53:33 AM To: [EMAIL PROTECTED] Cc: [EMAIL PROTECTED] Subject: Re: [R] disappointed with x-axes in hist and density plots On Thu, 17 Jun 2004, Rishi Ganti wrote: I've got a few issues with the x-axes in the histogram and density plots. First, often the default x-axis doesn't even extend to the length of my data. R often draws histogram bars (or density lines) farther than the drawn x-axis extends. For example, I might have a histogram bar at -15,000. But I wouldn't know that, because the most negative number on the x-axis is -10,000. The second issue is the use of scientific notation. Yes I can read it, but I don't prefer it. Is there any way for R just to print out 100 and not 1e+6 on these charts? Thanks for your help. Rishi You can use the axis() function to draw axes with any set of labels you want. -thomas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] disappointed with x-axes in hist and density plots
x - rnorm(1000) y - rnorm(1000) plot(x,y) axis(1,-4,4) (speculation that attempting above ... not what you want to do...) rather do plot(x,y,xlim = c(-4,4)) Rishi Ganti wrote: Thanks, but even with axis() I can't get the x-axis to extend to the sides. Try, e.g., x = rnorm(1000) you should have some values in excess of 3 (or below -3). I want to draw the x-axis from -4 to 4, thus encapsulating all points. axis(1,-4:4) but it won't draw. It TRIES to draw it, but I don't see a -4 or 4 on the plot. - Original Message - From: Thomas Lumley Sent: 6/17/2004 9:53:33 AM To: [EMAIL PROTECTED] Cc: [EMAIL PROTECTED] Subject: Re: [R] disappointed with x-axes in hist and density plots On Thu, 17 Jun 2004, Rishi Ganti wrote: I've got a few issues with the x-axes in the histogram and density plots. First, often the default x-axis doesn't even extend to the length of my data. R often draws histogram bars (or density lines) farther than the drawn x-axis extends. For example, I might have a histogram bar at -15,000. But I wouldn't know that, because the most negative number on the x-axis is -10,000. The second issue is the use of scientific notation. Yes I can read it, but I don't prefer it. Is there any way for R just to print out 100 and not 1e+6 on these charts? Thanks for your help. Rishi You can use the axis() function to draw axes with any set of labels you want. -thomas __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R help in Firefox on Windows XP
Works for me with Firefox 0.8. James W. MacDonald Affymetrix and cDNA Microarray Core University of Michigan Cancer Center 1500 E. Medical Center Drive 7410 CCGC Ann Arbor MI 48109 734-647-5623 Prof Brian Ripley [EMAIL PROTECTED] 06/17/04 02:10PM Works for me. Did it work with Firefox 0.8? (I upgraded from 0.8.) I would try 0.8 and then upgrade. On Thu, 17 Jun 2004, Erich Neuwirth wrote: I had to reinstall my machine, so I installed Firefox 0.9 as browser I am using WinXP and R 1.9.1 beta. Now search in R html help does not work. I checked that the Java VM is working correctlt, Sun's test site says my installation is OK. Firefoxalso tells me that Applet Searchengine loaded Applet Searchengine started it just does not find anything. Does anybody know how to solve this? -- 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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] disappointed with x-axes in hist and density plots
Rishi Ganti [EMAIL PROTECTED] writes: Thanks, but even with axis() I can't get the x-axis to extend to the sides. Try, e.g., x = rnorm(1000) you should have some values in excess of 3 (or below -3). I want to draw the x-axis from -4 to 4, thus encapsulating all points. axis(1,-4:4) but it won't draw. It TRIES to draw it, but I don't see a -4 or 4 on the plot. Well, the x limits likely don't extend that far. (You have noticed by now that the axis extents are not the same as the sides of the bounding box, right?) So set them: hist(x,xlim=c(-4,4)) Or generically: hist(x,xlim=range(pretty(x))) tends to give the kind of external axis that you're requesting. Less efficient use of space though (since the axis is always longer than it needed to be). -- O__ Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] nlme graphics in a loop problem
Could it be that you need print(plot(model)) inside the loop? I believe
plot() methods in nlme are mostly lattice graphics, which needs to be
explicitly print()ed inside functions and loops.
Andy
From: Chris Knight
Hi, I'm fitting mixed effects models using the lme function
of the nlme
package. This involves using the various associated plot functions.
However, when I attempt to fit a number of models using an
loop, whilst
the models work, the plot functions fail. As a trivial example, the
following works:
library(nlme)
graphics.off()
x-c(1:10)
y-c(1:4,7:12)
b-as.factor(c(rep(A,5),rep(B,5)))
model-lme(y~x, random=~1|b)
plot(model)
however the following, identical code in a loop, doesn't:
graphics.off()
for(i in 1:2){
x-c(1:10)
y-c(1:4,7:12)
b-as.factor(c(rep(A,5),rep(B,5)))
model-lme(y~x, random=~1|b)
plot(model)
}
Mostly this is only inconvenient, since a similar plot may be created
successfully within the loop using plot(fitted(model),resid(model)),
however, I'd be keen to find out whether this is a general
problem/sign
of something deeper or I'm just missing something easy that could sort
it out.
Thanks,
Chris
--
~~
Dr. Christopher Knight Tel:+44
1865 275111
Dept. Plant Sciences +44
1865 275790
South Parks Road
Oxford OX1 3RB Fax:+44
1865 275074
` · . , ,(((º
~~
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] nlme graphics in a loop problem
yes, use the print command, and if you can save each plot in the loop,
that will work.
I had to do this for a loop and tables I was making. I had to save the
first table to a matrix, then bind each successive table to that
matrix, then print the final matrix. Not sure how to do this for a
plot, though...
On Jun 17, 2004, at 12:27 PM, Liaw, Andy wrote:
Could it be that you need print(plot(model)) inside the loop? I
believe
plot() methods in nlme are mostly lattice graphics, which needs to be
explicitly print()ed inside functions and loops.
Andy
From: Chris Knight
Hi, I'm fitting mixed effects models using the lme function
of the nlme
package. This involves using the various associated plot functions.
However, when I attempt to fit a number of models using an
loop, whilst
the models work, the plot functions fail. As a trivial example, the
following works:
library(nlme)
graphics.off()
x-c(1:10)
y-c(1:4,7:12)
b-as.factor(c(rep(A,5),rep(B,5)))
model-lme(y~x, random=~1|b)
plot(model)
however the following, identical code in a loop, doesn't:
graphics.off()
for(i in 1:2){
x-c(1:10)
y-c(1:4,7:12)
b-as.factor(c(rep(A,5),rep(B,5)))
model-lme(y~x, random=~1|b)
plot(model)
}
Mostly this is only inconvenient, since a similar plot may be created
successfully within the loop using plot(fitted(model),resid(model)),
however, I'd be keen to find out whether this is a general
problem/sign
of something deeper or I'm just missing something easy that could sort
it out.
Thanks,
Chris
--
~~
Dr. Christopher Knight Tel:+44
1865 275111
Dept. Plant Sciences +44
1865 275790
South Parks Road
Oxford OX1 3RB Fax:+44
1865 275074
` · . , ,(((º
~~
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
Jennifer A. Emond, MS
Statistician
Department of Biostatistics, UCSD
9500 Gilman Drive
M/C 0949
La Jolla, CA 92093-0949
(858) 622-5877
[EMAIL PROTECTED]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] nlme graphics in a loop problem
Use array for your tables and list for plots; e.g.,
allTable - array(NA, dim=c(numRow, numCol, numIter))
allPlot - vector(mode=list, length=numIter)
for (i in 1:numIter) {
...
allPlot[[i]] - plot(...)
allTable[,,i] - however you generate the table
}
HTH,
Andy
From: Jennifer Emond
yes, use the print command, and if you can save each plot in
the loop,
that will work.
I had to do this for a loop and tables I was making. I had
to save the
first table to a matrix, then bind each successive table to that
matrix, then print the final matrix. Not sure how to do this for a
plot, though...
On Jun 17, 2004, at 12:27 PM, Liaw, Andy wrote:
Could it be that you need print(plot(model)) inside the loop? I
believe
plot() methods in nlme are mostly lattice graphics, which
needs to be
explicitly print()ed inside functions and loops.
Andy
From: Chris Knight
Hi, I'm fitting mixed effects models using the lme function
of the nlme
package. This involves using the various associated plot functions.
However, when I attempt to fit a number of models using an
loop, whilst
the models work, the plot functions fail. As a trivial example, the
following works:
library(nlme)
graphics.off()
x-c(1:10)
y-c(1:4,7:12)
b-as.factor(c(rep(A,5),rep(B,5)))
model-lme(y~x, random=~1|b)
plot(model)
however the following, identical code in a loop, doesn't:
graphics.off()
for(i in 1:2){
x-c(1:10)
y-c(1:4,7:12)
b-as.factor(c(rep(A,5),rep(B,5)))
model-lme(y~x, random=~1|b)
plot(model)
}
Mostly this is only inconvenient, since a similar plot may
be created
successfully within the loop using
plot(fitted(model),resid(model)),
however, I'd be keen to find out whether this is a general
problem/sign
of something deeper or I'm just missing something easy
that could sort
it out.
Thanks,
Chris
--
~~
Dr. Christopher Knight Tel:+44
1865 275111
Dept. Plant Sciences +44
1865 275790
South Parks Road
Oxford OX1 3RB Fax:+44
1865 275074
` · . , ,(((º
~~
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
Jennifer A. Emond, MS
Statistician
Department of Biostatistics, UCSD
9500 Gilman Drive
M/C 0949
La Jolla, CA 92093-0949
(858) 622-5877
[EMAIL PROTECTED]
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] problem with restore and some .RData
The short answer is to add importFrom(stats, family) to MASS's NAMESPACE and reinstall. Another way is to use library(stats) in your .Rprofile in that directory. The real question is why your .RData is loading namespace MASS. Presumably some object in your workspace has MASS in its environment. Please find out which and tell us exactly how you created it. (Have you copied a MASS object, for example?) On Thu, 17 Jun 2004, Ronaldo Reis Jr. wrote: Hi, I have problem with the restore function in some .RData using R 1.9.0 Look the error: [EMAIL PROTECTED] RAnalise]$ R ... Error: object 'family' not found whilst loading namespace 'MASS' Fatal error: unable to restore saved data in .RData But if I load this .RData with the load() function all objects are recovered. [EMAIL PROTECTED] RAnalise]$ R --no-restore-data ... ls() character(0) load(.RData) [1] .Traceback RespY [3] Trat1Trat2 ... [63] varied.371.84210116.9varied.371.842101169.5 [65] varied.371.8421011695.12 All recovered objects work fine. I can use q() and save workspace image to quit session. But to load this .RData I need to use R with --no-restore-data option and load the .RData manually. I try to find the problem, but I dont find this. -- 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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Re: Clustering in R
Thank you for all the responses. I've downloaded TIGR MeV, it seems to do
everything I need it to do. The only problem is that none of the algorithms
seem to work (K-means, hierarchical) giving errors that just say -2. I think
probably the reason for this is that I'm running linux (since there aren't any
available windows machines), and MeV was debugged in windows. It seems that
the lab is getting a new windows machine soon so at that time I'll be able to
try it out.
In the meantime, I'm going to see if I can get workable data by removing
unexpressed genes and using hclust with cutree to get into each of the
branches. From the replies I've gotten, this should probably work, I just
need to sit down and try it. If I get stuck on either of these, I know where
to ask for more help. At least now I have viable directions. Thanks again
for everyone who responded, it was enormously helpful.
Wayne
Quoting Martin Maechler [EMAIL PROTECTED]:
Thanks a lot, Michael!
I cc to R-help, where this question really belongs {as the
'Subject' suggests itself...} -- please drop 'bioconductor' from
CC'ing further replies.
michael == michael watson (IAH-C) [EMAIL PROTECTED]
on Thu, 17 Jun 2004 09:16:59 +0100 writes:
michael OK, admittedly it is not incredibly simple, but it
michael is not *that* difficult.
michael If you are familiar with R, it should take you an
michael hour or two; if unfamiliar, perhaps a day or two.
michael The commands you want (and need to read the help on) are:
michael hclust
michael plclust
michael cutree
and I would add identify.hclust() {and rect.hclust()}
a very neat but not known / used enough function
a link to which I have just added to the help(hclust) page.
Look at its examples {not with example() since they are
dontrun} correcting the extraneous . in the last (and
coolest!) example!
michael dendrogram
michael as.dendrogram
michael heatmap
where you use dendrograms produced from hclust objects via
as.dendrogram(hc-obj) or also twins objects produced from
package cluster's agnes() or diana() via
as.dendrogram(as.hclust( twins-obj ) )
help(dendrogram) also mentions
[[ (and shows examples) and cut() for cutting dendrograms and shows
how you can depict dendrograms into its parts.
michael With intelligent use of hclust - cutree - subsetting -
hclust
michael (in that order) you will be able to drill down
michael into your dendrogram and create sub-trees - until
michael you get to the level where you can see your gene
michael names.
or also
hclust - as.dendrogram - cut - ..
- [[ -
Note that there also is reorder.dendrogram() for reordering
dendrogram nodes ``sensibly'' --- something that heatmap() does,
but you can play with quite a bit.
Further, note Catherine Hurley's gclus package which
orders/reorders 'hclust' objects directly, but with a more
interesting algorithm.
Note that I'd strongly recommend to use R 1.9.1 beta for these,
since I know which bugs in the dendrogram code I have fixed
since R 1.9.0...
michael An important message to take home here is that if
michael you have 14000 genes and therefore 14000 labels,
michael it's going to be difficult to display your tree in
michael ANY software, including the expensive commercial products.
not showing the labels and using identify.hclust() and the
command line to extract the indices of observations in
clusters (and subclusters) and visualize them in other, non-dendrogram
plots,
might well be feasible.
michael Let me know how you get on
michael Thanks
michael Mick
michael -Original Message-
michael From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
michael Sent: 16 June 2004 21:26
michael To: [EMAIL PROTECTED]
michael Subject: [BioC] Clustering in R
Dear list members,
I'm an undergrad and I work in a lab at Brandeis.
I am trying to cluster around 14,000 genes across 6
microarray experiments. Two of these experiments
are replicates. I have decided to use R since it
seems to be the most complete and flexible software
package for normalization and clustering of
microarray data.
The problem is that I am new to clustering and to
R. Just to mention of a few of the problems I'm
having: the dendrogram that is drawn by R from the
agnes object is far too dense to see any of the
gene names; kmeans won't work, returning an error
saying that my data has NAs in it (there weren't
any missing values in the original table though);
I'd like to be able to see a heatmap or a
cumulative plot of expression profiles for genes
that are clustered together or are on the same
branch of the dendrogram.
I know that these questions are probably very
simple, but I can't
Re: [R] Resolution of plots
On Thu, 17-Jun-2004 at 08:00AM -0700, Thomas Lumley wrote: | On Wed, 16 Jun 2004, Prof Brian Ripley wrote: | | You will have to tell us more. Exporting how: to what format using what | device and what exact command on what operating system? | | The only device I know of that even knows about dpi is bitmap() and that | has no such limit unless imposed by your implementation of ghostscript. | | There is an issue with PNG. libpng provides png_set_pHYs to set resolution | (in pixels/metre) but provides a default if it is not set. We don't set | it, and so get the default resolution. I don't follow. Is that in relation to the function bitmap()0 or to png()? -- Patrick Connolly HortResearch Mt Albert Auckland New Zealand Ph: +64-9 815 4200 x 7188 ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~ I have the world`s largest collection of seashells. I keep it on all the beaches of the world ... Perhaps you`ve seen it. ---Steven Wright ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~ __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Resolution of plots
On Fri, 18 Jun 2004, Patrick Connolly wrote: On Thu, 17-Jun-2004 at 08:00AM -0700, Thomas Lumley wrote: | On Wed, 16 Jun 2004, Prof Brian Ripley wrote: | | You will have to tell us more. Exporting how: to what format using what | device and what exact command on what operating system? | | The only device I know of that even knows about dpi is bitmap() and that | has no such limit unless imposed by your implementation of ghostscript. | | There is an issue with PNG. libpng provides png_set_pHYs to set resolution | (in pixels/metre) but provides a default if it is not set. We don't set | it, and so get the default resolution. I don't follow. Is that in relation to the function bitmap() or to png()? png(). That records a nominal 72dpi in the info, although it isn't much used even by PhotoShop/GIMP etc. (Almost all digital cameras record JPEGs set to 72dpi, for example. If you want to use such a program to resize to a physical size it is very easy to edit their idea of the dpi.) The original question was about a 96dpi limit and I still await elucidation of what that was about. -- 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 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Since package.contents is deprecated, what is replacement function?
I did a ?package.contents in version 1.9.0. The help file says this is deprecated, and says to see also Deprecated and Defunct, but what is the new function that replaces packge.contents? The function still seems to work, but it's just tagged as deprecated. efg [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] [TIP] Apropos resolution of plots
If you need to make JPEG plots for the web or suchlike, here's a method for poor mans antialiasing that seems to turn out rather nice: [Requires ghostscript, netpbm] In R (lifted out of example(contour)) bitmap(out.ppm,ppmraw,res=4*72, pointsize=12) data(volcano) rx - range(x - 10*1:nrow(volcano)) ry - range(y - 10*1:ncol(volcano)) ry - ry + c(-1,1) * (diff(rx) - diff(ry))/2 tcol - terrain.colors(12) par(pty = s, bg = lightcyan) plot(x = 0, y = 0,type = n, xlim = rx, ylim = ry, xlab = , ylab = ) u - par(usr) rect(u[1], u[3], u[2], u[4], col = tcol[8], border = red) contour(x, y, volcano, col = tcol[2], lty = solid, add = TRUE, vfont = c(sans serif, plain)) title(A Topographic Map of Maunga Whau, font = 4) abline(h = 200*0:4, v = 200*0:4, col = white, lty = 2, lwd= 0.2) dev.off() and then from the command line pnmsmooth -size 5 5 out.ppm smooth.ppm pnmscale .25 smooth.ppm aa.ppm pnmtojpeg -quality 100 aa.ppm aa.jpg Notice that the intermediate files tend to get rather large even though the end result is quite compact, so don't forget to clean up: $ ls -l *.ppm *.jpg -rw-rw-r--1 pd pd 75393 Jun 18 00:34 aa.jpg -rw-rw-r--1 pd pd 559887 Jun 18 00:34 aa.ppm -rw-rw-r--1 pd pd8958022 Jun 18 00:34 out.ppm -rw-rw-r--1 pd pd8957969 Jun 18 00:34 smooth.ppm You could in principle use a lower -quality in the final jpeg conversion, but at least to my eyes some unsightly artifacts creep in rather quickly. -- O__ Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Manova question.
On Jun 15, 2004, at 10:11 AM, Douglas Bates wrote: Jens Schumacher wrote: knussear wrote: Hi list, I'm attempting to re-create a Repeated Measures Compositional Analysis as described in the work by Aebischer et. al. (Ecology. 1993. 74(5): 1313-1325). In this paper they describe transitions of data into a log ratio difference matrix, from which they obtain two matrices using a monova routine. I am able to produce the second of the two matrices, but I'm having trouble with the first. the difference matrix going in is given here. AnimalScrubBl woodCon wood Grass 10.970-2.380-5.154-9.408 21.217-0.173-4.955-5.521 31.178-0.248-4.0890.338 40.5200.466-4.801-1.946 58.4459.31910.7538.171 68.6549.32710.7328.152 78.4299.35010.8188.141 89.1209.5653.8138.127 99.2279.8823.8137.779 109.4238.0863.8138.539 119.6269.3923.8138.135 129.2348.3023.8138.537 138.6728.9089.8328.416 And the first of the matrices is given here, and is matrix of mean-corrected sums of squares and cross products calculated from the difference matrix. ScrubBl woodCon wood Grass Scrub179.52214.59244.58273.75 Bl wood214.59268.44314.35343.86 Con wood 244.58314.35471.09400.22 Grass273.75343.86400.22477.78 From manova on the data set I can get the diagonal of the matrix, but not the others. manova(y ~ NULL) Terms: Residuals Scrub179.5273 Bl.wood 268.4347 Con.wood 471.0845 Grass477.8014 Deg. of Freedom12 Could anyone offer a suggestion ? Thanks Ken __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Let data.matrixbe the above difference matrix. You obtain the raw sums of squares and cross-products matrix by R2 - t(data.matrix) %*% data.matrix or even R2 - crossprod(data.matrix, data.matrix) or, the preferred form, R2 - crossprod(data.matrix) I'm not sure if this helps but I get the correct matrices if I use JMP. It refers to them as EH matrices. Is there a way to get R to produce these? Thanks Ken __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Is there an easy way to generate linearly independent vectors
On 06/17/04 19:04, Fred wrote:
Dear R-listers:
I am trying to test an algorithm on a set of linearly independent vectors
{x1,x2,...,xn}.
Well, here's an idea, for 10 vectors of length 10,
as columns of a matrix m1. The 11th seems to be needed.
m1 - matrix(rnorm(110),10,11)
for (i in 2:11) {
m1[,i] - resid(lm(m1[,i] ~ m1[, 1:(i-1)]))
}
Test it with cor(m1[,-11])
I'm sure there are better ways.
Of course
m1 - matrix(rnorm(100),10,10)
is ALMOST what you want.
Jon
--
Jonathan Baron, Professor of Psychology, University of Pennsylvania
Home page:http://www.sas.upenn.edu/~baron
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Is there an easy way to generate linearly independent vec tors
I believe eigen(), svd() and qr() can all do it.
Andy
From: Jonathan Baron
On 06/17/04 19:04, Fred wrote:
Dear R-listers:
I am trying to test an algorithm on a set of linearly
independent vectors
{x1,x2,...,xn}.
Well, here's an idea, for 10 vectors of length 10,
as columns of a matrix m1. The 11th seems to be needed.
m1 - matrix(rnorm(110),10,11)
for (i in 2:11) {
m1[,i] - resid(lm(m1[,i] ~ m1[, 1:(i-1)]))
}
Test it with cor(m1[,-11])
I'm sure there are better ways.
Of course
m1 - matrix(rnorm(100),10,10)
is ALMOST what you want.
Jon
--
Jonathan Baron, Professor of Psychology, University of Pennsylvania
Home page:http://www.sas.upenn.edu/~baron
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] can't get text to appear over individual panels in multi-panel plot
On Thursday 17 June 2004 19:37, Patrick Bennett wrote: I'm trying to learn how to create Trellis multi-panel plots, but I'm having some trouble reproducing the graphs shown in Venables Ripley (2002) (e.g., Figs 4.14 4.15). Actually, everything looks fine except for the fact that I can't see any text above the individual panels. I'm using R 1.9.0 for OS-X running on Mac OSX 10.3.3, and I'm drawing the graphs into the Quartz device. Sorry for bothering the list with what is probably a trivial question, but I've run out of ideas... Any help would be appreciated. I don't have that edition, and it's difficult to guess without code, but have you looked at the scripts that come with the MASS package ? They are supposed to reproduce the analysis done in the book (I'm assuming that the book you are referring to is in fact MASS). It might be as simple as replacing text() calls by ltext(). Deepayan __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Is there an easy way to generate linearly independent vec tors
I want to get linearly independent vectors, not orthogonal ones.
The functions eigen, svd, I think it may provide orthogonal
vectors which are not what I expect.
- Original Message -
From: Liaw, Andy [EMAIL PROTECTED]
To: 'Jonathan Baron' [EMAIL PROTECTED]; Fred
[EMAIL PROTECTED]
Cc: R-help mailing list [EMAIL PROTECTED]
Sent: Thursday, June 17, 2004 7:50 PM
Subject: RE: [R] Is there an easy way to generate linearly independent vec
tors
I believe eigen(), svd() and qr() can all do it.
Andy
From: Jonathan Baron
On 06/17/04 19:04, Fred wrote:
Dear R-listers:
I am trying to test an algorithm on a set of linearly
independent vectors
{x1,x2,...,xn}.
Well, here's an idea, for 10 vectors of length 10,
as columns of a matrix m1. The 11th seems to be needed.
m1 - matrix(rnorm(110),10,11)
for (i in 2:11) {
m1[,i] - resid(lm(m1[,i] ~ m1[, 1:(i-1)]))
}
Test it with cor(m1[,-11])
I'm sure there are better ways.
Of course
m1 - matrix(rnorm(100),10,10)
is ALMOST what you want.
Jon
--
Jonathan Baron, Professor of Psychology, University of Pennsylvania
Home page:http://www.sas.upenn.edu/~baron
--
Notice: This e-mail message, together with any attachments, contains
information of Merck Co., Inc. (One Merck Drive, Whitehouse Station, New
Jersey, USA 08889), and/or its affiliates (which may be known outside the
United States as Merck Frosst, Merck Sharp Dohme or MSD and in Japan, as
Banyu) that may be confidential, proprietary copyrighted and/or legally
privileged. It is intended solely for the use of the individual or entity
named on this message. If you are not the intended recipient, and have
received this message in error, please notify us immediately by reply e-mail
and then delete it from your system.
--
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] can't get text to appear over individual panels in multi-panel plot
I neglected to say that I am using the R-Aqua interface and the MASS, grid, lattice packages. Here is one specific example where I'm having trouble. After loading the crabs data set, I create the figure with the following code (which is taken from MASS): lcrabs.pc-predict(princomp(log(crabs[,4:8]))) sex-crabs$sex;levels(sex)-c(Female,Male) sp-crabs$sp;levels(sp)-c(Blue,Orange) splom(~lcrabs.pc[,1:3] | sp*sex,cex=0.5,pscales=0) The figure is plotted in the Quartz window: everything looks OK except for the lack of text above the individual panels. Patrick J. Bennett, Ph.D. Professor Canada Research Chair Department of Psychology McMaster University 1280 Main St. West Hamilton, ON L8S 4K1 Canada Office: 905-525-9140 ext. 23012 FAX: 905-529-6225 URL: www.psychology.mcmaster.ca/bennett/www.chairs.gc.ca Verbosity leads to unclear, inarticulate things. -- Dan Quayle On Jun 17, 2004, at 8:59 PM, Deepayan Sarkar wrote: On Thursday 17 June 2004 19:37, Patrick Bennett wrote: I'm trying to learn how to create Trellis multi-panel plots, but I'm having some trouble reproducing the graphs shown in Venables Ripley (2002) (e.g., Figs 4.14 4.15). Actually, everything looks fine except for the fact that I can't see any text above the individual panels. I'm using R 1.9.0 for OS-X running on Mac OSX 10.3.3, and I'm drawing the graphs into the Quartz device. Sorry for bothering the list with what is probably a trivial question, but I've run out of ideas... Any help would be appreciated. I don't have that edition, and it's difficult to guess without code, but have you looked at the scripts that come with the MASS package ? They are supposed to reproduce the analysis done in the book (I'm assuming that the book you are referring to is in fact MASS). It might be as simple as replacing text() calls by ltext(). Deepayan __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] can't get text to appear over individual panels in multi-panel plot
On Thursday 17 June 2004 22:24, Patrick Bennett wrote: I neglected to say that I am using the R-Aqua interface and the MASS, grid, lattice packages. Here is one specific example where I'm having trouble. After loading the crabs data set, I create the figure with the following code (which is taken from MASS): lcrabs.pc-predict(princomp(log(crabs[,4:8]))) sex-crabs$sex;levels(sex)-c(Female,Male) sp-crabs$sp;levels(sp)-c(Blue,Orange) splom(~lcrabs.pc[,1:3] | sp*sex,cex=0.5,pscales=0) The figure is plotted in the Quartz window: everything looks OK except for the lack of text above the individual panels. Here's what I get (on a pdf device), and it looks OK to me. http://www.stat.wisc.edu/~deepayan/R/crabs.pdf If this is not what you see, then there might be a problem with your device driver. Deepayan __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] msm
Hello, I'm writing about msm. It may be that consistent users of Markov models have a good idea as to what constitutes workable data for a model. I think of general rules, in basic statistical studies where n is limited to exclude fairly precise figures in the lower range. On the other hand Markov models don't seem to be often enough used for parameters to be as well laid out. I also get the feeling that msm is organized to work optimally with certain sizes and shapes of data. Is there a source that anyone is aware of on this? (I have the Nelder text on optimization, and also have a feeling that what's possible is pretty closely connected with optimization questions) Thanks again, Russell [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] can't get text to appear over individual panels in multi-panel plot
On Thursday 17 June 2004 22:57, Patrick Bennett wrote: yes, i can reproduce that same graph when i print to the pdf-device. but the panel titles do not appear when I print to the Quartz-device. Hmm. I won't be able to help you then, let's hope someone else can. Deepayan On Jun 17, 2004, at 11:47 PM, Deepayan Sarkar wrote: On Thursday 17 June 2004 22:24, Patrick Bennett wrote: I neglected to say that I am using the R-Aqua interface and the MASS, grid, lattice packages. Here is one specific example where I'm having trouble. After loading the crabs data set, I create the figure with the following code (which is taken from MASS): lcrabs.pc-predict(princomp(log(crabs[,4:8]))) sex-crabs$sex;levels(sex)-c(Female,Male) sp-crabs$sp;levels(sp)-c(Blue,Orange) splom(~lcrabs.pc[,1:3] | sp*sex,cex=0.5,pscales=0) The figure is plotted in the Quartz window: everything looks OK except for the lack of text above the individual panels. Here's what I get (on a pdf device), and it looks OK to me. http://www.stat.wisc.edu/~deepayan/R/crabs.pdf If this is not what you see, then there might be a problem with your device driver. Deepayan __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
