[R] combine words and dataframe

[Ravi S. Shankar]

Hi R,

 

I am trying to get an output like this

 

Hi

Hello

1   a   b

2   a   b

3   a   b

4   a   b

5   a   b

 

And write it as a text file

 

cat(paste("Hi",sep='\n',"Hello")) gives me

Hi

Hello

 

And when I try
cat(paste("Hi",sep='\n',"Hello"),sep='\n',data.frame(rep("a",3),rep("b",
3))) I get an error!

 

Thanks in advance for your help!

 

Regards

Ravi

 

 

 

 

This e-mail may contain confidential and/or privileged i...{{dropped:13}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Multiple Hexbinplots in 2 columns with a Single Categorical Variable

[Dr Stuart Reece]

Dear Ladies and Gentlemen,

 

I have a fairly large database (N=13,000) and a single main categorical
discriminator between the groups.

 

I want to look at the time course of a number of continuous biochemical
variables over chronologic age.

 

Therefore I believe I need to prepare hexbinplots in two columns with simple
regression lines in them (with useOuterStrips (in library(latticeExtra) if
possible) and also single hexbinplots with two regression lines
(panel.lmlines) corresponding to the two groups to facilitate data
comparison.

 

Tick marks and labels should be off in most panels, and the scale in each
row will be different.

 

I have some code written for the NHANES dataset in library(hexbin) which
almost does this job, but leaves blank paper up the middle of the two
columns of panels.

 

I have been exploring lattice for this which almost gets me there but not
quite.

 

The code I have written is as follows.

 

I would be ever so grateful for any advice anyone may be able to offer.

 

With best wishes,

 

Stuart Reece,

Australia.

 

 

useOuterStrips(hexbinplot(Transferin ~ Age | factor(Race) + factor(Sex),
data = NHANES, type = "r",

   aspect = 1, 

   scales = 

   list(x =

 list(relation = "free", rot = 0,

 at=list(TRUE, TRUE, NULL, NULL)),

y = 

 list(relation = "free", rot = 0,

 at=list(TRUE, NULL, TRUE, NULL))),

   par.settings =
   list(layout.heights = list(axis.panel = rep(c(1, 0),c(1, 1,
   par.strip.text = list(cex = 1)))  
 

 


[[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] Anova interaction not tested

[Gabriel Murray]

Playing around with this some more, it looks like it may purely be due to
the way that R handles unbalanced data. I will investigate lme for this.

Best,
Gabe

On Wed, Apr 8, 2009 at 3:24 PM, Gabriel Murray  wrote:

> Here is some toy data to give an example. Imagine someone is test-riding
> several bicycles to decide which to buy, and measures their average speed
> during each test-ride. They want to assess the main effects and interaction
> of brand and material on speed.
>
> SpeedSizeBrandMaterial
> 20 smallgiantsteel
> 25.1  mediumgiantsteel
> 25.2  medium  trekaluminum
> 26 smallcannondalecarbon
> 19.9  largegiantaluminum
> 21 mediumtrekcarbon
> 30 mediumcannondalecarbon
>
> This is the output I get:
>
> > fit = aov(Speed~Brand*Material, data=myData)
> > summary(fit)
> Df Sum Sq Mean Sq F value Pr(>F)
> Brand2 49.862  24.931  2.3738 0.2964
> Material 2 13.502   6.751  0.6428 0.6087
> Residuals2 21.005  10.503
>
>
> I suppose I am missing something very obvious here as to why the
> interaction is not reported, but I don't see what it is.
> Cheers,
>
> Gabe
>
>
>
> On Wed, Apr 8, 2009 at 1:34 PM, Patrick Connolly <
> p_conno...@slingshot.co.nz> wrote:
>
>> On Wed, 08-Apr-2009 at 12:59PM -0700, Gabriel Murray wrote:
>>
>> |> I've noticed with certain datasets that when I try to do an anova and
>> test
>> |> for main effects and interaction for two explanatory variables,
>> sometimes
>> |> the main effect results are given but not the interaction results. For
>> |> example,
>> |>
>> |> ex1 = aov(Score ~ var1*var2, data=myData)
>> |> summary(ex1)
>>
>> Could it be that there are no degrees of freedom left?
>>
>> That's about my best guess without enough information.
>>
>>
>> |> __
>> |> R-help@r-project.org mailing list
>> |> https://stat.ethz.ch/mailman/listinfo/r-help
>> |> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> |> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>> Nothing you sent is reproducible.
>>
>>
>> --
>> ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
>>___Patrick Connolly
>>  {~._.~}   Great minds discuss ideas
>>  _( Y )_ Average minds discuss events
>> (:_~*~_:)  Small minds discuss people
>>  (_)-(_)  . Eleanor Roosevelt
>>
>> ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
>>
>
>

[[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] Anova interaction not tested

[Gabriel Murray]

Here is some toy data to give an example. Imagine someone is test-riding
several bicycles to decide which to buy, and measures their average speed
during each test-ride. They want to assess the main effects and interaction
of brand and material on speed.

SpeedSizeBrandMaterial
20 smallgiantsteel
25.1  mediumgiantsteel
25.2  medium  trekaluminum
26 smallcannondalecarbon
19.9  largegiantaluminum
21 mediumtrekcarbon
30 mediumcannondalecarbon

This is the output I get:

> fit = aov(Speed~Brand*Material, data=myData)
> summary(fit)
Df Sum Sq Mean Sq F value Pr(>F)
Brand2 49.862  24.931  2.3738 0.2964
Material 2 13.502   6.751  0.6428 0.6087
Residuals2 21.005  10.503


I suppose I am missing something very obvious here as to why the interaction
is not reported, but I don't see what it is.
Cheers,

Gabe


On Wed, Apr 8, 2009 at 1:34 PM, Patrick Connolly  wrote:

> On Wed, 08-Apr-2009 at 12:59PM -0700, Gabriel Murray wrote:
>
> |> I've noticed with certain datasets that when I try to do an anova and
> test
> |> for main effects and interaction for two explanatory variables,
> sometimes
> |> the main effect results are given but not the interaction results. For
> |> example,
> |>
> |> ex1 = aov(Score ~ var1*var2, data=myData)
> |> summary(ex1)
>
> Could it be that there are no degrees of freedom left?
>
> That's about my best guess without enough information.
>
>
> |> __
> |> R-help@r-project.org mailing list
> |> https://stat.ethz.ch/mailman/listinfo/r-help
> |> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> |> and provide commented, minimal, self-contained, reproducible code.
>
>
> Nothing you sent is reproducible.
>
>
> --
> ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
>___Patrick Connolly
>  {~._.~}   Great minds discuss ideas
>  _( Y )_ Average minds discuss events
> (:_~*~_:)  Small minds discuss people
>  (_)-(_)  . Eleanor Roosevelt
>
> ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
>

[[alternative HTML version deleted]]

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

Re: [R] R File I/O Capability - Writing output to specific lines of existing file

[Jason Rupert]

Darn.

I was afraid of this.  Always kind of weak when it comes to file I/O 
approaches, so I guess I will have to stretch and try to put something 
together.   

Yeah.  The input flat text file is about 600 lines long.  I will be replacing 
about 200.  I already have the template for those 200 lines.  I guess I will 
just need to read down to line 400 and replace until line 600.  

Thanks again for all the insights regarding this. 



--- On Wed, 4/8/09, jim holtman  wrote:

> From: jim holtman 
> Subject: Re: [R] R File I/O Capability - Writing output to specific lines of  
> existing file
> To: jasonkrup...@yahoo.com
> Cc: R-help@r-project.org
> Date: Wednesday, April 8, 2009, 8:02 PM
> You can always read in the initialization file, make the
> updates to it
> and then write it back out.  If it is a text file, it would
> be very
> hard to write into the middle of it since there is no
> structure to the
> file.  You can read it in as a table (read.table) or just
> as lines
> (readLines) and the make any changes you want.  You can use
> regular
> expressions if you want to put something in the middle of
> one of the
> lines you have read in.
> 
> On Wed, Apr 8, 2009 at 10:13 AM, Jason Rupert
>  wrote:
> >
> > Currently I am using the R "write" command
> to output results to a *.txt file and then copying those
> results into an initialization file.  In an attempt to
> continue to automate the process I would like to have R
> write to the location in the existing initialization file,
> instead of me copying the data over.
> >
> > By any chance are there any R commands to help?
> >
> > Primarily, I will be using Windows, and the
> initialization file is a simple flat file, i.e. not XML or
> binary.
> >
> > I looked at:
> > (a)
> http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/R.io/html/00Index.html
> >
> > (b)
> http://search.r-project.org/cgi-bin/namazu.cgi?query=File+I%2FO&max=100&result=normal&sort=score&idxname=functions&idxname=Rhelp08
> >
> > But, this did not appear to provide the functionality
> to edit an existing file by adding information to the middle
> of the file.
> >
> > Thank you again for any help and insight.
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/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 that 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.

Re: [R] Help plotting image to window without margin

[Kingsford Jones]

See 'mar' under ?par

e.g.,

x <- y <-  -100:100
z <- outer(x, y, function(x,y) sqrt(x^2 + y^2))
par(mar=c(0,0,0,0))
image(x,y,z)

hth,

Kingsford Jones


On Wed, Apr 8, 2009 at 9:07 PM, Bob Meglen  wrote:
> I am using several scripts that employ various packages to process images 
> generated from multispectral data sets.
> I have used rimage and other packages to try to "plot" or "image" to place 
> images in a window. I have not been able to find a way to suppress the 
> margins that surround the image in the window. I would like to be able to 
> have the image fill the window without any white-space .
>
> I would appreciate some suggestions about how this can be done.
>
> Thanks,
> Bob Meglen
> Boulder, CO
>        [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

[R] Help plotting image to window without margin

[Bob Meglen]

I am using several scripts that employ various packages to process images 
generated from multispectral data sets.
I have used rimage and other packages to try to "plot" or "image" to place 
images in a window. I have not been able to find a way to suppress the margins 
that surround the image in the window. I would like to be able to have the 
image fill the window without any white-space .

I would appreciate some suggestions about how this can be done.

Thanks,
Bob Meglen
Boulder, CO
[[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 with biOps loading

[Bob Meglen]

I have tried several times to load biOps package after reading the posts in 
this archive regarding the necessity of placing these two libs in the PATH 
(libjpeg62.dll and libtiff3.dll).  I have tried locating the libs in several 
directories that should have worked, but I still get the following message.
--
 Error in inDL(x, as.logical(local), as.logical(now), ...) : 
 unable to load shared library 'C:/Program 
Files/R/R-2.8.1/library/biOps/libs/biOps.dll':
 LoadLibrary failure:  The specified module could not be found.

 Error: package/namespace load failed for 'biOps'
-
I am using R version 2.8.1 (2008-12-22) with Win-XP and Jaguar GUI

Can anyone give me a hint about where to place the libs and how to get the 
biOps package to work. It looks like the package could solve several of my 
image processing problems.

Thanks
[[alternative HTML version deleted]]

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

Re: [R] MLE for bimodal distribution

[Ravi Varadhan]

Dear Ted,

Thanks for your comments on the profilie-likelihood approach for ratio of 
sigmas.

I would like to look at your R code and the note on Aitken acceleration. I 
would appreciate if you could share this with me.

I am glad you nibbled at my "bait" on EM acceleration.  This is one of my 
favourite topics.  EM is indeed slow, but it pretty much guarantees convergence 
to a local maximum.  Acceleration methods, on the other hand, do not have this 
guarantee, as they forsake the ascent property.  This trade-off between rate of 
convergence and monotonicity is the crux of the acceleration problem.  I 
recently wrote a paper on this (Scand J Stats, June  2008).  

I have developed a class of acceleration methods, called SQUAREM, which strikes 
a good balance between speed and stability.  They are monotone, yet fast.  They 
will not be as fast as unbridled Aitken acceleration, which are susceptible to 
numerical instabilities.  SQUAREM, on the other hand,  is guarenteed to 
converge like the original EM algorithm, and will provide significant 
improvements in speed.  In other words, you can have your cake and eat it too!

I have written an R function to implement this.  I am planning to release an R 
package soon (as soon as I can free-up some time).  This package can be used to 
accelerate "any" EM algorithm.  The user is only required to supply the EM 
update function, i.e. the function that computes a single E and M step.  This 
function can be used as an "off-the-shelf" accelerator without needing any 
problem-specific input.  


Best,
Ravi.



Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvarad...@jhmi.edu


- Original Message -
From: ted.hard...@manchester.ac.uk (Ted Harding)
Date: Wednesday, April 8, 2009 7:43 pm
Subject: Re: [R] MLE for bimodal distribution
To: r-h...@stat.math.ethz.ch


> On 08-Apr-09 22:10:26, Ravi Varadhan wrote:
>  > EM algorithm is a better approach for maximum likelihood estimation
>  > of finite-mixture models than direct maximization of the mixture
>  > log-likelihood.  Due to its ascent properties, it is guaranteed to
>  > converge to a local maximum.  By theoretical construction, it also
>  > automatically ensures that all the parameter constraints are satisfied.
>  > 
>  > The problem with mixture estimation, however, is that the local maximum
>  > may not be a good solution.  Even global maximum may not be a good
>  > solution.  For example, it is well known that in a Gaussian mixture,
>  > you can get a log-likelihood of +infinity by letting mu = one of the
>  > data points, and sigma = 0 for one of the mixture components.  You 
> can
>  > detect trouble in MLE if you observe one of the following: (1) a
>  > degenerate solution, i.e. one of the mixing proportions is close to 
> 0,
>  > (2) one of the sigmas is too small.
>  > 
>  > EM convergence is sensitive to starting values.  Although, there are
>  > several starting strategies (see MacLachlan and Basford's book on
>  > Finite Mixtures), there is no universally sound strategy for 
> ensuring a
>  > good estimate (e.g. global maximum, when it makes sense). It is always
>  > a good practice to run EM with multiple starting values.  Be warned
>  > that EM convergence can be excruciatingly slow.  Acceleration methods
>  > can be of help in large simulation studies or for bootstrapping.
>  > 
>  > Best,
>  > Ravi.
>  
>  Well put! I've done a fair bit of mixture-fitting in my time,
>  and one approach which can be worth trying (and for which there
>  is often a reasonably objective justification) is to do a series
>  of fits (assuming a given number of components, e.g. 2 for the
>  present purpose) with the sigma's in a constant ratio in each one:
>  
>sigma2 = lambda*sigma1
>  
>  Then you won't get those singularities where it tries to climb
>  up a single data-point. If you do this for a range of values of
>  lambda, and keep track of the log-likelihood, then you get in
>  effect a profile likelihood for lambda. Once you see what that
>  looks like, you can then set about penalising ranges you don't
>  want to see.
>  
>  The reason I say "there is often a reasonably objective justification"
>  is that often you can be pretty sure, if there is a mixture present,
>  that lambda does not go outside say 1/10 - 10 (or whatever),
>  since you expect that the latent groups are fairly similar,
>  and unlikely to have markedly different sigma's.
>  
>  As to acceleration: agreed, EM can be slow. Aitken acceleration
>  can be dramatically faster. Several outlines of the Aitken procedure
>  can be found by googling on "aitken acceleration".
>  
>  I recently wrote a short note, describing its general principle
>  and outlining its application to the EM algorithm, using the Probit
>  model as illustration (with R code). For fitting th

Re: [R] R File I/O Capability - Writing output to specific lines of existing file

[jim holtman]

You can always read in the initialization file, make the updates to it
and then write it back out.  If it is a text file, it would be very
hard to write into the middle of it since there is no structure to the
file.  You can read it in as a table (read.table) or just as lines
(readLines) and the make any changes you want.  You can use regular
expressions if you want to put something in the middle of one of the
lines you have read in.

On Wed, Apr 8, 2009 at 10:13 AM, Jason Rupert  wrote:
>
> Currently I am using the R "write" command to output results to a *.txt file 
> and then copying those results into an initialization file.  In an attempt to 
> continue to automate the process I would like to have R write to the location 
> in the existing initialization file, instead of me copying the data over.
>
> By any chance are there any R commands to help?
>
> Primarily, I will be using Windows, and the initialization file is a simple 
> flat file, i.e. not XML or binary.
>
> I looked at:
> (a) 
> http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/R.io/html/00Index.html
>
> (b) 
> http://search.r-project.org/cgi-bin/namazu.cgi?query=File+I%2FO&max=100&result=normal&sort=score&idxname=functions&idxname=Rhelp08
>
> But, this did not appear to provide the functionality to edit an existing 
> file by adding information to the middle of the file.
>
> Thank you again for any help and insight.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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 that 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.

Re: [R] Need help in calculating studentized residuals/leverage values of non-linear model [nls()]

[Giam Xingli]

Hi Dieter,

Yes, I understand the definition of studentized residuals and I also know that 
studentizing the residuals is an easy exercise. However, I need the hat 
('leverage')-values to studentize the residuals. I tried to use hatvalues() to 
obtain the leverage values but while it works on glm and lm objects, it doesn't 
work on nls() objects.

I also tried to do this:

X <- model.matrix(model) # where model is an nls() object
H <- X %*% solve( t(X) %*% X ) %*% t(X)
hatvalues <- diag(H)

But model.matrix didn't work either. 

That is why I posted the question of getting the studentized residuals or the 
leverage values of nls() objects. I hope you don't misunderstand me, I'm sorry 
the original post came across as a 'homework' question.

Many thanks,
Xingli
(P.S. I saw your reply on CRAN. However I didn't get the daily digest, so I am 
replying to your previous post. Sorry about that!)  

--

Message: 72
Date: Mon, 6 Apr 2009 12:39:07 -0700
From: 
Subject: Re: [R] Need help in calculating studentized
residuals/leverage values   of non-linear model [nls()]
To: Dieter Menne ,
r-h...@stat.math.ethz.ch
Message-ID: <20090406153907.97vpo.402395.r...@mp16>
Content-Type: text/plain; charset=utf-8

Is the output of residuals() the studentized residuals or just the residuals?

 Dieter Menne  wrote: 
> Giam Xingli  nus.edu.sg> writes:
> 
> >I hope I can get advice regarding the calculation of leverage values or
> >studentized residual values of a non-linear regression model. It seems 
> > like
> >rstudent() does not work on a nls object.
> 
> residuals() should work for nls.
> 
> Dieter
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] MLE for bimodal distribution

[Ted Harding]

On 08-Apr-09 22:10:26, Ravi Varadhan wrote:
> EM algorithm is a better approach for maximum likelihood estimation
> of finite-mixture models than direct maximization of the mixture
> log-likelihood.  Due to its ascent properties, it is guaranteed to
> converge to a local maximum.  By theoretical construction, it also
> automatically ensures that all the parameter constraints are satisfied.
> 
> The problem with mixture estimation, however, is that the local maximum
> may not be a good solution.  Even global maximum may not be a good
> solution.  For example, it is well known that in a Gaussian mixture,
> you can get a log-likelihood of +infinity by letting mu = one of the
> data points, and sigma = 0 for one of the mixture components.  You can
> detect trouble in MLE if you observe one of the following: (1) a
> degenerate solution, i.e. one of the mixing proportions is close to 0,
> (2) one of the sigmas is too small.
> 
> EM convergence is sensitive to starting values.  Although, there are
> several starting strategies (see MacLachlan and Basford's book on
> Finite Mixtures), there is no universally sound strategy for ensuring a
> good estimate (e.g. global maximum, when it makes sense). It is always
> a good practice to run EM with multiple starting values.  Be warned
> that EM convergence can be excruciatingly slow.  Acceleration methods
> can be of help in large simulation studies or for bootstrapping.
> 
> Best,
> Ravi.

Well put! I've done a fair bit of mixture-fitting in my time,
and one approach which can be worth trying (and for which there
is often a reasonably objective justification) is to do a series
of fits (assuming a given number of components, e.g. 2 for the
present purpose) with the sigma's in a constant ratio in each one:

  sigma2 = lambda*sigma1

Then you won't get those singularities where it tries to climb
up a single data-point. If you do this for a range of values of
lambda, and keep track of the log-likelihood, then you get in
effect a profile likelihood for lambda. Once you see what that
looks like, you can then set about penalising ranges you don't
want to see.

The reason I say "there is often a reasonably objective justification"
is that often you can be pretty sure, if there is a mixture present,
that lambda does not go outside say 1/10 - 10 (or whatever),
since you expect that the latent groups are fairly similar,
and unlikely to have markedly different sigma's.

As to acceleration: agreed, EM can be slow. Aitken acceleration
can be dramatically faster. Several outlines of the Aitken procedure
can be found by googling on "aitken acceleration".

I recently wrote a short note, describing its general principle
and outlining its application to the EM algorithm, using the Probit
model as illustration (with R code). For fitting the location
parameter alone, Raw EM took 59 iterations, Aitken-accelerated EM
took 3. For fitting the location and scale paramaters, Raw EM took 108,
and Aitken took 4.

If anyone would like a copy (PS or PDF) of this, drop me a line.
Or is there some "repository" for R-help (like the "Files" area
in Google Groups) where one can place files for others to download?

[And, if not, do people think it might be a good idea?]

Best wishes tgo all,
Ted.


E-Mail: (Ted Harding) 
Fax-to-email: +44 (0)870 094 0861
Date: 09-Apr-09   Time: 00:33:02
-- XFMail --

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

Re: [R] Doubt about aov and lm function... bug?

[William Dunlap]

I'm also suspicious of code using parse to construct
expressions.  However this problem arose because the
terms class is a subclass of the formula class and
AlgDesign defines a model.matrix method for the formula
class but not for the terms class.  Hence
  model.matrix(terms(...))
invokes model.matrix.formula, which causes the problem
because AlgDesign::model.matrix.formula isn't general enough to
handle a terms object.

If AlgDesign had a model.matrix.terms function that
called model.matrix.default or was a copy of
model.matrix.default then the problem would go away.

Bill Dunlap
TIBCO Software Inc - Spotfire Division
wdunlap tibco.com 

>> library(AlgDesign)
> 
>> aov(Sepal.Length ~ Species, data=iris)
> 
> Error in parse(text = x) :
>   unexpected symbol in "Sepal(Sepal.Length+Species)Length"

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] MLE for bimodal distribution

[Ravi Varadhan]

EM algorithm is a better approach for maximum likelihood estimation of 
finite-mixture models than direct maximization of the mixture log-likelihood.  
Due to its ascent properties, it is guaranteed to converge to a local maximum.  
By theoretical construction, it also automatically ensures that all the 
parameter constraints are satisfied.  

The problem with mixture estimation, however, is that the local maximum may not 
be a good solution.  Even global maximum may not be a good solution.  For 
example, it is well known that in a Gaussian mixture, you can get a 
log-likelihood of +infinity by letting mu = one of the data points, and sigma = 
0 for one of the mixture components.  You can detect trouble in MLE if you 
observe one of the following: (1) a degenerate solution, i.e. one of the mixing 
proportions is close to 0, (2) one of the sigmas is too small.

EM convergence is sensitive to starting values.  Although, there are several 
starting strategies (see MacLachlan and Basford's book on Finite Mixtures), 
there is no universally sound strategy for ensuring a good estimate (e.g. 
global maximum, when it makes sense). It is always a good practice to run EM 
with multiple starting values.  Be warned that EM convergence can be 
excruciatingly slow.  Acceleration methods can be of help in large simulation 
studies or for bootstrapping.

Best,
Ravi.


Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvarad...@jhmi.edu


- Original Message -
From: Bert Gunter 
Date: Wednesday, April 8, 2009 4:14 pm
Subject: Re: [R] MLE for bimodal distribution
To: 'Ben Bolker' , r-help@r-project.org


> The problem is that fitting mixtures models is hard -- (non)identifiability
>  is a serious problem: very large sample sizes may be necessary to clearly
>  distinguish the modes. As V&R say in MASS, 4th edition, p. 320: " ...
>  fitting normal mixtures is a difficult problem, and the results 
> obtained are
>  often heavily dependent on the initial configuration ([i.e. starting 
> values.
>  BG] supplied"
>  
>  The EM algorithm is quite commonly used: you might have a look at 
> em() and
>  related functions in the mclust package if Ben's straightforward approach
>  fails to do the job for you. No guarantee em will work either, of course.
>  
>  -- Bert
>  
>  
>  Bert Gunter
>  Genentech Nonclinical Biostatistics
>  650-467-7374
>  
>  -Original Message-
>  From: r-help-boun...@r-project.org [ On
>  Behalf Of Ben Bolker
>  Sent: Wednesday, April 08, 2009 12:47 PM
>  To: r-help@r-project.org
>  Subject: Re: [R] MLE for bimodal distribution
>  
>  
>  
>  
>  _nico_ wrote:
>  > 
>  > Hello everyone,
>  > 
>  > I'm trying to use mle from package stats4 to fit a bi/multi-modal
>  > distribution to some data, but I have some problems with it.
>  > Here's what I'm doing (for a bimodal distribution):
>  > 
>  > # Build some fake binormally distributed data, the procedure fails 
> also
>  > with real data, so the problem isn't here
>  > data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))
>  > # Just to check it's bimodal
>  > plot(density(data))
>  > f = function(m, s, m2, s2, w)
>  > {
>  > -log( sum(w*dnorm(data, mean=m, sd=s)) + sum((1-w)*dnorm(data, mean=m2,
>  > sd=s2)) )
>  > }
>  > 
>  > res = mle(f, start=list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6))
>  > 
>  > This gives an error:
>  > Error in optim(start, f, method = method, hessian = TRUE, ...) : 
>  >   non-finite finite-difference value [2]
>  > In addition: There were 50 or more warnings (use warnings() to see 
> the
>  > first 50)
>  > And the warnings are about dnorm producing NaNs
>  > 
>  > So, my questions are:
>  > 1) What does "non-finite finite-difference value" mean? My guess 
> would be
>  > that an Inf was passed somewhere where a finite number was expected...?
>  > I'm not sure how optim works though...
>  > 2) Is the log likelihood function I wrote correct? Can I use the 
> same type
>  > of function for 3 modes?
>  > 3) What should I do to avoid function failure? I tried by changing 
> the
>  > parameters, but it did not work.
>  > 4) Can I put constraints to the parameters? In particular, I would 
> like w
>  > to be between 0 and 1.
>  > 5) Is there a way to get the log-likelihood value, so that I can compare
>  > different extimations?
>  > 6) Do you have any (possibly a bit "pratically oriented") readings 
> about
>  > MLE to suggest?
>  > 
>  > Thanks in advance
>  > Nico
>  > 
>  
>Here's some tweaked code that works.
>  Read about the "L-BFGS-B" method in ?optim to constrain parameters,
>  although you will have some difficulty making this work for more than
>  two components.  I think there's also a mixture model problem in
>  Venables & Ripley (MASS).
>  
>There are almost certainly more efficient approaches for mixture
>  model fitting, al

Re: [R] Minimum Spanning Tree

[Gábor Csárdi]

On Wed, Apr 8, 2009 at 10:20 PM, jpearl01  wrote:
>
> That's like a miracle!  The only thing that would make this graph perfect is
> if the lengths of the edges were in the same ratio as the actual edge
> lengths from the matrix.  Is it possible to alter that?

Not really. The thing is that the nodes are ordered into layers based
on the distance (in steps) from the root. This is how the
Reingold-Tilford layouting algorithm works. If you start changing the
lengths of the edges, then the nodes could move on top of each other.

You would need to modify the algorithm itself for this, I don't know
how difficult that would be.

What you can do is modifying the width of the edges, I know that is
not as good, but it is still something:

library(igraph)

tab <- read.csv("http://www.nabble.com/file/p22957493/sp_matrix.csv";)
tab <- tab[,-1]

g <- graph.adjacency(as.matrix(tab), weighted=TRUE)
V(g)$label <- V(g)$name

mst <- minimum.spanning.tree(g)

rescale <- function(x,from,to) {# linearly rescale a vector
  r <- range(x)
  (x-r[1]) / (r[2]-r[1]) * (to-from) + from
}

E(mst)$width <- rescale(E(mst)$weight,1,5)   # width between 1 and 5

lay <- layout.reingold.tilford(mst, root=which.max(degree(mst))-1)
lay <- cbind(lay[,2], lay[,1])# rotate
x11(width=15, height=8)
plot(mst, layout=lay, vertex.size=25, vertex.size2=10,
 vertex.shape="crectangle", asp=FALSE,
 vertex.label.cex=0.7, vertex.color="white", edge.arrow.mode=0)

G.

> Thank you!!
> ~josh
>
[...]

-- 
Gabor Csardi  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] Doubt about aov and lm function... bug?

[Peter Dalgaard]

Jose Claudio Faria wrote:


library(AlgDesign)



aov(Sepal.Length ~ Species, data=iris)


Error in parse(text = x) :
  unexpected symbol in "Sepal(Sepal.Length+Species)Length"

I wonder it is really a small bug. Do you agree?


Yes, in AlgDesign.

I'm always suspicious of things involving parse() - see also fortune(106).



--
   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
~~ - (p.dalga...@biostat.ku.dk)  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] Doubt about aov and lm function... bug?

[Rolf Turner]

On 9/04/2009, at 9:17 AM, Jose Claudio Faria wrote:


Hi Peter and Rolf,

Many thanks for all!

The "bug" was just found: it is related with the package AlgDesign!
(and I've been loading it from Rprofile.site already for a long  
long time). :-(





I wonder it is really a small bug. Do you agree?


It seems to me that there is a conflict between what lm() does and what
the functions in AlgDesign do.  Presumably they are used for completely
different purposes.  So this isn't a bug as such, it is simply a matter
of arranging things so that the right model.matrix() function gets  
called

in the right circumstances.

Sounds like a job for namespaces (about which I was *complaining*  
recently;

:-) ) although I don't understand these well enough to be sure that they
could help here.  Peter will know, I have every confidence!

cheers,

Rolf Turner

##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

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

Re: [R] does grid.rect() not accept NULL argument from gpar(col=NULL) ?

[Paul Murrell]

Hi


Mark Heckmann wrote:
> I want to draw a grid rectangle without a border.
> 
> ?gpar says:
> 
> "Specifying the value NULL for a parameter is the same as not specifying any
> value for that parameter, except for col and fill, where NULL indicates not
> to draw a border or not to fill an area (respectively)." 
> 
> pushViewport(viewport(height=unit(.8, "npc")))
>   grid.rect(gp=gpar(col=NULL, fill="green"))
> popViewport()
> 
> Still a border is drawn. What am I doing wrong?


Your mistake was to assume that the documentation is correct! :)

