Re: [R] question on function arguments

2008-02-18 Thread baptiste Auguié
any input?

(and also, why would all my posts get bounced off while I clearly  
subscribed to this list?)

Best regards,


On 18 Feb 2008, at 09:53, baptiste Auguié wrote:

> Hi,
>
> I have two small issues with my R code, no big deal but curiosity
> really. Here is a sample code,
>
>
>>
>> x <- rnorm(1:10)
>>
>> foo <- function(a = 1, b = list(x = c(1:10), y = c(1:10))){
>>  
>>  for (ii in seq(along=b$y)){
>>  
>>  print(x[ii] + b$x[ii])
>>  }
>>  
>>  
>> }
>>
>> foo() # default OK
>>
>> foo(b=list(x=1, y=c(1:10))) # only the first term works
>>
>> foo(b$y = 1:5) # error
>>
>
> In the last call, i wish to use all default arguments but b$y (that
> is, using the default a and b$x). Is this possible?
>
> The second call is more related to indices: i would like my argument b
> $x to be either a vector (default), or a scalar. In the latter, the
> loop b$x[ii] breaks when i would like it to recycle the single value.
> I can check the length of b$x with a "if" statement, but it becomes
> intricate when several arguments have this option of being vector or
> scalar. Is there a way to use b$x that would work in both cases?
>
>
> Many thanks,
>
> baptiste
>
> _
>
> Baptiste Auguié
>
> Physics Department
> University of Exeter
> Stocker Road,
> Exeter, Devon,
> EX4 4QL, UK
>
> Phone: +44 1392 264187
>
> http://newton.ex.ac.uk/research/emag
> http://projects.ex.ac.uk/atto
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting- 
> guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] R-problem

2008-02-18 Thread Theuser, Catrin
Dear all,
 
I encountered the following problem while using the "R" software: "unrecognised 
record type 7, subtype 16 encountered in system file".
Can you please help me with this?
 
Thanks in advance,
Catrin Theuser

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] rJava and matrices

2008-02-18 Thread Ben Woodcroft
Thanks for your prompt reply Martin. That should work, but still involves
changing the java code.

I was just hoping there might be an obvious solution I was missing. Sorry
for the mailing list wrongness. I didn't realise there was an rJava mailing
list.

Thanks,
ben

On 19/02/2008, Martin Morgan <[EMAIL PROTECTED]> wrote:
>
> Ben
>
> R represents matrices as a vectors (the result of as.vector(t) on
> your data below) with a 'dim' attribute, and I bet that this is how
> they're passed to Java. Write a method with sig
>
> public void data.structure.DistanceMatrix.setDistances(double[] data,
> int[] dim)
>
> invoked along the lines of
>
> .jcall(distanceMatrix, "V", "setDistances", as.vector(t), dim(t))
>
> (sorry, I'm a bit rusty on this so details might be incorrect)
>
> This probably belongs on R-devel, if not on the rJava mailing list.
>
> Martin
>
> "Ben Woodcroft" <[EMAIL PROTECTED]> writes:
>
> > Hi list,
> >
> > I'm wanting to use a matrix as input to some java code, but I seem to be
> > unable to do this (See code below). When searching online for a solution
> I
> > found that rJava 0.5-2 (the version under development not yet in CRAN)
> is
> > adding "direct support for raw vectors as method parameters". Is this
> > related?
> >
> > I also tried to use .jarray / .jcast in various combinations to no
> avail.
> >
> > The only way I can think of to solve this is to modify the java code to
> > allow me to upload a row of the matrix one at a time.
> >
> > Any other ideas?
> >
> >> require(rJava);
> > Loading required package: rJava
> >> t = rbind(c(1,2,3),c(2,3,4),c(1,2,10));
> >> .jinit()
> >> distanceMatrix = .jnew("data/structure/DistanceMatrix");
> >
> >> .jmethods(distanceMatrix);
> >  ... other methods omitted for brevity ...
> >  [4] "public void data.structure.DistanceMatrix.setDistances
> (double[][])"
> >
> >> .jcall(distanceMatrix, "V", "setDistances", as.matrix(t));
> > Error in .jcall(distanceMatrix, "V", "setDistances", as.matrix(t)) :
> >   method setDistances with signature ([D)V not found
> > Calls: .jcall -> .External
> > Execution halted
> >
> >
> >
> > Thanks in advance,
> > ben
> >
> >   [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> --
> Martin Morgan
> Computational Biology / Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N.
> PO Box 19024 Seattle, WA 98109
>
> Location: Arnold Building M2 B169
> Phone: (206) 667-2793
>



-- 
FYI: My email addresses at unimelb, uq and gmail all redirect to the same
place.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] rJava and matrices

2008-02-18 Thread Ben Woodcroft
Hi list,

I'm wanting to use a matrix as input to some java code, but I seem to be
unable to do this (See code below). When searching online for a solution I
found that rJava 0.5-2 (the version under development not yet in CRAN) is
adding "direct support for raw vectors as method parameters". Is this
related?

I also tried to use .jarray / .jcast in various combinations to no avail.

The only way I can think of to solve this is to modify the java code to
allow me to upload a row of the matrix one at a time.

Any other ideas?

> require(rJava);
Loading required package: rJava
> t = rbind(c(1,2,3),c(2,3,4),c(1,2,10));
> .jinit()
> distanceMatrix = .jnew("data/structure/DistanceMatrix");

> .jmethods(distanceMatrix);
 ... other methods omitted for brevity ...
 [4] "public void data.structure.DistanceMatrix.setDistances(double[][])"

> .jcall(distanceMatrix, "V", "setDistances", as.matrix(t));
Error in .jcall(distanceMatrix, "V", "setDistances", as.matrix(t)) :
  method setDistances with signature ([D)V not found
Calls: .jcall -> .External
Execution halted



Thanks in advance,
ben

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] tabulation on dataframe question

2008-02-18 Thread Moshe Olshansky
If d is your dataframe, how about

ind <- d$Dist == 2
aggregate(d,by=list(d$NameA,dNameA),FUN=length)

Regards,

Moshe.

--- Karin Lagesen <[EMAIL PROTECTED]>
wrote:

> 
> I have a data frame with data similar to this:
> 
> NameA  GrpA   NameB GrpB Dist
> A  Alpha  B Alpha   0.2
> A  Alpha  C Beta0.2
> A  Alpha  D Beta0.4
> B  Alpha  C Beta0.2
> B  Alpha  D Beta0.1
> C  Beta   D Beta0.3
> 
> Dist is a distance measure between two entities. The
> table displays
> all to all distances, but the distance between two
> entities only
> appears once. What I would like to get is a table
> where I get a count
> of entities per group where the distance satisfies a
> certain condition
> ( equal to zero for instance).
> 
> In this case, if the requirement was Distances ==
> 0.2
> 
> AlphaBeta
> Alpha   12
> Beta20
> 
> This resulting table would be symmetrical.
> 
> I hope I am able to convey what I would like, and
> TIA for your help!
> 
> Karin
> -- 
> Karin Lagesen, PhD student
> [EMAIL PROTECTED]
> http://folk.uio.no/karinlag
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained,
> reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Matrix inversion

2008-02-18 Thread Bill.Venables
Ben Domingue asks:

> I am trying to invert a matrix for the purposes of least squares.  I
> have tried a number of things, and the variety of results has me
> confused.

Don't be.

> 1. When I try solve() I get the following:
> >Error in solve.default(t(X) %*% X) : system is computationally
> singular: reciprocal condition number = 3.76391e-20

That seems clear enough.

> 2. When I try qr.solve(), I get:
> >Error in qr.solve(t(X) %*% X) : singular matrix 'a' in solve

That backs up what solve() said.

> 3. I can, however, use lm(y~X) to get coefficients.  This confuses me
> since I thought that lm() used qr().  So why did qr.solve() not work
> earlier?

It may use the QR decomposition, but that does not necessarily imply
that it uses the R function qr() in the same the way you did.  

If you look at the help information for lm() you will see that there 
is an argument singular.OK with the default value of "TRUE".  
Essentially this agrument allows you to say whether or not you 
anticipate the model matrix *could be* not of full column rank.  
If so, the regression coefficients are not unique and what you will 
see in the coefficients is one particular version.  The residuals and 
the residual sum of squares, however, are unique.

> 4. I have even tried using ginv().  The process works, but I end up
> with a different set of regression coefficients after I finish the
> process than what I had with lm().  To the best of my knowledge, this
> shouldn't happen.

Then you have learned something and it's a win.  If the model matrix
is not of full column rank, then the regression coefficients are not
unique.  This means that different solving methods can give different
answers.  Indeed you should expect this.  lm() uses the QR
decomposition; ginv() uses the singular value decomposition.  You
might like to see what aov() gives.  It could be the same as lm() or
it could be something different again.  All three, though, will give
you the same residuals up to computational accuracy.

> 
> I've been digging around all day and can't figure this out.  Thanks,
> 
> Ben Domingue
> PhD Student, School of Education
> University of Colorado at Boulder
> 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Softmax in nnet

2008-02-18 Thread G Ilhamto
Hi R help,

 I run my data in nnet with skip layer, factor response (with 0 & 1
 values) and explicitly put softmax=T to compare the result of the
 default nnet with no softmax specification. I assume this should give
 me the same result. I got the result the default one, but not the
 softmax version and I got the error message that I did not quite
 understand.

 test6.nn.skipT.softm.Yfac <- nnet(Yfac~ X1 +.. +X8, skip=T, size=0,
 softmax=T, data = train.set)

 Error in nnet.default(x, y, w, entropy = TRUE, ...) : no weights to fit

 In addition: Warning messages:
 1: In if (softmax) { :  the condition has length > 1 and only the
 first element will be used
 2: In if (skip) net <- add.net(net, seq(1, net$n[1]), seq(1 + net$n[1] +  :
  the condition has length > 1 and only the first element will be used

 When I specify the weights (0.1 just for a try)
 test6.nn.skipT.softm.Yfac <- nnet(Yfac~ X1 +.. +X8, skip=T, size=0,
 weights= 0.1, softmax=T, data = train.set)

 I got another error message:
 Error in model.frame.default(formula = Yfac ~ HusYEduc + AgeRespd +
 muslimat +  :
 variable lengths differ (found for '(weights)')
 -

 Q: Does softmax apply to two-category response?
 What is softmax require?

 Thank you,
 Ilham

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Matrix inversion

2008-02-18 Thread Ben Domingue
Howdy,
I am trying to invert a matrix for the purposes of least squares.  I
have tried a number of things, and the variety of results has me
confused.
1. When I try solve() I get the following:
>Error in solve.default(t(X) %*% X) : system is computationally
singular: reciprocal condition number = 3.76391e-20
2. When I try qr.solve(), I get:
>Error in qr.solve(t(X) %*% X) : singular matrix 'a' in solve
3. I can, however, use lm(y~X) to get coefficients.  This confuses me
since I thought that lm() used qr().  So why did qr.solve() not work
earlier?
4. I have even tried using ginv().  The process works, but I end up
with a different set of regression coefficients after I finish the
process than what I had with lm().  To the best of my knowledge, this
shouldn't happen.

I've been digging around all day and can't figure this out.  Thanks,

Ben Domingue
PhD Student, School of Education
University of Colorado at Boulder

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Hazard model with long-term survivor (cure model)

2008-02-18 Thread Gabor Grothendieck
Maybe:

http://post.queensu.ca/~pengp/software.html

On Feb 18, 2008 3:08 AM, 宋时歌 <[EMAIL PROTECTED]> wrote:
> Dear All,
>
> Are there R packages that can estimate survival model with long-term
> survivors? This is sometimes known as "cure" model or "split-population"
> model. Thanks.
>
> Shige
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] predicting memory usage

2008-02-18 Thread Moshe Olshansky
If your numbers are not integers then they are double
precision which means 8 bytes per number. In such a
case you will need 86000*2500*8 = 1.7Gb of memory
(plus a small amount for book keeping). You will need
more (and sometimes much more) if you intend to do
some operations on your arrays). So you probably will
need a 64 bit machine with enough RAM.
Integer arrays will take 4 bytes per entry, so about
half that amount.

Regards,

Moshe.

--- Federico Calboli <[EMAIL PROTECTED]> wrote:

> Hi All,
> 
> is there a way of predicting memory usage?
> 
> I need to build an array of 86000 by 2500 numbers
> (or I might create  
> a list of 2 by 2500 arrays 43000 long). How much
> memory should I  
> expect to use/need?
> 
> Cheers,
> 
> Fede
> 
> --
> Federico C. F. Calboli
> Department of Epidemiology and Public Health
> Imperial College, St. Mary's Campus
> Norfolk Place, London W2 1PG
> 
> Tel +44 (0)20 75941602   Fax +44 (0)20 75943193
> 
> f.calboli [.a.t] imperial.ac.uk
> f.calboli [.a.t] gmail.com
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained,
> reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Hazard model with long-term survivor (cure model)

2008-02-18 Thread David Winsemius
"=?UTF-8?B?5a6L5pe25q2M?=" <[EMAIL PROTECTED]> wrote in
news:[EMAIL PROTECTED]: 

> Dear All,
> 
> Are there R packages that can estimate survival model with long-term
> survivors? This is sometimes known as "cure" model or
> "split-population" model. Thanks.

The usual Cox model would certainly allow the analysis of observations 
with long-term survivors, but I am wondering if you want some sort of 
parametric model or one which compares a cohort's survival to some sort 
of external standard population expected survival. Therneau and Gramsch's 
text has a chapter on working with such expected survival estimates and 
the survival package can probably be considered an R atandard. If you set 
up a parametric Weibull model with a decreasing hazard, you get a cure, 
at least asymptotically. I think Harrell's Design package can conjure up 
accelerated time models that include the Weibull. 

Although the SRAB cancer methodologists assert that: "Neither SAS nor 
Splus can be used to fit survival models to relative survival data." (1)
I would be very dubious regarding such a claim. See for instance:





-- 
David Winsemius

1) 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] newbie (me) needs to model distribution as two overlapping gaussians

2008-02-18 Thread Wensui Liu
it could be a finite mixture.
take a look at flexmix package.

On Feb 18, 2008 11:03 AM,  <[EMAIL PROTECTED]> wrote:
>
> Recently, I have been working with some data that look like two overlapping 
> gaussian distributions.  I would like to either
>
> 1) determine the mean and SD for each of the two distributions
>
> OR
>
> 2) get some (bayesian ?) statistic that estimates how likely an observation 
> is to belong to the left-hand or right-hand distribution
>
> In case I'm using the wrong language, my data looks something like this:
>
> B <- rnorm(500,40,10)
> H <- rnorm(500,80,5 )
> N <- runif(200,0,99)
> D <- c(B,H,N)
>
> Where B=background, H=hits, N=noise, and D=my observed distribution
>
> I have seen analyses like this in the past, but I can't remember what it is 
> called.  If somebody out there can point me towards an R function, or even 
> the cannonical name for this kind of model, I think I can write the necessary 
> code.
>
> Thanks in advance,
> Mark
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
===
WenSui Liu
ChoicePoint Precision Marketing
Phone: 678-893-9457
Email : [EMAIL PROTECTED]
Blog   : statcompute.spaces.live.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Compare mean survival time

2008-02-18 Thread Ben Bolker
Xing Yuan  gmail.com> writes:

> 
> Dear List,
> 
> Does anybody no how to compare mean survival times for two (more) groups in
> R? What test statistics should I use?
> 
> Thank you very much!
> 
> Joe

  The question is probably a little bit too vague.  You should
certainly look at the survival package (which is 'recommended',
so included with R) -- if you really have individual-level survival
data, then survival analysis is more powerful than comparing mean
survival times (but not necessarily the same thing).  Try
Therneau and Grambsch (listed on the "books" page at www.r-project.org) ...

   Ben Bolker

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] index range

2008-02-18 Thread Ben Bolker

Or see the "Oarray" package (note that's
a letter O, not a zero!)

  Ben Bolker

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] index range

2008-02-18 Thread Gabor Grothendieck
On Feb 18, 2008 6:52 PM, Gabor Grothendieck <[EMAIL PROTECTED]> wrote:
> You can define origin 0 objects yourself if you like.
> Here is a partial implementation:
>
> "[.orig0" <- function(x, i)

The 0 should be a 1.

>if (is.numeric(i)) .subset(x, i+0) else .subset(x, i)
> orig0 <- function(x)
>structure(x, class = c("orig0", setdiff(class(x), "orig0")))
> x <- orig0(1:5)
> x[0:3] # 1:4
>
> Note that usually 0 means leave out that element in R and
> in this implementation -1 means leave it out.  Also normally
> -3 mean leave out 3rd element but in the implementation
> above -0 would be 0 so it would give the first element.
>
> Probably best to just get used to the R way.
>
>
> On Feb 18, 2008 6:31 PM, hill0093 <[EMAIL PROTECTED]> wrote:
> >
> > It looks to me like the index range starts at 1 in R.
> > Is this true?
> > If so, is there a way to change it to start at 0?
> > That way, I wouldn't have to make so many
> > changes when I translate a function from
> > another language.
> > --
> > View this message in context: 
> > http://www.nabble.com/index-range-tp15550797p15550797.html
> > Sent from the R help mailing list archive at Nabble.com.
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] index range

