Re: [R] date handling
d = as.POSIXlt(c(2005-07-01, 2005-07-02, 2005-07-03, 2005-07-04, 2005-07-05)) d$mon and d$year will get you part way there. That is assuming your dates are formated -mm-dd. strptime() might also be useful. Richard van Wingerden wrote: Hi, Given a frame with calendar date's: 2005-07-01, 2005-07-02,2005-07-03,2005-07-04,2005-07-05,etc. I want to extract the following from these dates: week number month number year number Any ideas how to accomplish this? Many thanks. Regards, Richard __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Bivariate Splines in R
Cal == Cal Stats [EMAIL PROTECTED] on Mon, 12 Dec 2005 10:25:38 -0800 (PST) writes: Cal Hi.., is there a function in R to fit bivariate splines Cal ? I came across 'polymars' (POLSPLINE) and 'mars' Cal (mda) packages. Are these the one to use or are there Cal other specific commands? I'd recommend to use gam(y ~ s(x1,x2)) from the recommended package 'mgcv'. help(gam) has many examples, some of which using bivariate splines. Cal Thanks. Cal Harsh you're welcome; Martin Maechler, ETH Zurich __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] date handling
D = as.POSIXlt(c(2005-07-01, 2005-07-02, 2005-07-03, 2005-07-04, 2005-07-05)) (Years = 1900+D$year) [1] 2005 2005 2005 2005 2005 (Months = D$mon+1) [1] 7 7 7 7 7 weekdays(D) [1] Friday Saturday Sunday Monday Tuesday # you better check this one carefully !!! (Weeknumbers = floor((D$yday+7)/7)) JeeBee On Mon, 12 Dec 2005 12:59:11 +0100, Richard van Wingerden wrote: Hi, Given a frame with calendar date's: 2005-07-01, 2005-07-02,2005-07-03,2005-07-04,2005-07-05,etc. I want to extract the following from these dates: week number month number year number Any ideas how to accomplish this? Many thanks. Regards, Richard __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Technique for reading large sparse fwf data file
Dear list: A datafile was sent to me that is very large (92890 x 1620) and is *very* sparse. Instead of leaving the entries with missing data blank, each cell with missing data contains a dot (.) The data are binary in almost all columns, with only a few columns containing whole numbers, which I believe requires 2 bytes for the binary and 4 for the others. So, by my calculations (assuming 4 bytes for all cells to create an upperbound) I should need around 92890 * 1620 * 4 = 574MB to read in these data and about twice that for analyses. My computer has 3GB. But, I am unable to read in the file even though I have allocated sufficient memory to R for this. My first question is do the dots in the empty cells consume additional memory? I am assuming the answer is yes and believe I should remove them before I do the read in. Because my data are in a fixed width format file, I can open the file in a text editor and find and replace all dots with nothing. Then, I should retry the read in process? Maybe this will work? I created a smaller data file (~ 14000 * 1620) in SAS and tried to import this subset (it still had the dots), but R still would not allow for me to do so. I could use a little guidance as I think I have allocated sufficient memory to read in a datafile assuming my calculations are right. Does anyone have any thoughts on a strategy? Harold [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Labeling a range of bars in barplot?
Hi, I am plotting a distribution of (ordered) values as a barplot. I would like to label groups of bars together to highlight aspects of the distribution. The label for the group should be the range of values in those bars. As this is hard to describe, here is an example; x - rlnorm(50)*2 barplot(sort(x,decreasing=T)) y - quantile(x, seq(0, 1, 0.2)) y plot(diff(y)) That last plot is to highlight that I want to label lots of the small columns together, and have a few more labels for the bigger columns (more densely labeled). I guess I will have to turn out my own labels using low level plotting functions, but I am stumped as to how to perform the calculation for label placement. I imagine drawing several line segments, one for each group of bars to be labeled together, and putting the range under each line segment as the label. Each line segment will sit under the group of bars that it covers. Thanks for any help with the above! Cheers, Dan. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Labeling a range of bars in barplot?
Is this what you want... op-par(xpd=TRUE) segments(0,-1,20,-1,col=2) text(10,-1,Interval,pos=1) Cheers, Hans Gardfjell Dept. of Ecology and Environmental science Umeå University, Sweden Dan Bolser wrote: Hi, I am plotting a distribution of (ordered) values as a barplot. I would like to label groups of bars together to highlight aspects of the distribution. The label for the group should be the range of values in those bars. As this is hard to describe, here is an example; x - rlnorm(50)*2 barplot(sort(x,decreasing=T)) y - quantile(x, seq(0, 1, 0.2)) y plot(diff(y)) That last plot is to highlight that I want to label lots of the small columns together, and have a few more labels for the bigger columns (more densely labeled). I guess I will have to turn out my own labels using low level plotting functions, but I am stumped as to how to perform the calculation for label placement. I imagine drawing several line segments, one for each group of bars to be labeled together, and putting the range under each line segment as the label. Each line segment will sit under the group of bars that it covers. Thanks for any help with the above! Cheers, Dan. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] about empirical sample size in partial correlations
Hello everyone I have estimated the partial correlations of a matrix (n=160, vars=7) with command pcor.shrink (library corpcor) and then i want to estimate their respective confidence intervals. So i have tried to use the command pcor.confint(pcor,kappa,alpha) from library GeneNT. My problem is that i don't know how to estimate kappa (empirical sample size). Is there any command with wich i can estimate kappa or is there any other command for estimating confidence intervals for partial correlations without using the kappa parameter? Thank you _ http://www.mailbox.gr ÁðïêôÞóôå äùñåÜí ôï ìïíáäéêü óáò e-mail. http://www.superweb.gr ÏéêïíïìéêÜ êáé áîéüðéóôá ðáêÝôá web hosting ìå áóöáëÝò Åëëçíéêü controlpanel http://www.domains.gr Ôï üíïìÜ óáò óôï internet ìüíï ìå 10 Åõñþ. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Question
Hello, I have a problem that I am trying to solve and I am not sure how to do it in R. Suppose, that 16 numbers are choosen at random from 0 to 9, what's the probability that their average will be between 4 and 6. I typed the following code: set.seed(100) sample(0:9, 16, replace =TRUE) [1] 3 2 5 0 4 4 8 3 5 1 6 8 2 3 7 6 Is what I got, however I realize the set.seed function locks in the number I get every time. My question is in order to run a true random sample, wouldn't I have to use the runif function? And then deliminate the sample to show the numbers that lie between 4 and 6? If that's the case, how do I do that? Davia S. Cox 517-575-8031 cell [EMAIL PROTECTED] Human potential, though not always apparent, is there waiting to be discovered and invited forth. -William W. Purkey __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Question
Please disregard this message and don't post it to the web. I found the answer. Thanks Davia S. Cox 517-575-8031 cell [EMAIL PROTECTED] Human potential, though not always apparent, is there waiting to be discovered and invited forth. -William W. Purkey On Dec 13, 2005, at 6:20 AM, Davia Cox wrote: Hello, I have a problem that I am trying to solve and I am not sure how to do it in R. Suppose, that 16 numbers are choosen at random from 0 to 9, what's the probability that their average will be between 4 and 6. I typed the following code: set.seed(100) sample(0:9, 16, replace =TRUE) [1] 3 2 5 0 4 4 8 3 5 1 6 8 2 3 7 6 Is what I got, however I realize the set.seed function locks in the number I get every time. My question is in order to run a true random sample, wouldn't I have to use the runif function? And then deliminate the sample to show the numbers that lie between 4 and 6? If that's the case, how do I do that? Davia S. Cox 517-575-8031 cell [EMAIL PROTECTED] Human potential, though not always apparent, is there waiting to be discovered and invited forth. -William W. Purkey __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Problem with understanding output of Cox model
Hi all, I am using a 'tricked' Cox Hazard regression model for discrete choice habitat modelling. However, I'm having a hard time understanding the meaning of the first line the following part of the summary() output: Rsquare= 0.307 (max possible= 0.475 ) Likelihood ratio test= 91.8 on 12 df, p=2.23e-14 Wald test = 26.3 on 12 df, p=0.00977 Score (logrank) test = 58.6 on 12 df, p=4.03e-08 Does anyone know how I can read the 'max possible' R-square? What does it signify? Is it some kind of quasi-correlation; in which case I read it that I have explained 65% (i.e., 0.307/0.475) of the variation? I hope someone can help on this. Thanks in advance, Cheers Roel Roel May Norwegian Wolverine Project Norwegian Institute for Nature Research (NINA) Tungasletta 2, N-7485 Trondheim, Norway Tlf. +47 73 80 14 65, Mob. +47 95 78 59 95 Email [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] Internett www.nina.no http://www.nina.no/ , www.jerv.info http://www.jerv.info/ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Time delay function or plot animation
Tobias, thanks to you and others, who made the same suggestion. This is exactly what I was looking for, Sys.sleep(0.1) does the job much more quietly than my previous solution ... Regards, Ulrike -- Original Message --- From: Tobias Verbeke [EMAIL PROTECTED] To: Ulrike Grömping [EMAIL PROTECTED] Cc: r-help@stat.math.ethz.ch Sent: Mon, 12 Dec 2005 21:14:26 +0100 Subject: Re: [R] Time delay function or plot animation Ulrike Grömping wrote: is it possible to specify a time delay for plotting the points in a curve? I would like to make the plotting process slow enough to show the development of the graph, and therefore I am looking either for the possibility within the plot function to specify a plotting speed or (if that doesn't exist) for a function like pause or wait that allows to specify a time delay until the next statement is executed. I have searched help and mailing list archives, but I don't seem to look for the right keywords. My current solution is very crude: I add the following to an existing plot: for (i in 1:steps) { lines(x[i:(i+1)], y[i:(i+1)]) #calculations with the sole purpose of generating a time delay zoeger-1 for (j in 1:85000) {zoeger-zoeger*j} } Is there a better way to achieve this? See ?Sys.sleep HTH, Tobias Regards, Ulrike __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html --- End of Original Message --- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Comments please, how to set your mailer to read R-help in digest format
2005/12/12, Michael Dewey [EMAIL PROTECTED]: There are occasional comments on this list about how difficult it is to read the digest format. [snip] Any other comments welcome of course. I use google mail with a R-Label. In this way, all my R-Mails are available in the R-Label-view and are searchable. The conversations are instanteous updated and threadwise ordered. So it feels quasi the same as a real newsgroup (which I would prefer btw). [Well I think, it's a major space waste if everybody stores just the same mails on his/her email account but since Google offers 1 GB space, who cares]. Best regards, Hans-Peter [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] help for multivariate analysis
dear R users, I need some help for multivariate analysis. I have 2 anaesthetic treatment groups (20 patients/group) where I register heart frequency and pressure for 60 min (repeated measures every 5 minutes). I would like to perform a test to check if treatments are different in controlling freq and pressures during the anaesthesia, but i would like to have also an overall measure and not only multiple p for different time intervals. I also think I should choose a test in which time is meaningful since the measures are not simple repeated measurements but measurements taken at specific time points. 1 million dollar question how to do in R? thanks in advance Dr Matteo Vidali Dep. of Medical Sciences University of East Piedmont A. Avogadro ITALY -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] help with writing function
I'm trying to write a function that takes a vector of length n and then takes the first value of the vector i.e j=1 and forms a new vector of length n (i.e replicate the first value n times). This function will then calculate the absoulte difference of the original vector and the new vector and store the results omitting the difference between the value and itself. This function should be able to repeat the procedure for each of the j's i.e j=2 to n. The results should all be stored together. Below is what I've tried so far but it seems to work only for j=1 . Your help will be highly appreciated. IED-function(risk){ n-length(risk) i-c(1:n) Diff-numeric() for(j in 1:n){ relrisk-risk relrisk[i]-relrisk[j] Difference-abs(risk-relrisk) Difference-Difference[-c(1:j)] Difference-append(Diff,Difference) return(Difference) } } Oarabile __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Question
On 12/13/05 6:20 AM, Davia Cox [EMAIL PROTECTED] wrote: Hello, I have a problem that I am trying to solve and I am not sure how to do it in R. Suppose, that 16 numbers are choosen at random from 0 to 9, what's the probability that their average will be between 4 and 6. I typed the following code: set.seed(100) sample(0:9, 16, replace =TRUE) [1] 3 2 5 0 4 4 8 3 5 1 6 8 2 3 7 6 Is what I got, however I realize the set.seed function locks in the number I get every time. Just don't use set.seed() before every run (unless you want to always get the same answers). Set.seed() is available to allow you to generate reproducible results, so not using it means that you will get a different set of random numbers every time you run your sample from above. Sean __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Question
Davia Cox wrote: Hello, I have a problem that I am trying to solve and I am not sure how to do it in R. Suppose, that 16 numbers are choosen at random from 0 to 9, what's the probability that their average will be between 4 and 6. I typed the following code: set.seed(100) sample(0:9, 16, replace =TRUE) [1] 3 2 5 0 4 4 8 3 5 1 6 8 2 3 7 6 Is what I got, however I realize the set.seed function locks in the number I get every time. Yes, that's what set.seed is intended to do... otherwise don't use set.seed (and make your work unreproducible). My question is in order to run a true random sample, We have to disappoint you: Your computer cannot generate true random samples. wouldn't I have to use the runif function? And then deliminate the sample to show the No, if you want integers, the sample above fits perfectly well. numbers that lie between 4 and 6? If that's the case, how do I do that? Is this a homework question? Uwe Ligges Davia S. Cox 517-575-8031 cell [EMAIL PROTECTED] Human potential, though not always apparent, is there waiting to be discovered and invited forth. -William W. Purkey __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] sample matrix
Please, I´d like to store this sample matrix as a new object. How can I do this ? pulse - c(67, 67, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 71, 71, 72, 72, 73, 74) m - NULL x - 0 for (i in 1:5) { x - sample(pulse,3) m - mean(x) cat(x,m,\n) } Thanks, Mauricio __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Generation of missiing values in a time serie...
Hi, First, thank you. It 'almost' works: I can convert a matrix to a 'zoo' object, but it can not convert it to a 'ts'. The sitruation is this: I load a set of samples with h_types - list (0, 0, character(0), character(0), 0, 0, 0, 0, 0) h_names - list (time, flow, type, dir, seq, ts, x, rtt, size) pcks_file - pipe (grep ' P ' server.dat), r) pcks- scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) Read 1429 records mat - data.frame (pcks) colnames (mat) - h_names mat time flow type dir seq ts x rtt size 1 1.0008930P +0 1.000893 1472 0.00 1472 2 1.5144540P +1 1.514454 2944 0.513142 1472 3 2.0150930P +2 2.015093 2944 0.513142 1472 4 2.5150250P +3 2.515025 4806 0.504488 1472 5 2.8219760P +4 2.821976 5730 0.496728 1472 6 3.0789310P +5 3.078931 5832 0.489744 1472 7 3.3318970P +6 3.331897 5832 0.489744 1472 [...] 1425 176.9255040P + 1424 176.925504 12141 0.764699 1472 1426 177.0394890P + 1425 177.039489 12141 0.764699 1472 1427 177.1534690P + 1426 177.153469 12141 0.764699 1472 1428 177.2674640P + 1427 177.267464 12141 0.764699 1472 1429 177.3814340P + 1428 177.381434 12141 0.764699 1472 Then I convert it to a 'zoo', removing the 'time' column from the input last_col- ncol (mat) range_cols - 2:last_col matrix_values - mat [,range_cols] z - zoo (matrix_values, mat $ time) So far, everything is fine. But then I need to apply a function to samples taken every t seconds, ie, the mean every 10 seconds. And, AFAIK, rapply needs a constant 'width', something that can only be obtained with regular time series... But if I try to get a 'ts' object from my 'zoo' I get as.ts (z) Error in if (del == 0 to == 0) return(to) : missing value where TRUE/FALSE needed Ummm, do you know if there is any error in my 'zoo' object? Has it the right format? Or, could I avoid the conversion and use the 'rapply' function but on time intervals (instead of points intervals)? Thank you very much. Alvaro On 12 Dec 2005, at 16:31, Gabor Grothendieck wrote: First we generate some sample data x and its times tt. Then using the zoo package we create an irregularly spaced time series. Now if you want a regularly spaced time series convert it to ts class. After loading zoo as shown below, the R command vignette(zoo) gives more info. x - 1:4 tt - c(1, 3, 4, 6) library(zoo) x.zoo - zoo(x, tt) # irregularly spaced time series x.ts - as.ts(x.zoo) # regular time series with NAs On 12/12/05, Alvaro Saurin [EMAIL PROTECTED] wrote: Hi, I am a R begginer and I have a small problem with time series. I was wondering if someone could help me I am collecting data from packets going through a network, and using R for obtaining some simple statistics of connections. However, my data is not collected at a constant frequency, so I would like to create a evenly spaced TS from my traces, using the minimum time difference between two samples as the period for my new TS, and filling the gaps with NA (or 0s). I think ther must be some simple solution for this... Anyone could help me? Thanks in advance. -- Alvaro Saurin [EMAIL PROTECTED] [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html -- Alvaro Saurin [EMAIL PROTECTED] [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Incomplete Beta
Is there any function available in R for computing the incomplete Beta function? I'll appreciate any suggestion -- Albert Sorribas Grup de Bioestadística i Biomatematica Departament de Ciències Mèdiques Bàsiques Universitat de Lledia tel: +34 973 702 406 FAX: +34 973 702 426 Home page: http://www.udl.es/Biomath/Group __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Generation of missiing values in a time serie...
Your variable mat is not a matrix; its a data frame. Check it with: class(mat) Here is an example: x - cbind(A = 1:4, B = 5:8) tt - c(1, 3:4, 6) library(zoo) x.zoo - zoo(x, tt) x.ts - as.ts(x.zoo) On 12/13/05, Alvaro Saurin [EMAIL PROTECTED] wrote: Hi, First, thank you. It 'almost' works: I can convert a matrix to a 'zoo' object, but it can not convert it to a 'ts'. The sitruation is this: I load a set of samples with h_types - list (0, 0, character(0), character(0), 0, 0, 0, 0, 0) h_names - list (time, flow, type, dir, seq, ts, x, rtt, size) pcks_file - pipe (grep ' P ' server.dat), r) pcks- scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) Read 1429 records mat - data.frame (pcks) colnames (mat) - h_names mat time flow type dir seq ts x rtt size 1 1.0008930P +0 1.000893 1472 0.00 1472 2 1.5144540P +1 1.514454 2944 0.513142 1472 3 2.0150930P +2 2.015093 2944 0.513142 1472 4 2.5150250P +3 2.515025 4806 0.504488 1472 5 2.8219760P +4 2.821976 5730 0.496728 1472 6 3.0789310P +5 3.078931 5832 0.489744 1472 7 3.3318970P +6 3.331897 5832 0.489744 1472 [...] 1425 176.9255040P + 1424 176.925504 12141 0.764699 1472 1426 177.0394890P + 1425 177.039489 12141 0.764699 1472 1427 177.1534690P + 1426 177.153469 12141 0.764699 1472 1428 177.2674640P + 1427 177.267464 12141 0.764699 1472 1429 177.3814340P + 1428 177.381434 12141 0.764699 1472 Then I convert it to a 'zoo', removing the 'time' column from the input last_col- ncol (mat) range_cols - 2:last_col matrix_values - mat [,range_cols] z - zoo (matrix_values, mat $ time) So far, everything is fine. But then I need to apply a function to samples taken every t seconds, ie, the mean every 10 seconds. And, AFAIK, rapply needs a constant 'width', something that can only be obtained with regular time series... But if I try to get a 'ts' object from my 'zoo' I get as.ts (z) Error in if (del == 0 to == 0) return(to) : missing value where TRUE/FALSE needed Ummm, do you know if there is any error in my 'zoo' object? Has it the right format? Or, could I avoid the conversion and use the 'rapply' function but on time intervals (instead of points intervals)? Thank you very much. Alvaro On 12 Dec 2005, at 16:31, Gabor Grothendieck wrote: First we generate some sample data x and its times tt. Then using the zoo package we create an irregularly spaced time series. Now if you want a regularly spaced time series convert it to ts class. After loading zoo as shown below, the R command vignette(zoo) gives more info. x - 1:4 tt - c(1, 3, 4, 6) library(zoo) x.zoo - zoo(x, tt) # irregularly spaced time series x.ts - as.ts(x.zoo) # regular time series with NAs On 12/12/05, Alvaro Saurin [EMAIL PROTECTED] wrote: Hi, I am a R begginer and I have a small problem with time series. I was wondering if someone could help me I am collecting data from packets going through a network, and using R for obtaining some simple statistics of connections. However, my data is not collected at a constant frequency, so I would like to create a evenly spaced TS from my traces, using the minimum time difference between two samples as the period for my new TS, and filling the gaps with NA (or 0s). I think ther must be some simple solution for this... Anyone could help me? Thanks in advance. -- Alvaro Saurin [EMAIL PROTECTED] [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html -- Alvaro Saurin [EMAIL PROTECTED] [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Constrained Log-Likelihood with SQP Solver
Dear R-Users, I'm searching for somebody who can support me or even likes to collaborate with me in setting up an R-package for constrained maximim log-likelihood parameter estimation. For example fitting the parameters of a MA(1)-APARCH(1,1) model for a time series of 17'000 points (e.g. the famous Ding-Granger-Engle mode) takes about 10 minutes with the existing optimization algorithms available under R. Modern state of the art algorithms, like SQP algorithms as implemented in Gauss, Matlab, Ox, take about a few seconds. I tested this finding with a free constrained SQP solver written in FORTRAN under R and found these results confirmed. I got the results in a few seconds instead of a few minutes! Now I'm looking for a collegue who has the experience in implementing FORTRAN Optimization Code in R, calling the objective function and optionally gradient and hessian from R functions. I have already inspected a lot of Fortran, C, and R sources from the base package, but I didn't succeed so far with a reasonable effort. Many thanks in advance Diethelm Wuertz __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] superimpose density line over hist
Hi all, I'm trying to superimpose a rchisq density line over a histogram with something like: hist(alnlength) lines(density(rchisq(length(alnlength), 4)),col=red) But the rchisq line won't appear anywhere, Anyone knows what I am missing here? Thanks in advance, Albert. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] How to make a plot?
Does anyone have an idea of how to make a chart in R like the ones here that use a graphic: http://bigpicture.typepad.com/comments/images/slide1.gif __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] sample matrix as a new object
Please, I´d like to store this sample matrix as a new object. How can I do this ? pulse - c(67, 67, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 71, 71, 72, 72, 73, 74) m - NULL x - 0 for (i in 1:5) { x - sample(pulse,3) m - mean(x) cat(x,m,\n) } Thanks, Mauricio __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Technique for reading large sparse fwf data file
I should have also noted in this email how I have allocated memory and an error that appears. I'm using Windows, so as in FAQ 2.2 I did C:\Program Files\R\R-2.2.0\bin\Rgui.exe --sdi --max-mem-size=2Gb # Check memory size in R example(memory.size) mmry.s memory.size() [1] 11894064 mmry.s memory.size(TRUE) [1] 12500992 mmry.s round(memory.limit()/1048576, 2) [1] 2048 An interesting issue appears after trying to import the subset of the larger file (which is a csv file 75,238 KB). R indicates it has run out of memory as: Error: vector memory exhausted (limit reached?) Error: vector memory exhausted (limit reached?) So, when I then try to quit R, it doesn't allow me to. Here is a copy and paste from my workspace. quit() Error: vector memory exhausted (limit reached?) quit() Error: recursive default argument reference quit() Error: vector memory exhausted (limit reached?) Clearly, enough memory is allocated to handle this file. But, I also wonder why R then locks and I need to do a forced shut down. -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Doran, Harold Sent: Tuesday, December 13, 2005 5:33 AM To: r-help@stat.math.ethz.ch Subject: [R] Technique for reading large sparse fwf data file Dear list: A datafile was sent to me that is very large (92890 x 1620) and is *very* sparse. Instead of leaving the entries with missing data blank, each cell with missing data contains a dot (.) The data are binary in almost all columns, with only a few columns containing whole numbers, which I believe requires 2 bytes for the binary and 4 for the others. So, by my calculations (assuming 4 bytes for all cells to create an upperbound) I should need around 92890 * 1620 * 4 = 574MB to read in these data and about twice that for analyses. My computer has 3GB. But, I am unable to read in the file even though I have allocated sufficient memory to R for this. My first question is do the dots in the empty cells consume additional memory? I am assuming the answer is yes and believe I should remove them before I do the read in. Because my data are in a fixed width format file, I can open the file in a text editor and find and replace all dots with nothing. Then, I should retry the read in process? Maybe this will work? I created a smaller data file (~ 14000 * 1620) in SAS and tried to import this subset (it still had the dots), but R still would not allow for me to do so. I could use a little guidance as I think I have allocated sufficient memory to read in a datafile assuming my calculations are right. Does anyone have any thoughts on a strategy? Harold [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] help with multivariate analysis
dear R users, I need some help for multivariate analysis. I have 2 anaesthetic treatment groups (20 patients/group) where I register heart frequency and pressure for 60 min (repeated measures every 5 minutes). I would like to perform a test to check if treatments are different in controlling freq and pressures during the anaesthesia, but i would like to have also an overall measure and not only multiple p for different time intervals. I also think I should choose a test in which time is meaningful since the measures are not simple repeated measurements but measurements taken at specific time points. 1 million dollar question how to do in R? thanks in advance Dr Matteo Vidali Dep. of Medical Sciences University of East Piedmont A. Avogadro ITALY -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Bivariate Splines in R
From: Martin Maechler Cal == Cal Stats [EMAIL PROTECTED] on Mon, 12 Dec 2005 10:25:38 -0800 (PST) writes: Cal Hi.., is there a function in R to fit bivariate splines Cal ? I came across 'polymars' (POLSPLINE) and 'mars' Cal (mda) packages. Are these the one to use or are there Cal other specific commands? I'd recommend to use gam(y ~ s(x1,x2)) from the recommended package 'mgcv'. help(gam) has many examples, some of which using bivariate splines. The `gss' package might be worth trying, too. Cheers, Andy Cal Thanks. Cal Harsh you're welcome; Martin Maechler, ETH Zurich __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] superimpose density line over hist
One thing to check: Are the data in the same range as chisq with 4 df? Andy From: Albert Vilella Hi all, I'm trying to superimpose a rchisq density line over a histogram with something like: hist(alnlength) lines(density(rchisq(length(alnlength), 4)),col=red) But the rchisq line won't appear anywhere, Anyone knows what I am missing here? Thanks in advance, Albert. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Question
On 13-Dec-05 Uwe Ligges wrote: Davia Cox wrote: [...] My question is in order to run a true random sample, We have to disappoint you: Your computer cannot generate true random samples. At least on Linux (and probably most Unix) systems, /dev/random must be a pretty good approximation (provided you don't over-work it ... ). Oh for the days of shot noise! Best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 13-Dec-05 Time: 14:16:59 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] QQ plot for deviance residuals in generalized linear models
Hello In an article in J. Computational and Graphical Statistics, v 13, p. 36-47, Ben and Yohai propose quantile quantile plots for deviance residuals of generalized linear models, Normal quantile plots of these residuals are unsatisfactory, but B and Y propose a QQ plot of the deviance residuals against F_D(d) = P(D leq d) = P(c(y,x,Beta) leq d) they note that there are functions in S-plus 6.0 to generate these QQ plots for logistic and Poisson regression. Has anyone implemented these in R? TIA Peter Peter L. Flom, PhD Assistant Director, Statistics and Data Analysis Core Center for Drug Use and HIV Research National Development and Research Institutes 71 W. 23rd St http://cduhr.ndri.org www.peterflom.com New York, NY 10010 (212) 845-4485 (voice) (917) 438-0894 (fax) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] superimpose density line over hist
try hist(alnlength, prob = TRUE) moreover I think that you do not need to sample from the Chi-squared to draw its density; you could use dchisq(seq(...), 4)) instead. I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://www.med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm - Original Message - From: Albert Vilella [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Tuesday, December 13, 2005 2:23 PM Subject: [R] superimpose density line over hist Hi all, I'm trying to superimpose a rchisq density line over a histogram with something like: hist(alnlength) lines(density(rchisq(length(alnlength), 4)),col=red) But the rchisq line won't appear anywhere, Anyone knows what I am missing here? Thanks in advance, Albert. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] superimpose density line over hist
Le 13.12.2005 14:23, Albert Vilella a écrit : Hi all, I'm trying to superimpose a rchisq density line over a histogram with something like: hist(alnlength) lines(density(rchisq(length(alnlength), 4)),col=red) But the rchisq line won't appear anywhere, Anyone knows what I am missing here? Thanks in advance, Albert. Hi Albert, A few comments : - your code should be reproductible, otherwise it is useless. (that recommandation is on the posting guide) - that question is a top ten question on that list, go to the archive and you will find answers. (also posting guide) BTW, it should be a FAQ and what about an example of overlaying in hist help page ? - then if you want to superimpose an histogram and a density curve, they should be on the same scale, it is not the case since you missed the argument freq in your hist call, it should be : R hist(alnlength, freq=FALSE) - Why do you simulate points to draw the density line ? Give a shot at : R curve(dchisq, df=4, col=red) - That might interrest you : http://addictedtor.free.fr/graphiques/search.php?q=hist Romain -- visit the R Graph Gallery : http://addictedtor.free.fr/graphiques mixmod 1.7 is released : http://www-math.univ-fcomte.fr/mixmod/index.php +---+ | Romain FRANCOIS - http://francoisromain.free.fr | | Doctorant INRIA Futurs / EDF | +---+ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to make a plot?
That can be done very easy using a program that can deal with layers, for example the Gimp. Just make one layer containing the barplot and one layer containing the background image. Then you can do all kind of nice things. Directly in R, I don't have a clue, but it wouldn't surprise me if someone will show you an example ;) JeeBee. On Tue, 13 Dec 2005 08:32:15 -0500, Gabor Grothendieck wrote: Does anyone have an idea of how to make a chart in R like the ones here that use a graphic: http://bigpicture.typepad.com/comments/images/slide1.gif __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Incomplete Beta
On Tue, 13 Dec 2005, Albert Sorribas wrote: Is there any function available in R for computing the incomplete Beta function? pbeta(). The incomplete Beta function is the cdf of the Beta distribution -thomas I'll appreciate any suggestion -- Albert Sorribas Grup de Bioestadística i Biomatematica Departament de Ciències Mèdiques Bàsiques Universitat de Lledia tel: +34 973 702 406 FAX: +34 973 702 426 Home page: http://www.udl.es/Biomath/Group __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Thomas Lumley Assoc. Professor, Biostatistics [EMAIL PROTECTED] University of Washington, Seattle__ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Problem with understanding output of Cox model
On Tue, 13 Dec 2005, May, Roel wrote: Hi all, I am using a 'tricked' Cox Hazard regression model for discrete choice habitat modelling. However, I'm having a hard time understanding the meaning of the first line the following part of the summary() output: Rsquare= 0.307 (max possible= 0.475 ) Likelihood ratio test= 91.8 on 12 df, p=2.23e-14 Wald test = 26.3 on 12 df, p=0.00977 Score (logrank) test = 58.6 on 12 df, p=4.03e-08 Does anyone know how I can read the 'max possible' R-square? What does it signify? It doesn't signify anything very useful. It is the proportional reduction in deviance. -thomas __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] interruption when pasting code into R under linux
hello, has anyone come across the following rather mysterious problem: when pasting large bits of code (100 and more lines) into the R console with the central mouse button (under linux), only part of the code is pasted, and the text interrupts somewhere arbitrarily. It does not happen when smaller bits are pasted subsequently. I use linux suse 9.3 with the latest version of R Sebastian Leuzinger Institute of Botany, University of Basel Schönbeinstr. 6 CH-4056 Basel ph0041 (0) 61 2673511 fax 0041 (0) 61 2673504 email [EMAIL PROTECTED] web http://pages.unibas.ch/botschoen/leuzinger __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] sample matrix as a new object
You mean writing to a file? Maybe you should read the R Data Import/Export Manual http://cran.r-project.org/doc/manuals/R-data.html or as pdf http://cran.r-project.org/doc/manuals/R-data.pdf You might also want to read the manual pages of these two commands: help('dump') help('write.table') JeeBee. On Tue, 13 Dec 2005 11:31:33 -0300, Carlos Mauricio Cardeal Mendes wrote: Please, I´d like to store this sample matrix as a new object. How can I do this ? pulse - c(67, 67, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 71, 71, 72, 72, 73, 74) m - NULL x - 0 for (i in 1:5) { x - sample(pulse,3) m - mean(x) cat(x,m,\n) } Thanks, Mauricio __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] superimpose density line over hist
On 13-Dec-05 Albert Vilella wrote: Hi all, I'm trying to superimpose a rchisq density line over a histogram with something like: hist(alnlength) lines(density(rchisq(length(alnlength), 4)),col=red) But the rchisq line won't appear anywhere, Anyone knows what I am missing here? Well, it's certainly somewhere, but not within the window of your graphic! What you're missing is a) Breaks in 'hist' b) Allowing for (1) the widths of the bins, (2) the size of the sample, when you drawn the curve. Example (this should work most of the time, but if it doesn't because you have a sample that goes out of range just try again): alnlength-rchisq(1000,4) x-0.25*(0:100) hist(alnlength,breaks=0.25*(0:100)) lines(x,1000*dchisq(x,4)*0.25) Note that to get the 'lines' matching the histogram you have to a) multiply the density by the binswidth to get probabilities; b) multiply by the sample size so that the vertical scale of the curve matches that of the histogram (which here is showing counts). If you don't want to set the breakpoints explicitly but leave it to 'hist', you can subsequently extract the breakpoints from the 'hist' object and carry on from there. Hoping this helps, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 13-Dec-05 Time: 15:50:21 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Question
On Tue, 13 Dec 2005 [EMAIL PROTECTED] wrote: On 13-Dec-05 Uwe Ligges wrote: Davia Cox wrote: [...] My question is in order to run a true random sample, We have to disappoint you: Your computer cannot generate true random samples. At least on Linux (and probably most Unix) systems, /dev/random must be a pretty good approximation (provided you don't over-work it ... ). See library(accuracy) ?trueRandom (for some definition of 'true'). -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] superimpose density line over hist
On 13-Dec-05 Albert Vilella wrote: Hi all, I'm trying to superimpose a rchisq density line over a histogram with something like: hist(alnlength) lines(density(rchisq(length(alnlength), 4)),col=red) But the rchisq line won't appear anywhere, Anyone knows what I am missing here? Following up my earlier reply, and more in line with your precise question above: alnlength-rchisq(1000,4) hist(alnlength,breaks=0.25*(0:100)) j-density(rchisq(1000,4)) lines(j$x,1000*0.25*j$y) Best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 13-Dec-05 Time: 15:59:50 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Forcing model parameters to a constant (in GLMs)
Hi, I'm looking for a way to bind specific model parameters to a constant (in my case: 0) before fitting the model. I'm fitting multinomial logistic regression models (via glm()). Assuming a model with a categorial response variable c {c1,c2,c3} that is a matrix like c2:b2 b2.1 b2.2 b2.3 c3:b3 b3.1 b3.2 b3.3, I would like to force, e.g. b2.1, b3.1 and b2.3 to 0 before estimating the other parameters. I'm aware of a similar - unresolved - question on list list back in July: https://stat.ethz.ch/pipermail/r-help/2005-July/075042.html The offset in an interaction term doesn't seem to do what I want - as far as I can see, it manipulates the intercept, and it's additive (i.e. an intercept is estimated anyways). Thanks David -- David Reitter - ICCS/HCRC, Informatics, University of Edinburgh __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] superimpose density line over hist
This certainly did the trick, Thanks Ted, Sean, Romain, Dimitris, Lian and Andy, alnlength-rchisq(1000,4) x-0.25*(0:100) hist(alnlength,breaks=0.25*(0:100)) lines(x,1000*dchisq(x,4)*0.25) And apologies for my newbieness in the posting, Albert. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Polytopic Vector Analysis (PVA)
I am curious whether anybody has or is developing this capability within R. I have not found any software yet that has this capability and I am not sure whether it is too new a method and nobody is actually using it, or if there are other means to get the same analysis that I do not know of. Any information regarding developments or use of this method would be helpful. Melanie Edwards [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] sample matrix
Hi Mauricio, Although you question is not clear, I'm supposing that you want to save as a new object (matrix) whatever is printed out from cat(x,m,\n). If this is the case, then your result is a matrix with 5 rows and 4 columns. For each row, the first 3 values are x and the last value is m. Maybe this will help you. pulse - c(67, 67, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 71, 71, 72, 72, 73, 74) m - NULL x - 0 result- matrix(0, 5,4) for (i in 1:5) { x - sample(pulse,3) result[i,1:3]-x m - mean(x) result[i,4]-m cat(x,m,\n) } Hopefully the object that you wanted is result HTH Thanks, Mauricio __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Ferdinand Alimadhi Programmer / Analyst Harvard University The Institute for Quantitative Social Science (617) 496-0187 [EMAIL PROTECTED] www.iq.harvard.edu __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] bug in geoR (?)
I've enconuntered this problem with the last cran version of geoR: library(geoR) day - rep(1:2, each=5) coords - matrix(rep(runif(10),2), 10, 2) data - rnorm(10) data[1] - NA as.geodata(cbind(coords, data, day), realisations=4) as.geodata: 1 points removed due to NA in the data Errore in as.geodata(cbind(coords, data, day), realisations = 4) : realisations and coords have incompatible dimensions The problem disappear if I remove the NA manually from the dataset before passing to as.geodata. I.e.: as.geodata(cbind(coords, data, day)[2:10,], realisations=4) works. Antonio, Fabio Di Narzo. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Incomplete Beta
On 13-Dec-05 Thomas Lumley wrote: On Tue, 13 Dec 2005, Albert Sorribas wrote: Is there any function available in R for computing the incomplete Beta function? pbeta(). The incomplete Beta function is the cdf of the Beta distribution But don't forget to multiply by beta(,): ibeta(x,a,b) - function(x,a,b){ pbeta(x,a,b)*beta(a,b) } ! Best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 13-Dec-05 Time: 16:18:35 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Generation of missiing values in a time serie...
On 13 Dec 2005, at 13:08, Gabor Grothendieck wrote: Your variable mat is not a matrix; its a data frame. Check it with: class(mat) Here is an example: x - cbind(A = 1:4, B = 5:8) tt - c(1, 3:4, 6) library(zoo) x.zoo - zoo(x, tt) x.ts - as.ts(x.zoo) Fixed, but anyway it fails: h_types - list (0, 0, NULL, NULL, 0, 0, 0, 0, 0) h_names - list (time, flow, seq, ts, x, rtt, size) pcks_file - pipe (grep ' P ' server.dat, r) pcks- scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) mat_df - data.frame (pcks[1:2], pcks[5:9]) mat - as.matrix (mat_df) colnames (mat) - h_names class (mat) [1] matrix z - zoo (mat, mat [,time]) z z time flow seq ts xrtt size 1.0009 1.000893 0.00 0.00 1.000893 1472.00 0.00 1472.00 1.5145 1.514454 0.00 1.00 1.514454 2944.00 0.513142 1472.00 2.0151 2.015093 0.00 2.00 2.015093 2944.00 0.513142 1472.00 2.5152.515025 0.00 3.00 2.515025 4806.00 0.504488 1472.00 2.8222.821976 0.00 4.00 2.821976 5730.00 0.496728 1472.00 [...] as.ts (z) Error in if (del == 0 to == 0) return(to) : missing value where TRUE/FALSE needed Any idea? Thanks for your help. Alvaro -- Alvaro Saurin [EMAIL PROTECTED] [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] convergence error (lme) which depends on the version of nlme (?)
On 12/12/05, Leo Gürtler [EMAIL PROTECTED] wrote: Dear list members, the following hlm was constructed: hlm - groupedData(laut ~ design | grpzugeh, data = imp.not.I) the grouped data object is located at and can be downloaded: www.anicca-vijja.de/lg/hlm_example.Rdata The following works: library(nlme) summary( fitlme - lme(hlm) ) with output: ... AIC BIClogLik 425.3768 465.6087 -197.6884 Random effects: Formula: ~design | grpzugeh Structure: General positive-definite StdDevCorr (Intercept) 0.3772478 (Intr) dsgn:8 dsgn:7 designmit:8 0.6776543 0.183 designohne:7 0.6619983 -0.964 0.086 designohne:8 1.0680576 -0.966 0.077 1.000 Residual 1.3468816 Notice that the estimated variance-covariance matrix for the random effects is singular (a correlation of +1.000). The estimates of the parameters in the model are on the boundary and it is not a proper linear mixed model. The definition of a linear mixed model (or at least my definition) requires that the variance-covariance matrix of the random effects be positive definite and this one is only positive semidefinite. Fixed effects: laut ~ design Value Std.Error DF t-value p-value (Intercept) 3.857143 0.2917529 102 13.220579 0. designmit:8 -0.285714 0.4417919 102 -0.646717 0.5193 designohne:7 -0.107143 0.4383878 102 -0.244402 0.8074 designohne:8 0.607143 0.5408713 102 1.122527 0.2643 Correlation: (Intr) dsgnm:8 dsgn:7 designmit:8 -0.451 designohne:7 -0.775 0.363 designohne:8 -0.763 0.304 0.699 Standardized Within-Group Residuals: Min Q1Med Q3Max -2.5074669 -0.4530573 0.1755326 0.5837670 2.374 Number of Observations: 112 Number of Groups: 7 The following does _not_ work and leads to a convergence error: fitlme1 - lme(laut ~ design, random = ~ design | grpzugeh, data = hlm) Fehler in lme.formula(laut ~ design, random = ~design | grpzugeh, data = hlm) : iteration limit reached without convergence (9) This was tried with R : Copyright 2005, The R Foundation for Statistical Computing Version 2.2.0 (2005-10-06 r35749) Using another R version (2.1.0, also windows with nlme version built under R 2.1.1) , it works. Thus, what's the problem then? I tried without the random effects, i.e. random = ~ 1 | grpzugeh This works. Comparing both calls on the version R2.1.0 that goes well, the following differences in the output of the random effects can be identified: summary( fitlme - lme(hlm) ) -- Random effects: ... Structure: General positive-definite /-- compared to summary(lme(laut ~ design, random = ~ design | grpzugeh, data = hlm)) -- Random effects: ... Structure: General positive-definite, Log-Cholesky parametrization /-- The estimates of the fixed effects are similar, the S.E.s not. The random effects are different, too. AIC/BIC/logLik are slightly different. Thus my question: 1) Do I have overseen a switch for the structure of the random effects? Is something wrong with the call/ formular? 2) What is the cause of the convergence error which seems to depend on the built of R/nlme? Thank you very much. Best wishes, leo gürtler As Dieter indicated in his response, the more current function lmer from the lme4 package (actually it's in the Matrix package but it would be in the lme4 package if a certain capability related to packages were available) is preferred to lme. Fitting your model with the control options for verbose output in both the EM and nlminb iterations produces (fm1 - lmer(laut ~ design + (design | grpzugeh), hlm, control = list(msV=1,EMv=1))) EM iterations 0 407.611 ( 6.0 1.5 1.5 1.5 0.0 0.0 0.0 0.0 0.0 0.0: -0.409-1.07-2.19 -0.969 -0.0472 -0.344 -0.0282 -0.491 -0.1630.941) 1 402.107 ( 10.4497 1.95422 3.22722 2.22340 0.196761 1.02069 0.00757874 1.13553 0.110538 -0.685820: -0.122 -0.550 -0.567 -0.181 0.0294 -0.112 -0.00789 -0.204 -0.01840.361) 2 399.890 ( 14.8865 2.30933 5.18627 2.99207 0.242029 2.06595 -0.0167045 2.18847 0.173349 -1.51318: -0.0497 -0.331 -0.209 0.00812 0.0311 -0.0667 -0.00119 -0.129 0.009420.222) 3 398.756 ( 19.0686 2.58783 7.19874 3.76967 0.147926 3.04342 -0.0686073 3.14563 0.190736 -2.40480: -0.0224 -0.217 -0.0877 0.0682 0.0250 -0.0508 0.00304 -0.0968 0.01780.166) 4 398.074 ( 23.0243 2.81061 9.22509 4.55494 -0.0495774 3.95755 -0.140106 4.03331 0.174045 -3.33077:-0.00975 -0.150 -0.0362 0.0864 0.0192 -0.0422 0.00605 -0.0784 0.02130.134) 5 397.620 ( 26.8048 2.99284 11.2543 5.34938 -0.321835 4.82191 -0.225236 4.87317 0.132590 -4.27703:-0.00344 -0.108 -0.0119 0.0876 0.0145 -0.0360 0.00810 -0.0653 0.02290.111) 6 397.297 ( 30.4530 3.14530 13.2827 6.15353 -0.648070 5.64798 -0.319808
Re: [R] Generation of missiing values in a time serie...
Please provide a reproducible example. Note that dput(x) will output an R object in a way that can be copied and pasted into another session. On 12/13/05, Alvaro Saurin [EMAIL PROTECTED] wrote: On 13 Dec 2005, at 13:08, Gabor Grothendieck wrote: Your variable mat is not a matrix; its a data frame. Check it with: class(mat) Here is an example: x - cbind(A = 1:4, B = 5:8) tt - c(1, 3:4, 6) library(zoo) x.zoo - zoo(x, tt) x.ts - as.ts(x.zoo) Fixed, but anyway it fails: h_types - list (0, 0, NULL, NULL, 0, 0, 0, 0, 0) h_names - list (time, flow, seq, ts, x, rtt, size) pcks_file - pipe (grep ' P ' server.dat, r) pcks- scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) mat_df - data.frame (pcks[1:2], pcks[5:9]) mat - as.matrix (mat_df) colnames (mat) - h_names class (mat) [1] matrix z - zoo (mat, mat [,time]) z z time flow seq ts xrtt size 1.0009 1.000893 0.00 0.00 1.000893 1472.00 0.00 1472.00 1.5145 1.514454 0.00 1.00 1.514454 2944.00 0.513142 1472.00 2.0151 2.015093 0.00 2.00 2.015093 2944.00 0.513142 1472.00 2.5152.515025 0.00 3.00 2.515025 4806.00 0.504488 1472.00 2.8222.821976 0.00 4.00 2.821976 5730.00 0.496728 1472.00 [...] as.ts (z) Error in if (del == 0 to == 0) return(to) : missing value where TRUE/FALSE needed Any idea? Thanks for your help. Alvaro -- Alvaro Saurin [EMAIL PROTECTED] [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] fSeries
I'm trying to use garchFit from fSeries, with Student or Skewed Student conditionnal distribution. Let's say that eps (vector) is my series of daily log-returns: data(EuStockMarkets) eps = diff(log(EuStockMarkets[,CAC])) library(fSeries) g = garchFit(series = eps, formula.var = ~garch(2,2), cond.dist = dstd) s = [EMAIL PROTECTED] All the coefficients are ok (checked with SAS 9.1) except nu (degrees of freedom of the student) and the log-likelyhood. I've really checked everything and can't find the estimated series sigma (volatility) and eta, such that eps = sigma * eta and eta is centered and reduced... I've tryed combinations of all s$x,s$h,s$z and nothing looks looks correct. Also, is it possible to fit EGARCH and TGARCH with R ? If anyone ever managed to make it work, i'd be grateful ;-) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Question
On 13-Dec-05 Prof Brian Ripley wrote: On Tue, 13 Dec 2005 [EMAIL PROTECTED] wrote: On 13-Dec-05 Uwe Ligges wrote: Davia Cox wrote: [...] My question is in order to run a true random sample, We have to disappoint you: Your computer cannot generate true random samples. At least on Linux (and probably most Unix) systems, /dev/random must be a pretty good approximation (provided you don't over-work it ... ). See library(accuracy) ?trueRandom (for some definition of 'true'). Spot on! Thanks for this -- I wasn't aware of it previously. The URL to HotBits in ?trueRandom is worth following, and the descriptions in How it Works: http://www.fourmilab.ch/hotbits/how.html are entertaining reading. Best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 13-Dec-05 Time: 16:47:37 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] help with writing function
On Tue, 13 Dec 2005, Oarabile Molaodi wrote: I'm trying to write a function that takes a vector of length n and then takes the first value of the vector i.e j=1 and forms a new vector of length n (i.e replicate the first value n times). This function will then calculate the absoulte difference of the original vector and the new vector and store the results omitting the difference between the value and itself. This function should be able to repeat the procedure for each of the j's i.e j=2 to n. The results should all be stored together. Below is what I've tried so far but it seems to work only for j=1 . Your help will be highly appreciated. IED-function(risk){ n-length(risk) i-c(1:n) Diff-numeric() for(j in 1:n){ relrisk-risk relrisk[i]-relrisk[j] Difference-abs(risk-relrisk) Difference-Difference[-c(1:j)] Difference-append(Diff,Difference) return(Difference) } } How about IED - function(risk){ n-length(risk) Diff-numeric(n) for(j in 1:n){ relrisk-rep(risk[j],n) Diff[j]-sum(abs(risk-relrisk)[-j]) } Diff } -- Robert Burrows, PhD New England Biometrics [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How to make a plot?
Le 13.12.2005 14:32, Gabor Grothendieck a écrit : Does anyone have an idea of how to make a chart in R like the ones here that use a graphic: http://bigpicture.typepad.com/comments/images/slide1.gif Hi Gabor, Here's a first draft. The following code will produce a pdf (i only found transparancy support in pdf devices) which will be converted afterwards into a png file using (thanks to imageMagick) : $ convert output.pdf output.png Everything is there : http://addictedtor.free.fr/misc/gabor/ Background image : http://addictedtor.free.fr/misc/gabor/cruise.pnm Code : http://addictedtor.free.fr/misc/gabor/gabor.R Pdf output : http://addictedtor.free.fr/misc/gabor/output.pdf Png output : http://addictedtor.free.fr/misc/gabor/output.png Romain Code : # require(pixmap) require(grid) x - read.pnm('cruise.pnm') pdf('output.pdf', width=6, height=4, version=1.4) par(mar=c(0,0,0,0)) plot(x) # base graphics y - c(6, 6.5, 7, 8, 8.5, 8.2, 10, 9.6, 9.7, 9) # some data like in the picture you gave # now the grid stuff pushViewport(viewport(xscale=c(0,10), yscale=c(0,10))) grid.rect(x=0:9, y=0, width=1, height=y, default.units=native, gp=gpar(fill=white, alpha=0.7, col=gray, lwd=2), just=c(left,bottom)) grid.rect(x=0:9, y=y, width=1, height= unit(1, npc) - unit(y, native) , default.units=native , gp=gpar(fill=white, col=white), just=c(left,bottom)) grid.lines(x=c(0,10), y=c(5, 5), default.units=native, gp=gpar(lwd=2, col=white, lty=dotted)) grid.lines(x=c(0,10), y=c(10, 10), default.units=native, gp=gpar(lwd=2, col=gray, lty=dotted)) popViewport() dev.off() # -- visit the R Graph Gallery : http://addictedtor.free.fr/graphiques mixmod 1.7 is released : http://www-math.univ-fcomte.fr/mixmod/index.php +---+ | Romain FRANCOIS - http://francoisromain.free.fr | | Doctorant INRIA Futurs / EDF | +---+ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Generation of missiing values in a time serie...
I think I have found the error. It appears when there are two entries with the same time. Using as input file: - CUT # Output format for PCKs: # TIME FLOW P [+-] SEQ TS X RTT SIZE # 123.125683 0 P + 967 123.125683 13394 0.798205 1472 123.241137 0 P + 968 123.241137 12680 0.796258 1472 123.241137 0 P + 969 123.241137 12680 0.796258 1472 123.472631 0 P + 970 123.472631 12680 0.796258 1472 123.588613 0 P + 971 123.588613 12680 0.796258 1472 123.704594 0 P + 972 123.704594 12680 0.796258 1472 - CUT I run fhe following code: - CUT h_types - list (0, 0, NULL, NULL, 0, 0, 0, 0, 0) h_names - list (time, flow, seq, ts, x, rtt, size) pcks_file- pipe (grep ' P ' data, r) pcks - scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) mat_df - data.frame (pcks[1:2], pcks[5:9]) mat - as.matrix (mat_df) colnames (mat) - h_names z - zoo (mat, mat [,time]) - CUT The dput of 'z' shows: - CUT structure(c(123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594, 0, 0, 0, 0, 0, 0, 967, 968, 969, 970, 971, 972, 123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594, 13394, 12680, 12680, 12680, 12680, 12680, 0.798205, 0.796258, 0.796258, 0.796258, 0.796258, 0.796258, 1472, 1472, 1472, 1472, 1472, 1472 ), .Dim = c(6, 7), .Dimnames = list(c(1, 2, 3, 4, 5, 6), c(time, flow, seq, ts, x, rtt, size)), index = structure(c(123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594), .Names = c(1, 2, 3, 4, 5, 6)), class = zoo) - CUT If I try a 'as.ts(z)', it fails. If I remove the duplicate entry, I can convert it to a TS with no problem. Is this made intentionally? Because then I have to filter the input matrix... But, anyway, the output matrix, after filtering, doesn't seem regular: - CUT as.ts (z) Time Series: Start = 1 End = 5 Frequency = 1 time flow seq ts x rtt size 1 123.12570 967 123.1257 13394 0.798205 1472 2 123.24110 969 123.2411 12680 0.796258 1472 3 123.47260 970 123.4726 12680 0.796258 1472 4 123.58860 971 123.5886 12680 0.796258 1472 5 123.70460 972 123.7046 12680 0.796258 1472 Warning message: ‘x’ does not have an underlying regularity in: as.ts.zoo(z) - CUT Weird... On 13 Dec 2005, at 16:33, Gabor Grothendieck wrote: Please provide a reproducible example. Note that dput(x) will output an R object in a way that can be copied and pasted into another session. On 12/13/05, Alvaro Saurin [EMAIL PROTECTED] wrote: On 13 Dec 2005, at 13:08, Gabor Grothendieck wrote: Your variable mat is not a matrix; its a data frame. Check it with: class(mat) Here is an example: x - cbind(A = 1:4, B = 5:8) tt - c(1, 3:4, 6) library(zoo) x.zoo - zoo(x, tt) x.ts - as.ts(x.zoo) Fixed, but anyway it fails: h_types - list (0, 0, NULL, NULL, 0, 0, 0, 0, 0) h_names - list (time, flow, seq, ts, x, rtt, size) pcks_file - pipe (grep ' P ' server.dat, r) pcks- scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) mat_df - data.frame (pcks[1:2], pcks[5:9]) mat - as.matrix (mat_df) colnames (mat) - h_names class (mat) [1] matrix z - zoo (mat, mat [,time]) z z time flow seq ts xrtt size 1.0009 1.000893 0.00 0.00 1.000893 1472.00 0.00 1472.00 1.5145 1.514454 0.00 1.00 1.514454 2944.00 0.513142 1472.00 2.0151 2.015093 0.00 2.00 2.015093 2944.00 0.513142 1472.00 2.5152.515025 0.00 3.00 2.515025 4806.00 0.504488 1472.00 2.8222.821976 0.00 4.00 2.821976 5730.00 0.496728 1472.00 [...] as.ts (z) Error in if (del == 0 to == 0) return(to) : missing value where TRUE/FALSE needed Any idea? Thanks for your help. Alvaro -- Alvaro Saurin [EMAIL PROTECTED] [EMAIL PROTECTED] -- Alvaro Saurin [EMAIL PROTECTED] [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Incomplete Beta
Thomas Lumley [EMAIL PROTECTED] writes: On Tue, 13 Dec 2005, Albert Sorribas wrote: Is there any function available in R for computing the incomplete Beta function? pbeta(). The incomplete Beta function is the cdf of the Beta distribution ..except for the normalizing constant, I think. Check your references, but I suspect that you get pbeta by multiplying the incomplete Beta function by Gamma(a+b)/(Gamma(a)Gamma(b)). -thomas I'll appreciate any suggestion -- Albert Sorribas Grup de Bioestadística i Biomatematica Departament de Ciències Mèdiques Bàsiques Universitat de Lledia tel: +34 973 702 406 FAX: +34 973 702 426 Home page: http://www.udl.es/Biomath/Group __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Thomas Lumley Assoc. Professor, Biostatistics [EMAIL PROTECTED] University of Washington, Seattle__ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] creating a subset of a dataset using ifelse statement?
I am using R in a Microsoft Windows environment. I have a dataset called mp1b. I have a variable called h. h can take a value from -1 to 5. If h 1, I want to create a new dataset called mp2 that is the same as mp1b: mp2-mp1b If h 0, I want to set create a dataset mp2, where I limit the original dataset to those where mp1b$group = =h. similar to: mp2-subset (mp1b, group= = h) I have tried this ifelse statement, but it does not seem to work as expected. mp2-ifelse(h1,mp1b,subset(mp1b,cluster_q==h)) Assistance is appreciated. - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] interruption when pasting code into R under linux
Sebastian Leuzinger wrote: hello, has anyone come across the following rather mysterious problem: when pasting large bits of code (100 and more lines) into the R console with the central mouse button (under linux), only part of the code is pasted, and the text interrupts somewhere arbitrarily. It does not happen when smaller bits are pasted subsequently. I use linux suse 9.3 with the latest version of R ... because of the limited buffersize of your clipboard. I'd suggest source()-ing a file rather than pasting those many lines of code. Uwe Ligges Sebastian Leuzinger Institute of Botany, University of Basel Schönbeinstr. 6 CH-4056 Basel ph0041 (0) 61 2673511 fax 0041 (0) 61 2673504 email [EMAIL PROTECTED] web http://pages.unibas.ch/botschoen/leuzinger __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Incomplete Beta
On Tue, 13 Dec 2005 [EMAIL PROTECTED] wrote: On 13-Dec-05 Thomas Lumley wrote: On Tue, 13 Dec 2005, Albert Sorribas wrote: Is there any function available in R for computing the incomplete Beta function? pbeta(). The incomplete Beta function is the cdf of the Beta distribution But don't forget to multiply by beta(,): ibeta(x,a,b) - function(x,a,b){ pbeta(x,a,b)*beta(a,b) } ! Depends on which definition you use, as ?pbeta explains. Thomas' advice was correct rather than yours for Abramowitz and Stegun's definition, for example. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] help with writing function
Does this do what you want: IED - function(risk){ tmp - outer(risk,risk,-) tmp - abs(tmp) return(tmp[lower.tri(tmp)]) } -- Gregory L. Snow Ph.D. Statistical Data Center, IHC [EMAIL PROTECTED] (801) 408-8111 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Oarabile Molaodi Sent: Tuesday, December 13, 2005 5:00 AM To: r-help@stat.math.ethz.ch Subject: [R] help with writing function I'm trying to write a function that takes a vector of length n and then takes the first value of the vector i.e j=1 and forms a new vector of length n (i.e replicate the first value n times). This function will then calculate the absoulte difference of the original vector and the new vector and store the results omitting the difference between the value and itself. This function should be able to repeat the procedure for each of the j's i.e j=2 to n. The results should all be stored together. Below is what I've tried so far but it seems to work only for j=1 . Your help will be highly appreciated. IED-function(risk){ n-length(risk) i-c(1:n) Diff-numeric() for(j in 1:n){ relrisk-risk relrisk[i]-relrisk[j] Difference-abs(risk-relrisk) Difference-Difference[-c(1:j)] Difference-append(Diff,Difference) return(Difference) } } Oarabile __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] interruption when pasting code into R under linux
On Tue, 2005-12-13 at 15:50 +0100, Sebastian Leuzinger wrote: hello, has anyone come across the following rather mysterious problem: when pasting large bits of code (100 and more lines) into the R console with the central mouse button (under linux), only part of the code is pasted, and the text interrupts somewhere arbitrarily. It does not happen when smaller bits are pasted subsequently. I use linux suse 9.3 with the latest version of R More than likely, you are getting an input buffer overflow. The arbitrary nature of when this occurs will be dependent upon a variety of factors, including the complexity of the R code you are pasting and how fast the R interpreter can process it. At some point, a bottleneck is created and the subsequent text in the input buffer is lost. This behavior is generally intentional to avoid security risks due to buffer overruns, which is a common method exploited by folks looking to compromise a system. See http://en.wikipedia.org/wiki/Buffer_overrun for more information. You can try to increase the input buffer size for the console to see if that helps. This will be dependent upon the console app you are using and the default buffer size in place. A better solution would be to save the R code in a text file and source() that file to bring the code into R. See ?source for more information. An even better solution, if you are comfortable with emacsen, is to use ESS. This provides for a more integrated development environment. See: http://stat.ethz.ch/ESS/ for more information. HTH, Marc Schwartz __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] batch mode problem
Hi, ok, I know I should be using a later version than 1.7.1 (64 bit) but it's not in my power. So here is the problem: In my R script I declare a data.frame that consists of 40 vectors, each having 125 numeric elements. This is no problem as long as I run the sript in interactive mode, but running it in batch mode I get strange error messages. Apparently, it has to do with the size of the data.frame because reducing the data.frame to 36 vector a 125 elements, I have no problems. How can I declare my large data.frame and still run the script in batch mode? Thanks, joerg __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Labeling a range of bars in barplot?
On Tue, 2005-12-13 at 10:53 +, Dan Bolser wrote: Hi, I am plotting a distribution of (ordered) values as a barplot. I would like to label groups of bars together to highlight aspects of the distribution. The label for the group should be the range of values in those bars. As this is hard to describe, here is an example; x - rlnorm(50)*2 barplot(sort(x,decreasing=T)) y - quantile(x, seq(0, 1, 0.2)) y plot(diff(y)) That last plot is to highlight that I want to label lots of the small columns together, and have a few more labels for the bigger columns (more densely labeled). I guess I will have to turn out my own labels using low level plotting functions, but I am stumped as to how to perform the calculation for label placement. I imagine drawing several line segments, one for each group of bars to be labeled together, and putting the range under each line segment as the label. Each line segment will sit under the group of bars that it covers. Thanks for any help with the above! Cheers, Dan. Dan, Here is a hint. barplot() returns the bar midpoints: mp - barplot(sort(x, decreasing = TRUE)) head(mp) [,1] [1,] 0.7 [2,] 1.9 [3,] 3.1 [4,] 4.3 [5,] 5.5 [6,] 6.7 There will be one value in 'mp' for each bar in your series. You can then use those values along the x axis to draw your line segments under the bars as you require, based upon the cut points you want to highlight. To get the center of a given group of bars, you can use: mean(mp[start:end]) where 'start' and 'end' are the extreme bars in each of your groups. Two other things that might be helpful. See ?cut and ?hist, noting the output in the latter when 'plot = FALSE'. HTH, Marc Schwartz __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Generation of missiing values in a time serie...
Yes, this is the definition of a time series and therefore of a zoo object. A time series is a mathematical function, i.e. it assigns a single element of the range to each element of the domain. This data does not describe a time series. Also it has no underlying regularity as the warning message states. To use as.ts one wants a series with an underlying regularity that has gaps and then as.ts will fill in the gaps with NAs. If we don't have an underlying regularity the question is not well posed but its likely we want to discretize time. The zoo command itself is somewhat forgiving, at least in this case, i.e. it allows one to specify this illegal zoo object with non-unique times for purposes of discretization; however, such a zoo object should not be used other than to get a legal zoo object out. For example, in the following we round the times to one decimal place and then within each set of values at the same discretized time take the last one. (Alternately specify mean instead of tail, 1 if the average is prefered.) Then we convert that to a ts object: as.ts(aggregate(z, round(time(z), 1), tail, 1)) Time Series: Start = c(123, 2) End = c(123, 8) Frequency = 10 time flow seq ts x rtt size 123.1 123.12570 967 123.1257 13394 0.798205 1472 123.2 123.24110 969 123.2411 12680 0.796258 1472 123.3 NA NA NA NANA NA NA 123.4 NA NA NA NANA NA NA 123.5 123.47260 970 123.4726 12680 0.796258 1472 123.6 123.58860 971 123.5886 12680 0.796258 1472 123.7 123.70460 972 123.7046 12680 0.796258 1472 On 12/13/05, Alvaro Saurin [EMAIL PROTECTED] wrote: I think I have found the error. It appears when there are two entries with the same time. Using as input file: - CUT # Output format for PCKs: # TIME FLOW P [+-] SEQ TS X RTT SIZE # 123.125683 0 P + 967 123.125683 13394 0.798205 1472 123.241137 0 P + 968 123.241137 12680 0.796258 1472 123.241137 0 P + 969 123.241137 12680 0.796258 1472 123.472631 0 P + 970 123.472631 12680 0.796258 1472 123.588613 0 P + 971 123.588613 12680 0.796258 1472 123.704594 0 P + 972 123.704594 12680 0.796258 1472 - CUT I run fhe following code: - CUT h_types - list (0, 0, NULL, NULL, 0, 0, 0, 0, 0) h_names - list (time, flow, seq, ts, x, rtt, size) pcks_file- pipe (grep ' P ' data, r) pcks - scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) mat_df - data.frame (pcks[1:2], pcks[5:9]) mat - as.matrix (mat_df) colnames (mat) - h_names z - zoo (mat, mat [,time]) - CUT The dput of 'z' shows: - CUT structure(c(123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594, 0, 0, 0, 0, 0, 0, 967, 968, 969, 970, 971, 972, 123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594, 13394, 12680, 12680, 12680, 12680, 12680, 0.798205, 0.796258, 0.796258, 0.796258, 0.796258, 0.796258, 1472, 1472, 1472, 1472, 1472, 1472 ), .Dim = c(6, 7), .Dimnames = list(c(1, 2, 3, 4, 5, 6), c(time, flow, seq, ts, x, rtt, size)), index = structure(c(123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594), .Names = c(1, 2, 3, 4, 5, 6)), class = zoo) - CUT If I try a 'as.ts(z)', it fails. If I remove the duplicate entry, I can convert it to a TS with no problem. Is this made intentionally? Because then I have to filter the input matrix... But, anyway, the output matrix, after filtering, doesn't seem regular: - CUT as.ts (z) Time Series: Start = 1 End = 5 Frequency = 1 time flow seq ts x rtt size 1 123.12570 967 123.1257 13394 0.798205 1472 2 123.24110 969 123.2411 12680 0.796258 1472 3 123.47260 970 123.4726 12680 0.796258 1472 4 123.58860 971 123.5886 12680 0.796258 1472 5 123.70460 972 123.7046 12680 0.796258 1472 Warning message: 'x' does not have an underlying regularity in: as.ts.zoo(z) - CUT Weird... On 13 Dec 2005, at 16:33, Gabor Grothendieck wrote: Please provide a reproducible example. Note that dput(x) will output an R object in a way that can be copied and pasted into another session. On 12/13/05, Alvaro Saurin [EMAIL PROTECTED] wrote: On 13 Dec 2005, at 13:08, Gabor Grothendieck wrote: Your variable mat is not a matrix; its a data frame. Check it with: class(mat) Here is an example: x - cbind(A = 1:4, B = 5:8) tt - c(1, 3:4, 6) library(zoo) x.zoo - zoo(x, tt) x.ts - as.ts(x.zoo) Fixed, but anyway it fails: h_types - list (0, 0, NULL, NULL, 0, 0, 0, 0, 0) h_names - list (time, flow, seq, ts, x, rtt, size) pcks_file - pipe (grep ' P ' server.dat, r) pcks- scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) mat_df
[R] Frailty with a Weibull PH Model
Hi, I am trying to fit a Weibull PH model with a gaussian frailty term (kidney data). Is the following approach correct? 1. Fit the model using survReg with the frailty term added in. Kidn1-survreg(formula = Surv(rtime, event)~ age+sex+gn+an+pkn+frailty.gaussian(patient), data=kidney, dist='weibull') 2. Based on Venables and Ripley (350-353), divide the negative of the coefficients by the scale parameter to obtain the Weibull PH model coefficients. Jindi Jatinder Singh Senior Manager, Analysis and Reporting PRA International 300-730 View Street Victoria, B.C. V8W 3Y7 Tel: 250-483-4416 Fax: 250 483 4588 http://www.prainternational.com e-mail: [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Ploting graphics using X tints from a color
Hi, I'm trying to draw a 2D plot using multiple tints of red. The (simplified) setup is the following: || year | x | y || My idea is that each year is plotted with a different tint of red. Older year (lightest) - Later year (darkest). I've managed to plot this with different scales of grays simply by doing: palette(gray(length(years):0/length(years))) before the plot and for each year the color used is a different tint of gray. So, is there any way to do this for any color? Any tip or advice? With this, I hope to visualize patterns in my dataset more easily. Thanks in advance for any help. Best regards, Sérgio Nunes __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Incomplete Beta
On 13-Dec-05 Prof Brian Ripley wrote: On Tue, 13 Dec 2005 [EMAIL PROTECTED] wrote: On 13-Dec-05 Thomas Lumley wrote: On Tue, 13 Dec 2005, Albert Sorribas wrote: Is there any function available in R for computing the incomplete Beta function? pbeta(). The incomplete Beta function is the cdf of the Beta distribution But don't forget to multiply by beta(,): ibeta(x,a,b) - function(x,a,b){ pbeta(x,a,b)*beta(a,b) } ! Depends on which definition you use, as ?pbeta explains. Thomas' advice was correct rather than yours for Abramowitz and Stegun's definition, for example. Hmmm ... In my edition (1964, Dover repr. 1966), Section 6.6 Incomplete Beta Function: 6.6.1 B_x(a,b) = the definition I was using 6.6.2 I_x(a,b) = B_x(a,b)/B(a,b) the latter referring on to Chapter 26 Probability Functions, Section 26.5 Incomplete Beta Function which reproduces the second (6.6.2) definition. There has clearly long been ambiguity here. AS use Incomplete Beta Function in 26.5 where I (and others) would prefer Beta Distribution. They do the same sort of thing for the Incomplete Gamma Function in 6.5, where their 6.5.1 is the analogue for Gamma of 6.6.2 for Beta, and their 6.5.2 the analogue of 6.6.1. Their use of it in Chap 26 Probability Functions is in relation to the Chi-Square Probability Function (see esp. 26.4.19). But the Father (or more accurately the Midwife) of the Incomplete Beta Function was Karl Pearson, whose Introduction (1933) to the Tables of the Incomplete Beta Function states: The function I proposed to have tabled was to be a *probability integral*; that is to say, if we represent by B(p,q) the complete B-function, = Int[0,1] x^(p-1) (1-x)^(q-1) dx, and by B_x(p,q) the incomplete B-function, or Int[0,x]...dx, [= AS 6.6.1] we tabled the ratio I_x(p,q) = B_x(p,q)/B(p,q) = ... [= AS 6.6.2] and the Table of Contents lists Table I: Incomplete Beta Function Ratio (though the title page of the Table section simply calls it Incomplete Beta Function). However, on balance it seems that Pearson meant to reserve Incomplete Beta Function for the simple integral, not normalised to the Ratio. My reasons for preferring the terminology Incomplete ... Function for the incomplete integral *not* divided by the normalising constant (for both Beta and Gamma), and using Distribution for the incomplete integral divided by the constant (i.e. Pearson's Ratio), are several, but in summary: 1. The Beta and Gamma functions (not normalised) are fundamental mathematical functions in their own right; likewise their incomplete versions. 2. When needed in probability applications, then of course they need to be normalised; but then why not simply call them distributions? 3. (1) and (2) encapsulate in the terminology an essential distinction, and using (2) instead of (1) could lead to interesting inferences (e.g. that the complete Beta function is identially 1). I.e. the Beta function should not change its definition as x passes from 1 - epsilon to 1. And similarly for the Gamma. Granted there is non-uniformity of usage; but this does lead to confusion, which could be avoided by simply sticking to the distinction between Incomplete ... Function and ... Distribution. Best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 13-Dec-05 Time: 18:40:17 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Manipulating matrices
Hi, I'm pretty new to R and I've been having some problems filtering data in matrices. I have the following initial dataset: || year | name | varA || I have multiple values for varA for the same year and the same name. Having this as the input I would like to obtain the following: || year | name | {varA mean} || Where I only have one line for each year and name with the mean of the values of varA in varA mean. Is there a simple way to achieve this without using control structures (for or while cycles)? Thanks in advance for any help. Best regards, Sérgio Nunes __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Fwd: Re: Wavelet reconstruction
I realized that I may not have answered the question you were asking and that no one else has responded. I can across a similar problem and may have an answer to your question now. If you have both the wavelet coefficients and the scaling coefficients then create a fake sequence of the same length as the original and decompose that sequence using wd form wavethersh with the same wavelet family and filter that was used to decompose the data. Then you can replace the wavelet coefficients and the scaling coefficients using putC and putD from wavethresh. This will leave you with a wd object that you can then reconstruct using wr from wavethresh. I hope this works and is still useful to you. Elizabeth Lawson Elizabeth Lawson [EMAIL PROTECTED] wrote: Date: Wed, 19 Oct 2005 08:04:02 -0700 (PDT) From: Elizabeth Lawson [EMAIL PROTECTED] Subject: Re: [R] Wavelet reconstruction To: Amir Safari [EMAIL PROTECTED] Using wavethresh you can recompose a decompose signal using wr. Here is an example of decomposing, thresholding and recomposing a signal. library(wavethresh) brain-c(0,0,0.5,1,-0.75,-0.25,1.8,-3,0.41667,1.08333,-1.8, -0.58333,2.16667,-4.08333,5.75,9.75,1.58333,0.75,15.8333, 16.6667,7.7,8.16667,1.3,-3.3,-2.75,-2.08333, -1.75,0.41667,1.25,8.58333,-4.58333,0.7,-6.41667,3.58333, 3.41667,-3.3,-7.25,-1.8,-1.5,-0.08333,-2.3,7.75,5, -2.3,12,10.5,-1.3,-3.3,-3.41667,14.0833,5.16667, 5.16667,2.25,-0.08333,-1.25,-1,0.08333,1.7,1,-1.3,0.41667, -0.16667,-0.25,-0.16667) x-( c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26, 27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51, 52,53,54,55,56,57,58,59,60,61,62,63,64)) x-x/64 par(mfrow=c(1,1)) plot(x,brain,xlab=Voxel,ylab=Activity,main=fMRI Data) wdbrain-wd(brain,4,family=DaubExPhase, bc=periodic) thres2-threshold(wdbrain,levels=3:(wdbrain$nlevels-1), type=soft, policy=manual, by.level=FALSE, value=7.32032, dev=var, boundary=FALSE, verbose = getOption(verbose), return.threshold=F) thr2 - wr(thres2) plot(x,brain, col = slateblue,xlab=Voxel,ylab=Activity,main=Wavelet Regression) mtext(N=4, Threshold=7.32032) lines(x, thr2, col= violetred , lwd=2,type=l) Good Luck!! Elizabeth Lawson Amir Safari [EMAIL PROTECTED] wrote: Hi There, I tried to find a function in {waveslim} or {wavethresh} in order to reconstruct the decomposed signals. As far as I found there is no function in {waveslim} to reconstruct decomposed data. The function wr{wavethresh} reconstructs the results of wd function. Apart from its limitations ( for example the length of vector must be power of 2 ) it apparently doesn't work with the functions and objects in waveslim. What could help ? So many thanks for your any idea. Amir Safari - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html - - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] creating a subset of a dataset using ifelse statement?
?ifelse ; it's clear that you can use it in your case. Where cluster_q is comming from ( I suppose is the column name you want to filter by!)? What if you use the simpler syntax : if(h1) mp2-mp1 else mp2-subset(mp1b,cluster_q==h) ? THT t c wrote: I am using R in a Microsoft Windows environment. I have a dataset called mp1b. I have a variable called h. h can take a value from -1 to 5. If h 1, I want to create a new dataset called mp2 that is the same as mp1b: mp2-mp1b If h 0, I want to set create a dataset mp2, where I limit the original dataset to those where mp1b$group = =h. similar to: mp2-subset (mp1b, group= = h) I have tried this ifelse statement, but it does not seem to work as expected. mp2-ifelse(h1,mp1b,subset(mp1b,cluster_q==h)) Assistance is appreciated. - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Ferdinand Alimadhi Programmer / Analyst Harvard University The Institute for Quantitative Social Science (617) 496-0187 [EMAIL PROTECTED] www.iq.harvard.edu [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] correct C function usage
Hello, I am not sure if I am interfacing with C correctly and _safely_ or if there is a better way esp. with regards to terminating the returned array. I am trying to fill an int array with values whose actual size is determined in the C function and is always maximally as large as length(values). I also don't understand the purpose of $ab in the example: conv - function(a, b) .C(convolve, -snip- ab = double(length(a) + length(b) - 1))$ab void testFill(int *values, int *newvalues, int* endposition ){ newvalues[0] = 1; newvalues[1] = 2; *endposition = 2; } dyn.load(../testFill.so) testTestFill - function(){ tempfilled - testFillC( c(30:40)) realfilled - tempfilled$newvalues[1:tempfilled$endposition] return(realfilled) } testFillC - function(a){ .C(testFill, as.integer(a), newvalues=integer(length(a)), endposition=integer(1)) } Thank you very much in advance Ido Tamir __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Generation of missiing values in a time serie...
In thinking about this some more, the trick I discussed is probably not the best way to do it since its possible that in the future zoo will completely disallow illegal zoo objects. I think a better way might be to construct it like this: aggregate(zoo(z.data), round(z.time, 1), tail, 1) where z.data is the matrix and z.time are the times. The variable z, which is an illegal zoo object, would not be created but in terms of z, since that is what I have reproducibly from your post, we have: z.data - coredata(z) z.time - time(z) On 12/13/05, Gabor Grothendieck [EMAIL PROTECTED] wrote: Yes, this is the definition of a time series and therefore of a zoo object. A time series is a mathematical function, i.e. it assigns a single element of the range to each element of the domain. This data does not describe a time series. Also it has no underlying regularity as the warning message states. To use as.ts one wants a series with an underlying regularity that has gaps and then as.ts will fill in the gaps with NAs. If we don't have an underlying regularity the question is not well posed but its likely we want to discretize time. The zoo command itself is somewhat forgiving, at least in this case, i.e. it allows one to specify this illegal zoo object with non-unique times for purposes of discretization; however, such a zoo object should not be used other than to get a legal zoo object out. For example, in the following we round the times to one decimal place and then within each set of values at the same discretized time take the last one. (Alternately specify mean instead of tail, 1 if the average is prefered.) Then we convert that to a ts object: as.ts(aggregate(z, round(time(z), 1), tail, 1)) Time Series: Start = c(123, 2) End = c(123, 8) Frequency = 10 time flow seq ts x rtt size 123.1 123.12570 967 123.1257 13394 0.798205 1472 123.2 123.24110 969 123.2411 12680 0.796258 1472 123.3 NA NA NA NANA NA NA 123.4 NA NA NA NANA NA NA 123.5 123.47260 970 123.4726 12680 0.796258 1472 123.6 123.58860 971 123.5886 12680 0.796258 1472 123.7 123.70460 972 123.7046 12680 0.796258 1472 On 12/13/05, Alvaro Saurin [EMAIL PROTECTED] wrote: I think I have found the error. It appears when there are two entries with the same time. Using as input file: - CUT # Output format for PCKs: # TIME FLOW P [+-] SEQ TS X RTT SIZE # 123.125683 0 P + 967 123.125683 13394 0.798205 1472 123.241137 0 P + 968 123.241137 12680 0.796258 1472 123.241137 0 P + 969 123.241137 12680 0.796258 1472 123.472631 0 P + 970 123.472631 12680 0.796258 1472 123.588613 0 P + 971 123.588613 12680 0.796258 1472 123.704594 0 P + 972 123.704594 12680 0.796258 1472 - CUT I run fhe following code: - CUT h_types - list (0, 0, NULL, NULL, 0, 0, 0, 0, 0) h_names - list (time, flow, seq, ts, x, rtt, size) pcks_file- pipe (grep ' P ' data, r) pcks - scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) mat_df - data.frame (pcks[1:2], pcks[5:9]) mat - as.matrix (mat_df) colnames (mat) - h_names z - zoo (mat, mat [,time]) - CUT The dput of 'z' shows: - CUT structure(c(123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594, 0, 0, 0, 0, 0, 0, 967, 968, 969, 970, 971, 972, 123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594, 13394, 12680, 12680, 12680, 12680, 12680, 0.798205, 0.796258, 0.796258, 0.796258, 0.796258, 0.796258, 1472, 1472, 1472, 1472, 1472, 1472 ), .Dim = c(6, 7), .Dimnames = list(c(1, 2, 3, 4, 5, 6), c(time, flow, seq, ts, x, rtt, size)), index = structure(c(123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594), .Names = c(1, 2, 3, 4, 5, 6)), class = zoo) - CUT If I try a 'as.ts(z)', it fails. If I remove the duplicate entry, I can convert it to a TS with no problem. Is this made intentionally? Because then I have to filter the input matrix... But, anyway, the output matrix, after filtering, doesn't seem regular: - CUT as.ts (z) Time Series: Start = 1 End = 5 Frequency = 1 time flow seq ts x rtt size 1 123.12570 967 123.1257 13394 0.798205 1472 2 123.24110 969 123.2411 12680 0.796258 1472 3 123.47260 970 123.4726 12680 0.796258 1472 4 123.58860 971 123.5886 12680 0.796258 1472 5 123.70460 972 123.7046 12680 0.796258 1472 Warning message: 'x' does not have an underlying regularity in: as.ts.zoo(z) - CUT Weird... On 13 Dec 2005, at 16:33, Gabor Grothendieck wrote: Please provide a reproducible example. Note that dput(x) will output an R object in a way that can be copied and pasted into another session.
[R] what does this warnings mean? and what should I do?
I use lmer to fit a mixed effect model.It give some warnings.what does this warnings mean? and what should I do? (fm2.mlm - lmer(qd ~ edu + jiankang + peixun +hunyin + cadcj + age + age2 + sex + dangyuan + Comp.1 + Comp.2+trust.cz1 +(trust.cz1|commid), data = individual,na.action = na.exclude,family=quasibinomial)) Generalized linear mixed model fit using PQL Formula: qd ~ edu + jiankang + peixun + hunyin + cadcj + age + age2 + sex + dangyuan + Comp.1 + Comp.2 + trust.cz1 + (trust.cz1 | commid) Data: individual Family: quasibinomial(logit link) AIC BIClogLik deviance 736.7059 821.8267 -349.3529 698.7059 Random effects: Groups NameVariance Std.Dev. Corr commid (Intercept) 1.56413 1.25065 trust.cz1 0.17922 0.42334 1.000 Residual 0.89728 0.94725 # of obs: 652, groups: commid, 39 Fixed effects: Estimate Std. Error DF t value Pr(|t|) (Intercept)-1.6115e-01 6.7997e-01 637 -0.2370 0.81274 edu-5.2585e-02 4.1048e-02 637 -1.2810 0.20064 jiankang -9.8243e-01 4.4645e-01 637 -2.2005 0.02813 * peixun -4.6307e-01 2.6397e-01 637 -1.7542 0.07988 . hunyin -1.2255e-02 2.8151e-01 637 -0.0435 0.96529 hunyin -2.7726e-01 1.3846e+00 637 -0.2002 0.84136 hunyin -2.9759e-01 8.7180e-01 637 -0.3414 0.73295 cadcj 2.2366e-01 7.6467e-01 637 0.2925 0.77000 age 9.3626e-02 4.0390e-02 637 2.3180 0.02076 * age2 -1.3095e-03 5.5104e-04 637 -2.3763 0.01778 * sex 3.9188e-01 1.9759e-01 637 1.9833 0.04776 * dangyuan -5.2558e-01 5.9091e-01 637 -0.8894 0.37410 Comp.1 5.2463e-02 1.0309e-01 637 0.5089 0.61100 Comp.2 -1.5048e-01 1.1435e-01 637 -1.3160 0.18863 trust.cz1 -8.0709e-01 4.4632e-01 637 -1.8083 0.07103 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 There were 11 warnings (use warnings() to see them) warnings() Warning messages: 1: optim or nlminb returned message false convergence (8) in: LMEopt(x = mer, value = cv) 2: optim or nlminb returned message false convergence (8) in: LMEopt(x = mer, value = cv) 3: optim or nlminb returned message false convergence (8) in: LMEopt(x = mer, value = cv) 4: optim or nlminb returned message false convergence (8) in: LMEopt(x = mer, value = cv) 5: optim or nlminb returned message false convergence (8) in: LMEopt(x = mer, value = cv) 6: optim or nlminb returned message false convergence (8) in: LMEopt(x = mer, value = cv) 7: optim or nlminb returned message false convergence (8) in: LMEopt(x = mer, value = cv) 8: optim or nlminb returned message false convergence (8) in: LMEopt(x = mer, value = cv) 9: optim or nlminb returned message singular convergence (7) in: LMEopt(x = mer, value = cv) 10: optim or nlminb returned message singular convergence (7) in: LMEopt(x = mer, value = cv) 11: optim or nlminb returned message singular convergence (7) in: LMEopt(x = mer, value = cv) version _ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major2 minor2.0 year 2005 month10 day 06 svn rev 35749 language R 2005-12-14 -- Deparment of Sociology Fudan University My new mail addres is [EMAIL PROTECTED] Blog:http://sociology.yculblog.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Ploting graphics using X tints from a color
Le 13.12.2005 19:34, Sérgio Nunes a écrit : Hi, I'm trying to draw a 2D plot using multiple tints of red. The (simplified) setup is the following: || year | x | y || My idea is that each year is plotted with a different tint of red. Older year (lightest) - Later year (darkest). I've managed to plot this with different scales of grays simply by doing: palette(gray(length(years):0/length(years))) before the plot and for each year the color used is a different tint of gray. So, is there any way to do this for any color? Any tip or advice? With this, I hope to visualize patterns in my dataset more easily. Thanks in advance for any help. Best regards, Sérgio Nunes Hi, You want to travel in the RGB space. See http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=107 plot(0, xlim=c(0,10), ylim=c(0,1)) pal - rgb(1, 0:10/10,0:10/10) rect(xleft=0:9, xright=1:10, ytop=1, ybottom=0, col=pal) Romain -- visit the R Graph Gallery : http://addictedtor.free.fr/graphiques mixmod 1.7 is released : http://www-math.univ-fcomte.fr/mixmod/index.php +---+ | Romain FRANCOIS - http://francoisromain.free.fr | | Doctorant INRIA Futurs / EDF | +---+ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Ploting graphics using X tints from a color
I'm trying to draw a 2D plot using multiple tints of red. The (simplified) setup is the following: || year | x | y || You might find this function useful: map_colour_gradient - function(x, low=red, mid=white,high=black, midpoint = 0) { x - as.numeric(x) low.rgb - col2rgb(low, TRUE)/256 mid.rgb - col2rgb(mid, TRUE)/256 high.rgb - col2rgb(high, TRUE)/256 interp_r - approxfun(c(min(x, na.rm=TRUE), midpoint, max(x, na.rm=TRUE)), c(low.rgb[1], mid.rgb[1], high.rgb[1])) interp_g - approxfun(c(min(x, na.rm=TRUE), midpoint, max(x, na.rm=TRUE)), c(low.rgb[2], mid.rgb[2], high.rgb[2])) interp_b - approxfun(c(min(x, na.rm=TRUE), midpoint, max(x, na.rm=TRUE)), c(low.rgb[3], mid.rgb[3], high.rgb[3])) rgb(interp_r(x), interp_g(x), interp_b(x)) } Given a numeric vector x, it will create a smooth gradient (linearly through RGB space). eg. x - rnorm(100) y - rnorm(100) plot(x, y, col=map_colour_gradient(x), pch=20) plot(x, y, col=map_colour_gradient(x, low=black, high=black, mid=yellow), pch=20) Note, however, using colour is only likely to find the most prominent of patterns. Hadley __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Age of an object?
It would be nice to have a date stamp on an object. In S/Splus this was always available, because objects were files. I have looked around, but I presume this information is not available. Trevor Hastie [EMAIL PROTECTED] Professor, Department of Statistics, Stanford University Phone: (650) 725-2231 (Statistics) Fax: (650) 725-8977 (650) 498-5233 (Biostatistics) Fax: (650) 725-6951 URL: http://www-stat.stanford.edu/~hastie address: room 104, Department of Statistics, Sequoia Hall 390 Serra Mall, Stanford University, CA 94305-4065 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Ploting graphics using X tints from a color
Le 13.12.2005 21:40, Romain Francois a écrit : Le 13.12.2005 19:34, Sérgio Nunes a écrit : Hi, I'm trying to draw a 2D plot using multiple tints of red. The (simplified) setup is the following: || year | x | y || My idea is that each year is plotted with a different tint of red. Older year (lightest) - Later year (darkest). I've managed to plot this with different scales of grays simply by doing: palette(gray(length(years):0/length(years))) before the plot and for each year the color used is a different tint of gray. So, is there any way to do this for any color? Any tip or advice? With this, I hope to visualize patterns in my dataset more easily. Thanks in advance for any help. Best regards, Sérgio Nunes Hi, You want to travel in the RGB space. See http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=107 plot(0, xlim=c(0,10), ylim=c(0,1)) pal - rgb(1, 0:10/10,0:10/10) rect(xleft=0:9, xright=1:10, ytop=1, ybottom=0, col=pal) Romain Ops ! s - seq(0,1, length=10) pal - rgb(1, s, s) is better. Sorry. And more generally, see ?colorRampPalette. f - colorRampPalette(c(royalblue, white)) pal - f(10) Romain -- visit the R Graph Gallery : http://addictedtor.free.fr/graphiques mixmod 1.7 is released : http://www-math.univ-fcomte.fr/mixmod/index.php +---+ | Romain FRANCOIS - http://francoisromain.free.fr | | Doctorant INRIA Futurs / EDF | +---+ __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] creating a subset of a dataset using ifelse statement?
t c wrote: I am using R in a Microsoft Windows environment. I have a dataset called “mp1b”. I have a variable called h. h can take a value from -1 to 5. If h 1, I want to create a new dataset called mp2 that is the same as mp1b: “mp2-mp1b” If h 0, I want to set create a dataset mp2, where I limit the original dataset to those where mp1b$group = =h. similar to: “mp2-subset (mp1b, group= = h)” I have tried this ifelse statement, but it does not seem to work as expected. “mp2-ifelse(h1,mp1b,subset(mp1b,cluster_q==h))” mp2 - if(h1) mp1b else subset(mp1b,cluster_q==h) Uwe Ligges Assistance is appreciated. - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] batch mode problem
Jörg Schaber wrote: Hi, ok, I know I should be using a later version than 1.7.1 (64 bit) but it's not in my power. So here is the problem: In my R script I declare a data.frame that consists of 40 vectors, each having 125 numeric elements. This is no problem as long as I run the sript in interactive mode, but running it in batch mode I get strange error messages. Apparently, it has to do with the size of the data.frame because reducing the data.frame to 36 vector a 125 elements, I have no problems. How can I declare my large data.frame and still run the script in batch mode? Put the data.frame in a seperate file and source that one. Uwe Ligges Thanks, joerg __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Ploting graphics using X tints from a color
Sérgio Nunes snunes at gmail.com writes: Hi, I'm trying to draw a 2D plot using multiple tints of red. The (simplified) setup is the following: || year | x | y || My idea is that each year is plotted with a different tint of red. Older year (lightest) - Later year (darkest). how about: x = runif(300) y = runif(300) year = factor(rep(1:30,each=10)) nyear = length(levels(year)) grays = gray((nyear:0)/nyear) reds = hsv(h=1,s=(0:nyear)/nyear,v=1) plot(x,y,col=grays[as.numeric(year)],pch=16) plot(x,y,col=reds[as.numeric(year)],pch=16) plot(x,y,bg=reds[as.numeric(year)],pch=21,fg=1) legend(0.8,1,levels(year),pch=21,pt.bg=reds,col=1,bg=white,ncol=2) __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] correct C function usage
On Tue, 13 Dec 2005, Ido M. Tamir wrote: Hello, I am not sure if I am interfacing with C correctly and _safely_ or if there is a better way esp. with regards to terminating the returned array. You need to pass the length to the C routine and check you do not overwrite it. (As in the parts you -snip-ed below.) I am trying to fill an int array with values whose actual size is determined in the C function and is always maximally as large as length(values). I also don't understand the purpose of $ab in the example: conv - function(a, b) .C(convolve, -snip- ab = double(length(a) + length(b) - 1))$ab .C returns a list, an element for each argument after the first, named if the arguments were named. So this selects the (copy) of the vector sent back as the last argument of the C function. void testFill(int *values, int *newvalues, int* endposition ){ newvalues[0] = 1; newvalues[1] = 2; *endposition = 2; } dyn.load(../testFill.so) testTestFill - function(){ tempfilled - testFillC( c(30:40)) realfilled - tempfilled$newvalues[1:tempfilled$endposition] return(realfilled) } testFillC - function(a){ .C(testFill, as.integer(a), newvalues=integer(length(a)), endposition=integer(1)) } What do testFillC(1) or testFillC(logical(0)) do? -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] export from R to MySQL
On Dec 12, 2005, at 6:50 PM, David James wrote: Prof Brian Ripley wrote: On Mon, 12 Dec 2005, Sean Davis wrote: On 12/12/05 9:21 AM, bogdan romocea [EMAIL PROTECTED] wrote: Sean Davis wrote: but you will have to create the table by hand There's no need for manual steps. To take advantage of MySQL's extremely fast 'load data infile' you could dump the data in CSV format, write a script for mysql (the command line tool), for example q - function(table,infile) { query - paste( create table ,table, (col1 float, col2 float); This is creating the table by hand, as opposed to using dbWriteTable. If your data.frame contains 67 columns, using dbWriteTable saves quite a bit of typing The RODBC equivalent creates the table for you, then fast imports the file. Might be worthwhile contribution to RMySQL for someone. That's what RMySQL's dbWriteTable() does. The original posting mentioned problems associated with speed of data.frame and dbWriteTable, which seems plausible (but I haven't quantified it myself) given the fact that dbWriteTable outputs a data.frame to an intermediate file via write.table and then uses the LOAD DATA for fast loading that intermediate file. Thanks at all! As write.table and read.table itself are to some degree slow, for matrizes which are only numeric cat() and scan() could be faster. however it's a special case. Just be careful with client-server systems to have the file in the right place (if indeed you are allowed to have files on the server). -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- David __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] fSeries
Braesch Lucas wrote: I'm trying to use garchFit from fSeries, with Student or Skewed Student conditionnal distribution. Let's say that eps (vector) is my series of daily log-returns: data(EuStockMarkets) eps = diff(log(EuStockMarkets[,CAC])) library(fSeries) g = garchFit(series = eps, formula.var = ~garch(2,2), cond.dist = dstd) s = [EMAIL PROTECTED] All the coefficients are ok (checked with SAS 9.1) except nu (degrees of freedom of the student) and the log-likelyhood. I've really checked everything and can't find the estimated series sigma (volatility) and eta, such that eps = sigma * eta and eta is centered and reduced... I've tryed combinations of all s$x,s$h,s$z and nothing looks looks correct. Also, is it possible to fit EGARCH and TGARCH with R ? If anyone ever managed to make it work, i'd be grateful ;-) Do you think, that SAS is right? - Please can you post the results from SAS? This is a good example which shows what can go wrong in GARCH Modelling!!! First simulate with Rmetrics: data(EuStockMarkets) eps = as.vector(diff(log(EuStockMarkets[,CAC]))) var(x) # Important - Maybe you have a scale problem in optimization because # your variance paramater is so small compared with the other parameters! # Thus, multiply with 100: x = 100* as.vector(eps) # Rmetrics: garchFit(formula.mean = ~arma(0,0), formula.var = ~garch(2,2), cond.dist = dstd) #mu omega alpha1 alpha2 beta1 beta2 shape # 0.0523284 0.0421556 0.0455789 0.010 0.8678519 0.0523520 7.9870453 # Now I simulated with Ox and S-Plus, in both cases I found convergence problems. # The reason may be that your model is not a GARCH(2,2) it's a GARCH(1,2) model! # Now Try: garchFit(formula.mean = ~arma(0,0), formula.var = ~garch(1,2), cond.dist = dstd) #mu omega alpha1 beta1 beta2 shape # 0.0523284 0.0421547 0.0455790 0.8678688 0.0523368 7.9870458 # Great, we get the same result! # Now, try Ox/[EMAIL PROTECTED], the result is: Coefficient Std.Error t-value t-prob Cst(M) 0.052328 0.0237722.201 0.0278 Cst(V) 0.042139 0.0275971.527 0.1269 ARCH(Alpha1) 0.045604 0.0253771.797 0.0725 GARCH(Beta1) 0.8676640.648081.339 0.1808 GARCH(Beta2) 0.0525550.60865 0.08635 0.9312 Student(DF) 7.983069 1.15536.910 0. # Now try S-Plus/Finmetrics, the result is: Conditional Distribution: t with estimated parameter 7.937377 and standard error 1.148712 Value Std.Error t value Pr(|t|) C 0.053110.02377 2.2344 0.01279 A 0.043550.02818 1.5455 0.06120 ARCH(1) 0.046530.02553 1.8230 0.03423 GARCH(1) 0.855120.64209 1.3318 0.09155 GARCH(2) 0.063030.60239 0.1046 0.45834 # So Rmetrics, Ox, and S-Plus are in agreement!!! # What is with SAS? Please give us the results for GARCH(1,2) # and GARCH(2,2)! # Please note, garchFit() from Rmetrics is still in # testing phase. An updated version is just under preparation. Diethelm Wuertz __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Manipulating matrices
aggregate(x$varA, by=list(x$year, x$name), mean) On 12/13/05, Sérgio Nunes [EMAIL PROTECTED] wrote: Hi, I'm pretty new to R and I've been having some problems filtering data in matrices. I have the following initial dataset: || year | name | varA || I have multiple values for varA for the same year and the same name. Having this as the input I would like to obtain the following: || year | name | {varA mean} || Where I only have one line for each year and name with the mean of the values of varA in varA mean. Is there a simple way to achieve this without using control structures (for or while cycles)? Thanks in advance for any help. Best regards, Sérgio Nunes __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Jim Holtman Cincinnati, OH +1 513 247 0281 What the problem you are trying to solve? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] getting faster results
Hey, Can anyone answer this question. I am working with really large datasets and most of the programs I have been running take quite some time. I heard that R may be faster in Unix. I sthis true and if so can anyone reccomend which system and requirements may allow things to go faster for? Thanks!! Elizabeth Lawson - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Polytopic Vector Analysis (PVA)
Melanie Edwards wrote: I am curious whether anybody has or is developing this capability within R. I have not found any software yet that has this capability and I am not sure whether it is too new a method and nobody is actually using it, googling around I found that this (PVA) is a variant of factor analysis restricted to find non-negative factors. It is not a new method, although maybe the name is. This has been/is used for instance in air quality monitoring to identify sources of pollution, and if you have some prior information about possible sources that maybe could be used to. I guess this could be called 'source unmixing' or something similar, which indicatyes a similarity with independent component analysis. Maybe enough to restrict optimization to non-negative values? help.search(factor analysis) shows that factor analysis is multiply implemented for R, so maybe there is something, and if not, maybe simple to adapt. Kjetil Halvorsen or if there are other means to get the same analysis that I do not know of. Any information regarding developments or use of this method would be helpful. Melanie Edwards [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Generation of missiing values in a time serie...
Gabor Grothendieck wrote: Yes, this is the definition of a time series and therefore of a zoo object. A time series is a mathematical function, i.e. it assigns a single element of the range to each element of the domain. This data does not describe a time series. Since nobody else has mentiones it on this thread: Tha CRAN package pastecs has function `regul' to regularize irregular time series. maybe that is what the original poster want. Kjetil Also it has no underlying regularity as the warning message states. To use as.ts one wants a series with an underlying regularity that has gaps and then as.ts will fill in the gaps with NAs. If we don't have an underlying regularity the question is not well posed but its likely we want to discretize time. The zoo command itself is somewhat forgiving, at least in this case, i.e. it allows one to specify this illegal zoo object with non-unique times for purposes of discretization; however, such a zoo object should not be used other than to get a legal zoo object out. For example, in the following we round the times to one decimal place and then within each set of values at the same discretized time take the last one. (Alternately specify mean instead of tail, 1 if the average is prefered.) Then we convert that to a ts object: as.ts(aggregate(z, round(time(z), 1), tail, 1)) Time Series: Start = c(123, 2) End = c(123, 8) Frequency = 10 time flow seq ts x rtt size 123.1 123.12570 967 123.1257 13394 0.798205 1472 123.2 123.24110 969 123.2411 12680 0.796258 1472 123.3 NA NA NA NANA NA NA 123.4 NA NA NA NANA NA NA 123.5 123.47260 970 123.4726 12680 0.796258 1472 123.6 123.58860 971 123.5886 12680 0.796258 1472 123.7 123.70460 972 123.7046 12680 0.796258 1472 On 12/13/05, Alvaro Saurin [EMAIL PROTECTED] wrote: I think I have found the error. It appears when there are two entries with the same time. Using as input file: - CUT # Output format for PCKs: # TIME FLOW P [+-] SEQ TS X RTT SIZE # 123.125683 0 P + 967 123.125683 13394 0.798205 1472 123.241137 0 P + 968 123.241137 12680 0.796258 1472 123.241137 0 P + 969 123.241137 12680 0.796258 1472 123.472631 0 P + 970 123.472631 12680 0.796258 1472 123.588613 0 P + 971 123.588613 12680 0.796258 1472 123.704594 0 P + 972 123.704594 12680 0.796258 1472 - CUT I run fhe following code: - CUT h_types - list (0, 0, NULL, NULL, 0, 0, 0, 0, 0) h_names - list (time, flow, seq, ts, x, rtt, size) pcks_file- pipe (grep ' P ' data, r) pcks - scan (pcks_file, what = h_types, comment.char = '#', fill = TRUE) mat_df - data.frame (pcks[1:2], pcks[5:9]) mat - as.matrix (mat_df) colnames (mat) - h_names z - zoo (mat, mat [,time]) - CUT The dput of 'z' shows: - CUT structure(c(123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594, 0, 0, 0, 0, 0, 0, 967, 968, 969, 970, 971, 972, 123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594, 13394, 12680, 12680, 12680, 12680, 12680, 0.798205, 0.796258, 0.796258, 0.796258, 0.796258, 0.796258, 1472, 1472, 1472, 1472, 1472, 1472 ), .Dim = c(6, 7), .Dimnames = list(c(1, 2, 3, 4, 5, 6), c(time, flow, seq, ts, x, rtt, size)), index = structure(c(123.125683, 123.241137, 123.241137, 123.472631, 123.588613, 123.704594), .Names = c(1, 2, 3, 4, 5, 6)), class = zoo) - CUT If I try a 'as.ts(z)', it fails. If I remove the duplicate entry, I can convert it to a TS with no problem. Is this made intentionally? Because then I have to filter the input matrix... But, anyway, the output matrix, after filtering, doesn't seem regular: - CUT as.ts (z) Time Series: Start = 1 End = 5 Frequency = 1 time flow seq ts x rtt size 1 123.12570 967 123.1257 13394 0.798205 1472 2 123.24110 969 123.2411 12680 0.796258 1472 3 123.47260 970 123.4726 12680 0.796258 1472 4 123.58860 971 123.5886 12680 0.796258 1472 5 123.70460 972 123.7046 12680 0.796258 1472 Warning message: 'x' does not have an underlying regularity in: as.ts.zoo(z) - CUT Weird... On 13 Dec 2005, at 16:33, Gabor Grothendieck wrote: Please provide a reproducible example. Note that dput(x) will output an R object in a way that can be copied and pasted into another session. On 12/13/05, Alvaro Saurin [EMAIL PROTECTED] wrote: On 13 Dec 2005, at 13:08, Gabor Grothendieck wrote: Your variable mat is not a matrix; its a data frame. Check it with: class(mat) Here is an example: x - cbind(A = 1:4, B = 5:8) tt - c(1, 3:4, 6) library(zoo) x.zoo - zoo(x, tt) x.ts - as.ts(x.zoo) Fixed, but anyway it fails: h_types - list (0, 0, NULL, NULL, 0, 0, 0, 0, 0) h_names - list (time, flow, seq, ts, x,
Re: [R] getting faster results
Please read the posting guide and repost. Your question is too vague to respond to sensibly. Speed in R depends on two things: How you program and what sort of system resources R has available (especially memory). OS's can certainly make a difference, but this is a 2nd order effect, I think. -- Bert Gunter Genentech Non-Clinical Statistics South San Francisco, CA The business of the statistician is to catalyze the scientific learning process. - George E. P. Box -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Elizabeth Lawson Sent: Tuesday, December 13, 2005 2:23 PM To: r-help@stat.math.ethz.ch Subject: [R] getting faster results Hey, Can anyone answer this question. I am working with really large datasets and most of the programs I have been running take quite some time. I heard that R may be faster in Unix. I sthis true and if so can anyone reccomend which system and requirements may allow things to go faster for? Thanks!! Elizabeth Lawson - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Simplifying matrices ?
Hi again, the previous answers were great. I was able to do what was planned. Now, I would like to do the following to a matrix: || year | event | team | total || where I can have multiple event per team, but each team only has a year and a total. Thus, this table has multiple lines for the same team where only the event changes. Considering this, how can I output this: || year | team | total || where each team occurs only once, and the event was discarded. Hope I made myself clear. Thanks in advance. This is a great mailing list. Regards, Sérgio Nunes __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] bug in geoR (?)
Dear Fabio PLease download the latest geoR from: www.est.ufpr.br/geoR and let me know in case the problem persists P.J. On Tue, 13 Dec 2005, Antonio, Fabio Di Narzo wrote: Date: Tue, 13 Dec 2005 17:18:14 +0100 From: Antonio, Fabio Di Narzo [EMAIL PROTECTED] To: [EMAIL PROTECTED] Cc: R-help@stat.math.ethz.ch Subject: [R] bug in geoR (?) I've enconuntered this problem with the last cran version of geoR: library(geoR) day - rep(1:2, each=5) coords - matrix(rep(runif(10),2), 10, 2) data - rnorm(10) data[1] - NA as.geodata(cbind(coords, data, day), realisations=4) as.geodata: 1 points removed due to NA in the data Errore in as.geodata(cbind(coords, data, day), realisations = 4) : realisations and coords have incompatible dimensions The problem disappear if I remove the NA manually from the dataset before passing to as.geodata. I.e.: as.geodata(cbind(coords, data, day)[2:10,], realisations=4) works. Antonio, Fabio Di Narzo. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Paulo Justiniano Ribeiro Jr LEG (Laboratório de Estatística e Geoinformação) Departamento de Estatística Universidade Federal do Paraná Caixa Postal 19.081 CEP 81.531-990 Curitiba, PR - Brasil Tel: (+55) 41 3361 3573 Fax: (+55) 41 3361 3141 e-mail: [EMAIL PROTECTED] http://www.est.ufpr.br/~paulojus__ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Sorting vectors according to a desired correlation coef
Dear all First time posting in this list. This is my question: I have three vectors (x, y, z) containing random numbers from three different distributions (normal, random, exponential). I want to sort the vectors so the pairwise correlation coefficient between them is 0.5. How can I do this. I can't use mvrnorm because I'm using different distributions. Any help will be greatly appreciated. Carlos [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lme intervals
Have you tried ?intervals and ?intervals.lme? The following is an example under ?intervals.lme: fm1 - lme(distance ~ age, Orthodont, random = ~ age | Subject) intervals(fm1) Also, have you checked Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer)? If this does not meet your needs and you don't find the answer in the documentation for predict.lme, please provide self-contained, toy example of what you want, as suggested in the posting guide! www.R-project.org/posting-guide.html. hope this helps. spencer graves Michael Kubovy wrote: Hi Dieter, No, because I'm looking for the CIs on the means of baLO of an additive model which has 20 cells, exactly as stated. Essentially I want to have the values of 'lower' and 'upper' to plug into xYplot when used in the form xYplot(Cbind(baLO,lower,upper) ~ bar | sub, groups = delta, data=e7). Thanks. On Dec 12, 2005, at 4:53 AM, Dieter Menne wrote: Michael Kubovy kubovy at virginia.edu writes: I run e7.lmeb3 - lme(baLO ~ bar + factor( delta), data = e7, random = ~ 1 | sub, method = ML) ... cut how can I get the CIs for the fixed effects in the 20 cells of the bar * delta design? A typo, maybe? Your design is bar+factor(delta), but you want bar*delta? Dieter __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Simplifying matrices ?
unique(x[,-2]) Sérgio Nunes a écrit : Hi again, the previous answers were great. I was able to do what was planned. Now, I would like to do the following to a matrix: || year | event | team | total || where I can have multiple event per team, but each team only has a year and a total. Thus, this table has multiple lines for the same team where only the event changes. Considering this, how can I output this: || year | team | total || where each team occurs only once, and the event was discarded. Hope I made myself clear. Thanks in advance. This is a great mailing list. Regards, Sérgio Nunes __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Glitch when creating online help
___ Hi, I'm writing up the online help for a package I'm developing (in-house only, sorry), and I've come across an odd glitch when trying to nest a list inside the arguments section of the .Rd file. I was just wondering if anyone could provide some insights. I'm using R 2.2.0 on Windows XP, along with ActivePerl 5.8.7 (build 815), MikTeX 2.4, and the tools downloaded from http://www.murdoch-sutherland.com/Rtools/ . Here is some code to reproduce the glitch. First, in R: f - function(x) x package.skeleton(foo, list=f) This creates the package skeleton, with a template f.Rd provided. Edit f.Rd to contain \name{f} \alias{f} \title{ ~~function to do ... ~~ } \description{ ~~ A concise (1-5 lines) description of what the function does. ~~ } \usage{f(x)} \arguments{ \item{item1}{ This is item 1. } \item{itemlist}{ Here is a list. \describe{ \item{subitem1}{Item 1 of the list.} \item{subitem2}{Item 2 of the list.} } } \item{item3}{ This is the item after the list. } } Then at the command prompt: R CMD INSTALL --build foo Once the package has been created, in R type: library(foo) ?f The result looks like fpackage:fooR Documentation ~~function to do ... ~~ Description: ~~ A concise (1-5 lines) description of what the function does. ~~ Usage: f(x) Arguments: item1: This is item 1. itemlist: Here is a list. .in +5 subitem1 Item 1 of the list. subitem2 Item 2 of the list. item3: This is the item after the list. Note the .in +5 at the top of the nested list. This is only in the online help within R, not the html version. -- Hong Ooi Senior Research Analyst, IAG Limited 388 George St, Sydney NSW 2000 +61 (2) 9292 1566 ___ The information transmitted in this message and its attachme...{{dropped}} __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] broken package?
I've just tried to install the tseries package (which seems to include quite a bit) and got the following results: package 'dynlm' successfully unpacked and MD5 sums checked package 'leaps' successfully unpacked and MD5 sums checked package 'chron' successfully unpacked and MD5 sums checked package 'fCalendar' successfully unpacked and MD5 sums checked package 'strucchange' successfully unpacked and MD5 sums checked package 'DAAG' successfully unpacked and MD5 sums checked package 'quadprog' successfully unpacked and MD5 sums checked Error in sprintf(gettext(unable to move temp installation '%d' to '%s'), : use format %s for character objects I'm running R-2.1.0 Is this updated in 2.2? --- Jeff D. Hamann Forest Informatics, Inc. PO Box 1421 Corvallis, Oregon USA 97339-1421 541-754-1428 [EMAIL PROTECTED] www.forestinformatics.com __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] About help on 'mahalanobis'
Hi, help on 'mahalanobis' (in the stats package in Rv2.2.0) now says: Description: Returns the Mahalanobis distance of all rows in 'x' and the vector mu='center' with respect to Sigma='cov'. This is (for vector 'x') defined as D^2 = (x - mu)' Sigma^{-1} (x - mu) It does return D^2 as written. However, would the text be more clear if it says Returns the _squared_ Mahalanobis distance D^2... instead? If so, then text in the example code, e.g. ##- Here, D^2 = usual Euclidean distances and the title of the first plot will also have to be updated. Compare this with what dist() in the same package returns. When asking for the Equlidean distance (matrix) between rows in a matrix, we get D not D^2, e.g. dist(c(1,3)) == 2. Cheers Henrik __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Memory shortage running Repeated Measures (nlme)
Dear group, I tried to run a Repeated Mesures Anova for Mixed effects model and I got a warnning after entering the model specification saying: Reached total allocation of 254Mb: see help(memory.size). here is part of the log: *** aphids-read.table(aphid.txt,header=T) attach(aphids) names(aphids) [1] site time shade treatp aphid library(nlme) aphids-groupedData(aphid~time|site/shade/treatp,data=aphids,outer=~shade*treatp) model-lme(aphid~shade*treatp,random=~time|site) Error in logLik.lmeStructInt(lmeSt, lmePars) : Calloc could not allocate (89359209 of 8) memory In addition: Warning message: Reached total allocation of 254Mb: see help(memory.size) *** My response variable is number of aphids/colony on shaded cocoa trees, measured over 15 days (which I treat as a random factor) in four sites (random factor), two levels of shade (fixed) within each site, and three aphid predation treatmens (fixed) within each shade level. My total N = 2,256 observations (actually counts). The comand lines were the same as in the example in: Crawley, M. J. 2002. Statistical Computing: an intro to data analysis using S-Plus. Wiley. pages: 702-704. My questions are: 1-Did I use the appropriate analysis to my data? 2-Is the lack of memory caused by my large N and many factor interactions? 3-Since I have counts, should I specify (family=poisson) in the model? If so, where in the command line this term must be located? I am a complete beginner with R, so any help will be wellcome. Evandro Silva Assistant Professor Universidade Estadual de Feira de Santana Laboratory of Insect Ecology Feira de Santana, BA, Brazil __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Glitch when creating online help
___ Ah, I think I've solved it. Apparently \describe needs to have a newline before it, or things get funny. Thus \arguments{ \item{item1}{ This is item 1. } \item{itemlist}{ Here is a list. \describe{ \item{subitem1}{Item 1 of the list.} \item{subitem2}{Item 2 of the list.} } } } works fine. OTOH, \details{ Here is another list. \describe{ \item{subitem1}{Item 1 of the list.} \item{subitem2}{Item 2 of the list.} } } (no nesting of \describe within an \item) gives the same glitch as described below. (Maybe I need to brush up on my TeX) -- Hong Ooi Senior Research Analyst, IAG Limited 388 George St, Sydney NSW 2000 +61 (2) 9292 1566 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Hong Ooi Sent: Wednesday, 14 December 2005 3:24 PM To: r-help@stat.math.ethz.ch Subject: [R] Glitch when creating online help ___ Hi, I'm writing up the online help for a package I'm developing (in-house only, sorry), and I've come across an odd glitch when trying to nest a list inside the arguments section of the .Rd file. I was just wondering if anyone could provide some insights. I'm using R 2.2.0 on Windows XP, along with ActivePerl 5.8.7 (build 815), MikTeX 2.4, and the tools downloaded from http://www.murdoch-sutherland.com/Rtools/ . Here is some code to reproduce the glitch. First, in R: f - function(x) x package.skeleton(foo, list=f) This creates the package skeleton, with a template f.Rd provided. Edit f.Rd to contain \name{f} \alias{f} \title{ ~~function to do ... ~~ } \description{ ~~ A concise (1-5 lines) description of what the function does. ~~ } \usage{f(x)} \arguments{ \item{item1}{ This is item 1. } \item{itemlist}{ Here is a list. \describe{ \item{subitem1}{Item 1 of the list.} \item{subitem2}{Item 2 of the list.} } } \item{item3}{ This is the item after the list. } } Then at the command prompt: R CMD INSTALL --build foo Once the package has been created, in R type: library(foo) ?f The result looks like fpackage:fooR Documentation ~~function to do ... ~~ Description: ~~ A concise (1-5 lines) description of what the function does. ~~ Usage: f(x) Arguments: item1: This is item 1. itemlist: Here is a list. .in +5 subitem1 Item 1 of the list. subitem2 Item 2 of the list. item3: This is the item after the list. Note the .in +5 at the top of the nested list. This is only in the online help within R, not the html version. -- Hong Ooi Senior Research Analyst, IAG Limited 388 George St, Sydney NSW 2000 +61 (2) 9292 1566 ___ The information transmitted in this message and its attachme...{{dropped}} __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] broken package?
On Tue, 13 Dec 2005, Jeff D. Hamann wrote: I've just tried to install the tseries package (which seems to include quite a bit) and got the following results: package 'dynlm' successfully unpacked and MD5 sums checked package 'leaps' successfully unpacked and MD5 sums checked package 'chron' successfully unpacked and MD5 sums checked package 'fCalendar' successfully unpacked and MD5 sums checked package 'strucchange' successfully unpacked and MD5 sums checked package 'DAAG' successfully unpacked and MD5 sums checked package 'quadprog' successfully unpacked and MD5 sums checked Error in sprintf(gettext(unable to move temp installation '%d' to '%s'), : use format %s for character objects I'm running R-2.1.0 Is this updated in 2.2? Aren't you supposed to find out if a problem is already fixed before posting, to use proper R version numbers and to tell us your OS? My copy of the posting guide says so in all three cases. What did yours say? The error message was changed very soon after 2.1.0 was released. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Ploting graphics using X tints from a color
On Tue, 13 Dec 2005, Sérgio Nunes wrote: Hi, I'm trying to draw a 2D plot using multiple tints of red. The (simplified) setup is the following: || year | x | y || My idea is that each year is plotted with a different tint of red. Older year (lightest) - Later year (darkest). I've managed to plot this with different scales of grays simply by doing: There have been other replies using rgb and hsv colour spaces directly, but the new colorRamp and colorRampPalette functions iare worth mentioning because they pack those solutions into a general framework: years - 1:20 plot(years, pch=19, col=grey(length(years):0/length(years))) redPalette - colorRampPalette(c(white, red)) plot(years, pch=19, col=redPalette(length(years))) and the grey values can be reconstructed by: greyPalette - colorRampPalette(c(white, black)) plot(years, pch=19, col=greyPalette(length(years))) with some values differing by 1 RGB unit. I've found colorRampPalette() very flexible, and that it supplements RColorBrewer well, so that intermediate colours can be interpolated for those palettes: library(RColorBrewer) redPalette - colorRampPalette(brewer.pal(5, Reds)) plot(years, pch=19, col=redPalette(length(years))) for a touch more tomato (subjectively) in the red. palette(gray(length(years):0/length(years))) before the plot and for each year the color used is a different tint of gray. So, is there any way to do this for any color? Any tip or advice? With this, I hope to visualize patterns in my dataset more easily. Thanks in advance for any help. Best regards, Sérgio Nunes __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html