That behaviour is actually WAY out of date.  You can do what you want
using NA ...

pushViewport(viewport(height=unit(.8, "npc")))
grid.rect(gp=gpar(col=NA, fill="green"))
popViewport()

Thanks for reporting the problem (the docs are now fixed).

Paul


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

-- 
Dr Paul Murrell
Department of Statistics
The University of Auckland
Private Bag 92019
Auckland
New Zealand
64 9 3737599 x85392
p...@stat.auckland.ac.nz
http://www.stat.auckland.ac.nz/~paul/

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

Re: [R] Minimum Spanning Tree

[jpearl01]

That's like a miracle!  The only thing that would make this graph perfect is
if the lengths of the edges were in the same ratio as the actual edge
lengths from the matrix.  Is it possible to alter that?

Thank you!!
~josh



Actually,

library(igraph)

tab <- read.csv("http://www.nabble.com/file/p22957493/sp_matrix.csv";)
tab <- tab[,-1]

g <- graph.adjacency(as.matrix(tab), weighted=TRUE)
V(g)$label <- V(g)$name

mst <- as.undirected(minimum.spanning.tree(g))

lay <- layout.reingold.tilford(mst, root=which.max(degree(mst))-1)
lay <- cbind(lay[,2], lay[,1])# rotate
x11(width=15, height=8)
plot(mst, layout=lay, vertex.size=25, vertex.size2=10,
 vertex.shape="rectangle", asp=FALSE,
 vertex.label.cex=0.7, vertex.color="white")

works relatively well for me on your graph. It doesn't for you?

Best,
Gabor

-- 
View this message in context: 
http://www.nabble.com/Minimum-Spanning-Tree-tp22934813p22958820.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] Doubt about aov and lm function... bug?

[Jose Claudio Faria]

Hi Peter and Rolf,

Many thanks for all!

The "bug" was just found: it is related with the package AlgDesign!
(and I've been loading it from Rprofile.site already for a long long time). :-(

See below:

R version 2.8.1 Patched (2009-01-22 r47680)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

Design library by Frank E Harrell Jr

Type library(help='Design'), ?Overview, or ?Design.Overview')
to see overall documentation.

deldir 0.0-7

> aov(Sepal.Length ~ Species, data=iris)
Call:
   aov(formula = Sepal.Length ~ Species, data = iris)

Terms:
Species Residuals
Sum of Squares   6339
Deg. of Freedom   2   147

Residual standard error: 0.51
Estimated effects may be unbalanced

> library(AlgDesign)

> aov(Sepal.Length ~ Species, data=iris)

Error in parse(text = x) :
  unexpected symbol in "Sepal(Sepal.Length+Species)Length"

I wonder it is really a small bug. Do you agree?

All the best,
-- 
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\
Jose Claudio Faria
Estatistica - prof. Titular
UESC/DCET/Brasil
joseclaudio.fa...@gmail.com
joseclaudio.fa...@terra.com.br
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Doubt about aov and lm function... bug?

[Peter Dalgaard]

Jose Claudio Faria wrote:

# a) aov

av <- aov(Sepal.Length ~ Species, data=iris)


Error in parse(text = x) :
  unexpected symbol in "Sepal(Sepal.Length+Species)Length"


traceback()

13: parse(text = x)
12: eval(parse(text = x)[[1]])
11: formula(eval(parse(text = x)[[1]]))
10: formula.character(object, env = baseenv())
9: formula(object, env = baseenv())
8: as.formula(frml)
7: expand.formula(frml, colnames(data))
6: model.matrix.formula(mt, mf, contrasts)


Er, I have

>  getAnywhere(model.matrix.formula)
no object named 'model.matrix.formula' was found
> methods(model.matrix)
[1] model.matrix.default model.matrix.lm

Where did your model.matrix.formula come from???

-p



5: model.matrix(mt, mf, contrasts)
4: lm(formula = Sepal.Length ~ Species, data = iris, singular.ok = TRUE)
3: eval(expr, envir, enclos)
2: eval(lmcall, parent.frame())
1: aov(Sepal.Length ~ Species, data = iris)

# b) lm

lm1 <- lm(Sepal.Length ~ Sepal.Width, data=iris)


Error in parse(text = x) :
  unexpected symbol in "Sepal(Sepal.Length+Sepal.Width)Length"


traceback()

10: parse(text = x)
9: eval(parse(text = x)[[1]])
8: formula(eval(parse(text = x)[[1]]))
7: formula.character(object, env = baseenv())
6: formula(object, env = baseenv())
5: as.formula(frml)
4: expand.formula(frml, colnames(data))
3: model.matrix.formula(mt, mf, contrasts)
2: model.matrix(mt, mf, contrasts)
1: lm(Sepal.Length ~ Sepal.Width, data = iris)

HTH,



--
   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
~~ - (p.dalga...@biostat.ku.dk)  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] Doubt about aov and lm function... bug?

[Jose Claudio Faria]

# a) aov
> av <- aov(Sepal.Length ~ Species, data=iris)

Error in parse(text = x) :
  unexpected symbol in "Sepal(Sepal.Length+Species)Length"

> traceback()
13: parse(text = x)
12: eval(parse(text = x)[[1]])
11: formula(eval(parse(text = x)[[1]]))
10: formula.character(object, env = baseenv())
9: formula(object, env = baseenv())
8: as.formula(frml)
7: expand.formula(frml, colnames(data))
6: model.matrix.formula(mt, mf, contrasts)
5: model.matrix(mt, mf, contrasts)
4: lm(formula = Sepal.Length ~ Species, data = iris, singular.ok = TRUE)
3: eval(expr, envir, enclos)
2: eval(lmcall, parent.frame())
1: aov(Sepal.Length ~ Species, data = iris)

# b) lm
> lm1 <- lm(Sepal.Length ~ Sepal.Width, data=iris)

Error in parse(text = x) :
  unexpected symbol in "Sepal(Sepal.Length+Sepal.Width)Length"

> traceback()
10: parse(text = x)
9: eval(parse(text = x)[[1]])
8: formula(eval(parse(text = x)[[1]]))
7: formula.character(object, env = baseenv())
6: formula(object, env = baseenv())
5: as.formula(frml)
4: expand.formula(frml, colnames(data))
3: model.matrix.formula(mt, mf, contrasts)
2: model.matrix(mt, mf, contrasts)
1: lm(Sepal.Length ~ Sepal.Width, data = iris)

HTH,
-- 
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\
Jose Claudio Faria
Estatistica - prof. Titular
UESC/DCET/Brasil
joseclaudio.fa...@gmail.com
joseclaudio.fa...@terra.com.br
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Doubt about aov and lm function... bug?

[Peter Dalgaard]

Rolf Turner wrote:


There is no bug.

Both of your examples work fine for me.  


Ditto

You may have a corrupt

version of ``iris'' somewhere in your search path before ``datasets''.

What does find("iris") tell you?  What does names("iris") tell you?


Could also be something interfering with the model.frame/model.matrix 
bits. The error message suggests that something did a substitute of ".", 
but I can't offhand see a parse(text=x) construct in any of the usual 
suspects. A traceback() could be informative.


--
   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
~~ - (p.dalga...@biostat.ku.dk)  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] Doubt about aov and lm function... bug?

[Jose Claudio Faria]

I've been finding the same problem in all R versions I'm using after
2.6.0 (I think).

> find("iris")
[1] ".GlobalEnv"   "package:datasets"

It is really strange...
-- 
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\
Jose Claudio Faria
Estatistica - prof. Titular
UESC/DCET/Brasil
joseclaudio.fa...@gmail.com
joseclaudio.fa...@terra.com.br
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Anova interaction not tested

[Patrick Connolly]

On Wed, 08-Apr-2009 at 12:59PM -0700, Gabriel Murray wrote:

|> I've noticed with certain datasets that when I try to do an anova and test
|> for main effects and interaction for two explanatory variables, sometimes
|> the main effect results are given but not the interaction results. For
|> example,
|> 
|> ex1 = aov(Score ~ var1*var2, data=myData)
|> summary(ex1)

Could it be that there are no degrees of freedom left?

That's about my best guess without enough information.


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


Nothing you sent is reproducible.


-- 
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.   
   ___Patrick Connolly   
 {~._.~}   Great minds discuss ideas
 _( Y )_ Average minds discuss events 
(:_~*~_:)  Small minds discuss people  
 (_)-(_)  . Eleanor Roosevelt
  
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Doubt about aov and lm function... bug?

[Rolf Turner]

There is no bug.

Both of your examples work fine for me.  You may have a corrupt
version of ``iris'' somewhere in your search path before ``datasets''.

What does find("iris") tell you?  What does names("iris") tell you?

cheers,

Rolf Turner

On 9/04/2009, at 8:16 AM, Jose Claudio Faria wrote:


Hi,

The below very strange:

#  a) aov function
av <- aov(Sepal.Length ~ Species, data=iris)

#  Error in parse(text = x) :
#   unexpected symbol in "Sepal(Sepal.Length+Species)Length"

av <- aov(iris[, 1] ~ iris[, 5])
#  summary(av)
#   Df Sum Sq Mean Sq F value Pr(>F)
#  iris[, 5] 2   63.231.6 119 <2e-16 ***
#  Residuals   147   39.0 0.3
#  ---
#  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



#  b) lm function
lm1 <- lm(Sepal.Length ~ Sepal.Width, data=iris)

#  Error in parse(text = x) :
#   unexpected symbol in "Sepal(Sepal.Length+Sepal.Width)Length"

lm1 <- lm(iris[, 1] ~ iris[, 2])
summary(lm1)

#  Call:
#  lm(formula = iris[, 1] ~ iris[, 2])
#
#  Residuals:
# Min 1Q Median 3QMax
#  -1.556 -0.633 -0.112  0.558  2.223
#
#  Coefficients:
#  Estimate Std. Error t value Pr(>|t|)
#  (Intercept)6.526  0.479   13.63   <2e-16 ***
#  iris[, 2] -0.223  0.155   -1.44 0.15
#  ---
#  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
#  Residual standard error: 0.82 on 148 degrees of freedom
#  Multiple R-squared: 0.0138,  Adjusted R-squared: 0.00716
#  F-statistic: 2.07 on 1 and 148 DF,  p-value: 0.152

I'm wrong or is it a bug?

platform   i386-pc-mingw32
arch   i386
os mingw32
system i386, mingw32
status Patched
major  2
minor  8.1
year   2009
month  01
day22
svn rev47680
language   R
version.string R version 2.8.1 Patched (2009-01-22 r47680)

Thanks,
--  
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\

Jose Claudio Faria
Estatistica - prof. Titular
UESC/DCET/Brasil
joseclaudio.fa...@gmail.com
joseclaudio.fa...@terra.com.br
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting- 
guide.html

and provide commented, minimal, self-contained, reproducible code.



##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

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

Re: [R] MLE for bimodal distribution

[Rubén Roa-Ureta]

_nico_ wrote:

Hello everyone,

I'm trying to use mle from package stats4 to fit a bi/multi-modal
distribution to some data, but I have some problems with it.
Here's what I'm doing (for a bimodal distribution):

# Build some fake binormally distributed data, the procedure fails also with
real data, so the problem isn't here
data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))
# Just to check it's bimodal
plot(density(data))
f = function(m, s, m2, s2, w)
{
-log( sum(w*dnorm(data, mean=m, sd=s)) + sum((1-w)*dnorm(data, mean=m2,
sd=s2)) )
}

res = mle(f, start=list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6))

This gives an error:
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  non-finite finite-difference value [2]

In addition: There were 50 or more warnings (use warnings() to see the first
50)
And the warnings are about dnorm producing NaNs

So, my questions are:
1) What does "non-finite finite-difference value" mean? My guess would be
that an Inf was passed somewhere where a finite number was expected...? I'm
not sure how optim works though...
2) Is the log likelihood function I wrote correct? Can I use the same type
of function for 3 modes?
  

I think it is correct but it is difficult to fit as others have pointed out.
As an alternative, you may discretise your variable so the mixture is 
fit to counts on a histogram using a multinomial likelihood.

HTH
Rubén

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] MLE for bimodal distribution

[Ravi Varadhan]

Ben,

You can also do this nicely with the function spg() in the "BB" package

set.seed(1001)

data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))

f = function(m, s, m2, s2, w)
{
  -sum(log(w*dnorm(data, mean=m, sd=s)+
   (1-w)*dnorm(data, mean=m2, sd=s2)))
}


start0 <- list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6)

fcn <- function(x) f(x[1], x[2], x[3], x[4], x[5])  

spg(par=as.numeric(start0), fn=fcn, lower=c(-Inf, 0, -Inf, 0, 0), upper=c(Inf, 
Inf, Inf, Inf, 1))