2008-02-18 Thread Gabor Grothendieck
You can define origin 0 objects yourself if you like.
Here is a partial implementation:

"[.orig0" <- function(x, i)
if (is.numeric(i)) .subset(x, i+0) else .subset(x, i)
orig0 <- function(x)
structure(x, class = c("orig0", setdiff(class(x), "orig0")))
x <- orig0(1:5)
x[0:3] # 1:4

Note that usually 0 means leave out that element in R and
in this implementation -1 means leave it out.  Also normally
-3 mean leave out 3rd element but in the implementation
above -0 would be 0 so it would give the first element.

Probably best to just get used to the R way.

On Feb 18, 2008 6:31 PM, hill0093 <[EMAIL PROTECTED]> wrote:
>
> It looks to me like the index range starts at 1 in R.
> Is this true?
> If so, is there a way to change it to start at 0?
> That way, I wouldn't have to make so many
> changes when I translate a function from
> another language.
> --
> View this message in context: 
> http://www.nabble.com/index-range-tp15550797p15550797.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Huge number

2008-02-18 Thread jim holtman
If you want to compute (157+221)! then sum up the log:

> a <- 1:(157+221)
> sum(log10(a))
[1] 811.8165

This is about 6.55e811 which exceeds the range of floating point
numbers (1.797693e+308).  You might check out the Brobdingnag package.

On Feb 18, 2008 6:23 PM, Hyojin Lee <[EMAIL PROTECTED]> wrote:
> Hi,
> I'm trying to calculate p-value to findout definitely expressed genes
> compare A to B situation.
> I got this data(this is a part of data) from whole organism , and each
> number means each expression values (that means, we could think 'a' gene
> is 13 in A situation, and it turns 30 in B situation)
> To findout probability, I'm going to use Audic - Claverie Method. ( The
> significance of digital gene expression profiles. 1997)
>
> But using this equation p(x|y), I have to calculate (x+y)!  first. but I
> can't  calculate (157+221)! or (666+1387)! in R.
> That's probabily the handling problem of huge number, How could I
> calculate p value in this data with R?
>
>
> A B
> Total5874641 6295980
> a13  30
> b36  39
> c0   5
> d40  61
> e16  20
> f13  11
> g3   3
> h9   5
> i12  35
> j157 221
> k17  39
> l6   17
> m666 1387
> n2   5
>
>
>
>
>
>
>
>
>
> The significance of digital gene expression profiles.
>
> Audic S
>  udic%20S%22%5BAuthor%5D&itool=EntrezSystem2.PEntrez.Pubmed.Pubmed_Result
> sPanel.Pubmed_RVAbstractPlusDrugs1> , Claverie JM
>  laverie%20JM%22%5BAuthor%5D&itool=EntrezSystem2.PEntrez.Pubmed.Pubmed_Re
> sultsPanel.Pubmed_RVAbstractPlusDrugs1> .
>
>[[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem you are trying to solve?

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] index range

2008-02-18 Thread hill0093

It looks to me like the index range starts at 1 in R.
Is this true?
If so, is there a way to change it to start at 0?
That way, I wouldn't have to make so many
changes when I translate a function from
another language.
-- 
View this message in context: 
http://www.nabble.com/index-range-tp15550797p15550797.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Huge number

2008-02-18 Thread Hyojin Lee
Hi, 
I'm trying to calculate p-value to findout definitely expressed genes
compare A to B situation.
I got this data(this is a part of data) from whole organism , and each
number means each expression values (that means, we could think 'a' gene
is 13 in A situation, and it turns 30 in B situation) 
To findout probability, I'm going to use Audic - Claverie Method. ( The
significance of digital gene expression profiles. 1997)
 
But using this equation p(x|y), I have to calculate (x+y)!  first. but I
can't  calculate (157+221)! or (666+1387)! in R.
That's probabily the handling problem of huge number, How could I
calculate p value in this data with R? 
 
 
 A B 
Total5874641 6295980
a13  30 
b36  39 
c0   5  
d40  61 
e16  20 
f13  11 
g3   3  
h9   5  
i12  35 
j157 221
k17  39 
l6   17 
m666 1387   
n2   5  
 
 
 
 
 
 
 
 

The significance of digital gene expression profiles.

Audic S
 , Claverie JM
 .

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] mean and variance of ratio

2008-02-18 Thread Moshe Olshansky
Hi Irene,

The result depends on the joint distribution of
x1,x2,y1,y2. The mean (or the variance) not always
exist (if, for example, all of them are independent
and normally distributed).
One possibility to estimate mean and variance is to
use simulation (once you verify that they exist!!!).

Regards,

Moshe.

--- Irene Mantzouni <[EMAIL PROTECTED]> wrote:

> Hi all!
> 
>  
> 
> I try to estimate a statistic of the form:
> (x1-x2)/(y1-y2), where
> x1,x2,y1,y2 represent variable means, so each has an
> estimate and
> standard error associated with it. 
> 
> How is it possible to estimate the mean and the
> variance of this ratio? 
> 
>  
> 
> Thank you!
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained,
> reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] newbie (me) needs to model distribution as two overlapping gaussians

2008-02-18 Thread David Winsemius
<[EMAIL PROTECTED]> wrote in
news:[EMAIL PROTECTED]: 

> 
> Recently, I have been working with some data that look like two
> overlapping gaussian distributions.  I would like to either 
> 
> 1) determine the mean and SD for each of the two distributions
> 
> OR
> 
> 2) get some (bayesian ?) statistic that estimates how likely an
> observation is to belong to the left-hand or right-hand distribution 
> 
> In case I'm using the wrong language, my data looks something like
> this: 
> 
> B <- rnorm(500,40,10)
> H <- rnorm(500,80,5 )
> N <- runif(200,0,99)
> D <- c(B,H,N)

First hit on  with: functions mixture gaussians



Task Views:




nor1mix: Normal (1-d) Mixture Models (S3 Classes and Methods)
  

In the nor1mix documentation, Mächler recommends the mclust package 
for estimation.

-- 
David Winsemius

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Huge RDCOM.err file with Bloomberg download

2008-02-18 Thread Jozef Thomas Kohut
Hello,

 

I use Rbloomberg package to download bloomberg data. Other required packages
are RDCOMClient, zoo and chron. The problem I have is  an error file
(RDCOM.err) that is created by starting the download process. The growth of
the file increases with each downloaded ticker and series by 2 Kb.  After
downloading 1200 tickers, each 3 series (3600 datasets), the growth of the
file is 7.2 Gb.

The file is on the  C: folder located. 

 

The contents of RDCOM.err (for each download process) is: 

 

Type of argument 14

Type of argument 14

Type of argument 14

Type of argument 14

Type of argument 16

Type of argument 16

Type of argument 14

Converting VARIANT to R 8204

Finishing convertDCOMObjectToR - convert array

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 7

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

Converting VARIANT to R 5

Finished convertDCOMObjectToR

 

Does somebody have the same problem ? Could somebody give me an advice, how
to prevent creating this file?

 

Thank you for your advice.

Best regards

 

Jozef Thomas Kohut

Email:   [EMAIL PROTECTED]

 


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] remove column names from a data frame

2008-02-18 Thread Charilaos Skiadas
I can't resist asking, why would you want to remove the variable  
names? If it is for printing purposes, then you can probably work  
around it on the printing side, depending on what you want to print.  
I can't think of another reason for wanting to remove the column  
names altogether.

Haris Skiadas
Department of Mathematics and Computer Science
Hanover College

On Feb 18, 2008, at 3:08 PM, Benilton Carvalho wrote:

> Well, Henrique's solution is nicer than mine. :-)
> b
>
> On Feb 18, 2008, at 3:05 PM, joseph wrote:
>
>> Hi Benilton
>> I actually tried "names(df) <- NULL" before I asked for help.  
>> Please see below how it looks like
>>
>> > names(df) <- NULL
>> > df
>>   structure(c("chr1", "chr2", "chr3"), class = "AsIs") structure(c 
>> ("1", "2", "3"), class = "AsIs")
>> 1  
>> chr1   1
>> 2  
>> chr2   2
>> 3  
>> chr3   3
>>   structure(c("4", "5", "6"), class = "AsIs") structure(c("7",  
>> "8", "9"), class = "AsIs")
>> 1
>> 4   7
>> 2
>> 5   8
>> 3
>> 6   9
>> >
>>
>>
>> - Original Message 
>> From: Benilton Carvalho <[EMAIL PROTECTED]>
>> To: joseph <[EMAIL PROTECTED]>
>> Cc: r-help@r-project.org
>> Sent: Monday, February 18, 2008 11:47:50 AM
>> Subject: Re: [R] remove column names from a data frame
>>
>> apparently you want to check the "Introduction to R" document I
>> found it very useful when I started working with R:
>>
>> http://cran.r-project.org/doc/manuals/R-intro.pdf
>>
>> try:
>>
>> names(df) <- NULL
>>
>> b
>>
>> ps: "df" is the name of the function to get the density for an F
>> distribution...
>>
>> On Feb 18, 2008, at 2:42 PM, joseph wrote:
>>
>> >
>> >
>> > I want to remove the column names from a data frame. I do
>> > it the long way, can any body show me a better way ?
>> >
>> >
>> > df= data.frame(chrN= c(“chr1”, “chr2”, “chr3”), start= c(1,
>> > 2, 3), end= c(4, 5, 6), score= c(7, 8, 9))
>> >
>> >
>> > df
>> >
>> >
>> > #I write a txt file without row or column names
>> >
>> >
>> > write.table 
>> (df,"df1.txt",sep='\t',quote=FALSE,row.names=F,col.names=F)
>> >
>> >
>> > #then I read it with the header = F to obtain the format
>> > I want
>> >
>> >
>> > df1=read.table("df1.txt", header = F)
>> >
>> >
>> > df1
>> >
>> >
>> >
>> >  
>> _ 
>> ___
>> > Looking for last minute shopping deals?
>> >
>> > /category.php?category=shopping
>> > [[alternative HTML version deleted]]
>> >
>> > __
>> > R-help@r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-help
>> > PLEASE do read the posting guide http://www.R-project.org/ 
>> posting-guide.html
>> > and provide commented, minimal, self-contained, reproducible code.
>>
>>
>>
>> Never miss a thing. Make Yahoo your homepage.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Compare mean survival time

2008-02-18 Thread Xing Yuan
Dear List,

Does anybody no how to compare mean survival times for two (more) groups in
R? What test statistics should I use?

Thank you very much!

Joe

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] remove column names from a data frame

2008-02-18 Thread Benilton Carvalho

Well, Henrique's solution is nicer than mine. :-)
b

On Feb 18, 2008, at 3:05 PM, joseph wrote:


Hi Benilton
I actually tried "names(df) <- NULL" before I asked for help. Please  
see below how it looks like


> names(df) <- NULL
> df
  structure(c("chr1", "chr2", "chr3"), class = "AsIs")  
structure(c("1", "2", "3"), class = "AsIs")
1  
chr1   1
2  
chr2   2
3  
chr3   3
  structure(c("4", "5", "6"), class = "AsIs") structure(c("7", "8",  
"9"), class = "AsIs")
1
4   7
2
5   8
3
6   9

>


- Original Message 
From: Benilton Carvalho <[EMAIL PROTECTED]>
To: joseph <[EMAIL PROTECTED]>
Cc: r-help@r-project.org
Sent: Monday, February 18, 2008 11:47:50 AM
Subject: Re: [R] remove column names from a data frame

apparently you want to check the "Introduction to R" document I
found it very useful when I started working with R:

http://cran.r-project.org/doc/manuals/R-intro.pdf

try:

names(df) <- NULL

b

ps: "df" is the name of the function to get the density for an F
distribution...

On Feb 18, 2008, at 2:42 PM, joseph wrote:

>
>
> I want to remove the column names from a data frame. I do
> it the long way, can any body show me a better way ?
>
>
> df= data.frame(chrN= c(“chr1”, “chr2”, “chr3”), start= c(1,
> 2, 3), end= c(4, 5, 6), score= c(7, 8, 9))
>
>
> df
>
>
> #I write a txt file without row or column names
>
>
>  
write.table(df,"df1.txt",sep='\t',quote=FALSE,row.names=F,col.names=F)

>
>
> #then I read it with the header = F to obtain the format
> I want
>
>
> df1=read.table("df1.txt", header = F)
>
>
> df1
>
>
>
>  


> Looking for last minute shopping deals?
>
> /category.php?category=shopping
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



Never miss a thing. Make Yahoo your homepage.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] remove column names from a data frame

2008-02-18 Thread joseph
Hi Benilton
I actually tried "names(df) <- NULL" before I asked for help. Please see below 
how it looks like

> names(df) <- NULL
> df
  structure(c("chr1", "chr2", "chr3"), class = "AsIs") structure(c("1", "2", 
"3"), class = "AsIs")
1 chr1  
 1
2 chr2  
 2
3 chr3  
 3
  structure(c("4", "5", "6"), class = "AsIs") structure(c("7", "8", "9"), class 
= "AsIs")
1   4   
7
2   5   
8
3   6   
9
> 


- Original Message 
From: Benilton Carvalho <[EMAIL PROTECTED]>
To: joseph <[EMAIL PROTECTED]>
Cc: r-help@r-project.org
Sent: Monday, February 18, 2008 11:47:50 AM
Subject: Re: [R] remove column names from a data frame


apparently 
you 
want 
to 
check 
the 
"Introduction 
to 
R" 
document 
I  
found 
it 
very 
useful 
when 
I 
started 
working 
with 
R:

http://cran.r-project.org/doc/manuals/R-intro.pdf

try:

names(df) 
<- 
NULL

b

ps: 
"df" 
is 
the 
name 
of 
the 
function 
to 
get 
the 
density 
for 
an 
F  
distribution...

On 
Feb 
18, 
2008, 
at 
2:42 
PM, 
joseph 
wrote:

>
>
> 
I 
want 
to 
remove 
the 
column 
names 
from 
a 
data 
frame. 
I 
do
> 
it 
the 
long 
way, 
can 
any 
body 
show 
me 
a 
better 
way 
?
>
>
> 
df= 
data.frame(chrN= 
c(“chr1”, 
“chr2”, 
“chr3”), 
start= 
c(1,
> 
2, 
3), 
end= 
c(4, 
5, 
6), 
score= 
c(7, 
8, 
9))
>
>
> 
df
>
>
> 
#I 
write 
a 
txt 
file 
without 
row 
or 
column 
names
>
>
> 
write.table(df,"df1.txt",sep='\t',quote=FALSE,row.names=F,col.names=F)
>
>
> 
#then 
I 
read 
it 
with 
the 
header 
= 
F 
to 
obtain 
the 
format
> 
I 
want
>
>
> 
df1=read.table("df1.txt", 
header 
= 
F)
>
>
> 
df1
>
>
>  
  
  
 
> 

> 
Looking 
for 
last 
minute 
shopping 
deals?
>
> 
/category.php?category=shopping
> 

[[alternative 
HTML 
version 
deleted]]
>
> 
__
> 
R-help@r-project.org 
mailing 
list
> 
https://stat.ethz.ch/mailman/listinfo/r-help
> 
PLEASE 
do 
read 
the 
posting 
guide 
http://www.R-project.org/posting-guide.html
> 
and 
provide 
commented, 
minimal, 
self-contained, 
reproducible 
code.







  

Never miss a thing.  Make Yahoo your home page. 

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] newbie (me) needs to model distribution as two overlapping gaussians

2008-02-18 Thread tcornulier
Is a finite mixture of 2 gaussians the name you are looking for?
This specific model will not deal with your N component however.

you can find some functions here:
http://cran.r-project.org/web/views/Cluster.html

Thomas

[EMAIL PROTECTED] wrote:
> Recently, I have been working with some data that look like two overlapping 
> gaussian distributions.  I would like to either
> 
> 1) determine the mean and SD for each of the two distributions
> 
> OR
> 
> 2) get some (bayesian ?) statistic that estimates how likely an observation 
> is to belong to the left-hand or right-hand distribution
> 
> In case I'm using the wrong language, my data looks something like this:
> 
> B <- rnorm(500,40,10)
> H <- rnorm(500,80,5 )
> N <- runif(200,0,99)
> D <- c(B,H,N)
> 
> Where B=background, H=hits, N=noise, and D=my observed distribution
> 
> I have seen analyses like this in the past, but I can't remember what it is 
> called.  If somebody out there can point me towards an R function, or even 
> the cannonical name for this kind of model, I think I can write the necessary 
> code.
> 
> Thanks in advance,
> Mark
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 

