Re: [R] main title x title and y title with ggplot2
Guillaume, Have a look at the ggplot book on p. 29 (http://had.co.nz/ggplot2/book.pdf). HTH, Thierry ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 [EMAIL PROTECTED] www.inbo.be Do not put your faith in what statistics say until you have carefully considered what they do not say. ~William W. Watt A statistical analysis, properly conducted, is a delicate dissection of uncertainties, a surgery of suppositions. ~M.J.Moroney -Oorspronkelijk bericht- Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Namens guillaume chaumet Verzonden: woensdag 5 maart 2008 8:51 Aan: r-help@r-project.org Onderwerp: [R] main title x title and y title with ggplot2 Hi R people, I'm a R newbie and I'm trying to put main title, x title and y title in my graph with no success. Any idea? I'm sorry for this newbie question.. Guillaume [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Marginal Effects in a Logit Regression
Hi I am carrying out some logit regressions, so have a (0,1) dependent variable and various dependent variables, each in the (-inf, inf) range. As well as the usual output, I need to report the marginal effects of each the explanatory variables. How can this be done? I've seen a couple of similar questions on the mailing list but no answers. Thanks! Paul [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] main title x title and y title with ggplot2
Thierry First thank you for the celerity of your response. Second I use ggplot2 like this : ggplot(data, aes(x,y,fill)) + geom_point() + etc. Where did you your xlab and ylab when using ggplot2 like that? Guillaume 2008/3/5, ONKELINX, Thierry [EMAIL PROTECTED]: Guillaume, Have a look at the ggplot book on p. 29 (http://had.co.nz/ggplot2/book.pdf). HTH, Thierry ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 [EMAIL PROTECTED] www.inbo.be Do not put your faith in what statistics say until you have carefully considered what they do not say. ~William W. Watt A statistical analysis, properly conducted, is a delicate dissection of uncertainties, a surgery of suppositions. ~M.J.Moroney -Oorspronkelijk bericht- Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Namens guillaume chaumet Verzonden: woensdag 5 maart 2008 8:51 Aan: r-help@r-project.org Onderwerp: [R] main title x title and y title with ggplot2 Hi R people, I'm a R newbie and I'm trying to put main title, x title and y title in my graph with no success. Any idea? I'm sorry for this newbie question.. Guillaume [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] main title x title and y title with ggplot2
Guillaume, You'll have to add the appropriate scales. ggplot(data, aes(x,y,fill)) + geom_point() + scale_x_continuous(your xlabel) + scale_y_continuous(your ylabel) I suppose you can add a main title in a similar way, but I haven't found that yet. But I shure that Hadley will answer this. Thierry ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 [EMAIL PROTECTED] www.inbo.be Do not put your faith in what statistics say until you have carefully considered what they do not say. ~William W. Watt A statistical analysis, properly conducted, is a delicate dissection of uncertainties, a surgery of suppositions. ~M.J.Moroney -Oorspronkelijk bericht- Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Namens guillaume chaumet Verzonden: woensdag 5 maart 2008 9:45 Aan: r-help@r-project.org Onderwerp: Re: [R] main title x title and y title with ggplot2 Thierry First thank you for the celerity of your response. Second I use ggplot2 like this : ggplot(data, aes(x,y,fill)) + geom_point() + etc. Where did you your xlab and ylab when using ggplot2 like that? Guillaume 2008/3/5, ONKELINX, Thierry [EMAIL PROTECTED]: Guillaume, Have a look at the ggplot book on p. 29 (http://had.co.nz/ggplot2/book.pdf). HTH, Thierry ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 [EMAIL PROTECTED] www.inbo.be Do not put your faith in what statistics say until you have carefully considered what they do not say. ~William W. Watt A statistical analysis, properly conducted, is a delicate dissection of uncertainties, a surgery of suppositions. ~M.J.Moroney -Oorspronkelijk bericht- Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Namens guillaume chaumet Verzonden: woensdag 5 maart 2008 8:51 Aan: r-help@r-project.org Onderwerp: [R] main title x title and y title with ggplot2 Hi R people, I'm a R newbie and I'm trying to put main title, x title and y title in my graph with no success. Any idea? I'm sorry for this newbie question.. Guillaume [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] vertex labels in igraph from adjacency matrix
On Wed, Mar 05, 2008 at 02:27:21AM -0500, Charilaos Skiadas wrote: [...] Btw, you will likely want to take the betweenness call out, and call it once and store the result, instead of calling it twice (well, assuming the graph is largish). Or even better, use which.max: which.max(betweenness(graph = my.graph, v=V(my.graph), directed = FALSE)) This is almost good, but there is a catch, in igraph vertices are numbered from zero. So if you want an igraph vertex id, then you need to subtract one from this, i.e.: maxb - which.max(betweennness(my.graph, directed=FALSE))-1 You can double check it: betweenness(my.graph, maxb, directed=FALSE) Gabor PS. there is also an igraph mailing list, see the igraph homepage at igraph.sf.net Haris Skiadas Department of Mathematics and Computer Science Hanover College [...] -- Csardi Gabor [EMAIL PROTECTED]UNIL DGM __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] ipf function in R
Hi I have a 3 x 2 contingency table: 10 20 30 40 50 60 I want to update the frequencies to new marginal totals: 100 130 40 80 110 I want to use the ipf (iterative proportional fitting) function which is apparently in the cat package. Can somebody please advice me how to input this data and invoke ipf in R to obtain an updated contingency table? Thanks. By the way I am quite new to R. -- Dr Chandra Shah Senior Research Fellow Monash University-ACER Centre for the Economics of Education and Training Faculty of Education, Building 6, Monash University Victoria Australia 3800 Tel. +61 3 9905 2787 Fax +61 3 9905 9184 www.education.monash.edu.au/centres/ceet __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] coxme - fitting random treatment effect nested within centre
Dear all, I am using coxme function in Kinship library to fit random treatment effect nested within centre. I got 3 treatments (0,1,2) and 3 centres. I used following commands, but got an error. ugroup=paste(rep(1:3,each=3),rep(0:2,3),sep='/') mat1=bdsmatrix(rep(c(1,1,1,1,1,1,1,1,1),3),blocksize=rep(3,3),dimnames=list(ugroup,ugroup)) mat2=bdsmatrix(rep(c(0,0,0,0,0,0,0,0,1),3),blocksize=rep(3,3),dimnames=list(ugroup,ugroup)) group=paste(dat1$centre,dat1$treat,sep='/') coxme(Surv(time,status) ~ as.factor(treat), data=dat1,random= ~1|group,varlist=list(mat1,mat2),rescale=F,pdcheck=FALSE) Error in coxme.fit(X, Y, strats, offset, init, control, weights = weights, : Random effects variance is not spd Could anyone help me correcting this error? Many thanks in advance. Ruwanthi Looking for last minute shopping deals? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] vertex labels in igraph from adjacency matrix
Mark, graph.adjacency always preserves the order of the vertices, so the vertex at row/column 1 will be vertex #0 in the igraph graph, etc. I'll document this in a minute. This means that you can always do g - graph.adjacency(A) V(g)$name - colnames(A) But i completely agree that this should be done by default, i'll do that in the near future. Btw, weighted shortest path calculation (= where the edges have weights or capacities) is only implemented in the development version of igraph. Just in case you're looking for that. Best, Gabor On Wed, Mar 05, 2008 at 01:39:41AM -0500, Mark W Kimpel wrote: I am getting some unexpected results from some functions of igraph and it is possible that I am misinterpreting the vertex numbers. Eg., the max betweenness measure seems to be from a vertex that is not connected to a single other vertex. Below if my code snippet: require(igraph) my.graph - graph.adjacency(adjmatrix = my.adj.matrix, mode=c(undirected)) most.between.vert - colnames(my.adj.matrix)[which(betweenness(graph = my.graph, v=V(my.graph), directed = FALSE) == max(betweenness(graph = my.graph, v=V(my.graph), directed = FALSE))) + 1] sum(my.adj.matrix[most.between.vertext,]) 0 Is there a way to automatically set the vertex name attributes from the row and colnames of the inputted adjacency matrix to graph.adjacency? This would seem like an intuitive feature, but I can't find it documented. In the absence of this feature, how can I make sure that I am setting the vertex name attributes correctly? -- Mark W. Kimpel MD ** Neuroinformatics ** Dept. of Psychiatry Indiana University School of Medicine 15032 Hunter Court, Westfield, IN 46074 (317) 490-5129 Work, Mobile VoiceMail (317) 204-4202 Home (no voice mail please) mwkimpelatgmaildotcom ** -- Csardi Gabor [EMAIL PROTECTED]UNIL DGM __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] main title x title and y title with ggplot2
ONKELINX, Thierry Thierry.ONKELINX at inbo.be writes: Guillaume, You'll have to add the appropriate scales. ggplot(data, aes(x,y,fill)) + geom_point() + scale_x_continuous(your xlabel) + scale_y_continuous(your ylabel) I suppose you can add a main title in a similar way, but I haven't found that yet. But I shure that Hadley will answer this. Thierry Guillaume, I had the same problem and found a solution in some forums. Try this: p-ggplot(data, aes(x,y,fill)) + geom_point() + scale_x_continuous(your xlabel) + scale_y_continuous(your ylabel) p$title-your title print(p) Ingo. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] questions about p.adjust
Hi, I have some questions about p.adjust. The false discovery rate is a less stringent condition than the family wise error rate, so these methods are more powerful than the others., these methods refer to FDR methods or FWER methods. Simply what are the differences/pros/cons of both classes of methods. The 'BH' and 'BY' methods control the false discovery rate, are they FDR method ? fdr stands for false discovery rate, how different is it from BH or BY ? Lastly, which method should I use concerning microarray analysis for selecting differentially expressed genes. TIA [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] queries about bitmap()
Hi, There are different types of tiff methods in bitmap(), which one should be used for publication-quality pictures ? 'tiffcrle', 'tiffg3', 'tiffg32d', 'tiffg4', 'tifflzw', 'tiffpack', 'tiff12nc', 'tiff24nc', Thanks [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] R_alloc with structures with flexible array members
Dear All, In a package, I want to use some C code where I am using a structure (as the basic element of a linked list) with flexible array members. Basically, this is a structure where the last component is an incomplete array type (e.g., Harbison Steel, C, a reference manual, 5th ed., p. 159) such as: struct Sequence { struct Sequence *next; int len; unsigned int state_count[]; }; To create one such sequence, I allocate storage (following Harbison and Steel) in a C program as follows: struct Sequence *A; int n = 4; A = malloc( sizeof(struct Sequence) + n * sizeof(unsigned int)); If I understand correctly, however, it would be better to use R_alloc instead of malloc (memory automagically freed on exit and error; error-checking). But I do not know how to make the call to R_alloc here, since R_alloc allocates n units of size bytes each. I've tried, without success, the following two: int to_add_for_R_alloc = (int) ceil((float) sizeof(struct sequence) / sizeof(unsigned int)); A = (struct sequence *) R_alloc(to_add_for_R_alloc + n, sizeof(unsigned int)); or even a brute force attempt as: A = (struct sequence *) R_alloc( 100, sizeof(struct sequence)); but both result in segmentation faults. Should I just keep using malloc (and free at end)? Thanks, R. -- Ramon Diaz-Uriarte Statistical Computing Team Structural Biology and Biocomputing Programme Spanish National Cancer Centre (CNIO) http://ligarto.org/rdiaz __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] ipf function in R
On 05-Mar-08 07:14:28, Chandra Shah wrote: Hi I have a 3 x 2 contingency table: 10 20 30 40 50 60 I want to update the frequencies to new marginal totals: 100 130 40 80 110 I want to use the ipf (iterative proportional fitting) function which is apparently in the cat package. Can somebody please advice me how to input this data and invoke ipf in R to obtain an updated contingency table? Thanks. By the way I am quite new to R. It is not clear from your description what you really want to do. Hence, not clear whether you should be using ipf() to do it! The purpose of ipf() in the 'cat' package is very simple: You suppply a complete contingency table, and specify the order of interactions you want to include in the fit (using the margins parameter in ipf()); and then the value of ipf() is a table of the fitted (expected) values corresponding to a model of the type you have specified, fitted by maximum likelihood. Your task seems to be different. If I understand your statement. it looks as though you want to find a,b,c,d such that 10+a 20+c | 40 30+b 40+d | 80 5060 |110 -+--- 100 130 |230 or, since there are row and column constraints on a,b,c,d: 10+a 30-a | 40 40-a 40+a | 80 5060 |110 -+--- 100 130 |230 so then the task would be to find a. But now it depends on what you would mean by find a. You could treat 'a' as an unknown parameter, and fit it by maximum likelihood under various assumptions about the dependence between rows and columns (e.g. that the log odds-ratios have particular values). Or perhaps (and here is where ipf() may come in) for a given choice of 'a' you could submit the resulting table to ipf() and compute chisq from the resulting fitted values and the given observed values; then seek the value of 'a' which minimises chisq. However, it's not really possible to advise in more detail without being sure of what you are trying to do. Above, I am only guessing! Best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 05-Mar-08 Time: 10:26:57 -- XFMail -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] queries about bitmap()
Ng Stanley wrote: Hi, There are different types of tiff methods in bitmap(), which one should be used for publication-quality pictures ? I'd avoid bitmap at all and try vector formats such as PostScript or PDF or similar formats depending on your publication. Uwe Ligges 'tiffcrle', 'tiffg3', 'tiffg32d', 'tiffg4', 'tifflzw', 'tiffpack', 'tiff12nc', 'tiff24nc', Thanks [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Marginal Effects in a Logit Regression
On Wednesday 05 March 2008 09:25:11 am Paul Sweeting wrote: PS I am carrying out some logit regressions, so have a (0,1) dependent PS , I need to report the marginal effects of Have a look at lrm of the design package. It reports Dxy which is maybe what you want... Stefan __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] problem
Goodmorning Jim, My file has not only more than a million values, but more than a million rows and moreless 30 columns (it is a productive dataset for cows), infact with read.table i'm not able to import it. It is an xls file. How do you import your million rows and 4-5 columns files? thank you regards Dr.ssa Erika Frigo Università degli Studi di Milano Facoltà di Medicina Veterinaria Dipartimento di Scienze e Tecnologie Veterinarie per la Sicurezza Alimentare (VSA) Via Grasselli, 7 20137 Milano Tel. 02/50318515 Fax 02/50318501 - Original Message - From: jim holtman [EMAIL PROTECTED] To: Erika Frigo [EMAIL PROTECTED] Cc: r-help@r-project.org Sent: Tuesday, March 04, 2008 6:13 PM Subject: Re: [R] problem Is it just a file with a million values or is it some type of a structure with a million rows of indeterinent columns? If it is just a million numbers, you can easily read with is 'scan' or 'read.table' with no problem. I work with data structures that have several million rows and 4-5 columns without any problems. What is the format of the input? On 3/4/08, Erika Frigo [EMAIL PROTECTED] wrote: Good evening to everybody, I have problems to import in R a really big dataset (more than 100 values). Which is the best package to install? Is there someone who works with this kind of dataset and can help me, please? Thank you very much, Regards Dr.ssa Erika Frigo Department of Veterinary Sciences and Technology for Food Safety University of Milan Via Grasselli, 7 20137 Milano Tel. +39 0250318515 Fax +39 0250318501 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem you are trying to solve? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Principle component analysis
Thanks to Mr.Liviu Androvic and Mr.Richard Rowe helped me in PCA. Because I have just learn R language in a few day so I have many problem. 1) I don't know why PCA rotation function not run although I try many times. Would you please hepl me and explain how to read the PCA map (both of rotated and unrotated) in a concrete example. 2) Where I can find document relate: Plan S(A), S(A*B), S(A)*B? Thanks alot. -- View this message in context: http://www.nabble.com/Principle-component-analysis-tp15846902p15846902.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Spidergram
A parallel coordinate plot would do fine. Load the package iplots and then use the command ipcp(x1, x2,...) Antony Unwin __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R-Terminal
Thats really great, but now the sink() uses also this width. Is it possible to make sink using another width (i.e. 80 characters)? Thanks, Martin Am Dienstag, den 04.03.2008, 13:23 +0100 schrieb Martin Elff: On Tuesday 04 March 2008 (12:34:47), Peter Dalgaard wrote: Martin Kaffanke wrote: Hi there! I use an gnome-terminal for using R. When I resize the termial to the maximum size, R uses only the left side of the window. Can I tell R to use the whole window somehow? This seems to do it: options(width=Sys.getenv(COLUMNS)) (Surprisingly, at least to me, the COLUMNS environment variable appears to be automagically updated by gnome-terminal even after R is started. I was expecting to have to play with system(tset) and suchlike.) That seems to work also with KDE's Konsole. It may be worthwile considering to implement some automatic width-adjustment like: .adjustwidth - function(...){ options(width=Sys.getenv(COLUMNS)) TRUE } addTaskCallback(.adjustwidth) After that, the output width is adapted to the actual terminal width each time one hits return. Just my 20 cents. Best, Martin __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Ihr Partner für Webdesign, Webapplikationen und Webspace. http://www.roomandspace.com/ Martin Kaffanke +43 650 4514224 signature.asc Description: Dies ist ein digital signierter Nachrichtenteil __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Reversed but positive axis in trellis plots?
Hi, In my discpipline, it is common to plot one acoustic property on a positive scale but from top to bottom on the ordinate and the same for another measurement on the abscissa. So, the origin of the plot is on the top right of the plot, with increasing values to the left /down. This is to highlight the correlation between the acoustic measurement and the position of the forming structure, for instance when teaching it to students. The grouping ability of the trellis plot is quite handy whan plotting many instances of the same thing, so I was wondering if it is possible to make a trellis xyplot behave this way? Converting all values to negative and changing labels to the negative of the negative seems one solution to the reverseness of the axes, but how do I change the position? Is it possible? /Fredrik -- Give up learning, and put an end to your troubles. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] non-linear correlation
Hi all! This is a rather statistical question; Which effect sizes (parametric or not) could I use in order to estimate the amount of non-linear correlation between 2 variables? Is it possible to correct for auto-correlation within the correlated times series? Any suggestions for the appropriate packages/ functions are more than welcome!! Thank you! [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] t.test p-Value
Hello list, I am trying to apply the paired t.test between diseased and not diseased patients to identify genes that are more expressed in the one situation under the other. In order to retrieve the genes that are more expressed in the positive disease state I do: p.values-c() for(i in 1:length(Significant[,1])){ p.values[i]-try(t.test(positive[i,],negative[i,],alternative =greater)$p.value) } which(p.values0.01) where Significant is my matrix of genes and their expression in tumors and positive, negative are subsets of thes matrix. Whn p0.01, I reject the null hypothesis and I accept the alternative one, that I have greater gene expression in positive than in negative. I assume I must be doing sth wrong because the heatmap that I get with the genes that pass the filter of p-value is wrong. Could anyone help me with this? thanks a lot, Eleni [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] problem
On Wed, Mar 05, 2008 at 12:32:19PM +0100, Erika Frigo wrote: My file has not only more than a million values, but more than a million rows and moreless 30 columns (it is a productive dataset for cows), infact with read.table i'm not able to import it. It is an xls file. read.table() expects clear text -- e.g. csv or tab separated in the case of read.delim(). If your file is in xls format the simplest option would be to export the data to CSV format from Excel. If for some reason that is not an option please have a look at the R Data Import/Export manual. Of course neither will solve the problem of not enough memory if your file is simply too large. In that case you will may want to put your data into a database and have R connect to it and retrieve the data in smaller chunks as required. cu Philipp -- Dr. Philipp Pagel Tel. +49-8161-71 2131 Lehrstuhl für Genomorientierte Bioinformatik Fax. +49-8161-71 2186 Technische Universität München Wissenschaftszentrum Weihenstephan 85350 Freising, Germany and Institut für Bioinformatik und Systembiologie / MIPS Helmholtz Zentrum München - Deutsches Forschungszentrum für Gesundheit und Umwelt Ingolstädter Landstrasse 1 85764 Neuherberg, Germany http://mips.gsf.de/staff/pagel __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R-Terminal
On Wednesday 05 March 2008 (12:56:17), Martin Kaffanke wrote: Thats really great, but now the sink() uses also this width. Is it possible to make sink using another width (i.e. 80 characters)? # auto width adjustment .adjustWidth - function(...){ options(width=Sys.getenv(COLUMNS)) TRUE } .adjustWidthCallBack - addTaskCallback(.adjustWidth) # a wrapper for 'sink' Sink - function(...,width=80){ if(nargs()){ removeTaskCallback(.adjustWidthCallBack) rm(.adjustWidthCallBack,envir=globalenv()) options(width=width) sink(...) } else{ sink() assign(.adjustWidthCallBack, addTaskCallback(.adjustWidth), envir=globalenv() ) } } # test Sink(test.out,width=20) # nonsensically narrow, but shows that # width=something has an effect... iris[1:10,] Sink() as.matrix(readLines(test.out)) iris[1:10,] Best, Martin - Dr. Martin Elff Faculty of Social Sciences LSPWIVS (van Deth) University of Mannheim A5, 6 68131 Mannheim Germany Phone: +49-621-181-2093 Fax: +49-621-181-2099 E-Mail: [EMAIL PROTECTED] Web: http://webrum.uni-mannheim.de/sowi/elff/ http://www.sowi.uni-mannheim.de/lspwivs/ - __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] t.test p-Value
I am sorry, the test is unpaired...But my question remains Thanks, Eleni On Wed, Mar 5, 2008 at 2:33 PM, Eleni Christodoulou [EMAIL PROTECTED] wrote: Hello list, I am trying to apply the paired t.test between diseased and not diseased patients to identify genes that are more expressed in the one situation under the other. In order to retrieve the genes that are more expressed in the positive disease state I do: p.values-c() for(i in 1:length(Significant[,1])){ p.values[i]-try(t.test(positive[i,],negative[i,],alternative =greater)$p.value) } which(p.values0.01) where Significant is my matrix of genes and their expression in tumors and positive, negative are subsets of thes matrix. Whn p0.01, I reject the null hypothesis and I accept the alternative one, that I have greater gene expression in positive than in negative. I assume I must be doing sth wrong because the heatmap that I get with the genes that pass the filter of p-value is wrong. Could anyone help me with this? thanks a lot, Eleni [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] matrix inversion using solve() and matrices containing large/small values
Hello I've stumbled upon a problem for inversion of a matrix with large values, and I haven't found a solution yet... I wondered if someone could give a hand. (It is about automatic optimisation of a calibration process, which involves the inverse of the information matrix) code: * macht=0.8698965 coeff=1.106836*10^(-8) echtecoeff=c(481.46,19919.23,-93.41188,0.5939589,-0.002846272,8.030726e-6 ,-1.155094e-8,6.357603e-12)/1000 dosis=c(0,29,70,128,201,290,396) dfdb - array(c(1,1,1,1,1,1,1,dosis,dosis^2,dosis^3,dosis^4,dosis^5,dosis^6,dosis^7),dim=c(7,8)) dfdbtrans = aperm(dfdb) sigerr=sqrt(coeff*dosis^macht) sigmadosis = c(1:7) for(i in 1:7){ sigmadosis[i]=ifelse(sigerr[i]2.257786084*10^(-4),2.257786084*10^(-4),sigerr[i]) } omega = diag(sigmadosis) infomatrix = dfdbtrans%*%omega%*%dfdb ** I need the inverse of this information matrix, and infomatrix_inv = solve(infomatrix, tol = 10^(-43)) does not deliver adequate results (matrixproduct of infomatrix_inv and infomatrix is not 1). Regular use of solve() delivers the error 'system is computationally singular: reciprocal condition number: 2.949.10^(-41)' So I went over to an eigendecomposition using eigen(): (so that infomatrix = V D V^(-1) == infomatrix^(-1)= V D^(-1) V^(-1) ) in the hope this would deliver better results.) *** infomatrix_eigen = eigen(infomatrix) infomatrix_eigen_D = diag(infomatrix_eigen$values) infomatrix_eigen_V = infomatrix_eigen$vectors infomatrix_eigen_V_inv = solve(infomatrix_eigen_V) *** however, the matrix product of these are not the same as the infomatrix itself, only in certain parts: infomatrix_eigen_V %*% infomatrix_eigen_D %*% infomatrix_eigen_V_inv infomatrix Therefore, I reckon the inverse of eigendecomposition won't be correct either. As far as I understand, the problem is due to the very large range of data, and therefore results in numerical problems, but I can't come up with a way to do it otherwise. Would anyone know how I could solve this problem? (PS, i'm running under linux suse 10.0, latest R version with MASS libraries (RV package)) F. Crop UGent -- Medical Physics [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] main title x title and y title with ggplot2
I had the same problem and found a solution in some forums. Try this: p-ggplot(data, aes(x,y,fill)) + geom_point() + scale_x_continuous(your xlabel) + scale_y_continuous(your ylabel) The new (more ggplot-like way) is to do: ggplot(data, aes(x,y,fill)) + ... + opts(title = my title) Hadley -- http://had.co.nz/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] matrix inversion using solve() and matrices containing large/small values
On 3/5/2008 8:21 AM, gerardus vanneste wrote: Hello I've stumbled upon a problem for inversion of a matrix with large values, and I haven't found a solution yet... I wondered if someone could give a hand. (It is about automatic optimisation of a calibration process, which involves the inverse of the information matrix) If the matrix actually isn't singular, then a rescaling of the parameters should help a lot. I see the diagonal of infomatrix as diag(infomatrix) [1] 5.930720e-03 3.872339e+02 4.562529e+07 6.281634e+12 9.228140e+17 [6] 1.398687e+23 2.154124e+28 3.345598e+33 so multiplying the parameters by numbers on the order of the square roots of these entries (e.g. 10^c(-1, 1, 4, 6, 9, 12, 14, 17)), and redoing the rest of the calculations on that scale, should work. Duncan Murdoch code: * macht=0.8698965 coeff=1.106836*10^(-8) echtecoeff=c(481.46,19919.23,-93.41188,0.5939589,-0.002846272,8.030726e-6 ,-1.155094e-8,6.357603e-12)/1000 dosis=c(0,29,70,128,201,290,396) dfdb - array(c(1,1,1,1,1,1,1,dosis,dosis^2,dosis^3,dosis^4,dosis^5,dosis^6,dosis^7),dim=c(7,8)) dfdbtrans = aperm(dfdb) sigerr=sqrt(coeff*dosis^macht) sigmadosis = c(1:7) for(i in 1:7){ sigmadosis[i]=ifelse(sigerr[i]2.257786084*10^(-4),2.257786084*10^(-4),sigerr[i]) } omega = diag(sigmadosis) infomatrix = dfdbtrans%*%omega%*%dfdb ** I need the inverse of this information matrix, and infomatrix_inv = solve(infomatrix, tol = 10^(-43)) does not deliver adequate results (matrixproduct of infomatrix_inv and infomatrix is not 1). Regular use of solve() delivers the error 'system is computationally singular: reciprocal condition number: 2.949.10^(-41)' So I went over to an eigendecomposition using eigen(): (so that infomatrix = V D V^(-1) == infomatrix^(-1)= V D^(-1) V^(-1) ) in the hope this would deliver better results.) *** infomatrix_eigen = eigen(infomatrix) infomatrix_eigen_D = diag(infomatrix_eigen$values) infomatrix_eigen_V = infomatrix_eigen$vectors infomatrix_eigen_V_inv = solve(infomatrix_eigen_V) *** however, the matrix product of these are not the same as the infomatrix itself, only in certain parts: infomatrix_eigen_V %*% infomatrix_eigen_D %*% infomatrix_eigen_V_inv infomatrix Therefore, I reckon the inverse of eigendecomposition won't be correct either. As far as I understand, the problem is due to the very large range of data, and therefore results in numerical problems, but I can't come up with a way to do it otherwise. Would anyone know how I could solve this problem? (PS, i'm running under linux suse 10.0, latest R version with MASS libraries (RV package)) F. Crop UGent -- Medical Physics [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] matrix inversion using solve() and matrices containing large/small values
Sorry, I meant to send this to the whole list. On Mar 5, 2008, at 8:46 AM, Charilaos Skiadas wrote: The problem doesn't necessarily have to do with the range of data. At first level, it has to do with the simple fact that dfdb has rank 6 at most, (7 at most in general, though in your case since 290=10*29, it is at most 6). Since matrix rank only goes down when multiplying, your infomatrix is an 8x8 matrix, with rank at most 6 (7 if you were more lucky with that 290, still not good enough), so it is certainly not invertible, even forgetting the computational issues of computing the inverse. You would need either power smaller than 7, or a longer dosis vector, I think. If you manage to make your problem in a case where dfdb is square, then you just have to invert that, which might be easier, seeing how it's a Vandermonde matrix. Haris Skiadas Department of Mathematics and Computer Science Hanover College On Mar 5, 2008, at 8:21 AM, gerardus vanneste wrote: Hello I've stumbled upon a problem for inversion of a matrix with large values, and I haven't found a solution yet... I wondered if someone could give a hand. (It is about automatic optimisation of a calibration process, which involves the inverse of the information matrix) code: * macht=0.8698965 coeff=1.106836*10^(-8) echtecoeff=c (481.46,19919.23,-93.41188,0.5939589,-0.002846272,8.030726e-6 ,-1.155094e-8,6.357603e-12)/1000 dosis=c(0,29,70,128,201,290,396) dfdb - array(c (1,1,1,1,1,1,1,dosis,dosis^2,dosis^3,dosis^4,dosis^5,dosis^6,dosis^7) ,dim=c(7,8)) dfdbtrans = aperm(dfdb) sigerr=sqrt(coeff*dosis^macht) sigmadosis = c(1:7) for(i in 1:7){ sigmadosis[i]=ifelse(sigerr[i]2.257786084*10^(-4),2.257786084*10^ (-4),sigerr[i]) } omega = diag(sigmadosis) infomatrix = dfdbtrans%*%omega%*%dfdb ** I need the inverse of this information matrix, and infomatrix_inv = solve(infomatrix, tol = 10^(-43)) does not deliver adequate results (matrixproduct of infomatrix_inv and infomatrix is not 1). Regular use of solve() delivers the error 'system is computationally singular: reciprocal condition number: 2.949.10^ (-41)' So I went over to an eigendecomposition using eigen(): (so that infomatrix = V D V^(-1) == infomatrix^(-1)= V D^(-1) V^(-1) ) in the hope this would deliver better results.) *** infomatrix_eigen = eigen(infomatrix) infomatrix_eigen_D = diag(infomatrix_eigen$values) infomatrix_eigen_V = infomatrix_eigen$vectors infomatrix_eigen_V_inv = solve(infomatrix_eigen_V) *** however, the matrix product of these are not the same as the infomatrix itself, only in certain parts: infomatrix_eigen_V %*% infomatrix_eigen_D %*% infomatrix_eigen_V_inv infomatrix Therefore, I reckon the inverse of eigendecomposition won't be correct either. As far as I understand, the problem is due to the very large range of data, and therefore results in numerical problems, but I can't come up with a way to do it otherwise. Would anyone know how I could solve this problem? (PS, i'm running under linux suse 10.0, latest R version with MASS libraries (RV package)) F. Crop UGent -- Medical Physics [[alternative HTML version deleted]] Haris Skiadas Department of Mathematics and Computer Science Hanover College __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] t.test p-Value
On Wed, Mar 5, 2008 at 2:05 PM, ian white [EMAIL PROTECTED] wrote: Don't you need to make some allowance for multiple testing? E.g. to get a experiment-wise significance level of 0.01 you need which(p.values very small number) where the very small number is approximately 0.01/(total number of genes). On Wed, 2008-03-05 at 14:38 +0200, Eleni Christodoulou wrote: I am sorry, the test is unpaired...But my question remains Thanks, Eleni On Wed, Mar 5, 2008 at 2:33 PM, Eleni Christodoulou [EMAIL PROTECTED] wrote: Hello list, I am trying to apply the paired t.test between diseased and not diseased patients to identify genes that are more expressed in the one situation under the other. In order to retrieve the genes that are more expressed in the positive disease state I do: p.values-c() for(i in 1:length(Significant[,1])){ p.values[i]-try(t.test(positive[i,],negative[i,],alternative =greater)$p.value) } which(p.values0.01) where Significant is my matrix of genes and their expression in tumors and positive, negative are subsets of thes matrix. Whn p0.01, I reject the null hypothesis and I accept the alternative one, that I have greater gene expression in positive than in negative. I assume I must be doing sth wrong because the heatmap that I get with the genes that pass the filter of p-value is wrong. Could anyone help me with this? thanks a lot, Eleni [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] nls: different results if applied to normal or linearized data
Dear all, I did a non-linear least square model fit y ~ a * x^b (a) nls(y ~ a * x^b, start=list(a=1,b=1)) to obtain the coefficients a b. I did the same with the linearized formula, including a linear model log(y) ~ log(a) + b * log(x) (b) nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1)) (c) lm(log10(y) ~ log10(x)) I expected coefficient b to be identical for all three cases. Hoever, using my dataset, coefficient b was: (a) 0.912 (b) 0.9794 (c) 0.9794 Coefficient a also varied between option (a) and (b), 107.2 and 94.7, respectively. Is this supposed to happen? Which is the correct coefficient b? Regards, Wolfgang -- Laboratory of Animal Physiology Department of Biology University of Turku FIN-20014 Turku Finland __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Asking, are simple effects different from 0
On 3/4/2008 2:45 PM, Jarrett Byrnes wrote: Hello, R-i-zens. I'm working on an data set with a factorial ANOVA that has a significant interaction. I'm interested in seeing whether the simple effects are different from 0, and I'm pondering how to do this. So, I have my.anova-lm(response ~ trtA*trtB) The output for which gives me a number of coefficients and whether they are different from 0. However, I want the simple effects only, incorporating their intercepts, with their error estimates. Can I somehow manipulate this object to get that? Or, would a good shortcut be my.simple.anova-lm(response ~ trtA:trtB + 0) and then use those coefficients and their error estimates? One approach would be to use glht() in the multcomp package. You need to work out how to formulate the matrix of coefficients that give the desired contrasts. Here is an example using the warpbreaks data frame: fm - lm(breaks ~ tension*wool, data=warpbreaks) # names(coef(fm)) # (Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB cm - rbind( A vs. B at L = c(0, 0, 0,-1, 0, 0), A vs. B at M = c(0, 0, 0,-1,-1, 0), A vs. B at H = c(0, 0, 0,-1, 0,-1), M vs. L at A = c(0, 1, 0, 0, 0, 0), M vs. H at A = c(0, 1,-1, 0, 0, 0), L vs. H at A = c(0, 0,-1, 0, 0, 0), M vs. L at B = c(0, 1, 0, 0, 1, 0), M vs. H at B = c(0, 1,-1, 0, 1,-1), L vs. H at B = c(0, 0,-1, 0, 0,-1)) library(multcomp) summary(glht(fm, linfct = cm), test = adjusted(type=none)) Simultaneous Tests for General Linear Hypotheses Fit: lm(formula = breaks ~ tension * wool, data = warpbreaks) Linear Hypotheses: Estimate Std. Error t value p value A vs. B at L == 0 16. 5.1573 3.167 0.002677 ** A vs. B at M == 0 -4.7778 5.1573 -0.926 0.358867 A vs. B at H == 0 5.7778 5.1573 1.120 0.268156 M vs. L at A == 0 -20.5556 5.1573 -3.986 0.000228 *** M vs. H at A == 0 -0.5556 5.1573 -0.108 0.914665 L vs. H at A == 0 20. 5.1573 3.878 0.000320 *** M vs. L at B == 0 0.5556 5.1573 0.108 0.914665 M vs. H at B == 0 10. 5.1573 1.939 0.058392 . L vs. H at B == 0 9. 5.1573 1.831 0.073270 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- none method) If so, I realize that R gives me t tests for each coefficient. Now, for those I know I'm using the residual degrees of freedom. Would it then be more appropriate to use those, all with the same residual DF and apply a bonferroni correction, or, use the mean and SE estimate with the sample size for that particular treatment and perform an uncorrected one sample t-test to see if the value is different from 0? I won't comment on whether to adjust and if so how, but you can implement various adjustments when summarizing. For example: summary(glht(fm, linfct = cm), test = adjusted(type=bonferroni)) Simultaneous Tests for General Linear Hypotheses Fit: lm(formula = breaks ~ tension * wool, data = warpbreaks) Linear Hypotheses: Estimate Std. Error t value p value A vs. B at L == 0 16. 5.1573 3.167 0.02409 * A vs. B at M == 0 -4.7778 5.1573 -0.926 1.0 A vs. B at H == 0 5.7778 5.1573 1.120 1.0 M vs. L at A == 0 -20.5556 5.1573 -3.986 0.00205 ** M vs. H at A == 0 -0.5556 5.1573 -0.108 1.0 L vs. H at A == 0 20. 5.1573 3.878 0.00288 ** M vs. L at B == 0 0.5556 5.1573 0.108 1.0 M vs. H at B == 0 10. 5.1573 1.939 0.52553 L vs. H at B == 0 9. 5.1573 1.831 0.65943 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- bonferroni method) Sorry for the bonehead question, but it's a situation I haven't seen before. -Jarrett __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Chuck Cleland, Ph.D. NDRI, Inc. 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] (no subject)
Dear all, I try to import a SPSS.por dataset with about 6000 cases and 650 variables with R commander and got error messages: 'Error: /temp/t.por use.value.labels=true, max.value.label=inf, 'Error: data is not data frame and cannot be attached' Any comment? Thanks in advance. wz http://geocities.com/wendyzhao78/HW.html [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls: different results if applied to normal or linearized data
Write out the objective functions that they are minimizing and it will be clear they are different so you can't expect the same results. On Wed, Mar 5, 2008 at 8:53 AM, Wolfgang Waser [EMAIL PROTECTED] wrote: Dear all, I did a non-linear least square model fit y ~ a * x^b (a) nls(y ~ a * x^b, start=list(a=1,b=1)) to obtain the coefficients a b. I did the same with the linearized formula, including a linear model log(y) ~ log(a) + b * log(x) (b) nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1)) (c) lm(log10(y) ~ log10(x)) I expected coefficient b to be identical for all three cases. Hoever, using my dataset, coefficient b was: (a) 0.912 (b) 0.9794 (c) 0.9794 Coefficient a also varied between option (a) and (b), 107.2 and 94.7, respectively. Is this supposed to happen? Which is the correct coefficient b? Regards, Wolfgang -- Laboratory of Animal Physiology Department of Biology University of Turku FIN-20014 Turku Finland __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Updating packages from one hard-drive to another, after upgrade of R
The CPU on my computer 'died' and I have had to purchase a new computer. I have just installed v 2.6.2. My previous computer had v 2.5.1 and a large number of files in the library folder. These files have been copied to a partition on my new hard drive, along with the old R installation. Can I just copy all the folders/files in the old library to the new v2.6.2 library or do the packages need to be downloaded again? There are quite a number of folders so it would be time-consuming to have to do a fresh download. Any assistance is appreciated, Bob __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Import SPSS.por file
Dear all, I try to import a SPSS.por dataset with about 6000 cases and 650 variables with R commander and got error messages: 'Error: /temp/t.por use.value.labels=true, max.value.label=inf, 'Error: data is not data frame and cannot be attached' Any comment? Thanks in advance. wz http://geocities.com/wendyzhao78/HW.html [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] coxme - fitting random treatment effect nested within centre
included message Dear all, I am using coxme function in Kinship library to fit random treatment effect nested within centre. I got 3 treatments (0,1,2) and 3 centres. I used following commands, but got an error. ugroup=paste(rep(1:3,each=3),rep(0:2,3),sep='/') mat1=bdsmatrix(rep(c(1,1,1,1,1,1,1,1,1),3),blocksize=rep(3,3),dimnames=list(ugro up,ugroup)) mat2=bdsmatrix(rep(c(0,0,0,0,0,0,0,0,1),3),blocksize=rep(3,3),dimnames=list(ugro up,ugroup)) group=paste(dat1$centre,dat1$treat,sep='/') coxme(Surv(time,status) ~ as.factor(treat), data=dat1,random= ~1|group,varlist=list(mat1,mat2),rescale=F,pdcheck=FALSE) Error in coxme.fit(X, Y, strats, offset, init, control, weights = weights, : Random effects variance is not spd Could anyone help me correcting this error? Many thanks in advance. Ruwanthi --- end inclusion - The build your own bdsmatrix style in coxme is for variance structures that are not built in, like pedigree data. Your problem is one that coxme can do directly, so it is easier to call the routine simply: fit - coxme(Surv(time, status) ~ factor(treat), data=data1, random = ~1 | centre/treat) Terry Therneau __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] legend for several graphics
Hi, I am trying to generate a figure of 9 plots that are contained in one device by using par(mfrow = c(3,3,)) I would like to have 1 common legend for all 9 plots somewhere outside of the plotting area (as opposed to one legend inside each of the 9 plots, which the function legend() seems to generate by default). Any hint how to do this? Best, Georg __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls: different results if applied to normal or linearized data
On Wednesday 05 March 2008 (14:53:27), Wolfgang Waser wrote: Dear all, I did a non-linear least square model fit y ~ a * x^b (a) nls(y ~ a * x^b, start=list(a=1,b=1)) to obtain the coefficients a b. I did the same with the linearized formula, including a linear model log(y) ~ log(a) + b * log(x) (b) nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1)) (c) lm(log10(y) ~ log10(x)) I expected coefficient b to be identical for all three cases. Hoever, using my dataset, coefficient b was: (a) 0.912 (b) 0.9794 (c) 0.9794 Coefficient a also varied between option (a) and (b), 107.2 and 94.7, respectively. Models (a) and (b) entail different distributions of the dependent variable y and different ranges of values that y may take. (a) implies that y has, conditionally on x, a normal distribution and has a range of feasible values from -Inf to +Inf. (b) and (c) imply that log(y) has a normal distribution, that is, y has a log-normal distribution and can take values from zero to +Inf. Is this supposed to happen? Given the above considerations, different results with respect to the intercept are definitely to be expected. Which is the correct coefficient b? That depends - is y strictly non-negative or not ... Just my 20 cents... __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Updating packages from one hard-drive to another, after upgrade of R
Bob, You can copy the files from the packages to your new computer. Then run update.packages(checkBuilt = TRUE). That should do. HTH, Thierry ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 [EMAIL PROTECTED] www.inbo.be Do not put your faith in what statistics say until you have carefully considered what they do not say. ~William W. Watt A statistical analysis, properly conducted, is a delicate dissection of uncertainties, a surgery of suppositions. ~M.J.Moroney -Oorspronkelijk bericht- Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Namens Bob Green Verzonden: woensdag 5 maart 2008 15:07 Aan: r-help@r-project.org Onderwerp: [R] Updating packages from one hard-drive to another,after upgrade of R The CPU on my computer 'died' and I have had to purchase a new computer. I have just installed v 2.6.2. My previous computer had v 2.5.1 and a large number of files in the library folder. These files have been copied to a partition on my new hard drive, along with the old R installation. Can I just copy all the folders/files in the old library to the new v2.6.2 library or do the packages need to be downloaded again? There are quite a number of folders so it would be time-consuming to have to do a fresh download. Any assistance is appreciated, Bob __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Updating packages from one hard-drive to another, after upgrade of R
See ?.libPaths for info on setting the location of your libraries. If you do want to copy them the copydir.bat utility in batchfiles.googlecode.com can do that. On Wed, Mar 5, 2008 at 9:07 AM, Bob Green [EMAIL PROTECTED] wrote: The CPU on my computer 'died' and I have had to purchase a new computer. I have just installed v 2.6.2. My previous computer had v 2.5.1 and a large number of files in the library folder. These files have been copied to a partition on my new hard drive, along with the old R installation. Can I just copy all the folders/files in the old library to the new v2.6.2 library or do the packages need to be downloaded again? There are quite a number of folders so it would be time-consuming to have to do a fresh download. Any assistance is appreciated, Bob __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] matrix inversion using solve() and matrices containing large/small values
On Wed, Mar 5, 2008 at 7:43 AM, Duncan Murdoch [EMAIL PROTECTED] wrote: On 3/5/2008 8:21 AM, gerardus vanneste wrote: Hello I've stumbled upon a problem for inversion of a matrix with large values, and I haven't found a solution yet... Someone with experience in numerical linear algebra will immediately focus on the words inversion of a matrix in that statement. There is a sort of a meta-result in numerical linear algebra that if the best way you can formulate an answer to your problem involves inversion of a matrix (without specifying a special structure like diagonal or bidiagonal or unit triangular or ...) then you need to think about your problem in more detail. Certainly a formula may be written out in terms of the inverse of an information matrix but it is a bad idea to try to calculate the result that way. If you have an information matrix is it positive definite? If so, you should be working with the Cholesky decomposition of the matrix or perhaps the QR decomposition of a model matrix that generates the information matrix. Think in terms of the factors of a matrix and how you can reexpress the calculation using them. I wondered if someone could give a hand. (It is about automatic optimisation of a calibration process, which involves the inverse of the information matrix) If the matrix actually isn't singular, then a rescaling of the parameters should help a lot. I see the diagonal of infomatrix as diag(infomatrix) [1] 5.930720e-03 3.872339e+02 4.562529e+07 6.281634e+12 9.228140e+17 [6] 1.398687e+23 2.154124e+28 3.345598e+33 so multiplying the parameters by numbers on the order of the square roots of these entries (e.g. 10^c(-1, 1, 4, 6, 9, 12, 14, 17)), and redoing the rest of the calculations on that scale, should work. Duncan Murdoch code: * macht=0.8698965 coeff=1.106836*10^(-8) echtecoeff=c(481.46,19919.23,-93.41188,0.5939589,-0.002846272,8.030726e-6 ,-1.155094e-8,6.357603e-12)/1000 dosis=c(0,29,70,128,201,290,396) dfdb - array(c(1,1,1,1,1,1,1,dosis,dosis^2,dosis^3,dosis^4,dosis^5,dosis^6,dosis^7),dim=c(7,8)) dfdbtrans = aperm(dfdb) sigerr=sqrt(coeff*dosis^macht) sigmadosis = c(1:7) for(i in 1:7){ sigmadosis[i]=ifelse(sigerr[i]2.257786084*10^(-4),2.257786084*10^(-4),sigerr[i]) } omega = diag(sigmadosis) infomatrix = dfdbtrans%*%omega%*%dfdb ** I need the inverse of this information matrix, and infomatrix_inv = solve(infomatrix, tol = 10^(-43)) does not deliver adequate results (matrixproduct of infomatrix_inv and infomatrix is not 1). Regular use of solve() delivers the error 'system is computationally singular: reciprocal condition number: 2.949.10^(-41)' So I went over to an eigendecomposition using eigen(): (so that infomatrix = V D V^(-1) == infomatrix^(-1)= V D^(-1) V^(-1) ) in the hope this would deliver better results.) *** infomatrix_eigen = eigen(infomatrix) infomatrix_eigen_D = diag(infomatrix_eigen$values) infomatrix_eigen_V = infomatrix_eigen$vectors infomatrix_eigen_V_inv = solve(infomatrix_eigen_V) *** however, the matrix product of these are not the same as the infomatrix itself, only in certain parts: infomatrix_eigen_V %*% infomatrix_eigen_D %*% infomatrix_eigen_V_inv infomatrix Therefore, I reckon the inverse of eigendecomposition won't be correct either. As far as I understand, the problem is due to the very large range of data, and therefore results in numerical problems, but I can't come up with a way to do it otherwise. Would anyone know how I could solve this problem? (PS, i'm running under linux suse 10.0, latest R version with MASS libraries (RV package)) F. Crop UGent -- Medical Physics [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] regex sulotion for seperating number and string
I have strings contain postcode and letters, some seperated with blank, some with comma, and some hasn't seperated. eg, 2324gz 2567 HK 3741,BF I want to seperate the number and letters into two new variables. I know this should be quite basic question, but searched on regex syntax and that seems a bit scarey to me, any one can shot me a quick solution on this particular question? thanks, Sun __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] vertex labels in igraph from adjacency matrix
Looks like I turned an off my one error into an off by two error by adding rather than subtracting. Clearly a logic error on my part. Also, which.max is clearly superior as it results in half as many function calls. Thanks guys! As an aside, although igraph may use the C indexing convention, R users and their code, is much more in tune with indexing beginning at 1. Adjusting this in the R interface to igraph may solve many headaches down the road. Mark Mark W. Kimpel MD ** Neuroinformatics ** Dept. of Psychiatry Indiana University School of Medicine 15032 Hunter Court, Westfield, IN 46074 (317) 490-5129 Work, Mobile VoiceMail (317) 204-4202 Home (no voice mail please) mwkimpelatgmaildotcom ** Gabor Csardi wrote: On Wed, Mar 05, 2008 at 02:27:21AM -0500, Charilaos Skiadas wrote: [...] Btw, you will likely want to take the betweenness call out, and call it once and store the result, instead of calling it twice (well, assuming the graph is largish). Or even better, use which.max: which.max(betweenness(graph = my.graph, v=V(my.graph), directed = FALSE)) This is almost good, but there is a catch, in igraph vertices are numbered from zero. So if you want an igraph vertex id, then you need to subtract one from this, i.e.: maxb - which.max(betweennness(my.graph, directed=FALSE))-1 You can double check it: betweenness(my.graph, maxb, directed=FALSE) Gabor PS. there is also an igraph mailing list, see the igraph homepage at igraph.sf.net Haris Skiadas Department of Mathematics and Computer Science Hanover College [...] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] problem with geepack
Hi all I am analyzing a data set containing information about the behaviour of marine molluscs on a vertical wall. Since I have replicate observations on the same individuals I was thinking to use the geepack library. The data are organised in a dataframe with the following variables Date = date of sampling, Size = dimensions (mm) Activity duration of activity (min) Water = duration of splashing by waves Hgt = resting eight of each specimen before activity begin Individual = a code indicating the id of the specimen. I have up to 12 replicate observations for individual. Some observation are missing and I organized the data frame to have exactly 12 rows for each specimen, with NAs where there is a missing observation. The following model worked fine: gee1-geese(Activity~Water, id=Individual, data=dataF, family=gaussian) but when I use other variables e.g gee2-geese(Activity~Hgt+Size+Water, id=Individual, data=dataF, family=gaussian) I get the error message Error in geese.fit(x, y, id, offset, soffset, w, waves, zsca, zcor, corp, : nrow(zsca) and length(y) not match which I am not able to understand. The same problem has been reported in the list in 2006, but I have not found any response to it. Any suggestion? Giacomo __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Constrained regression
I would like to acknowledge the answers I received from Tom Filloon, Mike Cheung and Berwyn Turlach. Berwyn's response was exactly what I needed. Use solve.QP from the quadprog package in R. S-Plus has the equivalent function solveQP in the NuOpt module. Berwyn's response is below G'day Carlos, On Mon, Mar 3, 2008 at 11:52 AM Carlos Alzola [EMAIL PROTECTED] wrote: I am trying to get information on how to fit a linear regression with constrained parameters. Specifically, I have 8 predictors , their coeffiecients should all be non-negative and add up to 1. I understand it is a quadratic programming problem but I have no experience in the subject. I searched the archives but the results were inconclusive. Could someone provide suggestions and references to the literature, please? A suggestion: library(MASS) ## to access the Boston data designmat - model.matrix(medv~., data=Boston) Dmat - crossprod(designmat, designmat) dvec - crossprod(designmat, Boston$medv) Amat - cbind(1, diag(NROW(Dmat))) bvec - c(1, rep(0,NROW(Dmat)) meq - 1 library(quadprog) res - solve.QP(Dmat, dvec, Amat, bvec, meq) The solution seems to contain values that are, for all practical purposes, actually zero: res$solution [1] 4.535581e-16 2.661931e-18 1.016929e-01 -1.850699e-17 [5] 1.458219e-16 -3.892418e-15 8.544939e-01 0.00e+00 [9] 2.410742e-16 2.905722e-17 -5.700600e-20 -4.227261e-17 [13] 4.381328e-02 -3.723065e-18 So perhaps better: zapsmall(res$solution) [1] 0.000 0.000 0.1016929 0.000 0.000 0.000 [7] 0.8544939 0.000 0.000 0.000 0.000 0.000 [13] 0.0438133 0.000 So the estimates seem to follow the constraints. And the unconstrained solution is: res$unconstrainted.solution [1] 3.645949e+01 -1.080114e-01 4.642046e-02 2.055863e-02 [5] 2.686734e+00 -1.776661e+01 3.809865e+00 6.922246e-04 [9] -1.475567e+00 3.060495e-01 -1.233459e-02 -9.527472e-01 [13] 9.311683e-03 -5.247584e-01 which seems to coincide with what lm() thinks it should be: coef(lm(medv~., Boston)) (Intercept) crimzn indus chas 3.645949e+01 -1.080114e-01 4.642046e-02 2.055863e-02 2.686734e+00 noxrm age dis rad -1.776661e+01 3.809865e+00 6.922246e-04 -1.475567e+00 3.060495e-01 tax ptratio black lstat -1.233459e-02 -9.527472e-01 9.311683e-03 -5.247584e-01 So there seem to be no numeric problems. Otherwise we could have done something else (e.g calculate the QR factorization of the design matrix, say X, and give the R factor to solve.QP, instead of calculating X'X and giving that one to solve.QP). If the intercept is not supposed to be included in the set of constrained estimates, then something like the following can be done: Amat[1,] - 0 res - solve.QP(Dmat, dvec, Amat, bvec, meq) zapsmall(res$solution) [1] 6.073972 0.00 0.109124 0.00 0.00 0.00 0.863421 [8] 0.00 0.00 0.00 0.00 0.00 0.027455 0.00 Of course, since after the first command in that last block the second column of Amat contains only zeros Amat[,2] [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 we might as well have removed it (and the corresponding entry in bvec) Amat - Amat[, -2] bvec - bvec[-2] before calling solve.QP(). Note, the Boston data set was only used to illustrate how to fit such models, I do not want to imply that these models are sensible for these data. :-) Hope this helps. Cheers, Berwin Carlos Alzola [EMAIL PROTECTED] (703) 242-6747 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] XLSolutions 9 Courses: Upcoming March-April 2008 R/S+ Course Schedule by XLSolutions Corp
Our March-April 2008 R/S+ course schedule is now available. Please check out this link for additional information and direct enquiries to Sue Turner [EMAIL PROTECTED] Phone: 206 686 1578 --Can't see your city? Please email us! -- Ask for Group Discount --- www.xlsolutions-corp.com/courselist.htm (1) R/S System: Advanced Programming *** San Francisco / April 28-29, 2008 *** *** Salt Lake City / April 21-22, 2008 *** (2) R/S Fundamentals and Programming Techniques *** San Francisco / March 27-28, 2008 *** *** New York City / March 25-26, 2008 *** *** Seattle / April 21-22 (3) Data Mining: Practical Tools and Techniques in R/Splus *** San Franciso / April 24-25, 2008 *** (4) R/S System: Complementing and Extending Statistical Computing for SAS Users *** Raleigh / April 28-29, 2008 *** (5) Introduction to R/S+ programming: Microarrays Analysis and Bioconductor. *** Los Angeles / April 25-26, 2008*** (6) Microarrays Data Analysis with R/S+ and GGobi *** New York City / April 28-29, 2008*** www.xlsolutions-corp.com/courselist.htm Payment due AFTER the class Email us for group discounts. Email Sue Turner: [EMAIL PROTECTED] Phone: 206-686-1578 Visit us: www.xlsolutions-corp.com/courselist.htm Please let us know if you and your colleagues are interested in this class to take advantage of group discount. Register now to secure your seat! Cheers, Elvis Miller, PhD Manager Training. XLSolutions Corporation 206 686 1578 www.xlsolutions-corp.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] regex sulotion for seperating number and string
Hi Sun, vec - c(2324gz,2567 HK,3741,BF) vec1 - gsub('[^[:digit:]]','',vec) vec2 - gsub('[^[:alpha:]]','',vec) vec1 [1] 2324 2567 3741 vec2 [1] gz HK BF Cheers Vincenzo --- Vincenzo Luca Di Iorio Consultant PME User support - GSK RD Limited --- sun [EMAIL PROTECTED] Sent by: [EMAIL PROTECTED] 05-Mar-2008 15:51 To [EMAIL PROTECTED] cc Subject [R] regex sulotion for seperating number and string I have strings contain postcode and letters, some seperated with blank, some with comma, and some hasn't seperated. eg, 2324gz 2567 HK 3741,BF I want to seperate the number and letters into two new variables. I know this should be quite basic question, but searched on regex syntax and that seems a bit scarey to me, any one can shot me a quick solution on this particular question? thanks, Sun __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] CROSSOVER TRIALS IN R (Binary Outcomes)
I will like to analyse a binary cross over design using the random effects model. The probability of success is assumed to be logistic. Suppose as an example, we have 4 subjects undergoing a crossover design, where the outcome is either success or failure. The first two subjects receive treatment A first followed by treatment B. The remaining two subjects receive treatments in the reverse order. The outcomes for the subjects is the sequence AB are as follows: (0,1) and (0,0). While the outcomes for the subjects in the sequence BA are (1,1) and (1,0). How can i analyse this using R. I have done the problem with PROC NLMIXED in SAS, I simply want to compare the results from SAS with those from R. Please help. R-Codes will be highly appreciated. Boikanyo. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Asking, are simple effects different from 0
Huh. Very interesting. I haven't really worked with manipulating contrast matrices before, save to do a prior contrasts. Could you explain the matrix you laid out just a bit more so that I can generalize it to my case? Chuck Cleland wrote: One approach would be to use glht() in the multcomp package. You need to work out how to formulate the matrix of coefficients that give the desired contrasts. Here is an example using the warpbreaks data frame: fm - lm(breaks ~ tension*wool, data=warpbreaks) # names(coef(fm)) # (Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB cm - rbind( A vs. B at L = c(0, 0, 0,-1, 0, 0), A vs. B at M = c(0, 0, 0,-1,-1, 0), A vs. B at H = c(0, 0, 0,-1, 0,-1), M vs. L at A = c(0, 1, 0, 0, 0, 0), M vs. H at A = c(0, 1,-1, 0, 0, 0), L vs. H at A = c(0, 0,-1, 0, 0, 0), M vs. L at B = c(0, 1, 0, 0, 1, 0), M vs. H at B = c(0, 1,-1, 0, 1,-1), L vs. H at B = c(0, 0,-1, 0, 0,-1)) library(multcomp) summary(glht(fm, linfct = cm), test = adjusted(type=none)) Simultaneous Tests for General Linear Hypotheses Fit: lm(formula = breaks ~ tension * wool, data = warpbreaks) Linear Hypotheses: Estimate Std. Error t value p value A vs. B at L == 0 16. 5.1573 3.167 0.002677 ** A vs. B at M == 0 -4.7778 5.1573 -0.926 0.358867 A vs. B at H == 0 5.7778 5.1573 1.120 0.268156 M vs. L at A == 0 -20.5556 5.1573 -3.986 0.000228 *** M vs. H at A == 0 -0.5556 5.1573 -0.108 0.914665 L vs. H at A == 0 20. 5.1573 3.878 0.000320 *** M vs. L at B == 0 0.5556 5.1573 0.108 0.914665 M vs. H at B == 0 10. 5.1573 1.939 0.058392 . L vs. H at B == 0 9. 5.1573 1.831 0.073270 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- none method) -- View this message in context: http://www.nabble.com/Asking%2C-are-simple-effects-different-from-0-tp15835552p15852362.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] loop
Hello all, I am trying to use m - seq(-1,1,0.1) x1 - vector() x2 - vector() for(i in m){ x1[i] - i x2[i] - i^2 } dat - data.frame(x1,x2) But, I have false result dat x1 x2 1 1 1 could some tell me how it is possible to do this? Thank you! - [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] rrp.impute: for what sizes does it work?
Hi, I have a survey dataset of about 2 observations where for 2 factor variables I have about 200 missing values each. I want to impute these using 10 possibly explanatory variables which are a mixture of integers and factors. Since I was quite intrigued by the concept of rrp I wanted to use it but it takes ages and terminates with an error. First time it stopped complaining about too little memory which I increased then from 1.5 to 2GB. Now it terminates after a long time with this: Error in nn[[i]] : subscript out of bounds Has anybody encountered this problem before, is it due to my data? Could you recommend another package for such imputation? I looked at the aregImpute in Hmisc but I don't really understand what it is doing. Thanks a lot, Werner E-Mails jetzt auf Ihrem Handy. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Need help for calculating cross-correlation between 4 multivariate time series data
Hi all, I would like to know whether there is any function in R were i can find the cross-correlation of two or more multivariate (time series) data. I tried the function ccf() but it seems like to have two univariate datasets. Please let me know. sincerely, sandeep -- Sandeep Joseph PhD Post Doctoral Associate Center for Tropical Emerging Global Diseases Paul D. Coverdell Center, Rm 335A 500 D. W. Brooks Drive Athens, GA 30602-7394 ph 404-578-6790 2091 south Miledge Ave Apt A6, Athens, GA-30605. phone: 706-850-0056 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls: different results if applied to normal or linearized data
On Wed, 5 Mar 2008, Wolfgang Waser wrote: Dear all, I did a non-linear least square model fit y ~ a * x^b (a) nls(y ~ a * x^b, start=list(a=1,b=1)) to obtain the coefficients a b. I did the same with the linearized formula, including a linear model log(y) ~ log(a) + b * log(x) (b) nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1)) (c) lm(log10(y) ~ log10(x)) I expected coefficient b to be identical for all three cases. Hoever, using my dataset, coefficient b was: (a) 0.912 (b) 0.9794 (c) 0.9794 Coefficient a also varied between option (a) and (b), 107.2 and 94.7, respectively. Is this supposed to happen? Which is the correct coefficient b? Yes. You are fitting by least-squares on two different scales: differences in y and differences in log(y) are not comparable. Both are correct solutions to different problems. Since we have no idea what 'x' and 'y' are, we cannot even guess which is more appropriate in your context. Regards, Wolfgang -- Laboratory of Animal Physiology Department of Biology University of Turku FIN-20014 Turku Finland __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- 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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] loop
Date: Wed, 05 Mar 2008 15:59:59 +0100 (CET) From: Neuer Arkadasch [EMAIL PROTECTED] Sender: [EMAIL PROTECTED] Precedence: list Hello all, I am trying to use m - seq(-1,1,0.1) x1 - vector() x2 - vector() for(i in m){ x1[i] - i x2[i] - i^2 } dat - data.frame(x1,x2) But, I have false result dat x1 x2 1 1 1 It does not seem to be false: I am getting the same result. You should probably explain what you expected and what you are trying to achieve. Giovanni could some tell me how it is possible to do this? Thank you! - [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Giovanni Petris [EMAIL PROTECTED] Associate Professor Department of Mathematical Sciences University of Arkansas - Fayetteville, AR 72701 Ph: (479) 575-6324, 575-8630 (fax) http://definetti.uark.edu/~gpetris/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] loop
Try this: m - seq(-1,1,0.1) x1 - vector(length=length(m)) x2 - vector(length=length(m)) for(i in m){ x1[i] - i x2[i] - i^2 } dat - data.frame(x1,x2) Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Neuer Arkadasch Sent: Wednesday, March 05, 2008 10:00 AM To: [EMAIL PROTECTED] Subject: [R] loop Hello all, I am trying to use m - seq(-1,1,0.1) x1 - vector() x2 - vector() for(i in m){ x1[i] - i x2[i] - i^2 } dat - data.frame(x1,x2) But, I have false result dat x1 x2 1 1 1 could some tell me how it is possible to do this? Thank you! - [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] loop
m - seq(-1,1,0.1) x1 - c() x2 - c() for(i in 1:length(m)){ x1[i] - m[i] x2[i] - m[i]^2 } dat - data.frame(x1,x2) -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Neuer Arkadasch Sent: March 5, 2008 10:00 AM To: [EMAIL PROTECTED] Subject: [R] loop Hello all, I am trying to use m - seq(-1,1,0.1) x1 - vector() x2 - vector() for(i in m){ x1[i] - i x2[i] - i^2 } dat - data.frame(x1,x2) But, I have false result dat x1 x2 1 1 1 could some tell me how it is possible to do this? Thank you! - [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] regex sulotion for seperating number and string
Thanks all for the prompt answers!!! All works perfectly! up and running! Thanks! jim holtman [EMAIL PROTECTED] wrote in message news:[EMAIL PROTECTED] This should do it for you: x - c(2564gc, 2367,GH, 2134 JHG) x.sep - gsub(([[:digit:]]+)[ ,]*([[:alpha:]]+), \\1 \\2, x) # now create separate values strsplit(x.sep, ) [[1]] [1] 2564 gc [[2]] [1] 2367 GH [[3]] [1] 2134 JHG On 3/5/08, sun [EMAIL PROTECTED] wrote: I have strings contain postcode and letters, some seperated with blank, some with comma, and some hasn't seperated. eg, 2324gz 2567 HK 3741,BF I want to seperate the number and letters into two new variables. I know this should be quite basic question, but searched on regex syntax and that seems a bit scarey to me, any one can shot me a quick solution on this particular question? thanks, Sun __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem you are trying to solve? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] loop
Why not simply? m - seq(-1, 1, by = 0.1) dat - data.frame(m, m^2) - Erik Iverson Neuer Arkadasch wrote: Hello all, I am trying to use m - seq(-1,1,0.1) x1 - vector() x2 - vector() for(i in m){ x1[i] - i x2[i] - i^2 } dat - data.frame(x1,x2) But, I have false result dat x1 x2 1 1 1 could some tell me how it is possible to do this? Thank you! - [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R_alloc with structures with flexible array members
Ramon Diaz-Uriarte wrote on 03/05/2008 04:25 AM: Dear All, In a package, I want to use some C code where I am using a structure (as the basic element of a linked list) with flexible array members. Basically, this is a structure where the last component is an incomplete array type (e.g., Harbison Steel, C, a reference manual, 5th ed., p. 159) such as: struct Sequence { struct Sequence *next; int len; unsigned int state_count[]; }; To create one such sequence, I allocate storage (following Harbison and Steel) in a C program as follows: struct Sequence *A; int n = 4; A = malloc( sizeof(struct Sequence) + n * sizeof(unsigned int)); If I understand correctly, however, it would be better to use R_alloc instead of malloc (memory automagically freed on exit and error; error-checking). But I do not know how to make the call to R_alloc here, since R_alloc allocates n units of size bytes each. I've tried, without success, the following two: int to_add_for_R_alloc = (int) ceil((float) sizeof(struct sequence) / sizeof(unsigned int)); A = (struct sequence *) R_alloc(to_add_for_R_alloc + n, sizeof(unsigned int)); or even a brute force attempt as: A = (struct sequence *) R_alloc( 100, sizeof(struct sequence)); but both result in segmentation faults. Should I just keep using malloc (and free at end)? Hi Ramon, You should be able to use R_alloc without seg faults, so there's something wrong with your code somewhere. R_alloc multiplies its arguments together to come up with the total number of bytes to allocate then it allocates a raw vector and returns the data portion. So you can just treat R_alloc similarly to malloc by calling R_alloc(1,sizeof(struct Sequence) + n * sizeof(unsigned int)). Best, Jeff -- http://biostat.mc.vanderbilt.edu/JeffreyHorner __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] rrp.impute: for what sizes does it work?
rrp is working! Sorry, it was my mistake... fiddling around to find out what the problem is I forgot to re-include the variables which are to be imputed. It seems like this case is not caught but the algorithm finishes with the mentioned error. Anyway, I am still a little fuzzy about imputation and browsing the web I found no clear recommendations regarding this topic. Which R function would you recommend for imputing categorical variables / factors? Many thanks, Werner Hi, I have a survey dataset of about 2 observations where for 2 factor variables I have about 200 missing values each. I want to impute these using 10 possibly explanatory variables which are a mixture of integers and factors. Since I was quite intrigued by the concept of rrp I wanted to use it but it takes ages and terminates with an error. First time it stopped complaining about too little memory which I increased then from 1.5 to 2GB. Now it terminates after a long time with this: Error in nn[[i]] : subscript out of bounds Has anybody encountered this problem before, is it due to my data? Could you recommend another package for such imputation? I looked at the aregImpute in Hmisc but I don't really understand what it is doing. Thanks a lot, Werner Lesen Sie Ihre E-Mails jetzt einfach von unterwegs. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] ipf function in R
On Wed, 5 Mar 2008, Chandra Shah wrote: Hi I have a 3 x 2 contingency table: 10 20 30 40 50 60 I want to update the frequencies to new marginal totals: 100 130 40 80 110 I want to use the ipf (iterative proportional fitting) function which is apparently in the cat package. Can somebody please advice me how to input this data and invoke ipf in R to obtain an updated contingency table? I'd use loglin() newtab - loglin( rowmarg%o%colmarg/sum(colmarg), margin=list(1,2), start=tab, fit=TRUE )$fit with rowmarg and colmarg set to your updated marginals. As for inputting the data, if this is all you have, type it in at the command line. See ?matrix ?c and note this: PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. HTH, Chuck Thanks. By the way I am quite new to R. -- Dr Chandra Shah Senior Research Fellow Monash University-ACER Centre for the Economics of Education and Training Faculty of Education, Building 6, Monash University Victoria Australia 3800 Tel. +61 3 9905 2787 Fax +61 3 9905 9184 www.education.monash.edu.au/centres/ceet __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:[EMAIL PROTECTED] UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] package for repeated measures ANOVA with various link functions REDUX
On Tue, Mar 4, 2008 at 9:48 PM, John Sorkin [EMAIL PROTECTED] wrote: Prof. Bates was correct to point out the lack of specifics in my original posting. I am looking for a package that will allow we to choose among link functions and account for repeated measures in a repeated measures ANOVA. My question is what package should I use to facilitate estimating rates of illegal drug use at three centers, and the effect two interventions have on usage. At each center data describing the rate of drug use was obtained once a month. For the first six-months, there was no intervention at any of the three centers. For months seven through 13 intervention one was applied at each of the three centers. For months 14 through 24 intervention two was applied at each center. The question I am trying to answer is did intervention one or two change drug usage at any of the three centers. I am treating center as a repeated measure, i.e. the rate of drug use at month one will be correlated with the rate of drug use at center one at months two, three, etc. I have accounted for repeated measures several ways in the past. (1) I have used SAS proc MIXED with a REPEATED statement. The REPEATED statement allows for the specification of the within-subject correlation of repeated measures by specifying the structure of the within-subject variance-covariance matrix of the repeated measures. The matrix is block diagonal with one block for each subject. But does such a structure extend to models with binary or count responses? You have mentioned that you want to use an arbitrary link function such as quasibinomial. What I understand the effect of the REPEATED statement to be is to specify a parameterized form of the marginal variance-covariance matrix of the responses. If the response variable has a multivariate normal distribution it is possible to independently specify the mean (determined by the fixed-effects parameters) and the marginal variance-covariance. However, in the case of generalized linear models the mean response is determined by a linear predictor and a link function while the variance-covariance of the response is determined by prior weights and a variance function. The same is true for generalized linear mixed models except that this description applies to the conditional distribution of the response given the random effects. The link and the variance functions must agree so, for example, using a logit or probit link which restricts the value of mu to the interval [0,1] would imply a variance function (up to prior weights) of mu(1-mu). At least I think so - others may feel that it is possible to specify an arbitrary variance function but I don't see how that can make sense. To me the whole point of generalized linear models is to transform the linear predictor to the desired range for the mean and to take into account what this implies about the variance. Even if you feel that it is possible to relax the ties between the link function and the variance function I don't see how it would be possible to specify an arbitrary structure for the marginal variance-covariance of the response. If you say that the marginal variance-covariance must have a block-wise compound symmetry structure but you are going to restrict the mean to the range [0,1] I think you have painted yourself into a corner. I don't think it is possible to specify a mean on a restricted range and separately specify an arbitrary variance-covariance structure. In particular, when the mean is on the range [0,1] then you better have the variance going to zero as the mean goes to 0 or to 1. You can't arbitrarily say that the variance within a block must be constant, regardless of the values of the means in those blocks. (2) I have used SAS proc GENMOD which uses GEE to adjust the parameter estimates and their standard errors for the fact that a repeated measurements of a parameter are obtained from a given subjects. Is there any package in R that will allow me to perform a repeated measures ANOVA with a selection of link functions that will allow me to account for repeated measures by either specifying the correlation structure of the repeated measures from a subject a la SAS proc mixed or by adjusting the parameter estimates using GEE a la proc GENMOD? Perhaps R has a package that accounts for repeated measures in some other manner. Thank you, John Sorkin John Sorkin M.D., Ph.D. Chief, Biostatistics and Informatics University of Maryland School of Medicine Division of Gerontology Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call phone number above prior to faxing) Douglas Bates [EMAIL PROTECTED] 3/4/2008 5:13 PM On Tue, Mar 4, 2008 at 10:52 AM, John Sorkin [EMAIL PROTECTED] wrote: R 2.6.0 Windows XP At the risk of raising the ire of the R gods . .
[R] Correlation matrix one side with significance
Hi there! In my case, cor(d[1:20]) makes me a good correlation matrix. Now I'd like to have it one sided, means only the left bottom side to be printed (the others are the same) and I'd like to have * where the p-value is lower than 0.05 and ** lower than 0.01. How can I do this? And another thing: Is there a way to output that table as a latex table? Thanks, Martin -- Ihr Partner für Webdesign, Webapplikationen und Webspace. http://www.roomandspace.com/ Martin Kaffanke +43 650 4514224 signature.asc Description: Dies ist ein digital signierter Nachrichtenteil __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Asking, are simple effects different from 0
On 3/5/2008 10:09 AM, jebyrnes wrote: Huh. Very interesting. I haven't really worked with manipulating contrast matrices before, save to do a prior contrasts. Could you explain the matrix you laid out just a bit more so that I can generalize it to my case? Each column corresponds to one of the coefficients in the model, and each row specifies a particular contrast. The numbers in the matrix indicate how the model coefficients are combined to indicate a particular difference in means. For example, the first row indicates that the third coefficient (woolB) is multiplied by -1. The baseline categories are A and L for the wool and tension factors, so the woolB effect in fm is the simple effect of B vs. A in the baseline category of the tension factor. Multiplying this coefficient by -1 produces an A vs. B comparison in the baseline category of the tension factor. When I want to check that contrasts are as intended, I use contrast() in the contrast package (by Steve Weston, Jed Wing, Max Kuhn). That function allows you to specify factor levels by name to construct the contrast. For example: library(contrast) # M vs. H at B contrast(fm, a=list(tension = M, wool = B), b=list(tension = H, wool = B)) lm model parameter contrast Contrast S.E. LowerUppert df Pr(|t|) 10 5.157299 -0.3694453 20.36945 1.94 48 0.0584 It also allows you to print the design matrix for a contrast: contrast(fm, a=list(tension = M, wool = B), b=list(tension = H, wool = B))$X (Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB 1 01 -1 0 1 -1 Chuck Cleland wrote: One approach would be to use glht() in the multcomp package. You need to work out how to formulate the matrix of coefficients that give the desired contrasts. Here is an example using the warpbreaks data frame: fm - lm(breaks ~ tension*wool, data=warpbreaks) # names(coef(fm)) # (Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB cm - rbind( A vs. B at L = c(0, 0, 0,-1, 0, 0), A vs. B at M = c(0, 0, 0,-1,-1, 0), A vs. B at H = c(0, 0, 0,-1, 0,-1), M vs. L at A = c(0, 1, 0, 0, 0, 0), M vs. H at A = c(0, 1,-1, 0, 0, 0), L vs. H at A = c(0, 0,-1, 0, 0, 0), M vs. L at B = c(0, 1, 0, 0, 1, 0), M vs. H at B = c(0, 1,-1, 0, 1,-1), L vs. H at B = c(0, 0,-1, 0, 0,-1)) library(multcomp) summary(glht(fm, linfct = cm), test = adjusted(type=none)) Simultaneous Tests for General Linear Hypotheses Fit: lm(formula = breaks ~ tension * wool, data = warpbreaks) Linear Hypotheses: Estimate Std. Error t value p value A vs. B at L == 0 16. 5.1573 3.167 0.002677 ** A vs. B at M == 0 -4.7778 5.1573 -0.926 0.358867 A vs. B at H == 0 5.7778 5.1573 1.120 0.268156 M vs. L at A == 0 -20.5556 5.1573 -3.986 0.000228 *** M vs. H at A == 0 -0.5556 5.1573 -0.108 0.914665 L vs. H at A == 0 20. 5.1573 3.878 0.000320 *** M vs. L at B == 0 0.5556 5.1573 0.108 0.914665 M vs. H at B == 0 10. 5.1573 1.939 0.058392 . L vs. H at B == 0 9. 5.1573 1.831 0.073270 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- none method) -- Chuck Cleland, Ph.D. NDRI, Inc. 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] inheritence in S4
Well well well... To summarize : let assume that A is a class (slot x) and C is a class containing A (A and slot y) - as(c,A) calls new(A). So new(A) HAS TO works, you can not decide to forbid empty object (unless you define setAs(C,A) ?) - In addition, any test that you would like to set in initialize or validity have to first check is some field are empty (because 'if([EMAIL PROTECTED] 0)' will fail if x=numerical(0)) - When you call new(C), the neither new(A) nor intialize(A) are called. () So, all the security test you write in initialize A, you have to rewrite them on C ? I start to undestand why not that much people use S4... Christophe [EMAIL PROTECTED] writes: Hi Martin, thanks for your answer But a couple of other quick points. I would have written setMethod(initialize, A, function(.Object, ..., xValue=numeric(0)){ callNextMethod(.Object, ..., x=xValue) }) I am not that much familiar with S3... In our way of writing this method, 'initialize' for 'A' will call the next method for A, which is initialise' for 'numeric', is that right ? from below, setClass(A, representation(x=numeric)) The 'next' method refers to arguments in the signature of the generic, i.e., the 'next' method relevant to .Object. For A, the next method is the default initialize method, which will use all named arguments to fill the slots of .Object (and perhaps call validObject, among other things). You can see the code with showMethods(initialize, classes=ANY, includeDef=TRUE) S3 is not relevant here. And finally, the position of 'xValue' and 'yValue' means that the arugment has to be named, e.g., new(B, yValue=12). This seems a little awkward at first, but seems like a best practice I agree with you. But I do not like the use of ... , it lets us to make many mistake like in : print(3.5165,digitts=2) There is a typo in digitts, but since print works with ... , R does not detect the mistake. So I agree with the importance of naming argument, I always do it but without 'forcing' it If this were initialize, and you had provided an invalid named argument, the default method would have failed because there is no slot of that name. And finally, in Jim's thread I mention using a constructor. So in practice for a case like the above I would not define any initialize methods, Interesting. Will you still define a validity method or not even ? For a simple case like this there is no extra validity to check beyond that ensured by S4 (e.g., slots of the correct type). If there were constraints, e.g., only positive values, or length one vectors, then I would define a validity function. Martin Christophe Christophe Genolini [EMAIL PROTECTED] writes: Thanks Martin Well it works except that as seems to not like the initialize method : the following code (that is the same than yours with some initialize for A B and C) does not compile. It seems that as(c,A) does not work if we definie a initialize for A... --- 8 -- setClass(A, representation(x=numeric)) setMethod(initialize,A,function(.Object,value)[EMAIL PROTECTED] - value;return(.Object)}) a - new(A,4) setClass(B, representation(y=numeric)) setMethod(initialize,B,function(.Object,value)[EMAIL PROTECTED] - value;return(.Object)}) b - new(B,5) setClass(C, contains=c(A, B)) setMethod(initialize,C,function(.Object,valueA, valueB){ [EMAIL PROTECTED] - valueA [EMAIL PROTECTED] - valueB return(.Object) }) c - new(C,valueA=10,valueB=12) setMethod(show, A, function(object) cat(A\n)) setMethod(show, B, function(object) cat(B\n)) setMethod(show, C, function(object) { callGeneric(as(object, A)) callGeneric(as(object, B)) cat(C\n) }) c --- 8 Is there something wrong with the use of 'as' between class and father class? Christophe Hi Christophe -- I don't know whether there's a particularly elegant way. This works setClass(A, representation(x=numeric)) setClass(B, representation(y=numeric)) setClass(C, contains=c(A, B)) setMethod(show, A, function(object) cat(A\n)) setMethod(show, B, function(object) cat(B\n)) setMethod(show, C, function(object) { callGeneric(as(object, A)) callGeneric(as(object, B)) cat(C\n) }) new(C) A B C but obviously involves the developer in making explicit decisions about method dispatch when there is multiple inheritance. Martin [EMAIL PROTECTED] writes: Hi the list I define a class A (slot a and b), a class C (slot c and d) and a class E that inherit from A and B. I define print(A) and print(B). For print(C), I would like to use both of them, but I do not see how... Thanks for your help... Christophe Ce message a ete envoye par IMP, grace a l'Universite Paris 10 Nanterre __ R-help@r-project.org mailing list
Re: [R] Reversed but positive axis in trellis plots?
On 3/5/08, Fredrik Karlsson [EMAIL PROTECTED] wrote: Hi, In my discpipline, it is common to plot one acoustic property on a positive scale but from top to bottom on the ordinate and the same for another measurement on the abscissa. So, the origin of the plot is on the top right of the plot, with increasing values to the left /down. This is to highlight the correlation between the acoustic measurement and the position of the forming structure, for instance when teaching it to students. The grouping ability of the trellis plot is quite handy whan plotting many instances of the same thing, so I was wondering if it is possible to make a trellis xyplot behave this way? Converting all values to negative and changing labels to the negative of the negative seems one solution to the reverseness of the axes, but how do I change the position? Is it possible? You can always specify explicit limits that are reversed, e.g., ylim=c(100, 0). If you want this as a general behaviour but don't know your data's range beforehand, one option might have been prepanel = function(x, y, ...) { list(ylim = rev(range(y))) } Unfortunately, this doesn't work when relation=same, as the step of combining the per-panel limits to obtain a common range disregards the order. -Deepayan __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] loop
Thank you Yinghai, that's what I need :-)! Yinghai Deng [EMAIL PROTECTED] schrieb: m - seq(-1,1,0.1) x1 - c() x2 - c() for(i in 1:length(m)){ x1[i] - m[i] x2[i] - m[i]^2 } dat - data.frame(x1,x2) -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Neuer Arkadasch Sent: March 5, 2008 10:00 AM To: [EMAIL PROTECTED] Subject: [R] loop Hello all, I am trying to use m - seq(-1,1,0.1) x1 - vector() x2 - vector() for(i in m){ x1[i] - i x2[i] - i^2 } dat - data.frame(x1,x2) But, I have false result dat x1 x2 1 1 1 could some tell me how it is possible to do this? Thank you! - [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. - Lesen Sie Ihre E-Mails jetzt einfach von unterwegs.. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] density plot
Hi, I'm trying to create a density plot which I used to do in geneplotter using the following code. Unfortunately I can't find the combination of R release and geneplotter that works. Can anyone suggest a fix or an alternative to smoothScatter that will plot depth of one dive vs depth of the next as density? Many thanks, Tom library(geneplotter) require(RColorBrewer) layout(matrix(1:4, ncol=2, byrow=TRUE)) op - par(mar=rep(2,4)) x - cbind(depthcurrent1,depthcurrent2) smoothScatter(x,nrpoints=0, colramp=colorRampPalette(c(white, black)), ylab=Depth of Next Dive, xlab=Depth of Current Dive) colors - densCols(x) plot(x, col=black, pch=20, ylab=Depth of Next Dive, xlab=Depth of Current Dive) The Zoological Society of London is incorporated by Royal Charter Principal Office England. Company Number RC000749 Registered address: Regent's Park, London, England NW1 4RY Registered Charity in England and Wales no. 208728 _ This e-mail has been sent in confidence to the named add...{{dropped:20}} __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] CROSSOVER TRIALS IN R (Binary Outcomes)
On Wed, 5 Mar 2008, Boikanyo Makubate wrote: I will like to analyse a binary cross over design using the random effects model. The probability of success is assumed to be logistic. Suppose as an example, we have 4 subjects undergoing a crossover design, where the outcome is either success or failure. The first two subjects receive treatment A first followed by treatment B. The remaining two subjects receive treatments in the reverse order. The outcomes for the subjects is the sequence AB are as follows: (0,1) and (0,0). While the outcomes for the subjects in the sequence BA are (1,1) and (1,0). How can i analyse this using R. I have done the problem with PROC NLMIXED in SAS, I simply want to compare the results from SAS with those from R. Please help. R-Codes will be highly appreciated. You need to be clear about the statistical model. In crossover trials, there is usually a subject effect that one conditions out of the likelihood (often implicitly). In this case, clogit() in the survival package would be suitable for such an approach. Something like res - clogit( outcomes ~ rx = strata( subject ) ) However, the example you give is degenerate in the same way that an ordinary logistic regression is when the responses are linearly separable in the regressor space. So, a well crafted program will tell you that it cannot find an answer (and clogit does this) for your example data. If you knew all of this and had some other statistical model in mind (such as one that models the subject effects rather than conditioning on them), you should repost telling us what model you wanted to fit. HTH, Chuck Boikanyo. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:[EMAIL PROTECTED] UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Correlation matrix one side with significance
Try this: On 05/03/2008, Martin Kaffanke [EMAIL PROTECTED] wrote: Hi there! In my case, cor(d[1:20]) makes me a good correlation matrix. Now I'd like to have it one sided, means only the left bottom side to be printed (the others are the same) and I'd like to have * where the p-value is lower than 0.05 and ** lower than 0.01. How can I do this? d - matrix(rexp(16, 2), 4) corr - cor(d) sign - symnum(cor(d), cutpoints=c(0.05, 0.01), corr = T, symbols=c(***, **, *), abbr=T, diag=F) noquote(mapply(function(x, y)paste(x, format(y, dig=3), sep=''), as.data.frame(unclass(sign)), as.data.frame(corr))) And another thing: Is there a way to output that table as a latex table? See ?latex function in Hmisc package and also xtable package Thanks, Martin -- Ihr Partner für Webdesign, Webapplikationen und Webspace. http://www.roomandspace.com/ Martin Kaffanke +43 650 4514224 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Henrique Dallazuanna Curitiba-Paraná-Brasil 25° 25' 40 S 49° 16' 22 O __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] vector manipulations
Thank you everybody. Phil, your expand.grid works very nicely and I will use it for non-vectorized functions. Yet I am a bit confused about vectorization. For me it is synonymous of no loop. :-( I wrote a toy example (with a function which is not my log-likelihood). FIRST PART nir=1:10 logl=function(x,y,nir) sum(log(x*nir+y)) x=seq(0.1,0.3,by=10^(-1)) y=seq(0.1,0.3,by=10^(-1)) z=outer(x,y,logl,nir=nir) This does not work. Can you explain me why it is not vectorised ? SECOND PART nir=1:10 logl2=function(x,y,nir) { a=0 for (i in 1:10) { a=a+log(x*nir[i]+y) } return(a) } x=seq(0.1,0.3,by=10^(-1)) y=seq(0.1,0.3,by=10^(-1)) z2=outer(x,y,logl2,nir=nir) This seems to work, though the function does not seem to be vectorized. I am sorry for being such a noob. I'm ok in maths but I am bad at programming. I bought a book on R (Introductory Statistics with R by Dalgaard) ** on Amazon last week . I will read it when I receive it. Do you know other good books ? 2008/3/5, [EMAIL PROTECTED] [EMAIL PROTECTED]: No problems with it working. The main problem I have observed is unrealistic expectations. People write an *essentially* non-vectorized function and expect Vectorize to produce a version of it which will out-perform explicit loops every time. No magic bullets in this game. Bill. Bill Venables CSIRO Laboratories PO Box 120, Cleveland, 4163 AUSTRALIA Office Phone (email preferred): +61 7 3826 7251 Fax (if absolutely necessary): +61 7 3826 7304 Mobile: +61 4 8819 4402 Home Phone: +61 7 3286 7700 mailto:[EMAIL PROTECTED] http://www.cmis.csiro.au/bill.venables/ -Original Message- From: Duncan Murdoch [mailto:[EMAIL PROTECTED] Sent: Wednesday, 5 March 2008 9:36 AM To: Venables, Bill (CMIS, Cleveland) Cc: r-help@r-project.org Subject: Re: [R] vector manipulations On 3/4/2008 5:41 PM, [EMAIL PROTECTED] wrote: Your problem is that your function log1( , ) is not vectorized with respect to its arguments. For a function to work in outer(...) it must accept vectors for its first two arguments and it must produce a parallel vector of responses. To quote the help information for outer: FUN is called with these two extended vectors as arguments. Therefore, it must be a vectorized function (or the name of one), expecting at least two arguments. Sometimes Vectorize can be used to make a non-vectorized function into a vectorized one, but the results are not always entirely satisfactory in my experience. What problems have you seen? Duncan Murdoch __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R_alloc with structures with flexible array members
Dear Jeff, Thanks for the suggestion. However, something is still not working. This is a simple example: *** start C #include R.h struct Sequence { int len; unsigned int state_count[]; }; int main(void) { struct Sequence *A; int n = 4; // First line segfaults. Second doesn't A = (struct Sequence *) R_alloc(1, sizeof(struct Sequence) + n * sizeof(unsigned int)); // A = malloc(sizeof(struct Sequence) + n * sizeof(unsigned int)); return(0); } *** end C ** I then do gcc -std=gnu99 -Wall -I/usr/share/R/include -I/usr/share/R/include -L/usr/lib/R/lib -lR ex7.c and the ./a.out segfaults when I use R_alloc (not with malloc). Best, R. On Wed, Mar 5, 2008 at 5:23 PM, Jeffrey Horner [EMAIL PROTECTED] wrote: Ramon Diaz-Uriarte wrote on 03/05/2008 04:25 AM: Dear All, In a package, I want to use some C code where I am using a structure (as the basic element of a linked list) with flexible array members. Basically, this is a structure where the last component is an incomplete array type (e.g., Harbison Steel, C, a reference manual, 5th ed., p. 159) such as: struct Sequence { struct Sequence *next; int len; unsigned int state_count[]; }; To create one such sequence, I allocate storage (following Harbison and Steel) in a C program as follows: struct Sequence *A; int n = 4; A = malloc( sizeof(struct Sequence) + n * sizeof(unsigned int)); If I understand correctly, however, it would be better to use R_alloc instead of malloc (memory automagically freed on exit and error; error-checking). But I do not know how to make the call to R_alloc here, since R_alloc allocates n units of size bytes each. I've tried, without success, the following two: int to_add_for_R_alloc = (int) ceil((float) sizeof(struct sequence) / sizeof(unsigned int)); A = (struct sequence *) R_alloc(to_add_for_R_alloc + n, sizeof(unsigned int)); or even a brute force attempt as: A = (struct sequence *) R_alloc( 100, sizeof(struct sequence)); but both result in segmentation faults. Should I just keep using malloc (and free at end)? Hi Ramon, You should be able to use R_alloc without seg faults, so there's something wrong with your code somewhere. R_alloc multiplies its arguments together to come up with the total number of bytes to allocate then it allocates a raw vector and returns the data portion. So you can just treat R_alloc similarly to malloc by calling R_alloc(1,sizeof(struct Sequence) + n * sizeof(unsigned int)). Best, Jeff -- http://biostat.mc.vanderbilt.edu/JeffreyHorner -- Ramon Diaz-Uriarte Statistical Computing Team Structural Biology and Biocomputing Programme Spanish National Cancer Centre (CNIO) http://ligarto.org/rdiaz __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Replace values in data.frame conditional on another data.frame
Dear List, I am looking for an efficient method for replacing values in a data.frame conditional on the values of a separate data.frame. Here is my scenario: I have a data.frame (A) with say 1000 columns, and 365 rows. Each cell in the data.frame has either valid value, or NA. I have an additional data.frame (B) with the same number of rows and columns, with valid values in all cells. What I would like to do, is replace the cells in B with NA if and only if the corresponding cell in A is NA. I have search extensively for a method to do this, and I'm sure that in the end the solution will be embarrassingly simple, however, any help is greatly appreciated! Cheers, Carson Apologies if this has been posted twice, we are currently experiencing server problems... __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Replace values in data.frame conditional on another data.frame
Try bb[is.na(aa)] - NA It may be simple but it is not necessarily obvious :) --- Carson Farmer [EMAIL PROTECTED] wrote: Dear List, I am looking for an efficient method for replacing values in a data.frame conditional on the values of a separate data.frame. Here is my scenario: I have a data.frame (A) with say 1000 columns, and 365 rows. Each cell in the data.frame has either valid value, or NA. I have an additional data.frame (B) with the same number of rows and columns, with valid values in all cells. What I would like to do, is replace the cells in B with NA if and only if the corresponding cell in A is NA. I have search extensively for a method to do this, and I'm sure that in the end the solution will be embarrassingly simple, however, any help is greatly appreciated! Cheers, Carson Apologies if this has been posted twice, we are currently experiencing server problems... __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] New data source - now how do we build an R interface?
Folks, A nice new data resource has come up -- http://data.un.org/ I thought it would be wonderful to setup an R function like tseries::get.hist.quote() which would be able to pull in some or all of this data. I walked around a bit of it and I'm not able to map the resources to predictable URLs which can then be wget. There's some javascript going on that I'm not understanding. Perhaps someone on the mailing list can think about this? -- Ajay Shah http://www.mayin.org/ajayshah [EMAIL PROTECTED] http://ajayshahblog.blogspot.com *(:-? - wizard who doesn't know the answer. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] loop
m - seq(-1,1,0.1) x1 - vector() x2 - vector() # the loop statement was incorrect. for(i in 1:length(m)){ x1[i] - m[i] x2[i] - m[i]^2 } dat - data.frame(x1,x2) # But why not something like this? There is no need for a loop. x1 - seq(-1,1,0.1) mdat - data.frame(x1, x2=x1^2) --- Neuer Arkadasch [EMAIL PROTECTED] wrote: Hello all, I am trying to use m - seq(-1,1,0.1) x1 - vector() x2 - vector() for(i in m){ x1[i] - i x2[i] - i^2 } dat - data.frame(x1,x2) But, I have false result dat x1 x2 1 1 1 could some tell me how it is possible to do this? Thank you! - [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Asking, are simple effects different from 0
Ah. I see. So, if I want to test to see whether each simple effect is different from 0, I would do something like the following: cm2 - rbind( A:L = c(1, 0, 0, 0, 0, 0), A:M = c(1, 1, 0, 0, 0, 0), A:H = c(1, 0, 1, 0, 0, 0), B:L = c(1, 0, 0, 1, 0, 0), B:M = c(1, 1, 0, 1, 1, 0), B:H = c(1, 0, 1, 1, 0, 1)) summary(glht(fm, linfct = cm2), test = adjusted(type=none)) Correct? What is the df on those t-tests then? Is it 48? Interestingly, I find this produces results no different than fm2-lm(breaks ~ tension:wool+0, data=warpbreaks) summary(fm2) Also, here, it would seem each t-test was done with the full 48df. Hrm. Chuck Cleland wrote: Each column corresponds to one of the coefficients in the model, and each row specifies a particular contrast. The numbers in the matrix indicate how the model coefficients are combined to indicate a particular difference in means. For example, the first row indicates that the third coefficient (woolB) is multiplied by -1. The baseline categories are A and L for the wool and tension factors, so the woolB effect in fm is the simple effect of B vs. A in the baseline category of the tension factor. Multiplying this coefficient by -1 produces an A vs. B comparison in the baseline category of the tension factor. -- View this message in context: http://www.nabble.com/Asking%2C-are-simple-effects-different-from-0-tp15835552p15855630.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] box-constrained
Hello everybody, I have a question about box-constrained optimization. I've done some research and I found that optim could do that. Are there other ways in R ? Is the following correct if I have a function f of two parameters belonging for example to [0,1] and [0,Infinity] ? optim(par=param, fn=f, method=L-BFGS-B, lower=c(0,0), upper=c(1,Inf)) My other question is whether it is possible to add the derivatives of my function (like in nlm) and whether it is better to add them ? If there is no need to add the derivatives, then I guess I could wish to optimize the likelihood directly rather than to optimize the log-likelihood... Indeed one aspect of the log-likelihood is to make the derivatives tractable (in iid cases). Do you agree ? Thank you ! Gustave [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] legend for several graphics
?mtitle should do it. --- Georg Otto [EMAIL PROTECTED] wrote: Hi, I am trying to generate a figure of 9 plots that are contained in one device by using par(mfrow = c(3,3,)) I would like to have 1 common legend for all 9 plots somewhere outside of the plotting area (as opposed to one legend inside each of the 9 plots, which the function legend() seems to generate by default). Any hint how to do this? Best, Georg __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] simulation study using R
Thanks to All, The comments were very helpful; however, the the simulation is running very slow. I reduced the number of loops (conditions) so I have 36 loops, and the data-generation occurs 1000 times within each loop. At the end of each 1000 reps, I saved the summary (e.g., mean) of the reps to a single row vector, and then saved it in a fileshare database. When the simulation is finished, I rbind() the 36 rows of the database objects into a final simulation result. Thanks again, Davood Tofighi On Wed, Mar 5, 2008 at 8:44 AM, Paul Gilbert [EMAIL PROTECTED] wrote: Davood Tofighi wrote: Thanks for your reply. For each condition, I will have a matrix or data frames of 1000 rows and 4 columns. I also have a total of 64 conditions for now. So, in total, I will have 64 matrices or data frames of 1000 rows and 4 columns. The format of data I would like to store would be data frames or matrices. I also would like to store the data for later use, I generally find it is better to store the seed and other data you need to reproduce the experiment, rather than trying to store the data (see, for example, package setRNG). If you save only the summary statistics you need, then you can usually do it in memory. (Be sure to assign variables for the statistics to their full size first and then populate them, rather than extending them at each step.) If you write things to files then it will slow down your simulation a lot. In fact, in most cases it will be quicker to re-run the experiment than it will be to read the data from disk. Paul e.g., a plot of the empirical distribution of the chi^2, or to compute the power of Chi^2 across 1000 reps for each condition. On Mon, Mar 3, 2008 at 7:03 PM, jim holtman [EMAIL PROTECTED] wrote: What is the format of the data you are storing (single value, multivalued vector, matrix, dataframe, ...)? This will help formulate a solution. What do you plan to do with the data? Are you going to do further analysis, write it to flat files, store it in a data base, etc.? How big are the data objects you are manipulating? On Mon, Mar 3, 2008 at 7:05 PM, Davood Tofighi [EMAIL PROTECTED] wrote: Dear All, I am running a Monte Carlo simulation study and have some questions on how to manage data storage efficiently at the end of each 1000 replication loop. I have three conditions coded using the FOR {} loops and a FOR loop that generates data for each condition, performs analysis, and computes a statistic 1000 times. Therefore, for each condition, I will have 1000 statistic values. My question is what's the best way to store the 1000 statistic for each condition. Any suggestion on how to manage such simulation studies is greatly appreciated. Thanks, -- Davood Tofighi Department of Psychology Arizona State University [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem you are trying to solve? La version française suit le texte anglais. This email may contain privileged and/or confidential information, and the Bank of Canada does not waive any related rights. Any distribution, use, or copying of this email or the information it contains by other than the intended recipient is unauthorized. If you received this email in error please delete it immediately from your system and notify the sender promptly by email that you have done so. Le présent courriel peut contenir de l'information privilégiée ou confidentielle. La Banque du Canada ne renonce pas aux droits qui s'y rapportent. Toute diffusion, utilisation ou copie de ce courriel ou des renseignements qu'il contient par une personne autre que le ou les destinataires désignés est interdite. Si vous recevez ce courriel par erreur, veuillez le supprimer immédiatement et envoyer sans délai à l'expéditeur un message électronique pour l'aviser que vous avez éliminé de votre ordinateur toute copie du courriel reçu. -- Davood Tofighi Department of Psychology Arizona State University P.O. BOX 871104 Tempe, AZ 85287-1104 Tel.:480-727-7884 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and
Re: [R] Asking, are simple effects different from 0
On 3/5/2008 1:32 PM, jebyrnes wrote: Ah. I see. So, if I want to test to see whether each simple effect is different from 0, I would do something like the following: cm2 - rbind( A:L = c(1, 0, 0, 0, 0, 0), A:M = c(1, 1, 0, 0, 0, 0), A:H = c(1, 0, 1, 0, 0, 0), B:L = c(1, 0, 0, 1, 0, 0), B:M = c(1, 1, 0, 1, 1, 0), B:H = c(1, 0, 1, 1, 0, 1)) That does not corresponds to what I think of as the simple effects. That specifies the six cell means, but it does not *compare* any cell means. I think of a simple effect as the effect of one factor at a specific level of some other factor. summary(glht(fm, linfct = cm2), test = adjusted(type=none)) Correct? What is the df on those t-tests then? Is it 48? Yes, df = 48 for each contrast. Interestingly, I find this produces results no different than fm2-lm(breaks ~ tension:wool+0, data=warpbreaks) summary(fm2) Yes, but those are not what I would call the simple effects. Those are essentially one-sample t-tests for each of the 6 cell means. Also, here, it would seem each t-test was done with the full 48df. Hrm. The df are based on the whole model, not the 9 observations in one cell. Chuck Cleland wrote: Each column corresponds to one of the coefficients in the model, and each row specifies a particular contrast. The numbers in the matrix indicate how the model coefficients are combined to indicate a particular difference in means. For example, the first row indicates that the third coefficient (woolB) is multiplied by -1. The baseline categories are A and L for the wool and tension factors, so the woolB effect in fm is the simple effect of B vs. A in the baseline category of the tension factor. Multiplying this coefficient by -1 produces an A vs. B comparison in the baseline category of the tension factor. -- Chuck Cleland, Ph.D. NDRI, Inc. 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] N-way ANOVA
Dear All, Can R perform n-way ANOVA, i.e., with 3 or more factors? Thanks in advance, Paul __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] legend for several graphics
On Wed, 2008-03-05 at 15:28 +0100, Georg Otto wrote: Hi, I am trying to generate a figure of 9 plots that are contained in one device by using par(mfrow = c(3,3,)) I would like to have 1 common legend for all 9 plots somewhere outside of the plotting area (as opposed to one legend inside each of the 9 plots, which the function legend() seems to generate by default). Any hint how to do this? Here's one way: op - par(mfrow = c(3,3), ## split region oma = c(5,0,4,0) + 0.1, ## create outer margin mar = c(5,4,2,2) + 0.1) ## shrink some margins plot(1:10, main = a, pch = 1:2, col= 1:2) plot(1:10, main = b, pch = 1:2, col= 1:2) plot(1:10, main = c, pch = 1:2, col= 1:2) plot(1:10, main = d, pch = 1:2, col= 1:2) plot(1:10, main = e, pch = 1:2, col= 1:2) plot(1:10, main = f, pch = 1:2, col= 1:2) plot(1:10, main = g, pch = 1:2, col= 1:2) plot(1:10, main = h, pch = 1:2, col= 1:2) plot(1:10, main = i, pch = 1:2, col= 1:2) ## title mtext(My Plots, side = 3, outer = TRUE, font = 2, line = 1, cex = 1.2) ## draw legend legend(-12.5, -6, legend = c(Type 1, Type 2), pch = 1:2, col = 1:2, ncol = 2) par(op) I had to fiddle by hand with the legend x and y locations to get it roughly centred. There has to be better way - probably something to do with reseting the plot region, but I can't recall how to do that now. If there is, I'm sure someone will tell me what I overlooked. Is this what you were looking for? G Best, Georg __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Dr. Gavin Simpson [t] +44 (0)20 7679 0522 ECRC, UCL Geography, [f] +44 (0)20 7679 0565 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/ UK. WC1E 6BT. [w] http://www.freshwaters.org.uk %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] box-constrained
Hi, Let me make the following points in response to your questions: 1. Your call to optim() with L-BFGS-B as the method is correct. Just make sure that your function f is defined as negative log-likelihood, since optim is by default a minimizer. The other option is to define log-likelihood as usual, but set control=list(fnscale=-1). 2. You can add derivative (or gradient to be more precise) by defining that function and then using the gr argument in optim. Specifying exact gradient almost always improves the convergence of the iterative schemes, especially for ill-conditioned problems (flat region around the local minima). So, if it is not too much trouble, and you are confident of your differentiation skills, you should do that. However, in most cases the approximate finite-difference gradient used by optim() should be good enough. 3. Regardless of whether it is easy to compute the exact gradient or not, it is generally a bad idea to maximize the likelihood that involves the product of a large number of very small numbers. It is almost always better to maximize the log-likelihood. Since the objective function is additive rather than multiplicative, it has better numerical conditioning. Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Gustave Lefou Sent: Wednesday, March 05, 2008 1:34 PM To: r-help@r-project.org Subject: [R] box-constrained Hello everybody, I have a question about box-constrained optimization. I've done some research and I found that optim could do that. Are there other ways in R ? Is the following correct if I have a function f of two parameters belonging for example to [0,1] and [0,Infinity] ? optim(par=param, fn=f, method=L-BFGS-B, lower=c(0,0), upper=c(1,Inf)) My other question is whether it is possible to add the derivatives of my function (like in nlm) and whether it is better to add them ? If there is no need to add the derivatives, then I guess I could wish to optimize the likelihood directly rather than to optimize the log-likelihood... Indeed one aspect of the log-likelihood is to make the derivatives tractable (in iid cases). Do you agree ? Thank you ! Gustave [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] N-way ANOVA
On Wed, 5 Mar 2008, Paul Smith wrote: Dear All, Can R perform n-way ANOVA, i.e., with 3 or more factors? Yes. There are even examples on the help page! Thanks in advance, Paul __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- 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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] inheritence in S4
[EMAIL PROTECTED] writes: Well well well... You're partly misunderstanding... To summarize : let assume that A is a class (slot x) and C is a class containing A (A and slot y) - as(c,A) calls new(A). So new(A) HAS TO works, you can not decide to forbid empty object (unless you define setAs(C,A) ?) new(A) has to work (return a valid object), but you can still create arbitrary objects with a prototype (provided the prototype is consistent with your validity function) and / or with initialize methods (which do not have any named arguments, other than .Object, without default values). Creating valid objects seems like a good idea! - In addition, any test that you would like to set in initialize or validity have to first check is some field are empty (because if([EMAIL PROTECTED] 0)' will fail if x=numerical(0)) The validity and initialize methods have to be written in a way consistent with valid objects -- if you want your object to contain zero-length vectors, then the validity test has to accommodate that (e.g., if (length([EMAIL PROTECTED])==0 || all([EMAIL PROTECTED]0)) {}). If your initialize method is supposed to return a valid object, then you'd better construct one! An initialize method on A that expected an argument xValue whose default was 0 might be written as setMethod(A, initialize, function(.Object, ..., xValue=0) { callNextMethod(.Object, ..., x=xValue) }) For this simple case you could create the object with a suitable prototype and not have an initialize method at all setClass(A, representation(x=numeric), prototype=prototype(x=0)) - When you call new(C), the neither new(A) nor intialize(A) are called. () new(A) is not called (why would it be?). If you have an initialize method on A then it will be called. The problem you experienced before was that your initialize method for A REQUIRED additional arguments. So, all the security test you write in initialize A, you have to rewrite them on C ? That is not correct; you do not have to duplicate code to construct an object of class C. I start to undestand why not that much people use S4... Christophe [EMAIL PROTECTED] writes: Hi Martin, thanks for your answer But a couple of other quick points. I would have written setMethod(initialize, A, function(.Object, ..., xValue=numeric(0)){ callNextMethod(.Object, ..., x=xValue) }) I am not that much familiar with S3... In our way of writing this method, 'initialize' for 'A' will call the next method for A, which is initialise' for 'numeric', is that right ? from below, setClass(A, representation(x=numeric)) The 'next' method refers to arguments in the signature of the generic, i.e., the 'next' method relevant to .Object. For A, the next method is the default initialize method, which will use all named arguments to fill the slots of .Object (and perhaps call validObject, among other things). You can see the code with showMethods(initialize, classes=ANY, includeDef=TRUE) S3 is not relevant here. And finally, the position of 'xValue' and 'yValue' means that the arugment has to be named, e.g., new(B, yValue=12). This seems a little awkward at first, but seems like a best practice I agree with you. But I do not like the use of ... , it lets us to make many mistake like in : print(3.5165,digitts=2) There is a typo in digitts, but since print works with ... , R does not detect the mistake. So I agree with the importance of naming argument, I always do it but without 'forcing' it If this were initialize, and you had provided an invalid named argument, the default method would have failed because there is no slot of that name. And finally, in Jim's thread I mention using a constructor. So in practice for a case like the above I would not define any initialize methods, Interesting. Will you still define a validity method or not even ? For a simple case like this there is no extra validity to check beyond that ensured by S4 (e.g., slots of the correct type). If there were constraints, e.g., only positive values, or length one vectors, then I would define a validity function. Martin Christophe Christophe Genolini [EMAIL PROTECTED] writes: Thanks Martin Well it works except that as seems to not like the initialize method : the following code (that is the same than yours with some initialize for A B and C) does not compile. It seems that as(c,A) does not work if we definie a initialize for A... --- 8 -- setClass(A, representation(x=numeric)) setMethod(initialize,A,function(.Object,value)[EMAIL PROTECTED] - value;return(.Object)}) a - new(A,4) setClass(B, representation(y=numeric)) setMethod(initialize,B,function(.Object,value)[EMAIL PROTECTED] - value;return(.Object)}) b - new(B,5) setClass(C, contains=c(A, B)) setMethod(initialize,C,function(.Object,valueA, valueB){ [EMAIL PROTECTED] - valueA [EMAIL PROTECTED] - valueB
Re: [R] R_alloc with structures with flexible array members
On Wed, 5 Mar 2008, Ramon Diaz-Uriarte wrote: Dear Jeff, Thanks for the suggestion. However, something is still not working. This is a simple example: *** start C #include R.h struct Sequence { int len; unsigned int state_count[]; }; int main(void) { struct Sequence *A; int n = 4; // First line segfaults. Second doesn't A = (struct Sequence *) R_alloc(1, sizeof(struct Sequence) + n * sizeof(unsigned int)); // A = malloc(sizeof(struct Sequence) + n * sizeof(unsigned int)); return(0); } *** end C ** I then do gcc -std=gnu99 -Wall -I/usr/share/R/include -I/usr/share/R/include -L/usr/lib/R/lib -lR ex7.c and the ./a.out segfaults when I use R_alloc (not with malloc). You can't use R_alloc in a standalone program without initializing R, which has not been done here. You said 'in a package', but this is not in a package. Best, R. On Wed, Mar 5, 2008 at 5:23 PM, Jeffrey Horner [EMAIL PROTECTED] wrote: Ramon Diaz-Uriarte wrote on 03/05/2008 04:25 AM: Dear All, In a package, I want to use some C code where I am using a structure (as the basic element of a linked list) with flexible array members. Basically, this is a structure where the last component is an incomplete array type (e.g., Harbison Steel, C, a reference manual, 5th ed., p. 159) such as: struct Sequence { struct Sequence *next; int len; unsigned int state_count[]; }; To create one such sequence, I allocate storage (following Harbison and Steel) in a C program as follows: struct Sequence *A; int n = 4; A = malloc( sizeof(struct Sequence) + n * sizeof(unsigned int)); If I understand correctly, however, it would be better to use R_alloc instead of malloc (memory automagically freed on exit and error; error-checking). But I do not know how to make the call to R_alloc here, since R_alloc allocates n units of size bytes each. I've tried, without success, the following two: int to_add_for_R_alloc = (int) ceil((float) sizeof(struct sequence) / sizeof(unsigned int)); A = (struct sequence *) R_alloc(to_add_for_R_alloc + n, sizeof(unsigned int)); or even a brute force attempt as: A = (struct sequence *) R_alloc( 100, sizeof(struct sequence)); but both result in segmentation faults. Should I just keep using malloc (and free at end)? Hi Ramon, You should be able to use R_alloc without seg faults, so there's something wrong with your code somewhere. R_alloc multiplies its arguments together to come up with the total number of bytes to allocate then it allocates a raw vector and returns the data portion. So you can just treat R_alloc similarly to malloc by calling R_alloc(1,sizeof(struct Sequence) + n * sizeof(unsigned int)). Best, Jeff -- http://biostat.mc.vanderbilt.edu/JeffreyHorner -- Ramon Diaz-Uriarte Statistical Computing Team Structural Biology and Biocomputing Programme Spanish National Cancer Centre (CNIO) http://ligarto.org/rdiaz __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- 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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] types of vectors / lists
Hello, I am an advanced user of R. Recently I found out that apparently I do not fully understand vectors and lists fully Take this code snippet: T = c(02.03.2008 12:23, 03.03.2008 05:54) Times = strptime(T, %d.%m.%Y %H:%M) Times # OK class(Times) # OK is.list(Times)# sort of understand and not understand that length(Times) # 9 ??? why is it length(Times[1]) ?? Times[1] # OK Times[[1]] # Wrong so Times is a vector-style thing of a non-atomic type. What puzzles me is first that is.list returns true, and more importantly that the length-query returns the length of the first object of that apparent list and not how many elements are contained (i.e., I would have expected 2 as result - and I do not know what syntax to use in order to get the 2). Moreover, if the whole thing is part of a data.frame: DFTimes = as.data.frame(Times) dim(DFTimes) length(DFTimes$Times) # OK, 2 then everything is as expected. Could anyone please clearify why is.list returns true and why length in the first example returns 9 ? Is it that c() makes a list if the objects to be concatenated are not representable directly by atomic types ? thanks, Thomas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls: different results if applied to normal or linearized data
On 6/03/2008, at 2:53 AM, Wolfgang Waser wrote: Dear all, I did a non-linear least square model fit y ~ a * x^b (a) nls(y ~ a * x^b, start=list(a=1,b=1)) to obtain the coefficients a b. I did the same with the linearized formula, including a linear model log(y) ~ log(a) + b * log(x) (b) nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1)) (c) lm(log10(y) ~ log10(x)) I expected coefficient b to be identical for all three cases. Hoever, using my dataset, coefficient b was: (a) 0.912 (b) 0.9794 (c) 0.9794 Coefficient a also varied between option (a) and (b), 107.2 and 94.7, respectively. Is this supposed to happen? Which is the correct coefficient b? The two approaches assume two different models. Model (1) is y = a*x^b + E (where different errors are independent and identically --- usually normally --- distributed). Model (2) is y = a*(x^b)*E (and you are usually tacitly assuming that ln E is normally distributed). The point estimates of a and b will consequently be different --- although usually not hugely different. Their distributional properties will be substantially different. cheers, Rolf Turner ## Attention:\ This e-mail message is privileged and confid...{{dropped:9}} __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] collapsing list to data.frame
Hello, Given a list with all elements having identical layout, e.g.: l = NULL l[[1]] = list(4, hello) l[[2]] = list(7, world) l[[3]] = list(9, ) is there an easy way to collapse this list into a data.frame with each row being the elements of the list ? I.e. in this case I want to convert the list into a data.frame with 3 rows and 2 columns, where column 1 holds the integer values, and column 2 the character values. I can get it done by looping over all elements and rbind them together to the final result, but that is quite slow (for large sets) and ugly, so I was wondering if there's an easy syntax. thanks a lot in advance, Thomas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] N-way ANOVA
On Wed, Mar 5, 2008 at 7:41 PM, Prof Brian Ripley [EMAIL PROTECTED] wrote: Can R perform n-way ANOVA, i.e., with 3 or more factors? Yes. There are even examples on the help page! Thanks, Prof. Ripley. I will have a look at it. Paul __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] differentiating a numeric vector
What functions exist for differentiating a numeric vector (in my case spectral data)? That is, experimental data without an analytical function. ie, x - seq(1,10,0.1) y=x^3+rnorm(length(x),sd=0.01) #although the real function would be nothing simple like x^3... derivy - I know I could just use diff(y) but it would be nice to estimate derivatives at the endpoints, calculate higher-order derivatives, and maybe have some smoothing options ie by differentiating a spline or something like that. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] legend for several graphics
On Wed, Mar 5, 2008 at 8:28 AM, Georg Otto [EMAIL PROTECTED] wrote: Hi, I am trying to generate a figure of 9 plots that are contained in one device by using par(mfrow = c(3,3,)) I would like to have 1 common legend for all 9 plots somewhere outside of the plotting area (as opposed to one legend inside each of the 9 plots, which the function legend() seems to generate by default). If you provide more detail about your problem (what are the 9 plots?) I'm sure we can suggest other solutions using lattice or ggplot that will substantially simplify your code, as well as automatically drawing the legend. Hadley -- http://had.co.nz/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Asking, are simple effects different from 0
Indeed, but are not each of the cell means also evaluations of the effect of one factor at the specific level of another factor? Is this an issue of Tomato, tomahto. I guess my question is, if I want to know if each of those is different from 0, then should I use the 48df from the full model, or the 9 for each cell? Chuck Cleland wrote: That does not corresponds to what I think of as the simple effects. That specifies the six cell means, but it does not *compare* any cell means. I think of a simple effect as the effect of one factor at a specific level of some other factor. summary(glht(fm, linfct = cm2), test = adjusted(type=none)) Correct? What is the df on those t-tests then? Is it 48? Yes, df = 48 for each contrast. Interestingly, I find this produces results no different than fm2-lm(breaks ~ tension:wool+0, data=warpbreaks) summary(fm2) Yes, but those are not what I would call the simple effects. Those are essentially one-sample t-tests for each of the 6 cell means. Also, here, it would seem each t-test was done with the full 48df. Hrm. The df are based on the whole model, not the 9 observations in one cell. -- View this message in context: http://www.nabble.com/Asking%2C-are-simple-effects-different-from-0-tp15835552p15857771.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] collapsing list to data.frame
Perhaps data.frame(do.call(rbind, l)) ? - Erik Iverson [EMAIL PROTECTED] wrote: Hello, Given a list with all elements having identical layout, e.g.: l = NULL l[[1]] = list(4, hello) l[[2]] = list(7, world) l[[3]] = list(9, ) is there an easy way to collapse this list into a data.frame with each row being the elements of the list ? I.e. in this case I want to convert the list into a data.frame with 3 rows and 2 columns, where column 1 holds the integer values, and column 2 the character values. I can get it done by looping over all elements and rbind them together to the final result, but that is quite slow (for large sets) and ugly, so I was wondering if there's an easy syntax. thanks a lot in advance, Thomas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] read.table
Dear R community, I encounter a problem reading data into a dataframe. See attachment for the input. I use: data-read.table(test,fill=T,row.names=1) When you look at data you can see that some lines of my input were broken... I can avoid this problem by sorting the lines by length... Do I have it wrong or is there a bug? Thanks and wishing you a great remainder of the day, Georg. ** Georg Ehret JHU Baltimore - US __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] collapsing list to data.frame
try this: l - vector(list, 3) l[[1]] - list(4, hello) l[[2]] - list(7, world) l[[3]] - list(9, ) lis - lapply(l, names-, value = c(V1, V2)) do.call(rbind, lapply(lis, data.frame, stringsAsFactors = FALSE)) I hope it helps. Best, Dimitris Dimitris Rizopoulos Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm Quoting [EMAIL PROTECTED]: Hello, Given a list with all elements having identical layout, e.g.: l = NULL l[[1]] = list(4, hello) l[[2]] = list(7, world) l[[3]] = list(9, ) is there an easy way to collapse this list into a data.frame with each row being the elements of the list ? I.e. in this case I want to convert the list into a data.frame with 3 rows and 2 columns, where column 1 holds the integer values, and column 2 the character values. I can get it done by looping over all elements and rbind them together to the final result, but that is quite slow (for large sets) and ugly, so I was wondering if there's an easy syntax. thanks a lot in advance, Thomas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Axes in polymars
Hi All, I can't quite figure out how to change the parameters of the x and y axes when I plot a polymars object. I want to add a scatterplot of the data points, but the polymars plot seems to automatically set the parameters to fit the polymars line. I've tried using plot(poly,1,fig= c(x1,x2,y1,y2)) but have no luck. Any thoughts? Thanks, Ryan Please open the following attachment for important information regarding this e-mail communication. This email communication is privileged, confidential or otherwise protected by disclosure and is intended only for the individuals or entities named above and any others who have been specifically authorized to receive it. Any unauthorized dissemination, copying or use of the contents of this email is strictly prohibited and may be in violation of law. If you are not the intended recipient, please do not read, copy, use or disclose to others the contents of this communication. Please notify the sender that you have received this e-mail in error by replying to this e-mail or by phoning +1-212-583-5000 Monday to Friday between 7:30 am and 7:30 pm (EST). At any other time, please call +1-212-583-5221. Please then delete the e-mail from your system and any copies of it. Nothing contained in this disclaimer shall be construed in any way to grant permission to transmit confidential information via this firm's e-mail system or as a waiver of any confidentiality or privilege. This communication may contain highly confidential information regarding Blackstone's investments, strategy and organization. Your acceptance of this communication from Blackstone constitutes your agreement to (i) keep confidential all the confidential information contained in this communication, as well as any information derived by you from the confidential information contained in this communication (collectively, Confidential Information) and not disclose any such Confidential Information to any other person, (ii) not use any of the Confidential Information for any purpose other than pursuant to your relationship with Blackstone, (iii) not use the Confidential Information for purposes of trading any security, including, without limitation, securities of Blackstone or its portfolio companies, (iv) not copy any Confidential Information without Blackstone's prior consent and (v) promptly return any Confidential Information contained in this communication to Blackstone upon Blackstone's request. This communication is for informational purposes only and does not constitute an offer to sell or a solicitation of an offer to purchase any interest in any investment vehicles managed by any business unit of The Blackstone Group. Blackstone does not accept any responsibility or liability arising from the use of this communication. No representation is being made that the information presented is accurate, current or complete, and such information is at all times subject to change without notice. Opinions expressed may differ or be contrary to the opinions and recommendations of a Blackstone business unit. Blackstone does not provide legal, accounting or tax advice. Any statement regarding legal, accounting or tax matters was written in connection with the explanation of the matters described herein and was not intended or written to be relied upon by any person as definitive advice. Any discussion of U.S. tax matters contained within this communication is not intended to be used and cannot be used for the purpose of avoiding penalties that may be imposed under applicable Federal, state or local tax law or recommending to another party any transaction or matter addressed herein. Each person should seek advice based on its particular circumstances from independent legal, accounting and tax advisors regarding the matters discussed in this e-mail.__ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] types of vectors / lists
Hello - [EMAIL PROTECTED] wrote: Hello, I am an advanced user of R. Recently I found out that apparently I do not fully understand vectors and lists fully Take this code snippet: T = c(02.03.2008 12:23, 03.03.2008 05:54) Times = strptime(T, %d.%m.%Y %H:%M) Times # OK class(Times) # OK is.list(Times)# sort of understand and not understand that ?POSIXlt says Class 'POSIXlt' is a named list of vectors representing 'sec' 0-61: seconds 'min' 0-59: minutes 'hour' 0-23: hours 'mday' 1-31: day of the month 'mon' 0-11: months after the first of the year. 'year' Years since 1900. 'wday' 0-6 day of the week, starting on Sunday. 'yday' 0-365: day of the year. 'isdst' Daylight savings time flag. Positive if in force, zero if not, negative if unknown. length(Times) # 9 ??? why is it length(Times[1]) ?? Because Times is a list with 9 elements, try typing names(Times) to see what those elements are called. That will help you understand. Times[1] # OK Times[[1]] # Wrong Whether or not this is 'wrong' has been I believe debated before on the list. See this message and the rest of the thread associated with it. http://tolstoy.newcastle.edu.au/R/help/06/07/31508.html - Erik Iverson __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.