# of course, "L-BFGS-B" also works well

optim(par=as.numeric(start0), fn=fcn, method="L-BFGS-B", lower=c(-Inf, 0, -Inf, 
0, 0), upper=c(Inf, Inf, Inf, Inf, 1))

Ravi.


Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvarad...@jhmi.edu


- Original Message -
From: Ben Bolker 
Date: Wednesday, April 8, 2009 3:49 pm
Subject: Re: [R] MLE for bimodal distribution
To: r-help@r-project.org


>  
>  
>  _nico_ wrote:
>  > 
>  > Hello everyone,
>  > 
>  > I'm trying to use mle from package stats4 to fit a bi/multi-modal
>  > distribution to some data, but I have some problems with it.
>  > Here's what I'm doing (for a bimodal distribution):
>  > 
>  > # Build some fake binormally distributed data, the procedure fails 
> also
>  > with real data, so the problem isn't here
>  > data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))
>  > # Just to check it's bimodal
>  > plot(density(data))
>  > f = function(m, s, m2, s2, w)
>  > {
>  > -log( sum(w*dnorm(data, mean=m, sd=s)) + sum((1-w)*dnorm(data, mean=m2,
>  > sd=s2)) )
>  > }
>  > 
>  > res = mle(f, start=list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6))
>  > 
>  > This gives an error:
>  > Error in optim(start, f, method = method, hessian = TRUE, ...) : 
>  >   non-finite finite-difference value [2]
>  > In addition: There were 50 or more warnings (use warnings() to see 
> the
>  > first 50)
>  > And the warnings are about dnorm producing NaNs
>  > 
>  > So, my questions are:
>  > 1) What does "non-finite finite-difference value" mean? My guess 
> would be
>  > that an Inf was passed somewhere where a finite number was expected...?
>  > I'm not sure how optim works though...
>  > 2) Is the log likelihood function I wrote correct? Can I use the 
> same type
>  > of function for 3 modes?
>  > 3) What should I do to avoid function failure? I tried by changing 
> the
>  > parameters, but it did not work.
>  > 4) Can I put constraints to the parameters? In particular, I would 
> like w
>  > to be between 0 and 1.
>  > 5) Is there a way to get the log-likelihood value, so that I can compare
>  > different extimations?
>  > 6) Do you have any (possibly a bit "pratically oriented") readings 
> about
>  > MLE to suggest?
>  > 
>  > Thanks in advance
>  > Nico
>  > 
>  
>Here's some tweaked code that works.
>  Read about the "L-BFGS-B" method in ?optim to constrain parameters,
>  although you will have some difficulty making this work for more than
>  two components.  I think there's also a mixture model problem in
>  Venables & Ripley (MASS).
>  
>There are almost certainly more efficient approaches for mixture
>  model fitting, although this "brute force" approach should be fine
>  if you don't need to do anything too complicated.
>  (RSiteSearch("mixture model"))
>  
>  set.seed(1001)
>  data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))
>  # Just to check it's bimodal
>  plot(density(data))
>  f = function(m, s, m2, s2, w)
>  {
>-sum(log(w*dnorm(data, mean=m, sd=s)+
> (1-w)*dnorm(data, mean=m2, sd=s2)))
>  }
>  
>  library(stats4)
>  start0 <- list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6)
>  do.call("f",start0) ## OK
>  res = mle(f, start=start0)
>  
>  with(as.list(coef(res)),
>   curve(w*dnorm(x,m,s)+(1-w)*dnorm(x,m2,s2),add=TRUE,col=2))
>  
>  
>  -- 
>  View this message in context: 
>  Sent from the R help mailing list archive at Nabble.com.
>  
>  __
>  R-help@r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  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] Constrained, multiple response statistics

[Jonathan Greenberg]

Hmm, it LOOKS like mvpart may be along the lines of what I want, but is 
mvpart a nominal classification tree, or can it handle multiple, 
continuous response variables as well?


--j

Gene Leynes wrote:
This sounds very similar to what I've been working on, but I'm not 
sure without an example.


My solution has been to use an optimization that normalizes inside the 
objective function.  The betas that are provided by optim are not 
normalized, however since they were normalized inside the objective 
function, normalizing them after the fact mirrors the internal 
workings of the objective function.


see this example:
http://markmail.org/message/ze5237m6gbgvvvyf

Still, after looking at several statistical packages, and considering 
the thoughtful responses from my post, I think that there must be a 
better way using existing models, so I've been looking at other 
packages / models.


On Tue, Apr 7, 2009 at 6:10 PM, Jonathan Greenberg 
mailto:greenb...@ucdavis.edu>> wrote:


R'ers:

  I was hoping I could get some direction on this.  I have a
dataset of the form:

Y1,Y2,...,YM = f(X1,X2,...,XN), where N is >>> M

The response data (Y1,Y2,...,YM) is frequency data, such that the
sum of all Yi = 1.0.  Both Xj and Yi are continuous variables.

I'm trying to figure out the best approach(es) to solving for the
model f() -- any ideas?  I could solve each Y one at a time, but
the lack of constraint worries me, and I'm pretty sure that
normalizing the data afterwards to sum to 1.0 is not going to work
out properly.  Thoughts?  I've never worked with multiple response
statistics before, so I'm mostly trying to get some pointers on
where to begin investigating...

--j

-- 


Jonathan A. Greenberg, PhD
Postdoctoral Scholar
Center for Spatial Technologies and Remote Sensing (CSTARS)
University of California, Davis
One Shields Avenue
The Barn, Room 250N
Davis, CA 95616
Cell: 415-794-5043
AIM: jgrn307, MSN: jgrn...@hotmail.com
, Gchat: jgrn307

__
R-help@r-project.org  mailing list
https://stat.ethz.ch/mailman/listinfo/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] Doubt about aov and lm function... bug?

[Jose Claudio Faria]

Hi,

The below very strange:

#  a) aov function
av <- aov(Sepal.Length ~ Species, data=iris)

#  Error in parse(text = x) :
#   unexpected symbol in "Sepal(Sepal.Length+Species)Length"

av <- aov(iris[, 1] ~ iris[, 5])
#  summary(av)
#   Df Sum Sq Mean Sq F value Pr(>F)
#  iris[, 5] 2   63.231.6 119 <2e-16 ***
#  Residuals   147   39.0 0.3
#  ---
#  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



#  b) lm function
lm1 <- lm(Sepal.Length ~ Sepal.Width, data=iris)

#  Error in parse(text = x) :
#   unexpected symbol in "Sepal(Sepal.Length+Sepal.Width)Length"

lm1 <- lm(iris[, 1] ~ iris[, 2])
summary(lm1)

#  Call:
#  lm(formula = iris[, 1] ~ iris[, 2])
#
#  Residuals:
# Min 1Q Median 3QMax
#  -1.556 -0.633 -0.112  0.558  2.223
#
#  Coefficients:
#  Estimate Std. Error t value Pr(>|t|)
#  (Intercept)6.526  0.479   13.63   <2e-16 ***
#  iris[, 2] -0.223  0.155   -1.44 0.15
#  ---
#  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
#  Residual standard error: 0.82 on 148 degrees of freedom
#  Multiple R-squared: 0.0138,  Adjusted R-squared: 0.00716
#  F-statistic: 2.07 on 1 and 148 DF,  p-value: 0.152

I'm wrong or is it a bug?

platform   i386-pc-mingw32
arch   i386
os mingw32
system i386, mingw32
status Patched
major  2
minor  8.1
year   2009
month  01
day22
svn rev47680
language   R
version.string R version 2.8.1 Patched (2009-01-22 r47680)

Thanks,
-- 
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\
Jose Claudio Faria
Estatistica - prof. Titular
UESC/DCET/Brasil
joseclaudio.fa...@gmail.com
joseclaudio.fa...@terra.com.br
///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] MLE for bimodal distribution

[_nico_]

Ben Bolker wrote:
> 
> 
>   Here's some tweaked code that works.
> [cut]
> 

Thanks, that saved me a few headaches. I also find out the answer to my
(dumb) question #5, which is obviously to call f with the returned
parameters or use the logLik function.

I will have a look at the mixture model problem, as you suggested.
Looking at my data I don't think I really need more than bimodal
distributions anyway, yet I was going to have a go at it mostly out of
curiosity.

Well, thank you very much again for your help! It was really appreciated.

Nico
-- 
View this message in context: 
http://www.nabble.com/MLE-for-bimodal-distribution-tp22954970p22958318.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] MLE for bimodal distribution

[Bert Gunter]

The problem is that fitting mixtures models is hard -- (non)identifiability
is a serious problem: very large sample sizes may be necessary to clearly
distinguish the modes. As V&R say in MASS, 4th edition, p. 320: " ...
fitting normal mixtures is a difficult problem, and the results obtained are
often heavily dependent on the initial configuration ([i.e. starting values.
BG] supplied"

The EM algorithm is quite commonly used: you might have a look at em() and
related functions in the mclust package if Ben's straightforward approach
fails to do the job for you. No guarantee em will work either, of course.

-- Bert


Bert Gunter
Genentech Nonclinical Biostatistics
650-467-7374

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Ben Bolker
Sent: Wednesday, April 08, 2009 12:47 PM
To: r-help@r-project.org
Subject: Re: [R] MLE for bimodal distribution




_nico_ wrote:
> 
> Hello everyone,
> 
> I'm trying to use mle from package stats4 to fit a bi/multi-modal
> distribution to some data, but I have some problems with it.
> Here's what I'm doing (for a bimodal distribution):
> 
> # Build some fake binormally distributed data, the procedure fails also
> with real data, so the problem isn't here
> data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))
> # Just to check it's bimodal
> plot(density(data))
> f = function(m, s, m2, s2, w)
> {
> -log( sum(w*dnorm(data, mean=m, sd=s)) + sum((1-w)*dnorm(data, mean=m2,
> sd=s2)) )
> }
> 
> res = mle(f, start=list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6))
> 
> This gives an error:
> Error in optim(start, f, method = method, hessian = TRUE, ...) : 
>   non-finite finite-difference value [2]
> In addition: There were 50 or more warnings (use warnings() to see the
> first 50)
> And the warnings are about dnorm producing NaNs
> 
> So, my questions are:
> 1) What does "non-finite finite-difference value" mean? My guess would be
> that an Inf was passed somewhere where a finite number was expected...?
> I'm not sure how optim works though...
> 2) Is the log likelihood function I wrote correct? Can I use the same type
> of function for 3 modes?
> 3) What should I do to avoid function failure? I tried by changing the
> parameters, but it did not work.
> 4) Can I put constraints to the parameters? In particular, I would like w
> to be between 0 and 1.
> 5) Is there a way to get the log-likelihood value, so that I can compare
> different extimations?
> 6) Do you have any (possibly a bit "pratically oriented") readings about
> MLE to suggest?
> 
> Thanks in advance
> Nico
> 

  Here's some tweaked code that works.
Read about the "L-BFGS-B" method in ?optim to constrain parameters,
although you will have some difficulty making this work for more than
two components.  I think there's also a mixture model problem in
Venables & Ripley (MASS).

  There are almost certainly more efficient approaches for mixture
model fitting, although this "brute force" approach should be fine
if you don't need to do anything too complicated.
(RSiteSearch("mixture model"))

set.seed(1001)
data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))
# Just to check it's bimodal
plot(density(data))
f = function(m, s, m2, s2, w)
{
  -sum(log(w*dnorm(data, mean=m, sd=s)+
   (1-w)*dnorm(data, mean=m2, sd=s2)))
}

library(stats4)
start0 <- list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6)
do.call("f",start0) ## OK
res = mle(f, start=start0)

with(as.list(coef(res)),
 curve(w*dnorm(x,m,s)+(1-w)*dnorm(x,m2,s2),add=TRUE,col=2))


-- 
View this message in context:
http://www.nabble.com/MLE-for-bimodal-distribution-tp22954970p22957613.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] xyplot - show values of a series on graph

[alienz747]

This is great, thanks a lot!

On Wed, Apr 8, 2009 at 3:50 PM, baptiste auguie  wrote:

> not very clean, but perhaps,
>
> xyplot(cars+trucks~year, data=df2, type="o",
>  panel=function(x,y,subscripts,...){
>panel.xyplot(x,y,subscripts=subscripts,...)
>
>  grid.text(unit(df2$year,"native"),unit(df2$cars,"native"),label=df2$cars,
> just="top")}
>   )
>
> baptiste
>
> On 8 Apr 2009, at 20:25, taz9 wrote:
>
>
>> Thank you very much for your help. I tried to use lattice but I'm not sure
>> how to restrict it to display only the values of "cars".
>>
>> xyplot(cars+trucks~year, data=df2, type="o",
>> panel=function(x,y,...){
>>panel.xyplot(x,y,...)
>>grid.text(unit(x,"native"),unit(y,"native"),label=y, just="top")}
>>  )
>>
>>
>>
>> baptiste auguie-2 wrote:
>>
>>>
>>> with ggplot2,
>>>
>>> d <- melt(df2,id="year")
>>> qplot(year,value,data=d,colour=variable,geom=c("line","point")) +
>>> geom_text(data= subset(d, variable == "cars"), aes(label=value))
>>>
>>>
>>> with lattice, my best guess would be to use grid.text in a custom
>>> panel function.
>>>
>>>
>>> Hope this helps,
>>>
>>> baptiste
>>>
>>> On 8 Apr 2009, at 19:40, taz9 wrote:
>>>
>>>
 Hi All,

 I have a very simple graph:

 cars <- c(1, 3, 6, 4, 9)
 trucks <- c(2, 5, 4, 5, 12)
 year <- c(2004, 2005, 2006, 2007, 2008)
 df2<-data.frame(cars,trucks,year)
 xyplot(cars+trucks~year, data=df2, type="o")

 I need to show the values of "cars" on the graph. How can I do this?

 Thanks.

 --
 View this message in context:

 http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22956986.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.

>>>
>>> _
>>>
>>> Baptiste Auguié
>>>
>>> School of Physics
>>> University of Exeter
>>> Stocker Road,
>>> Exeter, Devon,
>>> EX4 4QL, UK
>>>
>>> Phone: +44 1392 264187
>>>
>>> http://newton.ex.ac.uk/research/emag
>>>
>>> __
>>> R-help@r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>>
>> --
>> View this message in context:
>> http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22957501.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.
>>
>
> _
>
> Baptiste Auguié
>
> School of Physics
> University of Exeter
> Stocker Road,
> Exeter, Devon,
> EX4 4QL, UK
>
> Phone: +44 1392 264187
>
> http://newton.ex.ac.uk/research/emag
> __
>
>

[[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] xyplot - show values of a series on graph

[alienz747]

Thank you so much Deepayan! My problem is solved.

On Wed, Apr 8, 2009 at 3:47 PM, Deepayan Sarkar
wrote:

> On Wed, Apr 8, 2009 at 12:25 PM, taz9  wrote:
> >
> > Thank you very much for your help. I tried to use lattice but I'm not
> sure
> > how to restrict it to display only the values of "cars".
> >
> > xyplot(cars+trucks~year, data=df2, type="o",
> > panel=function(x,y,...){
> >panel.xyplot(x,y,...)
> >grid.text(unit(x,"native"),unit(y,"native"),label=y, just="top")}
> >   )
>
> A custom 'panel.groups' will let you condition on group number (see
> ?panel.superpose):
>
> xyplot(cars+trucks~year, data=df2, type="o",
>   panel = panel.superpose,
>   panel.groups = function(x, y, ..., group.number) {
>   panel.xyplot(x,y,...)
>   if (group.number == 1) {
>   require(grid)
>   grid.text(unit(x,"native"),
> unit(y,"native"),
> label=y, just="top")
>   }
>   })
>
> -Deepayan
>
> >
> >
> >
> > baptiste auguie-2 wrote:
> >>
> >> with ggplot2,
> >>
> >> d <- melt(df2,id="year")
> >> qplot(year,value,data=d,colour=variable,geom=c("line","point")) +
> >> geom_text(data= subset(d, variable == "cars"), aes(label=value))
> >>
> >>
> >> with lattice, my best guess would be to use grid.text in a custom
> >> panel function.
> >>
> >>
> >> Hope this helps,
> >>
> >> baptiste
> >>
> >> On 8 Apr 2009, at 19:40, taz9 wrote:
> >>
> >>>
> >>> Hi All,
> >>>
> >>> I have a very simple graph:
> >>>
> >>> cars <- c(1, 3, 6, 4, 9)
> >>> trucks <- c(2, 5, 4, 5, 12)
> >>> year <- c(2004, 2005, 2006, 2007, 2008)
> >>> df2<-data.frame(cars,trucks,year)
> >>> xyplot(cars+trucks~year, data=df2, type="o")
> >>>
> >>> I need to show the values of "cars" on the graph. How can I do this?
> >>>
> >>> Thanks.
> >>>
> >>> --
> >>> View this message in context:
> >>>
> http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22956986.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.
> >>
> >> _
> >>
> >> Baptiste Auguié
> >>
> >> School of Physics
> >> University of Exeter
> >> Stocker Road,
> >> Exeter, Devon,
> >> EX4 4QL, UK
> >>
> >> Phone: +44 1392 264187
> >>
> >> http://newton.ex.ac.uk/research/emag
>

[[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] Anova interaction not tested

[Gabriel Murray]

I've noticed with certain datasets that when I try to do an anova and test
for main effects and interaction for two explanatory variables, sometimes
the main effect results are given but not the interaction results. For
example,

ex1 = aov(Score ~ var1*var2, data=myData)
summary(ex1)

gives me only the main effects for var1 and var2, but not the interaction. I
also tried

ex1 = aov(Score ~ var1+var2+var1:var2, data=myData)
summary(ex1)

but it still only gives the main effects. The only way I can get the
interaction results is if I test for ONLY that, e.g.

ex1 = aov(Score ~ var1:var2, data=myData)
summary(ex1)

Why would the interaction not be tested in the first two examples?

Thanks,
Gabe

[[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] Display a very low p-value

[hadley wickham]

>  pnorm(37:39,lower.tail=FALSE)
> [1] 5.725571e-300  0.00e+00  0.00e+00
>
>  This is just a limitation of double precision floating-point arithmetic
> ...
>
>  curve(pnorm(x,lower.tail=FALSE),from=30,to=40,log="y")
> .Machine$double.xmin

But note

curve(pnorm(x,lower.tail=FALSE, log=T),from=30,to=1000)

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] xyplot - show values of a series on graph

[baptiste auguie]

not very clean, but perhaps,

xyplot(cars+trucks~year, data=df2, type="o",
 panel=function(x,y,subscripts,...){
panel.xyplot(x,y,subscripts=subscripts,...)
 	 
grid 
.text(unit(df2$year,"native"),unit(df2$cars,"native"),label=df2$cars,  
just="top")}

   )

baptiste

On 8 Apr 2009, at 20:25, taz9 wrote:



Thank you very much for your help. I tried to use lattice but I'm  
not sure

how to restrict it to display only the values of "cars".

xyplot(cars+trucks~year, data=df2, type="o",
panel=function(x,y,...){
panel.xyplot(x,y,...)
grid.text(unit(x,"native"),unit(y,"native"),label=y, just="top")}
  )



baptiste auguie-2 wrote:


with ggplot2,

d <- melt(df2,id="year")
qplot(year,value,data=d,colour=variable,geom=c("line","point")) +
geom_text(data= subset(d, variable == "cars"), aes(label=value))


with lattice, my best guess would be to use grid.text in a custom
panel function.


Hope this helps,

baptiste

On 8 Apr 2009, at 19:40, taz9 wrote:



Hi All,

I have a very simple graph:

cars <- c(1, 3, 6, 4, 9)
trucks <- c(2, 5, 4, 5, 12)
year <- c(2004, 2005, 2006, 2007, 2008)
df2<-data.frame(cars,trucks,year)
xyplot(cars+trucks~year, data=df2, type="o")

I need to show the values of "cars" on the graph. How can I do this?

Thanks.

--
View this message in context:
http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22956986.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.


_

Baptiste Auguié

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

Phone: +44 1392 264187

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

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




--
View this message in context: 
http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22957501.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.


_

Baptiste Auguié

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

Phone: +44 1392 264187

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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Minimum Spanning Tree

[Gábor Csárdi]

Actually,

library(igraph)

tab <- read.csv("http://www.nabble.com/file/p22957493/sp_matrix.csv";)
tab <- tab[,-1]

g <- graph.adjacency(as.matrix(tab), weighted=TRUE)
V(g)$label <- V(g)$name

mst <- as.undirected(minimum.spanning.tree(g))

lay <- layout.reingold.tilford(mst, root=which.max(degree(mst))-1)
lay <- cbind(lay[,2], lay[,1])# rotate
x11(width=15, height=8)
plot(mst, layout=lay, vertex.size=25, vertex.size2=10,
 vertex.shape="rectangle", asp=FALSE,
 vertex.label.cex=0.7, vertex.color="white")

works relatively well for me on your graph. It doesn't for you?

Best,
Gabor

On Wed, Apr 8, 2009 at 9:14 PM, jpearl01  wrote:
>
>
> Make the graph undirected first and then choose the right plotting
> parameters. E.g. the following works fine for me:
>
> set.seed(2)
> g <- erdos.renyi.game(100, 300, type="gnm", directed=TRUE)
> E(g)$weight <- runif(ecount(g))
> mst <- minimum.spanning.tree(g)
>
> mst <- simplify(as.undirected(mst))
> lay <- layout.reingold.tilford(mst, root=which.max(degree(mst))-1)
> plot(mst, layout=lay, vertex.size=5, asp=FALSE, vertex.color=NA,
>       vertex.frame.color=NA)
>
> G.
>
>
> Thanks for all your help Gabor,  However, I'm unable to get a display in R
> which is at all readable.  I'm not sure why that is (if the data is
> generated randomly, like in your example the tree builds/reads fine...
> however on my dataset whichever node is picked as the root just completely
> over powers the rest of the nodes which all clump together in a manner that
> is unreadable).  My dataset is:
> http://www.nabble.com/file/p22957493/sp_matrix.csv sp_matrix.csv
> If you would like to take a look.  However, as it stands I think I'll have
> to find another option.  Thanks again for all your efforts.
>
> ~josh
> --
> View this message in context: 
> http://www.nabble.com/Minimum-Spanning-Tree-tp22934813p22957493.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.
>



-- 
Gabor Csardi  UNIL DGM

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

[R] Customize Legend in Juxtaposed Barplot

[Paul Jones]

Hello, I'm making barplots out of binary data,  typically 0=no 1=yes, 
but in some instances I have 2=Don't Know.


Anyway, when I create a juxtaposed barplot, what appears in the legend 
is 0,1,and 2.


I want it to read No, Yes, Don't Know.

I'm not sure how to do that.

In addition to this, the other variable that this juxtaposed barplot is 
grouped by is grade level, which is a string that takes the values of 
Elementary, Middle, and High Schools. When I create the Barplot it 
orders the groups in alphabetical order. Is there anyway I can change 
the order of the juxtaposed groupings to Elementary, Middle, and High 
School?


Also, is it possible to print in smaller text No,Yes, and Don't Know 
under each bar. If this is too hard, I'll just ignore it.



Here is an example of the code that I use.

## This is an R Script to analyze the data from the LVAC Survey

## set working directory and load data
setwd("/home/paul/LVAC/General/")
#load("/home/paul/LVAC/LVAC.RData")
General <- read.csv("/home/paul/LVAC/General.csv")

## This Data is split into three data frames for Elementary,Middle and 
High Schools


png(file="PTA.png")
   x <- table(General$PTA.Support, General$Grade.Level)
   #x <- table(x, exclude=NULL)
   #names(x) <- r ("Yes","No")
   yrange <- range(x,na.rm=TRUE)*1.2
   if (yrange[1]>0) yrange[1] <- 0
   if (yrange[2]<0) yrange[2] <- 0
   bplot <- barplot(x,col=rainbow(if(is.matrix(x)) dim(x) else 
length(x)),beside=TRUE,legend.text=TRUE,ylim=yrange,

   main ="Barplot of PTA Support in Grade Levels")
   text(bplot,x,labels=x,pos=3,offset=.5)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] xyplot - show values of a series on graph