-- 


**
Thomas Cornulier

Mammal Research Institute
Polish Academy of Sciences
ul. Waszkiewicza 1c
17-230 Bialowieza
http://www.zbs.bialowieza.pl/?en

tel (0048) 85 682 77 88

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] remove column names from a data frame

2008-02-18 Thread Henrique Dallazuanna
Try this:

names(df) <- NA

or

names(df) <- make.names(seq(ncol(df)))


On 18/02/2008, joseph <[EMAIL PROTECTED]> wrote:
>
>
> I want to remove the column names from a data frame. I do
> it the long way, can any body show me a better way ?
>
>
> df= data.frame(chrN= c("chr1", "chr2", "chr3"), start= c(1,
> 2, 3), end= c(4, 5, 6), score= c(7, 8, 9))
>
>
> df
>
>
> #I write a txt file without row or column names
>
>
> write.table(df,"df1.txt",sep='\t',quote=FALSE,row.names=F,col.names=F)
>
>
> #then I read it with the header = F to obtain the format
> I want
>
>
> df1=read.table("df1.txt", header = F)
>
>
> df1
>
>
>   
> 
> Looking for last minute shopping deals?
>
> /category.php?category=shopping
> [[alternative HTML version deleted]]
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>


-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] remove column names from a data frame

2008-02-18 Thread Benilton Carvalho
apparently you want to check the "Introduction to R" document I  
found it very useful when I started working with R:


http://cran.r-project.org/doc/manuals/R-intro.pdf

try:

names(df) <- NULL

b

ps: "df" is the name of the function to get the density for an F  
distribution...


On Feb 18, 2008, at 2:42 PM, joseph wrote:




I want to remove the column names from a data frame. I do
it the long way, can any body show me a better way ?


df= data.frame(chrN= c(“chr1”, “chr2”, “chr3”), start= c(1,
2, 3), end= c(4, 5, 6), score= c(7, 8, 9))


df


#I write a txt file without row or column names


write.table(df,"df1.txt",sep='\t',quote=FALSE,row.names=F,col.names=F)


#then I read it with the header = F to obtain the format
I want


df1=read.table("df1.txt", header = F)


df1


  


Looking for last minute shopping deals?

/category.php?category=shopping
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] remove column names from a data frame

2008-02-18 Thread joseph


I want to remove the column names from a data frame. I do
it the long way, can any body show me a better way ?


df= data.frame(chrN= c(“chr1”, “chr2”, “chr3”), start= c(1,
2, 3), end= c(4, 5, 6), score= c(7, 8, 9))


df


#I write a txt file without row or column names


write.table(df,"df1.txt",sep='\t',quote=FALSE,row.names=F,col.names=F)


#then I read it with the header = F to obtain the format
I want


df1=read.table("df1.txt", header = F)


df1


  

Looking for last minute shopping deals?  

/category.php?category=shopping
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] skip non-converging nls() in a list

2008-02-18 Thread Gabor Grothendieck
See warnOnly in ?nls.control

On Feb 18, 2008 1:45 PM, Marc Belisle <[EMAIL PROTECTED]> wrote:
> Howdee,
>
> My question appears at #6 below:
>
> 1. I want to model the growth of each of a large number of individuals using
> a 4-parameter logistic growth curve.
>
> 2. nlme does not converge with the random structure that I want to use.
>
> 3. nlsList does not converge for some individuals.
>
> 4. I decided to go around nlsList using:
>
> t(sapply(split(data, list(data$id)),
>function(subd){coef(nls(mass ~ SSfpl(age, A, B, xmid, scal), data =
> subd))}))
>
> 5. This does not converge either:
>
> 'Error in nls(y ~ cbind(1, 1/(1 + exp((xmid - x)/exp(lscal, data = xy,
> :
>singular gradient'
>
> 6. Would anyone know if I can get R to skip non-converging use of nls() so
> that I can at least obtain the parameters for the curves which R can fit?
>
> Thanks for your time,
>
> Marc
>
> ===
> Marc Bélisle
> Professeur adjoint
> Chaire de recherche du Canada en écologie spatiale et en écologie du paysage
> Département de biologie
> Université de Sherbrooke
> 2500 Boul. de l'Université
> Sherbrooke, Québec
> J1K 2R1 Canada
>
> Tél: +1-819-821-8000 poste 61313
> Fax: +1-819-821-8049
> Courriél: [EMAIL PROTECTED]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] skip non-converging nls() in a list

2008-02-18 Thread Andrew Robinson
Hi Marc,

use try()

Cheers

Andrew

On Mon, Feb 18, 2008 at 01:45:35PM -0500, Marc Belisle wrote:
> Howdee,
> 
> My question appears at #6 below:
> 
> 1. I want to model the growth of each of a large number of individuals using
> a 4-parameter logistic growth curve.
> 
> 2. nlme does not converge with the random structure that I want to use.
> 
> 3. nlsList does not converge for some individuals.
> 
> 4. I decided to go around nlsList using:
> 
> t(sapply(split(data, list(data$id)),
> function(subd){coef(nls(mass ~ SSfpl(age, A, B, xmid, scal), data =
> subd))}))
> 
> 5. This does not converge either:
> 
> 'Error in nls(y ~ cbind(1, 1/(1 + exp((xmid - x)/exp(lscal, data = xy,
> :
> singular gradient'
> 
> 6. Would anyone know if I can get R to skip non-converging use of nls() so
> that I can at least obtain the parameters for the curves which R can fit?
> 
> Thanks for your time,
> 
> Marc
> 
> ===
> Marc B?lisle
> Professeur adjoint
> Chaire de recherche du Canada en ?cologie spatiale et en ?cologie du paysage
> D?partement de biologie
> Universit? de Sherbrooke
> 2500 Boul. de l'Universit?
> Sherbrooke, Qu?bec
> J1K 2R1 Canada
> 
> T?l: +1-819-821-8000 poste 61313
> Fax: +1-819-821-8049
> Courri?l: [EMAIL PROTECTED]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> -- 
> This message has been scanned for viruses and
> dangerous content by MailScanner, and is
> believed to be clean.

-- 
Andrew Robinson  
Department of Mathematics and StatisticsTel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] skip non-converging nls() in a list

2008-02-18 Thread Marc Belisle
Howdee,

My question appears at #6 below:

1. I want to model the growth of each of a large number of individuals using
a 4-parameter logistic growth curve.

2. nlme does not converge with the random structure that I want to use.

3. nlsList does not converge for some individuals.

4. I decided to go around nlsList using:

t(sapply(split(data, list(data$id)),
function(subd){coef(nls(mass ~ SSfpl(age, A, B, xmid, scal), data =
subd))}))

5. This does not converge either:

'Error in nls(y ~ cbind(1, 1/(1 + exp((xmid - x)/exp(lscal, data = xy,
:
singular gradient'

6. Would anyone know if I can get R to skip non-converging use of nls() so
that I can at least obtain the parameters for the curves which R can fit?

Thanks for your time,

Marc

===
Marc Bélisle
Professeur adjoint
Chaire de recherche du Canada en écologie spatiale et en écologie du paysage
Département de biologie
Université de Sherbrooke
2500 Boul. de l'Université
Sherbrooke, Québec
J1K 2R1 Canada

Tél: +1-819-821-8000 poste 61313
Fax: +1-819-821-8049
Courriél: [EMAIL PROTECTED]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] paste("Mus., 10 ", expression(mu)," g", sep="")

2008-02-18 Thread Henrique Dallazuanna
Try this:

plot(rnorm(10), main=expression(paste("Mus., 10 ",mu," g", sep="")))

On 18/02/2008, Samor Gandhi <[EMAIL PROTECTED]> wrote:
> Dear all,
>
>   I am very thankful, if you could tell wheather it is possible to write
>
>   paste("Mus., 10 ", expression(mu)," g", sep="")
>
>   Thank you in advance,
>   Samor
>
>
> -
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] progress bar

2008-02-18 Thread Syed Abid Hussaini
Hi all, 
I have some functions which sometimes take longer to run and it would
 be useful to have a progress
bar showing how much time is left for the script to finish (something
 like a download progress
bar). I tried searching but could not find it. Is this possible with R?

Note: Using R version: 2.6.1 on windows xp.


  

Never miss a thing.  Make Yahoo your home page.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] paste("Mus., 10 ", expression(mu)," g", sep="")

2008-02-18 Thread Samor Gandhi
Dear all,
   
  I am very thankful, if you could tell wheather it is possible to write
   
  paste("Mus., 10 ", expression(mu)," g", sep="") 
   
  Thank you in advance,
  Samor

   
-

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] progress bar

2008-02-18 Thread Henrique Dallazuanna
See RSiteSearch("Progress Bar")

On 18/02/2008, Syed Abid Hussaini <[EMAIL PROTECTED]> wrote:
> I have some functions which sometimes take longer to run and it would be 
> useful to have a progress
> bar showing how much time is left for the script to finish (something like a 
> download progress
> bar). I tried searching but could not find it. Is this possible with R?
>
> Note: Using R version: 2.6.1 on windows xp.
>
>
>   
> 
> Be a better friend, newshound, and
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] View body of externally compiled functions

2008-02-18 Thread David Winsemius
"Michael Gormley" <[EMAIL PROTECTED]> wrote in
news:[EMAIL PROTECTED]: 

> This is a simple question.  With functions written and compiled in
> R, the body of the function can be seen by typing the function name
> without any input arguments or by using the body function.  How can
> I look at the body of code written and compiled externally and
> called using the .Call function? 

Look at the sources.

-- 
David Winsemius

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] progress bar

2008-02-18 Thread Syed Abid Hussaini
I have some functions which sometimes take longer to run and it would be useful 
to have a progress
bar showing how much time is left for the script to finish (something like a 
download progress
bar). I tried searching but could not find it. Is this possible with R? 

Note: Using R version: 2.6.1 on windows xp.


  

Be a better friend, newshound, and

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] View body of externally compiled functions

2008-02-18 Thread Michael Gormley
This is a simple question.  With functions written and compiled in R, the
body of the function can be seen by typing the function name without any
input arguments or by using the body function.  How can I look at the body
of code written and compiled externally and called using the .Call function?


Thanks,
Mike Gormley

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] predicting memory usage

2008-02-18 Thread jim holtman
If this is numeric, then for just storing one copy, you will require
86000 * 2500 * 8 = 1.7GB of memory.  You should have 3-4X that if you
want to analyze it, so you might need about 6GB of physical memory and
a 64-bit version of R.  Is there some other alternative?  Do you need
all the values at once, or can you use a database to access the
portions you want?

On 2/18/08, Federico Calboli <[EMAIL PROTECTED]> wrote:
> Hi All,
>
> is there a way of predicting memory usage?
>
> I need to build an array of 86000 by 2500 numbers (or I might create
> a list of 2 by 2500 arrays 43000 long). How much memory should I
> expect to use/need?
>
> Cheers,
>
> Fede
>
> --
> Federico C. F. Calboli
> Department of Epidemiology and Public Health
> Imperial College, St. Mary's Campus
> Norfolk Place, London W2 1PG
>
> Tel +44 (0)20 75941602   Fax +44 (0)20 75943193
>
> f.calboli [.a.t] imperial.ac.uk
> f.calboli [.a.t] gmail.com
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390

What is the problem you are trying to solve?

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] predicting memory usage

2008-02-18 Thread Federico Calboli
Hi All,

is there a way of predicting memory usage?

I need to build an array of 86000 by 2500 numbers (or I might create  
a list of 2 by 2500 arrays 43000 long). How much memory should I  
expect to use/need?

Cheers,

Fede

--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] "Error in package manager"

2008-02-18 Thread laure gallien
Hi all,

I am quite new to R, and I guess i made a mistake with
it... 
Now, when i want to open the package manager i got
this :

"Erreur dans package.manager(is.loaded, pkgs,
pkgs.desc, pkgs.url) : 
  invalid arguments (length mismatch)
De plus : Warning message:
In .find.package(pkgs) : aucun package nommé
'file56e509fe' n'est trouvé"

Well, it is in french... saying that it didn't find
the package call "ile56e509fe".

Can you help me ??
I already tried to download a new R version (for my
MacOS X) and the last R GUI...And the problem is still
there...


Thanks a lot,

Laure

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Number of digits of a value for problem 7.31 RFAQ SOLVED not really

2008-02-18 Thread Richard . Cotton
> Actually (after some trials) there is a little problem when faced with 
> zeros...
> 
>  >getndp(1.0)
> [1] 0

Are you sure this isn't what you want?  1.0 is just 1 in disguise, and 
round(1.0, 0) is the same as round(1.0, 1) anyway.

> Note that I thought on a very different way which was starting from the 
> point that as.character seems to work:
> 
> a<-as.character(1.23134)
> b<-strsplit(a,"\\.")
> 
> nchar(b)

I think the last line should be:

nchar(b[[1]][2])

If you choose this route, be careful with R clipping your values.
e.g.
 x <- 0.12345678901234567890123456789

format(x, nsmall=20)
[1] "0.12345678901234567737" #this is what is stored (actually some binary 
number, this is the closest decimal representation)

a <- as.character(x); a
[1] "0.123456789012346" #

b<-strsplit(a,"\\.")
nchar(b[[1]][2])
[1] 15

Regards,
Richie.

Mathematical Sciences Unit
HSL



ATTENTION:

This message contains privileged and confidential inform...{{dropped:20}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] can we include nonparametric component for survival regression?

2008-02-18 Thread Thomas Lumley
On Mon, 18 Feb 2008, gallon li wrote:

> i am trying to fit a survival regression model (cox model or parametric
> model) in R by including the covariate effects as a function m(x) instead of
> just beta*x. is it possible to fit such a model? can someone recommend some
> reference? I searched but only found a package called addreg where
> the hazard is actually modeled additively. That is not what i want.

The survival package has pspline() for this purpose (assuming that by 
'non-parametric' you mean a flexible smooth curve)

-thomas

Thomas Lumley   Assoc. Professor, Biostatistics
[EMAIL PROTECTED]   University of Washington, Seattle

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] predict.lm with matrix as newdata

2008-02-18 Thread Uwe Ligges


Marlin Keith Cox wrote:
> Thank you in advance for helping me with this.
> I included a part of the actual data.  The result of pred.est and
> pred.est1are different, but they should be identical.  For
> pred.est, I entered the slope and y intercept and received a value for each
> individual number in the matrix (Z).  For pred.est1, I was wanting R to do
> this for me and it seems that it combined values within each row of the
> matrix (Z), this was indicated by the both the row numbers and pred.est1 =
> 30.
> 
> rm(list=ls())
> 
> #Water
> w = c(9.037929, 8.273594,   9.253578,   8.039490,   8.381091,  11.610205,
>  11.261445,  11.185045,   9.903959,  12.307910,  10.307602,  13.950863,
> 13.028331,  13.677622,  13.051269,  13.289698,  14.391914,  21.411512,
> 75.568802,  78.168239,  60.829664, 112.558057, 127.296835, 108.739532,
> 133.516712, 130.492614, 129.783800, 168.289704, 152.732359, 168.405646)
> 
> #Resistance
> R=c(660, 548, 676, 763, 768, 692, 657, 630, 748, 680, 786, 645, 710, 677,
> 692, 732, 737, 651, 396,
> 601, 640, 448, 464, 472, 434, 487, 495, 426, 429, 456)
> 
> #Detector length
> Lend=c(37.0,  39.0,  39.0,  39.0,  40.0,  41.5,  44.0,  45.0,  46.0,  47.0,
> 47.0,  48.0,
> 48.5,  49.0,  51.0,  53.0, 53.0,  60.0,  89.0, 103.0, 108.5, 118.0, 118.0,
> 123.0,
> 126.0, 138.0, 139.0, 141.0, 141.0, 151.0)
> 
> #Errors to be multiplied by Restistance
> x=c(0,.05,.10,.15,.20,.25)
> 
> #Errors to be multiplied by Detector length
> y=c(0,.01,.02,.03,.04,.05)
> 
> 
> #equation to predict water weight in grams
> Rs<-(Lend^2)/R
> model.lm<-lm(w~Rs)
> a=predict.lm(model.lm)
> 
> 
> 
> X=(R%o%x+R)
> 
> Y=((Lend%o%y+Lend)^2)
> X
> Y
> num.x.col <- length(X[1,])
> num.y.col <- length(Y[1,])
> num.rows <- length(X[,1])
> 
> Z <- matrix(nrow=num.rows, ncol=num.x.col*num.y.col)
> 
> for( i in 1:num.rows)  {
>Z[i,] <- as.vector( Y[i,] %*% t(X[i,])^-1 )
> }
> Z
> pred.est <- 3.453*Z+1.994


