Re: [R] Structural equation modeling in R(lavaan,sem)
On 03/29/2011 10:49 PM, jouba wrote: Dear all I have an error mesage « error message : the MLM estimator can not be used when data are incomplete » When i use the function sem(package lavaan) and when i fill the paramters estimator and missing estimator="MLM", missing="ml" i understand by this that i am not allowed to use estimator= ``mlm `` in parralle to the parameter missing You understood correctly. Estimator 'MLM' (=ML estimation + Satorra-Bentler scaled test statistic) is only defined for complete data. The equivalent for incomplete data is the estimator 'MLR' (=ML estimation + Yuan-Bentler scaled test statistic). Yves. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem with Snowball & RWeka
Dear Forum,
I have also problems Snowball (macos, fresh install of rJava, RWeka, RWekajars,
Snowball
from http://cran.ch.r-project.org/bin/macosx/leopard/ ):
> library(Snowball)
> example(SnowballStemmer)
SnwblS> ## Test the supplied vocabulary for the default stemmer ('porter'):
SnwblS> source <- readLines(system.file("words", "porter","voc.txt",
SnwblS+ package = "Snowball"))
SnwblS> result <- SnowballStemmer(source)
Erreur dans .jnew(name) :
java.lang.InternalError: Can't start the AWT because Java was started on the
first thread. Make sure StartOnFirstThread is not specified in your
application's Info.plist or on the command line
Trying to add database driver (JDBC): RmiJdbc.RJDriver - Warning, not in
CLASSPATH?
Trying to add database driver (JDBC): jdbc.idbDriver - Warning, not in
CLASSPATH?
Trying to add database driver (JDBC): org.gjt.mm.mysql.Driver - Warning, not in
CLASSPATH?
Trying to add database driver (JDBC): com.mckoi.JDBCDriver - Warning, not in
CLASSPATH?
Trying to add database driver (JDBC): org.hsqldb.jdbcDriver - Warning, not in
CLASSPATH?
> example(SnowballStemmer)
SnwblS> ## Test the supplied vocabulary for the default stemmer ('porter'):
SnwblS> source <- readLines(system.file("words", "porter","voc.txt",
SnwblS+ package = "Snowball"))
SnwblS> result <- SnowballStemmer(source)
SnwblS> target <- readLines(system.file("words", "porter", "output.txt",
SnwblS+ package = "Snowball"))
SnwblS> ## Any differences?
SnwblS> any(result != target)
[1] TRUE
Stemmer 'porter' unknown!
> sessionInfo()
R version 2.12.1 (2010-12-16)
Platform: i386-apple-darwin9.8.0/i386 (32-bit)
locale:
[1] fr_FR/fr_FR/fr_FR/C/fr_FR/fr_FR
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] Snowball_0.0-7
loaded via a namespace (and not attached):
[1] grid_2.12.1 rJava_0.8-8 RWeka_0.4-3 RWekajars_3.7.3-1
tools_2.12.1
Any hints?
Thanks!
Le 24 mars 2011 à 12:31, Mike Marchywka a écrit :
>
>
>
>> Date: Thu, 24 Mar 2011 03:35:31 -0700
>> From: kont...@alexanderbachmann.de
>> To: r-help@r-project.org
>> Subject: [R] Problem with Snowball & RWeka
>>
>> Dear Forum,
>>
>> when I try to use SnowballStemmer() I get the following error message:
>>
>> "Could not initialize the GenericPropertiesCreator. This exception was
>> produced: java.lang.NullPointerException"
>>
>> It seems to have something to do with either Snowball or RWeka, however I
>> can't figure out, what to do myself. If you could spend 5 minutes of your
>> valuable time, to help me or give me a hint where to look for, it would be
>> very much appreciated.
>>
>> Thank you very much.
>>
> If you only want answers from people who have encountered this exact problem
> before then that's great but you are more likely to get a useful response
> if you include reproducible code and some data to produce the error you have
> seen. Sometimes I investigate these things because they involve a package
> or objective I wanted to look at anyway. It could be that the only problem
> is that the OP
> missed something in documentation or had typo etc.
>
> In this case, to pursue it from the perspective of debuggin the code,
> you probably want to find some way to get a stack trace and then find out
> which java variable was null and relate it back to how you invoked it.
> This likely points to a missing object
> in your call or maybe the installation lacked a dependency as this
> occured during init, but hard to speculate with what you have provided.
> You could try reinstalling and check for errors.
>
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Jean-Pierre Müller
SSP / Anthropole 4123 / UNIL / CH - 1015 Lausanne
Voice:+41 21 692 3116 / Fax:+41 21 692 3115
Please avoid sending me Word or PowerPoint attachments.
See http://www.gnu.org/philosophy/no-word-attachments.html
S'il vous plaît, évitez de m'envoyer des attachements au format Word ou
PowerPoint.
Voir http://www.gnu.org/philosophy/no-word-attachments.fr.html
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] a for loop to lapply
Hi Alaios,
there are some problems with your code:
1. In the expression "Shadowlist[,,i]<-i", i corresponds to dimmaps, not
dimx. The for loop should be either
for (i in c(1:dimmaps)){
Shadowlist[,,i]<-i
}
or
for (i in c(1:dimx)){
Shadowlist[i,,]<-i
}
2. in lapply, the function must be second argument, your expression
lapply(seq(1:dimx),Shadowlist[,,seq(1:dimx)],returni)
fails because lapply wants to call Shadowlist[...] as a function. You
probably meant (including the correction above):
lapply(seq(1:dimmaps),returni,Shadowlist[,,seq(1:dimmaps)])
or
lapply(seq(1:dimx),returni,Shadowlist[seq(1:dimx),,])
3. lapply always returns a list, there is no way to make it return an array.
apply is better to deal with arrays, but as far as I know it cannot return
arrays with dimension greater than 2.
Finally a general remark for future posts: Please make sure your examples are
self-contained and executable (i.e. they don't throw an error and they do not
depend on undefined variables, such as dimx etc. in this case)
Best regards,
Andreas
--
Andreas Borg
Medizinische Informatik
UNIVERSITÄTSMEDIZIN
der Johannes Gutenberg-Universität
Institut für Medizinische Biometrie, Epidemiologie und Informatik
Obere Zahlbacher Straße 69, 55131 Mainz
www.imbei.uni-mainz.de
Telefon +49 (0) 6131 175062
E-Mail: b...@imbei.uni-mainz.de
Diese E-Mail enthält vertrauliche und/oder rechtlich geschützte Informationen.
Wenn Sie nicht der
richtige Adressat sind oder diese E-Mail irrtümlich erhalten haben, informieren
Sie bitte sofort den
Absender und löschen Sie diese Mail. Das unerlaubte Kopieren sowie die
unbefugte Weitergabe
dieser Mail und der darin enthaltenen Informationen ist nicht gestattet.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] a for loop to lapply
Hello Alex.
A few issues:
* you want seq(dimx) instead of seq(1:dimx) (d'oh)
* I think you have problems with your dimensions in the original code as
well: you use i, which runs up to dimx as an indexer for your third
dimension, of size dimmaps. If dimx > dimmaps, you're in for unexpected
results.
* basic idea of the apply-style functions (nicked *apply below):
- first argument = a collection of items to run over. Could be a
list or a vector
- second argument a function, that could take any of the items in
the collection as its first argument
- other arguments: either tuning parameters (like simplify) for
*apply or passed on as more arguments to the function
- each item from the collection is sequentially fed as the first
argument, the extra arguments (always the same) are also passed to *apply.
- normally, the results of each call are collected into a list,
where the names of the list items refers to your original collection. In
more elaborate versions (sapply) and under some circumstances, this list is
transformed into a simpler structure.
* your test case is rather complicated: I don't think there is a way to make
lapply or one of its cousins to return a threedimensional array just like
that. With sapply (and simplify=TRUE, the default), if the result for each
item of your collection has the same length, the result is coerced into a
twodimensional array with one column for each item in your collection.
* on the other hand, for your example, you probably don't want to use *apply
functions nor loops: it can be done with some clever use of seq and rep and
dim, for sure.
All in all, it seems you may need to get your basics up to speed first, then
shift to *apply (and use a simpler example to get started, like: given a
matrix with two columns, create a vector holding the differences and the
sums of the columns - I know this can be done without *apply as well, but
apart from that it is a more attainable exercise).
Good luck to you on that!
HTH,
Nick Sabbe
--
ping: nick.sa...@ugent.be
link: http://biomath.ugent.be
wink: A1.056, Coupure Links 653, 9000 Gent
ring: 09/264.59.36
-- Do Not Disapprove
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Alaios
Sent: woensdag 30 maart 2011 8:31
To: R-help@r-project.org
Subject: [R] a for loop to lapply
Dear all,
I am trying to learn lapply.
I would like, as a test case, to try the lapply alternative for the
Shadowlist<-array(data=NA,dim=c(dimx,dimy,dimmaps))
for (i in c(1:dimx)){
Shadowlist[,,i]<-i
}
---so I wrote the following---
returni <-function(i,ShadowMatrix) {ShadowMatrix<-i}
lapply(seq(1:dimx),Shadowlist[,,seq(1:dimx)],returni)
So far I do not get same results with both ways.
Could you please help me understand what might be wrong?
Regards
Alex
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] a for loop to lapply
Hi Alex,
lapply is not a substitute for for, so it not only does things
differenly, it does a different thing.
> Shadowlist<-array(data=NA,dim=c(dimx,dimy,dimmaps))
> for (i in c(1:dimx)){
> Shadowlist[,,i]<-i
> }
Note that your test case is not reproducible as you haven't defined
dimx, dimy, dimmaps.
> returni <-function(i,ShadowMatrix) {ShadowMatrix<-i}
> lapply(seq(1:dimx),Shadowlist[,,seq(1:dimx)],returni)
> So far I do not get same results with both ways.
> Could you please help me understand what might be wrong?
I don't suppose you're getting any results with lapply this way at all.
1. The second argument to lapply should be a function but here you
have something else.
2. Lapply returns a list of results and does *not* modify its
arguments. If you really want to, you can use <<- within lapply, so
your example (slightly modified) could look something like this:
A<-B<- array(data=NA,dim=c(5,5,5))
for (i in c(1:5))A[,,i]<-i
lapply(1:5, function(i) B[,,i][]<<-i)
all.equal(A,B)
Notice that in this case, what lapply returns is not identical to what
it does to B, and this is not nice. The normal use of lapply is
applying a function to each element of a list (or a vector)...
C<- lapply(1:5, diag) # a list with 5 elements, matrices with increasing size
lapply(C, dim) # dimensions of each matrix
# sapply would give it in a more convenient form
# or you could print it more concisely as str(lapply(C, dim))
the equivalent with for would be:
for(i in 1:5) print(dim(C[[i]]))
The _for_ version is less short and there are much more ['s and ('s to
take care about so lapply/sapply is easier here. Notice also that it
creates a new variable (i) or potentially overwrites anything named i
in your workspace. Lapply doesn't do it unless you explicitly tell it
to (using things like <<- , not recommended).
Good luck,
Kenn
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] a for loop to lapply
alaios wrote:
>
> I am trying to learn lapply.
> I would like, as a test case, to try the lapply alternative for the
> Shadowlist<-array(data=NA,dim=c(dimx,dimy,dimmaps))
> for (i in c(1:dimx)){
> Shadowlist[,,i]<-i
> }
> ---so I wrote the following---
> returni <-function(i,ShadowMatrix) {ShadowMatrix<-i}
> lapply(seq(1:dimx),Shadowlist[,,seq(1:dimx)],returni)
>
>
The basic problem is that you have an array, but you name it "list". This is
valid code, but by calling it so you shot a list-bullet into your
vector-foot. An array should be treated like an array, and the best way to
fill it uses simple vectoring. No loop, no xapply needed.
Use lapply if you have a data.frame or (more generic) a real list, which can
contain rather heterogeneous components. And try to forget what you learned
in your c++ course: Seeing
lapply(,Shadowlist[,,seq(1:dimx)],)
to fill Shadowlist tell me that it must be wrong, because R (in all
recommended cases) never "fills" a parameter passed to it like c(++) can do.
You should always have something like:
Shadowlist <- lapply()
Dieter
#--
dimx=2 # don't forget to make the example self-contained
dimy=3
dimmaps=3
ShadowArray<-array(data=NA,dim=c(dimx,dimy,dimmaps))
# Are you shure you mean "dimx" ? But let's assume it is correct
for (i in c(1:dimx)){
ShadowArray[,,i]<-i
}
ShadowArray
ShadowArray<-array(data=NA,dim=c(dimx,dimy,dimmaps))
ShadowArray[,,1:dimx]<-1:dimx
ShadowArray
--
View this message in context:
http://r.789695.n4.nabble.com/a-for-loop-to-lapply-tp3417169p3417377.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] new syntax: bash-like pipe operator
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1
On 29/03/11 18:00, Gabor Grothendieck wrote:
> On Tue, Mar 29, 2011 at 9:56 AM, Robert Sugar wrote:
>> Dear R Community,
>>
>> One thing that always bugged me about R is the abundance of multi-level
>> nested statements, like:
>>
>> cumsum(ifelse(c(1,diff(v)),1,0))
>>
>> because:
>>
>> a) you have to read them inside out as opposed to left-to-right
>> b) in the console you always have to go back and type the parenthesis if you
>> want to nest your statement
>> I kind of like the UNIX pipe operator as a very good abstraction of a data
>> flow and would fancy to have something like this in R:
>> v %|% diff %|% c(1, value) %|% ifelse(value, 1, 0) %|% cumsum
>> so I went ahead and wrote it:
>>
>> "%|%" <- function(x,y)
>> {
>> # Operator similar to the UNIX pipe operator. Allows to set up a chain of
>> functions
>> # where the result of one function will be the input of the next.
>> #
>> # x: any R statement. will be passed to y as a parameter
>> # y: either a) the name of a single-parameter function or b) a function
>> call, where "value"
>> # will be replaced with the x parameter
>> #
>> # Example use 1: > c(1,2,3) %|% sum
>> # results in [1] 6
>> # Example use 2: > c(1,2,3) %|% rep(value, 10)
>> # results in [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2
>> 3
>>
>>thecall<-match.call() #get call
>>
>>#for functions with one parameter call it with x. Example: "c(1,2,3)
>> %|% sum" -> 6
>>if(is.name(thecall$y) || is.function(thecall$y)) return (y(x))
>>
>>#instead of value use x
>>value <- eval(thecall$x)
>>
>>return(eval(thecall$y)) #evaluate
>> }
>>
>> would be happy to receive some feedback on
>>
>> 1. Would you find this useful?
Yes - I definitaly like it - although I will in programming cases still
use the nested calls, but in an "ad-hoc" analysis situation, this would
be really usefull.
>> 2. Is there any bugs/features that you'd like to mention?
While we are at it:
separate pipes for output and error stream which can be send to
different targets (displs, file, console, ...). But I guess that would
require some major changes.
>> 3. Should I go ahead and make an R package on that?
Yes please - I would use it and give feedback.
>
> That is interesting although note that a simple way to get that order
> without introducing new features is to use intermediate variables.
> Here we use the variable . in order to save space and use .= to assign
> to it.
>
> # Example 1 above
> .= c(1, 2, 3); sum(.)
>
> # Example 2 above
> .= c(1, 2, 3); rep(., 10)
>
> # example in preamble
> v <- c(1, 2, 3)
> .= diff(v); .= c(1, .); .= ifelse(., 1, 0); cumsum(.)
Definitely true, but not nearly as nice and elegant (for Linux user
eyes) as the | piping operator.
Looking forward to the package,
Rainer
>
>
- --
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
Biology, UCT), Dipl. Phys. (Germany)
Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa
Tel:+33 - (0)9 53 10 27 44
Cell: +27 - (0)8 39 47 90 42
Fax (SA): +27 - (0)8 65 16 27 82
Fax (D) : +49 - (0)3 21 21 25 22 44
Fax (FR): +33 - (0)9 58 10 27 44
email: rai...@krugs.de
Skype: RMkrug
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
iEYEARECAAYFAk2S5k4ACgkQoYgNqgF2egqaJQCghjFqor1sDCdwwwMXIoLEzUhZ
EZAAnA4Pt0uTCKURUWvclQ2aXXARFK2Y
=+2uT
-END PGP 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] Using graphics straight from R into published articles
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 On 29/03/11 19:52, Philipp Pagel wrote: > On Tue, Mar 29, 2011 at 09:31:18AM -0700, blanco wrote: >> I was just wondering if people use graphics from R straight into articles or >> are they always edited in some way; fonts, headers, axis, color etc? Using >> photoshop or some other programs? >> >> I would like to think it is possible, better and more profession to do it >> all in R. Absolutely. You can get one paper = one file (if you want to, including your data - but usually that one is separate) when you use e.g. sweave or org-mode babel (http://orgmode.org/, http://orgmode.org/worg/org-contrib/babel/) - all graphs are generated when you compile your document. I use org-mode. >> I tried google and the search option but found nothing on the topic. >> >> What are the experiences for all the professionals out there that use R? >> Are there any articles on this specific subject? > > I'm not aware of any articles on the topic but I can share what I do: > > 95% of the time I tweak various graphics parameters in R and see no > necessity for postprocessing in other applications. Up to now, even 100% > > In about 5% I do some manual editing for a "camera ready" figure. > These are usually the result of exotic request from referees. But > under no circumstances would I use Photoshop or any other pixel > graphics software for this. My R graphics are always created as eps or > pdf vector graphics and any editing is done with a proper vector > graphics software (Illustrator or Inkscape). Absolutely vector - no jpeg, png, ... although it takes sometimes some convincing of co-authors... > > I share your feeling that it is better to do as much as possible in R > because it means that I won't have to do it again if I need to produce > another revision of the figure - all it takes is anoother run of my > script. And I can re-use good solutions in the future. Any manual > touch-ups have to be done manually every single time => not my idea of > efficiency. Absolutely - reproducible research. And if you have several graphs, they all share the same layout, font, fontsize, margins, ... I rather spend an hour fiddling on my graphs in R (and sometimes swearing about the layout of multiple graphs in one figure) then spend half an hour "painting" my graphs. Cheers, Rainer > > cu > Philipp > - -- Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, UCT), Dipl. Phys. (Germany) Centre of Excellence for Invasion Biology Natural Sciences Building Office Suite 2039 Stellenbosch University Main Campus, Merriman Avenue Stellenbosch South Africa Tel:+33 - (0)9 53 10 27 44 Cell: +27 - (0)8 39 47 90 42 Fax (SA): +27 - (0)8 65 16 27 82 Fax (D) : +49 - (0)3 21 21 25 22 44 Fax (FR): +33 - (0)9 58 10 27 44 email: rai...@krugs.de Skype: RMkrug -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.10 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/ iEYEARECAAYFAk2S6CwACgkQoYgNqgF2egpvFgCfewRlqSfz+RSGmxcn15I4H9Yh riUAn2KzUCS8wD55HFQy+M6O5QRf0yIr =v8nN -END PGP 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] calculating the mode in R...
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 On 30/03/11 02:47, Fernando Marmolejo Ramos wrote: > Dear R users > > I?m aware that the package ?modest? is useful to find the mode in an array. > > However, I?d like to know if someone has translated the ?mode? function > built-in > in MATLAB into R language. This function finds the most frequent value in an > array (http://www.mathworks.com/help/techdoc/ref/mode.html). This sounds like a combination of the table() which tabulates your data, and the max() function to identify the max - but then you have to specify what to do when two values have the same count. Rainer > > Best > > Fer > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. - -- Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, UCT), Dipl. Phys. (Germany) Centre of Excellence for Invasion Biology Natural Sciences Building Office Suite 2039 Stellenbosch University Main Campus, Merriman Avenue Stellenbosch South Africa Tel:+33 - (0)9 53 10 27 44 Cell: +27 - (0)8 39 47 90 42 Fax (SA): +27 - (0)8 65 16 27 82 Fax (D) : +49 - (0)3 21 21 25 22 44 Fax (FR): +33 - (0)9 58 10 27 44 email: rai...@krugs.de Skype: RMkrug -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.10 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/ iEYEARECAAYFAk2S6TIACgkQoYgNqgF2ego48ACdHfCL+BdGA6wZ4bHrjq2wCXJW vIoAnREiHFeSbJy9vYQPEnRhpV6nLbNB =suk2 -END PGP 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: wireframe "eats up" points; how to make points on wireframe visible?
On Fri, Mar 4, 2011 at 1:47 PM, Marius Hofert wrote:
> Dear expeRts,
>
> I would like to add two points to a wireframe plot. The points have (x,y,z)
> coordinates
> where z is determined to be on the wireframe [same z-value]. Now something
> strange
> happens. One point is perfectly plotted, the other isn't shown at all. It only
> appears if I move it upwards in z-direction by adding a positive number. So
> somehow
> it disappears in the wireframe-surface *although* the plot symbol [the cross]
> has
> a positive length in each dimension [I also chose cex=5 to make it large
> enough so
> that it should (theoretically) be visible].
>
> My wireframe plot is a complicated function which I cannot post here. Below
> is a minimal
> example, however, it didn't show the same problem [the surface is too nice I
> guess].
> I therefore *artifically* create the problem in the example below so that you
> know
> what I mean. For one of the points, I subtract an epsilon [=0.25] in
> z-direction and
> suddenly the point completely disappears. The strange thing is that the point
> is
> not even "under" the surface [use the screen-argument to rotate the wireframe
> plot to check this],
> it's simply gone, eaten up by the surface.
>
> How can I make the two points visible?
> I also tried to use the alpha-argument to make the wireframe transparent, but
> I couldn't
> solve the problem.
>
> Cheers,
>
> Marius
>
> PS: One also faces this problem for example if one wants to make points
> visible that are on "opposite sides" of the wireframe.
>
> library(lattice)
>
> f <- function(x) 1/((1-x[1])*(1-x[2])+1)
>
> u <- seq(0, 1, length.out=20)
> grid <- expand.grid(x=u, y=u)
> x <- grid[,1]
> y <- grid[,2]
> z <- apply(grid, 1, f)
>
> pt.x <- c(0.2, 0.5)
> pt.y <- c(0.6, 0.8)
> eps <- 0.25
> pts <- rbind(c(pt.x, f(pt.x)-eps), c(pt.y, f(pt.y))) # points to add to the
> wireframe
The reason in this case is fairly obvious: you have
> pts
[,1] [,2] [,3]
[1,] 0.2 0.5 0.4642857
[2,] 0.6 0.8 0.9259259
So the z-value for Point 1 is 0.4642857, which is less than 0.5, the
minimum of the z-axis.panel.3dscatter() "clips" any points outside the
range of the bounding box, so this point is not plotted.
I can't say what the problem was in your original example without
looking at it, but I would guess it's caused by something similar.
Hope this is not too late to be useful.
-Deepayan
>
> wireframe(z~x*y, pts=pts, aspect=1, scales=list(col=1, arrows=FALSE),
> panel.3d.wireframe = function(x,y,z,xlim,ylim,zlim,xlim.scaled,
> ylim.scaled,zlim.scaled,pts,...){
> panel.3dwire(x=x, y=y, z=z, xlim=xlim, ylim=ylim, zlim=zlim,
> xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled,
> zlim.scaled=zlim.scaled, ...)
> panel.3dscatter(x=pts[,1], y=pts[,2], z=pts[,3],
> xlim=xlim, ylim=ylim, zlim=zlim,
> xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled,
> zlim.scaled=zlim.scaled, type="p", col=c(2,3),
> cex=1.8, .scale=TRUE, ...)
> }, key=list(x=0.5, y=0.95, points=list(col=c(2,3)),
> text=list(c("Point 1", "Point 2")),
> cex=1, align=TRUE, transparent=TRUE))
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Using graphics straight from R into published articles
Large snip. > Absolutely vector - no jpeg, png, ... although it takes > sometimes some convincing of co-authors... > That depends on the kind of graph. I aggree that you should try vector at first. But when it generates very larges files (e.g. scatterplots with thousands of points) then you better switch to bitmaps like tiff or png. Jpeg can create artefacts, so is not very good for graphics. Best regards, Thierry __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] calculating the mode in R...
Hello Rainer and Fernando,
Actually, I think this function should involve the which.max (not max):
Here is a tiny function to perform this (with smarter handeling of multiple
modes and giving proper warning in such cases)
The.mode <- function(x, show_all_modes = F)
{
x_freq <- table(x)
mode_location <- which.max(x_freq)
The_mode <- names(x_freq)[mode_location]
Number_of_modes <- length(mode_location)
#
if(show_all_modes) {
if(Number_of_modes >1) {
warning(paste("Multiple modes exist - returning all",Number_of_modes,"of
them"))}
return(The_mode)
} else {
if(Number_of_modes >1) {
warning(paste("Multiple modes exist - returning only the first one out
of", Number_of_modes))}
return(The_mode[1])
}
}
Cheers,
Tal
Contact
Details:---
Contact me: tal.gal...@gmail.com | 972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
--
On Wed, Mar 30, 2011 at 11:26 AM, Rainer M Krug wrote:
> -BEGIN PGP SIGNED MESSAGE-
> Hash: SHA1
>
> On 30/03/11 02:47, Fernando Marmolejo Ramos wrote:
> > Dear R users
> >
> > I?m aware that the package ?modest? is useful to find the mode in an
> array.
> >
> > However, I?d like to know if someone has translated the ?mode? function
> built-in
> > in MATLAB into R language. This function finds the most frequent value in
> an
> > array (http://www.mathworks.com/help/techdoc/ref/mode.html).
>
> This sounds like a combination of the table() which tabulates your data,
> and the max() function to identify the max - but then you have to
> specify what to do when two values have the same count.
>
> Rainer
>
> >
> > Best
> >
> > Fer
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
>
> - --
> Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
> Biology, UCT), Dipl. Phys. (Germany)
>
> Centre of Excellence for Invasion Biology
> Natural Sciences Building
> Office Suite 2039
> Stellenbosch University
> Main Campus, Merriman Avenue
> Stellenbosch
> South Africa
>
> Tel:+33 - (0)9 53 10 27 44
> Cell: +27 - (0)8 39 47 90 42
> Fax (SA): +27 - (0)8 65 16 27 82
> Fax (D) : +49 - (0)3 21 21 25 22 44
> Fax (FR): +33 - (0)9 58 10 27 44
> email: rai...@krugs.de
>
> Skype: RMkrug
> -BEGIN PGP SIGNATURE-
> Version: GnuPG v1.4.10 (GNU/Linux)
> Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
>
> iEYEARECAAYFAk2S6TIACgkQoYgNqgF2ego48ACdHfCL+BdGA6wZ4bHrjq2wCXJW
> vIoAnREiHFeSbJy9vYQPEnRhpV6nLbNB
> =suk2
> -END PGP 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.
>
[[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] Not all rows are being read-in
On Tue, Mar 29, 2011 at 06:58:59PM -0400, Dimitri Liakhovitski wrote: > I have a tab-delimited .txt file (size 800MB) with about 3.4 million > rows and 41 columns. About 15 columns contain strings. > Tried to read it in in R 2.12.2 on a laptop that has Windows XP: > mydata<-read.delim(file="FileName.TXT",sep="\t") > R did not complain (!) and I got: dim(mydata) 1692063 41. My guess would be that there are (unexpected) quotes and/or double quotes in your file and so R thinks that rather large blocks of your file are actually very long strings. This routinely happens in situations like this: ID x description 1 0.4 my first measurement 2 1.6 Normal 5" object 3 0.4 Some measuremetn 4 0.7 A 4" long sample R thinks that the description in row 2 ends in row 4 and you loose data. Try read.delim(..., quote=""). cu Philipp -- Dr. Philipp Pagel Lehrstuhl für Genomorientierte Bioinformatik Technische Universität München Wissenschaftszentrum Weihenstephan Maximus-von-Imhof-Forum 3 85354 Freising, Germany http://webclu.bio.wzw.tum.de/~pagel/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Class noquote to matrix
On 2011-03-29 19:12, Mark Ebbert wrote:
Hi,
I apologize if the solution is right in front of me, but I can't find anything on
how to convert a class of 'noquote' to 'matrix'. I've tried as.matrix and I've
tried coercing it by 'class(x)<-matrix', but of course that didn't work. I've
been using the function 'symnum' which returns an object of class 'noquote'.
Here is an example from the 'symnum' help page:
pval<- rev(sort(c(outer(1:6, 10^-(1:3)
symp<- symnum(pval, corr=FALSE,
cutpoints = c(0, .001,.01,.05, .1, 1),
symbols = c("***","**","*","."," "))
One way would be to unclass the object first:
as.matrix(unclass(symp))
Peter Ehlers
Thanks!
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Using graphics straight from R into published articles
On Wed, Mar 30, 2011 at 08:48:55AM +, ONKELINX, Thierry wrote: > Large snip. > > > Absolutely vector - no jpeg, png, ... although it takes > > That depends on the kind of graph. I aggree that you should try > vector at first. But when it generates very larges files (e.g. > scatterplots with thousands of points) then you better switch to > bitmaps like tiff or png. Jpeg can create artefacts, so is not very > good for graphics. True. Sometimes one can get away with switching from a normal scatterplot to hexbin or something like this but if that is not anoption a high resolution tiff or png is the way out. And of course, I agree that jpeg should never be used for graphs. cu Philipp -- Dr. Philipp Pagel Lehrstuhl für Genomorientierte Bioinformatik Technische Universität München Wissenschaftszentrum Weihenstephan Maximus-von-Imhof-Forum 3 85354 Freising, Germany http://webclu.bio.wzw.tum.de/~pagel/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using graphics straight from R into published articles
On Wed, 30 Mar 2011, Philipp Pagel wrote: On Wed, Mar 30, 2011 at 08:48:55AM +, ONKELINX, Thierry wrote: Large snip. Absolutely vector - no jpeg, png, ... although it takes That depends on the kind of graph. I aggree that you should try vector at first. But when it generates very larges files (e.g. scatterplots with thousands of points) then you better switch to bitmaps like tiff or png. Jpeg can create artefacts, so is not very good for graphics. True. Sometimes one can get away with switching from a normal scatterplot to hexbin or something like this but if that is not anoption a high resolution tiff or png is the way out. And of course, I agree that jpeg should never be used for graphs. That depends on the definition of 'graphs'. One reason that I added jpeg support to Sweave is that some of us do work in areas where figures (and it is 'figures' you put in papers in my native language) are images, even photographic images. cu Philipp -- Dr. Philipp Pagel Lehrstuhl für Genomorientierte Bioinformatik Technische Universität München Wissenschaftszentrum Weihenstephan Maximus-von-Imhof-Forum 3 85354 Freising, Germany http://webclu.bio.wzw.tum.de/~pagel/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Brian D. Ripley, rip...@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595__ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using graphics straight from R into published articles
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 On 30/03/11 11:12, Philipp Pagel wrote: > On Wed, Mar 30, 2011 at 08:48:55AM +, ONKELINX, Thierry wrote: >> Large snip. >> >>> Absolutely vector - no jpeg, png, ... although it takes >> >> That depends on the kind of graph. I aggree that you should try >> vector at first. But when it generates very larges files (e.g. >> scatterplots with thousands of points) then you better switch to >> bitmaps like tiff or png. Jpeg can create artefacts, so is not very >> good for graphics. > > True. Sometimes one can get away with switching from a normal > scatterplot to hexbin or something like this but if that is not > anoption a high resolution tiff or png is the way out. Well - is there a rule without an exception? Of course - you are right. Cheers, Rainer > > And of course, I agree that jpeg should never be used for graphs. > > cu > Philipp > - -- Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, UCT), Dipl. Phys. (Germany) Centre of Excellence for Invasion Biology Natural Sciences Building Office Suite 2039 Stellenbosch University Main Campus, Merriman Avenue Stellenbosch South Africa Tel:+33 - (0)9 53 10 27 44 Cell: +27 - (0)8 39 47 90 42 Fax (SA): +27 - (0)8 65 16 27 82 Fax (D) : +49 - (0)3 21 21 25 22 44 Fax (FR): +33 - (0)9 58 10 27 44 email: rai...@krugs.de Skype: RMkrug -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.10 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/ iEUEARECAAYFAk2S9rkACgkQoYgNqgF2egpp2ACcCrToeRM2OWPBVAs99Q/u0c/z Vl8AmKoP0l2RR+EyL+VVlYM2S+hQdhw= =hKPu -END PGP 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] How to put line linking two plots
Hello!
Suppose I have three charts like below. The top chart is a general
overview and the bottom charts are related so some point of this chart.
To make clear this relationship I want to draw a line between (4,0.9) in
the top chart and (10,1) in the bottom-left one.
Currently I add it manually using Inkscape on the resulting pdf file.
Is it possible to add it inside R? Should I switch to other charting
packages?
Thanks for the advice!
mario
set.seed(123)
pdf("test.pdf", width=14, height=7)
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(runif(10), type='b')
plot(runif(20), type='l')
plot(runif(20), type='l')
dev.off()
R 2.12.2 on Windows 7 (32bits)
--
Ing. Mario Valle
Data Analysis and Visualization Group| http://www.cscs.ch/~mvalle
Swiss National Supercomputing Centre (CSCS) | Tel: +41 (91) 610.82.60
v. Cantonale Galleria 2, 6928 Manno, Switzerland | Fax: +41 (91) 610.82.82
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] save ordinary numerical calculations as pdf
I'd like to save some calculation outputs as a pdf, to incorporate with others
in a document. I've tried
pdf("filename")
name_of_object_to_output
dev.off()
but it doesn't seem to work, appears that this pdf function is for graphics?
Is there a way to output numerical objects to pdf?
Thanks
J
Dr. Jim Maas
University of East Anglia
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Creating 3 vectors that sum to 1
On Tue, Mar 29, 2011 at 01:42:18PM -0600, Greg Snow wrote:
> Or we could expand a bit more:
>
> require(TeachingDemos)
> require(gtools)
>
> n <- 1000
> rtrg <- matrix(NA, n, 3)
> for (i in 1:n) rtrg[i,] <- diff(c(0, sort(runif(2)), 1))
>
> rtrg2 <- matrix(NA, n, 3)
> for (i in 1:n) {
> tmp <- runif(3)
> rtrg2[i, ] <- tmp/sum(tmp)
> }
>
> rtrg3 <- matrix( rexp(n*3), ncol=3 )
> rtrg3 <- rtrg3/rowSums(rtrg3)
>
> rtrg4 <- rdirichlet(n, rep(1,3))
If i understand correctly, this is an efficient way to generate the
uniform distribution over the triangle. Thank you for pointing this
out.
Generating the uniform distribution is a natural question related
to the original request.
> par(mfrow=c(2,2))
> triplot(rtrg, pch='.') # Looks more uniformly distributed
This distribution is also exactly uniform.
If x is generated as sort(runif(2)), then it is uniformly distributed
over the two dimensional triangle 0 <= x[1] <= x[2] <= 1.
The transformation, which maps x to c(x[1], x[2] - x[1], 1 - x[2])
is linear, so it preserves the uniform distribution and its output
ranges over the triangle with corners [1, 0, 0], [0, 1, 0], [0, 0, 1].
> triplot(rtrg2, col=2, pch='.') # Corners are sparsely populated
The ratio of the density in the center and in the corners seems
to be 27.
Consider a small area A on the triangle between the points (1, 0, 0),
(0, 1, 0), (0, 0, 1). The points x generated as runif(3), which are mapped
to A by the transformation x/sum(x), consist of two types of points.
(1) The cone (or a pyramid) between A and the point (0, 0, 0).
(2) Points on the extension of this cone away from (0, 0, 0),
which are still inside the cube [0, 1]^3.
The volume of (1) depends on the area of A, but not on its location
within the triangle, since the height of the cone is always the same.
If A is close to a corner, then (2) is small compared to (1) (the ratio
converges to 0, if A is limited to small neighborhood of the corner).
If A is close to the center of the triangle, then the union of (1) and (2)
contains a cone approximately 3 times larger than the cone (1). So,
the volume of (1) and (2) together is about 27 times larger than (1) alone.
(The ratio converges to 27, if A is limited to small neighborhood of the
center).
Petr Savicky.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] offset in gam
Hi, we are modelling dusky grouper distribution using underwater visual census. Our response variable is the abundance of groupers but we want to see differences in habitat related to maturity because dusky groupers use different habitats when they're juveniles versus adults. So we're using the % of maturity. ie 1 individual that is 50cm TL gets aprox 0.5 (50% chances to be mature and 50% to be inmature) and so on. If we have let's say another grouper that is 40 cm TL it will have a probability for example 0.7 to be inmature so for our transect we will have a total of 0.7 +0.5= 1.2 for 2 individuals. We're using the Poisson distribution because our data is not overdispersed and it's not normal. But we don't know what offset we should include. Would it be the area of the transect? all transects have the same area, 90 m2. Or should it be the total abundance of groupers per transect? in the example offset=2 Thanks a lot in adavance Lucia -- View this message in context: http://r.789695.n4.nabble.com/offset-in-gam-tp3417491p3417491.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] calculating the mode in R...
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1
On 30/03/11 10:59, Tal Galili wrote:
> Hello Rainer and Fernando,
>
> Actually, I think this function should involve the which.max (not max):
Right - was not awarw of which.max - I would have used which(max(...)) -
which.max is much more elegant.
Thanks,
Rainer
>
> Here is a tiny function to perform this (with smarter handeling of
> multiple modes and giving proper warning in such cases)
> The.mode <- function(x, show_all_modes = F)
It is advisable to use FaALSE instead of F, especially in a function:
you can assign values to F and T, but you can not assign values to FALSE
and TRUE
Cheers,
Rainer
> {
> x_freq <- table(x)
> mode_location <- which.max(x_freq)
> The_mode <- names(x_freq)[mode_location]
> Number_of_modes <- length(mode_location)
> #
> if(show_all_modes) {
> if(Number_of_modes >1) {
>warning(paste("Multiple modes exist - returning
> all",Number_of_modes,"of them"))}
> return(The_mode)
> } else {
> if(Number_of_modes >1) {
>warning(paste("Multiple modes exist - returning only the first one
> out of", Number_of_modes))}
> return(The_mode[1])
> }
> }
>
>
> Cheers,
> Tal
>
> Contact
> Details:---
> Contact me: tal.gal...@gmail.com |
> 972-52-7275845
> Read me: www.talgalili.com (Hebrew) |
> www.biostatistics.co.il (Hebrew) |
> www.r-statistics.com (English)
> --
>
>
>
>
> On Wed, Mar 30, 2011 at 11:26 AM, Rainer M Krug > wrote:
>
> On 30/03/11 02:47, Fernando Marmolejo Ramos wrote:
>> Dear R users
>
>> I?m aware that the package ?modest? is useful to find the mode in
> an array.
>
>> However, I?d like to know if someone has translated the ?mode?
> function built-in
>> in MATLAB into R language. This function finds the most frequent
> value in an
>> array (http://www.mathworks.com/help/techdoc/ref/mode.html).
>
> This sounds like a combination of the table() which tabulates your data,
> and the max() function to identify the max - but then you have to
> specify what to do when two values have the same count.
>
> Rainer
>
>
>> Best
>
>> Fer
>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/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.
- --
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
Biology, UCT), Dipl. Phys. (Germany)
Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa
Tel:+33 - (0)9 53 10 27 44
Cell: +27 - (0)8 39 47 90 42
Fax (SA): +27 - (0)8 65 16 27 82
Fax (D) : +49 - (0)3 21 21 25 22 44
Fax (FR): +33 - (0)9 58 10 27 44
email: rai...@krugs.de
Skype: RMkrug
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
iEYEARECAAYFAk2S9ccACgkQoYgNqgF2egqLtQCeM5ntWTHOy2wyk4CbOxwnUfJG
VZUAnRcOpRImBm/AbhsK14O54+pzCP6H
=QW4d
-END PGP 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] Violin Plot - using ggplot2
Hi there,
I am trying to create a violin plot with ggplot2 (which looks awesome
by the way). I am probably not understanding the code correctly by
this is what I am trying (this works if you want to copy and paste it
in R)
library(ggplot2)
p <- rep(c(rep("condition_a", 4), rep("condition_b", 4)), 2)
q <- c(rep("grp_1", 8), rep("grp_2", 8))
r <- sample(1:5, 16, rep = T)
d <- data.frame(p, q, r)
gg1 <- ggplot(d, aes(y=r, x=factor(p)))
gg2 <- gg1 + geom_ribbon(aes(ymax = ..density.., ymin = -..density..),
stat = "density") + facet_grid(q ~ ., as.table = F, scales = "free_y")
+ labs(x = "some x desc", y = "some y desc")
print(gg2)
The result I am getting is posted as an image at
http://awesomescreenshot.com/0e8ae8c52 (I tried to draw what I would
like to achieve.)
What I would like to achieve instead is the equivalent of this using ggplot2
library(lattice)
bwplot( r ~ p | q,
data=d,
panel = panel.violin
)
What I want is to draw traditional violin plots grouping them by
grp_1, grp_2 etc (so q above). For each group I have two conditions
(p above) which I want to plot violin plots for their values. I want
the plots to be vertical (as drawn in red in my image). Any ideas on
how to change the above ggplot2 code to do just that? I have been
playing with facet_grid all day... Also the values on the axis are
incorrect ... they should be in the range 1:5 (as r above) and a plot
should show the probability density of that value...
Any help? I am going mad...
Many Thanks
JP
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] glm: modelling zeros as binary and non-zeroes as coming from a continuous distribution
Hello, I'd like to implement a regression model for extremely zero-inflated continuous data using a conditional approach, whereby zeroes are modelled as coming from a binary distribution, while non-zero values are modelled as log-normal. So far, I've come across two solutions for this: one, in R, is described in the book by Gelman & Hill (http://www.amazon.com/dp/052168689X), where they just model zeros and non-zeros separately and then bring them together by simulation. I can do this, but it makes it difficult to assess the significance of regression coefficients wrt to zero and each other. Another solution I have been pointed at is in SAS: http://listserv.uga.edu/cgi-bin/wa?A2=ind0805A&L=sas-l&P=R20779, where they use NLMIXED (with only fixed effects) to specify their own log-likelihood function. I'm wondering if there's any way to do the same in R (lme can't deal with this, as far as I'm aware). Finally, I'm wondering whether anyone has experience with the COZIGAM package - does it do something like this? Many thanks, Mikhail __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Loading a package without having to install it?
Dear UseRs, If I build an R package, is there a way to load it without having to install it? Kind regards, -- David Lindelöf, Ph.D. +41 (0)79 415 66 41 or skype:david.lindelof Better Software for Tomorrow's Cities: http://computersandbuildings.com Follow me on Twitter: http://twitter.com/dlindelof __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Violin Plot - using ggplot2
Dear JP,
Please do not cross post between lists.
Here is a possible solution.
library(ggplot2)
p <- rep(c(rep("condition_a", 4), rep("condition_b", 4)), 2)
q <- c(rep("grp_1", 8), rep("grp_2", 8))
r <- sample(1:5, 16, rep = T)
d <- data.frame(p, q, r)
ggplot(d, aes(x = r)) + geom_ribbon(aes(ymax = ..density.., ymin =
-..density..), stat = "density") + facet_grid(q ~ p) + coord_flip()
ggplot(d, aes(x = r, fill = p)) + geom_ribbon(alpha = 0.2, aes(ymax =
..density.., ymin = -..density..), stat = "density") + facet_wrap( ~ q) +
coord_flip()
Best regards,
Thierry
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium
Research Institute for Nature and Forest
team Biometrics & 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 JP
> Verzonden: woensdag 30 maart 2011 11:26
> Aan: r-help@r-project.org
> Onderwerp: [R] Violin Plot - using ggplot2
>
> Hi there,
>
> I am trying to create a violin plot with ggplot2 (which looks awesome
> by the way). I am probably not understanding the code correctly by
> this is what I am trying (this works if you want to copy and paste it
> in R)
>
> library(ggplot2)
> p <- rep(c(rep("condition_a", 4), rep("condition_b", 4)), 2)
> q <- c(rep("grp_1", 8), rep("grp_2", 8))
> r <- sample(1:5, 16, rep = T)
> d <- data.frame(p, q, r)
> gg1 <- ggplot(d, aes(y=r, x=factor(p)))
> gg2 <- gg1 + geom_ribbon(aes(ymax = ..density.., ymin = -..density..),
> stat = "density") + facet_grid(q ~ ., as.table = F, scales = "free_y")
> + labs(x = "some x desc", y = "some y desc")
> print(gg2)
>
> The result I am getting is posted as an image at
> http://awesomescreenshot.com/0e8ae8c52 (I tried to draw what I would
> like to achieve.)
>
> What I would like to achieve instead is the equivalent of
> this using ggplot2
>
> library(lattice)
> bwplot( r ~ p | q,
>data=d,
>panel = panel.violin
> )
>
> What I want is to draw traditional violin plots grouping them by
> grp_1, grp_2 etc (so q above). For each group I have two conditions
> (p above) which I want to plot violin plots for their values. I want
> the plots to be vertical (as drawn in red in my image). Any ideas on
> how to change the above ggplot2 code to do just that? I have been
> playing with facet_grid all day... Also the values on the axis are
> incorrect ... they should be in the range 1:5 (as r above) and a plot
> should show the probability density of that value...
>
> Any help? I am going mad...
>
> Many Thanks
> JP
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Connect observed values by a smooth curve, is there any method to get coordinate value of arbitrary point on the curve?
On 11-03-29 9:32 PM, Zhipeng Wang wrote: Dear all, I have the problem as the title shows. It is time series of measurements. I want to estimate the value at time point when we have no actual measured value (Only in the time range, not for prediction in future time). Hope you could give me some suggestion or hint. Thank you so much. If you want to interpolate the data, use splinefun. For example, x <- 1:5 y <- rnorm(5) f <- splinefun(x,y) curve(f, from=1, to=5) points(x,y) If you want to draw a smooth curve through the data, then use one of the model fitting functions, and then predict to compute values at new x values. For example, library(mgcv) x <- 1:50 # Smoothing needs more data y <- rnorm(50) fit <- gam(y ~ s(x)) newx <- seq(1,50, len=1000) newy <- predict(fit, newdata=data.frame(x=newx)) plot(x,y) lines(newx, newy) 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.
Re: [R] Violin Plot - using ggplot2
Sorry for the cross post - but I will post your answer in the other place
too for completeness
And thanks for your answer... that is brilliant especially the second
transparent plot
Fantastic
On 30 March 2011 12:21, ONKELINX, Thierry wrote:
> Dear JP,
>
> Please do not cross post between lists.
>
> Here is a possible solution.
>
> library(ggplot2)
> p <- rep(c(rep("condition_a", 4), rep("condition_b", 4)), 2)
> q <- c(rep("grp_1", 8), rep("grp_2", 8))
> r <- sample(1:5, 16, rep = T)
> d <- data.frame(p, q, r)
> ggplot(d, aes(x = r)) + geom_ribbon(aes(ymax = ..density.., ymin =
> -..density..), stat = "density") + facet_grid(q ~ p) + coord_flip()
> ggplot(d, aes(x = r, fill = p)) + geom_ribbon(alpha = 0.2, aes(ymax =
> ..density.., ymin = -..density..), stat = "density") + facet_wrap( ~ q) +
> coord_flip()
>
> Best regards,
>
> Thierry
>
>
>
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek
> team Biometrie & Kwaliteitszorg
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium
>
> Research Institute for Nature and Forest
> team Biometrics & 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 JP
> > Verzonden: woensdag 30 maart 2011 11:26
> > Aan: r-help@r-project.org
> > Onderwerp: [R] Violin Plot - using ggplot2
> >
> > Hi there,
> >
> > I am trying to create a violin plot with ggplot2 (which looks awesome
> > by the way). I am probably not understanding the code correctly by
> > this is what I am trying (this works if you want to copy and paste it
> > in R)
> >
> > library(ggplot2)
> > p <- rep(c(rep("condition_a", 4), rep("condition_b", 4)), 2)
> > q <- c(rep("grp_1", 8), rep("grp_2", 8))
> > r <- sample(1:5, 16, rep = T)
> > d <- data.frame(p, q, r)
> > gg1 <- ggplot(d, aes(y=r, x=factor(p)))
> > gg2 <- gg1 + geom_ribbon(aes(ymax = ..density.., ymin = -..density..),
> > stat = "density") + facet_grid(q ~ ., as.table = F, scales = "free_y")
> > + labs(x = "some x desc", y = "some y desc")
> > print(gg2)
> >
> > The result I am getting is posted as an image at
> > http://awesomescreenshot.com/0e8ae8c52 (I tried to draw what I would
> > like to achieve.)
> >
> > What I would like to achieve instead is the equivalent of
> > this using ggplot2
> >
> > library(lattice)
> > bwplot( r ~ p | q,
> >data=d,
> >panel = panel.violin
> > )
> >
> > What I want is to draw traditional violin plots grouping them by
> > grp_1, grp_2 etc (so q above). For each group I have two conditions
> > (p above) which I want to plot violin plots for their values. I want
> > the plots to be vertical (as drawn in red in my image). Any ideas on
> > how to change the above ggplot2 code to do just that? I have been
> > playing with facet_grid all day... Also the values on the axis are
> > incorrect ... they should be in the range 1:5 (as r above) and a plot
> > should show the probability density of that value...
> >
> > Any help? I am going mad...
> >
> > Many Thanks
> > JP
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
--
Jean-Paul Ebejer
Early Stage Researcher
InhibOx Ltd
Pembroke House
36-37 Pembroke Street
Oxford
OX1 1BP
UK
(+44 / 0) 1865 262 034
This email and any files transmitted with it are confide...{{dropped:22}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Dirichlet surface
On Mar 29, 2011, at 4:45 PM, Kehl Dániel wrote:
Dear list members,
I want to draw surfaces of Dirichlet distributions with different
parameter settings.
My code is the following:
#
a1 <- a2 <- a3 <- 2
#a2 <- .5
#a3 <- .5
x1 <- x2 <- seq(0.01, .99, by=.01)
f <- function(x1, x2){
term1 <- gamma(a1+a2+a3)/(gamma(a1)*gamma(a2)*gamma(a3))
term2 <- x1^(a1-1)*x2^(a2-1)*(1-x1-x2)^(a3-1)
term3 <- (x1 + x2 < 1)
term1*term2*term3
}
z <- outer(x1, x2, f)
z[z<=0] <- NA
persp(x1, x2, z,
main = "Dirichlet Distribution",
col = "lightblue",
theta = 50,
phi = 20,
r = 50,
d = 0.1,
expand = 0.5,
ltheta = 90,
lphi = 180,
shade = 0.75,
ticktype = "detailed",
nticks = 5)
#
It works fine (I guess), except for a1=a2=a3=1. In that case I get
the error: Error in persp.default... invalid 'z' limits.
The z matrix has only elements 2 and NA.
Might be a buglet in persp. If you set the zlim argument to c(0,2.1),
you get the expected constant plane at z=2 for those arguments.
Any ideas are appreciated.
Thank you:
Daniel
University of Pécs
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Connect observed values by a smooth curve, is there any method to get coordinate value of arbitrary point on the curve?
On Mar 29, 2011, at 9:32 PM, Zhipeng Wang wrote: Dear all, I have the problem as the title shows. It is time series of measurements. I want to estimate the value at time point when we have no actual measured value (Only in the time range, not for prediction in future time). Hope you could give me some suggestion or hint. Thank you so much. ?approxfun -- David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] glm: modelling zeros as binary and non-zeroes as coming from a continuous distribution
Hi: You might want to consider hurdle models in the pscl package. HTH, Dennis On Wed, Mar 30, 2011 at 2:41 AM, a11msp wrote: > Hello, > > I'd like to implement a regression model for extremely zero-inflated > continuous data using a conditional approach, whereby zeroes are > modelled as coming from a binary distribution, while non-zero values > are modelled as log-normal. > > So far, I've come across two solutions for this: one, in R, is > described in the book by Gelman & Hill > (http://www.amazon.com/dp/052168689X), where they just model zeros and > non-zeros separately and then bring them together by simulation. I can > do this, but it makes it difficult to assess the significance of > regression coefficients wrt to zero and each other. > > Another solution I have been pointed at is in SAS: > http://listserv.uga.edu/cgi-bin/wa?A2=ind0805A&L=sas-l&P=R20779, > where they use NLMIXED (with only fixed effects) to specify their own > log-likelihood function. > I'm wondering if there's any way to do the same in R (lme can't deal > with this, as far as I'm aware). > > Finally, I'm wondering whether anyone has experience with the COZIGAM > package - does it do something like this? > > Many thanks, > Mikhail > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/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] glm: modelling zeros as binary and non-zeroes as coming from a continuous distribution
Hi Dennis, Thanks - these were the first things I tried, but the problem is that they refuse to work with non-count data... Mikhail On Wed, Mar 30, 2011 at 12:56 PM, Dennis Murphy wrote: > Hi: > > You might want to consider hurdle models in the pscl package. > > HTH, > Dennis > > On Wed, Mar 30, 2011 at 2:41 AM, a11msp wrote: >> >> Hello, >> >> I'd like to implement a regression model for extremely zero-inflated >> continuous data using a conditional approach, whereby zeroes are >> modelled as coming from a binary distribution, while non-zero values >> are modelled as log-normal. >> >> So far, I've come across two solutions for this: one, in R, is >> described in the book by Gelman & Hill >> (http://www.amazon.com/dp/052168689X), where they just model zeros and >> non-zeros separately and then bring them together by simulation. I can >> do this, but it makes it difficult to assess the significance of >> regression coefficients wrt to zero and each other. >> >> Another solution I have been pointed at is in SAS: >> http://listserv.uga.edu/cgi-bin/wa?A2=ind0805A&L=sas-l&P=R20779, >> where they use NLMIXED (with only fixed effects) to specify their own >> log-likelihood function. >> I'm wondering if there's any way to do the same in R (lme can't deal >> with this, as far as I'm aware). >> >> Finally, I'm wondering whether anyone has experience with the COZIGAM >> package - does it do something like this? >> >> Many thanks, >> Mikhail >> >> __ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/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] glm: modelling zeros as binary and non-zeroes as coming from a continuous distribution
On Mar 30, 2011; 11:41am Mikhail wrote: >> I'm wondering if there's any way to do the same in R (lme can't deal >> with this, as far as I'm aware). You can do this using the pscl package. Regards, Mark. -- View this message in context: http://r.789695.n4.nabble.com/glm-modelling-zeros-as-binary-and-non-zeroes-as-coming-from-a-continuous-distribution-tp3417718p3417857.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] wrong calculation of R-squared in summary.lm
Dear all, I just stumbled upon the fact, that when I perform a regression on multivariate responses, that are not centred, I get a ricilulously high R-squared value. After reading the code of summary.lm, I found a bug in the function summary.lm: mss is calculated by: mss <-sum((f - mean(f))^2) - where f are the fitted values. This works only for a single response variable, because otherwise the matrix containing the fitted values is scaled by a scalar. replacing this with: mss <- scale(f,scale=FALSE) solved the problem for me. Example for bug: response<-cbind(rnorm(500)+1000,rnorm(500)+300) predictor<-rnorm(500) fit<-lm(response ~ predictor) summary.lm(fit) now reports a R^2 value of 1!!! Please correct me, if I'm wrong Stefan Stefan Schlager M.A. Anthropologie Medizinische Fakultät der der Albert Ludwigs- Universität Freiburg Hebelstr. 29 79104 Freiburg Anthropology Faculty of Medicine, Albert-Ludwigs-University Freiburg Hebelstr. 29 D- 79104 Freiburg phone +49 (0)761 203-5522 fax +49 (0)761 203-6898 On 30/03/11 13:54, r-help-requ...@r-project.org wrote: Mailing list subscription confirmation notice for mailing list R-help We have received a request from 129.132.148.130 for subscription of your email address, "stefan.schla...@uniklinik-freiburg.de", to the r-help@r-project.org mailing list. To confirm that you want to be added to this mailing list, simply reply to this message, keeping the Subject: header intact. Or visit this web page: https://stat.ethz.ch/mailman/confirm/r-help/4aa823dc72684b2ac88f8b630c5669c89da37c7c Or include the following line -- and only the following line -- in a message to r-help-requ...@r-project.org: confirm 4aa823dc72684b2ac88f8b630c5669c89da37c7c Note that simply sending a `reply' to this message should work from most mail readers, since that usually leaves the Subject: line in the right form (additional "Re:" text in the Subject: is okay). If you do not wish to be subscribed to this list, please simply disregard this message. If you think you are being maliciously subscribed to the list, or have any other questions, send them to r-help-ow...@r-project.org. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] glm: modelling zeros as binary and non-zeroes as coming from a continuous distribution
Update: turns out there was a sister posting to mine two years ago: http://r.789695.n4.nabble.com/Zinb-for-Non-interger-data-td898206.html It was then suggested to use a zero-inflated distribution from the gamlss package. It turns out that they do have a zero-adjusted (albeit not strictly speaking inflated) logarithmic distribution. Will try it out and see how it goes... Mikhail __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] fonts in mosaic
I need to change the font(s) used in mosaic from package vcd. The help file and the vignette do not give very explicit examples for doing that. The easiest solution would be changing the font for everything on the graph: var labels, var names, title, subtitle, and cell labels. What is the easiest way of accomplishing this? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] maximum likelihood accuracy - comparison with Stata
On 28 March 2011 17:08, Peter Ehlers wrote: > On 2011-03-27 21:37, Alex Olssen wrote: >> >> Hi everyone, >> >> I am looking to do some manual maximum likelihood estimation in R. I >> have done a lot of work in Stata and so I have been using output >> comparisons to get a handle on what is happening. >> >> I estimated a simple linear model in R with lm() and also my own >> maximum likelihood program. I then compared the output with Stata. >> Two things jumped out at me. >> >> Firstly, in Stata my coefficient estimates are identical in both the >> OLS estimation by -reg- and the maximum likelihood estimation using >> Stata's ml lf command. >> In R my coefficient estimates differ slightly between lm() and my >> own maximum likelihood estimation. >> >> Secondly, the estimates for sigma2 are very different between R >> and Stata. 3.14 in R compared to 1.78 in Stata. >> >> I have copied my maximum likelihood program below. It is copied from >> a great intro to MLE in R by Macro Steenbergen >> http://artsci.wustl.edu/~jmonogan/computing/r/MLE_in_R.pdf >> >> Any comments are welcome. In particular I would like to know why the >> estimate of sigma2 is so different. I would also like to know >> about the accuracy of the coefficient estimates. > > Some comments: > > 1. > Since Kmenta is not in the datasets package, you should mention > where you're getting it. (I suppose that for economists, it's > obvious, but we're not all economists.) I used the version in > the systemfit package. > > 2. > I don't believe your Stata value for sigma2 (by which I assume > you mean sigma^2). Are you quoting sigma? > > 3. > I can't reproduce your R value of 3.14 for sigma2. I get 3.725391. > > 4. > (most important) There's a difference between the simple ML > estimate (which is biased) and R's unbiased estimate for sigma^2. > This dataset has 20 cases, so try multiplying your result by 20/17. > > 5. > As to any difference in coefficients, I would guess that the > difference is slight (you don't show what it is); it > may be due to the fact that optim() produces approximate > solutions, whereas in this simple case, an 'exact' solution > is possible (and found by R). > > 6. > In case you weren't aware of it: the stats4 package has an > mle() function. ... and there is the "maxLik" package. http://cran.r-project.org/package=maxLik http://www.springerlink.com/content/973512476428614p/ Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] save ordinary numerical calculations as pdf
I'm not sure about .pdf, but look at ?sink to create a text file. You could
then print this as a .pdf if you so desired.
Rob
--
Robert W. Baer, Ph.D.
Professor of Physiology
Kirksville College of Osteopathic Medicine
A. T. Still University of Health Sciences
800 W. Jefferson St.
Kirksville, MO 63501
660-626-2322
FAX 660-626-2965
--
From: "Maas James Dr (MED)"
Sent: Wednesday, March 30, 2011 4:51 AM
To:
Subject: [R] save ordinary numerical calculations as pdf
I'd like to save some calculation outputs as a pdf, to incorporate with
others in a document. I've tried
pdf("filename")
name_of_object_to_output
dev.off()
but it doesn't seem to work, appears that this pdf function is for
graphics?
Is there a way to output numerical objects to pdf?
Thanks
J
Dr. Jim Maas
University of East Anglia
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] wrong calculation of R-squared in summary.lm
This is user error: that fit has class c("mlm", "lm") and you should
be calling summary() on it, which will dispatch the mlm method.
For multiple-response linear models, there is also summary.manova.
On Wed, 30 Mar 2011, Stefan Schlager wrote:
Dear all,
I just stumbled upon the fact, that when I perform a regression on
multivariate responses, that are not centred, I get a ricilulously high
R-squared value. After reading the code of summary.lm, I found a bug in the
function summary.lm:
mss is calculated by:
mss <-sum((f - mean(f))^2) - where f are the fitted values.
This works only for a single response variable, because otherwise the matrix
containing the fitted values is scaled by a scalar. replacing
this with:
mss <- scale(f,scale=FALSE) solved the problem for me.
Example for bug:
response<-cbind(rnorm(500)+1000,rnorm(500)+300)
predictor<-rnorm(500)
fit<-lm(response ~ predictor)
summary.lm(fit) now reports a R^2 value of 1!!!
Please correct me, if I'm wrong
Stefan
Stefan Schlager M.A.
Anthropologie
Medizinische Fakultät der der Albert Ludwigs- Universität Freiburg
Hebelstr. 29
79104 Freiburg
Anthropology
Faculty of Medicine, Albert-Ludwigs-University Freiburg
Hebelstr. 29
D- 79104 Freiburg
phone +49 (0)761 203-5522
fax +49 (0)761 203-6898
On 30/03/11 13:54, r-help-requ...@r-project.org wrote:
Mailing list subscription confirmation notice for mailing list R-help
We have received a request from 129.132.148.130 for subscription of
your email address, "stefan.schla...@uniklinik-freiburg.de", to the
r-help@r-project.org mailing list. To confirm that you want to be
added to this mailing list, simply reply to this message, keeping the
Subject: header intact. Or visit this web page:
https://stat.ethz.ch/mailman/confirm/r-help/4aa823dc72684b2ac88f8b630c5669c89da37c7c
Or include the following line -- and only the following line -- in a
message to r-help-requ...@r-project.org:
confirm 4aa823dc72684b2ac88f8b630c5669c89da37c7c
Note that simply sending a `reply' to this message should work from
most mail readers, since that usually leaves the Subject: line in the
right form (additional "Re:" text in the Subject: is okay).
If you do not wish to be subscribed to this list, please simply
disregard this message. If you think you are being maliciously
subscribed to the list, or have any other questions, send them to
r-help-ow...@r-project.org.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Brian D. Ripley, rip...@stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Quick recode of -999 to NA in R
Hi,
I am trying to write a loop to recode my data from -999 to NA in R. What's
the most efficient way to do this? Below is what I'm presently doing, which
is inefficient. Thanks,
Chris
dat0 <- read.table("time1.dat")
colnames(dat0) <- c("e1dq", "e1arcp", "e1dev", "s1prcp", "s1nrcp", "s1ints",
"a1gpar", "a1pias", "a1devt")
dat0[dat0$e1dq==-999.,"e1dq"] <- NA
dat0[dat0$e1arcp==-999.,"e1arcp"] <- NA
dat0[dat0$e1dev==-999.,"e1dev"] <- NA
dat0[dat0$s1prcp==-999.,"s1prcp"] <- NA
dat0[dat0$s1nrcp==-999.,"s1nrcp"] <- NA
dat0[dat0$s1ints==-999.,"s1ints"] <- NA
dat0[dat0$a1gpar==-999.,"a1gpar"] <- NA
dat0[dat0$a1pias==-999.,"a1pias"] <- NA
dat0[dat0$a1devt==-999.,"a1devt"] <- NA
[[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] Quick recode of -999 to NA in R
Try:
dat0 <- read.table('tim1.dat', na = -999)
On Wed, Mar 30, 2011 at 10:15 AM, Christopher Desjardins
wrote:
> Hi,
> I am trying to write a loop to recode my data from -999 to NA in R. What's
> the most efficient way to do this? Below is what I'm presently doing, which
> is inefficient. Thanks,
> Chris
>
>
> dat0 <- read.table("time1.dat")
>
> colnames(dat0) <- c("e1dq", "e1arcp", "e1dev", "s1prcp", "s1nrcp", "s1ints",
> "a1gpar", "a1pias", "a1devt")
>
> dat0[dat0$e1dq==-999.,"e1dq"] <- NA
>
> dat0[dat0$e1arcp==-999.,"e1arcp"] <- NA
>
> dat0[dat0$e1dev==-999.,"e1dev"] <- NA
>
> dat0[dat0$s1prcp==-999.,"s1prcp"] <- NA
>
> dat0[dat0$s1nrcp==-999.,"s1nrcp"] <- NA
>
> dat0[dat0$s1ints==-999.,"s1ints"] <- NA
>
> dat0[dat0$a1gpar==-999.,"a1gpar"] <- NA
>
> dat0[dat0$a1pias==-999.,"a1pias"] <- NA
>
> dat0[dat0$a1devt==-999.,"a1devt"] <- NA
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Quick recode of -999 to NA in R
On 30.03.2011 15:15, Christopher Desjardins wrote:
Hi,
I am trying to write a loop to recode my data from -999 to NA in R. What's
the most efficient way to do this? Below is what I'm presently doing, which
is inefficient. Thanks,
Chris
I think read.table(na.string="-999.") is.
Uwe Ligges
dat0<- read.table("time1.dat")
colnames(dat0)<- c("e1dq", "e1arcp", "e1dev", "s1prcp", "s1nrcp", "s1ints",
"a1gpar", "a1pias", "a1devt")
dat0[dat0$e1dq==-999.,"e1dq"]<- NA
dat0[dat0$e1arcp==-999.,"e1arcp"]<- NA
dat0[dat0$e1dev==-999.,"e1dev"]<- NA
dat0[dat0$s1prcp==-999.,"s1prcp"]<- NA
dat0[dat0$s1nrcp==-999.,"s1nrcp"]<- NA
dat0[dat0$s1ints==-999.,"s1ints"]<- NA
dat0[dat0$a1gpar==-999.,"a1gpar"]<- NA
dat0[dat0$a1pias==-999.,"a1pias"]<- NA
dat0[dat0$a1devt==-999.,"a1devt"]<- NA
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Quick recode of -999 to NA in R
Am 30.03.2011 09:15, schrieb Christopher Desjardins:
Hi,
I am trying to write a loop to recode my data from -999 to NA in R. What's
the most efficient way to do this? Below is what I'm presently doing, which
is inefficient. Thanks,
Chris
dat0<- read.table("time1.dat")
colnames(dat0)<- c("e1dq", "e1arcp", "e1dev", "s1prcp", "s1nrcp", "s1ints",
"a1gpar", "a1pias", "a1devt")
dat0[dat0$e1dq==-999.,"e1dq"]<- NA
dat0[dat0$e1arcp==-999.,"e1arcp"]<- NA
dat0[dat0$e1dev==-999.,"e1dev"]<- NA
dat0[dat0$s1prcp==-999.,"s1prcp"]<- NA
dat0[dat0$s1nrcp==-999.,"s1nrcp"]<- NA
dat0[dat0$s1ints==-999.,"s1ints"]<- NA
dat0[dat0$a1gpar==-999.,"a1gpar"]<- NA
dat0[dat0$a1pias==-999.,"a1pias"]<- NA
dat0[dat0$a1devt==-999.,"a1devt"]<- NA
Given that all your variables share the same code for missing values,
e.g. -99, you could do the following:
df <- data.frame(a = c(1, 2, -99),
b = c(33, -99, 22),
c = c(-99, -99, 101))
df
df[df == -99] <- NA
df
HTH,
Bernd
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Quick recode of -999 to NA in R
On Wed, 2011-03-30 at 08:15 -0500, Christopher Desjardins wrote:
> Hi,
> I am trying to write a loop to recode my data from -999 to NA in R. What's
> the most efficient way to do this? Below is what I'm presently doing, which
> is inefficient. Thanks,
> Chris
>
>
>dat0 <- read.table("time1.dat")
dat0 <- read.table("time1.dat", na.strings = c("NA","-999"))
would seem the easiest option.
Do read ?read.table for details.
HTH
G
> colnames(dat0) <- c("e1dq", "e1arcp", "e1dev", "s1prcp", "s1nrcp", "s1ints",
> "a1gpar", "a1pias", "a1devt")
>
> dat0[dat0$e1dq==-999.,"e1dq"] <- NA
>
> dat0[dat0$e1arcp==-999.,"e1arcp"] <- NA
>
> dat0[dat0$e1dev==-999.,"e1dev"] <- NA
>
> dat0[dat0$s1prcp==-999.,"s1prcp"] <- NA
>
> dat0[dat0$s1nrcp==-999.,"s1nrcp"] <- NA
>
> dat0[dat0$s1ints==-999.,"s1ints"] <- NA
>
> dat0[dat0$a1gpar==-999.,"a1gpar"] <- NA
>
> dat0[dat0$a1pias==-999.,"a1pias"] <- NA
>
> dat0[dat0$a1devt==-999.,"a1devt"] <- NA
>
> [[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.
--
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Dr. Gavin Simpson [t] +44 (0)20 7679 0522
ECRC, UCL Geography, [f] +44 (0)20 7679 0565
Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/
UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Quick recode of -999 to NA in R
It could be done in a large number of ways depending on how often you
need it etc.
You might take a look at defmacro in package gtools:
# library(gtools)
setNA <- macro(df, var, values)
{
df$var[df$var %in% values] <- NA
}
then instead of
> dat0[dat0$e1dq==-999.,"e1dq"] <- NA
you could have
setNA(dat0, e1dq, -999)
Otherwise, with lapply, you could have
dat0.na_ized <- lapply(dat0, function(x) x[x%in% 999.0] <- NA)
(but I didn't try this so I can't tell if you need to use
as.data.frame on the result or not)
Kenn
On Wed, Mar 30, 2011 at 4:15 PM, Christopher Desjardins
wrote:
> Hi,
> I am trying to write a loop to recode my data from -999 to NA in R. What's
> the most efficient way to do this? Below is what I'm presently doing, which
> is inefficient. Thanks,
> Chris
>
>
> dat0 <- read.table("time1.dat")
>
> colnames(dat0) <- c("e1dq", "e1arcp", "e1dev", "s1prcp", "s1nrcp", "s1ints",
> "a1gpar", "a1pias", "a1devt")
>
> dat0[dat0$e1dq==-999.,"e1dq"] <- NA
>
> dat0[dat0$e1arcp==-999.,"e1arcp"] <- NA
>
> dat0[dat0$e1dev==-999.,"e1dev"] <- NA
>
> dat0[dat0$s1prcp==-999.,"s1prcp"] <- NA
>
> dat0[dat0$s1nrcp==-999.,"s1nrcp"] <- NA
>
> dat0[dat0$s1ints==-999.,"s1ints"] <- NA
>
> dat0[dat0$a1gpar==-999.,"a1gpar"] <- NA
>
> dat0[dat0$a1pias==-999.,"a1pias"] <- NA
>
> dat0[dat0$a1devt==-999.,"a1devt"] <- NA
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] save ordinary numerical calculations as pdf
Take a look in ?Sweave
On Wed, Mar 30, 2011 at 6:51 AM, Maas James Dr (MED) wrote:
> I'd like to save some calculation outputs as a pdf, to incorporate with
> others in a document. I've tried
>
> pdf("filename")
> name_of_object_to_output
> dev.off()
>
> but it doesn't seem to work, appears that this pdf function is for graphics?
>
> Is there a way to output numerical objects to pdf?
>
> Thanks
>
> J
>
> Dr. Jim Maas
> University of East Anglia
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to regress data into coefficients for a gamma function
Regression for the gamma distribution can be expressed as a generalized linear model. Check Chapter 8 of McCullagh, P. & Nelder, J. A. (1989), Generalized linear models, Chapman & Hall, London, UK. Walter Anderson Sent by: r-help-boun...@r-project.org 03/29/2011 09:57 AM To r-help@r-project.org cc Subject [R] How to regress data into coefficients for a gamma function Hello, I need to regress data like the example below. The data points represent friction factors derived from observed trip length data. The function used to describe this data is a gamma function of the form, f(t) = a * t^b * e^(c*t) and I need to regress the data to obtain the a,b, and c coefficients. The gamma function can also be expressed in the log-linear form, ln[f(t)] = ln[a] + b * ln[t] + c*t, which may be easier to perform a regression on. I have performed a search for information on the subject, and have found a few possibly related sites, I can not figure out how to perform the regression in R. Any help/guidance would be appreciated. Walter tf ln(f) 1 195289313 2 295107713 3 394599113 4 494235813 5 593945213 6 689549413 7 789186113 8 888241613 9 987369713 101086788513 111179777113 121279159513 131377924313 141477306813 151576580213 161663538313 171762848113 181862048813 191961285913 202060959013 212150859713 222250351113 232349842513 242449006913 252548607313 262642286212 272741886612 282841523312 292940833112 303040542412 313126338012 323226120112 333325720412 343425284512 353525066512 363620997812 373720598112 383820343812 393919907912 404019726312 414115475811 424215294211 434314930911 444414858311 454514749311 46469808611 47479772311 48489590611 49499409011 50509336311 51518173811 52528101211 53538064911 54547955911 55557919511 56567047711 57577047711 58587047711 59597047711 60607047711 61614286710 62624214010 63634214010 64644214010 65654177710 66663669110 67673632810 68683596510 69693596510 70703596510 71713196810 72723196810 73733196810 74743196810 75753160510 7676210709 7777207079 7878207079 7979203439 8080199809 8181159849 8282159849 8383159849 8484159849 8585159849 8686141689 8787141689 8888141689 8989141689 9090141689 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Not all rows are being read-in
Philipp, you are a savior! That's exactly what has been happening - and it was driving me crazy. quote="" fixed things. Thank you very much! Dimitri On Wed, Mar 30, 2011 at 5:01 AM, Philipp Pagel wrote: > On Tue, Mar 29, 2011 at 06:58:59PM -0400, Dimitri Liakhovitski wrote: >> I have a tab-delimited .txt file (size 800MB) with about 3.4 million >> rows and 41 columns. About 15 columns contain strings. >> Tried to read it in in R 2.12.2 on a laptop that has Windows XP: >> mydata<-read.delim(file="FileName.TXT",sep="\t") >> R did not complain (!) and I got: dim(mydata) 1692063 41. > > My guess would be that there are (unexpected) quotes and/or double quotes in > your > file and so R thinks that rather large blocks of your file are > actually very long strings. This routinely happens in situations like > this: > > ID x description > 1 0.4 my first measurement > 2 1.6 Normal 5" object > 3 0.4 Some measuremetn > 4 0.7 A 4" long sample > > R thinks that the description in row 2 ends in row 4 and you loose > data. > > Try read.delim(..., quote=""). > > cu > Philipp > > -- > Dr. Philipp Pagel > Lehrstuhl für Genomorientierte Bioinformatik > Technische Universität München > Wissenschaftszentrum Weihenstephan > Maximus-von-Imhof-Forum 3 > 85354 Freising, Germany > http://webclu.bio.wzw.tum.de/~pagel/ > > __ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > -- Dimitri Liakhovitski Ninah Consulting www.ninah.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] save ordinary numerical calculations as pdf
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1
On 30/03/11 11:51, Maas James Dr (MED) wrote:
> I'd like to save some calculation outputs as a pdf, to incorporate with
> others in a document. I've tried
>
> pdf("filename")
> name_of_object_to_output
> dev.off()
>
> but it doesn't seem to work, appears that this pdf function is for graphics?
>
> Is there a way to output numerical objects to pdf?
Yes - there is an equivalent to the sink() which pute the results in a
pdf, but I can not remember the name of the function.
I know it does not help much, but at least I can tell you that there is
such a function.
Cheers,
Rainer
>
> Thanks
>
> J
>
> Dr. Jim Maas
> University of East Anglia
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
- --
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
Biology, UCT), Dipl. Phys. (Germany)
Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa
Tel:+33 - (0)9 53 10 27 44
Cell: +27 - (0)8 39 47 90 42
Fax (SA): +27 - (0)8 65 16 27 82
Fax (D) : +49 - (0)3 21 21 25 22 44
Fax (FR): +33 - (0)9 58 10 27 44
email: rai...@krugs.de
Skype: RMkrug
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
iEYEARECAAYFAk2TMzAACgkQoYgNqgF2egrATwCcCrcCfRBkj+Q8Cb+RddIlyfO6
Fg4An2RKwh0GoCRJJ0VJPt8GIQPTcl3B
=GNAp
-END PGP 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] Quick recode of -999 to NA in R
Try using a loop like the following
dat0 <- read.table("time1.dat")
id <- c("e1dq", "e1arcp", "e1dev", "s1prcp", "s1nrcp","s1ints","a1gpar", "a1pias",
"a1devt")
for (a in 1:length(id)) {
dat0[dat0$id[a]==-999.,as.character(id[a])] <- NA
}
--
Muhammad Rahiz
Researcher & DPhil Candidate (Climate Systems & Policy)
School of Geography & the Environment
University of Oxford
On Wed, 30 Mar 2011, Christopher Desjardins wrote:
Hi,
I am trying to write a loop to recode my data from -999 to NA in R. What's
the most efficient way to do this? Below is what I'm presently doing, which
is inefficient. Thanks,
Chris
dat0 <- read.table("time1.dat")
colnames(dat0) <- c("e1dq", "e1arcp", "e1dev", "s1prcp", "s1nrcp", "s1ints",
"a1gpar", "a1pias", "a1devt")
dat0[dat0$e1dq==-999.,"e1dq"] <- NA
dat0[dat0$e1arcp==-999.,"e1arcp"] <- NA
dat0[dat0$e1dev==-999.,"e1dev"] <- NA
dat0[dat0$s1prcp==-999.,"s1prcp"] <- NA
dat0[dat0$s1nrcp==-999.,"s1nrcp"] <- NA
dat0[dat0$s1ints==-999.,"s1ints"] <- NA
dat0[dat0$a1gpar==-999.,"a1gpar"] <- NA
dat0[dat0$a1pias==-999.,"a1pias"] <- NA
dat0[dat0$a1devt==-999.,"a1devt"] <- NA
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] save ordinary numerical calculations as pdf
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1
On 30/03/11 15:42, Rainer M Krug wrote:
> On 30/03/11 11:51, Maas James Dr (MED) wrote:
>> I'd like to save some calculation outputs as a pdf, to incorporate with
>> others in a document. I've tried
>
>> pdf("filename")
>> name_of_object_to_output
>> dev.off()
>
>> but it doesn't seem to work, appears that this pdf function is for graphics?
>
>> Is there a way to output numerical objects to pdf?
>
> Yes - there is an equivalent to the sink() which pute the results in a
> pdf, but I can not remember the name of the function.
>
> I know it does not help much, but at least I can tell you that there is
> such a function.
textplot in the gplots package:
This function displays text output in a graphics window. It is
the equivalent of 'print' except that the output is displayed as a
plot.
Thanks google and Greg Snow
(http://www.mail-archive.com/r-help@r-project.org/msg02187.html)
Rainer
>
> Cheers,
>
> Rainer
>
>
>> Thanks
>
>> J
>
>> Dr. Jim Maas
>> University of East Anglia
>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>
- --
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
Biology, UCT), Dipl. Phys. (Germany)
Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa
Tel:+33 - (0)9 53 10 27 44
Cell: +27 - (0)8 39 47 90 42
Fax (SA): +27 - (0)8 65 16 27 82
Fax (D) : +49 - (0)3 21 21 25 22 44
Fax (FR): +33 - (0)9 58 10 27 44
email: rai...@krugs.de
Skype: RMkrug
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
iEYEARECAAYFAk2TNEQACgkQoYgNqgF2egp2qwCfbVJfyetIkxyzqcRuTASrDSxt
qtUAnjJVuCYzjozwOfhCqYcZzz3ZBJow
=9/4z
-END PGP 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] Quick recode of -999 to NA in R
Ah ... yes. I knew that but clearly didn't at the time of my question or
script writing.
Thanks,
Chris
On Wed, Mar 30, 2011 at 8:22 AM, Henrique Dallazuanna wrote:
> Try:
>
> dat0 <- read.table('tim1.dat', na = -999)
>
Ah ... yes. I knew that but clearly didn't at the time of my question or
script writing.
Thanks,
Chris
>
> On Wed, Mar 30, 2011 at 10:15 AM, Christopher Desjardins
> wrote:
> > Hi,
> > I am trying to write a loop to recode my data from -999 to NA in R.
> What's
> > the most efficient way to do this? Below is what I'm presently doing,
> which
> > is inefficient. Thanks,
> > Chris
> >
> >
> > dat0 <- read.table("time1.dat")
> >
> > colnames(dat0) <- c("e1dq", "e1arcp", "e1dev", "s1prcp", "s1nrcp",
> "s1ints",
> > "a1gpar", "a1pias", "a1devt")
> >
> > dat0[dat0$e1dq==-999.,"e1dq"] <- NA
> >
> > dat0[dat0$e1arcp==-999.,"e1arcp"] <- NA
> >
> > dat0[dat0$e1dev==-999.,"e1dev"] <- NA
> >
> > dat0[dat0$s1prcp==-999.,"s1prcp"] <- NA
> >
> > dat0[dat0$s1nrcp==-999.,"s1nrcp"] <- NA
> >
> > dat0[dat0$s1ints==-999.,"s1ints"] <- NA
> >
> > dat0[dat0$a1gpar==-999.,"a1gpar"] <- NA
> >
> > dat0[dat0$a1pias==-999.,"a1pias"] <- NA
> >
> > dat0[dat0$a1devt==-999.,"a1devt"] <- NA
> >
> >[[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
>
>
> --
> Henrique Dallazuanna
> Curitiba-Paraná-Brasil
> 25° 25' 40" S 49° 16' 22" O
>
[[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] save ordinary numerical calculations as pdf
The latex() function in the Hmisc package will typeset your
objects. Embed that in a tex document and run pdflatex.
Rich
Sent from my iPhone
On Mar 30, 2011, at 5:51, "Maas James Dr (MED)" wrote:
> I'd like to save some calculation outputs as a pdf, to incorporate with
> others in a document. I've tried
>
> pdf("filename")
> name_of_object_to_output
> dev.off()
>
> but it doesn't seem to work, appears that this pdf function is for graphics?
>
> Is there a way to output numerical objects to pdf?
>
> Thanks
>
> J
>
> Dr. Jim Maas
> University of East Anglia
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] summing values by week - based on daily dates - but with some dates missing
Dear everybody,
I have the following challenge. I have a data set with 2 subgroups,
dates (days), and corresponding values (see example code below).
Within each subgroup: I need to aggregate (sum) the values by week -
for weeks that start on a Monday (for example, 2008-12-29 was a
Monday).
I find it difficult because I have missing dates in my data - so that
sometimes I don't even have the date for some Mondays. So, I can't
write a proper loop.
I want my output to look something like this:
group dates value
group.1 2008-12-29 3.0937
group.1 2009-01-05 3.8833
group.1 2009-01-12 1.362
...
group.2 2008-12-29 2.250
group.2 2009-01-05 1.4057
group.2 2009-01-12 3.4411
...
Thanks a lot for your suggestions! The code is below:
Dimitri
### Creating example data set:
mydates<-rep(seq(as.Date("2008-12-29"), length = 43, by = "day"),2)
myfactor<-c(rep("group.1",43),rep("group.2",43))
set.seed(123)
myvalues<-runif(86,0,1)
myframe<-data.frame(dates=mydates,group=myfactor,value=myvalues)
(myframe)
dim(myframe)
## Removing same rows (dates) unsystematically:
set.seed(123)
removed.group1<-sample(1:43,size=11,replace=F)
set.seed(456)
removed.group2<-sample(44:86,size=11,replace=F)
to.remove<-c(removed.group1,removed.group2);length(to.remove)
to.remove<-to.remove[order(to.remove)]
myframe<-myframe[-to.remove,]
(myframe)
--
Dimitri Liakhovitski
Ninah Consulting
www.ninah.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] Using xlevels
I'm working on predict.survreg and am confused about xlevels. The model.frame method has the argument, but none of the standard methods (model.frame.lm, model.frame.glm) appear to make use of it. The documentation for model.matrix states: xlev: to be used as argument of model.frame if data has no "terms" attribute. But the terms attribute has no xlevels information in it, so I find this statement completely confusing. Any insight is appreciated. Terry Therneau __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] fonts in mosaic
On Wed, 30 Mar 2011, Erich Neuwirth wrote:
I need to change the font(s) used in mosaic from package vcd.
The help file and the vignette do not give very explicit examples for
doing that.
The easiest solution would be changing the font for everything
on the graph: var labels, var names, title, subtitle, and cell labels.
What is the easiest way of accomplishing this?
Personally, I simply change the size of the device I'm plotting on. When I
plot on a large device, the fonts will be relatively smaller, and vice
versa. This is what I do when including graphics in PDF files (papers,
slides, reports, etc.).
For fine control, you can set the arguments of the labeling function
employed. ?strucplot shows that the default is ?labeling_border which has
several arguments. For example you can set the graphical parameters of the
labels (gp_labels) or the graphical parameters of the variable names
(gp_varnames). Both arguments take ?gpar lists ("grid" graphical
parameters). For example you may do
mosaic(UCBAdmissions)
and
mosaic(UCBAdmissions, labeling_args = list(
gp_labels = gpar(fontsize = 12, fontface = 3),
gp_varnames = gpar(fontsize = 16, fontface = 2)
))
etc.
Hope that helps,
Z
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] sampling design runs with no errors but returns empty data set
Dear colleagues,
I'm working with the 2008 Canada Election Studies
(http://www.queensu.ca/cora/_files/_CES/CES2008.sav.zip), trying to construct a
weighted national sample using the survey package.
Three weights are included in the national survey (a household weight, a
provincial weight and a national weight which is a product of the first two).
In the following code I removed variables with missing national weights and
tried to construct the sample from advice I've gleaned from the documentation
for the survey package and other help requests.
There are no errors, but the data frame (weight_test) contains no
What am I missing?
Yours, Simon Kiss
P.S. The code is only reproducible if the data set is downloadable. I'm nt sure
ces<-read.spss(file.choose(), to.data.frame=TRUE, use.value.labels=FALSE)
missing_data<-subset(ces1, !is.na(ces08_NATWGT))
weight_test<-svydesign(id=~0, weights=~ces08_NATWGT, data=missing_data)
Note: this is some reproducible code that creates a data set that is a very
stripped down version of what I'm working with, but with this, the surveydesign
function appears to work properly.
mydat<-data.frame(ces08_HHWGT=runif(3000, 0.5, 5), ces08_PROVWGT=runif(3000,
0.6, 1.2), party=sample(c("NDP", "BQ", "Lib", "Con"), 3000, replace=TRUE),
age=sample(seq(18, 72,1), 3000, replace=TRUE), income=sample(seq(21,121,1),
3000, replace=TRUE))
mydat$ces08_NATWGT<-mydat$ces08_HHWGT*mydat$ces08_PROVWGT
weight_test<-svydesign(id=~1, weights=~ces08_NATWGT, data=mydat)
*
Simon J. Kiss, PhD
Assistant Professor, Wilfrid Laurier University
73 George Street
Brantford, Ontario, Canada
N3T 2C9
Cell: +1 519 761 7606
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] VECM with UNRESTRICTED TREND
Dear All, My question is: how can I estimate VECM system with "unrestricted trend" (aka "case 5") option as a deterministic term? As far as I know, ca.jo in urca package allows for "restricted trend" only [vecm <- ca.jo(data, type = "trace"/"eigen", ecdet = "trend", K = n, spec = "transitory"/"longrun")]. Obviously, I don't have to do this in urca, so if another package gives the possibility, please let me know too! Thanks in advance! Greg [[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] Multiple area plots to share the same x-axis
I'm going to go ahead and shamelessly bump this question up the list.
I saw that out of the 34 posts yesterday, only 9 did not receive an
answer. Perhaps someone who finds this question trivial is checking
his e-mail right now :)
Werner
On Tue, Mar 29, 2011 at 10:20 AM, jovian wrote:
> Hello,
>
> I asked a similar question before but in an existing thread. I am not sure
> if it is proper etiquette to repost a similar question as a new tread, I
> think in this case, it might be because this way more people can see it and
> perhaps learn from it. (Also, part of the existing thread became private)
>
> I want to know how to plot multiple ggplot area plots on top of one another
> so that the same x-axis is shared?
>
> This solution simply stitches multiple plots on top of each other:
>
> vp.layout <- function(x, y) viewport(layout.pos.row=x, layout.pos.col=y)
> arrange <- function(..., nrow=NULL, ncol=NULL, as.table=FALSE) {
> dots <- list(...)
> n <- length(dots)
> if(is.null(nrow) & is.null(ncol)) { nrow = floor(n/2) ; ncol =
> ceiling(n/nrow)}
> if(is.null(nrow)) { nrow = ceiling(n/ncol)}
> if(is.null(ncol)) { ncol = ceiling(n/nrow)}
> ## NOTE see n2mfrow in grDevices for possible alternative
> grid.newpage()
> pushViewport(viewport(layout=grid.layout(nrow,ncol) ) )
> ii.p <- 1
> for(ii.row in seq(1, nrow)){
> ii.table.row <- ii.row
> if(as.table) {ii.table.row <- nrow - ii.table.row + 1}
> for(ii.col in seq(1, ncol)){
> ii.table <- ii.p
> if(ii.p > n) break
> print(dots[[ii.table]], vp=vp.layout(ii.table.row, ii.col))
> ii.p <- ii.p + 1
> }
> }
> }
>
> set <- read.table(file="http://www.jovian.nl/set.csv";, head=1, sep=",")
> set2 <- read.table(file="http://www.jovian.nl/set2.csv";, head=1, sep=",")
> library(ggplot2)
> s <- ggplot(set, aes(x = time, y = hours)) + geom_area(colour = 'red', fill
> = 'red', alpha = 0.5) +
> geom_area(stat = 'smooth', span = 0.2, alpha = 0.3) + ylim(0,40)
> s1 <- ggplot(set2, aes(x = time, y = hours)) + geom_area(colour = 'red',
> fill = 'red', alpha = 0.5) +
> geom_area(stat = 'smooth', span = 0.2, alpha = 0.3) + ylim(0,40)
> arrange(s,s1,ncol=1)
>
>
> The arrange() function was taken from
> http://gettinggeneticsdone.blogspot.com/2010/03/arrange-multiple-ggplot2-plots-in-same.html.
> In this example, the x-axes are only similar because the data sets have the
> same range. In effect nothing more happens than that two images are plotted
> on top of one another. Now how to "merge" these two (and later more) area
> plots on top of each other so that they share the same x-axis (so that only
> one x-axis would be necessary on the bottom of the plot)?
>
> Thanks,
>
> Werner
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Multiple-area-plots-to-share-the-same-x-axis-tp3414050p3414050.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] Using xlevels
On Wed, 30 Mar 2011, Terry Therneau wrote: I'm working on predict.survreg and am confused about xlevels. The model.frame method has the argument, but none of the standard methods (model.frame.lm, model.frame.glm) appear to make use of it. But I see this in predict.lm: m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$xlevels) It is used to remap levels in newdata to those used in the fit. The documentation for model.matrix states: xlev: to be used as argument of model.frame if data has no "terms" attribute. Well, the code says if (is.null(attr(data, "terms"))) data <- model.frame(object, data, xlev=xlev) But the terms attribute has no xlevels information in it, so I find this statement completely confusing. Any insight is appreciated. It means exactly what it says: a 'data' argument with a terms attribute is considered to be a model frame. Terry Therneau __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Brian D. Ripley, rip...@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multiple area plots to share the same x-axis
Use facetting:
sets <- rbind(cbind(set="set1", set), cbind(set="set2", set2))
ggplot(sets, aes(x = time, y = hours)) +
geom_area(colour = 'red', fill = 'red', alpha = 0.5) +
geom_area(stat = 'smooth', span = 0.2, alpha = 0.3) +
ylim(0,40) +
facet_grid(set ~ .)
Best,
Ista
On Wed, Mar 30, 2011 at 10:54 AM, Werner Heijstek wrote:
> I'm going to go ahead and shamelessly bump this question up the list.
> I saw that out of the 34 posts yesterday, only 9 did not receive an
> answer. Perhaps someone who finds this question trivial is checking
> his e-mail right now :)
>
> Werner
>
> On Tue, Mar 29, 2011 at 10:20 AM, jovian wrote:
>> Hello,
>>
>> I asked a similar question before but in an existing thread. I am not sure
>> if it is proper etiquette to repost a similar question as a new tread, I
>> think in this case, it might be because this way more people can see it and
>> perhaps learn from it. (Also, part of the existing thread became private)
>>
>> I want to know how to plot multiple ggplot area plots on top of one another
>> so that the same x-axis is shared?
>>
>> This solution simply stitches multiple plots on top of each other:
>>
>> vp.layout <- function(x, y) viewport(layout.pos.row=x, layout.pos.col=y)
>> arrange <- function(..., nrow=NULL, ncol=NULL, as.table=FALSE) {
>> dots <- list(...)
>> n <- length(dots)
>> if(is.null(nrow) & is.null(ncol)) { nrow = floor(n/2) ; ncol =
>> ceiling(n/nrow)}
>> if(is.null(nrow)) { nrow = ceiling(n/ncol)}
>> if(is.null(ncol)) { ncol = ceiling(n/nrow)}
>> ## NOTE see n2mfrow in grDevices for possible alternative
>> grid.newpage()
>> pushViewport(viewport(layout=grid.layout(nrow,ncol) ) )
>> ii.p <- 1
>> for(ii.row in seq(1, nrow)){
>> ii.table.row <- ii.row
>> if(as.table) {ii.table.row <- nrow - ii.table.row + 1}
>> for(ii.col in seq(1, ncol)){
>> ii.table <- ii.p
>> if(ii.p > n) break
>> print(dots[[ii.table]], vp=vp.layout(ii.table.row, ii.col))
>> ii.p <- ii.p + 1
>> }
>> }
>> }
>>
>> set <- read.table(file="http://www.jovian.nl/set.csv";, head=1, sep=",")
>> set2 <- read.table(file="http://www.jovian.nl/set2.csv";, head=1, sep=",")
>> library(ggplot2)
>> s <- ggplot(set, aes(x = time, y = hours)) + geom_area(colour = 'red', fill
>> = 'red', alpha = 0.5) +
>> geom_area(stat = 'smooth', span = 0.2, alpha = 0.3) + ylim(0,40)
>> s1 <- ggplot(set2, aes(x = time, y = hours)) + geom_area(colour = 'red',
>> fill = 'red', alpha = 0.5) +
>> geom_area(stat = 'smooth', span = 0.2, alpha = 0.3) + ylim(0,40)
>> arrange(s,s1,ncol=1)
>>
>>
>> The arrange() function was taken from
>> http://gettinggeneticsdone.blogspot.com/2010/03/arrange-multiple-ggplot2-plots-in-same.html.
>> In this example, the x-axes are only similar because the data sets have the
>> same range. In effect nothing more happens than that two images are plotted
>> on top of one another. Now how to "merge" these two (and later more) area
>> plots on top of each other so that they share the same x-axis (so that only
>> one x-axis would be necessary on the bottom of the plot)?
>>
>> Thanks,
>>
>> Werner
>>
>>
>> --
>> View this message in context:
>> http://r.789695.n4.nabble.com/Multiple-area-plots-to-share-the-same-x-axis-tp3414050p3414050.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.
>
--
Ista Zahn
Graduate student
University of Rochester
Department of Clinical and Social Psychology
http://yourpsyche.org
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Degrees of freedom for lm in logLik and AIC
On Mar 28, 2011, at 16:53 , Ben Bolker wrote: > Rubén Roa azti.es> writes: > >> >> >> However, shouldn't _free parameters_ only be counted for degrees of >> freedom and for calculation of AIC? >> The sigma parameter is profiled out in a least-squares >> linear regression, so it's not free, it's not a >> dimension of the likelihood. >> Just wondering ... >> > > For AIC I think this distinction should only matter for the purposes > of consistency between computations in different packages/languages/contexts. > Only the differences between AIC values matter for inference. (If you were > talking about AICc then I would tend to agree with you -- nuisance > parameters should not affect 'residual degrees of freedom'/finite-size > corrections.) > You're having me confused... AIC and friends are not within my "core competences" (i.e., I know what they are about, but I don't use them intensively on a daily basis), so I may be completely off base, but AFAICS _residual_ degrees of freedom never enters in the logLik computations. Also, I don't get the point about profiling -- you're just maximizing in two steps: first over sigma then over everything else, how is that different from just maximizing? On the other hand, I suppose it might be argued that the REML variant really only makes sense as a likelihood for sigma, so should have df=1. After all, it is by definition based on a transformation which makes the mean value parameters disappear. -- Peter Dalgaard Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] VECM with UNRESTRICTED TREND
Hello Greg,
you can exploit the argument 'dumvar' for this. See ?ca.jo
Best,
Bernhard
> -Ursprüngliche Nachricht-
> Von: r-help-boun...@r-project.org
> [mailto:r-help-boun...@r-project.org] Im Auftrag von Grzegorz Konat
> Gesendet: Mittwoch, 30. März 2011 16:46
> An: r-help@r-project.org
> Betreff: [R] VECM with UNRESTRICTED TREND
>
> Dear All,
>
> My question is:
>
> how can I estimate VECM system with "unrestricted trend" (aka
> "case 5") option as a deterministic term?
>
> As far as I know, ca.jo in urca package allows for "restricted trend"
> only [vecm
> <- ca.jo(data, type = "trace"/"eigen", ecdet = "trend", K =
> n, spec = "transitory"/"longrun")].
> Obviously, I don't have to do this in urca, so if another
> package gives the possibility, please let me know too!
>
> Thanks in advance!
>
> Greg
>
> [[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.
>
*
Confidentiality Note: The information contained in this ...{{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] Quick recode of -999 to NA in R
On Mar 30, 2011, at 16:05 , Christopher Desjardins wrote:
>>
>> dat0 <- read.table('tim1.dat', na = -999)
>>
>
> Ah ... yes. I knew that but clearly didn't at the time of my question or
> script writing.
> Thanks,
> Chris
Depending on where your data came from, you could get caught by the fact that
the above is really ...na.strings="-999"... and that is not going to work if
the actual code is (say) -999.00.
Another straightforward option is dat0[dat0 == -999] <- NA
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk Priv: pda...@gmail.com
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Quick recode of -999 to NA in R
Amen. Ditto for "-999.000", "-999.00" and all of the other ones
that various (usually Fortran) programmers have used. Has the most
recent Fortran standard come around to understanding NA?
--
Clint BowmanINTERNET: cl...@ecy.wa.gov
Air Quality Modeler INTERNET: cl...@math.utah.edu
Department of Ecology VOICE: (360) 407-6815
PO Box 47600FAX:(360) 407-7534
Olympia, WA 98504-7600
USPS: PO Box 47600, Olympia, WA 98504-7600
Parcels:300 Desmond Drive, Lacey, WA 98503-1274
On Wed, 30 Mar 2011, peter dalgaard wrote:
On Mar 30, 2011, at 16:05 , Christopher Desjardins wrote:
dat0 <- read.table('tim1.dat', na = -999)
Ah ... yes. I knew that but clearly didn't at the time of my question or
script writing.
Thanks,
Chris
Depending on where your data came from, you could get caught by the fact that the above
is really ...na.strings="-999"... and that is not going to work if the actual
code is (say) -999.00.
Another straightforward option is dat0[dat0 == -999] <- NA
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Multistep forecasting with neural networks
HI! Would anybody suggest what function/package/method can I use to obtain forecasts for horizons longer than 1 period ahead? I want to do multistep forecasts with multilayer feedforward networks. I have tried several packages, but it turned out that the functions yield only one-step forecasts. To my best knowledge, *neuralne*t and *AMORE* packages don;t allow forecasting (you need to supply new input data for the forecast period and a function sim() from AMORE package gives you 1-step aheads). *nnet* has only 1 hidden layer, while in *RSNNS* I have tried predict.rsnns(), but having used mlp for training I get only binary result (e.g. a vector of "1"s). Regards, [[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] optim and optimize are not finding the right parameter
Dear all,
I have a function that predicts DV based on one predictor pred:
pred<-c(0,300,780,1560,2340,13120)
DV<-c(0,500,1000,1400,1700,1900)
## I define Function 1 that computes the predicted value based on pred
values and parameters a and b:
calc_DV_pred <- function(a,b) {
DV_pred <- rep(0,(length(pred)))
for(i in 1:length(DV_pred)){
DV_pred[i] <- a * (1- exp( (0-b)*pred[i] ))
}
return(DV_pred)
}
## I define Function 2 that computes the sum of squared deviations for
predicted DV from actual DV:
my.function <- function(param){
b<-param
a<-max(DV)
mysum <- sum((DV - calc_DV_pred(a,b))^2)
return(mysum)
}
If I test my function for parameter b =0.001, then I get a very good fit:
pred<-c(0,300,780,1560,2340,13120)
DV<-c(0,500,1000,1400,1700,1900)
test<-my.function(0.001)
(test) # I get 11,336.81
However, when I try to optimize my.function with optim - trying to
find the value of b that minimizes the output of my.function - I get a
wrong solution:
opt1 <- optim(fn=my.function,par=c(b=0.0001),
method="L-BFGS-B", lower=0,upper=1)
(opt1)
test2<-my.function(opt1$par) # is much larger than the value of "test" above.
# And getting the same - wrong - result with optimize:
opt2 <- optimize(f=my.function,interval=c(0,0.1))
test3<-my.function(opt2$minimum)
What am I doing wrong? Thanks a lot for your recomendations!
--
Dimitri Liakhovitski
Ninah Consulting
www.ninah.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] Quick recode of -999 to NA in R
On Wed, Mar 30, 2011 at 11:51 AM, peter dalgaard wrote:
>
> On Mar 30, 2011, at 16:05 , Christopher Desjardins wrote:
>
>>>
>>> dat0 <- read.table('tim1.dat', na = -999)
>>>
>>
>> Ah ... yes. I knew that but clearly didn't at the time of my question or
>> script writing.
>> Thanks,
>> Chris
>
> Depending on where your data came from, you could get caught by the fact that
> the above is really ...na.strings="-999"... and that is not going to work if
> the actual code is (say) -999.00.
>
> Another straightforward option is dat0[dat0 == -999] <- NA
That's a good point about -999.00. If we knew all the variations in
the input then a workaround for that would be to cover all of them:
read.table("myfile", ...whatever..., na.strings = c("-999",
"-999.0", "-999.00"))
--
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] summing values by week - based on daily dates - but with somedates missing
Hi,
How about something like:
sum.by.day <- function(ff) {
by.day <- split(ff$value,weekdays(ff$dates))
lapply(by.day,sum)
}
by.grp <- split(myframe,myframe$group)
lapply(by.grp,sum.by.day)
Martyn
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On Behalf Of Dimitri Liakhovitski
Sent: 30 March 2011 15:23
To: r-help
Subject: [R] summing values by week - based on daily dates - but with
somedates missing
Dear everybody,
I have the following challenge. I have a data set with 2 subgroups,
dates (days), and corresponding values (see example code below).
Within each subgroup: I need to aggregate (sum) the values by week -
for weeks that start on a Monday (for example, 2008-12-29 was a
Monday).
I find it difficult because I have missing dates in my data - so that
sometimes I don't even have the date for some Mondays. So, I can't
write a proper loop.
I want my output to look something like this:
group dates value
group.1 2008-12-29 3.0937
group.1 2009-01-05 3.8833
group.1 2009-01-12 1.362
...
group.2 2008-12-29 2.250
group.2 2009-01-05 1.4057
group.2 2009-01-12 3.4411
...
Thanks a lot for your suggestions! The code is below:
Dimitri
### Creating example data set:
mydates<-rep(seq(as.Date("2008-12-29"), length = 43, by = "day"),2)
myfactor<-c(rep("group.1",43),rep("group.2",43))
set.seed(123)
myvalues<-runif(86,0,1)
myframe<-data.frame(dates=mydates,group=myfactor,value=myvalues)
(myframe)
dim(myframe)
## Removing same rows (dates) unsystematically:
set.seed(123)
removed.group1<-sample(1:43,size=11,replace=F)
set.seed(456)
removed.group2<-sample(44:86,size=11,replace=F)
to.remove<-c(removed.group1,removed.group2);length(to.remove)
to.remove<-to.remove[order(to.remove)]
myframe<-myframe[-to.remove,]
(myframe)
--
Dimitri Liakhovitski
Ninah Consulting
www.ninah.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.
This e-mail has been scanned for all viruses by Star.\ _...{{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.
Re: [R] Dirichlet surface
Dear David,
I think that is a small bug too, maybe because the function is constant?
is there a nice way to put the c(0,2.1) argument optionally, only if all
the parameters are 1?
Should I post the problem somewhere else (developers maybe?)
thanks:
Daniel
2011-03-30 04:42 keltezéssel, David Winsemius írta:
On Mar 29, 2011, at 4:45 PM, Kehl Dániel wrote:
Dear list members,
I want to draw surfaces of Dirichlet distributions with different
parameter settings.
My code is the following:
#
a1 <- a2 <- a3 <- 2
#a2 <- .5
#a3 <- .5
x1 <- x2 <- seq(0.01, .99, by=.01)
f <- function(x1, x2){
term1 <- gamma(a1+a2+a3)/(gamma(a1)*gamma(a2)*gamma(a3))
term2 <- x1^(a1-1)*x2^(a2-1)*(1-x1-x2)^(a3-1)
term3 <- (x1 + x2 < 1)
term1*term2*term3
}
z <- outer(x1, x2, f)
z[z<=0] <- NA
persp(x1, x2, z,
main = "Dirichlet Distribution",
col = "lightblue",
theta = 50,
phi = 20,
r = 50,
d = 0.1,
expand = 0.5,
ltheta = 90,
lphi = 180,
shade = 0.75,
ticktype = "detailed",
nticks = 5)
#
It works fine (I guess), except for a1=a2=a3=1. In that case I get
the error: Error in persp.default... invalid 'z' limits.
The z matrix has only elements 2 and NA.
Might be a buglet in persp. If you set the zlim argument to c(0,2.1),
you get the expected constant plane at z=2 for those arguments.
Any ideas are appreciated.
Thank you:
Daniel
University of Pécs
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] summing values by week - based on daily dates - but with somedates missing
Thank you, Martyn.
But it looks like this way we are getting sums by day - i.e., across
all Mondays, all Tuesdays, etc.
Maybe I did not explain well, sorry! The desired output would contain
sums for each WHOLE week - across all days that comprise that week -
Monday through Sunday.
Makes sense?
Dimitri
On Wed, Mar 30, 2011 at 12:53 PM, Martyn Byng wrote:
> Hi,
>
> How about something like:
>
> sum.by.day <- function(ff) {
> by.day <- split(ff$value,weekdays(ff$dates))
> lapply(by.day,sum)
> }
>
> by.grp <- split(myframe,myframe$group)
>
> lapply(by.grp,sum.by.day)
>
>
> Martyn
>
> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
> On Behalf Of Dimitri Liakhovitski
> Sent: 30 March 2011 15:23
> To: r-help
> Subject: [R] summing values by week - based on daily dates - but with
> somedates missing
>
> Dear everybody,
>
> I have the following challenge. I have a data set with 2 subgroups,
> dates (days), and corresponding values (see example code below).
> Within each subgroup: I need to aggregate (sum) the values by week -
> for weeks that start on a Monday (for example, 2008-12-29 was a
> Monday).
> I find it difficult because I have missing dates in my data - so that
> sometimes I don't even have the date for some Mondays. So, I can't
> write a proper loop.
> I want my output to look something like this:
> group dates value
> group.1 2008-12-29 3.0937
> group.1 2009-01-05 3.8833
> group.1 2009-01-12 1.362
> ...
> group.2 2008-12-29 2.250
> group.2 2009-01-05 1.4057
> group.2 2009-01-12 3.4411
> ...
>
> Thanks a lot for your suggestions! The code is below:
> Dimitri
>
> ### Creating example data set:
> mydates<-rep(seq(as.Date("2008-12-29"), length = 43, by = "day"),2)
> myfactor<-c(rep("group.1",43),rep("group.2",43))
> set.seed(123)
> myvalues<-runif(86,0,1)
> myframe<-data.frame(dates=mydates,group=myfactor,value=myvalues)
> (myframe)
> dim(myframe)
>
> ## Removing same rows (dates) unsystematically:
> set.seed(123)
> removed.group1<-sample(1:43,size=11,replace=F)
> set.seed(456)
> removed.group2<-sample(44:86,size=11,replace=F)
> to.remove<-c(removed.group1,removed.group2);length(to.remove)
> to.remove<-to.remove[order(to.remove)]
> myframe<-myframe[-to.remove,]
> (myframe)
>
>
>
> --
> Dimitri Liakhovitski
> Ninah Consulting
> www.ninah.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.
>
>
> This e-mail has been scanned for all viruses by Star.
>
>
>
> The Numerical Algorithms Group Ltd is a company registered in England
> and Wales with company number 1249803. The registered office is:
> Wilkinson House, Jordan Hill Road, Oxford OX2 8DR, United Kingdom.
>
> This e-mail has been scanned for all viruses by Star. The service is
> powered by MessageLabs.
>
>
--
Dimitri Liakhovitski
Ninah Consulting
www.ninah.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] optim and optimize are not finding the right parameter
Not sure it's the case here, but numeric optimizers are well-known to
be subject to scaling issues.
-- Bert
On Wed, Mar 30, 2011 at 9:45 AM, Dimitri Liakhovitski
wrote:
> Dear all,
>
> I have a function that predicts DV based on one predictor pred:
>
> pred<-c(0,300,780,1560,2340,13120)
> DV<-c(0,500,1000,1400,1700,1900)
>
> ## I define Function 1 that computes the predicted value based on pred
> values and parameters a and b:
> calc_DV_pred <- function(a,b) {
> DV_pred <- rep(0,(length(pred)))
> for(i in 1:length(DV_pred)){
> DV_pred[i] <- a * (1- exp( (0-b)*pred[i] ))
> }
> return(DV_pred)
> }
>
> ## I define Function 2 that computes the sum of squared deviations for
> predicted DV from actual DV:
> my.function <- function(param){
> b<-param
> a<-max(DV)
> mysum <- sum((DV - calc_DV_pred(a,b))^2)
> return(mysum)
> }
>
> If I test my function for parameter b =0.001, then I get a very good fit:
> pred<-c(0,300,780,1560,2340,13120)
> DV<-c(0,500,1000,1400,1700,1900)
> test<-my.function(0.001)
> (test) # I get 11,336.81
>
> However, when I try to optimize my.function with optim - trying to
> find the value of b that minimizes the output of my.function - I get a
> wrong solution:
> opt1 <- optim(fn=my.function,par=c(b=0.0001),
> method="L-BFGS-B", lower=0,upper=1)
> (opt1)
> test2<-my.function(opt1$par) # is much larger than the value of "test" above.
>
> # And getting the same - wrong - result with optimize:
> opt2 <- optimize(f=my.function,interval=c(0,0.1))
> test3<-my.function(opt2$minimum)
>
> What am I doing wrong? Thanks a lot for your recomendations!
>
> --
> Dimitri Liakhovitski
> Ninah Consulting
> www.ninah.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.
>
--
"Men by nature long to get on to the ultimate truths, and will often
be impatient with elementary studies or fight shy of them. If it were
possible to reach the ultimate truths without the elementary studies
usually prefixed to them, these would not be preparatory studies but
superfluous diversions."
-- Maimonides (1135-1204)
Bert Gunter
Genentech Nonclinical Biostatistics
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Quick recode of -999 to NA in R
On Wed, Mar 30, 2011 at 10:51 AM, peter dalgaard wrote:
>
> On Mar 30, 2011, at 16:05 , Christopher Desjardins wrote:
>
> >>
> >> dat0 <- read.table('tim1.dat', na = -999)
> >>
> >
> > Ah ... yes. I knew that but clearly didn't at the time of my question or
> > script writing.
> > Thanks,
> > Chris
>
> Depending on where your data came from, you could get caught by the fact
> that the above is really ...na.strings="-999"... and that is not going to
> work if the actual code is (say) -999.00.
>
> Another straightforward option is dat0[dat0 == -999] <- NA
>
>
Thanks Peter. That way does work well too and it does what I need.
Chris
>
> --
> Peter Dalgaard
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: pd@cbs.dk Priv: pda...@gmail.com
>
>
[[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] Using xlevels
I see the logic now. I think that more sentences in the document would
be very helpful, however. What is written is very subtle.
I suggest the following small expansion for model.matrix.Rd:
\item{data}{a data frame. If the object has a \code{terms} attribute
then it is assumed to be the result of a call to \code{model.frame},
otherwise \code{model.frame} will be called first.}
I often forget that model.frames are not a class, but an "implied"
class based on the presence of a terms component. Many users, I
suspect, do not even have this starting knowledge.
Off to make changes to model.frame.coxph and model.matrix.coxph...
Thanks for the feeback.
Terry
On Wed, 2011-03-30 at 16:36 +0100, Prof Brian Ripley wrote:
> On Wed, 30 Mar 2011, Terry Therneau wrote:
>
> > I'm working on predict.survreg and am confused about xlevels.
> > The model.frame method has the argument, but none of the standard
> > methods (model.frame.lm, model.frame.glm) appear to make use of it.
>
> But I see this in predict.lm:
>
> m <- model.frame(Terms, newdata, na.action = na.action,
> xlev = object$xlevels)
>
> It is used to remap levels in newdata to those used in the fit.
>
> >
> > The documentation for model.matrix states:
> > xlev: to be used as argument of model.frame if data has no "terms"
> > attribute.
>
> Well, the code says
>
> if (is.null(attr(data, "terms")))
> data <- model.frame(object, data, xlev=xlev)
>
> > But the terms attribute has no xlevels information in it, so I find this
> > statement completely confusing. Any insight is appreciated.
>
> It means exactly what it says: a 'data' argument with a terms
> attribute is considered to be a model frame.
>
> >
> > Terry Therneau
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Quick recode of -999 to NA in R
On Mar 30, 2011, at 11:52 AM, Gabor Grothendieck wrote:
> On Wed, Mar 30, 2011 at 11:51 AM, peter dalgaard wrote:
>>
>> On Mar 30, 2011, at 16:05 , Christopher Desjardins wrote:
>>
dat0 <- read.table('tim1.dat', na = -999)
>>>
>>> Ah ... yes. I knew that but clearly didn't at the time of my question or
>>> script writing.
>>> Thanks,
>>> Chris
>>
>> Depending on where your data came from, you could get caught by the fact
>> that the above is really ...na.strings="-999"... and that is not going to
>> work if the actual code is (say) -999.00.
>>
>> Another straightforward option is dat0[dat0 == -999] <- NA
>
> That's a good point about -999.00. If we knew all the variations in
> the input then a workaround for that would be to cover all of them:
>
> read.table("myfile", ...whatever..., na.strings = c("-999",
> "-999.0", "-999.00"))
There is also the seminal:
is.na(dat0) <- dat0 == -999
See ?is.na, specifically the assignment variant.
Regards,
Marc Schwartz
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Dirichlet surface
On Wed, 2011-03-30 at 09:55 -0700, Kehl Dániel wrote:
> Dear David,
>
> I think that is a small bug too, maybe because the function is constant?
> is there a nice way to put the c(0,2.1) argument optionally, only if all
> the parameters are 1?
> Should I post the problem somewhere else (developers maybe?)
I don't think this is a bug really; the code is having to compute limits
of the z axis and you supplied it one bit of information: a 2. If your
data are so degenerate then it is not unreasonable to expect some user
intervention. Admittedly, persp doesn't seem to work like other R
plotting functions.
You could do something like:
persp(x1, x2, z,
zlim = if(length(na.omit(unique(as.vector(z < 2){ c(0,2.1) }
else { NULL})
in your call to persp so it only uses user-defined limits if the number
of numeric values in z is less than 2.
HTH
G
> thanks:
> Daniel
>
> 2011-03-30 04:42 keltezéssel, David Winsemius írta:
> >
> > On Mar 29, 2011, at 4:45 PM, Kehl Dániel wrote:
> >
> >> Dear list members,
> >>
> >> I want to draw surfaces of Dirichlet distributions with different
> >> parameter settings.
> >> My code is the following:
> >> #
> >> a1 <- a2 <- a3 <- 2
> >> #a2 <- .5
> >> #a3 <- .5
> >> x1 <- x2 <- seq(0.01, .99, by=.01)
> >>
> >> f <- function(x1, x2){
> >> term1 <- gamma(a1+a2+a3)/(gamma(a1)*gamma(a2)*gamma(a3))
> >> term2 <- x1^(a1-1)*x2^(a2-1)*(1-x1-x2)^(a3-1)
> >> term3 <- (x1 + x2 < 1)
> >> term1*term2*term3
> >> }
> >>
> >> z <- outer(x1, x2, f)
> >> z[z<=0] <- NA
> >>
> >> persp(x1, x2, z,
> >> main = "Dirichlet Distribution",
> >> col = "lightblue",
> >> theta = 50,
> >> phi = 20,
> >> r = 50,
> >> d = 0.1,
> >> expand = 0.5,
> >> ltheta = 90,
> >> lphi = 180,
> >> shade = 0.75,
> >> ticktype = "detailed",
> >> nticks = 5)
> >> #
> >>
> >> It works fine (I guess), except for a1=a2=a3=1. In that case I get
> >> the error: Error in persp.default... invalid 'z' limits.
> >> The z matrix has only elements 2 and NA.
> >
> > Might be a buglet in persp. If you set the zlim argument to c(0,2.1),
> > you get the expected constant plane at z=2 for those arguments.
> >>
> >> Any ideas are appreciated.
> >>
> >> Thank you:
> >> Daniel
> >> University of Pécs
> >>
> >> __
> >> R-help@r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
> >
> > David Winsemius, MD
> > Heritage Laboratories
> > West Hartford, CT
> >
> >
> >
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Dr. Gavin Simpson [t] +44 (0)20 7679 0522
ECRC, UCL Geography, [f] +44 (0)20 7679 0565
Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/
UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Dirichlet surface
It helped a lot indeed, thank you very much!
Now I understand why it was a problem for persp!
Daniel
2011-03-30 10:31 keltezéssel, Gavin Simpson írta:
On Wed, 2011-03-30 at 09:55 -0700, Kehl Dániel wrote:
Dear David,
I think that is a small bug too, maybe because the function is constant?
is there a nice way to put the c(0,2.1) argument optionally, only if all
the parameters are 1?
Should I post the problem somewhere else (developers maybe?)
I don't think this is a bug really; the code is having to compute limits
of the z axis and you supplied it one bit of information: a 2. If your
data are so degenerate then it is not unreasonable to expect some user
intervention. Admittedly, persp doesn't seem to work like other R
plotting functions.
You could do something like:
persp(x1, x2, z,
zlim = if(length(na.omit(unique(as.vector(z< 2){ c(0,2.1) }
else { NULL})
in your call to persp so it only uses user-defined limits if the number
of numeric values in z is less than 2.
HTH
G
thanks:
Daniel
2011-03-30 04:42 keltezéssel, David Winsemius írta:
On Mar 29, 2011, at 4:45 PM, Kehl Dániel wrote:
Dear list members,
I want to draw surfaces of Dirichlet distributions with different
parameter settings.
My code is the following:
#
a1<- a2<- a3<- 2
#a2<- .5
#a3<- .5
x1<- x2<- seq(0.01, .99, by=.01)
f<- function(x1, x2){
term1<- gamma(a1+a2+a3)/(gamma(a1)*gamma(a2)*gamma(a3))
term2<- x1^(a1-1)*x2^(a2-1)*(1-x1-x2)^(a3-1)
term3<- (x1 + x2< 1)
term1*term2*term3
}
z<- outer(x1, x2, f)
z[z<=0]<- NA
persp(x1, x2, z,
main = "Dirichlet Distribution",
col = "lightblue",
theta = 50,
phi = 20,
r = 50,
d = 0.1,
expand = 0.5,
ltheta = 90,
lphi = 180,
shade = 0.75,
ticktype = "detailed",
nticks = 5)
#
It works fine (I guess), except for a1=a2=a3=1. In that case I get
the error: Error in persp.default... invalid 'z' limits.
The z matrix has only elements 2 and NA.
Might be a buglet in persp. If you set the zlim argument to c(0,2.1),
you get the expected constant plane at z=2 for those arguments.
Any ideas are appreciated.
Thank you:
Daniel
University of Pécs
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] optim and optimize are not finding the right parameter
Try this:
pred <- pred/1e06
DV <- DV/1e03
opt1 <- optim(fn=my.function, par=1.0)
opt2 <- optim(fn=my.function, par=1.0, method="BFGS")
opt3 <- optim(fn=my.function, par=1.0, method="L-BFGS-B", lower=0, upper=1)
opt1
opt2
opt3
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: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Dimitri Liakhovitski
Sent: Wednesday, March 30, 2011 12:45 PM
To: r-help
Subject: [R] optim and optimize are not finding the right parameter
Dear all,
I have a function that predicts DV based on one predictor pred:
pred<-c(0,300,780,1560,2340,13120)
DV<-c(0,500,1000,1400,1700,1900)
## I define Function 1 that computes the predicted value based on pred
values and parameters a and b:
calc_DV_pred <- function(a,b) {
DV_pred <- rep(0,(length(pred)))
for(i in 1:length(DV_pred)){
DV_pred[i] <- a * (1- exp( (0-b)*pred[i] ))
}
return(DV_pred)
}
## I define Function 2 that computes the sum of squared deviations for
predicted DV from actual DV:
my.function <- function(param){
b<-param
a<-max(DV)
mysum <- sum((DV - calc_DV_pred(a,b))^2)
return(mysum)
}
If I test my function for parameter b =0.001, then I get a very good
fit:
pred<-c(0,300,780,1560,2340,13120)
DV<-c(0,500,1000,1400,1700,1900)
test<-my.function(0.001)
(test) # I get 11,336.81
However, when I try to optimize my.function with optim - trying to
find the value of b that minimizes the output of my.function - I get a
wrong solution:
opt1 <- optim(fn=my.function,par=c(b=0.0001),
method="L-BFGS-B", lower=0,upper=1)
(opt1)
test2<-my.function(opt1$par) # is much larger than the value of "test"
above.
# And getting the same - wrong - result with optimize:
opt2 <- optimize(f=my.function,interval=c(0,0.1))
test3<-my.function(opt2$minimum)
What am I doing wrong? Thanks a lot for your recomendations!
--
Dimitri Liakhovitski
Ninah Consulting
www.ninah.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] Dirichlet surface
Actually, it works for the a=1 case, not for the others. It still gives
the invalid 'zlim argument' error.
I'll try to work it out maybe instead of NULL giving a c which is
dependent on the max(z).
Daniel
2011-03-30 10:42 keltezéssel, Kehl Dániel írta:
It helped a lot indeed, thank you very much!
Now I understand why it was a problem for persp!
Daniel
2011-03-30 10:31 keltezéssel, Gavin Simpson írta:
On Wed, 2011-03-30 at 09:55 -0700, Kehl Dániel wrote:
Dear David,
I think that is a small bug too, maybe because the function is
constant?
is there a nice way to put the c(0,2.1) argument optionally, only if
all
the parameters are 1?
Should I post the problem somewhere else (developers maybe?)
I don't think this is a bug really; the code is having to compute limits
of the z axis and you supplied it one bit of information: a 2. If your
data are so degenerate then it is not unreasonable to expect some user
intervention. Admittedly, persp doesn't seem to work like other R
plotting functions.
You could do something like:
persp(x1, x2, z,
zlim = if(length(na.omit(unique(as.vector(z< 2){ c(0,2.1) }
else { NULL})
in your call to persp so it only uses user-defined limits if the number
of numeric values in z is less than 2.
HTH
G
thanks:
Daniel
2011-03-30 04:42 keltezéssel, David Winsemius írta:
On Mar 29, 2011, at 4:45 PM, Kehl Dániel wrote:
Dear list members,
I want to draw surfaces of Dirichlet distributions with different
parameter settings.
My code is the following:
#
a1<- a2<- a3<- 2
#a2<- .5
#a3<- .5
x1<- x2<- seq(0.01, .99, by=.01)
f<- function(x1, x2){
term1<- gamma(a1+a2+a3)/(gamma(a1)*gamma(a2)*gamma(a3))
term2<- x1^(a1-1)*x2^(a2-1)*(1-x1-x2)^(a3-1)
term3<- (x1 + x2< 1)
term1*term2*term3
}
z<- outer(x1, x2, f)
z[z<=0]<- NA
persp(x1, x2, z,
main = "Dirichlet Distribution",
col = "lightblue",
theta = 50,
phi = 20,
r = 50,
d = 0.1,
expand = 0.5,
ltheta = 90,
lphi = 180,
shade = 0.75,
ticktype = "detailed",
nticks = 5)
#
It works fine (I guess), except for a1=a2=a3=1. In that case I get
the error: Error in persp.default... invalid 'z' limits.
The z matrix has only elements 2 and NA.
Might be a buglet in persp. If you set the zlim argument to c(0,2.1),
you get the expected constant plane at z=2 for those arguments.
Any ideas are appreciated.
Thank you:
Daniel
University of Pécs
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using xlevels
Terry,
The fact that model.frame attaches xlevels to
the terms based on factors in the input data.frame
(and attaches dataClass based on the input data.frame),
but the subsequent call to model.matrix is responsible
for turning character vectors in the data.frame into
factors (and then into contrasts) is part of the reason
that you cannot use predict() on an lmObject created
using a data.frame with character vectors in it.
> d <- data.frame(y=1:10, x=rep(LETTERS[1:3],c(3,3,4)),
stringsAsFactors=FALSE)
> fit <- lm(data=d, y~x)
Warning message:
In model.matrix.default(mt, mf, contrasts) :
variable 'x' converted to a factor
> predict(fit, newdata=data.frame(x=c("A","C"))) # expect c(2.0, 8.5)
Error: variable 'x' was fitted with type "other" but type "factor" was
supplied
This is one way that changing the default stringsAsFactors=TRUE
can cause problems.
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
> -Original Message-
> From: r-help-boun...@r-project.org
> [mailto:r-help-boun...@r-project.org] On Behalf Of Terry Therneau
> Sent: Wednesday, March 30, 2011 10:29 AM
> To: Prof Brian Ripley
> Cc: r-help@r-project.org
> Subject: Re: [R] Using xlevels
>
> I see the logic now. I think that more sentences in the
> document would
> be very helpful, however. What is written is very subtle.
> I suggest the following small expansion for model.matrix.Rd:
>
> \item{data}{a data frame. If the object has a \code{terms}
> attribute
> then it is assumed to be the result of a call to \code{model.frame},
> otherwise \code{model.frame} will be called first.}
>
> I often forget that model.frames are not a class, but an "implied"
> class based on the presence of a terms component. Many users, I
> suspect, do not even have this starting knowledge.
>
> Off to make changes to model.frame.coxph and model.matrix.coxph...
>
> Thanks for the feeback.
>
> Terry
>
>
> On Wed, 2011-03-30 at 16:36 +0100, Prof Brian Ripley wrote:
> > On Wed, 30 Mar 2011, Terry Therneau wrote:
> >
> > > I'm working on predict.survreg and am confused about xlevels.
> > > The model.frame method has the argument, but none of the standard
> > > methods (model.frame.lm, model.frame.glm) appear to make
> use of it.
> >
> > But I see this in predict.lm:
> >
> > m <- model.frame(Terms, newdata, na.action = na.action,
> > xlev = object$xlevels)
> >
> > It is used to remap levels in newdata to those used in the fit.
> >
> > >
> > > The documentation for model.matrix states:
> > > xlev: to be used as argument of model.frame if data has
> no "terms"
> > > attribute.
> >
> > Well, the code says
> >
> > if (is.null(attr(data, "terms")))
> > data <- model.frame(object, data, xlev=xlev)
> >
> > > But the terms attribute has no xlevels information in it,
> so I find this
> > > statement completely confusing. Any insight is appreciated.
> >
> > It means exactly what it says: a 'data' argument with a terms
> > attribute is considered to be a model frame.
> >
> > >
> > > Terry Therneau
> > >
> > > __
> > > R-help@r-project.org mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
> > >
> >
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] fonts in mosaic
Achim
I simply want to replace the font R uses on mosaic (whatever it is)
by a font of my choice (say Calibri or Arial)
because I need to embed the R charts in a PowerPoint
presentation and want the fonts to match.
And I want the most simple way of accomplishing this.
I worked my way through the strucplot vignette,
but I could not extract enough information there.
Is there some information about the proper font names
to use in R?
> Personally, I simply change the size of the device I'm plotting on. When
> I plot on a large device, the fonts will be relatively smaller, and vice
> > versa. This is what I do when including graphics in PDF files (papers,
> > slides, reports, etc.).
> >
> > For fine control, you can set the arguments of the labeling function
> > employed. ?strucplot shows that the default is ?labeling_border which
> > has several arguments. For example you can set the graphical parameters
> > of the labels (gp_labels) or the graphical parameters of the variable
> > names (gp_varnames). Both arguments take ?gpar lists ("grid" graphical
> > parameters). For example you may do
> >
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] optim and optimize are not finding the right parameter
Thank you, Ravi - definitely better!
On Wed, Mar 30, 2011 at 2:06 PM, Ravi Varadhan wrote:
> Try this:
>
> pred <- pred/1e06
>
> DV <- DV/1e03
>
> opt1 <- optim(fn=my.function, par=1.0)
> opt2 <- optim(fn=my.function, par=1.0, method="BFGS")
> opt3 <- optim(fn=my.function, par=1.0, method="L-BFGS-B", lower=0, upper=1)
> opt1
> opt2
> opt3
>
> 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: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
> Behalf Of Dimitri Liakhovitski
> Sent: Wednesday, March 30, 2011 12:45 PM
> To: r-help
> Subject: [R] optim and optimize are not finding the right parameter
>
> Dear all,
>
> I have a function that predicts DV based on one predictor pred:
>
> pred<-c(0,300,780,1560,2340,13120)
> DV<-c(0,500,1000,1400,1700,1900)
>
> ## I define Function 1 that computes the predicted value based on pred
> values and parameters a and b:
> calc_DV_pred <- function(a,b) {
> DV_pred <- rep(0,(length(pred)))
> for(i in 1:length(DV_pred)){
> DV_pred[i] <- a * (1- exp( (0-b)*pred[i] ))
> }
> return(DV_pred)
> }
>
> ## I define Function 2 that computes the sum of squared deviations for
> predicted DV from actual DV:
> my.function <- function(param){
> b<-param
> a<-max(DV)
> mysum <- sum((DV - calc_DV_pred(a,b))^2)
> return(mysum)
> }
>
> If I test my function for parameter b =0.001, then I get a very good
> fit:
> pred<-c(0,300,780,1560,2340,13120)
> DV<-c(0,500,1000,1400,1700,1900)
> test<-my.function(0.001)
> (test) # I get 11,336.81
>
> However, when I try to optimize my.function with optim - trying to
> find the value of b that minimizes the output of my.function - I get a
> wrong solution:
> opt1 <- optim(fn=my.function,par=c(b=0.0001),
> method="L-BFGS-B", lower=0,upper=1)
> (opt1)
> test2<-my.function(opt1$par) # is much larger than the value of "test"
> above.
>
> # And getting the same - wrong - result with optimize:
> opt2 <- optimize(f=my.function,interval=c(0,0.1))
> test3<-my.function(opt2$minimum)
>
> What am I doing wrong? Thanks a lot for your recomendations!
>
> --
> Dimitri Liakhovitski
> Ninah Consulting
> www.ninah.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.
>
>
--
Dimitri Liakhovitski
Ninah Consulting
www.ninah.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] fonts in mosaic
Try this:
windowsFonts(calibri = windowsFont("Calibri"))
mosaic(UCBAdmissions, labeling_args = list(
gp_labels = gpar(fontsize = 12, fontfamily = "calibri"),
gp_varnames = gpar(fontsize = 16, fontfamily = "calibri")
))
On Wed, Mar 30, 2011 at 3:25 PM, Erich Neuwirth
wrote:
> Achim
> I simply want to replace the font R uses on mosaic (whatever it is)
> by a font of my choice (say Calibri or Arial)
> because I need to embed the R charts in a PowerPoint
> presentation and want the fonts to match.
> And I want the most simple way of accomplishing this.
> I worked my way through the strucplot vignette,
> but I could not extract enough information there.
> Is there some information about the proper font names
> to use in R?
>
>
>> Personally, I simply change the size of the device I'm plotting on. When
>> I plot on a large device, the fonts will be relatively smaller, and vice
>> > versa. This is what I do when including graphics in PDF files (papers,
>> > slides, reports, etc.).
>> >
>> > For fine control, you can set the arguments of the labeling function
>> > employed. ?strucplot shows that the default is ?labeling_border which
>> > has several arguments. For example you can set the graphical parameters
>> > of the labels (gp_labels) or the graphical parameters of the variable
>> > names (gp_varnames). Both arguments take ?gpar lists ("grid" graphical
>> > parameters). For example you may do
>> >
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] save ordinary numerical calculations as pdf
Hi,
If you want grid graphics:
For data.frames and matrices, gridExtra has a grid.table() function.
For strings (paragraph), Rgraphics has a function too, whose name i
forget. It could be possible to combine the two and define a method to
display lists as well.
HTH,
baptiste
On 30 March 2011 22:51, Maas James Dr (MED) wrote:
> I'd like to save some calculation outputs as a pdf, to incorporate with
> others in a document. I've tried
>
> pdf("filename")
> name_of_object_to_output
> dev.off()
>
> but it doesn't seem to work, appears that this pdf function is for graphics?
>
> Is there a way to output numerical objects to pdf?
>
> Thanks
>
> J
>
> Dr. Jim Maas
> University of East Anglia
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Dirichlet surface
On Wed, 2011-03-30 at 11:12 -0700, Kehl Dániel wrote:
> Actually, it works for the a=1 case, not for the others. It still gives
> the invalid 'zlim argument' error.
> I'll try to work it out maybe instead of NULL giving a c which is
> dependent on the max(z).
Sorry, I misread the helpfile - the other lims can be NULL but not zlim.
This works:
persp(x1, x2, z,
zlim = if(length(na.omit(unique(as.vector(z < 2) {
c(0,2.1)
} else {
range(z, na.rm = TRUE)
})
HTH
G
> Daniel
>
> 2011-03-30 10:42 keltezéssel, Kehl Dániel írta:
> > It helped a lot indeed, thank you very much!
> > Now I understand why it was a problem for persp!
> >
> > Daniel
> >
> > 2011-03-30 10:31 keltezéssel, Gavin Simpson írta:
> >> On Wed, 2011-03-30 at 09:55 -0700, Kehl Dániel wrote:
> >>> Dear David,
> >>>
> >>> I think that is a small bug too, maybe because the function is
> >>> constant?
> >>> is there a nice way to put the c(0,2.1) argument optionally, only if
> >>> all
> >>> the parameters are 1?
> >>> Should I post the problem somewhere else (developers maybe?)
> >> I don't think this is a bug really; the code is having to compute limits
> >> of the z axis and you supplied it one bit of information: a 2. If your
> >> data are so degenerate then it is not unreasonable to expect some user
> >> intervention. Admittedly, persp doesn't seem to work like other R
> >> plotting functions.
> >>
> >> You could do something like:
> >>
> >> persp(x1, x2, z,
> >>zlim = if(length(na.omit(unique(as.vector(z< 2){ c(0,2.1) }
> >> else { NULL})
> >>
> >> in your call to persp so it only uses user-defined limits if the number
> >> of numeric values in z is less than 2.
> >>
> >> HTH
> >>
> >> G
> >>
> >>> thanks:
> >>> Daniel
> >>>
> >>> 2011-03-30 04:42 keltezéssel, David Winsemius írta:
> On Mar 29, 2011, at 4:45 PM, Kehl Dániel wrote:
>
> > Dear list members,
> >
> > I want to draw surfaces of Dirichlet distributions with different
> > parameter settings.
> > My code is the following:
> > #
> > a1<- a2<- a3<- 2
> > #a2<- .5
> > #a3<- .5
> > x1<- x2<- seq(0.01, .99, by=.01)
> >
> > f<- function(x1, x2){
> > term1<- gamma(a1+a2+a3)/(gamma(a1)*gamma(a2)*gamma(a3))
> > term2<- x1^(a1-1)*x2^(a2-1)*(1-x1-x2)^(a3-1)
> > term3<- (x1 + x2< 1)
> > term1*term2*term3
> > }
> >
> > z<- outer(x1, x2, f)
> > z[z<=0]<- NA
> >
> > persp(x1, x2, z,
> > main = "Dirichlet Distribution",
> > col = "lightblue",
> > theta = 50,
> > phi = 20,
> > r = 50,
> > d = 0.1,
> > expand = 0.5,
> > ltheta = 90,
> > lphi = 180,
> > shade = 0.75,
> > ticktype = "detailed",
> > nticks = 5)
> > #
> >
> > It works fine (I guess), except for a1=a2=a3=1. In that case I get
> > the error: Error in persp.default... invalid 'z' limits.
> > The z matrix has only elements 2 and NA.
> Might be a buglet in persp. If you set the zlim argument to c(0,2.1),
> you get the expected constant plane at z=2 for those arguments.
> > Any ideas are appreciated.
> >
> > Thank you:
> > Daniel
> > University of Pécs
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> David Winsemius, MD
> Heritage Laboratories
> West Hartford, CT
>
>
>
> >>> __
> >>> R-help@r-project.org mailing list
> >>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>> PLEASE do read the posting guide
> >>> http://www.R-project.org/posting-guide.html
> >>> and provide commented, minimal, self-contained, reproducible code.
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Dr. Gavin Simpson [t] +44 (0)20 7679 0522
ECRC, UCL Geography, [f] +44 (0)20 7679 0565
Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/
UK. W
Re: [R] fonts in mosaic
On Wed, 30 Mar 2011, Erich Neuwirth wrote:
Achim
I simply want to replace the font R uses on mosaic (whatever it is)
by a font of my choice (say Calibri or Arial)
because I need to embed the R charts in a PowerPoint
presentation and want the fonts to match.
Ah, ok, sorry I misread your mail.
And I want the most simple way of accomplishing this.
I worked my way through the strucplot vignette,
but I could not extract enough information there.
Is there some information about the proper font names
to use in R?
Hmm, if you look at
help("gpar", package = "grid")
it allows specification of a "fontfamily" argument. But I'm not sure how
this can be utilized to specify certain TTFs. But I remember seeing posts
on the mailing list related to this.
In any case, I don't think that this is specific to "vcd" but pertains to
"grid" in general.
hth,
Z
Personally, I simply change the size of the device I'm plotting on. When
I plot on a large device, the fonts will be relatively smaller, and vice
versa. This is what I do when including graphics in PDF files (papers,
slides, reports, etc.).
For fine control, you can set the arguments of the labeling function
employed. ?strucplot shows that the default is ?labeling_border which
has several arguments. For example you can set the graphical parameters
of the labels (gp_labels) or the graphical parameters of the variable
names (gp_varnames). Both arguments take ?gpar lists ("grid" graphical
parameters). For example you may do
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] summing values by week - based on daily dates - but with some dates missing
Try this:
aggregate(value ~ group + format(dates, "%Y.%W"), myframe, FUN = sum)
On Wed, Mar 30, 2011 at 11:23 AM, Dimitri Liakhovitski
wrote:
> Dear everybody,
>
> I have the following challenge. I have a data set with 2 subgroups,
> dates (days), and corresponding values (see example code below).
> Within each subgroup: I need to aggregate (sum) the values by week -
> for weeks that start on a Monday (for example, 2008-12-29 was a
> Monday).
> I find it difficult because I have missing dates in my data - so that
> sometimes I don't even have the date for some Mondays. So, I can't
> write a proper loop.
> I want my output to look something like this:
> group dates value
> group.1 2008-12-29 3.0937
> group.1 2009-01-05 3.8833
> group.1 2009-01-12 1.362
> ...
> group.2 2008-12-29 2.250
> group.2 2009-01-05 1.4057
> group.2 2009-01-12 3.4411
> ...
>
> Thanks a lot for your suggestions! The code is below:
> Dimitri
>
> ### Creating example data set:
> mydates<-rep(seq(as.Date("2008-12-29"), length = 43, by = "day"),2)
> myfactor<-c(rep("group.1",43),rep("group.2",43))
> set.seed(123)
> myvalues<-runif(86,0,1)
> myframe<-data.frame(dates=mydates,group=myfactor,value=myvalues)
> (myframe)
> dim(myframe)
>
> ## Removing same rows (dates) unsystematically:
> set.seed(123)
> removed.group1<-sample(1:43,size=11,replace=F)
> set.seed(456)
> removed.group2<-sample(44:86,size=11,replace=F)
> to.remove<-c(removed.group1,removed.group2);length(to.remove)
> to.remove<-to.remove[order(to.remove)]
> myframe<-myframe[-to.remove,]
> (myframe)
>
>
>
> --
> Dimitri Liakhovitski
> Ninah Consulting
> www.ninah.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.
>
--
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using graphics straight from R into published articles
Wow - thanks all for your helpful replies. Awesome forum. Am I right to assume that you use the postscript function to create .ps and .pdf files from R? blanco -- View this message in context: http://r.789695.n4.nabble.com/Using-graphics-straight-from-R-into-published-articles-tp3415401p3418682.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] How to define new operator in R?
Hello, everyone! Does anyone know how make some symbols have special means in R? For example, we know that "+" in R means the sum of the two operand on its left and right. I want to define some operators in R by myself. Is this possible? Regards! -- Chuanlong Du Department of Statistcis Iowa State University Ames, IA, US 50011 [[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] How can we modify attributes of files in Windows Systems using R?
Hello! Does any one know how to modify attributes of files (hide, read only and etc) in Windows Systems using R? I tried to use shell to call system command, but it seems that it doesn't work well. Sometimes it works but sometimes not. Regards! -- Chuanlong Du Department of Statistcis Iowa State University Ames, IA, US 50011 [[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] Noble America Announcement for Summer Internship 2011
# Noble America Announcement for Summer Internship 2011 ## Organization Name: Noble America, http://www.thisisnoble.com/ Location: Stamford, Connecticut, USA Period: 10-12 Weeks for Summer 2011 starting May/June 2011 Brief Description of work: The summer intern will work on building tools for the Trading/Risk groups, and would require: 1. Strong knowledge of R and C# programming 2. Strong knowledge of time series (more generally mathematical and statistical knowledge is welcome but not strictly necessary) 3. Financial engineering knowledge, specifically in the area of derivative pricing, is a great plus. General Finance knowledge is _not_ required. Typical tasks include: 1. Helping with the testing and building of new systems, tools, models, etc for commodity trading desks 2. Documentation of existing tools, models, etc The typical intern may expect to work closely with Risk and Trading and gain some valuable experience in the makings of tools and analysis pertaining to commodities trading (especially Natural Gas and Power) == Application: Please send CV by email For information please address : Wang Yu (wan...@thisisnoble.com) and Yannis Tzamouranis (y...@thisisnoble.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] Is there a function to save the content in R console to some file?
Hello, everyone! Does anyone know whether there's a function in R which can save the content in R console to some file? I'm using Windows system and the usual R console instead of Rstudio. Regards! -- Chuanlong Du Department of Statistcis Iowa State University Ames, IA, US 50011 [[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] postscript rotation (bug?)
Dear Marc, Your answer on post https://stat.ethz.ch/pipermail/r-help/2005-March/067634.html has a broken link. Thank you, W __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] VECM with UNRESTRICTED TREND
Hello Bernhard, Thank You very much. Unfortunately I'm still not really sure how should I use dummy vars in this context... If I have a system of three variables (x, y, z), lag order = 2 and 1 cointegrating relation, what should I do? I mean, what kind of 'pattern' should be used to create those dummy variables, what should they represent and how many of them do I need? Many thanks in advance! Best, Greg 2011/3/30 Pfaff, Bernhard Dr. > Hello Greg, > > you can exploit the argument 'dumvar' for this. See ?ca.jo > > Best, > Bernhard > > > -Ursprüngliche Nachricht- > > Von: r-help-boun...@r-project.org > > [mailto:r-help-boun...@r-project.org] Im Auftrag von Grzegorz Konat > > Gesendet: Mittwoch, 30. März 2011 16:46 > > An: r-help@r-project.org > > Betreff: [R] VECM with UNRESTRICTED TREND > > > > Dear All, > > > > My question is: > > > > how can I estimate VECM system with "unrestricted trend" (aka > > "case 5") option as a deterministic term? > > > > As far as I know, ca.jo in urca package allows for "restricted trend" > > only [vecm > > <- ca.jo(data, type = "trace"/"eigen", ecdet = "trend", K = > > n, spec = "transitory"/"longrun")]. > > Obviously, I don't have to do this in urca, so if another > > package gives the possibility, please let me know too! > > > > Thanks in advance! > > > > Greg > > > > [[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. > > > * > Confidentiality Note: The information contained in this message, > and any attachments, may contain confidential and/or privileged > material. It is intended solely for the person(s) or entity to > which it is addressed. Any review, retransmission, dissemination, > or taking of any action in reliance upon this information by > persons or entities other than the intended recipient(s) is > prohibited. If you received this in error, please contact the > sender and delete the material from any computer. > * > > [[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] Class noquote to matrix
I knew there had to be a simple solution. Thank you!
On Mar 30, 2011, at 3:04 AM, Peter Ehlers wrote:
> On 2011-03-29 19:12, Mark Ebbert wrote:
>> Hi,
>>
>> I apologize if the solution is right in front of me, but I can't find
>> anything on how to convert a class of 'noquote' to 'matrix'. I've tried
>> as.matrix and I've tried coercing it by 'class(x)<-matrix', but of course
>> that didn't work. I've been using the function 'symnum' which returns an
>> object of class 'noquote'.
>>
>> Here is an example from the 'symnum' help page:
>> pval<- rev(sort(c(outer(1:6, 10^-(1:3)
>> symp<- symnum(pval, corr=FALSE,
>> cutpoints = c(0, .001,.01,.05, .1, 1),
>> symbols = c("***","**","*","."," "))
>
> One way would be to unclass the object first:
>
> as.matrix(unclass(symp))
>
> Peter Ehlers
>
>>
>> Thanks!
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to define new operator in R?
On Mar 30, 2011, at 11:04 AM, Chuanlong Du wrote: Hello, everyone! Does anyone know how make some symbols have special means in R? For example, we know that "+" in R means the sum of the two operand on its left and right. I want to define some operators in R by myself. Is this possible? You can create your own infix operators by surrounding them with matching percent signs: `%` ?"%in%" -- David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] postscript rotation (bug?)
On Mar 30, 2011, at 11:51 AM, Wayne Lee wrote: > Dear Marc, > > Your answer on post > https://stat.ethz.ch/pipermail/r-help/2005-March/067634.html > has a broken link. > > Thank you, > W Wayne, here is the original thread being referenced from March of 2005: https://stat.ethz.ch/pipermail/r-help/2005-March/067030.html Possible that some reorg of the archives happened since then resulting in the URL being lost, or the particular post that I was referencing in my reply is gone. I would focus on Prof. Ripley's reply to Roger Peng in that thread: https://stat.ethz.ch/pipermail/r-help/2005-March/067051.html HTH, Marc Schwartz __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Using graphics straight from R into published articles
On Wed, Mar 30, 2011 at 09:56:09AM -0700, blanco wrote: > Wow - thanks all for your helpful replies. Awesome forum. > > Am I right to assume that you use the postscript function to create .ps and > .pdf files from R? almost: postscript(..., onefile=FALSE) # for eps pdf() # for PDF And don't forget to close the device with dev.off() after the plot. cu Philipp -- Dr. Philipp Pagel Lehrstuhl für Genomorientierte Bioinformatik Technische Universität München Wissenschaftszentrum Weihenstephan Maximus-von-Imhof-Forum 3 85354 Freising, Germany http://webclu.bio.wzw.tum.de/~pagel/ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Creating 3 vectors that sum to 1
Here is one comparison:
library(TeachingDemos)
library(gtools)
dirfun1 <- function(n, pch='.',...,orig=TRUE) {
if(orig) {
tmp <- matrix( runif(n*2), n, 2 )
rtmp <- cbind( pmin( tmp[,1], tmp[,2] ), abs( tmp[,1]-tmp[,2]
), 1-pmax( tmp[,1], tmp[,2] ) )
triplot(rtmp, pch=pch, ...)
} else {
triplot( rdirichlet(n, rep(1,3)), pch=pch, ... )
}
}
vis.test( 1000, FUN=dirfun1 )
Now on each set of plots click (or have an independent person click) on the one
that looks most different from the others (singing/humming the song from Sesame
Street is purely optional).
I did it twice and got a p-value of 1 each time.
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111
> -Original Message-
> From: Ravi Varadhan [mailto:rvarad...@jhmi.edu]
> Sent: Tuesday, March 29, 2011 4:41 PM
> To: Greg Snow
> Cc: r-help@r-project.org
> Subject: Re: RE: [R] Creating 3 vectors that sum to 1
>
> Hi Greg,
>
> Thanks.
>
> Here is one approach to speeding up the 1st method that I had
> suggested:
>
> n <- 1
> set.seed(123)
> rtrg <- matrix(NA, n, 3)
> system.time(for (i in 1:n) rtrg[i,] <- diff(c(0, sort(runif(2)), 1)))
>
> set.seed(123)
> system.time({
> tmp <- matrix(runif(n*2), n, 2, byrow=TRUE)
> rtrg.1 <- cbind(pmin(tmp[,1], tmp[,2]), abs(tmp[,1] - tmp[,2]),1 -
> pmax(tmp[,1], tmp[,2]))
> })
>
> all.equal(rtrg, rtrg.1)
>
> Now, how can we use vis.test to test differences between these?
>
> 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: Greg Snow
> Date: Tuesday, March 29, 2011 3:42 pm
> Subject: RE: [R] Creating 3 vectors that sum to 1
> To: Ravi Varadhan
> Cc: "r-help@r-project.org"
>
>
> > Or we could expand a bit more:
> >
> > require(TeachingDemos)
> > require(gtools)
> >
> > n <- 1000
> > rtrg <- matrix(NA, n, 3)
> > for (i in 1:n) rtrg[i,] <- diff(c(0, sort(runif(2)), 1))
> >
> > rtrg2 <- matrix(NA, n, 3)
> > for (i in 1:n) {
> > tmp <- runif(3)
> > rtrg2[i, ] <- tmp/sum(tmp)
> > }
> >
> > rtrg3 <- matrix( rexp(n*3), ncol=3 )
> > rtrg3 <- rtrg3/rowSums(rtrg3)
> >
> > rtrg4 <- rdirichlet(n, rep(1,3))
> >
> > par(mfrow=c(2,2))
> > triplot(rtrg, pch='.') # Looks more uniformly distributed
> > triplot(rtrg2, col=2, pch='.') # Corners are sparsely populated
> > triplot(rtrg3, col=3, pch='.')
> > triplot(rtrg4, col=4, pch='.')
> >
> >
> >
> >
> > What could also be interesting in using vis.test (also TeachingDemos)
> > to see which can be told apart from each other. My guess is that
> > rtrg2 method will be visible different from the other 3, but the
> other
> > 3 will be indistinguishable from each other.
> >
> > The last 2 have the advantage (the 2nd could be rewritten to have the
> > same advantage) of being much quicker, not sure how to speed up the
> > 1st noticibly.
> >
> >
> >
> > --
> > Gregory (Greg) L. Snow Ph.D.
> > Statistical Data Center
> > Intermountain Healthcare
> > greg.s...@imail.org
> > 801.408.8111
> >
> >
> > > -Original Message-
> > > From: Ravi Varadhan [
> > > Sent: Tuesday, March 29, 2011 12:59 PM
> > > To: Ravi Varadhan
> > > Cc: Greg Snow; r-help@r-project.org
> > > Subject: Re: [R] Creating 3 vectors that sum to 1
> > >
> > >
> > > Here is an exploration of two different 3-tuple generators (that
> sum
> > to
> > > 1) based on Greg's triplot function:
> > >
> > > require(TeachingDemos)
> > >
> > > n <- 1000
> > > rtrg <- matrix(NA, n, 3)
> > > for (i in 1:n) rtrg[i,] <- diff(c(0, sort(runif(2)), 1))
> > >
> > > rtrg2 <- matrix(NA, n, 3)
> > > for (i in 1:n) {
> > > tmp <- runif(3)
> > > rtrg2[i, ] <- tmp/sum(tmp)
> > > }
> > >
> > > par(mfrow=c(2,1))
> > > triplot(rtrg) # Looks more uniformly distributed
> > > triplot(rtrg2, col=2) # Corners are sparsely populated
> > >
> > > 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: Ravi Varadhan
> > > Date: Tuesday, March 29, 2011 2:33 pm
> > > Subject: Re: [R] Creating 3 vectors that sum to 1
> > > To: Greg Snow
> > > Cc: "r-help@r-project.org"
> > >
> > >
> > > > The following one-liner generates uniformly distributed 3-tuples
> that
> > > > sum to 1:
> > > >
> > > > diff(c(0, sort(runif(2)), 1))
> > > >
> > > > More, generally you can generate n-tuples that sum to unity as:
> > > >
> > > > diff(c(0, sort(runif(n-1)), 1))
> > > >
> > > >
> > > > Ravi.
> > > >
> > > >
> __
Re: [R] Using graphics straight from R into published articles
On Mar 30, 2011, at 11:56 AM, blanco wrote: Wow - thanks all for your helpful replies. Awesome forum. Am I right to assume that you use the postscript function to create .ps and .pdf files from R? No, just .ps and .eps files. The pdf() functon is for the obvious purposes. ?Devices ?capabilities -- David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