[Deepayan Sarkar]

On Wed, Apr 8, 2009 at 12:25 PM, taz9  wrote:
>
> Thank you very much for your help. I tried to use lattice but I'm not sure
> how to restrict it to display only the values of "cars".
>
> xyplot(cars+trucks~year, data=df2, type="o",
> panel=function(x,y,...){
>        panel.xyplot(x,y,...)
>        grid.text(unit(x,"native"),unit(y,"native"),label=y, just="top")}
>   )

A custom 'panel.groups' will let you condition on group number (see
?panel.superpose):

xyplot(cars+trucks~year, data=df2, type="o",
   panel = panel.superpose,
   panel.groups = function(x, y, ..., group.number) {
   panel.xyplot(x,y,...)
   if (group.number == 1) {
   require(grid)
   grid.text(unit(x,"native"),
 unit(y,"native"),
 label=y, just="top")
   }
   })

-Deepayan

>
>
>
> baptiste auguie-2 wrote:
>>
>> with ggplot2,
>>
>> d <- melt(df2,id="year")
>> qplot(year,value,data=d,colour=variable,geom=c("line","point")) +
>> geom_text(data= subset(d, variable == "cars"), aes(label=value))
>>
>>
>> with lattice, my best guess would be to use grid.text in a custom
>> panel function.
>>
>>
>> Hope this helps,
>>
>> baptiste
>>
>> On 8 Apr 2009, at 19:40, taz9 wrote:
>>
>>>
>>> Hi All,
>>>
>>> I have a very simple graph:
>>>
>>> cars <- c(1, 3, 6, 4, 9)
>>> trucks <- c(2, 5, 4, 5, 12)
>>> year <- c(2004, 2005, 2006, 2007, 2008)
>>> df2<-data.frame(cars,trucks,year)
>>> xyplot(cars+trucks~year, data=df2, type="o")
>>>
>>> I need to show the values of "cars" on the graph. How can I do this?
>>>
>>> Thanks.
>>>
>>> --
>>> View this message in context:
>>> http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22956986.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.
>>
>> _
>>
>> Baptiste Auguié
>>
>> School of Physics
>> University of Exeter
>> Stocker Road,
>> Exeter, Devon,
>> EX4 4QL, UK
>>
>> Phone: +44 1392 264187
>>
>> http://newton.ex.ac.uk/research/emag

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] MLE for bimodal distribution

[Ben Bolker]

_nico_ wrote:
> 
> Hello everyone,
> 
> I'm trying to use mle from package stats4 to fit a bi/multi-modal
> distribution to some data, but I have some problems with it.
> Here's what I'm doing (for a bimodal distribution):
> 
> # Build some fake binormally distributed data, the procedure fails also
> with real data, so the problem isn't here
> data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))
> # Just to check it's bimodal
> plot(density(data))
> f = function(m, s, m2, s2, w)
> {
> -log( sum(w*dnorm(data, mean=m, sd=s)) + sum((1-w)*dnorm(data, mean=m2,
> sd=s2)) )
> }
> 
> res = mle(f, start=list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6))
> 
> This gives an error:
> Error in optim(start, f, method = method, hessian = TRUE, ...) : 
>   non-finite finite-difference value [2]
> In addition: There were 50 or more warnings (use warnings() to see the
> first 50)
> And the warnings are about dnorm producing NaNs
> 
> So, my questions are:
> 1) What does "non-finite finite-difference value" mean? My guess would be
> that an Inf was passed somewhere where a finite number was expected...?
> I'm not sure how optim works though...
> 2) Is the log likelihood function I wrote correct? Can I use the same type
> of function for 3 modes?
> 3) What should I do to avoid function failure? I tried by changing the
> parameters, but it did not work.
> 4) Can I put constraints to the parameters? In particular, I would like w
> to be between 0 and 1.
> 5) Is there a way to get the log-likelihood value, so that I can compare
> different extimations?
> 6) Do you have any (possibly a bit "pratically oriented") readings about
> MLE to suggest?
> 
> Thanks in advance
> Nico
> 

  Here's some tweaked code that works.
Read about the "L-BFGS-B" method in ?optim to constrain parameters,
although you will have some difficulty making this work for more than
two components.  I think there's also a mixture model problem in
Venables & Ripley (MASS).

  There are almost certainly more efficient approaches for mixture
model fitting, although this "brute force" approach should be fine
if you don't need to do anything too complicated.
(RSiteSearch("mixture model"))

set.seed(1001)
data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))
# Just to check it's bimodal
plot(density(data))
f = function(m, s, m2, s2, w)
{
  -sum(log(w*dnorm(data, mean=m, sd=s)+
   (1-w)*dnorm(data, mean=m2, sd=s2)))
}

library(stats4)
start0 <- list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6)
do.call("f",start0) ## OK
res = mle(f, start=start0)

with(as.list(coef(res)),
 curve(w*dnorm(x,m,s)+(1-w)*dnorm(x,m2,s2),add=TRUE,col=2))


-- 
View this message in context: 
http://www.nabble.com/MLE-for-bimodal-distribution-tp22954970p22957613.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] rhipe v0.1

[Saptarshi Guha]

Greetings,
I would like to announce the 0.1 release of RHIPE:R and Hadoop
Integrated Processing Environment.  The website is located at :
http://www.stat.purdue.edu/~sguha/rhipe.
The download link is the bottom most
link on left side of the page.

RHIPE works on top of Hadoop, providing the R user a way to distribute
commands
over the Hadoop computing framework.

The 0.1 release has one main command - rhlapply which is a parallel version
of
lapply.  rhlapply outputs the results to Sequence files on the  Hadoop
Distributed Filesystem, and
RHIPE comes with commands to read R objects fromthese Sequence files(a
Hadoop file
format). rhlapply has features to share files/load libraries/execute code on
the cluster machines
and collect side effect files.

Since RHIPE uses Hadoop for distributing computation, it also benefits from
Hadoops stability: load balancing and machine failure recovery being two
important features, scheduling of jobs etc.

RHIPE also implements a basic Shared Associative Space via IBM's TSpaces -
with
commands like rhput,rhtake,rhread etc. This is optional.

See the website for details and performance results. Unfortunately, I
haven't
had the chance to compare with rmpi and snow.

TODO(will appear soon ): another function: rhmr - mapreduce using R.

Regards
Saptarshi Guha

[[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] Lattice Groups

[Lyman, Mark]

I guess I didn't look too closely. I didn't even notice the points were gone. 
Thanks, Deepayan.

Mark Lyman


-Original Message-
From: Deepayan Sarkar [mailto:deepayan.sar...@gmail.com] 
Sent: Wednesday, April 08, 2009 1:07 PM
To: Lyman, Mark
Cc: r-help@r-project.org; dickywe...@hotmail.com
Subject: Re: [R] Lattice Groups

On Wed, Apr 8, 2009 at 10:36 AM, Lyman, Mark  wrote:
> I don't understand your first question, but, since no one else has
> responded I can answer your second question. panel.bwplot, unlike
> panel.xyplot doesn't use panel.superpose when groups is not NULL. In
> order to get an analogous result you need to specify that you want to
> use panel.superpose.
>
> cols <- c("Sepal.Width", "Petal.Length", "Petal.Width")
> stackedData <- stack(iris[, cols])
> df <- data.frame(y = stackedData$values, x = rep(iris$Species, 3), which
> = gl(3, nrow(iris)))
>
> bwplot(y ~ x:which, data = df, groups = which, panel=panel.superpose,
> panel.groups = panel.bwplot)
>
> If you don't like the default colors, you can set the fill colors with
> par.settings like:
>
> bwplot(y ~ x:which, data = df, groups = which, panel=panel.superpose,
> panel.groups = panel.bwplot,
> par.settings=list(superpose.symbol=list(fill=2:4)))

And to answer the first question: using panel.superpose hijacks the
parameters of the median spot, but they can be supplied explicity:

bwplot(y ~ x:which, data = df, groups = which, panel=panel.superpose,
panel.groups = panel.bwplot,
par.settings=list(superpose.symbol=list(fill=2:4)), col = "black", pch = 16)

-Deepayan

>
> Without the groups, the fill colors are controlled like this
> bwplot(y~x:which, data = df,
> par.settings=list(box.rectangle=list(fill=2:4)))
>
> Although if you have groups, using the groups argument is probably
> better.
>
> Mark Lyman
>
>
> Message: 41
> Date: Tue, 7 Apr 2009 10:50:33 +0100
> From: Richard Weeks 
> Subject: [R] Lattice Groups
> To: 
> Message-ID: 
> Content-Type: text/plain
>
>
> Hi all,
>
>
>
> I'm trying to achieve a few things using the lattice package but am
> failing miserably.
>
> I am plotting side by side box plots and using a grouping variable, e.g.
>
>
>
> cols <- c("Sepal.Width", "Petal.Length", "Petal.Width")
> stackedData <- stack(iris[, cols])
> df <- data.frame(y = stackedData$values, x = rep(iris$Species, 3), which
> = gl(3, nrow(iris)))
>
> bwplot(y ~ x:which, data = df, group = which, panel.groups =
> panel.bwplot)
>
>
>
> My questions are
>
> 1) How am I able to retain the median spot in the boxes?
>
> 2) How can I change the fill using the par.settings argument rather than
> fill =1:3 say?
>
>
>
> Best wishes,
>
>
>
> Biff
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Display a very low p-value

[Ben Bolker]

Dimitris Rizopoulos-4 wrote:
> 
> in this case you need to use the 'lower.tail' argument, e.g.,
> 
> pnorm(8:15, lower.tail = FALSE)
> 
> 
> Bhoom Suktitipat wrote:
>> Hello R-Help,
>> 
>> I ran some analysis and were hit with some low Z-score. I tried to
>> convert
>> it to a p-value, however, it seems like the ceiling is around 1e-16.
>> 
>>> 1-pnorm(8)
>> [1] 6.661338e-16
>>> 1-pnorm(9)
>> [1] 0
>> 
>> Do you have any suggestion how I can display a very low p-value in the
>> form
>> of scientific number format?
> 
> 

For what it's worth, the pnorm calculation will still underflow for Z values
of
greater than about 37 -- 

 pnorm(37:39,lower.tail=FALSE)
[1] 5.725571e-300  0.00e+00  0.00e+00

  This is just a limitation of double precision floating-point arithmetic
...

 curve(pnorm(x,lower.tail=FALSE),from=30,to=40,log="y")
.Machine$double.xmin


-- 
View this message in context: 
http://www.nabble.com/Display-a-very-low-p-value-tp22956997p22957574.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] Minimum Spanning Tree

[jpearl01]

Make the graph undirected first and then choose the right plotting
parameters. E.g. the following works fine for me:

set.seed(2)
g <- erdos.renyi.game(100, 300, type="gnm", directed=TRUE)
E(g)$weight <- runif(ecount(g))
mst <- minimum.spanning.tree(g)

mst <- simplify(as.undirected(mst))
lay <- layout.reingold.tilford(mst, root=which.max(degree(mst))-1)
plot(mst, layout=lay, vertex.size=5, asp=FALSE, vertex.color=NA,
   vertex.frame.color=NA)

G.


Thanks for all your help Gabor,  However, I'm unable to get a display in R
which is at all readable.  I'm not sure why that is (if the data is
generated randomly, like in your example the tree builds/reads fine...
however on my dataset whichever node is picked as the root just completely
over powers the rest of the nodes which all clump together in a manner that
is unreadable).  My dataset is:
http://www.nabble.com/file/p22957493/sp_matrix.csv sp_matrix.csv 
If you would like to take a look.  However, as it stands I think I'll have
to find another option.  Thanks again for all your efforts.

~josh
-- 
View this message in context: 
http://www.nabble.com/Minimum-Spanning-Tree-tp22934813p22957493.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] xyplot - show values of a series on graph

[taz9]

Thank you very much for your help. I tried to use lattice but I'm not sure
how to restrict it to display only the values of "cars".

xyplot(cars+trucks~year, data=df2, type="o",
panel=function(x,y,...){
panel.xyplot(x,y,...)
grid.text(unit(x,"native"),unit(y,"native"),label=y, just="top")}
   )



baptiste auguie-2 wrote:
> 
> with ggplot2,
> 
> d <- melt(df2,id="year")
> qplot(year,value,data=d,colour=variable,geom=c("line","point")) +
> geom_text(data= subset(d, variable == "cars"), aes(label=value))
> 
> 
> with lattice, my best guess would be to use grid.text in a custom  
> panel function.
> 
> 
> Hope this helps,
> 
> baptiste
> 
> On 8 Apr 2009, at 19:40, taz9 wrote:
> 
>>
>> Hi All,
>>
>> I have a very simple graph:
>>
>> cars <- c(1, 3, 6, 4, 9)
>> trucks <- c(2, 5, 4, 5, 12)
>> year <- c(2004, 2005, 2006, 2007, 2008)
>> df2<-data.frame(cars,trucks,year)
>> xyplot(cars+trucks~year, data=df2, type="o")
>>
>> I need to show the values of "cars" on the graph. How can I do this?
>>
>> Thanks.
>>
>> -- 
>> View this message in context:
>> http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22956986.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.
> 
> _
> 
> Baptiste Auguié
> 
> School of Physics
> University of Exeter
> Stocker Road,
> Exeter, Devon,
> EX4 4QL, UK
> 
> Phone: +44 1392 264187
> 
> http://newton.ex.ac.uk/research/emag
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22957501.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] predict "interval" for lmRob?

[Galkowski, Jan]

[snip]


>Discarding actual data points always makes me nervous.  Sometimes the points 
>we want to discard are actually the most interesting.

No doubt this is true, and there's a lot of information in those outliers, a 
lot of structure.  For instance, in this case, one part of the outlier 
population is actually identifiable as a valid part of the primary dataset 
having the abscissa shifted by a known constant.  The mechanism for that is 
known, so it could be defended that this portion of the outliers could be added 
back into the main population by removing the shift.  Still, not much else is 
known about its surround, so I/we wonder what else we'd be picking up if that 
were done. 

But for the primary application, which is a calibration, I think going after 
the main population is what's wanted right now.

Thanks again.

 - Jan

[snip]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 "interval" for lmRob?

[Greg Snow]

> -Original Message-
> From: Galkowski, Jan [mailto:jgalk...@akamai.com]
> Sent: Wednesday, April 08, 2009 12:28 PM
> To: Greg Snow; r-help@r-project.org
> Subject: RE: predict "interval" for lmRob?

[snip]

> It sounds to me like I might use the robust regression to decide what
> to discard and then apply standard linear "lm" to the remainder,
> minding the diagnostics. Should they prove favorable, I'll proceed with
> the result of "lm".

Discarding actual data points always makes me nervous.  Sometimes the points we 
want to discard are actually the most interesting.

> Thanks for pointing out the limitations of "robust" and its kin for me.

I consider anything that encourages me to ask questions, contemplate the 
answers, and really think about my data and scientific question to be a benefit 
rather than a limitation (one of the reasons I like R so much).

> BTW, if "robust" does not adopt a normal model for the y variable,
> what's the proper interpretation of the standard errors for slope and
> intercept it yields?  A reference?

Well there are several references on the help page for lmRob, there is also a 
section in MASS (the book).  But I think that while some of the techniques may 
have been developed for one particular distribution, it has been found that 
they work for a larger set of distributions and the theory does not depend on a 
particular distribution (you have to decide which makes the most sense for your 
data/application area).  For simulations to show that they work I have seen:  
mixture of 2 normals, same mean but one with a much larger variance (giving the 
outliers), mixture of a normal and a t/cauchy, mixture of a normal and a gamma 
(some skewness/outliers), mixture of 2 normals with different means (outliers 
come from a different population mingled in with the population of interest and 
not easily distinguished), etc.



-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] April R/Splus Courses: Back2back (1) R/S+ Fundamentals and (2) R/S-Plus Advanced Programming. in San Francisco and New York City

[Sue Turner]

By popular demand, XLSolutions has scheduled first back2back courses 
in New York City  and San Francisco taught by top R/S+ gurus!

This is a first in our 9 years of teaching R/Splus Courses.

West Coast<---back2back--->  East Coast


(1) R/S-PLUS Fundamentals and Programming Techniques
http://www.xlsolutions-corp.com/coursedetail.asp?id=30

* San Francisco * April 23-24, 2009
* New York City * April 27-28, 2009

(2) R/S+ System: Advanced Programming 
http://www.xlsolutions-corp.com/coursedetail.asp?id=16

* San Francisco * April 27-28, 2009
* New York City * April 30-May 1st, 2009


Ask for group discount and reserve your seat Now - Earlybird Rates.
Payment due after the class! Email Sue Turner: s...@xlsolutions-corp.com

http://www.xlsolutions-corp.com/rplus.asp

Please let us know if you and your colleagues are interested in this
class to take advantage of group discount. Register now to secure your
seat!

Cheers,
Elvis Miller, PhD
Manager Training.
XLSolutions Corporation
206 686 1578
www.xlsolutions-corp.com
el...@xlsolutions-corp.com

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

Re: [R] Lattice Groups

[Deepayan Sarkar]

On Wed, Apr 8, 2009 at 10:36 AM, Lyman, Mark  wrote:
> I don't understand your first question, but, since no one else has
> responded I can answer your second question. panel.bwplot, unlike
> panel.xyplot doesn't use panel.superpose when groups is not NULL. In
> order to get an analogous result you need to specify that you want to
> use panel.superpose.
>
> cols <- c("Sepal.Width", "Petal.Length", "Petal.Width")
> stackedData <- stack(iris[, cols])
> df <- data.frame(y = stackedData$values, x = rep(iris$Species, 3), which
> = gl(3, nrow(iris)))
>
> bwplot(y ~ x:which, data = df, groups = which, panel=panel.superpose,
> panel.groups = panel.bwplot)
>
> If you don't like the default colors, you can set the fill colors with
> par.settings like:
>
> bwplot(y ~ x:which, data = df, groups = which, panel=panel.superpose,
> panel.groups = panel.bwplot,
> par.settings=list(superpose.symbol=list(fill=2:4)))

And to answer the first question: using panel.superpose hijacks the
parameters of the median spot, but they can be supplied explicity:

bwplot(y ~ x:which, data = df, groups = which, panel=panel.superpose,
panel.groups = panel.bwplot,
par.settings=list(superpose.symbol=list(fill=2:4)), col = "black", pch = 16)

-Deepayan

>
> Without the groups, the fill colors are controlled like this
> bwplot(y~x:which, data = df,
> par.settings=list(box.rectangle=list(fill=2:4)))
>
> Although if you have groups, using the groups argument is probably
> better.
>
> Mark Lyman
>
>
> Message: 41
> Date: Tue, 7 Apr 2009 10:50:33 +0100
> From: Richard Weeks 
> Subject: [R] Lattice Groups
> To: 
> Message-ID: 
> Content-Type: text/plain
>
>
> Hi all,
>
>
>
> I'm trying to achieve a few things using the lattice package but am
> failing miserably.
>
> I am plotting side by side box plots and using a grouping variable, e.g.
>
>
>
> cols <- c("Sepal.Width", "Petal.Length", "Petal.Width")
> stackedData <- stack(iris[, cols])
> df <- data.frame(y = stackedData$values, x = rep(iris$Species, 3), which
> = gl(3, nrow(iris)))
>
> bwplot(y ~ x:which, data = df, group = which, panel.groups =
> panel.bwplot)
>
>
>
> My questions are
>
> 1) How am I able to retain the median spot in the boxes?
>
> 2) How can I change the fill using the par.settings argument rather than
> fill =1:3 say?
>
>
>
> Best wishes,
>
>
>
> Biff
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] help with random forest package