What did you try above? What has this Z to do with the linear model?
Why does Z correspond to the former Rs - it is really far off the 
original values from Rs? Pay attention here.

I fear you are confusing some things here. You might want to explain 
what X, Y and Z are good for.


> pred.est1 <- predict.lm(model.lm, data.frame(Rs=Z))

If you want to apply predict() [please, not predict.lm()] on all columns 
of Z, try:

apply(Z, 2, function(x) predict(model.lm, newdata = data.frame(Rs=x)))

Uwe Ligges




> pred.est
> pred.est1
> 
> detach(z)
> detach(sen)
> 
> On Sat, Feb 16, 2008 at 2:10 AM, Uwe Ligges <[EMAIL PROTECTED]>
> wrote:
> 
>>
>> Marlin Keith Cox wrote:
>>> Z is a matrix and when I run the following line, it creates a prediction
>>> estimate using each column, how can I get it an estimate for each
>> individual
>>> number.  I have tried changing Z to a data.frame, but this does not do
>> it
>>> either.
>>>
>>> model.lm<-lm(w~x)
>>>
>>> pred.est <- predict.lm(model.lm, data.frame(x=Z))
>> That's what I do and it always worked for me. Can you please specify a
>> reproducible example that shows what you got and tells us exactly what
>> you expected?
>>
>> BTW: It is expected that you call the generic predict() rather than its
>> particular method.
>>
>> Uwe Ligges
>>
>>
>>> Thanks in advance,
>>> keith
>>>
> 
> 
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] predict.lm with matrix as newdata

2008-02-18 Thread Marlin Keith Cox
Thank you in advance for helping me with this.
I included a part of the actual data.  The result of pred.est and
pred.est1are different, but they should be identical.  For
pred.est, I entered the slope and y intercept and received a value for each
individual number in the matrix (Z).  For pred.est1, I was wanting R to do
this for me and it seems that it combined values within each row of the
matrix (Z), this was indicated by the both the row numbers and pred.est1 =
30.

rm(list=ls())

#Water
w = c(9.037929, 8.273594,   9.253578,   8.039490,   8.381091,  11.610205,
 11.261445,  11.185045,   9.903959,  12.307910,  10.307602,  13.950863,
13.028331,  13.677622,  13.051269,  13.289698,  14.391914,  21.411512,
75.568802,  78.168239,  60.829664, 112.558057, 127.296835, 108.739532,
133.516712, 130.492614, 129.783800, 168.289704, 152.732359, 168.405646)

#Resistance
R=c(660, 548, 676, 763, 768, 692, 657, 630, 748, 680, 786, 645, 710, 677,
692, 732, 737, 651, 396,
601, 640, 448, 464, 472, 434, 487, 495, 426, 429, 456)

#Detector length
Lend=c(37.0,  39.0,  39.0,  39.0,  40.0,  41.5,  44.0,  45.0,  46.0,  47.0,
47.0,  48.0,
48.5,  49.0,  51.0,  53.0, 53.0,  60.0,  89.0, 103.0, 108.5, 118.0, 118.0,
123.0,
126.0, 138.0, 139.0, 141.0, 141.0, 151.0)

#Errors to be multiplied by Restistance
x=c(0,.05,.10,.15,.20,.25)

#Errors to be multiplied by Detector length
y=c(0,.01,.02,.03,.04,.05)


#equation to predict water weight in grams
Rs<-(Lend^2)/R
model.lm<-lm(w~Rs)
a=predict.lm(model.lm)



X=(R%o%x+R)

Y=((Lend%o%y+Lend)^2)
X
Y
num.x.col <- length(X[1,])
num.y.col <- length(Y[1,])
num.rows <- length(X[,1])

Z <- matrix(nrow=num.rows, ncol=num.x.col*num.y.col)

for( i in 1:num.rows)  {
   Z[i,] <- as.vector( Y[i,] %*% t(X[i,])^-1 )
}
Z
pred.est <- 3.453*Z+1.994
pred.est1 <- predict.lm(model.lm, data.frame(Rs=Z))
pred.est
pred.est1

detach(z)
detach(sen)

On Sat, Feb 16, 2008 at 2:10 AM, Uwe Ligges <[EMAIL PROTECTED]>
wrote:

>
>
> Marlin Keith Cox wrote:
> > Z is a matrix and when I run the following line, it creates a prediction
> > estimate using each column, how can I get it an estimate for each
> individual
> > number.  I have tried changing Z to a data.frame, but this does not do
> it
> > either.
> >
> > model.lm<-lm(w~x)
> >
> > pred.est <- predict.lm(model.lm, data.frame(x=Z))
>
> That's what I do and it always worked for me. Can you please specify a
> reproducible example that shows what you got and tells us exactly what
> you expected?
>
> BTW: It is expected that you call the generic predict() rather than its
> particular method.
>
> Uwe Ligges
>
>
> > Thanks in advance,
> > keith
> >
>



-- 
Keith Cox, Ph.D.
Sitka Sound Science Center
Fisheries Biologist
P.O. Box 464
Sitka, Alaska, 99835

907 752-0563
[EMAIL PROTECTED]

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to change the axis labels of a wind rose plot?

2008-02-18 Thread David Winsemius
Huilin Chen <[EMAIL PROTECTED]> wrote in
news:[EMAIL PROTECTED]: 


> I am make plots using the function windrose() from the package
> "circular". I would like to have the north on the top instead of
> having 0 degree on the right side.
> Does anybody know how to change the labels in a windrose plot?

There are probably better ways, but this is the record of my experiments 
starting with the help file example:

set.seed(1)
dir <- circular(runif(100, 0, 360), units="degrees")
mag <-  rgamma(100, 15)
sample <- data.frame(dir=dir, mag=mag)
 
par(mfrow=c(2,2))
res <- windrose(sample)
#Rotates the data
sample$dir<-sample$dir+90
res <- windrose(sample)
axis.circular(labels=c("E","N","W","S"))
#Overlayed the letter labels but leaves the degree labels

#Read the help file again and now suppress the axes 
res <- windrose(sample,axes=FALSE)
# and correctly(?) position labels
axis.circular(labels=c("E","N","W","S"),tcl.text = -0.1)

-- 
David Winsemius

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] newbie (me) needs to model distribution as two overlapping gaussians

2008-02-18 Thread mmiller

Recently, I have been working with some data that look like two overlapping 
gaussian distributions.  I would like to either

1) determine the mean and SD for each of the two distributions

OR

2) get some (bayesian ?) statistic that estimates how likely an observation is 
to belong to the left-hand or right-hand distribution

In case I'm using the wrong language, my data looks something like this:

B <- rnorm(500,40,10)
H <- rnorm(500,80,5 )
N <- runif(200,0,99)
D <- c(B,H,N)

Where B=background, H=hits, N=noise, and D=my observed distribution

I have seen analyses like this in the past, but I can't remember what it is 
called.  If somebody out there can point me towards an R function, or even the 
cannonical name for this kind of model, I think I can write the necessary code.

Thanks in advance,
Mark

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] how to plot image() without painting a map (the background)

2008-02-18 Thread hadley wickham
> I would like to plot akilonlat03 and then akilonlat06 and keep the map of
> France in background.
> My script (see below) doesn't work as image "paints" the background as I
> read somewhere in this forum.

I don't see how you can plot both temperatures on the same plot -
won't they be on top of one another?  (And it seems obvious that you
should want to plot the map last, rather than first)

Hadley

-- 
http://had.co.nz/

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] tabulation on dataframe question

2008-02-18 Thread Henrique Dallazuanna
Try this:

with(x,
{tmp <- table(x[Dist==0.2,c('GrpA', 'GrpB')])
 tmp[lower.tri(tmp)] <- tmp[upper.tri(tmp)]
 tmp})

On 18/02/2008, Karin Lagesen <[EMAIL PROTECTED]> wrote:
>
> I have a data frame with data similar to this:
>
> NameA  GrpA   NameB GrpB Dist
> A  Alpha  B Alpha   0.2
> A  Alpha  C Beta0.2
> A  Alpha  D Beta0.4
> B  Alpha  C Beta0.2
> B  Alpha  D Beta0.1
> C  Beta   D Beta0.3
>
> Dist is a distance measure between two entities. The table displays
> all to all distances, but the distance between two entities only
> appears once. What I would like to get is a table where I get a count
> of entities per group where the distance satisfies a certain condition
> ( equal to zero for instance).
>
> In this case, if the requirement was Distances == 0.2
>
> AlphaBeta
> Alpha   12
> Beta20
>
> This resulting table would be symmetrical.
>
> I hope I am able to convey what I would like, and TIA for your help!
>
> Karin
> --
> Karin Lagesen, PhD student
> [EMAIL PROTECTED]
> http://folk.uio.no/karinlag
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] how to plot image() without painting a map (the background)

2008-02-18 Thread Ptit_Bleu

Hello,

I'm trying to plot dayly evolution of the temperature over France from
Global Forecast System files ("I'm trying" is the right expression...).

akilonlat03 is the temperature for different latitudes and longitudes à 3
o'clock.
akilonlat06 is the temperature for different latitudes and longitudes à 6
o'clock.

I would like to plot akilonlat03 and then akilonlat06 and keep the map of
France in background.
My script (see below) doesn't work as image "paints" the background as I
read somewhere in this forum.

So would have someone a solution to correct my script ?
Thanks in advance,
Ptit Bleu.




akilonlat03<-interp(lonlat03$Longitude, lonlat03$Latitude,
(5/9)*(lonlat03$TMP_200802150300-32))
akilonlat06<-interp(lonlat06$Longitude, lonlat06$Latitude,
(5/9)*(lonlat06$TMP_200802150600-32))

map('france')

image(akilonlat03, col=cm.colors(100), axes=FALSE, add=TRUE)
#contour(akilonlat03, col="blue", add=TRUE)
image(akilonlat06, col=cm.colors(100), axes=FALSE, add=TRUE)
#contour(akilonlat06, col="blue", add=TRUE)

-- 
View this message in context: 
http://www.nabble.com/how-to-plot-image%28%29-without-painting-a-map-%28the-background%29-tp15546906p15546906.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] tabulation on dataframe question

2008-02-18 Thread Gabor Grothendieck
Try this:

# read test data
Lines <- "NameA  GrpA   NameB GrpB Dist
A  Alpha  B Alpha   0.2
A  Alpha  C Beta0.2
A  Alpha  D Beta0.4
B  Alpha  C Beta0.2
B  Alpha  D Beta0.1
C  Beta   D Beta0.3
"
DF <- read.table(textConnection(Lines), header = TRUE)

# process
x <- xtabs(~ GrpA + GrpB, DF, Dist == 0.2)
x + t(x) - diag(diag(x))


On Feb 18, 2008 10:16 AM, Karin Lagesen <[EMAIL PROTECTED]> wrote:
>
> I have a data frame with data similar to this:
>
> NameA  GrpA   NameB GrpB Dist
> A  Alpha  B Alpha   0.2
> A  Alpha  C Beta0.2
> A  Alpha  D Beta0.4
> B  Alpha  C Beta0.2
> B  Alpha  D Beta0.1
> C  Beta   D Beta0.3
>
> Dist is a distance measure between two entities. The table displays
> all to all distances, but the distance between two entities only
> appears once. What I would like to get is a table where I get a count
> of entities per group where the distance satisfies a certain condition
> ( equal to zero for instance).
>
> In this case, if the requirement was Distances == 0.2
>
>AlphaBeta
> Alpha   12
> Beta20
>
> This resulting table would be symmetrical.
>
> I hope I am able to convey what I would like, and TIA for your help!
>
> Karin
> --
> Karin Lagesen, PhD student
> [EMAIL PROTECTED]
> http://folk.uio.no/karinlag
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] How to change the axis labels of a wind rose plot?

2008-02-18 Thread Huilin Chen
Dear all,

I am make plots using the function windrose() from the package "circular".
I would like to have the north on the top instead of having 0 degree on 
the right side.
Does anybody know how to change the labels in a windrose plot?


thanks,
Huilin

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] repeated measures ANOVA using a covariance matrix

2008-02-18 Thread sun71 sun
Dear all,
how can I perform a repeated measures ANOVA using a covariance matrix as input?
E.g., I have four repeated measures (N = 200) with mean vector tau (no
BS factor):

tau <- rbind(1.10, 2.51, 2.76, 3.52)

and covariance matrix (sigma):

sigma <- matrix(c(0.523, 0.268, 0.267, 0.211,
   0.268, 0.444, 0.492, 0.571,
   0.267, 0.492, 1.213, 1.112,
   0.211, 0.571, 1.112, 1.811), nrow = 4, ncol
= 4, byrow = TRUE)

Thank you very much in advance!

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Number of digits of a value for problem 7.31 RFAQ SOLVED not really

2008-02-18 Thread Matthieu Stigler
Matthieu Stigler a e'crit :
> [EMAIL PROTECTED] a e'crit :
 What I mean is if R shows 2.3456 I want to obtain the info that 
>> digits=4
 even if in facts the value has more (internal) digits.
>>> Try:
>>> x = 1.23456789
>>> format(x, nsmall=20)
>>> # [1] "1.2345678899989009"
>>
>> I've just re-read the question. I suspect what you really wanted was 
>> something like:
>>
>> getndp <- function(x, tol=2*.Machine$double.eps)
>> {
>> ndp <- 0
>> while(!isTRUE(all.equal(x, round(x, ndp), tol=tol))) ndp <- ndp+1 
>> if(ndp > -log10(tol)) warning("Tolerance reached, ndp possibly 
>> underestimated.")
>> ndp }
>>
>> x = 0.123456
>> getndp(x)
>> # [1] 6
>>
>> x2 <- 3
>> getndp(x2)
>> # [1] 0
>> x3 <- 0.1234567890123456789
>> getndp(x3)
>> # [1] 16
>> # Warning message:
>> # In getndp(x3) : Tolerance reached, ndp possibly underestimated.
>>
>> Regards,
>> Richie.
>>
>> Mathematical Sciences Unit
>> HSL
>>
>
> Nice! Thanks a lot!!
> Now I can solve my problem that:
>
> target<-0.223423874568437
> target<=target+0.1-0.1 #TRUE
> target>=target+0.1-0.1 #FALSE!!
> target==target+0.1-0.1 #FALSE!!
>
> by:
> target<=round(target+0.1-0.1, getndp(target)) #TRUE
> target>=round(target+0.1-0.1, getndp(target)) #TRUE!!
> target==round(target+0.1-0.1, getndp(target)) #TRUE!!
>
>
Actually (after some trials) there is a little problem when faced with 
zeros...

 >getndp(1.0)
[1] 0

Note that it is only a partial problem because this happens only when 
the value is taken with other and so 
max(apply(matrix(allvalues,1,getndp))) will give the expected answer.

But when faced with only one value...

Note that I thought on a very different way which was starting from the 
point that as.character seems to work:

a<-as.character(1.23134)
b<-strsplit(a,"\\.")

nchar(b)
But I don't see how to do so that I can extract only the second element 
(maybe otherwise with sub or grep?)... But the method you gave is really 
more reliable than my "makeshift repair"

Thanks again!

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] tabulation on dataframe question

2008-02-18 Thread Karin Lagesen

I have a data frame with data similar to this:

NameA  GrpA   NameB GrpB Dist
A  Alpha  B Alpha   0.2
A  Alpha  C Beta0.2
A  Alpha  D Beta0.4
B  Alpha  C Beta0.2
B  Alpha  D Beta0.1
C  Beta   D Beta0.3

Dist is a distance measure between two entities. The table displays
all to all distances, but the distance between two entities only
appears once. What I would like to get is a table where I get a count
of entities per group where the distance satisfies a certain condition
( equal to zero for instance).

In this case, if the requirement was Distances == 0.2

AlphaBeta
Alpha   12
Beta20

This resulting table would be symmetrical.

I hope I am able to convey what I would like, and TIA for your help!

Karin
-- 
Karin Lagesen, PhD student
[EMAIL PROTECTED]
http://folk.uio.no/karinlag

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] CRAN Taskviews returns 404

2008-02-18 Thread Achim Zeileis
On Mon, 18 Feb 2008 14:15:54 + Neil Shephard wrote:

> On Feb 18, 2008 2:02 PM, Peter Dalgaard <[EMAIL PROTECTED]>
> wrote:
> >
> > But what gave you the idea that
> > http://cran.r-project.org/web/packages/Views/ should work? Google
> > seems not to know it.
> >
>
> Its the target for the link to the TaskViews from
> http://cran.r-project.org/web/packages/index.html (which was reached
> itself by following the "packages" link at
> http://cran.r-project.org/).
>
> Google does know about http://cran.r-project.org/src/contrib/Views/
> but that simply redirects to the current location.
>
> Its a tangled web out there!

The "web" directory on CRAN was introduced two days ago. Hence the
problems with some old links and the old information on Google. I've
corrected the task view link and a few others. Let us know if there are
further problems or broken links.

