Re: [R] linear discriminant analysis in MASS
On Mon, 20 Feb 2006, Alain Paquette wrote: Hello R people I now know how to run my discriminant analysis with the lda function in MASS: lda.alain=lda(Groupes ~ Ht.D0 + Lc.Dc + Ram + IDF, gr, CV = FALSE) and it works fine. CV=FALSE is the default and so not needed. But I am missing a test and cannot find any help on how to get it, if it exist. The S equivalent: There is no such function in S, and I rather object as the S equivalent is lda() (and as the author of both I should know). Credit where credit is due: discrim() is an S-PLUS function, indebted to lda(). discrim(structure(.Data = Groupes ~ Ht.D0 + Lc.Dc + Ram + IDF, class = formula), data = gr, family = Canonical(cov.structure = homoscedastic), na.action = na.omit, prior = proportional) outputs a nice matrix of Mahalanobis distances between groups and even tests (Hotelling's T Squared) for significant distances. Well, it seems not to. That is part of the output of the summary() method, which itself calls the multicomp() method. Why don't I just take the S output you say? Because like you, I'd rather put in my paper that I did it using R of course! No `of course' applies. If you learnt of this output from S-PLUS, I urge you to credit it honestly and accurately. (If you refer to lda, you should credit that, not just R.) Does anyone know of a way to get this test out of lda? Or of another R package that does it? Mahalanobis distance between groups is easy, as this is just Euclidean distance between group centres in the scaled space. The test statistics can be produced, but - they are critically dependent on the unrealistic assumptions of multivariate normality and variance homogeneity and - there needs to be an adjustment for multiple comparisons. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Extracting error terms from an lme object
Hi all, I have a problem with extracting the right error terms from an lme-object. The fixed part of the model is quite simple (y~x1*x2, where x1 and x2 are factor variables with two levels), with random intercepts (random=~1|group). As I work with factor variables here, I would like to display this with a barplot (or something like that), with error bars representing the sem or sd of the mean values for each group (i.e. four bars). However, I have not figured out how to extract sem (sd) values for these means, or actually what would be the appropriate error terms to use here. The summary.lme function only provides me with the se for the parameter estimates (intercept and slopes) while I want to display error terms for the different levels. Could this be achieved by setting up an appropriate contrast matrix (or would any of the established contrasts methods work), or is there another way of getting error estimates that I could use in such a figure? All the best, Erlend Nilsen __ R-help@stat.math.ethz.ch 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 close all windows?
?graphics.off is for graphic windows only. You don't tell us, but I assume you are working with RGui under Windows. If you don't want many help windows, you can still change it in Edit - GUI preferences - Pager style - single window. Best, Philippe Grosjean Uwe Ligges wrote: Michael wrote: Hi all, I have 10+ graphic windows and 10+ R help window opening now... How do I close them all at once? Thanks a lot! [[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 See ?graphics.off Uwe Ligges __ R-help@stat.math.ethz.ch 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] Calculate R-Square for existing predictors
Hello everybody, in the past months I have developed a model to predict the mortality rate of Acute Myocardial Infarction Patients. My evaluation was mainly based on ROC Curves etc. I have now been asked to calculate the R-Square of my model, another existing model in regards to the actual mortality rate. So far I have only found a function calculating a regression, and then the R-Square of the calculated Regression. Is there a function to calculate the R-Square for existing models? Thanks for your help Katrin __ R-help@stat.math.ethz.ch 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 close all windows?
Or you can modify Rprofile.site with options(htmlhelp=TRUE) and get your help in HTML browser, which can be sometimes more convenient due to hyperlinks. Best regards Petr On 21 Feb 2006 at 9:58, Philippe Grosjean wrote: Date sent: Tue, 21 Feb 2006 09:58:24 +0100 From: Philippe Grosjean [EMAIL PROTECTED] Organization: SciViews UMH - EcoNum To: Uwe Ligges [EMAIL PROTECTED] Copies to: R-help@stat.math.ethz.ch Subject:Re: [R] how to close all windows? ?graphics.off is for graphic windows only. You don't tell us, but I assume you are working with RGui under Windows. If you don't want many help windows, you can still change it in Edit - GUI preferences - Pager style - single window. Best, Philippe Grosjean Uwe Ligges wrote: Michael wrote: Hi all, I have 10+ graphic windows and 10+ R help window opening now... How do I close them all at once? Thanks a lot! [[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 See ?graphics.off Uwe Ligges __ R-help@stat.math.ethz.ch 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 Petr Pikal [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] power law
See the methods here:
http://arxiv.org/abs/cond-mat/0412004
You can simply use the formula given here for continuous data.
Discrete data is a bit more tricky, this is the function I'm using:
# -
power.law.fit - function(x, xmin=NULL, start=2, ...) {
if (length(x) == 0) {
error(zero length vector)
}
if (length(x) == 1) {
error(vector should be at least of length two)
}
require(stats4)
if (is.null(xmin)) { xmin - min(x) }
n - length(x)
x - x[ x = xmin]
if (length(x) != n) {
warning(too small values eliminated from vector)
n - length(x)
}
mlogl - function(alpha) {
C - 1/sum( (xmin:1)^-alpha )
-n*log(C)+alpha*sum(log(x))
}
alpha - mle(mlogl, start=list(alpha=start), ...)
alpha
}
# -
The function (with is manual page) is also included in the igraph package,
see http://cneurocvs.rmki.kfki.hu/igraph
Hope this helps,
Gabor
On Wed, Feb 15, 2006 at 08:07:50AM -0500, Glazko, Galina wrote:
Dear list,
Does anyone know how to fit the power law distribution?
I have the empirical distribution and would like to check whether it fits
power law (with the power estimated from the data).
Any hints are appreciated
Best regards
Galina
[[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
--
Csardi Gabor [EMAIL PROTECTED]MTA RMKI, ELTE TTK
__
R-help@stat.math.ethz.ch 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] var-covar matrices comparison:
My first thought is to use a random permutation test. In this setting the main question you need to ask is what distance measure do you want to use between variance matrices -- there are lots of choices. One that I've found useful is the absolute value of the maximum eigenvalue of the difference of the matrices. If you have a hypothesis about how the variances may differ, then you should be able to come up with a more powerful statistic. Patrick Burns [EMAIL PROTECTED] +44 (0)20 8525 0696 http://www.burns-stat.com (home of S Poetry and A Guide for the Unwilling S User) Aldi Kraja wrote: Hi, Using package gclus in R, I have created some graphs that show the trends within subgroups of data and correlations among 9 variables (v1-v9). Being interested for more details on these data I have produced also the var-covar matrices. Question: From a pair of two subsets of data (with 9 variables each, I have two var-covar matrices for each subgroup, that differ for a treatment on one group (treatment A) vs (non-Treatment A). Is there a software that can compare if two var-covar matrices are statistically the same? Below are a pair of two matrices, from several others. Thank you in advance for any input. Aldi First group var-covar matrix (the data were under treatment a) v1v2 v3 v4 v5 v6 v7v8 v9 v1 730.87 3.406 -283.41-74.68 107.57 -1355.13 -112.46 14.000 5.776 v2 3.41 24.950 105.45 -121.31 -307.68 -285.4029.65 -2.500 -7.796 v3 -283.41 105.451 6292.19 -2676.46 -970.80 29296.23 10715.29 3.156 -66.313 v4 -74.68-121.307 -2676.46 124492.30 -2289.47 -20377.34 -409.71 183.500 563.102 v5 107.57-307.681 -970.80 -2289.47 7045.62 12118.09 954.51 38.258 96.355 v6 -1355.13-285.404 29296.23 -20377.34 12118.09 218555.93 70126.71 137.000-130.667 v7 -112.46 29.645 10715.29 -409.71 954.51 70126.71 28239.57 67.989 -26.370 v814.00 -2.500 3.16183.50 38.26137.0067.99 24.500 9.000 v9 5.78 -7.796 -66.31563.10 96.35 -130.67 -26.37 9.000 22.776 Second group var-covar matrix (the data were NOT under treatment a) v1v2 v3 v4 v5 v6 v7v8 v9 v1 2696.25 27.05 201.06 2745.54 -344.39540.48654.20 34.363 7.623 v227.05 86.37 -96.89 -497.28 -1185.10 -3108.71 -910.38 -4.254 -9.115 v3 201.06 -96.89 10647.26 8378.07 595.81 66122.43 26237.21 -65.093 -51.998 v4 2745.54 -497.28 8378.07 408391.25 -3887.28 40477.40 30652.01 450.539 50.311 v5 -344.39-1185.10 595.81 -3887.28 29204.00 65320.00 15238.41 -98.237 102.975 v6 540.48-3108.71 66122.43 40477.40 65320.00 549955.14 194691.90-555.552 -95.210 v7 654.20 -910.38 26237.21 30652.01 15238.41 194691.90 82698.88 -70.417 -75.585 v834.36 -4.25 -65.09450.54 -98.24 -555.55-70.42 79.689 8.164 v9 7.62 -9.11 -52.00 50.31 102.97-95.21-75.58 8.164 30.492 __ R-help@stat.math.ethz.ch 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-ERROR: unknown GUI none from Statistics::R
Hi, I download Statistics::R from http://search.cpan.org/~gmpassos/Statistics-R-0.02/ I am going to combine R with perl. but it keeps give the error msg: ERROR: unknown GUI none stop running script. I feel very strange, because previously it works. my perl version is 5.8.6 my R version is Version 2.2.1 Is there anybody know about it? Thanks a lot. Nina [[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] call row names
Hi R users. I have a table like that: table var A1 A2 A3 v1 41203 3.69 2.31 v2 20577 4.51 8.60 v3 20625 2.87 3.50 v4 6115 8.92 2.97 v5 3160 1.49 2.21 v6 2954 2.62 5.98 v7 4731 1.83 7.53 v8 2435 7.68 3.50 v9 2296 3.03 4.84 v10 6153 1.06 4.28 v11 3157 1.07 1.15 v12 2996 1.06 1.01 v13 6084 2.65 2.63 v14 3115 2.42 5.70 v15 2969 2.92 7.53 * If I want column A1 I do this: table$A1 * And if I want row v4 how can I do? (probably the problem happens because the column var is not considered as row names, but in the reality was with this purpose that was created by me) Thanks in advance Ana Quiterio [[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] call row names
1) It is not good practice to call your objects after existing R functions (e.g. table) 2) I think you are getting rows and columns confused. If you want to extract the rows/column of a matrix or dataframe, then try subsetting it by mat[A1, ] or mat[ , v4]. See help(subset) for more information. 3) It looks to me that your object is a list. Try doing class(table). Regards, Adai On Tue, 2006-02-21 at 11:56 +, Ana Quitério wrote: Hi R users. I have a table like that: table var A1 A2 A3 v1 41203 3.69 2.31 v2 20577 4.51 8.60 v3 20625 2.87 3.50 v4 6115 8.92 2.97 v5 3160 1.49 2.21 v6 2954 2.62 5.98 v7 4731 1.83 7.53 v8 2435 7.68 3.50 v9 2296 3.03 4.84 v10 6153 1.06 4.28 v11 3157 1.07 1.15 v12 2996 1.06 1.01 v13 6084 2.65 2.63 v14 3115 2.42 5.70 v15 2969 2.92 7.53 * If I want column A1 I do this: table$A1 * And if I want row v4 how can I do? (probably the problem happens because the column var is not considered as row names, but in the reality was with this purpose that was created by me) Thanks in advance Ana Quiterio [[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] R-help Digest, Vol 36, Issue 21
Hello, dear R users. I've already sent a question here, but I'm not sure that it had been read. I need to visualize classification of my numerical data based on 2-3 factors. As I suppose, the best way is a tree. With an orbitrary function at the ends (leaves), or at least with means of my data at the ends. What is the way to do it? As I found, ctree offers binary classification, but it that the only way? Of course, tree is not only way, may be you could offer other ways. 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] R and packages
Agreed. Hank On Feb 21, 2006, at 2:40 AM, J Dougherty wrote: You want to take the sysad by the hairy of his chinny-chin-chin and explain the issue in short sentences. He or she ought to be able produce a solution without difficulty, and should not have any problems about doing so. It IS their job. JD On Monday 20 February 2006 17:07, Erin Hodgess wrote: Dear R People: Here is yet another strange problem. I'm using R in one of my classes. However, the computer lab has something called Deep Freeze and the students cannot save anything to the hard drive. I had R installed and things were working well. They would save their .Rdata files to disks. Now, we need to add more packages. We can't download and we can't bring them in via zip files (already tried both!) When the zip files are expanded, of course they build new directories.. I'm completely annoyed because learning to download packages and installing them from local zips are actually important tasks! We're also losing good teaching/learning time! Anyhow, what I would like to do is produce a sort of combination R installation exe with our extra libraries as part of the package. Does anyone have any suggestions on how this could be done, please? This is R for Windows, Version 2.2.1 Thanks for any help! Sorry about the totally weird problem! Sincerely, Erin Hodgess Associate Professor Department of Computer and Mathematical Sciences University of Houston - Downtown mailto: [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch 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] R-help Digest, Vol 36, Issue 21
Evgeniy Kachalin wrote: Hello, dear R users. I've already sent a question here, but I'm not sure that it had been read. I need to visualize classification of my numerical data based on 2-3 factors. As I suppose, the best way is a tree. With an orbitrary function at the ends (leaves), or at least with means of my data at the ends. What is the way to do it? As I found, ctree offers binary classification, but it that the only way? Of course, tree is not only way, may be you could offer other ways. Or the best way of it is to do it with replacement, like a 'heatmap', but with means in the cells instead of colors, if it is possible. Sorry for the second letter. __ R-help@stat.math.ethz.ch 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] Nested AIC
-Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Aaron MacNeil Sent: 20 February 2006 15:17 To: r-help@stat.math.ethz.ch Subject: [R] Nested AIC Greetings, I have recently come into some confusion over weather or not AIC results for comparing among models requires that they be nested. Reading Burnham Anderson (2002) they are explicit that nested models are not required, but other respected statisticians have suggested that nesting is a pre-requisite for comparison. Could anyone who feels strongly regarding either position post their arguments for or against nested models and AIC? This would assist me greatly in some analysis I am currently conducting. Many thanks, Aaron Hi, Aaron, Burnham Anderson are explicit but they do not go into any depth regarding this issue. Akaike's colleagues Sakamoto, Ishiguro, and Kitagawa (Akaike Information Criterion Statistics, 1986, KTK Scientific Publishers) do no either, deal with it directly, and the examples they present that I have examined (not even half of the total in the book), are all of nested models. However, by reading some of Akaike's papers and the book quoted above it does not appear to me that there is any restriction on the use of the AIC related to nestedness. In fact, the theory does not preclude the comparison of models with different *probability densities (or mass)* as long as you keep all constants (like 1/sqrt(2pi) in the normal) in the calculation. Akaike (1973) wrote in the first sentence of his paper his general principle, which he called an extension of the maximum likelihood principle: Given a set of estimates theta_hat's of the vector of parameters theta of a probability distribution with density f(x|theta) we adopt as our final estimate the one which will give the maximum of the expected log-likelihood, which is by definition E(log f(X|theta_hat))=E(INTEGRAL f(x|theta)log f(x|theta_hat)dx) Where X is a random variable following the distribution with the density function f(x|theta) and is independent of theta_hat. All subsequent derivations in the paper, like the choice of distance measure, class of estimates, and elimination of the true parameter value, revolve around this principle. Now, nestedness is a mathematical property of what Burnham Anderson call the structural model, whereas Akaike's principle only concerns the probabilistic model f(x|theta) where the structural model is embedded. I reply to you even though I do not feel strongly about this issue and you asked for replies from people who feel strongly about this issue. Ruben __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] help-ERROR: unknown GUI none from Statistics::R
On Tue, 21 Feb 2006, jia ding wrote: Hi, I download Statistics::R from http://search.cpan.org/~gmpassos/Statistics-R-0.02/ I am going to combine R with perl. but it keeps give the error msg: ERROR: unknown GUI none stop running script. I feel very strange, because previously it works. my perl version is 5.8.6 my R version is Version 2.2.1 Is there anybody know about it? Yes. R --gui=none has not been needed for a long time, and it no longer works (i.e. has been removed). You need to remove it in lib/Statistics/R/Bridge/Linux.pm. 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
[R] glmm.admb - bug and possible solution??
Dear R list, A problem with the package glmmADMB was pointed out by R. Bagchi some time ago. We have now fixed the problem and a new version of the package can be downloaded from the usual place: http://otter-rsch.com/admbre/examples/glmmadmb/glmmADMB.html The likelihood values that previously were off by a constant have now been corrected. Please let me know if there are other problems. hans Dear Dr Skaug and R users, just discovered glmm.admb in R, and it seems a very useful tool. However, I ran into a problem when I compare two models: m1-glmm.admb(survival~light*species*damage, random=~1, group=table, data=bm, family=binomial, link=logit) m1.1-glmm.admb(survival~(light+species+damage)^2, random=~1, group=table, data=bm, family=binomial, link=logit) anova(m1, m1.1) I get the following output with the warning Analysis of Variance Table Model 1: survival ~ light * species * damage Model 2: survival ~ (light + damage + species)^2 NoPar LogLik Df -2logQ P.value 1 9.000 -103.307 2 9.000 -103.781 0 -0.948 Warning message: NaNs produced in: pchisq(q, df, lower.tail, log.p) _ Hans Julius Skaug Department of Mathematics University of Bergen Johannes Brunsgate 12 5008 Bergen Norway ph. (+47) 55 58 48 61 __ R-help@stat.math.ethz.ch 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 ROC for SVM
Message: 1 Date: Mon, 20 Feb 2006 17:19:41 +0530 From: karthi keyan [EMAIL PROTECTED] Subject: [R] how to get ROC for SVM To: r-help@stat.math.ethz.ch Message-ID: [EMAIL PROTECTED] Content-Type: text/plain Hi How can i get ROC for svm during modelling in SVM. how to plot the values thanks karthik Probably this will help: http://rocr.bioinf.mpi-sb.mpg.de/ Paolo __ R-help@stat.math.ethz.ch 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] rotated labels in barplot with beside=T and multiple groups
I have a data set that I display using barplot. I don't know what you call it, but when I look at it, it looks like this: lsu (0,0.1] (0.1,0.2] (0.2,0.3] (0.3,0.4] (0.4,0.5] (0.5,0.6] A 0.052631579 0.0 0.0 0.0 0.0 0.0 B 0.0 0.0 0.001007049 0.003021148 0.0 0.0 E 0.2 0.0 0.0 0.0 0.1 0.0 (0.6,0.7] (0.7,0.8] (0.8,0.9] (0.9,1] A 0.0 0.0 0.0 0.947368421 B 0.0 0.004028197 0.005035247 0.986908359 E 0.1 0.0 0.1 0.5 Now, trying the examples shown via the r-help mailing list I am trying to make a plot where each of the groups gets displayed in a histogram-like fashion upwards with the number 0.1, 0.2 and so forth underneath the group. What I do is the following: par(mar = c(6, 4, 4, 2) + 0.1) bplot = barplot(lsu, beside=TRUE, col=colors[1:length(lsu[,1])], ylim = c(0,1.0), xaxt = n, xlab = ) axis(side=1,at=bplot, labels=FALSE, tick=TRUE) NULL nam=rep(a,10) text(bplot, par(usr)[3] - 1.5, srt = 45, adj = 1, labels = nam, xpd = TRUE) NULL The result is the bars pointing upwards, like I want, but I get one tickmark per bar, and no labels underneath. I want no tickmark, and one label per group. Any ideas as to what I am doing wrong? TIA, Karin -- Karin Lagesen, PhD student [EMAIL PROTECTED] http://www.cmbn.no/rognes/ __ R-help@stat.math.ethz.ch 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] rotated labels in barplot with beside=T and multiple groups
why not: colnames(lsu) - seq(0.1,1,by=0.1) barplot(lsu, bes=T) Karin Lagesen a écrit : I have a data set that I display using barplot. I don't know what you call it, but when I look at it, it looks like this: lsu (0,0.1] (0.1,0.2] (0.2,0.3] (0.3,0.4] (0.4,0.5] (0.5,0.6] A 0.052631579 0.0 0.0 0.0 0.0 0.0 B 0.0 0.0 0.001007049 0.003021148 0.0 0.0 E 0.2 0.0 0.0 0.0 0.1 0.0 (0.6,0.7] (0.7,0.8] (0.8,0.9] (0.9,1] A 0.0 0.0 0.0 0.947368421 B 0.0 0.004028197 0.005035247 0.986908359 E 0.1 0.0 0.1 0.5 Now, trying the examples shown via the r-help mailing list I am trying to make a plot where each of the groups gets displayed in a histogram-like fashion upwards with the number 0.1, 0.2 and so forth underneath the group. What I do is the following: par(mar = c(6, 4, 4, 2) + 0.1) bplot = barplot(lsu, beside=TRUE, col=colors[1:length(lsu[,1])], ylim = c(0,1.0), xaxt = n, xlab = ) axis(side=1,at=bplot, labels=FALSE, tick=TRUE) NULL nam=rep(a,10) text(bplot, par(usr)[3] - 1.5, srt = 45, adj = 1, labels = nam, xpd = TRUE) NULL The result is the bars pointing upwards, like I want, but I get one tickmark per bar, and no labels underneath. I want no tickmark, and one label per group. Any ideas as to what I am doing wrong? TIA, Karin __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] How to Import Data
Hello, I am a very new user of R. I've spent several hours trying to import data, so I feel okay asking the list for help. I had an Excel file, then I turned it into a csv file, as instructed by directions. My filename is x111.csv. I then used the following commands to read this (fairly small) dataset in. x111 -read.table(file='x111.csv', sep=,header=T, quote=,comment.char=,as.is=T) I then get the following error message. Error in file(file, r) : unable to open connection In addition: Warning message: cannot open file 'x111.csv', reason 'No such file or directory' I would imagine I'm not putting my csv file in the right location for R to be able to read it. If that's the case, where should I put it? Or is there something else I need to do to it first? Thanks for your help, Carl __ R-help@stat.math.ethz.ch 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] rotated labels in barplot with beside=T and multiple groups
On Tue, 2006-02-21 at 14:07 +0100, Karin Lagesen wrote: I have a data set that I display using barplot. I don't know what you call it, but when I look at it, it looks like this: lsu (0,0.1] (0.1,0.2] (0.2,0.3] (0.3,0.4] (0.4,0.5] (0.5,0.6] A 0.052631579 0.0 0.0 0.0 0.0 0.0 B 0.0 0.0 0.001007049 0.003021148 0.0 0.0 E 0.2 0.0 0.0 0.0 0.1 0.0 (0.6,0.7] (0.7,0.8] (0.8,0.9] (0.9,1] A 0.0 0.0 0.0 0.947368421 B 0.0 0.004028197 0.005035247 0.986908359 E 0.1 0.0 0.1 0.5 It appears to be a matrix, comprised of the results of using cut() along with perhaps tapply() or similar on one or more continuous variables using a grouping variable. See ?cut and ?tapply. Now, trying the examples shown via the r-help mailing list I am trying to make a plot where each of the groups gets displayed in a histogram-like fashion upwards with the number 0.1, 0.2 and so forth underneath the group. What I do is the following: par(mar = c(6, 4, 4, 2) + 0.1) bplot = barplot(lsu, beside=TRUE, col=colors[1:length(lsu[,1])], ylim = c(0,1.0), xaxt = n, xlab = ) axis(side=1,at=bplot, labels=FALSE, tick=TRUE) NULL nam=rep(a,10) text(bplot, par(usr)[3] - 1.5, srt = 45, adj = 1, labels = nam, xpd = TRUE) NULL The result is the bars pointing upwards, like I want, but I get one tickmark per bar, and no labels underneath. I want no tickmark, and one label per group. Any ideas as to what I am doing wrong? TIA, Karin First, your code above has some errors. The use of col=colors[1:length(lsu[,1])] does not work. It can either be: col=colors()[1:length(lsu[,1])] or easier to read: colors()[1:nrow(lsu)] Second, the offset (-1.5) you have in the call to text(): par(usr)[3] - 1.5 puts the text labels well below the bottom of the plot, which is why they are not seen. You need to change it to: par(usr)[3] - .05 or a similarly reduced offset figure. Thus, you can now use: par(mar = c(6, 4, 4, 2) + 0.1) bplot - barplot(lsu, beside = TRUE, col = colors()[1:nrow(lsu)], ylim = c(0, 1.0), xaxt = n, xlab = ) axis(side = 1, at = bplot, labels = FALSE, tick = TRUE) nam - rep(a, 10) text(bplot, par(usr)[3] - .05, srt = 45, adj = 1, labels = nam, xpd = TRUE) In order to get group labels for each of the three bars, you need to pay attention to the Value section of ?barplot, which says: If 'beside' is true, use 'colMeans(mp)' for the midpoints of each _group_ of bars, see example. Thus, if you want to place group labels on the plot use: par(mar = c(6, 4, 4, 2) + 0.1) bplot - barplot(lsu, beside = TRUE, col = colors()[1:nrow(lsu)], ylim = c(0, 1.0), xaxt = n, xlab = ) nam - rep(a, 10) # Here is the change. Use colMeans(bplot) text(colMeans(bplot), par(usr)[3] - .05, srt = 45, adj = 1, labels = nam, xpd = TRUE) Do not use the call to axis() if you don't want the tick marks, since there is no other need for that function here. If you do want to put tick marks at the midpoint of each group, you could use: axis(1, at = colMeans(bplot), labels = FALSE) HTH, Marc Schwartz __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] visualise classification by factors (was Re: R-help Digest, Vol 36, Issue 21)
1) Please use a meaning subject line. Start a new thread instead of replying to another thread. 2) Please give a simple example (if possible reproducible) to help explain the problem. 3) Please read the posting guide. On Tue, 2006-02-21 at 15:12 +0300, Evgeniy Kachalin wrote: Hello, dear R users. I've already sent a question here, but I'm not sure that it had been read. I need to visualize classification of my numerical data based on 2-3 factors. As I suppose, the best way is a tree. With an orbitrary function at the ends (leaves), or at least with means of my data at the ends. What is the way to do it? As I found, ctree offers binary classification, but it that the only way? Of course, tree is not only way, may be you could offer other ways. 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 __ R-help@stat.math.ethz.ch 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 Import Data
carl, you wrote: I would imagine I'm not putting my csv file in the right location for R to be able to read it. If that's the case, where should I put it? Or is there something else I need to do to it first? getwd() gives you the working directory in which the datafile has to reside. johann __ R-help@stat.math.ethz.ch 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] wireframe, axis label-axis separation, xlab rotation
Greetings, A couple of questions: 1. I am using wireframe. It prints 3d plots nicely on screen, but when I use a postscript device the z-axis label merges with the z-axis. Is there any option to control their separation? 2. Again, using wireframe, I rotate the whole plot using the screen parameter. However, the x- and y-axis labels (xlab and ylab) remain horizontal. How could I rotate them? Thanks in advance, Joaquin PS: my apologizes if these questions are too basic. I am new to R and to Trellis. -- Joaquin Rapela PhD Student, Electrical Engineering University of Southern California Hedco Neuroscience Building 3641 Watt Way Los Angeles, CA 90089-2520 tel/fax: (213) 821-2070 http://www-scf.usc.edu/~rapela -- This is the origin of their (the Aborigines') belief in the omnipotence of thoughts, their unshakable confidence in the possibility of controlling the world and their inaccessibility to other experiences, so easily obtainable, which could teach them man's true position in the universe. Sigmund Freud Totem and Taboo __ R-help@stat.math.ethz.ch 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 Import Data
select the directory with setwd() and then import data: setwd(d:/.../yourdirectory) x111 - read.table(x111.csv,...) or indicate path behind filename: x111 - read.table(d:/.../yourdirectory/x111.csv,...) besides, there are other functions to import data. see ?read.table Carl Klarner a écrit : Hello, I am a very new user of R. I've spent several hours trying to import data, so I feel okay asking the list for help. I had an Excel file, then I turned it into a csv file, as instructed by directions. My filename is x111.csv. I then used the following commands to read this (fairly small) dataset in. x111 -read.table(file='x111.csv', sep=,header=T, quote=,comment.char=,as.is=T) I then get the following error message. Error in file(file, r) : unable to open connection In addition: Warning message: cannot open file 'x111.csv', reason 'No such file or directory' I would imagine I'm not putting my csv file in the right location for R to be able to read it. If that's the case, where should I put it? Or is there something else I need to do to it first? Thanks for your help, Carl __ R-help@stat.math.ethz.ch 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 Import Data
Carl Klarner cklarner at isugw.indstate.edu writes: I would imagine I'm not putting my csv file in the right location for R to be able to read it. If that's the case, where should I put it? Or is there something else I need to do to it first? Thanks for your help, Carl You're probably right. R opens by default in its installation directory (usually somewhere in Program Files\R\... on Windows) The easiest thing is to change the working directory to wherever you have the files using setwd() or (in a GUI) some menu option (Change dir in the File menu under Windows). Useful functions for working directory, file selection etc.: list.files() [list files in current directory], file.choose() [interactive file chooser], file.exists() [whether or not a file exists], getwd() [print current working directory], setwd() [set current working directory]. The other common problem, which you probably *aren't* having, is hidden file extensions under Windows. good luck Ben Bolker __ R-help@stat.math.ethz.ch 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] Resolution of plots?
Hello. I have a problem regarding the output from the R plot window. I have a quite dense dendrogram that I wish to visualize using a suitable software (like Illustrator), and to accomplish that I right click on the image from the plotting device and I save as metafile. The problem is when I scale up the image, some of the vectors in the dendrogram image are inseparable. My question is if there is any way of increasing the resolution of the meta-file that I get from the plotting device? Thanks for any help. /Marcus Marcus Gry Björklund Royal Institute of Technology AlbaNova University Center Department of Molecular Biotechnology 106 91 Stockholm, Sweden Phone (office): +46 8 553 783 44 Fax: + 46 8 553 784 81 Visiting address: Roslagstullsbacken 21, Floor 3 Delivery address: Roslagsvägen 30B Web: http://www.biotech.kth.se/molbio/microarray/index.html [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to Import Data
1) You need to use sep=, which is appropriate for a CSV file. 2) You need to specify the FULL path to the file. See http://cran.r-project.org/bin/windows/base/rw-FAQ.html#R-can_0027t-find-my-file 3) You can use read.csv which is the read.table variant for CSV files. For example a - read.csv( file=c:/Progra~1/Docume~1/ramasamy/x111.csv ) might work if you replace it with your full path. If you have the _unique_ rownames in the first column, you can add the argument row.names=1 in the call. Regards, Adai On Tue, 2006-02-21 at 08:52 -0500, Carl Klarner wrote: Hello, I am a very new user of R. I've spent several hours trying to import data, so I feel okay asking the list for help. I had an Excel file, then I turned it into a csv file, as instructed by directions. My filename is x111.csv. I then used the following commands to read this (fairly small) dataset in. x111 -read.table(file='x111.csv', sep=,header=T, quote=,comment.char=,as.is=T) I then get the following error message. Error in file(file, r) : unable to open connection In addition: Warning message: cannot open file 'x111.csv', reason 'No such file or directory' I would imagine I'm not putting my csv file in the right location for R to be able to read it. If that's the case, where should I put it? Or is there something else I need to do to it first? Thanks for your help, Carl __ R-help@stat.math.ethz.ch 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 Import Data
Or, you should put the whole path to that file, for example file='C:/My documents/x111.csv' Best, Ales Ziberna Johann Jacoby pravi: carl, you wrote: I would imagine I'm not putting my csv file in the right location for R to be able to read it. If that's the case, where should I put it? Or is there something else I need to do to it first? getwd() gives you the working directory in which the datafile has to reside. johann __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to Import Data
On Tue, Feb 21, 2006 at 08:52:04AM -0500, Carl Klarner wrote: x111 -read.table(file='x111.csv', sep=,header=T, quote=,comment.char=,as.is=T) to make things easier for you you could do x111 -read.table( file.choose(), ... ) Also: Why are you setting sep= ? Your filename suggests that you have some kind of separator in your file... cu Philipp -- Dr. Philipp PagelTel. +49-8161-71 2131 Dept. of Genome Oriented Bioinformatics Fax. +49-8161-71 2186 Technical University of Munich Science Center Weihenstephan 85350 Freising, Germany and Institute for Bioinformatics / MIPS Tel. +49-89-3187 3675 GSF - National Research Center Fax. +49-89-3187 3585 for Environment and Health Ingolstädter Landstrasse 1 85764 Neuherberg, Germany http://mips.gsf.de/staff/pagel __ R-help@stat.math.ethz.ch 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 Import Data
On 2/21/2006 9:13 AM, Adaikalavan Ramasamy wrote:
1) You need to use sep=, which is appropriate for a CSV file.
2) You need to specify the FULL path to the file. See
http://cran.r-project.org/bin/windows/base/rw-FAQ.html#R-can_0027t-find-my-file
The advice to look at the FAQ is good, but it's not true that you need
the full path. Like other programs, R maintains a current working
directory, and paths can be specified relative to that. I usually put a
line like
setwd('c:/my/working/dir')
at the start of scripts, so that the current directory is changed from
the beginning. Then I can use simple filenames to read files.
The other advice I usually give is to specify file.choose() instead of
an explicit filename; this will open a file selection dialog. (In my
opinion, this should be the default, but some people disagree quite
strongly.)
Duncan Murdoch
3) You can use read.csv which is the read.table variant for CSV files.
For example
a - read.csv( file=c:/Progra~1/Docume~1/ramasamy/x111.csv )
might work if you replace it with your full path. If you have the
_unique_ rownames in the first column, you can add the argument
row.names=1 in the call.
Regards, Adai
On Tue, 2006-02-21 at 08:52 -0500, Carl Klarner wrote:
Hello,
I am a very new user of R. I've spent several hours trying to import
data, so I feel okay asking the list for help. I had an Excel file,
then I turned it into a csv file, as instructed by directions. My
filename is x111.csv. I then used the following commands to read this
(fairly small) dataset in.
x111 -read.table(file='x111.csv',
sep=,header=T,
quote=,comment.char=,as.is=T)
I then get the following error message.
Error in file(file, r) : unable to open connection
In addition: Warning message:
cannot open file 'x111.csv', reason 'No such file or directory'
I would imagine I'm not putting my csv file in the right location for R
to be able to read it. If that's the case, where should I put it? Or
is there something else I need to do to it first?
Thanks for your help,
Carl
__
R-help@stat.math.ethz.ch 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] Resolution of plots?
On Tue, 21 Feb 2006, Marcus Gry Björklund wrote: I have a problem regarding the output from the R plot window. It seems you are using Windows, but you didn't say so! I will assume that from here on down. I have a quite dense dendrogram that I wish to visualize using a suitable software (like Illustrator), and to accomplish that I right click on the image from the plotting device and I save as metafile. The problem is when I scale up the image, some of the vectors in the dendrogram image are inseparable. My question is if there is any way of increasing the resolution of the meta-file that I get from the plotting device? No, as it is done to full metafile resolution (1/1200). Why not plot on a larger size directly to a metafile using win.metafile? -- 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] var-covar matrices comparison:
Thank you Patrick for your thoughts. I was expecting someone had written a function that does maybe a hierarchical comparison of var-covar matrices using the flury hierarchy as proposed by Philips and Arnold in Evolution 53(5), 1999;1506-1515. But your email has an interesting idea, so I will try to test it. Thank you, Aldi Patrick Burns wrote: My first thought is to use a random permutation test. In this setting the main question you need to ask is what distance measure do you want to use between variance matrices -- there are lots of choices. One that I've found useful is the absolute value of the maximum eigenvalue of the difference of the matrices. If you have a hypothesis about how the variances may differ, then you should be able to come up with a more powerful statistic. Patrick Burns [EMAIL PROTECTED] +44 (0)20 8525 0696 http://www.burns-stat.com (home of S Poetry and A Guide for the Unwilling S User) __ R-help@stat.math.ethz.ch 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] wireframe, axis label-axis separation, xlab rotation
On 2/21/06, Joaquin Rapela [EMAIL PROTECTED] wrote: Greetings, A couple of questions: 1. I am using wireframe. It prints 3d plots nicely on screen, but when I use a postscript device the z-axis label merges with the z-axis. Is there any option to control their separation? 2. Again, using wireframe, I rotate the whole plot using the screen parameter. However, the x- and y-axis labels (xlab and ylab) remain horizontal. How could I rotate them? Both are possible, e.g. wireframe(volcano, scales = list(z = list(distance = 3)), zlab = list(volcano, rot = 90), zoom = 0.4) (The documentation is a bit lacking of these details.) Deepayan __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Number of Days Between Dates: Incorrect Results For Date Calucations.
In some cases, incorrect results are produced by the code below intended to calculate the number of days between 2 dates. The year in question was a leap year. Note the results for 2004-04-04 and 2004-04-05 are the same! They should be 37 and 38 respectively. as.integer(as.POSIXct(2004-04-02) - as.POSIXct(2004-02-27)) [1] 35 as.integer(as.POSIXct(2004-04-03) - as.POSIXct(2004-02-27)) [1] 36 as.integer(as.POSIXct(2004-04-04) - as.POSIXct(2004-02-27)) [1] 37 as.integer(as.POSIXct(2004-04-05) - as.POSIXct(2004-02-27)) [1] 37 as.integer(as.POSIXct(2004-04-06) - as.POSIXct(2004-02-27)) [1] 38 as.integer(difftime(as.POSIXct(2004-04-06), as.POSIXct(2004-02-27),units=days)) [1] 38 as.integer(difftime(as.POSIXct(2004-04-04), as.POSIXct(2004-02-27),units=days)) [1] 37 as.integer(difftime(as.POSIXct(2004-04-05), as.POSIXct(2004-02-27),units=days)) [1] 37 It appears that difftime() and - are producing invalid results. Regards, Gerald Herbert __ R-help@stat.math.ethz.ch 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] mva.pairs
Hi Vijay, Vijay A Raghavan wrote: Hello, I am using the following code to plot an MVA plot. library(affy) library(Biobase) library(limma) library(gcrma) pd-read.phenoData(Clk.targets.2.txt,header=TRUE, row.names=1,as.is=TRUE,sep=\t) Data - ReadAffy(filenames=pData(pd)$FileName,phenoData=pd) Print(Data) eset - gcrma(Data) write.exprs(eset, file=clk.6-23-05.txt) bitmap(clk-1.mva.jpg,width=15, height=12) mva.pairs(exprs(eset)[,c(19,20,21,22,23,24)],log.it=FALSE) dev.off() but I am getting an subscript out of bound error. Any Idea ? The only reason I know for getting a subscript out of bound error is actually trying to subset an object using dimensions larger than the object itself. In other words, do you really have 24 columns (e.g., 24 chips) in your exprSet? Thanks, Vijay [[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 -- James W. MacDonald, M.S. Biostatistician Affymetrix and cDNA Microarray Core University of Michigan Cancer Center 1500 E. Medical Center Drive 7410 CCGC Ann Arbor MI 48109 734-647-5623 __ R-help@stat.math.ethz.ch 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 and packages
On Mon, 20 Feb 2006, Erin Hodgess wrote: Dear R People: Here is yet another strange problem. I'm using R in one of my classes. However, the computer lab has something called Deep Freeze and the students cannot save anything to the hard drive. I had R installed and things were working well. They would save their .Rdata files to disks. Now, we need to add more packages. We can't download and we can't bring them in via zip files (already tried both!) When the zip files are expanded, of course they build new directories.. I would suggest making this the system administrator's problem. Either your site really has a policy that would prohibit R packages (in which case trying to subvert it may not be a good idea) or it is exactly the sort of problem that system administrators are supposed to solve. Distrust for R binary (and even source) packages is not a completely unreasonable attitude. They can execute arbitrary code on your computer and, in contrast to R itself, you usually have no real knowledge of who wrote them. -thomas I'm completely annoyed because learning to download packages and installing them from local zips are actually important tasks! We're also losing good teaching/learning time! Anyhow, what I would like to do is produce a sort of combination R installation exe with our extra libraries as part of the package. Does anyone have any suggestions on how this could be done, please? This is R for Windows, Version 2.2.1 Thanks for any help! Sorry about the totally weird problem! Sincerely, Erin Hodgess Associate Professor Department of Computer and Mathematical Sciences University of Houston - Downtown mailto: [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 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] indexing within panels in xyplot
On 2/21/06, Sebastian Luque [EMAIL PROTECTED] wrote:
Dear R-helpers,
I need to show a linear fit through a subset of the data within each
combination of levels of two factors. So I prepared an xyplot with
different panels for each level of one of the factors, and different
symbols within each panel for the levels of the second factor. My problem
is selecting the subset of each combination through which the line should
be fit for subsequent plotting. This hopefully shows the idea:
--cut here---start-
toydf - expand.grid(1:100, c(A, B),
c(pop1, pop2, pop3, pop4, pop5))
toydf - data.frame(facA = toydf[[3]], facB = toydf[[2]],
x = toydf[[1]], y = rnorm(1000))
xyplot(y ~ x | facA, groups = facB, data = toydf,
panel.groups = function(x, y, subscripts, ...) {
panel.xyplot(x, y, ...)
lindx - which(y[subscripts] == max(y[subscripts], na.rm = TRUE))
xleft - mean(x[lindx], na.rm = TRUE)
fit - lm(y[x = xleft] ~ x[x = xleft])
panel.abline(fit)
})
--cut here---end---
i.e. the left limit for fitting the line is defined by the mean of x
values where y is equal to the maximum y values, *within* each combination
of levels of both factors. The above is giving me:
Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
0 (non-NA) cases
In addition: Warning message:
no finite arguments to max; returning -Inf
which shows I'm not understanding how the 'subscripts' argument works.
I'd appreciate some pointers on what I'm doing wrong, as I haven't been
able to find help in the help pages and List archives.
Well, there are exceptions to this rule, but generally x and y, when
they are passed on to the panel function, are _already_ subsetted, so
x[subscripts] makes absolutely no sense. Note how your panel function
calls
panel.xyplot(x, y, ...)
without referring to subscripts at all. The subscripts argument is
there for other variables (e.g. if you were drawing confidence
intervals, and had a separate vector in your data specifying the
interval lengths). In your case, there are no other variables
involved, so just get rid of the subscripts.
Deepayan
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Number of Days Between Dates: Incorrect Results For Date Calucations.
The results are actually correct if you consider daylight savings time. For example, try this and note that the difference is 23 hours, not 24 hours: as.POSIXct(2004-04-05) - as.POSIXct(2004-04-04) You can address this by either using Date or chron classes or adding the tz = GMT argument on your as.POSIXct calls as GMT does not have daylight savings time. See the Help Desk article in R News 4/1 for more on this. On 2/21/06, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: In some cases, incorrect results are produced by the code below intended to calculate the number of days between 2 dates. The year in question was a leap year. Note the results for 2004-04-04 and 2004-04-05 are the same! They should be 37 and 38 respectively. as.integer(as.POSIXct(2004-04-02) - as.POSIXct(2004-02-27)) [1] 35 as.integer(as.POSIXct(2004-04-03) - as.POSIXct(2004-02-27)) [1] 36 as.integer(as.POSIXct(2004-04-04) - as.POSIXct(2004-02-27)) [1] 37 as.integer(as.POSIXct(2004-04-05) - as.POSIXct(2004-02-27)) [1] 37 as.integer(as.POSIXct(2004-04-06) - as.POSIXct(2004-02-27)) [1] 38 as.integer(difftime(as.POSIXct(2004-04-06), as.POSIXct(2004-02-27),units=days)) [1] 38 as.integer(difftime(as.POSIXct(2004-04-04), as.POSIXct(2004-02-27),units=days)) [1] 37 as.integer(difftime(as.POSIXct(2004-04-05), as.POSIXct(2004-02-27),units=days)) [1] 37 It appears that difftime() and - are producing invalid results. Regards, Gerald Herbert __ R-help@stat.math.ethz.ch 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] Number of Days Between Dates: Incorrect Results For DateCalucations.
check again your results with as.integer() replaced by round(); check also ?as.integer for its usage. 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://www.med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm - Original Message - From: [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Tuesday, February 21, 2006 4:48 PM Subject: [R] Number of Days Between Dates: Incorrect Results For DateCalucations. In some cases, incorrect results are produced by the code below intended to calculate the number of days between 2 dates. The year in question was a leap year. Note the results for 2004-04-04 and 2004-04-05 are the same! They should be 37 and 38 respectively. as.integer(as.POSIXct(2004-04-02) - as.POSIXct(2004-02-27)) [1] 35 as.integer(as.POSIXct(2004-04-03) - as.POSIXct(2004-02-27)) [1] 36 as.integer(as.POSIXct(2004-04-04) - as.POSIXct(2004-02-27)) [1] 37 as.integer(as.POSIXct(2004-04-05) - as.POSIXct(2004-02-27)) [1] 37 as.integer(as.POSIXct(2004-04-06) - as.POSIXct(2004-02-27)) [1] 38 as.integer(difftime(as.POSIXct(2004-04-06), as.POSIXct(2004-02-27),units=days)) [1] 38 as.integer(difftime(as.POSIXct(2004-04-04), as.POSIXct(2004-02-27),units=days)) [1] 37 as.integer(difftime(as.POSIXct(2004-04-05), as.POSIXct(2004-02-27),units=days)) [1] 37 It appears that difftime() and - are producing invalid results. Regards, Gerald Herbert __ R-help@stat.math.ethz.ch 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] Nested AIC
This might be a more suitable message for eg the stats-discuss mailing list or one of the sci.stat.* newsgroups. It is more complicated that it looks, partly because of the Anna Karenina problem: all nested models are the same, but non-nested models can be non-nested in different ways Some notes: 1) Sometimes the AIC is clearly inappropriate: eg comparing the fit of a Poisson regression to a least-squares linear regression for count data. Here the likelihoods are not densities with respect to the same measure, so the likelihood ratio is meaningless. You could also argue that the linear model isn't really being fitted by maximum likelihood. 2) You need to be careful when fitting models with different R functions, since they may omit different constants in the likelihood. 3) Transformations of the outcome are a problem. You can frame this as a mathematical problem or just note the difficulty of saying what you mean when you decide that the multiplicative error in one model is smaller than the additive error in another model. 4) If you have two least-squares linear regression models with the same outcome variable and different predictors then the AIC is choosing based on a consistent estimate of the mean squared prediction error, and in that sense it is a valid way to choose the model that predicts best. This may or may not be the criterion you want, but if it isn't what you want then AIC isn't going to help. 5) If you have a large number of models then (nested or not) there is no guarantee that the estimate of prediction error is *uniformly* consistent, so the arguments behind AIC do not necessarily work. -thomas On Tue, 21 Feb 2006, Ruben Roa wrote: -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Aaron MacNeil Sent: 20 February 2006 15:17 To: r-help@stat.math.ethz.ch Subject: [R] Nested AIC Greetings, I have recently come into some confusion over weather or not AIC results for comparing among models requires that they be nested. Reading Burnham Anderson (2002) they are explicit that nested models are not required, but other respected statisticians have suggested that nesting is a pre-requisite for comparison. Could anyone who feels strongly regarding either position post their arguments for or against nested models and AIC? This would assist me greatly in some analysis I am currently conducting. Many thanks, Aaron Hi, Aaron, Burnham Anderson are explicit but they do not go into any depth regarding this issue. Akaike's colleagues Sakamoto, Ishiguro, and Kitagawa (Akaike Information Criterion Statistics, 1986, KTK Scientific Publishers) do no either, deal with it directly, and the examples they present that I have examined (not even half of the total in the book), are all of nested models. However, by reading some of Akaike's papers and the book quoted above it does not appear to me that there is any restriction on the use of the AIC related to nestedness. In fact, the theory does not preclude the comparison of models with different *probability densities (or mass)* as long as you keep all constants (like 1/sqrt(2pi) in the normal) in the calculation. Akaike (1973) wrote in the first sentence of his paper his general principle, which he called an extension of the maximum likelihood principle: Given a set of estimates theta_hat's of the vector of parameters theta of a probability distribution with density f(x|theta) we adopt as our final estimate the one which will give the maximum of the expected log-likelihood, which is by definition E(log f(X|theta_hat))=E(INTEGRAL f(x|theta)log f(x|theta_hat)dx) Where X is a random variable following the distribution with the density function f(x|theta) and is independent of theta_hat. All subsequent derivations in the paper, like the choice of distance measure, class of estimates, and elimination of the true parameter value, revolve around this principle. Now, nestedness is a mathematical property of what Burnham Anderson call the structural model, whereas Akaike's principle only concerns the probabilistic model f(x|theta) where the structural model is embedded. I reply to you even though I do not feel strongly about this issue and you asked for replies from people who feel strongly about this issue. Ruben __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html 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] Number of Days Between Dates: Incorrect Results For Date Calucations.
[EMAIL PROTECTED] writes: In some cases, incorrect results are produced by the code below intended to calculate the number of days between 2 dates. The year in question was a leap year. Note the results for 2004-04-04 and 2004-04-05 are the same! They should be 37 and 38 respectively. Nope. First, it depends on your timezone. Over here, they do actually differ. However, a few weeks earlier, we have the similar phenomenon as.POSIXct(2004-03-29) - as.POSIXct(2004-02-27) Time difference of 30.95833 days as.POSIXct(2004-03-28) - as.POSIXct(2004-02-27) Time difference of 30 days which is of course because March 28 was only 23 hours long which is in turn because as.POSIXct(2004-03-28) [1] 2004-03-28 CET as.POSIXct(2004-03-29) [1] 2004-03-29 CEST Get it? as.integer(as.POSIXct(2004-04-02) - as.POSIXct(2004-02-27)) [1] 35 as.integer(as.POSIXct(2004-04-03) - as.POSIXct(2004-02-27)) [1] 36 as.integer(as.POSIXct(2004-04-04) - as.POSIXct(2004-02-27)) [1] 37 as.integer(as.POSIXct(2004-04-05) - as.POSIXct(2004-02-27)) [1] 37 as.integer(as.POSIXct(2004-04-06) - as.POSIXct(2004-02-27)) [1] 38 as.integer(difftime(as.POSIXct(2004-04-06), as.POSIXct(2004-02-27),units=days)) [1] 38 as.integer(difftime(as.POSIXct(2004-04-04), as.POSIXct(2004-02-27),units=days)) [1] 37 as.integer(difftime(as.POSIXct(2004-04-05), as.POSIXct(2004-02-27),units=days)) [1] 37 It appears that difftime() and - are producing invalid results. Regards, Gerald Herbert __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- 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] Extracting variance components from lmer
Hi Rick There may be a better way, but the following should work: attributes(vc.fit)$sc Regards, Christoph Buser -- Christoph Buser [EMAIL PROTECTED] Seminar fuer Statistik, LEO C13 ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND phone: x-41-44-632-4673 fax: 632-1228 http://stat.ethz.ch/~buser/ -- Rick DeShon writes: Hi All. I need a bit of help extracting the residual error variance from the VarCorr structure from lmer. #Here's a 2-way random effects model lmer.1- lmer(rating ~ (1|person)+(1|rater), data = dat) #Get the structure vc.fit - VarCorr(lmer.1) #results in. $person 1 x 1 Matrix of class dpoMatrix (Intercept) (Intercept) 0.7755392 $rater 1 x 1 Matrix of class dpoMatrix (Intercept) (Intercept) 0.2054469 attr(,sc) [1] 0.5051518 #I can pull out the person and rater variance components easy enough. For example... vc.person - [EMAIL PROTECTED] I'm sure it's simple but I have not been able to grab the residual variance in the last matrix. I simply wish to asign the residual to a scalar variable. Any suggestions would be appreciated! Thanks! Rick DeShon -- Rick DeShon 306 Psychology Building Department of Psychology Michigan State University East Lansing, MI 48824-1116 [[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] indexing within panels in xyplot
Deepayan Sarkar [EMAIL PROTECTED] wrote:
[...]
Well, there are exceptions to this rule, but generally x and y, when
they are passed on to the panel function, are _already_ subsetted, so
x[subscripts] makes absolutely no sense. Note how your panel function
calls
panel.xyplot(x, y, ...)
without referring to subscripts at all. The subscripts argument is
there for other variables (e.g. if you were drawing confidence
intervals, and had a separate vector in your data specifying the
interval lengths). In your case, there are no other variables
involved, so just get rid of the subscripts.
Thanks Deepayan, I was indeed quite confused about this.
I realized I needed to limit the fitted line to the range of x values the
line is fit to, so I changed to panel.curve:
--cut here---start-
xyplot(y ~ x | facA, groups = facB, data = toydf,
panel.groups = function(x, y, ...) {
panel.xyplot(x, y, ...)
lindx - which(y == max(y, na.rm = TRUE))
xleft - mean(x[lindx], na.rm = TRUE)
fit - lm(y[x = xleft] ~ x[x = xleft])
panel.curve(coef(fit)[1] + (coef(fit)[2] * x),
xleft, max(x, na.rm = TRUE))
})
--cut here---end---
but can't find a way to color the line for each group differently. I
tried passing a length-2 vector as a 'col' argument to panel.curve.
Unfortunately it's only picking the first, so that both lines get colored
the same. I'm not sure, but it seems as if I need to use
'panel.superpose' directly to do this, as the help page suggests the above
would work?
Cheers,
--
Sebastian P. Luque
__
R-help@stat.math.ethz.ch 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] color quantization / binning a variable into levels
Hi all,
I'd like to quantize a variable to map it into a limited set of integers
for use with a colormap. image and filled.contour do this mapping
inside somewhere, but I'd like to choose the colors for plotting a set of
polygons. Is there a pre-existing function that does something like this
well? i.e., is capable of using 'breaks'?
quantize-function(x,n=10, breaks=NULL){
# bin the variable x into n levels
xmin-min(x)
xmax-max(x)
1+floor(n*(x-xmin)/(xmax-xmin)*.999)
}
x- -10:10
cbind(x,quantize(x,2),quantize(x),quantize(x,21))
quantize(x,breaks=c(5,7)) #
Thanks for your time,
Dave
--
Dr. David Forrest
[EMAIL PROTECTED](804)684-7900w
[EMAIL PROTECTED] (804)642-0662h
http://maplepark.com/~drf5n/
__
R-help@stat.math.ethz.ch 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] Compute a correlation matrix from an existing covariance matrix
Dear All, I am wondering if there is an R function to convert a covariance matrix to a correlation matrix. I have a covariance matrix sigma and I want to compute the corresponding correlation matrix R from sigma. Thank you very much, Bernard - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] R and packages
From: Thomas Lumley On Mon, 20 Feb 2006, Erin Hodgess wrote: Dear R People: Here is yet another strange problem. I'm using R in one of my classes. However, the computer lab has something called Deep Freeze and the students cannot save anything to the hard drive. I had R installed and things were working well. They would save their .Rdata files to disks. Now, we need to add more packages. We can't download and we can't bring them in via zip files (already tried both!) When the zip files are expanded, of course they build new directories.. I would suggest making this the system administrator's problem. Either your site really has a policy that would prohibit R packages (in which case trying to subvert it may not be a good idea) or it is exactly the sort of problem that system administrators are supposed to solve. Distrust for R binary (and even source) packages is not a completely unreasonable attitude. They can execute arbitrary code on your computer and, in contrast to R itself, you usually have no real knowledge of who wrote them. -thomas I'm completely annoyed because learning to download packages and installing them from local zips are actually important tasks! We're also losing good teaching/learning time! Anyhow, what I would like to do is produce a sort of combination R installation exe with our extra libraries as part of the package. Does anyone have any suggestions on how this could be done, please? This is R for Windows, Version 2.2.1 In case the sysadmin refuses to make it his/her problem, this would be a good alternative. It's not hard to build a custom R installer that includes packages of your choice. It's all documented. Look in the R-admin manual (e.g., http://cran.r-project.org/doc/manuals/R-admin.html) and search for EXTRA_PKGS. Andy Thanks for any help! Sorry about the totally weird problem! Sincerely, Erin Hodgess Associate Professor Department of Computer and Mathematical Sciences University of Houston - Downtown mailto: [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 Thomas Lumley Assoc. Professor, Biostatistics [EMAIL PROTECTED] University of Washington, Seattle __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch 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] Compute a correlation matrix from an existing covariance matrix
See ?cov2cor. Andy From: Marc Bernard Dear All, I am wondering if there is an R function to convert a covariance matrix to a correlation matrix. I have a covariance matrix sigma and I want to compute the corresponding correlation matrix R from sigma. Thank you very much, Bernard - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch 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] Compute a correlation matrix from an existing covariance matrix
Marc Bernard [EMAIL PROTECTED] writes: Dear All, I am wondering if there is an R function to convert a covariance matrix to a correlation matrix. I have a covariance matrix sigma and I want to compute the corresponding correlation matrix R from sigma. You mean something like cov2cor()? (BTW, why doesn't the page for cor/var/cov/cov2cor show up on help.search(covariance) or help.search(correlation) for me??) -- 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] Compute a correlation matrix from an existing covariance matrix
Marc Bernard wrote: Dear All, I am wondering if there is an R function to convert a covariance matrix to a correlation matrix. I have a covariance matrix sigma and I want to compute the corresponding correlation matrix R from sigma. Does cov2cor do what you want? Kjetil Thank you very much, Bernard - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch 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] Nested AIC
On Tue, 21 Feb 2006, Thomas Lumley wrote: This might be a more suitable message for eg the stats-discuss mailing list or one of the sci.stat.* newsgroups. It is more complicated that it looks, partly because of the Anna Karenina problem: all nested models are the same, but non-nested models can be non-nested in different ways And I am sure Akaike appreciated that, which may be why he only (AFAIK) derived a theoretical basis for AIC under strictly limited conditions including nesting. Some notes: 1) Sometimes the AIC is clearly inappropriate: eg comparing the fit of a Poisson regression to a least-squares linear regression for count data. Here the likelihoods are not densities with respect to the same measure, so the likelihood ratio is meaningless. You could also argue that the linear model isn't really being fitted by maximum likelihood. 2) You need to be careful when fitting models with different R functions, since they may omit different constants in the likelihood. 3) Transformations of the outcome are a problem. You can frame this as a mathematical problem or just note the difficulty of saying what you mean when you decide that the multiplicative error in one model is smaller than the additive error in another model. 4) If you have two least-squares linear regression models with the same outcome variable and different predictors then the AIC is choosing based on a consistent estimate of the mean squared prediction error, and in that sense it is a valid way to choose the model that predicts best. This may or may not be the criterion you want, but if it isn't what you want then AIC isn't going to help. 5) If you have a large number of models then (nested or not) there is no guarantee that the estimate of prediction error is *uniformly* consistent, so the arguments behind AIC do not necessarily work. (That only makes sense if the model class changes with 'n', suitably defined. You do get uniform consistency over a finite class of models, one of Akaike (1973)'s conditions. However, to use AIC you don't just need a consistent estimator, but to worry about the consistency of the O(1/n) term in the mean since AIC/n is effectively s^2 + 2p/n.) One other note. AIC/n is a consistent estimator but only if the model is true, and one with a lot of sampling error. Differences in AIC are much more precisely estimated for a pair of nested models than for some non-nested pairs. So sampling error can make comparisons of AIC meaningless unless the differences are large (and 'large' grows with 'n' for some appropriate 'n'). A recent talk of mine http://www.stats.ox.ac.uk/~ripley/Nelder80.pdf may be illuminating. There is a published paper version. On Tue, 21 Feb 2006, Ruben Roa wrote: -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Aaron MacNeil Sent: 20 February 2006 15:17 To: r-help@stat.math.ethz.ch Subject: [R] Nested AIC Greetings, I have recently come into some confusion over weather or not AIC results for comparing among models requires that they be nested. Reading Burnham Anderson (2002) they are explicit that nested models are not required, but other respected statisticians have suggested that nesting is a pre-requisite for comparison. Could anyone who feels strongly regarding either position post their arguments for or against nested models and AIC? This would assist me greatly in some analysis I am currently conducting. Many thanks, Aaron Hi, Aaron, Burnham Anderson are explicit but they do not go into any depth regarding this issue. Akaike's colleagues Sakamoto, Ishiguro, and Kitagawa (Akaike Information Criterion Statistics, 1986, KTK Scientific Publishers) do no either, deal with it directly, and the examples they present that I have examined (not even half of the total in the book), are all of nested models. However, by reading some of Akaike's papers and the book quoted above it does not appear to me that there is any restriction on the use of the AIC related to nestedness. In fact, the theory does not preclude the comparison of models with different *probability densities (or mass)* as long as you keep all constants (like 1/sqrt(2pi) in the normal) in the calculation. Akaike (1973) wrote in the first sentence of his paper his general principle, which he called an extension of the maximum likelihood principle: Given a set of estimates theta_hat's of the vector of parameters theta of a probability distribution with density f(x|theta) we adopt as our final estimate the one which will give the maximum of the expected log-likelihood, which is by definition E(log f(X|theta_hat))=E(INTEGRAL f(x|theta)log f(x|theta_hat)dx) Where X is a random variable following the distribution with the density function f(x|theta) and is independent of theta_hat. All subsequent derivations in the paper, like the choice of distance measure, class
Re: [R] lattice: calling functions
On Tue, 14 Feb 2006, Deepayan Sarkar wrote: On 2/14/06, Duncan Murdoch [EMAIL PROTECTED] wrote: On 2/14/2006 9:38 AM, Gabor Grothendieck wrote: On 2/14/06, Duncan Murdoch [EMAIL PROTECTED] wrote: On 2/14/2006 8:56 AM, Wolfram Fischer wrote: I defined three functions: fun0 - function( x=1:5, y=1:5, ... ) xyplot( y ~ x, ... ) fun1 - function( x=1:5, y=1:5, ... ) fun2( y ~ x, ... ) fun2 - function( ... ) xyplot( ... ) The call of fun0() works as expected. The call of fun1() causes the following error: 'Error in eval(expr, envir, enclos) : object y not found' How should I define fun2 to avoid the error? fun2 is fine, it's fun1 that has problems. It is passing a formula through fun2 to xyplot without telling xyplot where to evaluate the arguments. If you change it to fun1 - function( x=1:5, y=1:5, ... ) fun2( y ~ x, data=enviroment(), ... ) data=NULL works too, which is apparently what lm has. The point being that the environment of the formula is looked at, but the default data=parent.frame() subverts that because of the way eval works (using enclos only when envir is a list or data frame. What's wrong with environments?). The reason for not using enclos= when envir= is an environment is that an environment already consists of a frame plus enclosing environment. The point of enclos= is to supply the enclosing environment when envir= doesn't have one built in. -thomas Even the following works: fun1 - function( x=1:5, y=1:5, ... ) fun2( y ~ x, data = data.frame(x = x), ... ) I don't understand non-standard evaluation all that well, so I'll happily consider any suggestions. I'll try changing the defaults to NULL and see if there are any obvious problems. Deepayan it will tell xyplot to look in the current environment at the time of the call, i.e. the fun1 evaluation environment where x and y live. Although this does seem to be how xyplot works, I think it indicates there is a problem with it. The help file for xyplot indicates that for the xyplot formula method the default environment is the caller environment whereas it ought to be the environment of the formula: data: For the 'formula' method, a data frame containing values for any variables in the formula, as well as 'groups' and 'subset' if applicable. By default the environment where the function was called from is used. For example, if we replace xyplot with lm it does work as expected: fun1 - function( x=1:5, y=1:5, ... ) fun2( y ~ x, ... ) fun2 - function( ... ) lm( ... ) fun1() You're right, I forgot formulas have associated environments. I've added the lattice maintainer to the cc list. 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 Thomas Lumley Assoc. Professor, Biostatistics [EMAIL PROTECTED] University of Washington, Seattle __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] How to select only certain rows when making a new dataframe?
Dear R-users, I have two data frames. The FIRST data frame has 100 rows, the SECOND data frame has only 50 rows. The data frames have different variables in columns (VAR1, VAR2, etc) but they share a column called ID that contains a unique identifer linking the two data frames. I would like to make a THIRD data frame containing just the rows of the FIRST data frame that match the rows (on ID) in the SECOND data frame. THIRD=data.frame(FIRST$ID,FIRST$VAR1,FIRST$VAR2) How can I modify this line to include in THIRD just the rows in FIRST that match the rows in SECOND? (i.e., contain the same value in the shared ID column). Thanks! Mark __ R-help@stat.math.ethz.ch 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 select only certain rows when making a new dataframe?
THIRD - merge(FIRST, SECOND, by=ID, all.x=FALSE, all.y=TRUE) [EMAIL PROTECTED] wrote: Dear R-users, I have two data frames. The FIRST data frame has 100 rows, the SECOND data frame has only 50 rows. The data frames have different variables in columns (VAR1, VAR2, etc) but they share a column called ID that contains a unique identifer linking the two data frames. I would like to make a THIRD data frame containing just the rows of the FIRST data frame that match the rows (on ID) in the SECOND data frame. THIRD=data.frame(FIRST$ID,FIRST$VAR1,FIRST$VAR2) How can I modify this line to include in THIRD just the rows in FIRST that match the rows in SECOND? (i.e., contain the same value in the shared ID column). Thanks! Mark __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Chuck Cleland, Ph.D. NDRI, Inc. 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 452-1424 (M, W, F) fax: (917) 438-0894 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Convertin rows of a matrix to a list
Hello! I would like to convert rows of a matrix to a elements of a list. #For example, if I have mat-matrix(1:100,ncol=5, nrow=20) #I can do: list-apply(mat,1,list) list #however this is not quite what I want. To get what I want, I have to do: list-lapply(list,function(x)x[[1]]) list Is there a faster way? Best regards, Ales Ziberna __ R-help@stat.math.ethz.ch 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 select only certain rows when making a new dataframe?
-Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of [EMAIL PROTECTED] Sent: Tuesday, February 21, 2006 12:54 PM To: r-help@stat.math.ethz.ch Subject: [R] How to select only certain rows when making a new dataframe? Dear R-users, I have two data frames. The FIRST data frame has 100 rows, the SECOND data frame has only 50 rows. The data frames have different variables in columns (VAR1, VAR2, etc) but they share a column called ID that contains a unique identifer linking the two data frames. I would like to make a THIRD data frame containing just the rows of the FIRST data frame that match the rows (on ID) in the SECOND data frame. THIRD=data.frame(FIRST$ID,FIRST$VAR1,FIRST$VAR2) How can I modify this line to include in THIRD just the rows in FIRST that match the rows in SECOND? (i.e., contain the same value in the shared ID column). List this? one - data.frame(ID = 1:100, Var1 = runif(100), Var2 = runif(100)) two - data.frame(ID = 26:50, VarA = runif(50), VarB = runif(50)) three - one[one$ID == two$ID,] __ R-help@stat.math.ethz.ch 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] Unable to configure R 2.2.1 on Solaris 5.10 x86_64
Have we not seen such posts before from you? Perhaps you need to seek or employ a local Solaris guru (if they exist at `mac.com', or whatever your real but unstated affiliation is). The information about what is going wrong is in the file config.log, so please study it. But as it is not finding your readline headers, they appear not to be where you are saying to look. You certainly should not be having /foo/include/readline in your cpp path, just /foo/include since the headers looked for are of the form readline/history.h. On Mon, 20 Feb 2006, Aric Gregson wrote: Hello, Apologies for the post here. I have read the R-Admin (learned a lot!) and searched the web for days, but still fail at compiling R on my Ultra 20 running solaris 10 x86 1/06. This is the tail of './configure' output: . checking for dlopen in -ldl... yes checking readline/history.h usability... no checking readline/history.h presence... no checking for readline/history.h... no checking readline/readline.h usability... no checking readline/readline.h presence... no checking for readline/readline.h... no checking for rl_callback_read_char in -lreadline... no checking for main in -lncurses... no checking for main in -ltermcap... yes checking for rl_callback_read_char in -lreadline... no checking for history_truncate_file... no configure: error: --with-readline=yes (default) and headers/libs are not available My R-2.2.1/config.site has the following changes from default: CC=cc -xtarget=opteron -xarch=amd64 CFLAGS=-xO5 -xlibmil -dalign CPPFLAGS=-I/opt/sfw/include -I/opt/sfw/include/readline \ -I/opt/csw/include -I/opt/SUNWspro/include -I/opt/csw/include/readline F77=f95 -xarch=amd64 -xtarget=opteron FFLAGS=-xO5 -xlibmil -dalign LDFLAGS=-L/opt/SUNWspro/lib/amd64/ -L/opt/sfw/lib -L/opt/csw/lib CXX=CC -xarch=amd64 -xtarget=opteron CXXFLAGS=-xO5 -xlibmil -dalign R_BROWSER=mozilla MAKE=gmake My .zshrc file has the following (as I saw that some of the LDFLAGS should match my LD_LIBRARY_PATH and CONFIG_SHELL=/bin/ksh): LD_LIBRARY_PATH=/opt/sfw/include:/opt/sfw/include/readline:/usr/local/lib:/usr/X/lib:/usr/lib:/usr/ucblib:/lib:/usr/ccs/lib:/etc/lib:/usr/dt/lib:/opt/SUNWspro/lib/amd64/:/opt/sfw/lib:/opt/csw/lib You do not want include paths in the library path. Readline (version 4.2 from the Sun Companion CD) is located in /opt/sfw/include/readline. Both the readline.h and history.h files are there. Version 5.0 is located in /opt/csw/include/readline (from Blastwave). As you can see, I am new at compiling in general and Solaris specifically. Any help would be greatly appreciated. System information: System = SunOS Release = 5.10 KernelID = Generic_118844-26 Machine = i86pc OEM# = 0 Origin# = 1 NumCPU = 1 thanks, aric -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] How to get around heteroscedasticity with non-linear least squares in R?
I am using nls to fit dose-response curves but am not sure how to approach more robust regression in R to get around the problem of the my error showing increased variance with increasing dose. My understanding is that rlm or lqs would not be a good idea here. 'Fairly new to regression work, so apologies if I'm missing something obvious. [[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] Convertin rows of a matrix to a list
Perhaps something like: as.list(as.data.frame(t(mat))) Andy From: Aleš Žiberna Hello! I would like to convert rows of a matrix to a elements of a list. #For example, if I have mat-matrix(1:100,ncol=5, nrow=20) #I can do: list-apply(mat,1,list) list #however this is not quite what I want. To get what I want, I have to do: list-lapply(list,function(x)x[[1]]) list Is there a faster way? Best regards, Ales Ziberna __ R-help@stat.math.ethz.ch 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] I need your help
Hi,
Dear R users
I have problem with the following code. The matrix result must be a matrix
(3x3). But I have obtained a matrix(3x1) and I don't know why.
So, I need your help
Best regards
#
taille - function (delta, t, prob = 0.2)
{
niv.conf - c(0.90, 0.95, 0.99)
if(niv.conf - 0.90) {
t - 1.645
}
else {
if(niv.conf - 0.95) {
t - 1.96
}
t - 2.575
}
n - length(delta)
m - length(t)
result - matrix(nrow = n, ncol = m);
for (i in 1:n)
{
for(j in 1:m)
{
result[i,j]- prob*(1-prob)*((t[i])^2)/(delta[j])^2 ;
}
}
rownames(result) - delta
colnames(result) - niv.conf
round(result,2)
}
taille (delta - c( 0.01, 0.02, 0.03), niv.conf - c(0.90, 0.95, 0.99))
Lassana KOITA
Etudes de Sécurité et d'Exploitation aéroportuaires / Aerodrome Safety
Statistical analysis
Service Technique de l'Aviation Civile (STAC) / Civil Aviation Technical
Department
Direction Générale de l'Aviation Civile (DGAC) / French Civil Aviation
Authority
Tel: 01 49 56 80 60
Fax: 01 49 56 82 14
E-mail: [EMAIL PROTECTED]
http://www.stac.aviation-civile.gouv.fr/
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] indexing within panels in xyplot
On 2/21/06, Sebastian Luque [EMAIL PROTECTED] wrote:
Deepayan Sarkar [EMAIL PROTECTED] wrote:
[...]
Well, there are exceptions to this rule, but generally x and y, when
they are passed on to the panel function, are _already_ subsetted, so
x[subscripts] makes absolutely no sense. Note how your panel function
calls
panel.xyplot(x, y, ...)
without referring to subscripts at all. The subscripts argument is
there for other variables (e.g. if you were drawing confidence
intervals, and had a separate vector in your data specifying the
interval lengths). In your case, there are no other variables
involved, so just get rid of the subscripts.
Thanks Deepayan, I was indeed quite confused about this.
I realized I needed to limit the fitted line to the range of x values the
line is fit to, so I changed to panel.curve:
--cut here---start-
xyplot(y ~ x | facA, groups = facB, data = toydf,
panel.groups = function(x, y, ...) {
panel.xyplot(x, y, ...)
lindx - which(y == max(y, na.rm = TRUE))
xleft - mean(x[lindx], na.rm = TRUE)
fit - lm(y[x = xleft] ~ x[x = xleft])
panel.curve(coef(fit)[1] + (coef(fit)[2] * x),
xleft, max(x, na.rm = TRUE))
})
--cut here---end---
but can't find a way to color the line for each group differently. I
tried passing a length-2 vector as a 'col' argument to panel.curve.
Unfortunately it's only picking the first, so that both lines get colored
the same. I'm not sure, but it seems as if I need to use
'panel.superpose' directly to do this, as the help page suggests the above
would work?
The (somewhat mysterious) solution is the following:
xyplot(y ~ x | facA, groups = facB, data = toydf,
panel = panel.superpose,
panel.groups = function(x, y, col.line, ...) {
panel.xyplot(x, y, ...)
lindx - which(y == max(y, na.rm = TRUE))
xleft - mean(x[lindx], na.rm = TRUE)
fit - lm(y[x = xleft] ~ x[x = xleft])
panel.curve(coef(fit)[1] + (coef(fit)[2] * x),
col = col.line,
xleft, max(x, na.rm = TRUE))
})
This uses the fact that panel.groups is always supplied a 'col.line'
argument (along with many others) which has been suitably calculated
for each group (see panel.superpose for how this works).
You are in fact using 'panel.superpose' directly, as that's what the
panel function defaults to when there is a 'groups' argument. However,
this will change in R-2.3.0, and to use a panel.groups argument, you
will need to explicitly specify panel=panel.superpose. Sorry for the
confusion, but I believe this to be more sensible in the bigger scheme
of things.
If you now want a different set of line colors, you can
(1) either modify the superpose.line parameter, or
(2) specify col.line = c('red', 'blue') etc in the xyplot call.
Hope that makes things a bit clearer.
Deepayan
--
http://www.stat.wisc.edu/~deepayan/
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to get around heteroscedasticity with non-linear leas t squares in R?
Your understanding isn't similar to mine. Mine says robust/resistant methods are for data with heavy tails, not heteroscedasticity. The common ways to approach heteroscedasticity are transformation and weighting. The first is easy and usually quite effective for dose-response data. The second is not much harder. Both can be done in R with nls(). Andy From: Quin Wills I am using nls to fit dose-response curves but am not sure how to approach more robust regression in R to get around the problem of the my error showing increased variance with increasing dose. My understanding is that rlm or lqs would not be a good idea here. 'Fairly new to regression work, so apologies if I'm missing something obvious. [[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] Unable to configure R 2.2.1 on Solaris 5.10 x86_64
Which compiler are you using? I successfully built R for Solaris on SPARC using GCC, but it took a while to get it working the first time. The first step would be to see if you can build simple Hello World apps for C, C++ and Fortran in the directory in which you are going to install R. If you are unable to do this I would suggest that you build GCC from sources, and rather than use the Sun Software companion disk download the most recent versions of the various libraries and utilities from www.sunfreeware.com. HTH Phineas Campbell Birkbeck College University of London -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Aric Gregson Sent: Monday, February 20, 2006 11:22 PM To: r-help@stat.math.ethz.ch Subject: [R] Unable to configure R 2.2.1 on Solaris 5.10 x86_64 Hello, Apologies for the post here. I have read the R-Admin (learned a lot!) and searched the web for days, but still fail at compiling R on my Ultra 20 running solaris 10 x86 1/06. This is the tail of './configure' output: . checking for dlopen in -ldl... yes checking readline/history.h usability... no checking readline/history.h presence... no checking for readline/history.h... no checking readline/readline.h usability... no checking readline/readline.h presence... no checking for readline/readline.h... no checking for rl_callback_read_char in -lreadline... no checking for main in -lncurses... no checking for main in -ltermcap... yes checking for rl_callback_read_char in -lreadline... no checking for history_truncate_file... no configure: error: --with-readline=yes (default) and headers/libs are not available My R-2.2.1/config.site has the following changes from default: CC=cc -xtarget=opteron -xarch=amd64 CFLAGS=-xO5 -xlibmil -dalign CPPFLAGS=-I/opt/sfw/include -I/opt/sfw/include/readline \ -I/opt/csw/include -I/opt/SUNWspro/include -I/opt/csw/include/readline F77=f95 -xarch=amd64 -xtarget=opteron FFLAGS=-xO5 -xlibmil -dalign LDFLAGS=-L/opt/SUNWspro/lib/amd64/ -L/opt/sfw/lib -L/opt/csw/lib CXX=CC -xarch=amd64 -xtarget=opteron CXXFLAGS=-xO5 -xlibmil -dalign R_BROWSER=mozilla MAKE=gmake My .zshrc file has the following (as I saw that some of the LDFLAGS should match my LD_LIBRARY_PATH and CONFIG_SHELL=/bin/ksh): LD_LIBRARY_PATH=/opt/sfw/include:/opt/sfw/include/readline:/usr/local/lib:/u sr/X/lib:/usr/lib:/usr/ucblib:/lib:/usr/ccs/lib:/etc/lib:/usr/dt/lib:/opt/SU NWspro/lib/amd64/:/opt/sfw/lib:/opt/csw/lib Readline (version 4.2 from the Sun Companion CD) is located in /opt/sfw/include/readline. Both the readline.h and history.h files are there. Version 5.0 is located in /opt/csw/include/readline (from Blastwave). As you can see, I am new at compiling in general and Solaris specifically. Any help would be greatly appreciated. System information: System = SunOS Release = 5.10 KernelID = Generic_118844-26 Machine = i86pc OEM# = 0 Origin# = 1 NumCPU = 1 thanks, aric __ R-help@stat.math.ethz.ch 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] feature not available
Hi I am working with this data: my data summary is: summary(spi) open high low close volume Min. :4315 Min. :4365 Min. :4301 Min. :4352 Min. : 0 1st Qu.:4480 1st Qu.:4497 1st Qu.:4458 1st Qu.:4475 1st Qu.:11135 Median :4609 Median :4631 Median :4594 Median :4614 Median :14439 Mean :4620 Mean :4640 Mean :4599 Mean :4620 Mean :16590 3rd Qu.:4773 3rd Qu.:4796 3rd Qu.:4753 3rd Qu.:4766 3rd Qu.:18294 Max. :4944 Max. :4954 Max. :4912 Max. :4937 Max. :73559 openInterest direction volatilityvolumeXdirection Min. : 0 Min. :-62.0 Min. : 0.00 Min. :-2795242 1st Qu.:184685 1st Qu.:-19.25000 1st Qu.:30.00 1st Qu.: -248740 Median :193233 Median : -1.5 Median :38.50 Median : -15905 Mean :188825 Mean : 0.01563 Mean :41.58 Mean :6275 3rd Qu.:199800 3rd Qu.: 17.0 3rd Qu.:50.00 3rd Qu.: 206325 Max. :236759 Max. : 74.0 Max. :94.00 Max. : 2024470 volatilityXdirection upDown nextDay Min. : 0 Min. :0. Min. :0. 1st Qu.: 362816 1st Qu.:0. 1st Qu.:0. Median : 540187 Median :0. Median :0. Mean : 731595 Mean :0.4844 Mean :0.4844 3rd Qu.: 996650 3rd Qu.:1. 3rd Qu.:1. Max. :3604391 Max. :1. Max. :1. and trying to do a glm like this: logistic.model = glm(formula = as.factor(nextDay) ~ ., family=binomial, data=spi[1:50,]) summary(logistic.model) Call: glm(formula = as.factor(nextDay) ~ ., family = binomial, data = spi[1:50, ]) Deviance Residuals: Min 1Q Median 3Q Max -1.9994 -0.8575 -0.0330 0.7961 1.8376 Coefficients: (2 not defined because of singularities) Estimate Std. Error z value Pr(|z|) (Intercept) -3.016e+01 2.934e+01 -1.028 0.3040 open -2.514e-02 8.939e-02 -0.281 0.7785 high 7.598e-02 1.295e-01 0.587 0.5575 low -1.065e-01 1.176e-01 -0.905 0.3656 close 5.943e-02 7.995e-02 0.743 0.4573 volume1.960e-04 1.977e-04 0.991 0.3215 openInterest 5.728e-05 5.266e-05 1.088 0.2767 directionNA NA NA NA volatility NA NA NA NA volumeXdirection 3.605e-06 3.765e-06 0.958 0.3383 volatilityXdirection -5.815e-06 5.708e-06 -1.019 0.3082 upDown -2.561e+00 1.259e+00 -2.034 0.0419 * --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 69.235 on 49 degrees of freedom Residual deviance: 55.000 on 40 degrees of freedom AIC: 75 Number of Fisher Scoring iterations: 7 I am getting NA for direction and volatility. I got it both before and after casting them as.numeric. Can anyone tell me what I am doing wrong? Thanks Stephen -- 20/02/2006 [[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] feature not available
From: Stephen Choularton Hi I am working with this data: my data summary is: summary(spi) open high low close volume Min. :4315 Min. :4365 Min. :4301 Min. :4352 Min. : 0 1st Qu.:4480 1st Qu.:4497 1st Qu.:4458 1st Qu.:4475 1st Qu.:11135 Median :4609 Median :4631 Median :4594 Median :4614 Median :14439 Mean :4620 Mean :4640 Mean :4599 Mean :4620 Mean :16590 3rd Qu.:4773 3rd Qu.:4796 3rd Qu.:4753 3rd Qu.:4766 3rd Qu.:18294 Max. :4944 Max. :4954 Max. :4912 Max. :4937 Max. :73559 openInterest direction volatility volumeXdirection Min. : 0 Min. :-62.0 Min. : 0.00 Min. :-2795242 1st Qu.:184685 1st Qu.:-19.25000 1st Qu.:30.00 1st Qu.: -248740 Median :193233 Median : -1.5 Median :38.50 Median : -15905 Mean :188825 Mean : 0.01563 Mean :41.58 Mean :6275 3rd Qu.:199800 3rd Qu.: 17.0 3rd Qu.:50.00 3rd Qu.: 206325 Max. :236759 Max. : 74.0 Max. :94.00 Max. : 2024470 volatilityXdirection upDown nextDay Min. : 0 Min. :0. Min. :0. 1st Qu.: 362816 1st Qu.:0. 1st Qu.:0. Median : 540187 Median :0. Median :0. Mean : 731595 Mean :0.4844 Mean :0.4844 3rd Qu.: 996650 3rd Qu.:1. 3rd Qu.:1. Max. :3604391 Max. :1. Max. :1. and trying to do a glm like this: logistic.model = glm(formula = as.factor(nextDay) ~ ., family=binomial, data=spi[1:50,]) summary(logistic.model) Call: glm(formula = as.factor(nextDay) ~ ., family = binomial, data = spi[1:50, ]) Deviance Residuals: Min 1Q Median 3Q Max -1.9994 -0.8575 -0.0330 0.7961 1.8376 Coefficients: (2 not defined because of singularities) ^^ The answer is most likely here: You have multicollinearity. Andy Estimate Std. Error z value Pr(|z|) (Intercept) -3.016e+01 2.934e+01 -1.028 0.3040 open -2.514e-02 8.939e-02 -0.281 0.7785 high 7.598e-02 1.295e-01 0.587 0.5575 low -1.065e-01 1.176e-01 -0.905 0.3656 close 5.943e-02 7.995e-02 0.743 0.4573 volume1.960e-04 1.977e-04 0.991 0.3215 openInterest 5.728e-05 5.266e-05 1.088 0.2767 directionNA NA NA NA volatility NA NA NA NA volumeXdirection 3.605e-06 3.765e-06 0.958 0.3383 volatilityXdirection -5.815e-06 5.708e-06 -1.019 0.3082 upDown -2.561e+00 1.259e+00 -2.034 0.0419 * --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 69.235 on 49 degrees of freedom Residual deviance: 55.000 on 40 degrees of freedom AIC: 75 Number of Fisher Scoring iterations: 7 I am getting NA for direction and volatility. I got it both before and after casting them as.numeric. Can anyone tell me what I am doing wrong? Thanks Stephen -- 20/02/2006 [[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] Compute a correlation matrix from an existing covariance matrix
Peter Dalgaard [EMAIL PROTECTED] writes: (BTW, why doesn't the page for cor/var/cov/cov2cor show up on help.search(covariance) or help.search(correlation) for me??) It's there now. Must have been a temporary condition, possibly in the user... -- 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] How to get around heteroscedasticity with non-linear leas t squares in R?
Liaw, Andy [EMAIL PROTECTED] writes: Your understanding isn't similar to mine. Mine says robust/resistant methods are for data with heavy tails, not heteroscedasticity. The common ways to approach heteroscedasticity are transformation and weighting. The first is easy and usually quite effective for dose-response data. The second is not much harder. Both can be done in R with nls(). And there is gnls() which allows direct modelling of the variance. -p Andy From: Quin Wills I am using nls to fit dose-response curves but am not sure how to approach more robust regression in R to get around the problem of the my error showing increased variance with increasing dose. My understanding is that rlm or lqs would not be a good idea here. 'Fairly new to regression work, so apologies if I'm missing something obvious. [[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 -- 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] [OT] New book on software patents, etc.
I just saw a review of a book about software patents called Math you can't use http://www.brookings.edu/press/books/mathyoucantuse.htm which looks very interesting. The whole subject of software patents and licenses on code is important in the R Project. __ R-help@stat.math.ethz.ch 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 around heteroscedasticity with non-linear leas t squares in R?
On Tue, 21 Feb 2006, Peter Dalgaard wrote: Liaw, Andy [EMAIL PROTECTED] writes: Your understanding isn't similar to mine. Mine says robust/resistant methods are for data with heavy tails, not heteroscedasticity. The common ways to approach heteroscedasticity are transformation and weighting. The first is easy and usually quite effective for dose-response data. The second is not much harder. Both can be done in R with nls(). And there is gnls() which allows direct modelling of the variance. in package nlme, BTW. R-devel allows weights in nls, which makes it easier for those most familiar with that function. -p Andy From: Quin Wills I am using nls to fit dose-response curves but am not sure how to approach more robust regression in R to get around the problem of the my error showing increased variance with increasing dose. My understanding is that rlm or lqs would not be a good idea here. 'Fairly new to regression work, so apologies if I'm missing something obvious. [[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 -- 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] I need your help
Hello,
I made some adjustments to your code.
1) You used 't' twice in your code. When you assign t
- 1.645, the length m of 't' becomes 1 (that's why
your matrix was 3 x 1)
2) I assumed that the 'critical value' t is the
quantile of a Normal (0,1). By using 't -
qnorm(0.5*level + 0.5)' allows you to consider any
value between 0 and 1 for level. Otherwise, you must
write a 'neverending' list of 'if else'. If I assumed
wrong, put back in the code the 'if else' statement.
3) The 'kronecker' function avoids the loops for 'i'
and 'j'
I hope this is what you want.
Marco Geraci
##
taille - function (delta, level, prob = 0.2) {
t - qnorm(0.5*level + 0.5)
n - length(delta)
m - length(level)
result -
prob*(1-prob)*matrix(kronecker(t^2,delta^2,FUN=/),
n, m, byrow=T)
rownames(result) - delta
colnames(result) - level
round(result,2)
}
taille (delta = c(0.01, 0.02, 0.03), level = c(0.90,
0.95, 0.99))
#
--- KOITA Lassana - STAC/ACE
[EMAIL PROTECTED] wrote:
Hi,
Dear R users
I have problem with the following code. The matrix
result must be a matrix
(3x3). But I have obtained a matrix(3x1) and I don't
know why.
So, I need your help
Best regards
#
taille - function (delta, t, prob = 0.2)
{
niv.conf - c(0.90, 0.95, 0.99)
if(niv.conf - 0.90) {
t - 1.645
}
else {
if(niv.conf - 0.95) {
t - 1.96
}
t - 2.575
}
n - length(delta)
m - length(t)
result - matrix(nrow = n, ncol = m);
for (i in 1:n)
{
for(j in 1:m)
{
result[i,j]-
prob*(1-prob)*((t[i])^2)/(delta[j])^2 ;
}
}
rownames(result) - delta
colnames(result) - niv.conf
round(result,2)
}
taille (delta - c( 0.01, 0.02, 0.03), niv.conf -
c(0.90, 0.95, 0.99))
Lassana KOITA
Etudes de Sécurité et d'Exploitation aéroportuaires
/ Aerodrome Safety
Statistical analysis
Service Technique de l'Aviation Civile (STAC) /
Civil Aviation Technical
Department
Direction Générale de l'Aviation Civile (DGAC) /
French Civil Aviation
Authority
Tel: 01 49 56 80 60
Fax: 01 49 56 82 14
E-mail: [EMAIL PROTECTED]
http://www.stac.aviation-civile.gouv.fr/
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch 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 select only certain rows when making a new dataframe?
I think ?merge is what you need. The help page of this function is veryhelpful Marco --- [EMAIL PROTECTED] wrote: Dear R-users, I have two data frames. The FIRST data frame has 100 rows, the SECOND data frame has only 50 rows. The data frames have different variables in columns (VAR1, VAR2, etc) but they share a column called ID that contains a unique identifer linking the two data frames. I would like to make a THIRD data frame containing just the rows of the FIRST data frame that match the rows (on ID) in the SECOND data frame. THIRD=data.frame(FIRST$ID,FIRST$VAR1,FIRST$VAR2) How can I modify this line to include in THIRD just the rows in FIRST that match the rows in SECOND? (i.e., contain the same value in the shared ID column). Thanks! Mark __ R-help@stat.math.ethz.ch 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 get around heteroscedasticity with non-linear least squares in R?
Quin Wills wrote: I am using nls to fit dose-response curves but am not sure how to approach more robust regression in R to get around the problem of the my error showing increased variance with increasing dose. package sfsmisc has rnls (robust nls) which might be of use. Kjetil My understanding is that rlm or lqs would not be a good idea here. 'Fairly new to regression work, so apologies if I'm missing something obvious. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Composing a matrix for 'heatmap' from original dataset and genetic data analysis.
Hello.
I have following data:
dozeg1 g2 g3 g4 etc...
1200a d
1100b e
1158c f
1500a e
1403c f
...
etc.
Doze - numeric, levels(g1)=c('a','b','c'), levels(g2)=c('d','e','f').
This is the data of genetic analysis, where the factors are genes and
the numerical is some numerical value.
I want to build up a function, which will do a preliminary data
visualization based on factor levels combination and mean doze for that
combination.
The best way to do it, as I think, is to build a tree which divides at
the stage 1 into the levels of factor 1, then on the stage 2 into
the levels of factor two etc. And in the end it will be mean doze for
this branch. What I found was only 'ctree', this can only bifurcate, I'm
not sure it fits much. What's the convenient way to visualize such kind
of data in R?
I think 'heatmap' for any two factors I give to it with text mean labels
(how to add them?) would be nice solution. So that's the way I wrote
(ughly):
zar - function(g1,g2,param='doz') {
levels(g1)-lg1
levels(g2)-lg2
table(lg1,lg2)-f
#to make a matrix for 'heatmap', is it necessary?
#data labels are not preserved in this expression, badly!
for (i in lg1) {
for (j in lg2) {
f[i,j]-mean(subset(mydata,as.factor(g1)==ias.factor(g2)==j)[[param]],
na.rm=T)
}
}
heatmap(-f, Rowv=NA, Colv=NA, margins=c(10,10),main=param)
}
Please, make a hint to me for how to do this task easier and more
convenient,
how to add text labels to 'heatmap'
and is this way right?
--
Evgeniy
__
R-help@stat.math.ethz.ch 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] Limit of matrix + naming
Hi all, I have read a data matrix of 304 x 404 using read.table. When I am trying to name the colnames, with first try colnames(L5)-list(2:305), I keep getting error message such as Error in colnames-(`*tmp*`, value = list(c(2, 3, 4, 5, 6, 7, 8, 9, : length of 'dimnames' [2] not equal to array extent and I don't know why. But, if I look at the original names in terms of V1, V2, the rownames are repeated again after V256. Is there a limitation? Thanks. Yen Lin [[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] Selecting amongst large classes of models (Was: Nested AIC)
Professor Ripley, On Tue, Feb 21, 2006 at 05:36:15PM +, Prof Brian Ripley wrote: A recent talk of mine http://www.stats.ox.ac.uk/~ripley/Nelder80.pdf may be illuminating. There is a published paper version. Would you mind providing a citation for that published paper version? I do not find details on your website, and Current Contents does not provide any clues. Cheers, Andrew -- Andrew Robinson Department of Mathematics and StatisticsTel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 Email: [EMAIL PROTECTED] http://www.ms.unimelb.edu.au __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to get around heteroscedasticity with non-linear least squares in R?
Instead of thinking that the heteroscedasticity is a nuisance and something to get around, i.e, just wanting weighted estimates of the mean function, you might want to think about what heteroscedasticity is telling you and estimate some other quantities. Heteroscedasticity is telling you that the conditional distributions don't change at a constant rate across all portions of the distribution (think percentiles or more generally quantiles) and, therefore, a function for the mean (no matter how precisely estimated) cannot tell you all there is to know about your dose-response relation. Why not go after estimating the conditional quantile functions directly with nonlinear quantile regression, function nlrq() in the quantreg package? Brian Brian S. Cade U. S. Geological Survey Fort Collins Science Center 2150 Centre Ave., Bldg. C Fort Collins, CO 80526-8818 email: [EMAIL PROTECTED] tel: 970 226-9326 Kjetil Brinchmann Halvorsen [EMAIL PROTECTED] Sent by: [EMAIL PROTECTED] 02/21/2006 03:31 PM Please respond to [EMAIL PROTECTED] To Quin Wills [EMAIL PROTECTED] cc r-help@stat.math.ethz.ch Subject Re: [R] How to get around heteroscedasticity with non-linear least squares in R? Quin Wills wrote: I am using nls to fit dose-response curves but am not sure how to approach more robust regression in R to get around the problem of the my error showing increased variance with increasing dose. package sfsmisc has rnls (robust nls) which might be of use. Kjetil My understanding is that rlm or lqs would not be a good idea here. 'Fairly new to regression work, so apologies if I'm missing something obvious. [[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 [[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] R[1304] general protection rip - Error running R file in Batch mode
hi group, I have a list of 100 elements (genes), that are to be queried in postgres. The out of the file is processed for fisher test and written back as a file with pvalues. For one such operation there are 54 lines of code. I replicated the same code changing the query element and file name to write back. I created a jumbo 100*54 = 5400 lines of code by replication. I saved this file as .R and ran it as: R CMD BATCH myfile.R.. The program ran for ~40 mins. ( each operation takes 5 mins). after going home, I checked the stats and saw that the program got aborted leaving 8 output files in directory. I checked /var/log/messages and found : localhost kernel: R[1304] general protection rip:2a988c919e rsp:7fbfff87f0 error:0 postgres did not log any error message. This is my first time I am running my R code in BATCH mode. I am a new to programming in R and I have been working directly with R interface. Could any one guess what could have been the problem. Thanks Sri __ R-help@stat.math.ethz.ch 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] Limit of matrix + naming
Looking at length(list(2:305)) and str(list(2:305)) can be instructive. You don't need the list(), and you probably want to use `legal' names for variables. Andy From: Chia, Yen Lin Hi all, I have read a data matrix of 304 x 404 using read.table. When I am trying to name the colnames, with first try colnames(L5)-list(2:305), I keep getting error message such as Error in colnames-(`*tmp*`, value = list(c(2, 3, 4, 5, 6, 7, 8, 9, : length of 'dimnames' [2] not equal to array extent and I don't know why. But, if I look at the original names in terms of V1, V2, the rownames are repeated again after V256. Is there a limitation? Thanks. Yen Lin [[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] R[1304] general protection rip - Error running R file in Batch mode
Hi group,
I ran the file again and after 10 min. of run, the
process got interrupted.
Here is the output of the file
$bash: R CMD BATCH grand.R grand.log
R CMD BATCH grand.R
/usr/local/lib64/R/bin/BATCH: line 55: 1758 Broken
pipe ( echo invisible(options(echo =
TRUE)); cat ${in}; echo ''; echo proc.time() )
1760 Segmentation fault | ${R_HOME}/bin/R
${opts} ${out} 21
Is this because of the excess query code in one .R
file?
Any suggestions, will help me .
thanks
Sri
I created a jumbo 100*54 = 5400 lines of code by
replication.
I saved this file as .R and ran it as:
R CMD BATCH myfile.R..
The program ran for ~40 mins. ( each operation takes
5
mins). after going home, I checked the stats and saw
that the program got aborted leaving 8 output files
in
directory.
I checked /var/log/messages and found :
localhost kernel: R[1304] general protection
rip:2a988c919e rsp:7fbfff87f0 error:0
postgres did not log any error message.
This is my first time I am running my R code in
BATCH
mode. I am a new to programming in R and I have been
working directly with R interface.
Could any one guess what could have been the
problem.
Thanks
Sri
__
R-help@stat.math.ethz.ch 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] var-covar matrices comparison
Date: Mon, 20 Feb 2006 16:43:55 -0600 From: Aldi Kraja [EMAIL PROTECTED] Hi, Using package gclus in R, I have created some graphs that show the trends within subgroups of data and correlations among 9 variables (v1-v9). Being interested for more details on these data I have produced also the var-covar matrices. Question: From a pair of two subsets of data (with 9 variables each, I have two var-covar matrices for each subgroup, that differ for a treatment on one group (treatment A) vs (non-Treatment A). Is there a software that can compare if two var-covar matrices are statistically the same? This can be done in various structural equation modelling packages. I don't think it can be done automatically using the sem package, as that does not allow multiple groups. You can roll your own LR test (assuming MVN): f - (N-1) * (log(det(E)) - log(det(O)) + sum(diag((O %*% solve(E-p) N=size of group p=number of variables E=expected covariance matrix O=observed covariance matrix where in your example, E will be the observed covariance matrix for the pooled groups. There are GLS etc alternatives - see eg Bollen's book on SEM. | David Duffy (MBBS PhD) ,-_|\ | email: [EMAIL PROTECTED] ph: INT+61+7+3362-0217 fax: -0101 / * | Epidemiology Unit, Queensland Institute of Medical Research \_,-._/ | 300 Herston Rd, Brisbane, Queensland 4029, Australia GPG 4D0B994A v __ R-help@stat.math.ethz.ch 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] Ranking within factor subgroups
Hi,
I have a dataframe, x of the following form:
DateSymbol AB C
20041201 ABC 10 12 15
20041201 DEF 95 4
...
20050101 ABC 5 3 1
20050101 GHM 12 42
here A, B,C are properties of a set symbols recorded for a given date.
I wante to decile the symbols For each date and property and
create another set of columns bucketA,bucketB, bucketC containing the
decile rank
for each symbol. The following non-vectorized code does what I want,
bucket - function(data,nBuckets) {
q - quantile(data,seq(0,1,len=nBuckets+1),na.rm=T)
q[1] - q[1] - 0.1 # need to do this to ensure there are no extra NAs
cut(data,q,include.lowest=T,labels=F)
}
calcDeciles - function(x,colNames) {
nBuckets - 10
dates - unique(x$Date)
for ( date in dates) {
iVec - x$Date == date
xx - x[iVec,]
for (colName in colNames) {
data - xx[,colName]
bColName - paste(bucket,colName,sep=)
x[iVec,bColName] - bucket(data,nBuckets)
}
}
x
}
x - calcDeciles(x,c(A,B,C))
I was wondering if it is possible to vectorize the above function to make it
more efficient.
I tried,
rlist - tapply(x$A,x$Date,bucket)
but I am not sure how to assign the contents of rlist to their appropriate
slots in the original
dataframe.
Thanks,
Maneesh
__
R-help@stat.math.ethz.ch 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[1304] general protection rip - Error running R file in Batch mode
On 2/21/2006 8:24 PM, Srinivas Iyyer wrote: hi group, I have a list of 100 elements (genes), that are to be queried in postgres. The out of the file is processed for fisher test and written back as a file with pvalues. For one such operation there are 54 lines of code. I replicated the same code changing the query element and file name to write back. I created a jumbo 100*54 = 5400 lines of code by replication. That's a fairly strange way to write such a thing. Why not write a single function, and run it in a loop? I saved this file as .R and ran it as: R CMD BATCH myfile.R.. The program ran for ~40 mins. ( each operation takes 5 mins). after going home, I checked the stats and saw that the program got aborted leaving 8 output files in directory. I checked /var/log/messages and found : localhost kernel: R[1304] general protection rip:2a988c919e rsp:7fbfff87f0 error:0 postgres did not log any error message. This is my first time I am running my R code in BATCH mode. I am a new to programming in R and I have been working directly with R interface. Could any one guess what could have been the problem. You haven't really given us much to go on. You haven't stated the version, or what packages and functions you are using. I'd guess this is a bug with one of the packages you're using, but it could be a bug in R. Can you determine which line causes the error? Does it happen if you extract just that line, and enough previous ones to set things up for it? Duncan Murdoch __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to get around heteroscedasticity with non-linear leas t squares in R?
Thank you all, this has been a great help (including the methodological
advice). Very interesting - I'll be sure to read the lecture.
Quin
-Original Message-
From: Liaw, Andy [mailto:[EMAIL PROTECTED]
Sent: 22 February 2006 01:18
To: 'Brian S Cade'; [EMAIL PROTECTED]
Cc: Quin Wills; r-help@stat.math.ethz.ch; [EMAIL PROTECTED]
Subject: RE: [R] How to get around heteroscedasticity with non-linear leas t
squares in R?
From: Brian S Cade
Instead of thinking that the heteroscedasticity is a nuisance and
something to get around, i.e, just wanting weighted
estimates of the
mean function, you might want to think about what
heteroscedasticity is
telling you and estimate some other quantities.
Indeed! See Prof. Carroll's 2002 Fisher Lecture:
http://www.stat.tamu.edu/ftp/pub/rjcarroll/2003.papers.directory/published_F
isher_Lecture.pdf
(There's Powerpoint file on his web page, too.)
Andy
Heteroscedasticity is
telling you that the conditional distributions don't change
at a constant
rate across all portions of the distribution (think
percentiles or more
generally quantiles) and, therefore, a function for the mean
(no matter
how precisely estimated) cannot tell you all there is to know
about your
dose-response relation. Why not go after estimating the conditional
quantile functions directly with nonlinear quantile
regression, function
nlrq() in the quantreg package?
Brian
Brian S. Cade
U. S. Geological Survey
Fort Collins Science Center
2150 Centre Ave., Bldg. C
Fort Collins, CO 80526-8818
email: [EMAIL PROTECTED]
tel: 970 226-9326
Kjetil Brinchmann Halvorsen [EMAIL PROTECTED]
Sent by: [EMAIL PROTECTED]
02/21/2006 03:31 PM
Please respond to
[EMAIL PROTECTED]
To
Quin Wills [EMAIL PROTECTED]
cc
r-help@stat.math.ethz.ch
Subject
Re: [R] How to get around heteroscedasticity with non-linear
least squares
in R?
Quin Wills wrote:
I am using nls to fit dose-response curves but am not sure how to
approach
more robust regression in R to get around the problem of
the my error
showing increased variance with increasing dose.
package sfsmisc has rnls (robust nls)
which might be of use.
Kjetil
My understanding is that rlm or lqs would not be a good
idea here.
'Fairly new to regression work, so apologies if I'm missing
something
obvious.
[[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
[[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
--
Notice: This e-mail message, together with any attachments,...{{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] panel-VAR
Have you looked at the nlme package? If you provide more detail about your application (as suggested in the posting guide! www.R-project.org/posting-guide.html), you might receive more useful replies more quickly. hope this helps, spencer graves [EMAIL PROTECTED] wrote: Hi, Anybody knows if there is a package for panel VAR ? thanks __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch 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] OT Futility Analysis
What does this particular Steering Committee think a futility analysis is? Do they have any particular reference(s)? What do you find in your own literature review? If it were my problem, I think I'd start with questions like that. Your comments suggested to me a confounding of technical and political problems. The politics suggests the language you need to use in your response. Beyond that, I've never heard before of a futility analysis, but I think I could do one by just trying to be clear about the options the Steering Committee might consider plausible and then comparing them with appropriate simulations -- summarized as confidence intervals, as you suggest. And I hope that someone else will enlighten us both if there are better options available. Best Wishes, spencer graves p.s. For any attorneys who may read these comments, the suggestions are obviously warranteed up to the amount you paid for it, which is nothing. If you follow them and they turn out to be inappropriate, you will pay the price. I encourage you to share the problems with me, so I can learn from the experience. However, the limits of my liability are as already stated. Kevin E. Thorpe wrote: I beg your pardon if this is too off topic. I am posting here since I hope to find an R solution to my problem. Please indulge me while I give a little background about what I'm trying to do. I'm on a DSMB for a clinical trial. The Steering Committee for the trial has asked us to perform a futility analysis on their primary outcome which is a time-to-event endpoint. The trial was not designed with group sequential methods, nor was any futility analysis spelled out in the protocol. Another thing which may be relevant is that due to circumstances beyond the investigators' control, the trial will stop recruitment prematurely unless there is some compelling reason for them to find a way to continue the trial. Lastly, the trial has accrued not quite half of the planned sample size. Admittedly, I don't have a vast amount of experience implementing stopping rules. In other protocols I have seen where futility analyses have been planned but a group sequential design has not otherwise been employed, conditional power has been used for the futility rule. So naturally, that was my first thought (although I may well be wrong) in this case. I have done RSiteSearch() with the following terms (three different searches): futility analysis conditional power stochastic curtailment Nothing that looked relevant to my problem jumped out at me. I have read, somewhat recently, that there are problems with conditional power, although I don't remember the details at the moment. This has prompted me to consider other approaches to the problem. One simple thing that has occurred to me, although I don't know what the implications are is to simply look at a confidence interval around the hazard ratio for the treatment effect. In the event that the CI includes 1 and excludes any clinically important difference, I would take that as an indication of futility. I would appreciate your comments on this and to learn of any more formal methods, particularly of implementations in R. Thank you for reading. Kevin __ R-help@stat.math.ethz.ch 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] Ranking within factor subgroups
It might help to give a simple reproducible example in the future. For
example
df - cbind.data.frame( date=rep( 1:5, each=100 ), A=rpois(500, 100),
B=rpois(500, 50), C=rpois(500, 30) )
might generate something like
date A B C
11 93 51 32
21 95 51 30
31 102 59 28
41 105 52 32
51 105 53 26
61 99 59 37
.... ... .. ..
4955 100 57 19
4965 96 47 44
4975 111 56 35
4985 105 49 23
4995 105 61 30
5005 92 53 32
Here is my proposed solution. Can you double check with your existing
functions to see if they are correct.
decile.fn - function(x, nbreaks=10){
br - quantile( x, seq(0, 1, len=nbreaks+1), na.rm=T )
br[1] - -Inf
return( cut(x, br, labels=F) )
}
out - apply( df[ ,c(A, B, C)], 2,
function(v) unlist( tapply( v, df$date, decile.fn ) ) )
rownames(out) - rownames(df)
out - cbind(df$date, out)
Regards, Adai
On Tue, 2006-02-21 at 21:44 -0500, maneesh deshpande wrote:
Hi,
I have a dataframe, x of the following form:
DateSymbol AB C
20041201 ABC 10 12 15
20041201 DEF 95 4
...
20050101 ABC 5 3 1
20050101 GHM 12 42
here A, B,C are properties of a set symbols recorded for a given date.
I wante to decile the symbols For each date and property and
create another set of columns bucketA,bucketB, bucketC containing the
decile rank
for each symbol. The following non-vectorized code does what I want,
bucket - function(data,nBuckets) {
q - quantile(data,seq(0,1,len=nBuckets+1),na.rm=T)
q[1] - q[1] - 0.1 # need to do this to ensure there are no extra NAs
cut(data,q,include.lowest=T,labels=F)
}
calcDeciles - function(x,colNames) {
nBuckets - 10
dates - unique(x$Date)
for ( date in dates) {
iVec - x$Date == date
xx - x[iVec,]
for (colName in colNames) {
data - xx[,colName]
bColName - paste(bucket,colName,sep=)
x[iVec,bColName] - bucket(data,nBuckets)
}
}
x
}
x - calcDeciles(x,c(A,B,C))
I was wondering if it is possible to vectorize the above function to make it
more efficient.
I tried,
rlist - tapply(x$A,x$Date,bucket)
but I am not sure how to assign the contents of rlist to their appropriate
slots in the original
dataframe.
Thanks,
Maneesh
__
R-help@stat.math.ethz.ch 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] How do I tell it which directory to use?
From Tom: In R 2.2.0 under Windows, I want to be able to give it a filename such as myFile.txt without the quotes. But actually I mean: C:\Documents and Settings\Tom\My Documents\qpaper7\R Project Started 19 Dec 05\myFile.txt If I were to repeat this each time, my computer would get all bored and cranky and start to drop bits (only a joke, of course). I think I want to set the Home directory or the working directory or some directory or other to the above directory. I may or may not want to set some environmental variables. R 2.2.0; working directly from the console and copying and pasting code which I want to test into the console. Windows XP Home Edition. Administrator privileges are enabled. A curve ball: There are two accounts, Tom and Jones; the data are stored under Tom, whereas the computation is being done under the Jones account. I won't bore you with the details of why I am doing this. I was able to call Sys.getenv (R_USER) and get the home directory. I am a newbie to R and not familiar with the terminology. Tom Thomas L. Jones, Ph.D., Computer Science __ R-help@stat.math.ethz.ch 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] stripchart-y axis labels
Hello, I am trying to create a stripchart for my data, where y axis labels are characters (ie,names of cities). I would like to change the orientation of the y - axis labels, ie perpendicular to y axis. Below is the code i am using: par(srt=90) with(ozone.ne.trim, stripchart(Median~City,main = stripchart(ozone))) The par option doesn't seem to working. Kindly help. Thanks mathangi. __ R-help@stat.math.ethz.ch 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 do I tell it which directory to use?
I think you want getwd() and setwd(). HTH, Simon. Thomas L Jones wrote: From Tom: In R 2.2.0 under Windows, I want to be able to give it a filename such as myFile.txt without the quotes. But actually I mean: C:\Documents and Settings\Tom\My Documents\qpaper7\R Project Started 19 Dec 05\myFile.txt If I were to repeat this each time, my computer would get all bored and cranky and start to drop bits (only a joke, of course). I think I want to set the Home directory or the working directory or some directory or other to the above directory. I may or may not want to set some environmental variables. R 2.2.0; working directly from the console and copying and pasting code which I want to test into the console. Windows XP Home Edition. Administrator privileges are enabled. A curve ball: There are two accounts, Tom and Jones; the data are stored under Tom, whereas the computation is being done under the Jones account. I won't bore you with the details of why I am doing this. I was able to call Sys.getenv (R_USER) and get the home directory. I am a newbie to R and not familiar with the terminology. Tom Thomas L. Jones, Ph.D., Computer Science __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Simon Blomberg, B.Sc.(Hons.), Ph.D, M.App.Stat. Centre for Resource and Environmental Studies The Australian National University Canberra ACT 0200 Australia T: +61 2 6125 7800 email: Simon.Blomberg_at_anu.edu.au F: +61 2 6125 0757 CRICOS Provider # 00120C __ R-help@stat.math.ethz.ch 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 do I tell it which directory to use?
Tom, You can define your working directory by using: setwd(C:\Documents and Settings\Tom\My Documents\qpaper7\R Project Started 19 Dec 05) check that your file is there: list.files() and then use: source(myFile.txt) the machine should load myFile You can go to another directory: setwd(anotherdir) and repeat the procedure. Or even better if you define a number of directories in an external file: dir1 - c(C:\Documents and Settings\Tom\My Documents\qpaper7\) dir2 - c(C:\Documents and Settings\Tom\My Documents\) and after loading the file at the beginning of the sesion you can use: setwd(dir1) etc. Is it of any help to you? Cheers, Augusto Augusto Sanabria. MSc, PhD. Mathematical Modeller Risk Research Group Geospatial Earth Monitoring Division Geoscience Australia (www.ga.gov.au) Cnr. Jerrabomberra Av. Hindmarsh Dr. Symonston ACT 2609 Ph. (02) 6249-9155 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Thomas L Jones Sent: Wednesday, 22 February 2006 4:31 PM To: R-project help Subject: [R] How do I tell it which directory to use? From Tom: In R 2.2.0 under Windows, I want to be able to give it a filename such as myFile.txt without the quotes. But actually I mean: C:\Documents and Settings\Tom\My Documents\qpaper7\R Project Started 19 Dec 05\myFile.txt If I were to repeat this each time, my computer would get all bored and cranky and start to drop bits (only a joke, of course). I think I want to set the Home directory or the working directory or some directory or other to the above directory. I may or may not want to set some environmental variables. R 2.2.0; working directly from the console and copying and pasting code which I want to test into the console. Windows XP Home Edition. Administrator privileges are enabled. A curve ball: There are two accounts, Tom and Jones; the data are stored under Tom, whereas the computation is being done under the Jones account. I won't bore you with the details of why I am doing this. I was able to call Sys.getenv (R_USER) and get the home directory. I am a newbie to R and not familiar with the terminology. Tom Thomas L. Jones, Ph.D., Computer Science __ R-help@stat.math.ethz.ch 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] stripchart-y axis labels
Try par(las=2) Cheers, Simon. [EMAIL PROTECTED] wrote: Hello, I am trying to create a stripchart for my data, where y axis labels are characters (ie,names of cities). I would like to change the orientation of the y - axis labels, ie perpendicular to y axis. Below is the code i am using: par(srt=90) with(ozone.ne.trim, stripchart(Median~City,main = stripchart(ozone))) The par option doesn't seem to working. Kindly help. Thanks mathangi. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Simon Blomberg, B.Sc.(Hons.), Ph.D, M.App.Stat. Centre for Resource and Environmental Studies The Australian National University Canberra ACT 0200 Australia T: +61 2 6125 7800 email: Simon.Blomberg_at_anu.edu.au F: +61 2 6125 0757 CRICOS Provider # 00120C __ R-help@stat.math.ethz.ch 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] Gram-Charlier series
Good day everyone, I want to use the Gram-Charlier series expansion to model some data. To do that, I need functions to: 1) Calculate 'n' moments from given data 2) Transform 'n' moments to 'n' central moments, or 3) Transform 'n' moments to 'n' cumulants 4) Calculate a number of Hermite polynomials Are there R-functions to do any of the above? (mean, sd and cum3 are very limited) Thank you for your help, Augusto Augusto Sanabria. MSc, PhD. Mathematical Modeller Risk Research Group Geospatial Earth Monitoring Division Geoscience Australia (www.ga.gov.au) Cnr. Jerrabomberra Av. Hindmarsh Dr. Symonston ACT 2609 Ph. (02) 6249-9155 __ R-help@stat.math.ethz.ch 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 do I tell it which directory to use?
[EMAIL PROTECTED] wrote: Tom, You can define your working directory by using: setwd(C:\Documents and Settings\Tom\My Documents\qpaper7\R Project Started 19 Dec 05) Thanks for your contribution to R-help, but: No, see the R for Windows FAQs which tell you to use \\ or / rather than \! check that your file is there: list.files() and then use: source(myFile.txt) the machine should load myFile You can go to another directory: setwd(anotherdir) and repeat the procedure. Or even better if you define a number of directories in an external file: dir1 - c(C:\Documents and Settings\Tom\My Documents\qpaper7\) dir2 - c(C:\Documents and Settings\Tom\My Documents\) See above + don't forget the quotes! and after loading the file at the beginning of the sesion you can use: setwd(dir1) etc. In this case without quotes! setwd(dir1) Please don't confuse other R-help readers and try to be more precise in your answers. Uwe Ligges Is it of any help to you? Cheers, Augusto Augusto Sanabria. MSc, PhD. Mathematical Modeller Risk Research Group Geospatial Earth Monitoring Division Geoscience Australia (www.ga.gov.au) Cnr. Jerrabomberra Av. Hindmarsh Dr. Symonston ACT 2609 Ph. (02) 6249-9155 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Thomas L Jones Sent: Wednesday, 22 February 2006 4:31 PM To: R-project help Subject: [R] How do I tell it which directory to use? From Tom: In R 2.2.0 under Windows, I want to be able to give it a filename such as myFile.txt without the quotes. But actually I mean: C:\Documents and Settings\Tom\My Documents\qpaper7\R Project Started 19 Dec 05\myFile.txt If I were to repeat this each time, my computer would get all bored and cranky and start to drop bits (only a joke, of course). I think I want to set the Home directory or the working directory or some directory or other to the above directory. I may or may not want to set some environmental variables. R 2.2.0; working directly from the console and copying and pasting code which I want to test into the console. Windows XP Home Edition. Administrator privileges are enabled. A curve ball: There are two accounts, Tom and Jones; the data are stored under Tom, whereas the computation is being done under the Jones account. I won't bore you with the details of why I am doing this. I was able to call Sys.getenv (R_USER) and get the home directory. I am a newbie to R and not familiar with the terminology. Tom Thomas L. Jones, Ph.D., Computer Science __ R-help@stat.math.ethz.ch 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] How to get around heteroscedasticity with non-linear leas t squares in R?
From: Brian S Cade Instead of thinking that the heteroscedasticity is a nuisance and something to get around, i.e, just wanting weighted estimates of the mean function, you might want to think about what heteroscedasticity is telling you and estimate some other quantities. Indeed! See Prof. Carroll's 2002 Fisher Lecture: http://www.stat.tamu.edu/ftp/pub/rjcarroll/2003.papers.directory/published_F isher_Lecture.pdf (There's Powerpoint file on his web page, too.) Andy Heteroscedasticity is telling you that the conditional distributions don't change at a constant rate across all portions of the distribution (think percentiles or more generally quantiles) and, therefore, a function for the mean (no matter how precisely estimated) cannot tell you all there is to know about your dose-response relation. Why not go after estimating the conditional quantile functions directly with nonlinear quantile regression, function nlrq() in the quantreg package? Brian Brian S. Cade U. S. Geological Survey Fort Collins Science Center 2150 Centre Ave., Bldg. C Fort Collins, CO 80526-8818 email: [EMAIL PROTECTED] tel: 970 226-9326 Kjetil Brinchmann Halvorsen [EMAIL PROTECTED] Sent by: [EMAIL PROTECTED] 02/21/2006 03:31 PM Please respond to [EMAIL PROTECTED] To Quin Wills [EMAIL PROTECTED] cc r-help@stat.math.ethz.ch Subject Re: [R] How to get around heteroscedasticity with non-linear least squares in R? Quin Wills wrote: I am using nls to fit dose-response curves but am not sure how to approach more robust regression in R to get around the problem of the my error showing increased variance with increasing dose. package sfsmisc has rnls (robust nls) which might be of use. Kjetil My understanding is that rlm or lqs would not be a good idea here. 'Fairly new to regression work, so apologies if I'm missing something obvious. [[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 [[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