[Liaw, Andy]

I'm not quite sure what you're asking.  RF predicts by classifying the
new observation using all trees in the forest, and take plural vote.
The predict() method for randomForest objects does that for you.  The
getTree() function shows you what each individual tree is like (not
visually, just the underlying representation of the tree).
 
Andy




From: Chrysanthi A. [mailto:chrys...@gmail.com] 
Sent: Wednesday, April 08, 2009 2:56 PM
To: Liaw, Andy
Cc: r-help@r-project.org
Subject: Re: [R] help with random forest package


Many thanks for the reply.

So, extracting the votes, how can we clarify the classification
result? If I want to predict in which class will be included an unknown
sample, what is the rule that will give me that?

Thanks a lot,

Chrysanthi.




2009/4/8 Liaw, Andy 


The source code of the whole package is available on
CRAN.  All packages
are submitted to CRAN is source form.

There's no "rule" per se that gives the final
prediction, as the final
prediction is the result of plural vote by all trees in
the forest.

You may want to look at the varUsed() and getTree()
functions.

Andy

From:  Chrysanthi A.

> Hello,
>
> I am a phd student in Bioinformatics and I am using
the Random Forest
> package in order to classify my data, but I have some
questions.
> Is there a function in order to visualize the trees,
so as to
> get the rules?
> Also, could you please provide me with the code of
> "randomForest" function,
> as I would like to see how it works. I was wondering
if I can get the
> classification having the most votes over all the
trees in
> the forest (the
> final rules that will give me the final
classification).
> Also, is there a
> possibility to get a vector with the attributes that
are
> being selected for
> each node during the construction of each tree? I
mean, that
> I would like to
> know the m< the M input
> attributes.. Are they selected randomly? Is there a
> possibility to select
> the same variable in subsequent nodes?
>
> Thanks a lot,
>
> Chrysanthi.
>

>   [[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.
>
Notice:  This e-mail message, together with any
attachments, contains
information of Merck & Co., Inc. (One Merck Drive,
Whitehouse Station,
New Jersey, USA 08889), and/or its affiliates (which may
be known
outside the United States as Merck Frosst, Merck Sharp &
Dohme or
MSD and in Japan, as Banyu - direct contact information
for affiliates is
available at http://www.merck.com/contact/contacts.html)
that may be
confidential, proprietary copyrighted and/or legally
privileged. It is
intended solely for the use of the individual or entity
named on this
message. If you are not the intended recipient, and have
received this
message in error, please notify us immediately by reply
e-mail and
then delete it from your system.




Notice:  This e-mail message, together with any attachme...{{dropped:15}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Tinn-R pdf()

[Bert Gunter]

I **think** the explanation is:

When sending code to R through TINN-R, TINN-R **sources** the code into R.
Sourced code is **not** autoprinted. Indeed, ?source explicitly tells you
this.

-- Bert 


Bert Gunter
Genentech Nonclinical Biostatistics
650-467-7374

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Deepayan Sarkar
Sent: Wednesday, April 08, 2009 9:20 AM
To: Tobias Verbeke
Cc: r-help@r-project.org; ONKELINX, Thierry
Subject: Re: [R] Tinn-R pdf()

On 4/8/09, Tobias Verbeke  wrote:
> Hi Henning,
>
>
> > thanks for your help, with solved the problem, although i don't why,
> because when using the R editor accessible via the R console i created
many
> many lattice plots with the code i posted, i.e. without the print()
command.
> >
>
>  At the command line, R objects (including lattice plots which are
>  objects) get auto-printed, i.e. the print method is invoked automatically
> on these objects.
>
>  This is not the case when you write to a pdf file.

That claim is not accurate. Auto-printing should not be affected by
whether a pdf device is open at the moment.

I'm not familiar with TINN-R, but it's possible it executes R commands
through some sort of wrapper function, and that is failing to emulate
auto-printing correctly.

-Deepayan

>
>  Best,
>  Tobias
>
>
>
> >
> > >  Original-Nachricht 
> > > Datum: Wed, 8 Apr 2009 12:29:58 +0200
> > > Von: "ONKELINX, Thierry" 
> > > An: "Henning Wildhagen" , r-help@r-project.org
> > > Betreff: RE: [R] Tinn-R pdf()
> > >
> > >Dear Henning,
> > >
> > > You need to print() lattice plots when using a device:
> > >  library(lattice)
> > >
> > > pdf("plot1.pdf")
> > > PLOT<-(xyplot, ...)
> > > print(PLOT)
> > > dev.off()
> > >
> > > So this is not due to TINN-R.
> > >
> > > HTH,
> > >
> > > Thierry
> > >
> > >
> 
> > > 
> > > ir. Thierry Onkelinx
> > > Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> > > and Forest
> > > Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
> > > methodology and quality assurance
> > > Gaverstraat 4
> > > 9500 Geraardsbergen
> > > Belgium tel. + 32 54/436 185
> > > thierry.onkel...@inbo.be www.inbo.be
> > > To call in the statistician after the experiment is done may be no
more
> > > than asking him to perform a post-mortem examination: he may be able
to
> > > say what the experiment died of.
> > > ~ Sir Ronald Aylmer Fisher
> > >
> > > The plural of anecdote is not data.
> > > ~ Roger Brinner
> > >
> > > The combination of some data and an aching desire for an answer does
not
> > > ensure that a reasonable answer can be extracted from a given body of
> > > data.
> > > ~ John Tukey
> > >
> > > -Oorspronkelijk bericht-
> > > Van: r-help-boun...@r-project.org
> [mailto:r-help-boun...@r-project.org]
> > > Namens Henning Wildhagen
> > > Verzonden: woensdag 8 april 2009 11:25
> > > Aan: r-help@r-project.org
> > > Onderwerp: [R] Tinn-R pdf()
> > >
> > > Dear R and Tinn-R users,
> > >
> > > i recently switched to Tinn-R and sending code to R works fine (R
2.8.1,
> > >
> > > Tinn-R 2.2.0.2, OS Windows XP). However, i encountered a problem when
> trying to send plots to pdf files like this:
> > >
> > > library(lattice)
> > >
> > > pdf("plot1.pdf")
> > > PLOT<-(xyplot, ...)
> > > PLOT
> > > dev.off()
> > >
> > > The file "plot1.pdf" is created, but it is empty.
> > > If i paste the code above directly into the R console and run it, the
> > > file "plot1.pdf" is created and in this case contains "PLOT".
> > >
> > > I guess that some settings in Tinn-R are wrong, but i have no idea
> > > which. Maybe someone has a suggestion?
> > >
> > > Thanks,
> > >
> > > Henning
> > > --

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] help with random forest package

[Chrysanthi A.]

Many thanks for the reply.

So, extracting the votes, how can we clarify the classification result? If I
want to predict in which class will be included an unknown sample, what is
the rule that will give me that?

Thanks a lot,

Chrysanthi.



2009/4/8 Liaw, Andy 

> The source code of the whole package is available on CRAN.  All packages
> are submitted to CRAN is source form.
>
> There's no "rule" per se that gives the final prediction, as the final
> prediction is the result of plural vote by all trees in the forest.
>
> You may want to look at the varUsed() and getTree() functions.
>
> Andy
>
> From:  Chrysanthi A.
> > Hello,
> >
> > I am a phd student in Bioinformatics and I am using the Random Forest
> > package in order to classify my data, but I have some questions.
> > Is there a function in order to visualize the trees, so as to
> > get the rules?
> > Also, could you please provide me with the code of
> > "randomForest" function,
> > as I would like to see how it works. I was wondering if I can get the
> > classification having the most votes over all the trees in
> > the forest (the
> > final rules that will give me the final classification).
> > Also, is there a
> > possibility to get a vector with the attributes that are
> > being selected for
> > each node during the construction of each tree? I mean, that
> > I would like to
> > know the m< > the M input
> > attributes.. Are they selected randomly? Is there a
> > possibility to select
> > the same variable in subsequent nodes?
> >
> > Thanks a lot,
> >
> > Chrysanthi.
> >
> >   [[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.
> >
> Notice:  This e-mail message, together with any attach...{{dropped:17}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] xyplot - show values of a series on graph

[baptiste auguie]

with ggplot2,

d <- melt(df2,id="year")
qplot(year,value,data=d,colour=variable,geom=c("line","point")) +
geom_text(data= subset(d, variable == "cars"), aes(label=value))


with lattice, my best guess would be to use grid.text in a custom  
panel function.



Hope this helps,

baptiste

On 8 Apr 2009, at 19:40, taz9 wrote:



Hi All,

I have a very simple graph:

cars <- c(1, 3, 6, 4, 9)
trucks <- c(2, 5, 4, 5, 12)
year <- c(2004, 2005, 2006, 2007, 2008)
df2<-data.frame(cars,trucks,year)
xyplot(cars+trucks~year, data=df2, type="o")

I need to show the values of "cars" on the graph. How can I do this?

Thanks.

--
View this message in context: 
http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22956986.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.


_

Baptiste Auguié

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

Phone: +44 1392 264187

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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Display a very low p-value

[Dimitris Rizopoulos]

in this case you need to use the 'lower.tail' argument, e.g.,

pnorm(8:15, lower.tail = FALSE)


I hope it helps.

Best,
Dimitris


Bhoom Suktitipat wrote:

Hello R-Help,

I ran some analysis and were hit with some low Z-score. I tried to convert
it to a p-value, however, it seems like the ceiling is around 1e-16.


1-pnorm(8)

[1] 6.661338e-16

1-pnorm(9)

[1] 0

Do you have any suggestion how I can display a very low p-value in the form
of scientific number format?

Thank you very much.

Regards,
Bhoom

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



--
Dimitris Rizopoulos
Assistant Professor
Department of Biostatistics
Erasmus University Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Display a very low p-value

[stephen sefick]

?options look at scipen

On Wed, Apr 8, 2009 at 2:38 PM, Bhoom Suktitipat  wrote:
> Hello R-Help,
>
> I ran some analysis and were hit with some low Z-score. I tried to convert
> it to a p-value, however, it seems like the ceiling is around 1e-16.
>
>> 1-pnorm(8)
> [1] 6.661338e-16
>> 1-pnorm(9)
> [1] 0
>
> Do you have any suggestion how I can display a very low p-value in the form
> of scientific number format?
>
> Thank you very much.
>
> Regards,
> Bhoom
>
>        [[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.
>



-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] help with random forest package

[Liaw, Andy]

The source code of the whole package is available on CRAN.  All packages
are submitted to CRAN is source form.

There's no "rule" per se that gives the final prediction, as the final
prediction is the result of plural vote by all trees in the forest.

You may want to look at the varUsed() and getTree() functions.

Andy 

From:  Chrysanthi A.
> Hello,
> 
> I am a phd student in Bioinformatics and I am using the Random Forest
> package in order to classify my data, but I have some questions.
> Is there a function in order to visualize the trees, so as to 
> get the rules?
> Also, could you please provide me with the code of 
> "randomForest" function,
> as I would like to see how it works. I was wondering if I can get the
> classification having the most votes over all the trees in 
> the forest (the
> final rules that will give me the final classification). 
> Also, is there a
> possibility to get a vector with the attributes that are 
> being selected for
> each node during the construction of each tree? I mean, that 
> I would like to
> know the m< the M input
> attributes.. Are they selected randomly? Is there a 
> possibility to select
> the same variable in subsequent nodes?
> 
> Thanks a lot,
> 
> Chrysanthi.
> 
>   [[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.
> 
Notice:  This e-mail message, together with any attachme...{{dropped:12}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] xyplot - show values of a series on graph

[taz9]

Hi All,

I have a very simple graph:

cars <- c(1, 3, 6, 4, 9)
trucks <- c(2, 5, 4, 5, 12)
year <- c(2004, 2005, 2006, 2007, 2008)
df2<-data.frame(cars,trucks,year)
xyplot(cars+trucks~year, data=df2, type="o")

I need to show the values of "cars" on the graph. How can I do this?

Thanks.

-- 
View this message in context: 
http://www.nabble.com/xyplot---show-values-of-a-series-on-graph-tp22956986p22956986.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] Display a very low p-value

[Bhoom Suktitipat]

Hello R-Help,

I ran some analysis and were hit with some low Z-score. I tried to convert
it to a p-value, however, it seems like the ceiling is around 1e-16.

> 1-pnorm(8)
[1] 6.661338e-16
> 1-pnorm(9)
[1] 0

Do you have any suggestion how I can display a very low p-value in the form
of scientific number format?

Thank you very much.

Regards,
Bhoom

[[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] predict "interval" for lmRob?

[Galkowski, Jan]

Hi Greg,

Thanks for your guidance. 

In this case, the evidence is that the primary subpopulation of the data, 
accounting for  observes the standard statistical model in the sense that Rice 
uses the term.  It may by all accounts be normally distributed, and a Q-Q shows 
a large portion of the primary subpopulation behaves that way, out to 2 
theoretical quantiles. But, for the measurement ranges of interest, the 
complement of the "normal subpopulation", accounting for some 20% of the total 
two million data points, behaves in other ways, which are, as a matter of fact, 
poorly understood.  That's not likely to change soon.

The choice of a robust regression framework and of "robust" (and possibly 
"quantreg" as Prof Koenker suggested) was simply to automatically fit a line to 
the primary subpopulation, without having to make arbitrary choices as what to 
keep or what to discard. Also, use of any preexisting package was simply 
pursued as a timesaver, worksaver, and to have some conceptual framework within 
to proceed other than just throwing least squares at arbitrarily chosen 
subsets.  

It sounds to me like I might use the robust regression to decide what to 
discard and then apply standard linear "lm" to the remainder, minding the 
diagnostics. Should they prove favorable, I'll proceed with the result of "lm".

Thanks for pointing out the limitations of "robust" and its kin for me. 

BTW, if "robust" does not adopt a normal model for the y variable, what's the 
proper interpretation of the standard errors for slope and intercept it yields? 
 A reference?

 - Jan

-Original Message-
From: Greg Snow [mailto:greg.s...@imail.org] 
Sent: Wednesday, April 08, 2009 1:20 PM
To: Galkowski, Jan; r-help@r-project.org
Subject: RE: predict "interval" for lmRob?

Your problem is related to the theory underlying linear models (and is an 
example as to why it is important to understand the theory, not just know how 
to plug numbers into a computer).

The lm function is based on theory that assumes the y variable in normally 
distributed with the mean of that normal based on the model and the x values.  
This allows the predict function for lm to create prediction intervals based on 
the normal distribution, the predicted mean of that distribution, the estimated 
standard deviation, and the uncertainty in the predicted mean.  Note that if 
your y variable is not normally distributed, but the sample size is large 
enough for the Central Limit Theorem to hold, then the confidence intervals 
will be approximately correct, but the prediction intervals will probably not 
be.

When you switch to a robust regression approach, the assumption is that the y 
variable is not normal, so a prediction interval based on the normal 
distribution does not make sense.  To get an appropriate prediction interval 
you need some information on what the distribution of the y values is 
(conditional on the model), but most robust techniques are not based on a 
specific distribution, just some properties of the distribution.  Without some 
information (or at least an assumption) on that distribution, the predict 
method cannot create prediction intervals.

I know that this does not answer your question, but hopefully helps you to 
understand what is happening.  Think about what your actual scientific question 
is, it may be that you can answer the question without prediction intervals.

If you feel that you really need the prediction intervals, then you will need 
to do some additional background research into what distribution you think the 
data comes from, then you can proceed from there.  Some options include fitting 
a model based on that distribution, simulating data from the distribution given 
the model estimates and the uncertainty in those estimates, quantile 
regression, mixture of regressions, and others.

Hope this helps,

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Galkowski, Jan
> Sent: Wednesday, April 08, 2009 9:32 AM
> To: r-help@r-project.org
> Subject: [R] predict "interval" for lmRob?
> 
> lm's "predict" function offers an "interval" parameter to choose
> between 'confidence' and 'prediction' bands. In the package "robust"
> and for "lmRob", there is also a "predict" but it lacks such a
> parameter, and the documented "type" parameter has only "response"
> offerred.  Is there some way of obtaining prediction bands from lmRob?
> Is there an alternative robust (linear) regression package that offers
> such a capability?
> 
> Thanks for any and all help.
> 
>   - Jan Galkowski, Akamai Technologies, Cambridge, MA.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.htm

Re: [R] read.spss, locale and encodings

[Hans Ekbrand]

On Wed, Apr 08, 2009 at 07:12:23PM +0200, Peter Dalgaard wrote:
> Apparently, you can work around it like this
>
> lc <- Sys.setlocale("LC_CTYPE")
> Sys.setlocale("LC_CTYPE", "da_DK")
> x <- read.spss("~/Desktop/downloads/test.sav", reencode = "latin1")
> Sys.setlocale("LC_CTYPE", lc)
>
> -- which doesn't strike me as particularly logical, but whatever works

THANKS a lot Peter! This works perfectly! I had been struggling with
this problem way too long...

-- 
Hans Ekbrand (http://sociologi.cjb.net) 
GnuPG key: 1024D/7050614E
Fingerprint: 1408 C8D5 1E7D 4C9C C27E 014F 7C2C 872A 7050 614E
Learn about secure email at http://www.gnupg.org


signature.asc
Description: Digital signature
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Lattice Groups

[Lyman, Mark]

I don't understand your first question, but, since no one else has
responded I can answer your second question. panel.bwplot, unlike
panel.xyplot doesn't use panel.superpose when groups is not NULL. In
order to get an analogous result you need to specify that you want to
use panel.superpose.

cols <- c("Sepal.Width", "Petal.Length", "Petal.Width")
stackedData <- stack(iris[, cols])
df <- data.frame(y = stackedData$values, x = rep(iris$Species, 3), which
= gl(3, nrow(iris)))

bwplot(y ~ x:which, data = df, groups = which, panel=panel.superpose,
panel.groups = panel.bwplot)

If you don't like the default colors, you can set the fill colors with
par.settings like:

bwplot(y ~ x:which, data = df, groups = which, panel=panel.superpose,
panel.groups = panel.bwplot,
par.settings=list(superpose.symbol=list(fill=2:4)))

Without the groups, the fill colors are controlled like this
bwplot(y~x:which, data = df,
par.settings=list(box.rectangle=list(fill=2:4)))

Although if you have groups, using the groups argument is probably
better.

Mark Lyman


Message: 41
Date: Tue, 7 Apr 2009 10:50:33 +0100
From: Richard Weeks 
Subject: [R] Lattice Groups
To: 
Message-ID: 
Content-Type: text/plain


Hi all,

 

I'm trying to achieve a few things using the lattice package but am
failing miserably.

I am plotting side by side box plots and using a grouping variable, e.g.

 

cols <- c("Sepal.Width", "Petal.Length", "Petal.Width")
stackedData <- stack(iris[, cols])
df <- data.frame(y = stackedData$values, x = rep(iris$Species, 3), which
= gl(3, nrow(iris)))

bwplot(y ~ x:which, data = df, group = which, panel.groups =
panel.bwplot)

 

My questions are 

1) How am I able to retain the median spot in the boxes?

2) How can I change the fill using the par.settings argument rather than
fill =1:3 say?

 

Best wishes,

 

Biff

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Minimum Spanning Tree

[Gábor Csárdi]

On Wed, Apr 8, 2009 at 6:12 PM, jpearl01  wrote:
>
>
>>I am not sure what you mean. Of course you can plot it using different
>>layouts, e.g. with layout.reingold.tilford (after choosing the root
>>vertex in some way) and then it looks like a usual tree plot, but why
>>would that be any better?
>
> I'd like to be able to distinguish between the nodes better.  For example
> the image that I get looks like:
> http://www.nabble.com/file/p22954099/mst.gif mst.gif
>
> The tree nodes are just too close together to be able to read.  I think the
> problem might be that some distances between the nodes are *much* smaller
> than other edges.    The graph does not need to be directed, since the
> distance metric does not imply a direction. The most important part is just
> being able to see and differentiate between the nodes, and right now they
> seem to be all lumped together.  How can I make it more readable?

Make the graph undirected first and then choose the right plotting
parameters. E.g. the following works fine for me:

set.seed(2)
g <- erdos.renyi.game(100, 300, type="gnm", directed=TRUE)
E(g)$weight <- runif(ecount(g))
mst <- minimum.spanning.tree(g)

mst <- simplify(as.undirected(mst))
lay <- layout.reingold.tilford(mst, root=which.max(degree(mst))-1)
plot(mst, layout=lay, vertex.size=5, asp=FALSE, vertex.color=NA,
   vertex.frame.color=NA)

G.

> Thanks,
> ~josh
>
>
> --
> View this message in context: 
> http://www.nabble.com/Minimum-Spanning-Tree-tp22934813p22954099.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.
>



-- 
Gabor Csardi  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] predict "interval" for lmRob?

[Greg Snow]

Your problem is related to the theory underlying linear models (and is an 
example as to why it is important to understand the theory, not just know how 
to plug numbers into a computer).

The lm function is based on theory that assumes the y variable in normally 
distributed with the mean of that normal based on the model and the x values.  
This allows the predict function for lm to create prediction intervals based on 
the normal distribution, the predicted mean of that distribution, the estimated 
standard deviation, and the uncertainty in the predicted mean.  Note that if 
your y variable is not normally distributed, but the sample size is large 
enough for the Central Limit Theorem to hold, then the confidence intervals 
will be approximately correct, but the prediction intervals will probably not 
be.

When you switch to a robust regression approach, the assumption is that the y 
variable is not normal, so a prediction interval based on the normal 
distribution does not make sense.  To get an appropriate prediction interval 
you need some information on what the distribution of the y values is 
(conditional on the model), but most robust techniques are not based on a 
specific distribution, just some properties of the distribution.  Without some 
information (or at least an assumption) on that distribution, the predict 
method cannot create prediction intervals.

I know that this does not answer your question, but hopefully helps you to 
understand what is happening.  Think about what your actual scientific question 
is, it may be that you can answer the question without prediction intervals.

If you feel that you really need the prediction intervals, then you will need 
to do some additional background research into what distribution you think the 
data comes from, then you can proceed from there.  Some options include fitting 
a model based on that distribution, simulating data from the distribution given 
the model estimates and the uncertainty in those estimates, quantile 
regression, mixture of regressions, and others.

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Galkowski, Jan
> Sent: Wednesday, April 08, 2009 9:32 AM
> To: r-help@r-project.org
> Subject: [R] predict "interval" for lmRob?
> 
> lm's "predict" function offers an "interval" parameter to choose
> between 'confidence' and 'prediction' bands. In the package "robust"
> and for "lmRob", there is also a "predict" but it lacks such a
> parameter, and the documented "type" parameter has only "response"
> offerred.  Is there some way of obtaining prediction bands from lmRob?
> Is there an alternative robust (linear) regression package that offers
> such a capability?
> 
> Thanks for any and all help.
> 
>   - Jan Galkowski, Akamai Technologies, Cambridge, MA.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] read.spss, locale and encodings

[Peter Dalgaard]

Hans Ekbrand wrote:

On Wed, Apr 08, 2009 at 04:17:51PM +0200, Peter Dalgaard wrote:

Hans Ekbrand wrote:

Someone running foreign 8.34 that is willing to test my SPSS-file?

Someone with an SPSS file problem willing to help test the prereleases? :-)


http://sociologi.cjb.net/temp/test.sav


No joy.

> read.spss("~/Desktop/downloads/test.sav", reencode = "latin1")
Error in read.spss("~/Desktop/downloads/test.sav", reencode = "latin1") :
  error reading system-file header
In addition: Warning message:
In read.spss("~/Desktop/downloads/test.sav", reencode = "latin1") :
  ~/Desktop/downloads/test.sav: position 143: Variable name begins with 
invalid character


(I suppose the actual culprit could be number 144 which does indeed 
start with an A-ring ("ÅLDKAT"))


Apparently, you can work around it like this

lc <- Sys.setlocale("LC_CTYPE")
Sys.setlocale("LC_CTYPE", "da_DK")
x <- read.spss("~/Desktop/downloads/test.sav", reencode = "latin1")
Sys.setlocale("LC_CTYPE", lc)

-- which doesn't strike me as particularly logical, but whatever works

--
   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
~~ - (p.dalga...@biostat.ku.dk)  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] Reshape - strange outputs

[stephen sefick]

can you provide a reproducible example?

On Wed, Apr 8, 2009 at 11:59 AM, Steve Murray  wrote:
>
> Dear R Users,
>
> I am using the reshape package to reformat gridded data into column format 
> using the code shown below. However, when I display the resulting object, a 
> single column is fomed (instead of three) and all the latitude values (which 
> should be in either column one or two) are collected at the bottom. Also, the 
> NA values aren't removed, despite this being requested in the code.
>
> Code:
>
>
> # NetCDF file has been read in and is being processed...
>
> arunoff_1986_temp <- get.var.ncdf(netcdf_1036_temp, "arunoff")
>
>
> # Assign row and column names
>
> columnnames <- sprintf("%.2f", seq(from = -89.75, to = 89.75, length = 360))
> rnames <- sprintf("%.2f", seq(from = -179.75, to = 179.75, length = 720))
>
> colnames(arunoff_1986_temp) <- columnnames
> rownames(arunoff_1986_temp) <- rnames
>
>
> # Melt into columnar format
>
> arunoff_1986$Longitude <- rownames(arunoff_1986)
>    # Note: If I do: arunoff_1986$Latitude <- rownames(arunoff_1986)  (i.e. 
> change it to 'Latitude', I get the following: Warning message: In 
> arunoff_1986$Latitude <- rownames(arunoff_1986) : Coercing LHS to a list ), 
> thus proceed using 'Longitude' where no warning is apparent.
>
>
> arunoff_long_1986 <- melt(arunoff_1986, id.var="Longitude", na.rm=TRUE)
>
>
>> dim(arunoff_long_1986)
> [1] 259560      2
>
>>head(arunoff_long_1986, n=10)  # This displays what looks like one single 
>>column, but is in fact two: column entitled "L1" is empty until the end of 
>>the file, as shown by the 'tail' command below.
>   value L1
> 1
> 2
> 3
> 4
> 5
> 6
> 7
> 8
> 9
> 10
>
>> tail(arunoff_long_1986, n=10)
>             value       L1
> Latitude.351 85.25 Latitude
> Latitude.352 85.75 Latitude
> Latitude.353 86.25 Latitude
> Latitude.354 86.75 Latitude
> Latitude.355 87.25 Latitude
> Latitude.356 87.75 Latitude
> Latitude.357 88.25 Latitude
> Latitude.358 88.75 Latitude
> Latitude.359 89.25 Latitude
> Latitude.360 89.75 Latitude
>
>
>
> I'd be very grateful indeed if anyone is able to offer assistance by way of 
> pointing out what I've done wrong. I've spent a long time working on this and 
> trying various options, but am srill none the wiser! I'm aiming for three 
> columns: Latitude, Longitude, 'value'.
>
> Many thanks for any help offered,
>
> Steve
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Convert data frame containing time stamps to time series

[stephen sefick]

have you tried using zoo and then using the function as.ts()

On Wed, Apr 8, 2009 at 11:56 AM,   wrote:
> Converting dates is getting stranger still. I am coercing a data frame
> into a ts as follows:
>
>
> tst1<-as.POSIXct("1/21/09 5:01",format="%m/%d/%y %H:%M")
> tst2<-as.POSIXct("1/28/09 3:40",format="%m/%d/%y %H:%M")
> tsdat<-as.ts(dat,start=tst1,end=tst2,frequency=1)
>
> This generates a ts object. But strangely enough the first column of that
> matrix starts at the numeric value of 841 counts up to 1139 and then
> starts at 1 again, only to count up from there. The restart at 1 occurs at
> the first day "1/21/09" at 10:00:00.
>
> What is so special about that time? This phenomenon happens several times
> in the long file. But the restart count is always a different number.
> This creates a ramp with some bumps.
>
> Can anybody explain this?
> Thanks in advance,
> Alex van der Spek
>
>
>> I read records using scan:
>>
>> dat<-data.frame(scan(file="KDA.csv",what=list(t="%m/%d/%y
>> %H:%M",f=0,p=0,d=0,o=0,s=0,a=0,l=0,c=0),skip=2,sep=",",nmax=np,flush=TRUE,na.strings=c("I/OTimeout","ArcOff-line")))
>>
>> which results in:
>>
>>> dat[1:5,]
>>              t     f    p  d  o   s    a  l c
>> 1 1/21/09 5:01 16151  8.2 76 30 282 1060 53 7
>> 2 1/21/09 5:02 16256  8.3 76 23 282 1059 54 7
>> 3 1/21/09 5:03 16150  8.4 76 26 282 1059 55 7
>> 4 1/21/09 5:04 16150  9.0 76 25 282 1051 57 6
>> 5 1/21/09 5:05 15543 10.4 76  7 282 1024 58 6
>>
>> I have been unable to find a way to convert this into a time series. I did
>> read the manuals and came across a way to coerce a data frame to a ts
>> object: as.ts()
>>
>> Trouble is I do not know how to keep the timestamps in column t in the
>> data frame above. The t column is not strings. If I do:
>>
>> plot.ts(dat)
>>
>> I can see how the first graphics panel is indeed numbers not text. So I
>> think scan converted the text correctly per the format string I put in.
>>
>> Much more difficult still. The datafiles I have contain invalid data,
>> missing values and other none relevant information. I filter this out
>> using subset which works brilliantly. However, how can I filter using
>> subset and convert to a time series afterwards. Since after subsetting
>> there will be 'holes' i.e. missing records. Can a ts object deal with
>> missing records? If so, how? Just point me to a document. I can and will
>> put in the work to figure it out myself.
>>
>> Thank you!
>> Alex van der Spek
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/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.
>



-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] binary version of R 2.8.x

[stephen sefick]

Does solaris have a package managment system?  And do you have all of
the packages for building?  It looks like you do, but I am not sure.

On Wed, Apr 8, 2009 at 12:42 PM,   wrote:
>
> Yes, ./configure works with no error.
> I guess it's a problem of my environment setup.  It's my first attempt to 
> build it myself, and I have downloaded the binaries and installed it 
> previously.
>
>
> --- On Wed, 4/8/09, stephen sefick  wrote:
>
>> From: stephen sefick 
>> Subject: Re: [R] binary version of R 2.8.x
>> To: tomkur2006-takeh...@yahoo.com
>> Cc: r-help@r-project.org
>> Date: Wednesday, April 8, 2009, 7:35 AM
>> ./configure works with no errors ?
>> I have compiled on mac os x and debian so my advice may be
>> circumspect, but the errors I don't know what they are.
>>
>> On Tue, Apr 7, 2009 at 10:17 PM,
>>  wrote:
>> >
>> >
>> > Hello,
>> >
>> > Is there any binary version of R 2.8.x available for
>> x86 Solaris 10?  I need it so that I can use some R library
>> files.  I actually tried to compile the source myself, but
>> it failed.  I don't know much about how to compile
>> source files on Solaris, and it will be great if there is a
>> binary version of R 2.8.x available.  Thanks.
>> >
>> > /***error after running make***/
>> > util.c: In function `Rf_mbrtowc':
>> > util.c:1147: parse error before `char'
>> > util.c:1150: `p' undeclared (first use in this
>> function)
>> > util.c:1150: (Each undeclared identifier is reported
>> only once
>> > util.c:1150: for each function it appears in.)
>> > util.c:1150: `q' undeclared (first use in this
>> function)
>> > util.c:1150: `err' undeclared (first use in this
>> function)
>> > make[3]: *** [util.o] Error 1
>> > make[3]: Leaving directory
>> `/var/tmp/R-2.8.1/src/main'
>> > make[2]: *** [R] Error 2
>> > make[2]: Leaving directory
>> `/var/tmp/R-2.8.1/src/main'
>> > make[1]: *** [R] Error 1
>> > make[1]: Leaving directory `/var/tmp/R-2.8.1/src'
>> > make: *** [R] Error 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.
>> >
>>
>>
>>
>> --
>> Stephen Sefick
>>
>> Let's not spend our time and resources thinking about
>> things that are
>> so little or so large that all they really do for us is
>> puff us up and
>> make us feel like gods.  We are mammals, and have not
>> exhausted the
>> annoying little problems of being mammals.
>>
>>                                                               -K. Mullis
>



-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Genstat into R - Randomisation test

[Mike Lawrence]

I was taught that Fisher proposed the F-test as a computationally
simpler approximation to what he called a "Randomization test",
consisting of exhaustive permutations. I never looked at the original
Fisher reference myself, so this may be false.

However, I haven't observed a consistent nomenclature when I have seen
these tests discussed, so I typically ensure to mention whether what
I'm doing is exhaustive or non-exhaustive.

I do see the value in your interpretation, and think it makes sense to
drop "randomization"  as a name (despite it's possible historical
significance) and start using "exhaustive permutation test" (to
contrast with "non-exhaustive permutation test").


On Wed, Apr 8, 2009 at 1:18 PM, Peter Dalgaard  wrote:
> Mike Lawrence wrote:
>>
>> Looks like that code implements a non-exhaustive variant of the
>> randomization test, sometimes called a permutation test.
>
> Isn't it the other way around? (Permutation tests can be exhaustive by
> looking at all permutations, if a randomization test did that, then it
> wouldn't be random.)
>
>
> --
>   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
> ~~ - (p.dalga...@biostat.ku.dk)              FAX: (+45) 35327907
>



-- 
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University

Looking to arrange a meeting? Check my public calendar:
http://tinyurl.com/mikes-public-calendar

~ Certainty is folly... I think. ~

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Minimum Spanning Tree

[jpearl01]

>I am not sure what you mean. Of course you can plot it using different
>layouts, e.g. with layout.reingold.tilford (after choosing the root
>vertex in some way) and then it looks like a usual tree plot, but why
>would that be any better?

I'd like to be able to distinguish between the nodes better.  For example
the image that I get looks like:
http://www.nabble.com/file/p22954099/mst.gif mst.gif 

The tree nodes are just too close together to be able to read.  I think the
problem might be that some distances between the nodes are *much* smaller
than other edges.The graph does not need to be directed, since the
distance metric does not imply a direction. The most important part is just
being able to see and differentiate between the nodes, and right now they
seem to be all lumped together.  How can I make it more readable?

Thanks,
~josh


-- 
View this message in context: 
http://www.nabble.com/Minimum-Spanning-Tree-tp22934813p22954099.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] MLE for bimodal distribution

[_nico_]

Hello everyone,

I'm trying to use mle from package stats4 to fit a bi/multi-modal
distribution to some data, but I have some problems with it.
Here's what I'm doing (for a bimodal distribution):

# Build some fake binormally distributed data, the procedure fails also with
real data, so the problem isn't here
data = c(rnorm(1000, 3, 0.5), rnorm(500, 5, 0.3))
# Just to check it's bimodal
plot(density(data))
f = function(m, s, m2, s2, w)
{
-log( sum(w*dnorm(data, mean=m, sd=s)) + sum((1-w)*dnorm(data, mean=m2,
sd=s2)) )
}

res = mle(f, start=list("m"=3, "s"=0.5, "m2"=5, "s2"=0.35, "w"=0.6))

This gives an error:
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  non-finite finite-difference value [2]
In addition: There were 50 or more warnings (use warnings() to see the first
50)
And the warnings are about dnorm producing NaNs

So, my questions are:
1) What does "non-finite finite-difference value" mean? My guess would be
that an Inf was passed somewhere where a finite number was expected...? I'm
not sure how optim works though...
2) Is the log likelihood function I wrote correct? Can I use the same type
of function for 3 modes?
3) What should I do to avoid function failure? I tried by changing the
parameters, but it did not work.
4) Can I put constraints to the parameters? In particular, I would like w to
be between 0 and 1.
5) Is there a way to get the log-likelihood value, so that I can compare
different extimations?
6) Do you have any (possibly a bit "pratically oriented") readings about MLE
to suggest?

