Re: [R] date handling

2005-12-13 Thread paul sorenson
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

2005-12-13 Thread 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.

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

2005-12-13 Thread JeeBee

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

2005-12-13 Thread Doran, Harold
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?

2005-12-13 Thread Dan Bolser

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?

2005-12-13 Thread Hans Gardfjell
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

2005-12-13 Thread [EMAIL PROTECTED]
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

2005-12-13 Thread Davia Cox
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

2005-12-13 Thread Davia Cox
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

2005-12-13 Thread May, Roel
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

2005-12-13 Thread Ulrike Grömping
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-13 Thread Hans-Peter
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

2005-12-13 Thread Matteo Vidali
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

2005-12-13 Thread Oarabile Molaodi
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

2005-12-13 Thread Sean Davis



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

2005-12-13 Thread Uwe Ligges
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

2005-12-13 Thread Carlos Mauricio Cardeal Mendes
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...

2005-12-13 Thread Alvaro Saurin

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

2005-12-13 Thread Albert Sorribas

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...

2005-12-13 Thread Gabor Grothendieck
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

2005-12-13 Thread Diethelm Wuertz



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

2005-12-13 Thread 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] How to make a plot?

2005-12-13 Thread Gabor Grothendieck
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

2005-12-13 Thread Carlos Mauricio Cardeal Mendes
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

2005-12-13 Thread Doran, Harold
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

2005-12-13 Thread Matteo Vidali
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

2005-12-13 Thread Liaw, Andy
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

2005-12-13 Thread Liaw, Andy
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

2005-12-13 Thread Ted Harding
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

2005-12-13 Thread Peter Flom
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

2005-12-13 Thread Dimitris Rizopoulos
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

2005-12-13 Thread Romain Francois
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?

2005-12-13 Thread JeeBee

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

2005-12-13 Thread Thomas Lumley

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

2005-12-13 Thread Thomas Lumley
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

2005-12-13 Thread Sebastian Leuzinger
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

2005-12-13 Thread JeeBee
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

2005-12-13 Thread Ted Harding
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

2005-12-13 Thread Prof Brian Ripley
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

2005-12-13 Thread Ted Harding
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)

2005-12-13 Thread David Reitter
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

2005-12-13 Thread Albert Vilella
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)

2005-12-13 Thread Melanie Edwards
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

2005-12-13 Thread Ferdinand Alimadhi
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 (?)

2005-12-13 Thread Antonio, Fabio Di Narzo
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

2005-12-13 Thread Ted Harding
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...

2005-12-13 Thread Alvaro Saurin

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 (?)

2005-12-13 Thread Douglas Bates
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...

2005-12-13 Thread Gabor Grothendieck
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

2005-12-13 Thread Braesch Lucas
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

2005-12-13 Thread Ted Harding
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

2005-12-13 Thread Robert Burrows
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?

2005-12-13 Thread Romain Francois
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...

2005-12-13 Thread Alvaro Saurin

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

2005-12-13 Thread Peter Dalgaard
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?

2005-12-13 Thread t c
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

2005-12-13 Thread Uwe Ligges
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

2005-12-13 Thread Prof Brian Ripley
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

2005-12-13 Thread Gregory Snow
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

2005-12-13 Thread Marc Schwartz (via MN)
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

2005-12-13 Thread Jörg Schaber
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?

2005-12-13 Thread Marc Schwartz (via MN)
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...

2005-12-13 Thread Gabor Grothendieck
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

2005-12-13 Thread Singh, Jatinder
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

2005-12-13 Thread Sérgio Nunes
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

2005-12-13 Thread Ted Harding
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

2005-12-13 Thread Sérgio Nunes
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

2005-12-13 Thread Elizabeth Lawson
  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?

2005-12-13 Thread Ferdinand Alimadhi
?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

2005-12-13 Thread Ido M. Tamir
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...

2005-12-13 Thread Gabor Grothendieck
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?

2005-12-13 Thread ronggui
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

2005-12-13 Thread Romain Francois
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

2005-12-13 Thread hadley wickham
 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?

2005-12-13 Thread Trevor Hastie
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

2005-12-13 Thread Romain Francois
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?

2005-12-13 Thread Uwe Ligges
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

2005-12-13 Thread Uwe Ligges
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

2005-12-13 Thread Ben Bolker
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

2005-12-13 Thread Prof Brian Ripley
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

2005-12-13 Thread Meinhard Ploner

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

2005-12-13 Thread Diethelm Wuertz
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

2005-12-13 Thread jim holtman
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

2005-12-13 Thread Elizabeth Lawson
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)

2005-12-13 Thread Kjetil Brinchmann Halvorsen
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...

2005-12-13 Thread Kjetil Brinchmann Halvorsen
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

2005-12-13 Thread Berton Gunter
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 ?

2005-12-13 Thread Sérgio Nunes
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 (?)

2005-12-13 Thread Paulo Justiniano Ribeiro Jr

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

2005-12-13 Thread Sierra, Carlos
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

2005-12-13 Thread Spencer Graves
  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 ?

2005-12-13 Thread Jacques VESLOT
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

2005-12-13 Thread Hong Ooi

___


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?

2005-12-13 Thread Jeff D. Hamann
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'

2005-12-13 Thread Henrik Bengtsson
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)

2005-12-13 Thread Evandro do Nascimento Silva
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

2005-12-13 Thread Hong Ooi

___


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?

2005-12-13 Thread Prof Brian Ripley
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

2005-12-13 Thread Roger Bivand
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