thx,
Z

> Neil
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html and provide commented,
> minimal, self-contained, reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] R graphics question: "binary" bar chart

2008-02-18 Thread Henrique Dallazuanna
If I undestand your question you can try something about like this:

x <- matrix(rbinom(75, 1,.6), nc=3)
image(x, col=c("gray", "white"), axes=F)
axis(1, at=seq(0,1, l=nrow(x)), labels=paste("Day", 1:nrow(x)))
axis(2, at=seq(0,1, l=ncol(x)), labels=paste("P", 1:3, sep=""))


On 18/02/2008, Stephan Spat <[EMAIL PROTECTED]> wrote:
> Hello!
>
> I would like to visualize the hospitalization within one year of several
> patients using a bar chart. For each patient the stay in a hospital
> should be illustrated with a dark colour a if there is a stay at home
> between 2 hopital stays, it should be illustrated with a bright colour.
>
> e.g.
>
> P1   |//|  |///| ||
> P2 |//|
> P3  |//|   ||   |///||//|
> Day1Day100-Day365
>
>
> legend:
> Px   Patient with id x
> |///|Hospital stay (length of bar symbolizes length of stay)
>   Place between 2 hospital stays symbolizes stay at home
>  (not in hospital)
>
> Are there any package which can do that? If not, how can I implement
> this in R?
>
> Thanks a lot
>
> --
> Stephan Spat
>
> Institute of Medical Technologies and Health Management
> JOANNEUM RESEARCH Forschungsgesellschaft mbH
> Elisabethstrasse 11a, 8010 Graz, AUSTRIA
>
> phone:+43 (0) 316 876-2157
> fax:  +43 (0) 316 876-2130
> email:[EMAIL PROTECTED]
> homepage: http://www.joanneum.at/msg
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] mean and variance of ratio

2008-02-18 Thread Christian Ritz
Dear Irene,

if you have a vector of estimates of x1, x2, y1, y2 and the 
corresponding estimated variance-covariance matrix (accessible through 
the method "vcov"), then one possibility is to use the function 
delta.method() in the package alr3 (on CRAN).

This function returns:

1) an estimate of the ratio (which is simply obtained by plugging in 
the estimates of x1, x2, y1, y2) and

2) an approximate, large-sample standard error of the estimated ratio




Christian

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] R graphics question: "binary" bar chart

2008-02-18 Thread Stephan Spat
Hello!

I would like to visualize the hospitalization within one year of several
patients using a bar chart. For each patient the stay in a hospital
should be illustrated with a dark colour a if there is a stay at home
between 2 hopital stays, it should be illustrated with a bright colour.

e.g.

P1   |//|  |///| ||
P2 |//|
P3  |//|   ||   |///||//|
Day1Day100-Day365


legend:
Px   Patient with id x
|///|Hospital stay (length of bar symbolizes length of stay)
  Place between 2 hospital stays symbolizes stay at home
 (not in hospital)

Are there any package which can do that? If not, how can I implement
this in R?

Thanks a lot

-- 
Stephan Spat

Institute of Medical Technologies and Health Management
JOANNEUM RESEARCH Forschungsgesellschaft mbH
Elisabethstrasse 11a, 8010 Graz, AUSTRIA

phone:+43 (0) 316 876-2157
fax:  +43 (0) 316 876-2130
email:[EMAIL PROTECTED]
homepage: http://www.joanneum.at/msg

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] mean and variance of ratio

2008-02-18 Thread Dimitris Rizopoulos
two options are the Delta method (see, e.g., function deltamethod() in 
package 'msm'), and the Bootstrap method (check package 'boot').

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://med.kuleuven.be/biostat/
 http://www.student.kuleuven.be/~m0390867/dimitris.htm


- Original Message - 
From: "Irene Mantzouni" <[EMAIL PROTECTED]>
To: <[EMAIL PROTECTED]>
Sent: Monday, February 18, 2008 3:22 PM
Subject: [R] mean and variance of ratio


> Hi all!
>
>
>
> I try to estimate a statistic of the form: (x1-x2)/(y1-y2), where
> x1,x2,y1,y2 represent variable means, so each has an estimate and
> standard error associated with it.
>
> How is it possible to estimate the mean and the variance of this 
> ratio?
>
>
>
> Thank you!
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] mean and variance of ratio

2008-02-18 Thread Irene Mantzouni
Hi all!

 

I try to estimate a statistic of the form: (x1-x2)/(y1-y2), where
x1,x2,y1,y2 represent variable means, so each has an estimate and
standard error associated with it. 

How is it possible to estimate the mean and the variance of this ratio? 

 

Thank you!


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] CRAN Taskviews returns 404

2008-02-18 Thread Neil Shephard
On Feb 18, 2008 2:02 PM, Peter Dalgaard <[EMAIL PROTECTED]> wrote:
>
> But what gave you the idea that
> http://cran.r-project.org/web/packages/Views/ should work? Google seems
> not to know it.
>

Its the target for the link to the TaskViews from
http://cran.r-project.org/web/packages/index.html (which was reached
itself by following the "packages" link at
http://cran.r-project.org/).

Google does know about http://cran.r-project.org/src/contrib/Views/
but that simply redirects to the current location.

Its a tangled web out there!

Neil

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to make a vector/list/array of POSIXlt object?

2008-02-18 Thread Gabor Grothendieck
Your mapply creates a list so to do it that way convert result to a
POSIXct vector first.  Using your strptime2 and dd, tt & dat from mine:

dt <- do.call(c, mapply(strptime2, dd, tt, SIMPLIFY = FALSE, USE.NAMES = FALSE))
data.frame(dt, dat)

# or just
data.frame(dt = strptime2(dd, tt), dat)