Thanks in advance
Nico
-- 
View this message in context: 
http://www.nabble.com/MLE-for-bimodal-distribution-tp22954970p22954970.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] subscript out of bounds in eqscplot problem

[Pierre Moffard]

Dear R users,

I have the following problem. Suppose I have the following toy data set:
> data
 m1 m2 m3 m4 m5 state
[1,]  1  0  1 13 23   2
[2,]  0  1  0 23 94   2
[3,]  1  0  0 45 56   1
[4,]  0  1  0 35 84   2
[5,]  1  1  0 98 37   1
[6,]  1  1  0 68  1   2

where the last column is a categorical variable representing the outcome, 
columns m1 to m3 are binary and columns m4 and m5 are continuous. I wish to 
perform lda and then plot the results so I do:

sc.dat<-scale(data[,-6])
sc.dat<-cbind(sc.dat,data[,6])
try<-lda(sc.dat[,1:5],sc..dat[,6])
tr.pred<-predict(try)
eqscplot(tr.pred$x)

and then I get the following error (I also have collinearity which I ignore 
because it doesnt matter):

> eqscplot(tr.pred$x)
Error in eqscplot(tr.pred$x) : subscript out of bounds

Why am I getting this error?  What does it mean?

Any insight is highly appreciated.

Thanks in advance,
Pierre


  
[[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] binary version of R 2.8.x

[tomkur2006-takehome]

Yes, ./configure works with no error.
I guess it's a problem of my environment setup.  It's my first attempt to build 
it myself, and I have downloaded the binaries and installed it previously.


--- On Wed, 4/8/09, stephen sefick  wrote:

> From: stephen sefick 
> Subject: Re: [R] binary version of R 2.8.x
> To: tomkur2006-takeh...@yahoo.com
> Cc: r-help@r-project.org
> Date: Wednesday, April 8, 2009, 7:35 AM
> ./configure works with no errors ?
> I have compiled on mac os x and debian so my advice may be
> circumspect, but the errors I don't know what they are.
> 
> On Tue, Apr 7, 2009 at 10:17 PM, 
>  wrote:
> >
> >
> > Hello,
> >
> > Is there any binary version of R 2.8.x available for
> x86 Solaris 10?  I need it so that I can use some R library
> files.  I actually tried to compile the source myself, but
> it failed.  I don't know much about how to compile
> source files on Solaris, and it will be great if there is a
> binary version of R 2.8.x available.  Thanks.
> >
> > /***error after running make***/
> > util.c: In function `Rf_mbrtowc':
> > util.c:1147: parse error before `char'
> > util.c:1150: `p' undeclared (first use in this
> function)
> > util.c:1150: (Each undeclared identifier is reported
> only once
> > util.c:1150: for each function it appears in.)
> > util.c:1150: `q' undeclared (first use in this
> function)
> > util.c:1150: `err' undeclared (first use in this
> function)
> > make[3]: *** [util.o] Error 1
> > make[3]: Leaving directory
> `/var/tmp/R-2.8.1/src/main'
> > make[2]: *** [R] Error 2
> > make[2]: Leaving directory
> `/var/tmp/R-2.8.1/src/main'
> > make[1]: *** [R] Error 1
> > make[1]: Leaving directory `/var/tmp/R-2.8.1/src'
> > make: *** [R] Error 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.
> >
> 
> 
> 
> -- 
> Stephen Sefick
> 
> Let's not spend our time and resources thinking about
> things that are
> so little or so large that all they really do for us is
> puff us up and
> make us feel like gods.  We are mammals, and have not
> exhausted the
> annoying little problems of being mammals.
> 
>   -K. Mullis

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Constrained, multiple response statistics

[Gene Leynes]

This sounds very similar to what I've been working on, but I'm not sure
without an example.

My solution has been to use an optimization that normalizes inside the
objective function.  The betas that are provided by optim are not
normalized, however since they were normalized inside the objective
function, normalizing them after the fact mirrors the internal workings of
the objective function.

see this example:
http://markmail.org/message/ze5237m6gbgvvvyf

Still, after looking at several statistical packages, and considering the
thoughtful responses from my post, I think that there must be a better way
using existing models, so I've been looking at other packages / models.

On Tue, Apr 7, 2009 at 6:10 PM, Jonathan Greenberg wrote:

> R'ers:
>
>   I was hoping I could get some direction on this.  I have a dataset of the
> form:
>
> Y1,Y2,...,YM = f(X1,X2,...,XN), where N is >>> M
>
> The response data (Y1,Y2,...,YM) is frequency data, such that the sum of
> all Yi = 1.0.  Both Xj and Yi are continuous variables.
>
> I'm trying to figure out the best approach(es) to solving for the model f()
> -- any ideas?  I could solve each Y one at a time, but the lack of
> constraint worries me, and I'm pretty sure that normalizing the data
> afterwards to sum to 1.0 is not going to work out properly.  Thoughts?  I've
> never worked with multiple response statistics before, so I'm mostly trying
> to get some pointers on where to begin investigating...
>
> --j
>
> --
>
> Jonathan A. Greenberg, PhD
> Postdoctoral Scholar
> Center for Spatial Technologies and Remote Sensing (CSTARS)
> University of California, Davis
> One Shields Avenue
> The Barn, Room 250N
> Davis, CA 95616
> Cell: 415-794-5043
> AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

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

Re: [R] axis values on lattice log-scale plot

[Deepayan Sarkar]

On 4/7/09, Don McKenzie  wrote:
> I'm plotting the following (stripped of inessentials)
>
>  xyplot(sd ~ distance |
> wshed,data=sdvar.df,scales=list(x=list(log=TRUE),y=list(log=TRUE)))
>
>  sdvar.df  is a data frame, sd and distance are numeric, wshed is an ordered
> factor
>
>  trying to replicate the action of  log="xy" in plot()
>
>  The plot works fine but the axis values at the ticks are in scientific
> notation,
>  e.g. 10^1.5, 10^2.0, ... on the X axis and 10^-.5, 10^-.4, ... on the Y
> axis.
>
>  I've searched the lattice help (and I know it must be in there:-()  and
> the archive but can't seem to find a way to make
>  the axis values display as simply numeric  (i.e., 10,100,1000, ... on the X
> axis and proportions on the Y axis).
>
>  Suggestions much appreciated.

You need to provide your own axis annotation function.
?xscale.components.default gives you details and an example. Also, see
Figure 8.4 and 8.5 in

http://lmdvr.r-forge.r-project.org/

-Deepayan

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

[R] Comparing Proportions Among Groups

[Isabella Ghement]

Hi everyone,

I am trying to compare proportions among groups using the logistic
regression
approach as follows:

1) Fit the model log(p_i/(1-p_i)) = M + G_i, where p_i is the probability
   of success in group i and G_i is the effect of group i, i=1,..,I.
2) Test the hypotheses:
   Ho: G_1 = G_2 = ... = G_I (the probability of success is the same for all
groups)
   versus
   Ha: at least two G_i's are different (at least two groups have different
probabilities of success)

I am wondering if there any R functions available for accomplishing these
tasks.
Using dummy variables is one option, but is it possible to solve this
problem in R
without resorting to the use of dummy variables?

Many thanks and kind regards,

Isabella

Isabella R. Ghement, Ph.D.
Ghement Statistical Consulting Company
301-7031 Blundell Road, Richmond, B.C., Canada, V6Y 1J5
Tel: 604-767-1250
Fax: 604-270-3922
E-mail: isabe...@ghement.ca
Web: www.ghement.ca

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] vectors on top of contours, and lattice

[Deepayan Sarkar]

On 4/8/09, Cable, Samuel B Civ USAF AFMC AFRL/RVBXI
 wrote:
> OK, I needed to plot a set of vectors on top of a contour plot.  I
>  figured out a way to do this.  I create a panel function that calls
>  "larrows()" with arguments constructed from my vector data.  Then, when
>  I go to do the contour plot, I call contourplot() with the "panel"
>  argument set to point to my newly created panel function.
>
>
>
>  So far, so good.
>
>
>
>  Now, I need to do a series of contour plots, each overlaid with a
>  different set of vectors.  I can call contourplot(z~x*y,...)  but now
>  I'm stuck, because each contour plot will be overlaid with the *same*
>  vector plot.  This won't work because each set of contoured data has its
>  own particular set of associated vector data.  (BTW, I am not doing
>  gradients here.  The vector data is its own entity, not derivable from
>  the contoured data.)  Can I somehow define a whole series of panel
>  functions and coax contourplot() to choose the right one for each of its
>  plots?  Or do I have to revise my entire approach here?

See ?packet.number for one approach. This should allow you to extract
the correct component of the vector data inside the panel function.

-Deepayan

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

Re: [R] Colour each letter of a text string in a plot

[jim holtman]

Use 'text' to write out each one:

plot(0, type='n', xlim=c(0,1), ylim=c(0,1))
text(seq(0,1,length=10), rep(0.5,10), LETTERS[1:10], col=1:10)



On Wed, Apr 8, 2009 at 12:15 PM, Daren Tan  wrote:
> I am inserting a DNA sequence into a plot, and hope to colourize each
> of the four nucleotide of the DNA sequence with a unique colour i.e.,
> A ("red"), C ("green"), G ("blue", and T ("yellow"). I use the
> following codes, but the DNA sequence only shows as "red"
>
>
> DNA <- "ACGT"
>
> plot(1, xlim = c(0,1), ylim = c(0,1), axes=F, xlab="", ylab="",
> main="", xaxs="i", yaxs="i")
>
> text(x=0.5, y=0.5, labels=DNA, col=c("red", "green", "blue", "yellow"))
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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 that 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.

Re: [R] Tinn-R pdf()

[Deepayan Sarkar]

On 4/8/09, Tobias Verbeke  wrote:
> Hi Henning,
>
>
> > thanks for your help, with solved the problem, although i don't why,
> because when using the R editor accessible via the R console i created many
> many lattice plots with the code i posted, i.e. without the print() command.
> >
>
>  At the command line, R objects (including lattice plots which are
>  objects) get auto-printed, i.e. the print method is invoked automatically
> on these objects.
>
>  This is not the case when you write to a pdf file.

That claim is not accurate. Auto-printing should not be affected by
whether a pdf device is open at the moment.

I'm not familiar with TINN-R, but it's possible it executes R commands
through some sort of wrapper function, and that is failing to emulate
auto-printing correctly.

-Deepayan

>
>  Best,
>  Tobias
>
>
>
> >
> > >  Original-Nachricht 
> > > Datum: Wed, 8 Apr 2009 12:29:58 +0200
> > > Von: "ONKELINX, Thierry" 
> > > An: "Henning Wildhagen" , r-help@r-project.org
> > > Betreff: RE: [R] Tinn-R pdf()
> > >
> > >Dear Henning,
> > >
> > > You need to print() lattice plots when using a device:
> > >  library(lattice)
> > >
> > > pdf("plot1.pdf")
> > > PLOT<-(xyplot, ...)
> > > print(PLOT)
> > > dev.off()
> > >
> > > So this is not due to TINN-R.
> > >
> > > HTH,
> > >
> > > Thierry
> > >
> > >
> 
> > > 
> > > ir. Thierry Onkelinx
> > > Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> > > and Forest
> > > Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
> > > methodology and quality assurance
> > > Gaverstraat 4
> > > 9500 Geraardsbergen
> > > Belgium tel. + 32 54/436 185
> > > thierry.onkel...@inbo.be www.inbo.be
> > > To call in the statistician after the experiment is done may be no more
> > > than asking him to perform a post-mortem examination: he may be able to
> > > say what the experiment died of.
> > > ~ Sir Ronald Aylmer Fisher
> > >
> > > The plural of anecdote is not data.
> > > ~ Roger Brinner
> > >
> > > The combination of some data and an aching desire for an answer does not
> > > ensure that a reasonable answer can be extracted from a given body of
> > > data.
> > > ~ John Tukey
> > >
> > > -Oorspronkelijk bericht-
> > > Van: r-help-boun...@r-project.org
> [mailto:r-help-boun...@r-project.org]
> > > Namens Henning Wildhagen
> > > Verzonden: woensdag 8 april 2009 11:25
> > > Aan: r-help@r-project.org
> > > Onderwerp: [R] Tinn-R pdf()
> > >
> > > Dear R and Tinn-R users,
> > >
> > > i recently switched to Tinn-R and sending code to R works fine (R 2.8.1,
> > >
> > > Tinn-R 2.2.0.2, OS Windows XP). However, i encountered a problem when
> trying to send plots to pdf files like this:
> > >
> > > library(lattice)
> > >
> > > pdf("plot1.pdf")
> > > PLOT<-(xyplot, ...)
> > > PLOT
> > > dev.off()
> > >
> > > The file "plot1.pdf" is created, but it is empty.
> > > If i paste the code above directly into the R console and run it, the
> > > file "plot1.pdf" is created and in this case contains "PLOT".
> > >
> > > I guess that some settings in Tinn-R are wrong, but i have no idea
> > > which. Maybe someone has a suggestion?
> > >
> > > Thanks,
> > >
> > > Henning
> > > --

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Genstat into R - Randomisation test

[Peter Dalgaard]

Mike Lawrence wrote:

Looks like that code implements a non-exhaustive variant of the
randomization test, sometimes called a permutation test. 


Isn't it the other way around? (Permutation tests can be exhaustive by 
looking at all permutations, if a randomization test did that, then it 
wouldn't be random.)



--
   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
~~ - (p.dalga...@biostat.ku.dk)  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.

[R] Colour each letter of a text string in a plot

[Daren Tan]

I am inserting a DNA sequence into a plot, and hope to colourize each
of the four nucleotide of the DNA sequence with a unique colour i.e.,
A ("red"), C ("green"), G ("blue", and T ("yellow"). I use the
following codes, but the DNA sequence only shows as "red"


DNA <- "ACGT"

plot(1, xlim = c(0,1), ylim = c(0,1), axes=F, xlab="", ylab="",
main="", xaxs="i", yaxs="i")

text(x=0.5, y=0.5, labels=DNA, col=c("red", "green", "blue", "yellow"))

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Col Names

[jim holtman]

Is this what you are looking for:

> x
  Positions  DEV.MSCI
04/30/1980 00:00:00.000 -0.0150328542
05/31/1980 00:00:00.000  0.0005087752
06/30/1980 00:00:00.000  0.0586794492
07/31/1980 00:00:00.000  0.0458505592
08/31/1980 00:00:00.000  0.0350926728
> colnames(x)[2] <- "NewName"
> x
  Positions   NewName
04/30/1980 00:00:00.000 -0.0150328542
05/31/1980 00:00:00.000  0.0005087752
06/30/1980 00:00:00.000  0.0586794492
07/31/1980 00:00:00.000  0.0458505592
08/31/1980 00:00:00.000  0.0350926728


On Wed, Apr 8, 2009 at 10:45 AM, livia  wrote:
>
> Hello,
>
> I am working with S-plus TimeSeries Object. I wonder how can I change the
> column names of the variable instead of using the one default?
>
> i.e to change "DEV.MSCI" to other name
>              Positions            DEV.MSCI
>  04/30/1980 00:00:00.000 -0.0150328542
>  05/31/1980 00:00:00.000  0.0005087752
>  06/30/1980 00:00:00.000  0.0586794492
>  07/31/1980 00:00:00.000  0.0458505592
>  08/31/1980 00:00:00.000  0.0350926728
>
> --
> View this message in context: 
> http://www.nabble.com/Col-Names-tp22952055p22952055.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.
>



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

What is the problem that 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] Sweave problem, with multicolumn tables from R to LaTeX

[Christian Salas]

Hi there,

I have been using the example provided bellow for a while, and It was 
working without any problem. Nevertheless, just since 2-3 days is not 
working, probably because I did update.packages(). I have tried to 
re-install the older versions of the packages Hmisc() and xtable(), but 
still does not work. Can you run this example, and tell me if you got 
the same problems?


I use linux-ubuntu, and R 2.8.1

Here is my .Rnw file

%%%
\documentclass{article}
\usepackage[margin=1in]{geometry}
\usepackage{Sweave}
\SweaveOpts{echo=TRUE}

\title{Why is not working this? a table from R having multi columns
names}
\author{Myself}
\date{\today}

\begin{document}
\maketitle

<>=
#Calling some packages
library(Hmisc)
library(xtable)

#version of the packages
packageDescription('Hmisc')$Version #version of Hmisc
packageDescription('xtable')$Version  #version of xtable

#version and features of my PC and R
R.version$platform
R.version$version.string
@

\section{The table ready to be exported}
An R table as an example.
<<>>=
tab.coef.allmodels <-
structure(c(246.8809, 153.6952, 43.6259, 1823.686, -4.0138, -53.8741,
-13.9552, -699.7154, 0, -3.0084, -0.5188, -41.7249, 4.0898,
2.0342, 0.2798, 25.636, 235.5128, 160.0244, 30.0147, 1523.0723,
3.7601, -53.9238, -3.9789, -597.1446, 0, -3.2745, -0.0025, -28.7418,
4.4203, 2.1606, 0.0124, 17.9259), .Dim = c(4L, 8L), .Dimnames = list(
c("NHA", "GHA", "Hdom", "VHA"), c("Intercept", "stand2",
"band3", "band7", "Intercept", "stand2", "band3", "band7"
)))
tab.coef.allmodels
@

\section{Producing the \LaTeX\ table}
\subsection{Showing a simple output table, but not being
  the aimed one}

The following syntax produce Table \ref{tab:coefAllmod1}, which is not
what I want.
<>=
l <- latex(tab.coef.allmodels, file="tab.coefAllmod1.tex",
rowlabel="Model for",
  first.hline.double=FALSE,
  rowlabel.just="l",
caption="MLR and LME parameter estimates by stand variables",
  caption.loc="top", label="tab:coefAllmod1",where="!h")
@

\input{tab.coefAllmod1}


\subsection{What I want, but not getting it right}
The following syntax produce Table \ref{tab:coefAllmod2}, which
contains all the information that i want, but is not giving me
the correct format.
<>=
l <- latex(tab.coef.allmodels, file="tab.coefAllmod2.tex",
rowlabel="Model for",
  first.hline.double=FALSE,
  rowlabel.just="l",
 cgroup=c("MLR","LME"),
 n.cgroup=c(4,4),
caption="MLR and LME parameter estimates by stand variables",
  caption.loc="top", label="tab:coefAllmod2",where="!h")
@

\input{tab.coefAllmod2}

I have been using this
syntax for a while, and was working without any problem.
 Nevertheless, just since 2-3 days is not working, probably
  because I did {\tt update.packages()}.
 I have tried to re-install the older versions of the packages
 {\tt Hmisc()} and {\tt xtable()}, but still does not work.

\vspace{1in}
{\Large   {\bf What Am I missing here?}}

\end{document}

%%%

The the main .tex file, tables, the ouput pdf, and the .Rnw files are in 
the following url, if you want to check them


http://environment.yale.edu/salas/questions/shortTable/multcolTable_short.tex
http://environment.yale.edu/salas/questions/shortTable/tab.coefAllmod1.tex
http://environment.yale.edu/salas/questions/shortTable/tab.coefAllmod2.tex
http://environment.yale.edu/salas/questions/shortTable/multcolTable_short.pdf
http://environment.yale.edu/salas/questions/shortTable/multcolTable_short.Rnw

thanks in advance for your help

--
---
Christian Salas E-mail:christian.sa...@yale.edu
PhD candidate   http://environment.yale.edu/salas
School of Forestry and Environmental Studies
Yale University Tel: +1-(203)-432 5126
360 Prospect St Fax:+1-(203)-432 3809
New Haven, CT 06511-2189Office: Room 35, Marsh Hall
USA

Yale Biometrics Lab  http://environment.yale.edu/biometrics

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Is a point into an ellipse

[Duncan Murdoch]

On 4/8/2009 11:17 AM, Alain Guillet wrote:

Hi,

I drew an ellipse with the package ellipse. Now I would like to know  if 
a point is inside the ellipse. Is any R functions to do it without 
computing the equation of the ellipse manually? Thanks.


For example, if I do "plot(ellipse(0.8), type = 'l')", I would like to 
know if (0,1) belongs to the drawn ellipse.


It's probably too late by then, but if you saved the result of the call 
to ellipse(0.8), and treat it as a polygon, there is the function 
point.in.polygon in the sp package, or pip in splancs, or inside in 
RFOC, or inside in GEOmap, etc.


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] Reshape - strange outputs

[Steve Murray]

Dear R Users,

I am using the reshape package to reformat gridded data into column format 
using the code shown below. However, when I display the resulting object, a 
single column is fomed (instead of three) and all the latitude values (which 
should be in either column one or two) are collected at the bottom. Also, the 
NA values aren't removed, despite this being requested in the code.

Code:


# NetCDF file has been read in and is being processed...

arunoff_1986_temp <- get.var.ncdf(netcdf_1036_temp, "arunoff")


# Assign row and column names

columnnames <- sprintf("%.2f", seq(from = -89.75, to = 89.75, length = 360))
rnames <- sprintf("%.2f", seq(from = -179.75, to = 179.75, length = 720))

colnames(arunoff_1986_temp) <- columnnames
rownames(arunoff_1986_temp) <- rnames


# Melt into columnar format

arunoff_1986$Longitude <- rownames(arunoff_1986)
# Note: If I do: arunoff_1986$Latitude <- rownames(arunoff_1986)  (i.e. 
change it to 'Latitude', I get the following: Warning message: In 
arunoff_1986$Latitude <- rownames(arunoff_1986) : Coercing LHS to a list ), 
thus proceed using 'Longitude' where no warning is apparent.


arunoff_long_1986 <- melt(arunoff_1986, id.var="Longitude", na.rm=TRUE)


> dim(arunoff_long_1986)
[1] 259560  2

>head(arunoff_long_1986, n=10)  # This displays what looks like one single 
>column, but is in fact two: column entitled "L1" is empty until the end of the 
>file, as shown by the 'tail' command below.
   value L1
1  
2  
3  
4  
5  
6  
7  
8  
9  
10 

> tail(arunoff_long_1986, n=10)
 value   L1
Latitude.351 85.25 Latitude
Latitude.352 85.75 Latitude
Latitude.353 86.25 Latitude
Latitude.354 86.75 Latitude
Latitude.355 87.25 Latitude
Latitude.356 87.75 Latitude
Latitude.357 88.25 Latitude
Latitude.358 88.75 Latitude
Latitude.359 89.25 Latitude
Latitude.360 89.75 Latitude



I'd be very grateful indeed if anyone is able to offer assistance by way of 
pointing out what I've done wrong. I've spent a long time working on this and 
trying various options, but am srill none the wiser! I'm aiming for three 
columns: Latitude, Longitude, 'value'. 

Many thanks for any help offered,

Steve

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

[R] problem with creating a netcdf file under script sh

[tuatu]

Hi everyone, 
I try to make a netcdf file which disposes a difference between 2 variables
of 2netcdf files the same dimension 
When I programmed under R, everything is ok but when I put the code under
EOF of a sh script, an error occurs: 
"Error, passed variable has a dim that is NOT of class dim.ncdf!". My
programme is: 

x1<-open.ncdf("file1.nc") 
x2<-open.ncdf("file2.nc") 
y1<-get.var.ncdf(x1,"a") 
y2<-get.var.ncdf(x2,"a") 
y<- y1-y2 
t<-var.def.ncdf("t",units="c",dim=x1$dim,-1) 
nc<-create.ncdf("test.nc",vars=t) 
put.var.ncdf(nc,t,y) 
close.ncdf(nc) 

Thanks in advance, 
T.A 

-- 
View this message in context: 
http://www.nabble.com/problem-with-creating-a-netcdf-file-under-script-sh-tp22952125p22952125.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] Col Names

[livia]

Hello,

I am working with S-plus TimeSeries Object. I wonder how can I change the
column names of the variable instead of using the one default?

i.e to change "DEV.MSCI" to other name
  PositionsDEV.MSCI 
 04/30/1980 00:00:00.000 -0.0150328542
 05/31/1980 00:00:00.000  0.0005087752
 06/30/1980 00:00:00.000  0.0586794492
 07/31/1980 00:00:00.000  0.0458505592
 08/31/1980 00:00:00.000  0.0350926728

-- 
View this message in context: 
http://www.nabble.com/Col-Names-tp22952055p22952055.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] persp3d and rgl.viewpoint for rotating 3D plots

[Duncan Murdoch]

On 4/8/2009 11:47 AM, Camarda, Carlo Giovanni wrote:

Dear Duncan,

	thanks for your prompt reply. 
I was tempted to use movie3d, but below I copied what I get on my PC. Then, I have ImageMagick installed and, though I set convert equal to FALSE, I was not able to find where the function writes the .png files.
Maybe I lack experience in managing such complex functions, but I must confess that I have quite some hard time to grasp info from help(movie3d). 


Presumably ImageMagick is not on your PATH, which is why movie3d() 
failed to find it.  It writes the image files to the directory given by 
the dir parameter, which defaults to tempdir().  tempdir() returns a 
session-specific temporary directory, which R will try to clear at the 
end of a session.


Duncan Murdoch



Thanks in advance for any help,
Carlo Giovanni



library(rgl)
# set data
f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
x <- seq(-10, 10, length= 30); y <- x
z <- outer(x, y, f)
# set working directory
setwd("U:\\movie")
dir() # check

character(0)

# plot
persp3d(x, y, z)
movie3d(spin3d(), duration=5, convert=TRUE)
Writing movie000.png 
[...]
Writing movie050.png 

Error in movie3d(spin3d(), duration = 5, convert = TRUE) : 
  ImageMagick not found

dir() # check again

character(0)

movie3d(spin3d(), duration=5, convert=FALSE)
Writing movie000.png 
[...]
Writing movie050.png 


dir() # check again ??

character(0)





version
   _   
platform   i386-pc-mingw32 
arch   i386
os mingw32 
system i386, mingw32   
status 
major  2   
minor  8.1 
year   2008
month  12  
day22  
svn rev47281   
language   R   
version.string R version 2.8.1 (2008-12-22)



-Original Message-
From: Duncan Murdoch [mailto:murd...@stats.uwo.ca] 
Sent: Wednesday, April 08, 2009 5:13 PM

To: Camarda, Carlo Giovanni
Cc: r-h...@stat.math.ethz.ch
Subject: Re: [R] persp3d and rgl.viewpoint for rotating 3D plots

On 4/8/2009 11:01 AM, Camarda, Carlo Giovanni wrote:

Dear R-users,

within the rgl-package, I would have a question about the usage of persp3d in combination of rgl.viewpoint. 
I am not able to figure out how to let a 3D plot rotating around likewise the example in ?rgl.viewpoint. It seems that when I use persp3d(...) I see something on my screen, which is different from what I get when it's rotating. Is there any different behavior between shade3d and persp3d?

Please find below a simple example.

Actually, my final aim is to save each of the rotating graphs for creating a 
"clip". I commented below what I would use later on, but is there any way to do 
it better?


Take a look at ?movie3d.  I think it does what you want.

Duncan Murdoch




Thanks a lot for your help,
Carlo Giovanni Camarda


library(rgl)
# simplified version from help(rgl.viewpoint)
shade3d(oh3d())
coo <- 1:360
for(i in 1:length(coo)) {
rgl.viewpoint(coo[i])
#filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
#rgl.snapshot(filename)
}

# simple persp3d adaption
f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
x <- seq(-10, 10, length= 30); y <- x
z <- outer(x, y, f)
persp3d(x, y, z)
coo <- 1:360
for(i in 1:length(coo)) {
rgl.viewpoint(coo[i])
#filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
#rgl.snapshot(filename)
}



-
Dr. Carlo Giovanni Camarda
Research Scientist
Max Planck Institute for Demographic Research
Laboratory of Statistical Demography
Konrad-Zuse-Straße 1
18057 Rostock - Germany
Phone: +49 (0)381 2081 172
Fax: +49 (0)381 2081 472
cama...@demogr.mpg.de
http://www.demogr.mpg.de/en/staff/camarda/default.htm
-


--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}





__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 mail has been sent through the MPI for Demographic Research.  Should you 
receive a mail that is apparently from a MPI user without this text displayed, 
then the address has most likely been faked. If you are uncertain about the 
validity of this message, please check the mail header or ask your system 
administrator for assistance.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/lis