On Feb 18, 2008 8:37 AM, Bo Zhou <[EMAIL PROTECTED]> wrote:
>
> Ta. I will give that code a bash.
>
> Could you explain why my code didn't work?
>
>
>
>
>
> 
> > Date: Mon, 18 Feb 2008 00:25:44 -0500
>
>
> > From: [EMAIL PROTECTED]
> > To: [EMAIL PROTECTED]
> > Subject: Re: [R] How to make a vector/list/array of POSIXlt object?
> > CC: r-help@r-project.org
> >
> > If the problem is that you have a vector of dates, a vector of times
> > and a vector of data and you want to create a data frame with one
> > POSIXct column and one column of data then try this:
> >
> > dd <- c("01/22/2008", "02/13/2008")
> > tt <- c("01:01:00", "23:01:12")
> > dat <- 1:2
> >
> > data.frame(dt = strptime(paste(dd, tt), "%m/%d/%Y %H:%M:%S"), dat)
> >
> > # if you don't need subsecond data or time zones you could use chron
> >
> > library(chron)
> > data.frame(dt = chron(dd, tt), dat)
> >
> > If this is intended to be a time series you might want to look at the zoo
> > package. It has three vignettes that give more info.
> >
> >
> > On Feb 17, 2008 11:54 PM, Bo Zhou <[EMAIL PROTECTED]> wrote:
> > > Hi Gabor,
> > >
> > > I'm using this code but it doesn't work for me
> > >
> > > > strptime2<-function (date,time) as.POSIXct(strptime(paste(date,time),
> > > "%m/%d/%Y %H:%M:%S"))
> > > > dt=mapply(strptime2, "01/01/2008", "00:00:00", SIMPLIFY=FALSE,
> > > USE.NAMES=FALSE)
> > > > df=data.frame(X1=dt,X2=1)
> > > > dt
> > > [[1]]
> > > [1] "2008-01-01 Eastern Standard Time"
> > >
> > > > df
> > > structure.1199163600..class...c..POSIXtPOSIXcttzone.. X2
> > > 1 2008-01-01 1
> > >
> > > Here df looks very wrong to me.
> > >
> > > So I tested this code:
> > >
> > > > df2=data.frame(X1=as.POSIXct(Sys.time()),X2=1)
> > > > df2
> > > X1 X2
> > > 1 2008-02-17 23:43:08 1
> > > > class(df2$X1)
> > > [1] "POSIXt" "POSIXct"
> > >
> > > Ah this worked as I expected.
> > >
> > > So some tweaking here - SIMPLIFY is set TRUE now:
> > >
> > > > strptime2<-function (date,time) as.POSIXct(strptime(paste(date,time),
> > > "%m/%d/%Y %H:%M:%S"))
> > > > dt=mapply(strptime2, "01/01/2008", "00:00:00", SIMPLIFY=TRUE,
> > > USE.NAMES=FALSE)
> > > > df=data.frame(X1=dt,X2=1)
> > > > df
> > > X1 X2
> > > 1 1199163600 1
> > > > class(df$X1)
> > > [1] "numeric"
> > >
> > > Hmm... it worked, but not in a way I wanted. The class info is missing.
> > >
> > > So how to get the result like this below? I do need that mapply +
> > > strptime(paste), cos my CSV file is formatted in that way!
> > >
> > > > df2
> > > X1 X2
> > > 1 2008-02-17 23:43:08 1
> > > > class(df2$X1)
> > > [1] "POSIXt" "POSIXct"
> > >
> > >
> > > Any insight?
> > >
> > > Cheers,
> > >
> > > Bo
> > >
> > >
> > >
> > > > Date: Sun, 17 Feb 2008 15:53:28 -0500
> > > > From: [EMAIL PROTECTED]
> > > > To: [EMAIL PROTECTED]
> > > > Subject: Re: [R] How to make a vector/list/array of POSIXlt object?
> > > > CC: r-help@r-project.org
> > >
> > >
> > > >
> > > > Normally one uses POSIXct rather than POSIXlt for storage. See R News
> 4/1
> > > for
> > > > more info on date and time classes.
> > > >
> > > > On Feb 17, 2008 3:45 PM, Bo Zhou <[EMAIL PROTECTED]> wrote:
> > > > >
> > > > > Hi Guys,
> > > > >
> > > > > I'm cooking up my time series code. I want a data frame with first
> > > column as timestamp in POSIXlt format.
> > > > >
> > > > > I hit on this the problem of how to create an array/list/vector of
> > > POSIXlt objects. Code is as follows
> > > > >
> > > > >
> > > > >
> > > > > > dtt=array(dim = 2)
> > > > > > t=as.POSIXlt( strptime("07/12/07 13:20:01", "%m/%d/%Y
> > > %H:%M:%S",tz="GMT"))
> > > > > > dtt
> > > > > [1] NA NA
> > > > > > t
> > > > > [1] "0007-07-12 13:20:01 GMT"
> > > > > > dtt[1]=t
> > > > > Warning message:
> > > > > In dtt[1] = t :
> > > > > number of items to replace is not a multiple of replacement length
> > > > > > class(dtt)
> > > > > [1] "list"
> > > > > > class(t)
> > > > > [1] "POSIXt" "POSIXlt"
> > > > > > unclass(t)
> > > > > $sec
> > > > > [1] 1
> > > > >
> > > > > $min
> > > > > [1] 20
> > > > >
> > > > > $hour
> > > > > [1] 13
> > > > >
> > > > > $mday
> > > > > [1] 12
> > > > >
> > > > > $mon
> > > > > [1] 6
> > > > >
> > > > > $year
> > > > > [1] -1893
> > > > >
> > > > > $wday
> > > > > [1] 4
> > > > >
> > > > > $yday
> > > > > [1] 192
> > > > >
> > > > > $isdst
> > > > > [1] 0
> > > > >
> > > > > attr(,"tzone")
> > > > > [1] "GMT"
> > > > >
> > > > >
> > > > >
> > > > > Seems like POSIXlt is matrix in this case.
> > > > >
> > > > > Any suggestions?
> > > > >
> > > > > Cheers,
> > > > >
> > > > > B
> > > > > _
> > > > > [[elided Hotmail spam]]
> > > > >
> > > > > [[alter

Re: [R] Number of digits of a value for problem 7.31 RFAQ SOLVED

2008-02-18 Thread Matthieu Stigler
[EMAIL PROTECTED] a e'crit :
>>> What I mean is if R shows 2.3456 I want to obtain the info that 
>>>   
> digits=4 
>   
>>> even if in facts the value has more (internal) digits.
>>>   
>> Try:
>> x = 1.23456789
>> format(x, nsmall=20)
>> # [1] "1.2345678899989009"
>> 
>
> I've just re-read the question.  I suspect what you really wanted was 
> something like:
>
> getndp <- function(x, tol=2*.Machine$double.eps)
> {
>   ndp <- 0
>   while(!isTRUE(all.equal(x, round(x, ndp), tol=tol))) ndp <- ndp+1 
>   if(ndp > -log10(tol)) warning("Tolerance reached, ndp possibly 
> underestimated.")
>   ndp 
> }
>
> x = 0.123456
> getndp(x)
> # [1] 6
>
> x2 <- 3
> getndp(x2)
> # [1] 0 
>
> x3 <- 0.1234567890123456789
> getndp(x3)
> # [1] 16
> # Warning message:
> # In getndp(x3) : Tolerance reached, ndp possibly underestimated.
>
> Regards,
> Richie.
>
> Mathematical Sciences Unit
> HSL
>
>   

Nice! Thanks a lot!!
Now I can solve my problem that:

target<-0.223423874568437
target<=target+0.1-0.1 #TRUE
target>=target+0.1-0.1 #FALSE!!
target==target+0.1-0.1 #FALSE!!

by:
target<=round(target+0.1-0.1, getndp(target)) #TRUE
target>=round(target+0.1-0.1, getndp(target)) #TRUE!!
target==round(target+0.1-0.1, getndp(target)) #TRUE!!

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] CRAN Taskviews returns 404

2008-02-18 Thread Peter Dalgaard
Peter Dalgaard wrote:
> Neil Shephard wrote:
>   
>> The CRAN Taskviews (http://cran.r-project.org/web/packages/Views/) is
>> currently unavailable and returns 404 Object not found!
>>
>> This is possibly due to them now being linked from the side frame and the
>> URL now being http://cran.r-project.org/web/views/
>>
>> Unfortunately the webmaster link on the 404 page does not have any email
>> address associated with it, hence the mail to the list.
>>
>> Neil
>>   
>> 
> The proper address is likely [EMAIL PROTECTED], although the relevant 
> people will probably see it here as well.
>
> But what gave you the idea that 
> http://cran.r-project.org/web/packages/Views/ should work? Google seems 
> not to know it.
>   
Ah, it is the link used on http://cran.r-project.org/web/packages/ 

-- 
   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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] CRAN Taskviews returns 404

2008-02-18 Thread Peter Dalgaard
Neil Shephard wrote:
> The CRAN Taskviews (http://cran.r-project.org/web/packages/Views/) is
> currently unavailable and returns 404 Object not found!
>
> This is possibly due to them now being linked from the side frame and the
> URL now being http://cran.r-project.org/web/views/
>
> Unfortunately the webmaster link on the 404 page does not have any email
> address associated with it, hence the mail to the list.
>
> Neil
>   
The proper address is likely [EMAIL PROTECTED], although the relevant 
people will probably see it here as well.

But what gave you the idea that 
http://cran.r-project.org/web/packages/Views/ should work? Google seems 
not to know it.

-- 
   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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] working with weekdays

2008-02-18 Thread Gavin Simpson

On Mon, 2008-02-18 at 12:56 +, lamack lamack wrote:
> 
> Dear all, is it possible create, automatically, a variable with all the
> Monday and Thursday?
> 
> 18/02/2008
> 21/02/2008
> 25/02/2008
> 28/02/2008
> 03/03/2008
> .
> .
> for all months

Here is one way (for all months in 2008), if you know the first
Monday/Thursday

> first.mon <- as.Date("07-01-2008", format = "%d-%m-%Y")
> first.thu <- as.Date("03-01-2008", format = "%d-%m-%Y")
> mondays <- seq(from = first.mon, by = "weeks", length = 52)
> thursdays <- seq(from = first.thu, by = "weeks", length = 52)
> sort(c(mondays, thursdays))
  [1] "2008-01-03" "2008-01-07" "2008-01-10" "2008-01-14" "2008-01-17"
  [6] "2008-01-21".

And one if you don't:

> days <- seq(as.Date("01-01-2008", format = "%d-%m-%Y"), 
  as.Date("31-12-2008", format = "%d-%m-%Y"), by = "days") 
> monthu <- days[format(days, format = "%w") %in% c(1,4)]
> all.equal(monthu, sort(c(mondays, thursdays)))
[1] TRUE

HTH

G

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] working with weekdays

2008-02-18 Thread Henrique Dallazuanna
Try this:

values <- seq(as.Date("2008/01/01"), as.Date("2008/12/31"), by="days")
values[format(values, "%w") %in% c(1,4)]

On 18/02/2008, lamack lamack <[EMAIL PROTECTED]> wrote:
>
>
> Dear all, is it possible create, automatically, a variable with all the
> Monday and Thursday?
>
> 18/02/2008
> 21/02/2008
> 25/02/2008
> 28/02/2008
> 03/03/2008
> .
> .
> for all months
>
> Best regards
>
> JL
> _
> Confira vídeos com notícias do NY Times, gols direto do Lance, 
> videocas[[elided Hotmail spam]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to make a vector/list/array of POSIXlt object?

2008-02-18 Thread Bo Zhou
Ta. I will give that code a bash.

Could you explain why my code didn't work? 





> Date: Mon, 18 Feb 2008 00:25:44 -0500
> From: [EMAIL PROTECTED]
> To: [EMAIL PROTECTED]
> Subject: Re: [R] How to make a vector/list/array of POSIXlt object?
> CC: r-help@r-project.org
> 
> If the problem is that you have a vector of dates, a vector of times
> and a vector of data and you want to create a data frame with one
> POSIXct column and one column of data then try this:
> 
> dd <- c("01/22/2008", "02/13/2008")
> tt <- c("01:01:00", "23:01:12")
> dat <- 1:2
> 
> data.frame(dt = strptime(paste(dd, tt), "%m/%d/%Y %H:%M:%S"), dat)
> 
> # if you don't need subsecond data or time zones you could use chron
> 
> library(chron)
> data.frame(dt = chron(dd, tt), dat)
> 
> If this is intended to be a time series you might want to look at the zoo
> package.  It has three vignettes that give more info.
> 
> 
> On Feb 17, 2008 11:54 PM, Bo Zhou <[EMAIL PROTECTED]> wrote:
> > Hi Gabor,
> >
> > I'm using this code but it doesn't work for me
> >
> > > strptime2<-function (date,time) as.POSIXct(strptime(paste(date,time),
> > "%m/%d/%Y %H:%M:%S"))
> > > dt=mapply(strptime2, "01/01/2008", "00:00:00", SIMPLIFY=FALSE,
> > USE.NAMES=FALSE)
> > > df=data.frame(X1=dt,X2=1)
> > > dt
> > [[1]]
> > [1] "2008-01-01 Eastern Standard Time"
> >
> > > df
> >   structure.1199163600..class...c..POSIXtPOSIXcttzone.. X2
> > 12008-01-01  1
> >
> > Here df looks very wrong to me.
> >
> > So I tested this code:
> >
> > > df2=data.frame(X1=as.POSIXct(Sys.time()),X2=1)
> > > df2
> >X1 X2
> > 1 2008-02-17 23:43:08  1
> > > class(df2$X1)
> > [1] "POSIXt"  "POSIXct"
> >
> > Ah this worked as I expected.
> >
> > So some tweaking here - SIMPLIFY is set TRUE now:
> >
> > > strptime2<-function (date,time) as.POSIXct(strptime(paste(date,time),
> > "%m/%d/%Y %H:%M:%S"))
> > > dt=mapply(strptime2, "01/01/2008", "00:00:00", SIMPLIFY=TRUE,
> > USE.NAMES=FALSE)
> > > df=data.frame(X1=dt,X2=1)
> > > df
> >   X1 X2
> > 1 1199163600  1
> > > class(df$X1)
> > [1] "numeric"
> >
> > Hmm... it worked, but not in a way I wanted. The class info is missing.
> >
> > So how to get the result like this below? I do need that mapply +
> > strptime(paste), cos my CSV file is formatted in that way!
> >
> > > df2
> >X1 X2
> > 1 2008-02-17 23:43:08  1
> > > class(df2$X1)
> > [1] "POSIXt"  "POSIXct"
> >
> >
> > Any insight?
> >
> > Cheers,
> >
> > Bo
> >
> >
> >
> > > Date: Sun, 17 Feb 2008 15:53:28 -0500
> > > From: [EMAIL PROTECTED]
> > > To: [EMAIL PROTECTED]
> > > Subject: Re: [R] How to make a vector/list/array of POSIXlt object?
> > > CC: r-help@r-project.org
> >
> >
> > >
> > > Normally one uses POSIXct rather than POSIXlt for storage. See R News 4/1
> > for
> > > more info on date and time classes.
> > >
> > > On Feb 17, 2008 3:45 PM, Bo Zhou <[EMAIL PROTECTED]> wrote:
> > > >
> > > > Hi Guys,
> > > >
> > > > I'm cooking up my time series code. I want a data frame with first
> > column as timestamp in POSIXlt format.
> > > >
> > > > I hit on this the problem of how to create an array/list/vector of
> > POSIXlt objects. Code is as follows
> > > >
> > > >
> > > >
> > > > > dtt=array(dim = 2)
> > > > > t=as.POSIXlt( strptime("07/12/07 13:20:01", "%m/%d/%Y
> > %H:%M:%S",tz="GMT"))
> > > > > dtt
> > > > [1] NA NA
> > > > > t
> > > > [1] "0007-07-12 13:20:01 GMT"
> > > > > dtt[1]=t
> > > > Warning message:
> > > > In dtt[1] = t :
> > > > number of items to replace is not a multiple of replacement length
> > > > > class(dtt)
> > > > [1] "list"
> > > > > class(t)
> > > > [1] "POSIXt" "POSIXlt"
> > > > > unclass(t)
> > > > $sec
> > > > [1] 1
> > > >
> > > > $min
> > > > [1] 20
> > > >
> > > > $hour
> > > > [1] 13
> > > >
> > > > $mday
> > > > [1] 12
> > > >
> > > > $mon
> > > > [1] 6
> > > >
> > > > $year
> > > > [1] -1893
> > > >
> > > > $wday
> > > > [1] 4
> > > >
> > > > $yday
> > > > [1] 192
> > > >
> > > > $isdst
> > > > [1] 0
> > > >
> > > > attr(,"tzone")
> > > > [1] "GMT"
> > > >
> > > >
> > > >
> > > > Seems like POSIXlt is matrix in this case.
> > > >
> > > > Any suggestions?
> > > >
> > > > Cheers,
> > > >
> > > > B
> > > > _
> > > > [[elided Hotmail spam]]
> > > >
> > > > [[alternative HTML version deleted]]
> > > >
> > > > __
> > > > R-help@r-project.org mailing list
> > > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > > > and provide commented, minimal, self-contained, reproducible code.
> > > >
> >
> > 
> > Need to know the score, the latest news, or you need your Hotmail(R)-get 
> > your
> > "fix". Check it out.
> >

_


08
[[alternative HTML version deleted]]

Re: [R] Help wit barplot

2008-02-18 Thread Marc Schwartz
Tom.O wrote:
> Hi
>
> I have som problems with a barplot. It can be created easily in Excel but I
> cant manage to get it right in R.
>
> For example:
> MyData<-
> matrix(c(1,2,-1,-1,0,-2,-2,1,1,1,-2,3),ncol=3,dimnames=list(LETTERS[1:4],seq(3)))
>
> I want the barplot to stack the positive and negative values separatly so
> that positive values are stacked over 0 and negative is stacked under 0.
>
> I have added a jpg of the graph in excel
> http://www.nabble.com/file/p15542914/Excel.jpg Excel.jpg
>
> Regards Tom

Tom,

See this post from last year:

   http://finzi.psych.upenn.edu/R/Rhelp02a/archive/100161.html

HTH,

Marc Schwartz

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] fitted.values from zeroinfl (pscl package)

2008-02-18 Thread Achim Zeileis
On Mon, 18 Feb 2008, Sarah J Thomas wrote:

> Hello all:
>
> I have a question regarding the fitted.values returned from the
> zeroinfl() function. The values seem to be nearly identical to those
> fitted.values returned by the ordinary glm(). Why is this, shouldn't
> they be more "zero-inflated"?
>
> I construct a zero-inflated series of counts, called Y, like so:

To make this reproducible, I set the random seed to

set.seed(123)

in advance and then ran your source code

b= as.vector(c(1.5, -2))
g= as.vector(c(-3, 1))
x <- runif(100) # x is the covariate
X <- cbind(1,x)

p <- exp(X%*%g)/(1+exp(X%*%g))
m <- exp(X%*%b)   # log-link for the mean process
  # of the Poisson
Y <- rep(0, 100)

u <- runif(100)
for(i in 1:100) {
if( u[i] < p[i] ) { Y[i] = 0 }
else { Y[i] <- rpois(1, m[i]) }
}

# now let's compare the fitted.values from zeroinfl()
# and from glm()

z1 <- glm(Y ~ x, family=poisson)
z2 <- zeroinfl(Y ~ x|x) #poisson is the default

[snip]

> You can see that they are almost identical... and the fitted.values from
> zeroinfl don't seem to be zero-inflated at all! What is going on?

Well, let's see how zero inflated your observations are:

R> sum(u < p)
[1] 2

Wow, two (!) observations that have been zero-inflated. Let's see how much
the probability for observing a zero would have been

R> dpois(0, m[u < p])
[1] 0.3147816 0.1409670

which is not so low, in particular for the first one.

Overall, you've got

R> sum(Y < 1)
[1] 23

zeros in that data set and the expected number of zeros in a Poisson GLM
is

R> sum(dpois(0, fitted(z1)))
[1] 23.35615

So you have observed *less* zeros than expected by a Poisson GLM. Surely,
this is not the kind of data that zero-inflated models have been developed
for.

> Ultimately I want these fitted.values for a goodness of fit type of test
> to see if the zeroinfl model is needed or not for a given data series.
> With these fitted.values as they are, I am rejecting assumption of a
> zero-inflated model even when the data really are zero-inflated.

Maybe you ought to think about useful data-generating processes first
before designing tests or criticizing software...
Z

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Windows-GUI: "Install Packages from local zip files" and dependencies

2008-02-18 Thread Duncan Murdoch
On 18/02/2008 7:51 AM, Johannes Graumann wrote:
> Duncan Murdoch wrote:
> 
>>> 2) If that's not the default: is there a way to make it so?
>> Not simply, but of course it's possible with some work.  The problem is
>> that with repos = NULL, R doesn't know where to look for dependencies.
>> So you need to make two passes:  First, install the package, and second,
>> install its dependencies.
>>
>> You can use code like this to find the dependencies:
>>
>> library(tools)
>>
>> allpkgs <- installed.packages()
>> deps <- package.dependencies(allpkgs["test",])[[1]][,1]
>>
>> to find the dependencies of package test after you've installed it, and
>>
>> setdiff(deps, rownames(allpkgs))
>>
>> to find the ones that are not installed.  Where to find them to install
>> is the hard part:  are they also local, or on CRAN, or where?
> CRAN/BioC (both of which are set to Default == TRUE in the repository conf
> file.
> 
> This doesn't sound like something I can automate easily to appease the
> userbase of my little toolbox ...
> Is there a way to have a package upon installation run a post-install script
> with the functions you outline above?

We don't have post-install scripts, and you can't put this in the load 
hook, because that won't run without the dependencies.  So I think you 
need to offer something like biocLite, outside of your package.  Or set 
up a repository, and have your users install from there.

Duncan Murdoch

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] CRAN Taskviews returns 404

2008-02-18 Thread Neil Shephard

The CRAN Taskviews (http://cran.r-project.org/web/packages/Views/) is
currently unavailable and returns 404 Object not found!

This is possibly due to them now being linked from the side frame and the
URL now being http://cran.r-project.org/web/views/

Unfortunately the webmaster link on the 404 page does not have any email
address associated with it, hence the mail to the list.

Neil
-- 
View this message in context: 
http://www.nabble.com/CRAN-Taskviews-returns-404-tp15545185p15545185.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Number of digits of a value for problem 7.31 RFAQ

2008-02-18 Thread Richard . Cotton
> > What I mean is if R shows 2.3456 I want to obtain the info that 
digits=4 
> 
> > even if in facts the value has more (internal) digits.
> 
> Try:
> x = 1.23456789
> format(x, nsmall=20)
> # [1] "1.2345678899989009"

I've just re-read the question.  I suspect what you really wanted was 
something like:

getndp <- function(x, tol=2*.Machine$double.eps)
{
  ndp <- 0
  while(!isTRUE(all.equal(x, round(x, ndp), tol=tol))) ndp <- ndp+1 
  if(ndp > -log10(tol)) warning("Tolerance reached, ndp possibly 
underestimated.")
  ndp 
}

x = 0.123456
getndp(x)
# [1] 6

x2 <- 3
getndp(x2)
# [1] 0 

x3 <- 0.1234567890123456789
getndp(x3)
# [1] 16
# Warning message:
# In getndp(x3) : Tolerance reached, ndp possibly underestimated.

Regards,
Richie.

Mathematical Sciences Unit
HSL



ATTENTION:

This message contains privileged and confidential inform...{{dropped:20}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Number of digits of a value for problem 7.31 RFAQ

2008-02-18 Thread Richard . Cotton
> I need it to solve the problem (see RFAQ 7.31)that 0.2==0.2+0.1-0.1 
FALSE
> The solution suggested in RFAQ is to use 
isTRUE(all.equal(0.2,0.2+0.1-0.1))
> 
> But if I want to compare inequality:
> 0.2<=0.2 +0.1-0.1 TRUE
> but 0.2<=0.2 +0.1-0.1 FALSE
> bad!
> but in this case all.equal does not work I think... Unless to write a 
> double condition like
> 0.2>0.2 +0.1-0.1 | isTRUE(all.equal(0.2,0.2 +0.1-0.1))
> 
> The solution I found is to round the values, because 
> 0.2==round(0.2+0.1-0.1,2) TRUE
> However, one has to know the number of digits of the target value. How 
> can you do when it is unknown?

You can do this by using all.equal to check for cases close to the 
boundary:

if(isTRUE(all.equal(0.2, 0.2 + 0.1 - 0.1))) 
{
  message("Possible rounding problems.") 
  #investigate further
} else if(0.2 <= 0.2 + 0.1 - 0.1)
{
  #your code here
}

all.equal has a tolerance parameter that you can set to see how strict you 
want the equality to be; you may want to make this value smaller.

> What I mean is if R shows 2.3456 I want to obtain the info that digits=4 

> even if in facts the value has more (internal) digits.

Try:
x = 1.23456789
format(x, nsmall=20)
# [1] "1.2345678899989009"

Regards,
Richie.

Mathematical Sciences Unit
HSL




ATTENTION:

This message contains privileged and confidential inform...{{dropped:20}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] working with weekdays

2008-02-18 Thread lamack lamack


Dear all, is it possible create, automatically, a variable with all the
Monday and Thursday?

18/02/2008
21/02/2008
25/02/2008
28/02/2008
03/03/2008
.
.
for all months

Best regards

JL
_
Confira vídeos com notícias do NY Times, gols direto do Lance, videocas[[elided 
Hotmail spam]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Windows-GUI: "Install Packages from local zip files" and dependencies

2008-02-18 Thread Johannes Graumann
Duncan Murdoch wrote:

>> 2) If that's not the default: is there a way to make it so?
> 
> Not simply, but of course it's possible with some work.  The problem is
> that with repos = NULL, R doesn't know where to look for dependencies.
> So you need to make two passes:  First, install the package, and second,
> install its dependencies.
> 
> You can use code like this to find the dependencies:
> 
> library(tools)
> 
> allpkgs <- installed.packages()
> deps <- package.dependencies(allpkgs["test",])[[1]][,1]
> 
> to find the dependencies of package test after you've installed it, and
> 
> setdiff(deps, rownames(allpkgs))
> 
> to find the ones that are not installed.  Where to find them to install
> is the hard part:  are they also local, or on CRAN, or where?
CRAN/BioC (both of which are set to Default == TRUE in the repository conf
file.

This doesn't sound like something I can automate easily to appease the
userbase of my little toolbox ...
Is there a way to have a package upon installation run a post-install script
with the functions you outline above?

Thanks for any hints, Joh

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Windows-GUI: "Install Packages from local zip files" and dependencies

2008-02-18 Thread Prof Brian Ripley
On Mon, 18 Feb 2008, Johannes Graumann wrote:

> Hi,
>
> When using the "Packages --> Install packages from local zip files" menu
> item in the windows-gui:
> 1) is that supposed to automatically pull in dependencies (in that case I
> have to fix something in my package).

No (and ?install.packages tells you so).

> 2) If that's not the default: is there a way to make it so?

No.  There is no information available on dependencies unless you set up a 
local repository.  If you want the benfits of a repository, just use a 
repository.

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Windows-GUI: "Install Packages from local zip files" and dependencies

2008-02-18 Thread Duncan Murdoch
On 18/02/2008 7:14 AM, Johannes Graumann wrote:
> Hi,
> 
> When using the "Packages --> Install packages from local zip files" menu
> item in the windows-gui:
> 1) is that supposed to automatically pull in dependencies (in that case I
> have to fix something in my package).

No, it doesn't do that.  You can see the code it runs by printing

utils:::menuInstallLocal

It calls install.packages with repos = NULL.

> 2) If that's not the default: is there a way to make it so?

Not simply, but of course it's possible with some work.  The problem is 
that with repos = NULL, R doesn't know where to look for dependencies. 
So you need to make two passes:  First, install the package, and second, 
install its dependencies.

You can use code like this to find the dependencies:

library(tools)

allpkgs <- installed.packages()
deps <- package.dependencies(allpkgs["test",])[[1]][,1]

to find the dependencies of package test after you've installed it, and

setdiff(deps, rownames(allpkgs))

to find the ones that are not installed.  Where to find them to install 
is the hard part:  are they also local, or on CRAN, or where?

Duncan Murdoch

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Fw: r package

2008-02-18 Thread Alexandre Lerch Franco
help

- Original Message - 
From: Vito Ricci 
To: Alexandre Lerch Franco 
Sent: Monday, February 18, 2008 8:59 AM
Subject: Re: r package


Please send a message to R-help mailing list.
Regards.
Vito Ricci
 
Se non ora, quando?
Se non qui, dove?
Se non tu, chi? 

Personal Web Space: http://vr71.spaces.live.com/ 



- Messaggio originale -
Da: Alexandre Lerch Franco <[EMAIL PROTECTED]>
A: [EMAIL PROTECTED]
Inviato: Venerdì 15 febbraio 2008, 14:00:31
Oggetto: r package


Hi Vito:

my name is Alexandre Franco and I´m a doctoral student in finance and I´m 
looking for a package to run the JEP - Join Estimation Procedure - proposed by 
Chen and Liu (1993).

If you know the package, please send to me the directions, or name of the 
package.

Best regards.

Alexandre Franco







L'email della prossima generazione? Puoi averla con la nuova Yahoo! Mail



Esta mensagem foi verificada pelo E-mail Protegido Terra.
Scan engine: McAfee VirusScan / Atualizado em 15/02/2008 / Versão: 5.2.00/5231
Proteja o seu e-mail Terra: http://mail.terra.com.br/ 

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Solved (??) Behaviour of integrate (was 'Poisson-lognormal probab ility calculations')

2008-02-18 Thread k . jewell
Hi Again,

I think I've solved my problem, but please tell me if you think I'm wrong,
or you can see a better way!

A plot of the integrand showed a very sharp peak, so I was running into the
integrand "feature" mentioned in the note. I resolved it by limiting the
range of integration as shown here:
--
function (x, meanlog = 0, sdlog = 1, sdlim=6, rel.tol =
.Machine$double.eps^0.5,...) {
#+
# K. Jewell Feb 2008
# Based on VGAM::dpolono, modified to cover wider range of x
# Some simplification, but major change is to avoid integration problems by
#  limiting range of integration to sdlim (default = 6) standard deviations
each side 
#  of mean of Poisson or lognormal. Also increase default precision of
integration.
#  For x in 0:170 (working range of VGAM::dpolono) and default parameters,
max deviation 
#  from VGAM::dpolono is 3.243314e-05 of VGAM::dpolono result
#-
require(stats)
mapply(function(x, meanlog, sdlog, ...){
 lower <- min(max(0, x-sdlim*sqrt(x)), exp(meanlog-sdlim*sdlog))
 upper <- max(x+sdlim*sqrt(x), exp(meanlog+sdlim*sdlog))
 integrate(function(t, x, meanlog, sdlog) 
   exp(dpois(x,t,TRUE) + dlnorm(t, meanlog, sdlog, TRUE)),
  lower = lower, upper = upper, x = x, meanlog = meanlog, sdlog
= sdlog, 
  rel.tol = rel.tol, ...)$value
   },
   x, meanlog, sdlog, ...
   )
   }
---

Best regards,

Keith Jewell
mailto:[EMAIL PROTECTED] telephone (direct) +44 (0)1386 842055

> -Original Message-
> From: Jewell, Keith
> Sent: 15 February 2008 16:48
> To: 'r-help@r-project.org'
> Subject: Behaviour of integrate (was 'Poisson-lognormal probability
> calculations')
> 
> Hi again,
> 
> Adding further information to my own query, this function gets to the
> core of the problem, which I think lies in the behaviour of
> 'integrate'.
> -
> function (x, meanlog = 0, sdlog = 1, ...) {
> require(stats)
> integrand <- function(t, x, meanlog, sdlog) dpois(x,t)*dlnorm(t,
> meanlog, sdlog)
> mapply(function(x, meanlog, sdlog, ...)
> #(1/gamma(x+1))*
>  integrate(function(t, x, meanlog, sdlog)
> #   gamma(x+1)*
>   integrand(t, x, meanlog, sdlog),
>   lower = 0, upper = Inf, x = x, meanlog = meanlog, sdlog =
> sdlog, ...)$value,
>x, meanlog, sdlog, ...
>)
>}
> 
> Mathematically, the presence or not of the two commented lines should
> make no difference; they multiply the integrand by a constant (with
> respect to the integration), then divide the result by the same
> constant. In practice they make a big difference! I guess they're
> altering the behaviour of the 'integrate'.
> 
> I'd have thought the presence of the lines would worsen the behaviour.
> Without the lines the integrand is reasonably small, the integral is <
> 1. With the lines the limit on the integral is x!,  leading to "non-
> finite function values" for x much > 170, even if we use logs to get
> around the limit on gamma(x).
> 
> In fact with the lines the plot of function(x) v. x looks reasonable
> (but I don't know if the values are correct!!), but without the lines
> it looks silly, I just don't believe it!
> 
> I thought the problem might relate to the note in ?integrate "If the
> function is approximately constant (in particular, zero) over nearly
> all its range it is possible that the result and error estimate may be
> seriously wrong.". I wouldn't really expect multiplication by a large
> constant to fix such errors (??), but the lognormal distribution is
> skew so it might be considered "approximately constant ... over nearly
> all its range". Even
> though ...
> > integrate(dlnorm, 0, Inf)
> 1 with absolute error < 2.5e-07
> ... suggests that this is not the source of the problem,  I tried
> changing variables to integrate over a normally distributed variable:
> 
> function (x, meanlog = 0, sdlog = 1, ...) {
> require(stats)
> integrand <- function(t, x, meanlog, sdlog)
> dpois(x,exp(t))*dnorm(t, meanlog, sdlog)
> mapply(function(x, meanlog, sdlog, ...)
>  #  (1/gamma(x+1))*
>  integrate(function(t, x, meanlog, sdlog)
>  #  gamma(x+1)*
>   integrand(t, x, meanlog, sdlog),
>   lower = -Inf, upper = Inf, x = x, meanlog = meanlog,
> sdlog = sdlog, ...)$value,
>x, meanlog, sdlog, ...
>)
>}
> -
> Still no better; with the constants the values look reasonably smooth
> (but are they correct??), without the constants the values are silly.
> 
> I've tried reducing rel.tol and in

[R] Number of digits of a value for problem 7.31 RFAQ

2008-02-18 Thread Matthieu Stigler
Hello dear R users!

I did not find a function which gives information about the number of 
digits of a value shown by R.
Do you know one?

I need it to solve the problem (see RFAQ 7.31)that 0.2==0.2+0.1-0.1 FALSE
The solution suggested in RFAQ is to use isTRUE(all.equal(0.2,0.2+0.1-0.1))

But if I want to compare inequality:
0.2<=0.2 +0.1-0.1 TRUE
but 0.2<=0.2 +0.1-0.1 FALSE
bad!
but in this case all.equal does not work I think... Unless to write a 
double condition like
0.2>0.2 +0.1-0.1 | isTRUE(all.equal(0.2,0.2 +0.1-0.1))

The solution I found is to round the values, because 
0.2==round(0.2+0.1-0.1,2) TRUE
However, one has to know the number of digits of the target value. How 
can you do when it is unknown?

What I mean is if R shows 2.3456 I want to obtain the info that digits=4 
even if in facts the value has more (internal) digits.

I did not find a way to represent the number of digits but I found that 
as.numeric(as.character(x)) rounds the values to the number of digits 
shown. But is there a better way to make it?

Thanks!

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Windows-GUI: "Install Packages from local zip files" and dependencies

2008-02-18 Thread Johannes Graumann
Hi,

When using the "Packages --> Install packages from local zip files" menu
item in the windows-gui:
1) is that supposed to automatically pull in dependencies (in that case I
have to fix something in my package).
2) If that's not the default: is there a way to make it so?

Thanks, Joh

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] can we include nonparametric component for survival regression?

2008-02-18 Thread gallon li
i am trying to fit a survival regression model (cox model or parametric
model) in R by including the covariate effects as a function m(x) instead of
just beta*x. is it possible to fit such a model? can someone recommend some
reference? I searched but only found a package called addreg where
the hazard is actually modeled additively. That is not what i want.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Plotting a legend outside the figure - units

2008-02-18 Thread Duncan Murdoch
[EMAIL PROTECTED] wrote:
> Dear R gurus,
>
> I am trying to plot a legend in the margins of a figure. The basic idea is to 
> have two (or more) plots on the same figure, but then to have a common legend 
> at the bottom of the plot. The approach that I have taken is to setup the 
> figures, then do a "dummy" plot of the legend to determine the height of the 
> legend box (ie from the value returned by the legend() command). I then use 
> this height to adjust the outer margin of the plot to make sufficient space 
> for the legend. After plotting the actual plots, I then finally plot the 
> legend, moving it down to the bottom using the inset option in legend().
>
> The problem is a question of units. I have not been able to find what units 
> legend() returns anywhere in the documentation (ie in rect$h) - I've made the 
> asumption that its inches, but that may not necessarily be correct. 
> Similarly, the inset option to legend() requires units equal to a fraction of 
> the plot size, but I'm not exactly sure how to make that line up with what I 
> get out of legend().
>   
The rect values are returned in data coordinates, not inches.  So if you 
plot(1, 100), you'll get a legend with x coordinate near 1, and y 
coordinate near 100.  You specify the inset values in fractions of 
the plot region.
par("usr") can give the coordinates of the plot region in data 
coordinates if you want to translate.

Duncan Murdoch
> Suggestions and insights most welcome! Code follows at the bottom.
>
> Cheers,
>
> Mark
> -
>
> graphics.off()
> #First set number of plots
> par(mfrow=c(2,1))
>
> #Generate the legend once to get the height
> leg <-  legend("bottom",legend=c("Item A","Item B"),horiz=TRUE,pch=c(1,2),
>   col=c("red","blue"),plot=FALSE,trace=TRUE)
>
> par(oma=c(5,4,4,2)+0.1,mar=c(0,0,0,0),xpd=NA)
>
> #Add space at the bottom for the legend
> bottom.spacing  <-  par("omi") + c(leg$rect$h,0,0,0)
> par(omi=bottom.spacing)
>
> x <- runif(100)
> y <- runif(100)
> z <- rnorm(100)
> plot(x,y,xaxt="n",xlim=c(0,1),xlab="",pch=1,col="red")
> plot(y,z,xlim=c(0,1),pch=2,col="blue")
>
> #Actually plot legend this time - inset value is done manually by trial
> #and error, but obviosuly I'd rather have this done intelligently
> legend("bottom",legend=c("Item A","Item B"),horiz=TRUE,pch=c(1,2),
> col=c("red","blue"),inset=c(0,-0.475))
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Help wit barplot

2008-02-18 Thread Tom.O

Hi

I have som problems with a barplot. It can be created easily in Excel but I
cant manage to get it right in R.

For example:
MyData <-
matrix(c(1,2,-1,-1,0,-2,-2,1,1,1,-2,3),ncol=3,dimnames=list(LETTERS[1:4],seq(3)))

I want the barplot to stack the positive and negative values separatly so
that positive values are stacked over 0 and negative is stacked under 0. 

I have added a jpg of the graph in excel 
http://www.nabble.com/file/p15542914/Excel.jpg Excel.jpg 

Regards Tom
-- 
View this message in context: 
http://www.nabble.com/Help-wit-barplot-tp15542914p15542914.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] question on function arguments

2008-02-18 Thread baptiste Auguié
Hi,

I have two small issues with my R code, no big deal but curiosity  
really. Here is a sample code,


>
> x <- rnorm(1:10)
>
> foo <- function(a = 1, b = list(x = c(1:10), y = c(1:10))){
>   
>   for (ii in seq(along=b$y)){
>   
>   print(x[ii] + b$x[ii])
>   }
>   
>   
> }
>
> foo() # default OK
>
> foo(b=list(x=1, y=c(1:10))) # only the first term works
>
> foo(b$y = 1:5) # error
>

In the last call, i wish to use all default arguments but b$y (that  
is, using the default a and b$x). Is this possible?

The second call is more related to indices: i would like my argument b 
$x to be either a vector (default), or a scalar. In the latter, the  
loop b$x[ii] breaks when i would like it to recycle the single value.  
I can check the length of b$x with a "if" statement, but it becomes  
intricate when several arguments have this option of being vector or  
scalar. Is there a way to use b$x that would work in both cases?


Many thanks,

baptiste

_

Baptiste Auguié

Physics Department
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK

Phone: +44 1392 264187

http://newton.ex.ac.uk/research/emag
http://projects.ex.ac.uk/atto

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Plotting a legend outside the figure - units

2008-02-18 Thread rhelp . 20 . trevva
Dear R gurus,

I am trying to plot a legend in the margins of a figure. The basic idea is to 
have two (or more) plots on the same figure, but then to have a common legend 
at the bottom of the plot. The approach that I have taken is to setup the 
figures, then do a "dummy" plot of the legend to determine the height of the 
legend box (ie from the value returned by the legend() command). I then use 
this height to adjust the outer margin of the plot to make sufficient space for 
the legend. After plotting the actual plots, I then finally plot the legend, 
moving it down to the bottom using the inset option in legend().

The problem is a question of units. I have not been able to find what units 
legend() returns anywhere in the documentation (ie in rect$h) - I've made the 
asumption that its inches, but that may not necessarily be correct. Similarly, 
the inset option to legend() requires units equal to a fraction of the plot 
size, but I'm not exactly sure how to make that line up with what I get out of 
legend().

Suggestions and insights most welcome! Code follows at the bottom.

Cheers,

Mark
-

graphics.off()
#First set number of plots
par(mfrow=c(2,1))

#Generate the legend once to get the height
leg <-  legend("bottom",legend=c("Item A","Item B"),horiz=TRUE,pch=c(1,2),
  col=c("red","blue"),plot=FALSE,trace=TRUE)

par(oma=c(5,4,4,2)+0.1,mar=c(0,0,0,0),xpd=NA)

#Add space at the bottom for the legend
bottom.spacing  <-  par("omi") + c(leg$rect$h,0,0,0)
par(omi=bottom.spacing)

x <- runif(100)
y <- runif(100)
z <- rnorm(100)
plot(x,y,xaxt="n",xlim=c(0,1),xlab="",pch=1,col="red")
plot(y,z,xlim=c(0,1),pch=2,col="blue")

#Actually plot legend this time - inset value is done manually by trial
#and error, but obviosuly I'd rather have this done intelligently
legend("bottom",legend=c("Item A","Item B"),horiz=TRUE,pch=c(1,2),
col=c("red","blue"),inset=c(0,-0.475))

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Location of boxplots in bwplot

2008-02-18 Thread Sharon Kuhlmann-B

Hi,

I am trying to adjust the position of the boxplots when using bwplot(). In 
boxplot() there is the alternative to modify the position of the boxes by using 
boxplot(y~x, at=...).

However, I don't seem to find a comparable option in bwplot(). I would like in 
some instances to add more space between the boxplots, e.g. that I could 
specify at=c(1, 2, 5, 6) so that the first two boxplots are close to each 
other, and the other two a bit further apart. I can do something like this with 
the labels, but haven't managed with the actual boxplots. 

Grateful for any suggestions.

Sincerely,

Sharon



Sharon Kühlmann Berenzon, Ph.D.
Statistician
Dept. Epidemiology
Swedish Institute for Infectious Disease Control (SMI) 
 
[EMAIL PROTECTED]
tel. +46-8-457 2376; fax. +46-8-30 06 26

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] library(convert)

2008-02-18 Thread Prof Brian Ripley
On Mon, 18 Feb 2008, Gabor Csardi wrote:

> .. (Counting to ten.)
>
> The package is called 'convert'. It seems that this package is not
> on CRAN, however. I think you should ask the person whose code you're
> running.

It's on Bioconductor.  So select the BioC software repository on the 
Packages menu, then install.packages("convert") (or from the menu).
Or ask on the Bioconductor list.

> Gabor
>
> On Mon, Feb 18, 2008 at 10:17:41AM +0100, Schmitt, Corinna wrote:
>>
>> Hallo,
>>
>> I am running R-2.6. on Windows. I have a code which uses
>> library(convert). Can anyone tell me which package I need to install to
>> run this code. Everytime I receive the error message library (convert)
>> not found.
>>
>> Thanks, Corinna
>>
>>
>>
>>
>>
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> --
> Csardi Gabor <[EMAIL PROTECTED]>UNIL DGM
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] library(convert)

2008-02-18 Thread Gabor Csardi
Update: it is a Bioconductor package, so you need to do:

source("http://bioconductor.org/biocLite.R";)
biocLite("convert")

Gabor

On Mon, Feb 18, 2008 at 10:17:41AM +0100, Schmitt, Corinna wrote:
> 
> Hallo,
> 
> I am running R-2.6. on Windows. I have a code which uses
> library(convert). Can anyone tell me which package I need to install to
> run this code. Everytime I receive the error message library (convert)
> not found.
> 
> Thanks, Corinna
> 
> 
> 
> 
>  
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Csardi Gabor <[EMAIL PROTECTED]>UNIL DGM

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Weird SEs with effect()

2008-02-18 Thread Prof Brian Ripley
On Mon, 18 Feb 2008, Gustaf Granath wrote:

> Dear John & Brian,
>
> Ok. Now I start to understand. So basically  I cannot use the given SEs for 
> my purposes. They only make sense on the scale of log-counts. However in a 
> paper, the log-count scale is not very informative (the reader want to see 
> the effect on the scale you measure). If I understand you right, the 
> confidence intervals are fine though and maybe I will use them to illustrate 
> the reliability of the estimate. The problem is that showing SEs of the 
> adjusted means has become sort of the standard way to illustrate things in my 
> field (too many SAS users ?) and I might run into trouble with the reviewers. 
> I have several data sets with similar data and my plan was to use effect(). 
> That is why I want to figure this out. I hope I haven't been too annoying ;).
>
> Finally, is there a way to get correct SEs on the count scale (with adjusted 
> means)??
> I guess not, judging by your answers.

Yes, there is, at least approximately.  If X has mean mu and se s,
exp(X) has approximate mean exp(mu) and se s*exp(mu).  As in

> X <- rnorm(1000, 10, 0.1)
> sd(X)
[1] 0.09811928
> sd(exp(X))
[1] 2172.124
> exp(10)
[1] 22026.47

You will find slightly more accurate formulae on the help page ?plnorm
(they are exact if you assume normality of the estimated effects).


>
> Thanks again for your help,
>
> Gustaf
>
>
> John Fox wrote:
>> Dear Gustaf,
>>
>> 
>>> -Original Message-
>>> From: Gustaf Granath [mailto:[EMAIL PROTECTED]
>>> Sent: February-17-08 4:18 PM
>>> To: John Fox
>>> Cc: 'Prof Brian Ripley'; r-help@r-project.org
>>> Subject: RE: [R] Weird SEs with effect()
>>> 
>>> Dear John and Brian,
>>> Thank you for your help. I get the feeling that it is something
>>> fundamental that I do not understand here. Furthermore, a day of
>>> reading did not really help so maybe we have reached a dead end here.
>>> Nevertheless, here comes one last try.
>>> 
>>> I thought that the values produced by effect() were logs (e.g. in
>>> $fit). And then they were transformed (antilogged) with summary(). Was
>>> I wrong?
>>> 
>> 
>> I'm sorry that you're continuing to have problems with this.
>> 
>> Yes, there is a point that you don't understand: The SEs are on the scale 
>> of
>> the log-counts, but you can't get correct SEs on the scale of the counts by
>> exponentiating the SEs on the scale of the log-counts. What summary(), 
>> etc.,
>> do (and you can do) to produce confidence intervals on the count scale is
>> first to compute the intervals on the log-count scale and then to transform
>> the end-points.
>> 
>> I'm afraid that I can't make the point more clearly than that.
>> 
>> I hope this helps,
>>  John
>>
>> 
>>> What I want:
>>> I am trying to make a barplot with adjusted means with SEs (error
>>> bars), with the y axis labeled on the response scale.
>>> 
>>> #One of my GLM models (inf.level & def.level=factors, initial.size =
>>> covariate) #used as an example.
>>> #I was not able to make a reproducible example though. Sorry.
>>> 
>>> model <-
>>> glm(tot.fruit~initial.size+inf.level+def.level,family=quasipoisson)
>>> summary(model)
>>> Coefficients:
>>>  Estimate Std. Error t value Pr(>|t|)
>>> (Intercept)1.9368528  0.1057948  18.308  < 2e-16 ***
>>> initial.size   0.0015245  0.0001134  13.443  < 2e-16 ***
>>> inf.level50   -0.3142688  0.0908063  -3.461 0.000612 ***
>>> def.level12.5 -0.2329221  0.1236992  -1.883 0.060620 .
>>> def.level25   -0.1722354  0.1181993  -1.457 0.146062
>>> def.level50   -0.3543826  0.1212906  -2.922 0.003731 **
>>> 
>>> (Dispersion parameter for quasipoisson family taken to be 6.431139)
>>>  Null deviance: 2951.5  on 322  degrees of freedom
>>> Residual deviance: 1917.2  on 317  degrees of freedom
>>> 
>>> library(effects)
>>> def <- effect("def.level",model,se=TRUE)
>>> summary(def)
>>> $effect
>>> def.level
>>>  0  12.52550
>>> 11.145382  8.829541  9.381970  7.819672
>>> $lower
>>> def.level
>>> 0 12.5   25   50
>>> 9.495220 7.334297 7.867209 6.467627
>>> $upper
>>> def.level
>>> 0 12.5   25   50
>>> 13.08232 10.62962 11.18838  9.45436
>>> #Confidence intervals makes sense and are in line with the glm model
>>> result. Now #lets look at the standard errors. Btw, why aren't they
>>> given with summary?
>>> def$se
>>> 324325326327
>>> 0.08144281 0.09430438 0.08949864 0.09648573
>>> # As you can see, the SEs are very very very small.
>>> #In a graph it would look weird in combination with the glm result.
>>> #I thought that these values were logs. Thats why I used exp() which
>>> seems to be wrong.
>>> 
>>> Regards,
>>> 
>>> Gustaf
>>> 
>>>
>>> 
 Quoting John Fox <[EMAIL PROTECTED]>:
 Dear Brian and Gustaf,
 
 I too have a bit of trouble following what Gustaf is doing, but I
 
>>> think that
>>> 
 Brian's interpretation -- that Gustaf is trying to transform the
>

Re: [R] library(convert)

2008-02-18 Thread Gabor Csardi
.. (Counting to ten.)

The package is called 'convert'. It seems that this package is not 
on CRAN, however. I think you should ask the person whose code you're 
running.

Gabor

On Mon, Feb 18, 2008 at 10:17:41AM +0100, Schmitt, Corinna wrote:
> 
> Hallo,
> 
> I am running R-2.6. on Windows. I have a code which uses
> library(convert). Can anyone tell me which package I need to install to
> run this code. Everytime I receive the error message library (convert)
> not found.
> 
> Thanks, Corinna
> 
> 
> 
> 
>  
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Csardi Gabor <[EMAIL PROTECTED]>UNIL DGM

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] lineplot.CI problem