Re: [R] Convert data frame containing time stamps to time series

[amvds]

Converting dates is getting stranger still. I am coercing a data frame
into a ts as follows:


tst1<-as.POSIXct("1/21/09 5:01",format="%m/%d/%y %H:%M")
tst2<-as.POSIXct("1/28/09 3:40",format="%m/%d/%y %H:%M")
tsdat<-as.ts(dat,start=tst1,end=tst2,frequency=1)

This generates a ts object. But strangely enough the first column of that
matrix starts at the numeric value of 841 counts up to 1139 and then
starts at 1 again, only to count up from there. The restart at 1 occurs at
the first day "1/21/09" at 10:00:00.

What is so special about that time? This phenomenon happens several times
in the long file. But the restart count is always a different number.
This creates a ramp with some bumps.

Can anybody explain this?
Thanks in advance,
Alex van der Spek


> I read records using scan:
>
> dat<-data.frame(scan(file="KDA.csv",what=list(t="%m/%d/%y
> %H:%M",f=0,p=0,d=0,o=0,s=0,a=0,l=0,c=0),skip=2,sep=",",nmax=np,flush=TRUE,na.strings=c("I/OTimeout","ArcOff-line")))
>
> which results in:
>
>> dat[1:5,]
>  t fp  d  o   sa  l c
> 1 1/21/09 5:01 16151  8.2 76 30 282 1060 53 7
> 2 1/21/09 5:02 16256  8.3 76 23 282 1059 54 7
> 3 1/21/09 5:03 16150  8.4 76 26 282 1059 55 7
> 4 1/21/09 5:04 16150  9.0 76 25 282 1051 57 6
> 5 1/21/09 5:05 15543 10.4 76  7 282 1024 58 6
>
> I have been unable to find a way to convert this into a time series. I did
> read the manuals and came across a way to coerce a data frame to a ts
> object: as.ts()
>
> Trouble is I do not know how to keep the timestamps in column t in the
> data frame above. The t column is not strings. If I do:
>
> plot.ts(dat)
>
> I can see how the first graphics panel is indeed numbers not text. So I
> think scan converted the text correctly per the format string I put in.
>
> Much more difficult still. The datafiles I have contain invalid data,
> missing values and other none relevant information. I filter this out
> using subset which works brilliantly. However, how can I filter using
> subset and convert to a time series afterwards. Since after subsetting
> there will be 'holes' i.e. missing records. Can a ts object deal with
> missing records? If so, how? Just point me to a document. I can and will
> put in the work to figure it out myself.
>
> Thank you!
> Alex van der Spek
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] persp3d and rgl.viewpoint for rotating 3D plots

[Camarda, Carlo Giovanni]

Dear Duncan,

thanks for your prompt reply. 
I was tempted to use movie3d, but below I copied what I get on my PC. Then, I 
have ImageMagick installed and, though I set convert equal to FALSE, I was not 
able to find where the function writes the .png files.
Maybe I lack experience in managing such complex functions, but I must confess 
that I have quite some hard time to grasp info from help(movie3d). 

Thanks in advance for any help,
Carlo Giovanni


> library(rgl)
> # set data
> f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
> x <- seq(-10, 10, length= 30); y <- x
> z <- outer(x, y, f)
> # set working directory
> setwd("U:\\movie")
> dir() # check
character(0)
> # plot
> persp3d(x, y, z)
> movie3d(spin3d(), duration=5, convert=TRUE)
Writing movie000.png 
[...]
Writing movie050.png 

Error in movie3d(spin3d(), duration = 5, convert = TRUE) : 
  ImageMagick not found
> dir() # check again
character(0)
> movie3d(spin3d(), duration=5, convert=FALSE)
Writing movie000.png 
[...]
Writing movie050.png 

> dir() # check again ??
character(0)
> 

> version
   _   
platform   i386-pc-mingw32 
arch   i386
os mingw32 
system i386, mingw32   
status 
major  2   
minor  8.1 
year   2008
month  12  
day22  
svn rev47281   
language   R   
version.string R version 2.8.1 (2008-12-22)


-Original Message-
From: Duncan Murdoch [mailto:murd...@stats.uwo.ca] 
Sent: Wednesday, April 08, 2009 5:13 PM
To: Camarda, Carlo Giovanni
Cc: r-h...@stat.math.ethz.ch
Subject: Re: [R] persp3d and rgl.viewpoint for rotating 3D plots

On 4/8/2009 11:01 AM, Camarda, Carlo Giovanni wrote:
> Dear R-users,
> 
> within the rgl-package, I would have a question about the usage of 
> persp3d in combination of rgl.viewpoint. 
> I am not able to figure out how to let a 3D plot rotating around likewise the 
> example in ?rgl.viewpoint. It seems that when I use persp3d(...) I see 
> something on my screen, which is different from what I get when it's 
> rotating. Is there any different behavior between shade3d and persp3d?
> Please find below a simple example.
> 
> Actually, my final aim is to save each of the rotating graphs for creating a 
> "clip". I commented below what I would use later on, but is there any way to 
> do it better?

Take a look at ?movie3d.  I think it does what you want.

Duncan Murdoch


> 
> Thanks a lot for your help,
> Carlo Giovanni Camarda
> 
> 
> library(rgl)
> # simplified version from help(rgl.viewpoint)
> shade3d(oh3d())
> coo <- 1:360
> for(i in 1:length(coo)) {
> rgl.viewpoint(coo[i])
> #filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
> #rgl.snapshot(filename)
> }
> 
> # simple persp3d adaption
> f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
> x <- seq(-10, 10, length= 30); y <- x
> z <- outer(x, y, f)
> persp3d(x, y, z)
> coo <- 1:360
> for(i in 1:length(coo)) {
> rgl.viewpoint(coo[i])
> #filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
> #rgl.snapshot(filename)
> }
> 
> 
> 
> -
> Dr. Carlo Giovanni Camarda
> Research Scientist
> Max Planck Institute for Demographic Research
> Laboratory of Statistical Demography
> Konrad-Zuse-Straße 1
> 18057 Rostock - Germany
> Phone: +49 (0)381 2081 172
> Fax: +49 (0)381 2081 472
> cama...@demogr.mpg.de
> http://www.demogr.mpg.de/en/staff/camarda/default.htm
> -
> 
> 
> --
> This mail has been sent through the MPI for Demographic ...{{dropped:10}}
> 
> 
> 
> 
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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 mail has been sent through the MPI for Demographic Research.  Should you 
receive a mail that is apparently from a MPI user without this text displayed, 
then the address has most likely been faked. If you are uncertain about the 
validity of this message, please check the mail header or ask your system 
administrator for assistance.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Howto Plot With Transparent Background

[Christof Winter]

Eik Vettorazzi wrote, On 08.04.2009 15:08:
By default pdfs have a transparent background, see ?pdf and there the bg 
part,

which can be proven by a minimal R+ LaTeX example

#R code
pdf("test.pdf")
plot(1,1)
dev.off()

#minimal tex example
\documentclass[a4paper,12pt]{article}
\usepackage[svgnames]{xcolor}
\usepackage{graphicx}
\begin{document}
\pagecolor{green}


I really like those bizarre green pages. Thanks for sharing.
:-)

Christof


--
Christof Winter
Bioinformatics Group
Biotechnologisches Zentrum
Technische Universität Dresden
Tatzberg 47-51
01307 Dresden
Germany

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] vectors on top of contours, and lattice

[Cable, Samuel B Civ USAF AFMC AFRL/RVBXI]

OK, I needed to plot a set of vectors on top of a contour plot.  I
figured out a way to do this.  I create a panel function that calls
"larrows()" with arguments constructed from my vector data.  Then, when
I go to do the contour plot, I call contourplot() with the "panel"
argument set to point to my newly created panel function.

 

So far, so good.

 

Now, I need to do a series of contour plots, each overlaid with a
different set of vectors.  I can call contourplot(z~x*y,...)  but now
I'm stuck, because each contour plot will be overlaid with the *same*
vector plot.  This won't work because each set of contoured data has its
own particular set of associated vector data.  (BTW, I am not doing
gradients here.  The vector data is its own entity, not derivable from
the contoured data.)  Can I somehow define a whole series of panel
functions and coax contourplot() to choose the right one for each of its
plots?  Or do I have to revise my entire approach here?

 

Thanks!

 

--Sam

 


[[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] [OT ?] rant (was : Re: Conversions From standard to metric units)

[Michael Dewey]

At 00:48 04/04/2009, Duncan Murdoch wrote:

On 03/04/2009 5:37 PM, Emmanuel Charpentier wrote:

Le vendredi 03 avril 2009 à 14:17 -0400, stephen sefick a écrit :

I am starting to use R for almost any sort of calculation that I need.
 I am a biologist that works in the states, and there is often a need
to convert from standard units to metric units.


US/Imperial units are *not* standard units.


But they are fun:  you should see the arguments 
you can have about whether imperial fluid ounces 
are the same volume as US fluid ounces. (They're 
not: US ounces are bigger.  But not big enough so that their gallons catch up!)


Even late in the day I cannot resist sharing the 
astonishment of a colleague from the United 
States when I told him I had been taught as a child the rhyme

A pint of pure water
Weighs a pound and a quarter

whereas he had been taught
A pint's a pound
The world around




Duncan Murdoch


The former "metric system"

is now called "Système International" (International System) for a
reason, which is *not* gallocentrism of a "few" 6e7 frogs, but rather
laziness of about 5.6e9 losers who refuse to load their memories with
meaningless conversion factors...

Emmanuel Charpentier
who has served his time with
pounds per cubic feet, furlongs
per fortnight, BTU and other
figments of British/American
sadistic imagination, thank you
very much...
 # Again, didn't work the first time...


 Is there a package in
R for this already?  If not I believe that I am going to write some of
the most often used in function form.  My question is should I include
this in my StreamMetabolism package.  It is not along the same theme
lines, but could loosely fit.  The reason that I ask is that I don't
want to clutter CRAN with a small package containing some conversion
functions because I am to lazy to source them into R every time that I
use them, but I also don't want the StreamMetabolism package to turn
into StephenMisc Fuctions.  Thoughts, comments, or suggestions would
be appreciated.

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





Michael Dewey
http://www.aghmed.fsnet.co.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.

[R] predict "interval" for lmRob?

[Galkowski, Jan]

lm's "predict" function offers an "interval" parameter to choose between 
'confidence' and 'prediction' bands. In the package "robust" and for "lmRob", 
there is also a "predict" but it lacks such a parameter, and the documented 
"type" parameter has only "response" offerred.  Is there some way of obtaining 
prediction bands from lmRob?  Is there an alternative robust (linear) 
regression package that offers such a capability?

Thanks for any and all help.

  - Jan Galkowski, Akamai Technologies, Cambridge, MA.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Is a point into an ellipse

[Alain Guillet]

Hi,

I drew an ellipse with the package ellipse. Now I would like to know  if 
a point is inside the ellipse. Is any R functions to do it without 
computing the equation of the ellipse manually? Thanks.


For example, if I do "plot(ellipse(0.8), type = 'l')", I would like to 
know if (0,1) belongs to the drawn ellipse.



Regards,
Alain

--
Alain Guillet
Statistician and Computer Scientist

SMCS - Institut de statistique - Université catholique de Louvain
Bureau d.126
Voie du Roman Pays, 20
B-1348 Louvain-la-Neuve
Belgium

tel: +32 10 47 30 50

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] persp3d and rgl.viewpoint for rotating 3D plots

[Duncan Murdoch]

On 4/8/2009 11:01 AM, Camarda, Carlo Giovanni wrote:

Dear R-users,

within the rgl-package, I would have a question about the usage of persp3d in combination of rgl.viewpoint. 
I am not able to figure out how to let a 3D plot rotating around likewise the example in ?rgl.viewpoint. It seems that when I use persp3d(...) I see something on my screen, which is different from what I get when it's rotating. Is there any different behavior between shade3d and persp3d?

Please find below a simple example.

Actually, my final aim is to save each of the rotating graphs for creating a 
"clip". I commented below what I would use later on, but is there any way to do 
it better?


Take a look at ?movie3d.  I think it does what you want.

Duncan Murdoch




Thanks a lot for your help,
Carlo Giovanni Camarda


library(rgl)
# simplified version from help(rgl.viewpoint)
shade3d(oh3d())
coo <- 1:360
for(i in 1:length(coo)) {
rgl.viewpoint(coo[i])
#filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
#rgl.snapshot(filename)
}

# simple persp3d adaption
f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
x <- seq(-10, 10, length= 30); y <- x
z <- outer(x, y, f)
persp3d(x, y, z)
coo <- 1:360
for(i in 1:length(coo)) {
rgl.viewpoint(coo[i])
#filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
#rgl.snapshot(filename)
}



-
Dr. Carlo Giovanni Camarda
Research Scientist
Max Planck Institute for Demographic Research
Laboratory of Statistical Demography
Konrad-Zuse-Straße 1
18057 Rostock - Germany
Phone: +49 (0)381 2081 172
Fax: +49 (0)381 2081 472
cama...@demogr.mpg.de
http://www.demogr.mpg.de/en/staff/camarda/default.htm
-


--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}





__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] read.spss, locale and encodings

[Hans Ekbrand]

On Wed, Apr 08, 2009 at 04:17:51PM +0200, Peter Dalgaard wrote:
> Hans Ekbrand wrote:
>> Someone running foreign 8.34 that is willing to test my SPSS-file?
>
> Someone with an SPSS file problem willing to help test the prereleases? :-)

http://sociologi.cjb.net/temp/test.sav

-- 
Hans Ekbrand (http://sociologi.cjb.net) 
GPG Fingerprint: 1408 C8D5 1E7D 4C9C C27E 014F 7C2C 872A 7050 614E


signature.asc
Description: Digital signature
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] persp3d and rgl.viewpoint for rotating 3D plots

[Camarda, Carlo Giovanni]

Dear R-users,

within the rgl-package, I would have a question about the usage of persp3d 
in combination of rgl.viewpoint. 
I am not able to figure out how to let a 3D plot rotating around likewise the 
example in ?rgl.viewpoint. It seems that when I use persp3d(...) I see 
something on my screen, which is different from what I get when it's rotating. 
Is there any different behavior between shade3d and persp3d?
Please find below a simple example.

Actually, my final aim is to save each of the rotating graphs for creating a 
"clip". I commented below what I would use later on, but is there any way to do 
it better?

Thanks a lot for your help,
Carlo Giovanni Camarda


library(rgl)
# simplified version from help(rgl.viewpoint)
shade3d(oh3d())
coo <- 1:360
for(i in 1:length(coo)) {
rgl.viewpoint(coo[i])
#filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
#rgl.snapshot(filename)
}

# simple persp3d adaption
f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
x <- seq(-10, 10, length= 30); y <- x
z <- outer(x, y, f)
persp3d(x, y, z)
coo <- 1:360
for(i in 1:length(coo)) {
rgl.viewpoint(coo[i])
#filename <- paste("pic", formatC(i,digits=1,flag="00"), ".png", sep="")
#rgl.snapshot(filename)
}



-
Dr. Carlo Giovanni Camarda
Research Scientist
Max Planck Institute for Demographic Research
Laboratory of Statistical Demography
Konrad-Zuse-Straße 1
18057 Rostock - Germany
Phone: +49 (0)381 2081 172
Fax: +49 (0)381 2081 472
cama...@demogr.mpg.de
http://www.demogr.mpg.de/en/staff/camarda/default.htm
-


--
This mail has been sent through the MPI for Demographic ...{{dropped:10}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Genstat into R - Randomisation test

[Mike Lawrence]

*sigh*

and the elapsed time line should be:

cat('Elapsed time:',proc.time()[1]-start.time)


On Wed, Apr 8, 2009 at 11:35 AM, Mike Lawrence  wrote:
> *Hits head*
> Of course, the approach taken by your Genstat code of only shuffling
> one variable is sufficient and faster:
>
> n.obs = 100
> ID=rnorm(n.obs)
> CD=rnorm(n.obs)
>
> obs.cor = cor(ID,CD)
>
> num.permutations = 1e4
> perm.cor = rep(NA,num.permutations)
> start.time=proc.time()[1]
> index = 1:n.obs
> for(i in 1:num.permutations){
>        IDorder = sample(index)
>        perm.cor[i] = .Internal(cor(ID[IDorder], CD, 4, FALSE))
> }
> cat('Elapsed time:',start.time-proc.time(1))
> sum(perm.cor>obs.cor)/num.permutations
>
>
> On Wed, Apr 8, 2009 at 11:33 AM, Mike Lawrence  wrote:
>> "Assuming you have a data frame with columns ID & CD, this should do it"
>>
>> Oops, the code I sent doesn't assume this. It assumes that you have
>> two vectors, ID & CD, as generated in the first 3 lines.
>>
>> On Wed, Apr 8, 2009 at 11:31 AM, Mike Lawrence  wrote:
>>> Looks like that code implements a non-exhaustive variant of the
>>> randomization test, sometimes called a permutation test. Assuming you
>>> have a data frame with columns ID & CD, this should do it:
>>>
>>> n.obs = 100
>>> ID=rnorm(n.obs)
>>> CD=rnorm(n.obs)
>>>
>>> obs.cor = cor(ID,CD)
>>>
>>> num.permutations = 1e4
>>> perm.cor = rep(NA,num.permutations)
>>> start.time=proc.time()[1]
>>> index = 1:n.obs
>>> for(i in 1:num.permutations){
>>>        IDorder = sample(index)
>>>        CDorder = sample(index)
>>>        perm.cor[i] = .Internal(cor(ID[IDorder], CD[CDorder], 4, FALSE))
>>> }
>>> cat('Elapsed time:',start.time-proc.time(1))
>>> sum(perm.cor>obs.cor)/num.permutations
>>>
>>>
>>> On Wed, Apr 8, 2009 at 10:02 AM, Anne Kempel  wrote:
 Hello everybody,
 I have a question. I would like to get a correlation between constitutive
 and induced plant defence which I messured on 30 plant species. So I have
 table with Species, Induced defence (ID), and constitutive defence (CD).
 Since Induced and constitutive defence are not independant (so called
 spurious correlation) I should do a randomisation test. I have a syntax of
 my supervisor in Genstat, but I would really like to try this in R.


 "data from trade-off.IDCD"
 list
 variate [nval=1000] slope
 calc ID1=ID

 graph ID; CD
 calc b=corr(ID; CD)
 calc slope$[1]=b

 "slope$[1] is the correlation before permutating the data"

 for i=2...1000
    randomize ID1
    calc b=corr(CD1; ID1)
    calc slope$[i]=b
 endfor

 hist slope
 describe slope
 quantile [proportion=!(0.0005,0.005, 0.025, 0.975, 0.995,0.9995)]slope
 print slope$[1]
 corr[p=corr] ID,CD


 DHISTOGRAM [WINDOW=1; ORIENTATION=vertical; KEY=0; BARCOVERING=1.0] slope;
 PEN=2

 Does anybody have done something similar and has any idea how to make the
 randomisation part?
 I would be very grateful for any help!!
 Thanks in advance,
 Anne




 --
 Anne Kempel
 Institute of Plant Sciences
 University of Bern
 Altenbergrain 21
 CH-3103 Bern
 kem...@ips.unibe.ch

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

>>>
>>>
>>>
>>> --
>>> Mike Lawrence
>>> Graduate Student
>>> Department of Psychology
>>> Dalhousie University
>>>
>>> Looking to arrange a meeting? Check my public calendar:
>>> http://tinyurl.com/mikes-public-calendar
>>>
>>> ~ Certainty is folly... I think. ~
>>>
>>
>>
>>
>> --
>> Mike Lawrence
>> Graduate Student
>> Department of Psychology
>> Dalhousie University
>>
>> Looking to arrange a meeting? Check my public calendar:
>> http://tinyurl.com/mikes-public-calendar
>>
>> ~ Certainty is folly... I think. ~
>>
>
>
>
> --
> Mike Lawrence
> Graduate Student
> Department of Psychology
> Dalhousie University
>
> Looking to arrange a meeting? Check my public calendar:
> http://tinyurl.com/mikes-public-calendar
>
> ~ Certainty is folly... I think. ~
>



-- 
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University

Looking to arrange a meeting? Check my public calendar:
http://tinyurl.com/mikes-public-calendar

~ Certainty is folly... I think. ~

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] default print format for large numbers

[Richard . Cotton]

> Numbers like ``1239178547.653775" is inserted into a vector. I print
> the vector:
> 
> > route_9_80_end
>  [1] 1239178522 1239178526 1239178524 1239178524 1239178524 
> 1239178523 1239178524 1239178522 1239178521 1239178565 1239178566 
1239178566
> [13] 1239178565 1239178566 1239178566 1239178565 1239178566 1239178566
> 
> and the decimals are suppressed. While R is wonderful and 
> wonderfully documented, I can only find some quasi relavent 
> information on options(digits) or some such as might apply here. How
> can I force R to print some part of the decimal i.e. back to 
> 1239178547.653775?

When you create the vector, each number is stored as a double-precision 
floating point number.  When you type the name of the variable the print 
function is called, which displays the numbers on the console.

You can also manually call the print function, and set the digits argument 
to a high number to show all digits.

If you want to see all the digits every time any number is displayed, use 
options(digits).

For example:
x <- 1239178547.653775
x
# [1] 1239178548
print(x)
# [1] 1239178548
print(x, digits=20)
# [1] 1239178547.653775
options(digits=20)
x
# [1] 1239178547.653775

Note that the accuracy of the floating point numbers is limited to 
.Machine$double.eps, which is typically 2.2e-16.

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.