2008-02-18 Thread Dieter Vanderelst
Hi,

Thank you very much for this very rapid and helpful reply.

I'm giving the code a spin around the block. It looks to be working fine.

There is actually also something else on my lineplot.CI wish list. If I might 
be so bold to ask:

When plotting multiple data traces with relative high standard deviations, the 
whiskers tend to overlap. This hampers the interpretation of the data somewhat. 
Is there a way to let the function plot data points with some horizontal 
displacement to prevent this?

Regards,

Dieter

Manuel Morales wrote:
> Here's an updated version of lineplot.CI that will succeed even for
> cases where data are not present in all factor combinations. Also, this
> version has the option x.cont to specify that the x axis represents a
> continuous variable with proportional spacing. A new version of sciplot
> with these changes will be posted soon.
> 
> ## Examples:
> source("lineplot.CI.R")
> 
> ## Generate data
> time=c(rep(c(21:30),3),rep(c(1:10),3))
> y <- time+rnorm(60,0,1)
> factors <- rep(c(1:2),each=30)
> 
> ## Proportional spacing
> lineplot.CI(resp=y, x.factor=time, group=factors, x.cont=TRUE)
> 
> ## Factorial spacing
> lineplot.CI(resp=y, x.factor=time, group=factors)
> 
> Manuel
> 
> On Fri, 2008-02-15 at 15:18 +0100, Dieter Vanderelst wrote:
>> Hi List,
>>
>> I have a problem plotting data using the lineplot.CI command in the sciplot 
>> package.
>>
>> I want to plot the data of 2 experimental cases using different lines 
>> (traces). Time is on the X-axis. The tricky thing is that the data 
>> collection in the second case started later than for the first case. This is 
>> to say: the first n data points for the second case are missing.
>>
>> So far so good. However, when I plot the data using lineplot.CI, the 
>> standard error bars are not aligned correctly with the markers.
>>
>> I know that this might be difficult to imagine. Here you can find an 
>> example: http://i254.photobucket.com/albums/hh115/MarkerMe/example.png
>>
>> So, has anybody experienced this problem and solved it before? I think I 
>> could try padding the data of the second case with zeros to eliminate the 
>> missing data. But I hope there is a better solution.
>>
>> Regards,
>> Dieter
>> 
>> Dieter Vanderelst
>> dieter dot vanderelst at emailengine dot org
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] library(convert)

2008-02-18 Thread Schmitt, Corinna

Hallo,

I am running R-2.6. on Windows. I have a code which uses
library(convert). Can anyone tell me which package I need to install to
run this code. Everytime I receive the error message library (convert)
not found.

Thanks, Corinna




 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Possible overfitting of a GAM

2008-02-18 Thread Simon Wood
The figures don't obviously scream out `overfitting' to me, and the standard 
errors don't look excessively wide, given the data. Unless there is a strong 
reason for using `lo', you could also try the `gam' function in package 
`mgcv': it attempts to estimate the appropriate degree of smoothing 
automatically. If you get similar curves using mgcv::gam then you have some 
re-assurance that you don't have overfit here. 

On Saturday 16 February 2008 22:25, Thomas L Jones, PhD wrote:
> The subject is a Generalized Additive Model. Experts caution us against
> overfitting the data, which can cause inaccurate results. I am not a
> statistician (my background is in Computer Science). Perhaps some kind soul
> would take a look and vet the model for overfitting the data.
>
> The study estimated the ebb and flow of traffic through a voting place.
> Just one voting place was studied; the election was the U.S. mid-term
> election about a year ago. Procedure: The voting day was divided into
> five-minute bins, and the number of voters arriving in each bin was
> recorded. The voting day was 13 hours long, giving 156 bins.
>
> See http://tinyurl.com/36vzop for the scatterplot. There is a rather high
> random variation, due in part to the fact that the bin width was
> intentionally set to be narrow, in order to improve the amount of timing
> information gathered.
>
> http://tinyurl.com/3xjsyo displays the fitted curve. A GAM was used, with
> the loess smoothing algorithm (locally weighted regression). The default
> span was used. http://tinyurl.com/34av6l gives the scatterplot and the
> fitted curve. The two seem to match reasonably well.
>
> However, when I tried to generate the standard errors, things went awry.
> (Please see http://tinyurl.com/38ej2t ) There are three curves, seemingly
> the fitted curve and the curves for plus and minus two standard errors. The
> shapes seem okay, but there are large errors in the y values.
>
> Question: Have I overfitted the data?
>
> Feedback?
>
> Tom
> Thomas L. Jones, PhD, Computer Science
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html and provide commented, minimal,
> self-contained, reproducible code.

-- 
> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> +44 1225 386603  www.maths.bath.ac.uk/~sw283

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Weird SEs with effect()

2008-02-18 Thread Gustaf Granath
Dear John & Brian,

Ok. Now I start to understand. So basically  I cannot use the given SEs 
for my purposes. They only make sense on the scale of log-counts. 
However in a paper, the log-count scale is not very informative (the 
reader want to see the effect on the scale you measure). If I understand 
you right, the confidence intervals are fine though and maybe I will use 
them to illustrate the reliability of the estimate. The problem is that 
showing SEs of the adjusted means has become sort of the standard way to 
illustrate things in my field (too many SAS users ?) and I might run 
into trouble with the reviewers. I have several data sets with similar 
data and my plan was to use effect(). That is why I want to figure this 
out. I hope I haven't been too annoying ;).

Finally, is there a way to get correct SEs on the count scale (with 
adjusted means)??
I guess not, judging by your answers.

Thanks again for your help,

Gustaf


John Fox wrote:
> Dear Gustaf,
>
>   
>> -Original Message-
>> From: Gustaf Granath [mailto:[EMAIL PROTECTED]
>> Sent: February-17-08 4:18 PM
>> To: John Fox
>> Cc: 'Prof Brian Ripley'; r-help@r-project.org
>> Subject: RE: [R] Weird SEs with effect()
>>
>> Dear John and Brian,
>> Thank you for your help. I get the feeling that it is something
>> fundamental that I do not understand here. Furthermore, a day of
>> reading did not really help so maybe we have reached a dead end here.
>> Nevertheless, here comes one last try.
>>
>> I thought that the values produced by effect() were logs (e.g. in
>> $fit). And then they were transformed (antilogged) with summary(). Was
>> I wrong?
>> 
>
> I'm sorry that you're continuing to have problems with this.
>
> Yes, there is a point that you don't understand: The SEs are on the scale of
> the log-counts, but you can't get correct SEs on the scale of the counts by
> exponentiating the SEs on the scale of the log-counts. What summary(), etc.,
> do (and you can do) to produce confidence intervals on the count scale is
> first to compute the intervals on the log-count scale and then to transform
> the end-points.
>
> I'm afraid that I can't make the point more clearly than that.
>
> I hope this helps,
>  John
>
>   
>> What I want:
>> I am trying to make a barplot with adjusted means with SEs (error
>> bars), with the y axis labeled on the response scale.
>>
>> #One of my GLM models (inf.level & def.level=factors, initial.size =
>> covariate) #used as an example.
>> #I was not able to make a reproducible example though. Sorry.
>>
>> model <-
>> glm(tot.fruit~initial.size+inf.level+def.level,family=quasipoisson)
>> summary(model)
>> Coefficients:
>>  Estimate Std. Error t value Pr(>|t|)
>> (Intercept)1.9368528  0.1057948  18.308  < 2e-16 ***
>> initial.size   0.0015245  0.0001134  13.443  < 2e-16 ***
>> inf.level50   -0.3142688  0.0908063  -3.461 0.000612 ***
>> def.level12.5 -0.2329221  0.1236992  -1.883 0.060620 .
>> def.level25   -0.1722354  0.1181993  -1.457 0.146062
>> def.level50   -0.3543826  0.1212906  -2.922 0.003731 **
>>
>> (Dispersion parameter for quasipoisson family taken to be 6.431139)
>>  Null deviance: 2951.5  on 322  degrees of freedom
>> Residual deviance: 1917.2  on 317  degrees of freedom
>>
>> library(effects)
>> def <- effect("def.level",model,se=TRUE)
>> summary(def)
>> $effect
>> def.level
>>  0  12.52550
>> 11.145382  8.829541  9.381970  7.819672
>> $lower
>> def.level
>> 0 12.5   25   50
>> 9.495220 7.334297 7.867209 6.467627
>> $upper
>> def.level
>> 0 12.5   25   50
>> 13.08232 10.62962 11.18838  9.45436
>> #Confidence intervals makes sense and are in line with the glm model
>> result. Now #lets look at the standard errors. Btw, why aren't they
>> given with summary?
>> def$se
>> 324325326327
>> 0.08144281 0.09430438 0.08949864 0.09648573
>> # As you can see, the SEs are very very very small.
>> #In a graph it would look weird in combination with the glm result.
>> #I thought that these values were logs. Thats why I used exp() which
>> seems to be wrong.
>>
>> Regards,
>>
>> Gustaf
>>
>>
>> 
>>> Quoting John Fox <[EMAIL PROTECTED]>:
>>> Dear Brian and Gustaf,
>>>
>>> I too have a bit of trouble following what Gustaf is doing, but I
>>>   
>> think that
>> 
>>> Brian's interpretation -- that Gustaf is trying to transform the
>>>   
>> standard
>> 
>>> errors via the inverse link rather than transforming the ends of the
>>> confidence intervals -- is probably correct. If this is the case,
>>>   
>> then what
>> 
>>> Gustaf has done doesn't make sense.
>>>
>>> It is possible to get standard errors on the scale of the response
>>>   
>> (using,
>> 
>>> e.g., the delta method), but it's probably better to work on the
>>>   
>> scale of
>> 
>>> the linear predictor anyway. This is what the summary, print, and
>>>   
>> plot
>> 
>>> met

[R] Hazard model with long-term survivor (cure model)

2008-02-18 Thread 宋时歌
Dear All,

Are there R packages that can estimate survival model with long-term
survivors? This is sometimes known as "cure" model or "split-population"
model. Thanks.

Shige

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


  1   2   >