Re: [R] Plotting MDS (multidimensional scaling)
Daniel Malter daniel at umd.edu writes: Let's say I have an equilateral triangle. Then, the three Euclidean distances between points A, B, and C are all equal. That is, dist(AB)=dist(AC)=dist(BC). Let the points A, B, and C have (x,y)-coordinates (0,0), (2,0), and (1,sqrt(3)). Then, MDS should reproduce an equilateral triangle, which it does if there are only three points. require(MASS) x=c(0,2,1,0,0,sqrt(3)) dim(x)=c(3,2) d1=dist(x) fit1-isoMDS(d1) plot(fit1$points, xlab=Coordinate 1, ylab=Coordinate 2, main=Metric MDS,type=n) text(fit1$points, labels = c('A','B','C'), cex=1) So far so good, until I add more points. Now assume, I add a fourth point D at {0,2*sqrt(3)}. This produces the rectangular triangle ABD with hypothenuse BD that encompasses the smaller triangle ABC such that C lies in the middle between B and D. Then, MDS should reproduce the rectangular triangle ABD and the equilateral triangle ABC within it. However, even though distance matrix d2 below still indicates that ABC is an equilateral triangle, the plot of the MDS does not confirm this. x=c(0,2,1,0,0,0,sqrt(3),2*sqrt(3)) dim(x)=c(4,2) d2=dist(x) fit2-isoMDS(d2) plot(fit2$points, xlab=Coordinate 1, ylab=Coordinate 2, main=Metric MDS,type=n) text(fit2$points, labels = c('A','B','C','D'), cex=1) Daniel, Mario Valle already told you about asp=1 in plot() to force equal aspect ratio, and MASS also has eqscplot() function for plots with geometrically equal scaling. However, your example above hints that there is something else you should take care of: You label your plot as Metric MDS, but isoMDS does not do metric MDS. The title in its documentation reads Kruskal's Non-metric Multidimensional Scaling. In this case you happened to have metric MDS, because isoMDS uses metric scaling as its default starting configuration, and in this case that starting configuration is a perfect fit (stress = 0), and isoMDS() makes no iterations to change the starting configuration. If you want to work with metric MDS, use cmdscale() which does metric MDS. Cheers, jari Oksanen The reason for this is that the dimension of the plot is automatically scaled to fit the points. This distorts the visual impression of the distances, angular relationships, and relative locations. If you plot on a square pane, however, peace and order are restored in the galaxy. plot(fit2$points, xlab=Coordinate 1, ylab=Coordinate 2, main=Metric MDS,type=n,xlim=c(-3,3),ylim=c(-3,3)) text(fit2$points, labels = c('A','B','C','D'), cex=1) Best, Daniel -- View this message in context: http://r.789695.n4.nabble.com/Plotting-MDS-multidimensional- scaling-tp3422670p3422670.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] I think I just broke R
Am 03.04.2011 03:51, schrieb Daniel Malter: Check whether x, y, or glm have been redefined. If not, restart R. I wouldn't call my function 'glm'. However, I did call one 'binomial'. That was my mistake. Thanks :) A few weeks ago I asked how to set my error messages to english, and Richard Heiberger told me to use 'Sys.setenv(LANG=EN)'. He used this example, which did work for me at first, but doesn't work now anymore: Sys.setenv(LANG=DE) 2+a Fehler in 2 + a : nicht-numerisches Argument für binären Operator Sys.setenv(LANG=EN) 2+a Fehler in 2 + a : nicht-numerisches Argument für binären Operator Does someone have any idea why that could be the case? My sessionInfo() is here: sessionInfo() R version 2.10.1 (2009-12-14) i486-pc-linux-gnu locale: [1] LC_CTYPE=de_DE.utf8 LC_NUMERIC=C [3] LC_TIME=de_DE.utf8LC_COLLATE=de_DE.utf8 [5] LC_MONETARY=C LC_MESSAGES=de_DE.utf8 [7] LC_PAPER=de_DE.utf8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=de_DE.utf8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Unbalanced Anova: What is the best approach?
I have a three-way unbalanced ANOVA that I need to calculate (fixed effects plus interactions, no random effects). But word has it that aov() is good only for balanced designs. I have seen a number of different recommendations for working with unbalanced designs, but they seem to differ widely (car, nlme, lme4, etc.). So I would like to know what is the best or most usual way to go about working with unbalanced designs and extracting a reliable ANOVA table from them in R? [[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] Unbalanced Anova: What is the best approach?
On Apr 3, 2011, at 09:24 , Krishna Kirti Das wrote: I have a three-way unbalanced ANOVA that I need to calculate (fixed effects plus interactions, no random effects). But word has it that aov() is good only for balanced designs. I have seen a number of different recommendations for working with unbalanced designs, but they seem to differ widely (car, nlme, lme4, etc.). So I would like to know what is the best or most usual way to go about working with unbalanced designs and extracting a reliable ANOVA table from them in R? Actually, without random effects, aov() is not too crazy, but you might as well use plain lm(). In both cases, the main point is that you need to be aware that there is no such thing as the ANOVA table: Sums of squares will depend on the order of testing, and there is nothing to do about that (except getting balanced data). Pragmatically, I'd test the three-factor interaction, then use drop1() on a model with two-factor interactions, if nothing glaringly obvious pops up, try reduction to additive model and then use drop1() again. Obviously, if significant interactions appear, you cannot just remove them and need to investigate what they mean. -- 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] Discretizing data rows into regular intervals
On Sat, Apr 2, 2011 at 9:31 PM, Linh Tran tra...@berkeley.edu wrote: Hi guys, I'd like to thank you ahead of time for any help that you can offer me. I'm kind of stuck trying to do this. I have a data frame with dates and values (note: only two columns shown): head(test) date value stop 1 01/02/05 100 12/01/07 2 07/16/05 200 12/01/07 3 12/20/05 150 12/01/07 4 04/01/06 250 12/01/07 5 10/01/06 10 12/01/07 What I need to do is create regularly spaced 3-month intervals (starting with the first observed date) with values that are closest to but recorded after the date created. I would stop at the stop date. So the result would look like: new_date value 1 01/02/05 100 2 04/02/05 100 3 07/02/05 100 4 10/02/05 200 5 01/02/06 150 6 04/02/06 250 7 07/02/06 250 8 10/02/06 10 9 01/02/07 10 etc etc etc 10/02/07 --- ## Final obs since next one would be 1/2/08 (after stop date) See question #13 in the zoo-faq vignette: http://cran.r-project.org/web/packages/zoo/index.html and note the existence of zoo's yearqtr class. -- 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] I think I just broke R
On 03.04.2011 09:30, Alexander Engelhardt wrote: Am 03.04.2011 03:51, schrieb Daniel Malter: Check whether x, y, or glm have been redefined. If not, restart R. I wouldn't call my function 'glm'. However, I did call one 'binomial'. That was my mistake. Thanks :) A few weeks ago I asked how to set my error messages to english, and Richard Heiberger told me to use 'Sys.setenv(LANG=EN)'. He used this example, which did work for me at first, but doesn't work now anymore: Sys.setenv(LANG=DE) 2+a Fehler in 2 + a : nicht-numerisches Argument für binären Operator Sys.setenv(LANG=EN) 2+a Fehler in 2 + a : nicht-numerisches Argument für binären Operator Does someone have any idea why that could be the case? Use LANGUAGE rather than LANG as the environment variable. My sessionInfo() is here: sessionInfo() R version 2.10.1 (2009-12-14) and time to upgrade R Best, Uwe Ligges i486-pc-linux-gnu locale: [1] LC_CTYPE=de_DE.utf8 LC_NUMERIC=C [3] LC_TIME=de_DE.utf8 LC_COLLATE=de_DE.utf8 [5] LC_MONETARY=C LC_MESSAGES=de_DE.utf8 [7] LC_PAPER=de_DE.utf8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=de_DE.utf8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] kernel density plot
I am using the following commands for plotting kernel density for three kinds of crops density(s22$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-t plot(t, xlim=c(-3,4), main=Net Income Distribution, axes=F, ylim=c(0,0.00035). xlab=Value in Rupees) par(new=T) density(s33$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-u plot(u, xlim=c(-3,4), axes=F, main=, col=red, ylim=c(0,0.00035)) par(new=T) density(s44$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-v plot(v, xlim=c(-3,4), col=blue, axes=F, main=, ylim=c(0,0.00035)) the problem is that in the graph that is drawn 1. the xlab gets hidden with the [N= and the bandwidth=] values 2. when i do par(new=T) this N and bandwidth value appears multiple times..overlapping each time and making the graph look untidy.. Is there any way of making these N and Bandwidth values not appear in the graph? Thanks -- -- [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Plotting data on a US County Map
Hi, I have a data frame listing US counties and a quantity (number) per county and I have a shapefile of the US with county ID's. I would like to plot the number variable on a map (in the shapefile) using a color range per county (e.g. white = min(number) = 2, black = max(number) = 15). Can anyone help me actually plotting the data on the map? This is how far I got. Thanks! DF = data.frame(read.table(textConnection( A CNTY_FIPS number 1 US001 2 2 US002 8 3 US003 3 4 US004 5 5 US005 6 6 US006 7 7 US007 9 8 US008 9 9 US009 10 10 US010 11 11 US011 13 12 US012 15),head=TRUE,stringsAsFactors=FALSE)) library(maptools) library(ggplot2) library(gpclib) gpclibPermit() setwd(C:/Documents) us_counties.shp - readShapeSpatial(uscounties.shp) us_counties.shp.p - fortify.SpatialPolygonsDataFrame(us_counties.shp, region=CNTY_FIPS) us - merge(us_counties.shp.p, us_counties.shp, by.x=id, by.y=CNTY_FIPS) p - ggplot(data=us, aes(x=long, y=lat, group=group)) + geom_polygon(fill=#CFCFCF) p - p + geom_path(color=white) + coord_equal() ggsave(p, width=11.69, height=8.27, file=us_counties_a.jpg) -- View this message in context: http://r.789695.n4.nabble.com/Plotting-data-on-a-US-County-Map-tp3423342p3423342.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] I think I just broke R
On Sun, 3 Apr 2011, Uwe Ligges wrote: On 03.04.2011 09:30, Alexander Engelhardt wrote: Am 03.04.2011 03:51, schrieb Daniel Malter: Check whether x, y, or glm have been redefined. If not, restart R. I wouldn't call my function 'glm'. However, I did call one 'binomial'. That was my mistake. Thanks :) A few weeks ago I asked how to set my error messages to english, and Richard Heiberger told me to use 'Sys.setenv(LANG=EN)'. He used this example, which did work for me at first, but doesn't work now anymore: Sys.setenv(LANG=DE) 2+a Fehler in 2 + a : nicht-numerisches Argument für binären Operator Sys.setenv(LANG=EN) 2+a Fehler in 2 + a : nicht-numerisches Argument für binären Operator Does someone have any idea why that could be the case? Use LANGUAGE rather than LANG as the environment variable. Also, set it outside your R session, e.g. in your .Renviron file. You are supposed to be able to change this during an R session, but if you rely on OS facilities (as you probably do on Linux) rather than the gettext in the R sources, we have seen instances of the OS breaking this. My sessionInfo() is here: sessionInfo() R version 2.10.1 (2009-12-14) and time to upgrade R Best, Uwe Ligges i486-pc-linux-gnu locale: [1] LC_CTYPE=de_DE.utf8 LC_NUMERIC=C [3] LC_TIME=de_DE.utf8 LC_COLLATE=de_DE.utf8 [5] LC_MONETARY=C LC_MESSAGES=de_DE.utf8 [7] LC_PAPER=de_DE.utf8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=de_DE.utf8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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. -- 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] I think I just broke R
Am 03.04.2011 13:12, schrieb Prof Brian Ripley: On Sun, 3 Apr 2011, Uwe Ligges wrote: Use LANGUAGE rather than LANG as the environment variable. Also, set it outside your R session, e.g. in your .Renviron file. You are supposed to be able to change this during an R session, but if you rely on OS facilities (as you probably do on Linux) rather than the gettext in the R sources, we have seen instances of the OS breaking this. I did use LANGUAGE too, didn't work as well. Creating a ~/.Rprofile file with LANGUAGE=EN had no effect (weird..). When I edited /etc/R/Renviron.site to include LANGUAGE=EN, it worked. Thanks for the hints! sessionInfo() R version 2.10.1 (2009-12-14) and time to upgrade R I'm still fighting to find out how to upgrade stuff on Ubuntu. After a repository update the newest available version was still 2.10.1. I'll figure it out, sooner or later :) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] kernel density plot
On Apr 3, 2011, at 6:56 AM, Muzna Alvi wrote: I am using the following commands for plotting kernel density for three kinds of crops density(s22$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-t plot(t, xlim=c(-3,4), main=Net Income Distribution, axes=F, ylim=c(0,0.00035). xlab=Value in Rupees) par(new=T) density(s33$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-u plot(u, xlim=c(-3,4), axes=F, main=, col=red, ylim=c(0,0.00035)) par(new=T) density(s44$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-v plot(v, xlim=c(-3,4), col=blue, axes=F, main=, ylim=c(0,0.00035)) the problem is that in the graph that is drawn 1. the xlab gets hidden with the [N= and the bandwidth=] values 2. when i do par(new=T) this N and bandwidth value appears multiple times..overlapping each time and making the graph look untidy.. Is there any way of making these N and Bandwidth values not appear in the graph? Why not just set ylab= in subsequent calls to plot? -- David Winsemius, MD 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] kernel density plot
On Apr 3, 2011, at 7:55 AM, David Winsemius wrote: On Apr 3, 2011, at 6:56 AM, Muzna Alvi wrote: I am using the following commands for plotting kernel density for three kinds of crops density(s22$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-t plot(t, xlim=c(-3,4), main=Net Income Distribution, axes=F, ylim=c(0,0.00035). xlab=Value in Rupees) par(new=T) density(s33$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-u plot(u, xlim=c(-3,4), axes=F, main=, col=red, ylim=c(0,0.00035)) par(new=T) density(s44$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-v plot(v, xlim=c(-3,4), col=blue, axes=F, main=, ylim=c(0,0.00035)) the problem is that in the graph that is drawn 1. the xlab gets hidden with the [N= and the bandwidth=] values 2. when i do par(new=T) this N and bandwidth value appears multiple times..overlapping each time and making the graph look untidy.. Is there any way of making these N and Bandwidth values not appear in the graph? Why not just set ylab= in subsequent calls to plot? Sorry, I meant xlab=. -- David Winsemius, MD 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] help
Date: Sun, 3 Apr 2011 01:35:16 +0530 From: nandan.a...@gmail.com To: padmanabhan.vija...@gmail.com CC: r-help@r-project.org Subject: Re: [R] help One way that u might have thought of is to create plot in PDF in R and the use pdftools. Additionally one can also think of running R script using R CMD and then using pdftools in a .sh script file if u r in linux. I am not aware of pdftools capability in R. On 2 April 2011 23:01, Vijayan Padmanabhan wrote: Dear R Help group I need to run a command line script from within R session. I am not clear how i can acheive this. I tried shell and system function, but i am missing something critical.can someone provide help? My intention is to create a pdf file of a plot in R and then attach existing files from my system as attachment into the newly created pdf file. Any help would be greatly appreciated.. Here is the command line script i want to execute from within R. pdftools -S attachfiles=C:\test1.pdf -i C:\test2.pdf -o C:\test4.pdf Regards Vijayan Padmanabhan I just tried system(pdftk --help) and it appeared to work as I have pdftk from cygwin.I routinely do this the other way however and invoke R from a bash script and then use external tools like this from the bash script after R is done. If I'm generating various pieces, it seems to make sense to get them all first and release any resources R has accumulated as pdf manipulation itself can often require lots of memory etc. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Error in color2D.matplot : Error in plot.new() : figure margins too large
Hi, I am using color2D.matplot (...) function of plotrix package. I used a matrix of size around 20*20 However, apparently it failed to visualize the matrix and gave the following exception, which I don't have any idea about possible source of this error. Error in plot.new() : figure margins too large It would be appreciated if someone points me to the right origin of this error. best, /Shahab __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] I think I just broke R
On 11-04-03 7:50 AM, Alexander Engelhardt wrote: Am 03.04.2011 13:12, schrieb Prof Brian Ripley: On Sun, 3 Apr 2011, Uwe Ligges wrote: Use LANGUAGE rather than LANG as the environment variable. Also, set it outside your R session, e.g. in your .Renviron file. You are supposed to be able to change this during an R session, but if you rely on OS facilities (as you probably do on Linux) rather than the gettext in the R sources, we have seen instances of the OS breaking this. I did use LANGUAGE too, didn't work as well. Creating a ~/.Rprofile file with LANGUAGE=EN had no effect (weird..). That's not weird: you just created an R variable named LANGUAGE, not an environment variable. Duncan Murdoch When I edited /etc/R/Renviron.site to include LANGUAGE=EN, it worked. Thanks for the hints! sessionInfo() R version 2.10.1 (2009-12-14) and time to upgrade R I'm still fighting to find out how to upgrade stuff on Ubuntu. After a repository update the newest available version was still 2.10.1. I'll figure it out, sooner or later :) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] I think I just broke R
Am 03.04.2011 14:10, schrieb Duncan Murdoch: That's not weird: you just created an R variable named LANGUAGE, not an environment variable. Duncan Murdoch Silly me. It works now: alexx@derp:~$ cat ~/.Renviron LANGUAGE=EN 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.
Re: [R] Error in color2D.matplot : Error in plot.new() : figure margins too large
Totally guessing but did you make an earlier call to par() to adjust margins? Otherwise you may want to supply the matrix here for other people to experiment with. See ?dput as probably the best way to do this. --- On Sun, 4/3/11, shahab shahab.mok...@gmail.com wrote: From: shahab shahab.mok...@gmail.com Subject: [R] Error in color2D.matplot : Error in plot.new() : figure margins too large To: r-help@r-project.org Received: Sunday, April 3, 2011, 8:02 AM Hi, I am using color2D.matplot (...) function of plotrix package. I used a matrix of size around 20*20 However, apparently it failed to visualize the matrix and gave the following exception, which I don't have any idea about possible source of this error. Error in plot.new() : figure margins too large It would be appreciated if someone points me to the right origin of this error. best, /Shahab __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] recommendation on r scripting tutorial?
The documents accompanying the distribution can be a good start. And a suggestion for searching help on R over Google, use r-help as a basic keyword, coz a single letter of r hardly helps you find the desired topics. 2011/4/2 Wensui Liu liuwen...@gmail.com Good morning, dear listers I am wondering if you could recommend a good tutorial / book for r scripting. thank you so much in advance! WenSui Liu Credit Risk Manager, 53 Bancorp wensui@53.com 513-295-4370 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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.
[R] Inverse noncentral Beta
Hello I could not find whether there is any R-function that implements the inverse of a noncentral Beta. Could someone out there tell me where I can find it? Or how to implement it? Many thanks Ed [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R gui on windows how to force to always show the last line of output
CTRL+L helps. 2011/4/2 stan zimine szim...@gmail.com Hi. Googled but did not found the answer for the following little issue. how to force R gui on windows (maybe a specific setting) to always show the last line of output in the window console. My program in R makes measurements every 5 mins in indefinite loop and prints results in the console. The problem: last messages are not visible, The scrolling bar of the gui console gets shorter. I.e. you have to scroll for the last messages. Thanks if anybody knows the sol to this prob. SZ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] recommendation on r scripting tutorial?
On 11-04-03 9:01 AM, lcn wrote: The documents accompanying the distribution can be a good start. And a suggestion for searching help on R over Google, use r-help as a basic keyword, coz a single letter of r hardly helps you find the desired topics. When I google for R tutorial or r tutorial, the entire first page looks relevant (though I'm not familiar with most of the tutorials, so can't give a recommendation). I think it's a myth that Google doesn't know what you mean when you ask about R. Or perhaps it is tailoring its results to what it has seen me choose in the past. Duncan Murdoch 2011/4/2 Wensui Liuliuwen...@gmail.com Good morning, dear listers I am wondering if you could recommend a good tutorial / book for r scripting. thank you so much in advance! WenSui Liu Credit Risk Manager, 53 Bancorp wensui@53.com 513-295-4370 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Inverse noncentral Beta
On 03-Apr-11 13:03:30, Eduardo M. A. M.Mendes wrote: Hello I could not find whether there is any R-function that implements the inverse of a noncentral Beta. Could someone out there tell me where I can find it? Or how to implement it? Many thanks Ed Have a look at the 'ncp' paramater in ?qbeta -- this (default=0) is the non-centrality parameter for the beta distribution, and qbeta is the inverse function of the distribution. The noncentral Beta distribution (with ?ncp? = lambda) is defined (Johnson et al, 1995, pp. 502) as the distribution of X/(X+Y) where X ~ chi^2_2a(lambda) and Y ~ chi^2_2b. i.e. as in qbeta(p,a,b,ncp=lambda) Hoping this helps, Ted. E-Mail: (Ted Harding) ted.hard...@wlandres.net Fax-to-email: +44 (0)870 094 0861 Date: 03-Apr-11 Time: 14:17:58 -- XFMail -- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Unbalanced Anova: What is the best approach?
Dear Krishna, Although it's difficult to explain briefly, I'd argue that balanced and unbalanced ANOVA are not fundamentally different, in that the focus should be on the hypotheses that are tested, and these are naturally expressed as functions of cell means and marginal means. For example, in a two-way ANOVA, the null hypotheses of no interaction is equivalent to parallel profiles of cell means for one factor across levels of the other. What is different, though, is that in a balanced ANOVA all common approaches to constructing an ANOVA table coincide. Without getting into the explanation in detail (which you can find in a text like my Applied Regression Analysis and Generalized Linear Models), so-called type-I (or sequential) tests, such as those performed by the standard anova() function in R, test hypotheses that are rarely of substantive interest, and, even when they are, are of interest only by accident. So-called type-II tests, such as those performed by default by the Anova() function in the car package, test hypotheses that are almost always of interest. Type-III tests, which the Anova() function in car can perform optionally, require careful formulation of the model for the hypotheses tested to be sensible, and even then have less power than corresponding type-II tests in the circumstances in which a test would be of interest. Since you're addressing fixed-effects models, I'm not sure why you introduced nlme and lme4 into the discussion, but I note that Anova() in the car package has methods that can produce type-II and -III Wald tests for the fixed effects in mixed models fit by lme() and lmer(). Your question has been asked several times before on the r-help list. For example, if you enter terms like type-II or unbalanced ANOVA in the RSeek search engine and look under the Support Lists tab, you'll see many hits -- e.g., Mhttps://stat.ethz.ch/pipermail/r-help/2006-August/111927.html. I hope this helps, John John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Krishna Kirti Das Sent: April-03-11 3:25 AM To: r-help@r-project.org Subject: [R] Unbalanced Anova: What is the best approach? I have a three-way unbalanced ANOVA that I need to calculate (fixed effects plus interactions, no random effects). But word has it that aov() is good only for balanced designs. I have seen a number of different recommendations for working with unbalanced designs, but they seem to differ widely (car, nlme, lme4, etc.). So I would like to know what is the best or most usual way to go about working with unbalanced designs and extracting a reliable ANOVA table from them in R? [[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] Unbalanced Anova: What is the best approach?
Thank you, John. Yes, your answers do help. For me it's mainly about getting familiar with the R way of doing things. Thus your response also confirms what I suspected, that there is no explicit user-interface (at least one that is widely used) in terms of functions/packages that represents an unbalanced design in the same way that aov would represent a balanced one. Analyzing balanced and unbalanced data are obviously possible, but with balanced designs via aov what has to be done is intuitive within the language but unintuitive for unbalanced designs. I did notice that this question gets asked several times and in slightly different ways, and I think the lack of an interface that represents an unbalanced design in the same way aov represents balanced designs is why the question will probably keep getting asked again. I had mentioned nlme and lme4 because I saw in some of the discussions that using those were recommended for working with unbalanced designs. And specifying random effects with zero variance, for example, would probably serve my purposes. Thank you for your help. Sincerely, Krishna On Sun, Apr 3, 2011 at 7:28 AM, John Fox j...@mcmaster.ca wrote: Dear Krishna, Although it's difficult to explain briefly, I'd argue that balanced and unbalanced ANOVA are not fundamentally different, in that the focus should be on the hypotheses that are tested, and these are naturally expressed as functions of cell means and marginal means. For example, in a two-way ANOVA, the null hypotheses of no interaction is equivalent to parallel profiles of cell means for one factor across levels of the other. What is different, though, is that in a balanced ANOVA all common approaches to constructing an ANOVA table coincide. Without getting into the explanation in detail (which you can find in a text like my Applied Regression Analysis and Generalized Linear Models), so-called type-I (or sequential) tests, such as those performed by the standard anova() function in R, test hypotheses that are rarely of substantive interest, and, even when they are, are of interest only by accident. So-called type-II tests, such as those performed by default by the Anova() function in the car package, test hypotheses that are almost always of interest. Type-III tests, which the Anova() function in car can perform optionally, require careful formulation of the model for the hypotheses tested to be sensible, and even then have less power than corresponding type-II tests in the circumstances in which a test would be of interest. Since you're addressing fixed-effects models, I'm not sure why you introduced nlme and lme4 into the discussion, but I note that Anova() in the car package has methods that can produce type-II and -III Wald tests for the fixed effects in mixed models fit by lme() and lmer(). Your question has been asked several times before on the r-help list. For example, if you enter terms like type-II or unbalanced ANOVA in the RSeek search engine and look under the Support Lists tab, you'll see many hits -- e.g., Mhttps://stat.ethz.ch/pipermail/r-help/2006-August/111927.html. I hope this helps, John John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Krishna Kirti Das Sent: April-03-11 3:25 AM To: r-help@r-project.org Subject: [R] Unbalanced Anova: What is the best approach? I have a three-way unbalanced ANOVA that I need to calculate (fixed effects plus interactions, no random effects). But word has it that aov() is good only for balanced designs. I have seen a number of different recommendations for working with unbalanced designs, but they seem to differ widely (car, nlme, lme4, etc.). So I would like to know what is the best or most usual way to go about working with unbalanced designs and extracting a reliable ANOVA table from them in R? [[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. [[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] Unbalanced Anova: What is the best approach?
On Sun, Apr 3, 2011 at 3:10 AM, peter dalgaard pda...@gmail.com wrote: On Apr 3, 2011, at 09:24 , Krishna Kirti Das wrote: I have a three-way unbalanced ANOVA that I need to calculate (fixed effects plus interactions, no random effects). But word has it that aov() is good only for balanced designs. I have seen a number of different recommendations for working with unbalanced designs, but they seem to differ widely (car, nlme, lme4, etc.). So I would like to know what is the best or most usual way to go about working with unbalanced designs and extracting a reliable ANOVA table from them in R? Actually, without random effects, aov() is not too crazy, but you might as well use plain lm(). In both cases, the main point is that you need to be aware that there is no such thing as the ANOVA table: Sums of squares will depend on the order of testing, and there is nothing to do about that (except getting balanced data). Pragmatically, I'd test the three-factor interaction, then use drop1() on a model with two-factor interactions, if nothing glaringly obvious pops up, try reduction to additive model and then use drop1() again. Obviously, if significant interactions appear, you cannot just remove them and need to investigate what they mean. That helps. And I've been looking for something like drop1() for a while now. Thank you. Sincerely, Krishna [[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
Hello, nbsp; I want to sum first three terms of each column of matrix. But I don't calculate with apply function. nbsp; skwkrtlt;-function(N=1,mu=0,sigma=1,n=100, nboot=1000,alpha=0.05){ xlt;-rnorm(N,mu,sigma)#population samplexlt;-matrix(sample(x,n*nboot,replace=T),nrow=nboot) #... } nbsp; is that: suppose a is a 5x2 matrix. nbsp;a={1,2,3,4,5 6,7,8,9,10} nbsp; I want to sum first three terms. sm[1]=1+2+3 sm[2]=6+7+8 But I don't calculate. please help me!!! [[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] Unbalanced Anova: What is the best approach?
Hi, Krishna: in line On 4/3/2011 7:35 AM, Krishna Kirti Das wrote: Thank you, John. Yes, your answers do help. For me it's mainly about getting familiar with the R way of doing things. Thus your response also confirms what I suspected, that there is no explicit user-interface (at least one that is widely used) in terms of functions/packages that represents an unbalanced design in the same way that aov would represent a balanced one. Analyzing balanced and unbalanced data are obviously possible, but with balanced designs via aov what has to be done is intuitive within the language but unintuitive for unbalanced designs. Intuition is subject to one's background and expectations. If you think in terms of a series of nested hypotheses, then the standard R anova is very intuitive. I never use aov, because it's not intuitive to me and not very general. 'aov' is only useful for a balanced design with normal independent errors with constant variance. The real world is rarely so simple. The 'aov' algorithm was wonderful over half a century ago, when all computations were done by hand or using a mechanical calculator (e.g., an abacus or a calculator with gears). Unbalanced designs were largely impractical because of computational difficulties. There were many procedures for imputing missing values for a design that was almost balanced. I encourage you to think in terms of alternative sequences of nested hypotheses, including the implications of A being significant by itself, but not with B already present, except that the A:B interaction is or is not significant. I did notice that this question gets asked several times and in slightly different ways, and I think the lack of an interface that represents an unbalanced design in the same way aov represents balanced designs is why the question will probably keep getting asked again. I had mentioned nlme and lme4 because I saw in some of the discussions that using those were recommended for working with unbalanced designs. And specifying random effects with zero variance, for example, would probably serve my purposes. I'd be surprised if nlme or lme4 changes what I wrote above. Hope this helps. Spencer Thank you for your help. Sincerely, Krishna On Sun, Apr 3, 2011 at 7:28 AM, John Foxj...@mcmaster.ca wrote: Dear Krishna, Although it's difficult to explain briefly, I'd argue that balanced and unbalanced ANOVA are not fundamentally different, in that the focus should be on the hypotheses that are tested, and these are naturally expressed as functions of cell means and marginal means. For example, in a two-way ANOVA, the null hypotheses of no interaction is equivalent to parallel profiles of cell means for one factor across levels of the other. What is different, though, is that in a balanced ANOVA all common approaches to constructing an ANOVA table coincide. Without getting into the explanation in detail (which you can find in a text like my Applied Regression Analysis and Generalized Linear Models), so-called type-I (or sequential) tests, such as those performed by the standard anova() function in R, test hypotheses that are rarely of substantive interest, and, even when they are, are of interest only by accident. So-called type-II tests, such as those performed by default by the Anova() function in the car package, test hypotheses that are almost always of interest. Type-III tests, which the Anova() function in car can perform optionally, require careful formulation of the model for the hypotheses tested to be sensible, and even then have less power than corresponding type-II tests in the circumstances in which a test would be of interest. Since you're addressing fixed-effects models, I'm not sure why you introduced nlme and lme4 into the discussion, but I note that Anova() in the car package has methods that can produce type-II and -III Wald tests for the fixed effects in mixed models fit by lme() and lmer(). Your question has been asked several times before on the r-help list. For example, if you enter terms like type-II or unbalanced ANOVA in the RSeek search engine and look under the Support Lists tab, you'll see many hits -- e.g., Mhttps://stat.ethz.ch/pipermail/r-help/2006-August/111927.html. I hope this helps, John John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Krishna Kirti Das Sent: April-03-11 3:25 AM To: r-help@r-project.org Subject: [R] Unbalanced Anova: What is the best approach? I have a three-way unbalanced ANOVA that I need to calculate (fixed effects plus interactions, no random effects). But word has it that aov() is good only for balanced designs. I have seen a number of different recommendations
Re: [R] Unbalanced Anova: What is the best approach?
Dear Spencer, -Original Message- From: Spencer Graves [mailto:spencer.gra...@prodsyse.com] Sent: April-03-11 11:07 AM To: Krishna Kirti Das Cc: John Fox; r-help@r-project.org Subject: Re: [R] Unbalanced Anova: What is the best approach? Hi, Krishna: in line On 4/3/2011 7:35 AM, Krishna Kirti Das wrote: Thank you, John. Yes, your answers do help. For me it's mainly about getting familiar with the R way of doing things. Thus your response also confirms what I suspected, that there is no explicit user-interface (at least one that is widely used) in terms of functions/packages that represents an unbalanced design in the same way that aov would represent a balanced one. Analyzing balanced and unbalanced data are obviously possible, but with balanced designs via aov what has to be done is intuitive within the language but unintuitive for unbalanced designs. Intuition is subject to one's background and expectations. If you think in terms of a series of nested hypotheses, then the standard R anova is very intuitive. I never use aov, because it's not intuitive to me and not very general. 'aov' is only useful for a balanced design with normal independent errors with constant variance. The real world is rarely so simple. The 'aov' algorithm was wonderful over half a century ago, when all computations were done by hand or using a mechanical calculator (e.g., an abacus or a calculator with gears). Unbalanced designs were largely impractical because of computational difficulties. There were many procedures for imputing missing values for a design that was almost balanced. I encourage you to think in terms of alternative sequences of nested hypotheses, including the implications of A being significant by itself, but not with B already present, except that the A:B interaction is or is not significant. So-called type-II tests do exactly that -- that is, obey the principle of marginality; they are maximally powerful if the higher-order term(s) to which a particular term is marginal are 0. Best, John I did notice that this question gets asked several times and in slightly different ways, and I think the lack of an interface that represents an unbalanced design in the same way aov represents balanced designs is why the question will probably keep getting asked again. I had mentioned nlme and lme4 because I saw in some of the discussions that using those were recommended for working with unbalanced designs. And specifying random effects with zero variance, for example, would probably serve my purposes. I'd be surprised if nlme or lme4 changes what I wrote above. Hope this helps. Spencer Thank you for your help. Sincerely, Krishna On Sun, Apr 3, 2011 at 7:28 AM, John Foxj...@mcmaster.ca wrote: Dear Krishna, Although it's difficult to explain briefly, I'd argue that balanced and unbalanced ANOVA are not fundamentally different, in that the focus should be on the hypotheses that are tested, and these are naturally expressed as functions of cell means and marginal means. For example, in a two-way ANOVA, the null hypotheses of no interaction is equivalent to parallel profiles of cell means for one factor across levels of the other. What is different, though, is that in a balanced ANOVA all common approaches to constructing an ANOVA table coincide. Without getting into the explanation in detail (which you can find in a text like my Applied Regression Analysis and Generalized Linear Models), so-called type-I (or sequential) tests, such as those performed by the standard anova() function in R, test hypotheses that are rarely of substantive interest, and, even when they are, are of interest only by accident. So-called type-II tests, such as those performed by default by the Anova() function in the car package, test hypotheses that are almost always of interest. Type-III tests, which the Anova() function in car can perform optionally, require careful formulation of the model for the hypotheses tested to be sensible, and even then have less power than corresponding type-II tests in the circumstances in which a test would be of interest. Since you're addressing fixed-effects models, I'm not sure why you introduced nlme and lme4 into the discussion, but I note that Anova() in the car package has methods that can produce type-II and -III Wald tests for the fixed effects in mixed models fit by lme() and lmer(). Your question has been asked several times before on the r-help list. For example, if you enter terms like type-II or unbalanced ANOVA in the RSeek search engine and look under the Support Lists tab, you'll see many hits -- e.g., Mhttps://stat.ethz.ch/pipermail/r-help/2006-August/111927.html. I hope this helps, John John Fox Senator William McMaster
Re: [R] Unbalanced Anova: What is the best approach?
Dear Krishna, -Original Message- From: Krishna Kirti Das [mailto:krishnaki...@gmail.com] Sent: April-03-11 10:36 AM To: John Fox Cc: r-help@r-project.org Subject: Re: [R] Unbalanced Anova: What is the best approach? Thank you, John. Yes, your answers do help. For me it's mainly about getting familiar with the R way of doing things. Thus your response also confirms what I suspected, that there is no explicit user-interface (at least one that is widely used) in terms of functions/packages that represents an unbalanced design in the same way that aov would represent a balanced one. Analyzing balanced and unbalanced data are obviously possible, but with balanced designs via aov what has to be done is intuitive within the language but unintuitive for unbalanced designs. I don't agree with your characterization. For example, the representation of a two-way crossed ANOVA model as an R model formula is precisely the same for balanced and unbalanced data: for response Y and factors A and B, Y ~ A*B. Moreover, the issue of how to formulate tests is independent of the software you use. I did notice that this question gets asked several times and in slightly different ways, and I think the lack of an interface that represents an unbalanced design in the same way aov represents balanced designs is why the question will probably keep getting asked again. I suspect that the issue gets asked repeatedly for two reasons: (1) More fundamentally, I believe that the general level of understanding of hypothesis tests in unbalanced data is low; (2) people don't necessarily read previous posts to r-help. I had mentioned nlme and lme4 because I saw in some of the discussions that using those were recommended for working with unbalanced designs. And specifying random effects with zero variance, for example, would probably serve my purposes. I don't think that either lme() or lmer() will allow you to fit a model without random effects, but even if they did there wouldn't be much sense in doing so. You can compute a mean with lm() or glm(), but would you? Best, John Thank you for your help. Sincerely, Krishna On Sun, Apr 3, 2011 at 7:28 AM, John Fox j...@mcmaster.ca wrote: Dear Krishna, Although it's difficult to explain briefly, I'd argue that balanced and unbalanced ANOVA are not fundamentally different, in that the focus should be on the hypotheses that are tested, and these are naturally expressed as functions of cell means and marginal means. For example, in a two-way ANOVA, the null hypotheses of no interaction is equivalent to parallel profiles of cell means for one factor across levels of the other. What is different, though, is that in a balanced ANOVA all common approaches to constructing an ANOVA table coincide. Without getting into the explanation in detail (which you can find in a text like my Applied Regression Analysis and Generalized Linear Models), so-called type-I (or sequential) tests, such as those performed by the standard anova() function in R, test hypotheses that are rarely of substantive interest, and, even when they are, are of interest only by accident. So-called type-II tests, such as those performed by default by the Anova() function in the car package, test hypotheses that are almost always of interest. Type-III tests, which the Anova() function in car can perform optionally, require careful formulation of the model for the hypotheses tested to be sensible, and even then have less power than corresponding type-II tests in the circumstances in which a test would be of interest. Since you're addressing fixed-effects models, I'm not sure why you introduced nlme and lme4 into the discussion, but I note that Anova() in the car package has methods that can produce type-II and -III Wald tests for the fixed effects in mixed models fit by lme() and lmer(). Your question has been asked several times before on the r-help list. For example, if you enter terms like type-II or unbalanced ANOVA in the RSeek search engine and look under the Support Lists tab, you'll see many hits -- e.g., Mhttps://stat.ethz.ch/pipermail/r-help/2006-August/111927.html. I hope this helps, John John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of Krishna Kirti Das Sent: April-03-11 3:25 AM To: r-help@r-project.org Subject: [R] Unbalanced Anova: What is the
Re: [R] :HELP
Your message is severely garbled. However a={1,2,3,4,5 6,7,8,9,10} is not going to work since {} is incorrect. you need c() Anyway try this: a=matrix(1:10,ncol=2) sum(a[1:3,1]) sum(a[1:3,2]) --- On Sun, 4/3/11, mustafabinar mustafa_bi...@mynet.com wrote: From: mustafabinar mustafa_bi...@mynet.com Subject: [R] :HELP To: r-help@r-project.org Received: Sunday, April 3, 2011, 9:10 AM Hello, nbsp; I want to sum first three terms of each column of matrix. But I don't calculate with apply function. nbsp; skwkrtlt;-function(N=1,mu=0,sigma=1,n=100, nboot=1000,alpha=0.05){ xlt;-rnorm(N,mu,sigma)#population samplexlt;-matrix(sample(x,n*nboot,replace=T),nrow=nboot) #... } nbsp; is that: suppose a is a 5x2 matrix. nbsp;a={1,2,3,4,5 6,7,8,9,10} nbsp; I want to sum first three terms. sm[1]=1+2+3 sm[2]=6+7+8 But I don't calculate. please help me!!! [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] converting call objects into character
Dear all, I would like to log the calls to my functions. I am trying to do this using the function match.call(): fTest-function(x) { theCall-match.call() print(theCall) return(x) } fTest(2) fTest(x = 2) [1] 2 I can see theCall printed into the console, but I don't manage to convert it into a character to write it into a log file with other informations. Can anyone help? Many thanks, Samuel [[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] converting call objects into character
On Apr 3, 2011, at 12:14 PM, Samuel Le wrote: Dear all, I would like to log the calls to my functions. I am trying to do this using the function match.call(): fTest-function(x) { theCall-match.call() print(theCall) return(list(x=x, logf = theCall)) } fTest(x=2)$x [1] 2 fTest(x=2)$logf fTest(x = 2) str(fTest(x=2)$logf) language fTest(x = 2) You may want to convert that call component to a character object, since: cat(fTest(x=2)$logf) Error in cat(list(...), file, sep, fill, labels, append) : argument 1 (type 'language') cannot be handled by 'cat' I can see theCall printed into the console, but I don't manage to convert it into a character to write it into a log file with other informations. Can anyone help? David Winsemius, MD 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 in splitting ists into sub-lists
Dear List, Let's say I have a list whose components are 2 matrices (as exemplified in the mylist object below). I'd like to create a list with components being 4 matrices based on an logical index vector. is there a way to simplify what I'm doing to obtain the results in mylist2? I'd like something that would work on an arbitrary number of elements in mylist. mylist - list(matrix(1:9,3,3), matrix(10:18,3,3)) index1 - c(TRUE,FALSE,TRUE) index2- c(FALSE,TRUE,TRUE) mylist2 - list(mylist[[1]][,index1], mylist[[1]][,index1==FALSE], mylist[[2]][,index2], mylist[[2]][,index2==FALSE]) Thanks in advance, Axel. [[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] converting call objects into character
On Sun, Apr 3, 2011 at 11:42 AM, David Winsemius dwinsem...@comcast.net wrote: On Apr 3, 2011, at 12:14 PM, Samuel Le wrote: Dear all, I would like to log the calls to my functions. I am trying to do this using the function match.call(): fTest-function(x) { theCall-match.call() print(theCall) return(list(x=x, logf = theCall)) } fTest(x=2)$x [1] 2 fTest(x=2)$logf fTest(x = 2) str(fTest(x=2)$logf) language fTest(x = 2) You may want to convert that call component to a character object, since: cat(fTest(x=2)$logf) Error in cat(list(...), file, sep, fill, labels, append) : argument 1 (type 'language') cannot be handled by 'cat' If you want to examine a call object you need to ensure that it is not evaluated. Evaluating a number or a character string is not a problem because eval(4) is the same as 4 However, evaluating a function call should be different from the call itself. As David shows, the str function is careful not to evaluate the call object. (Martin and I found ourselves going around in circles when looking at the structure of a fitted model object that included a call and he kindly changed the behavior of str().) So you need to decide when a function, such as print(), evaluates its arguments or when it doesn't, which can get kind of complicated. An alternative is to use match.call() repeatedly instead of trying to save the value, as in fTest function(x) { print(match.call()) list(x=x, logf = match.call()) } fTest(x=2) fTest(x = 2) $x [1] 2 $logf fTest(x = 2) The trick there is that the value of match.call() is the unevaluated call whereas myCall - match.call() print(myCall) evaluates myCall in the call to print, thereby evaluating the function fTest again. Is this sufficiently confusing? :-) I can see theCall printed into the console, but I don't manage to convert it into a character to write it into a log file with other informations. Can anyone help? David Winsemius, MD 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] Function for finding NA's
Quick question, I tried to find a function in available packages to find NA's for an entire data set (or single variables) and report the row of missing values (NA's for each column). I searched the typical routes through the blogs and the help manuals for 15 minutes. Rather than spend any more time searching I created my own function to do this (probably in less time than it would have taken me to find the function). Now I still have the same question: Is this function (NAhunter I call it) already in existence? If so please direct me (because I'm sure they've written better code more efficiently). I highly doubt I'm this first person to want to find all the missing values in a data set so I assume there is a function for it but I just didn't spend enough time looking. If there is no existing function (big if here), is this something people feel is worthwhile for me to put into a package of some sort? Tyler Here's the code: NAhunter-function(dataset) { find.NA-function(variable) { if(is.numeric(variable)){ n-length(variable) mean-mean(variable, na.rm=T) median-median(variable, na.rm=T) sd-sd(variable, na.rm=T) NAs-is.na(variable) total.NA-sum(NAs) percent.missing-total.NA/n descriptives-data.frame(n,mean,median,sd,total.NA,percent.missing) rownames(descriptives)-c( ) Case.Number-1:n Missing.Values-ifelse(NAs0,Missing Value, ) missing.value-data.frame(Case.Number,Missing.Values) missing.values-missing.value[ which(Missing.Values=='Missing Value'),] list(NUMERIC DATA,DESCRIPTIVES=t(descriptives),CASE # OF MISSING VALUES=missing.values[,1]) } else{ n-length(variable) NAs-is.na(variable) total.NA-sum(NAs) percent.missing-total.NA/n descriptives-data.frame(n,total.NA,percent.missing) rownames(descriptives)-c( ) Case.Number-1:n Missing.Values-ifelse(NAs0,Missing Value, ) missing.value-data.frame(Case.Number,Missing.Values) missing.values-missing.value[ which(Missing.Values=='Missing Value'),] list(CATEGORICAL DATA,DESCRIPTIVES=t(descriptives),CASE # OF MISSING VALUES=missing.values[,1]) } } dataset-data.frame(dataset) options(scipen=100) options(digits=2) lapply(dataset,find.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] :HELP
Hi: Is this what you're after? x - matrix(1:10, nrow = 2, byrow = TRUE) rowSums(x[, 1:3]) HTH, Dennis PS: This list expects submissions to be in text, not HTML. See the Posting Guide, please. On Sun, Apr 3, 2011 at 6:10 AM, mustafabinar mustafa_bi...@mynet.comwrote: Hello, nbsp; I want to sum first three terms of each column of matrix. But I don't calculate with apply function. nbsp; skwkrtlt;-function(N=1,mu=0,sigma=1,n=100, nboot=1000,alpha=0.05){ xlt;-rnorm(N,mu,sigma)#population samplexlt;-matrix(sample(x,n*nboot,replace=T),nrow=nboot) #... } nbsp; is that: suppose a is a 5x2 matrix. nbsp;a={1,2,3,4,5 6,7,8,9,10} nbsp; I want to sum first three terms. sm[1]=1+2+3 sm[2]=6+7+8 But I don't calculate. please help me!!! [[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. [[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] converting call objects into character
On Apr 3, 2011, at 1:22 PM, Douglas Bates wrote: On Sun, Apr 3, 2011 at 11:42 AM, David Winsemius dwinsem...@comcast.net wrote: On Apr 3, 2011, at 12:14 PM, Samuel Le wrote: Dear all, I would like to log the calls to my functions. I am trying to do this using the function match.call(): fTest-function(x) { theCall-match.call() print(theCall) return(list(x=x, logf = theCall)) } fTest(x=2)$x [1] 2 fTest(x=2)$logf fTest(x = 2) str(fTest(x=2)$logf) language fTest(x = 2) You may want to convert that call component to a character object, since: cat(fTest(x=2)$logf) Error in cat(list(...), file, sep, fill, labels, append) : argument 1 (type 'language') cannot be handled by 'cat' If you want to examine a call object you need to ensure that it is not evaluated. Evaluating a number or a character string is not a problem because eval(4) is the same as 4 However, evaluating a function call should be different from the call itself. As David shows, the str function is careful not to evaluate the call object. (Martin and I found ourselves going around in circles when looking at the structure of a fitted model object that included a call and he kindly changed the behavior of str().) So you need to decide when a function, such as print(), evaluates its arguments or when it doesn't, which can get kind of complicated. An alternative is to use match.call() repeatedly instead of trying to save the value, as in fTest function(x) { print(match.call()) list(x=x, logf = match.call()) } fTest(x=2) fTest(x = 2) $x [1] 2 $logf fTest(x = 2) The trick there is that the value of match.call() is the unevaluated call whereas myCall - match.call() print(myCall) evaluates myCall in the call to print, thereby evaluating the function fTest again. Is this sufficiently confusing? :-) Yes, I am now sufficiently confused^W , ... er, motivated to look for another route. I think the way out of the confusion is to turn the call into text and since as.character doesn't do a very neat a job, I would suggest instead: deparse() fTest - function(x) { +print(match.call()) +list(x=x, logf = deparse(match.call())) + } fTest(x=3)$logf fTest(x = 3) [1] fTest(x = 3) cat(fTest(x=3)$logf) fTest(x = 3) fTest(x = 3) cat() is a convenient test of the capacity of an object to be written to a file. It has an append parameter that implies it could serve the logging function requested by the OP. I can see theCall printed into the console, but I don't manage to convert it into a character to write it into a log file with other informations. Can anyone help? David Winsemius, MD 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] Function for finding NA's
On Apr 3, 2011, at 1:44 PM, Tyler Rinker wrote: Quick question, I tried to find a function in available packages to find NA's for an entire data set (or single variables) and report the row of missing values (NA's for each column). I searched the typical routes through the blogs and the help manuals for 15 minutes. Rather than spend any more time searching I created my own function to do this (probably in less time than it would have taken me to find the function). Now I still have the same question: Is this function (NAhunter I call it) already in existence? If so please direct me (because I'm sure they've written better code more efficiently). I highly doubt I'm this first person to want to find all the missing values in a data set so I assume there is a function for it but I just didn't spend enough time looking. If there is no existing function (big if here), is this something people feel is worthwhile for me to put into a package of some sort? I'm not sure that it would have occurred to people to include it in a package. Consider: getNa - function(dfrm) lapply(dfrm, function(x) which(is.na(x) ) ) cities long lat city pop 1 -58.38194 -34.59972 Buenos Aires NA 2 14.25000 40.8 NA NA getNa(cities) $long integer(0) $lat integer(0) $city [1] 2 $pop [1] 1 2 There are several packages with functions by the name `describe` that do most or all of rest of what you have proposed. I happen to use Harrell's Hmisc but the other versions should also be reviewed if you want to avoid re-inventing the wheel. -- David. Tyler Here's the code: NAhunter-function(dataset) { find.NA-function(variable) { if(is.numeric(variable)){ n-length(variable) mean-mean(variable, na.rm=T) median-median(variable, na.rm=T) sd-sd(variable, na.rm=T) NAs-is.na(variable) total.NA-sum(NAs) percent.missing-total.NA/n descriptives-data.frame(n,mean,median,sd,total.NA,percent.missing) rownames(descriptives)-c( ) Case.Number-1:n Missing.Values-ifelse(NAs0,Missing Value, ) missing.value-data.frame(Case.Number,Missing.Values) missing.values-missing.value[ which(Missing.Values=='Missing Value'),] list(NUMERIC DATA,DESCRIPTIVES=t(descriptives),CASE # OF MISSING VALUES=missing.values[,1]) } else{ n-length(variable) NAs-is.na(variable) total.NA-sum(NAs) percent.missing-total.NA/n descriptives-data.frame(n,total.NA,percent.missing) rownames(descriptives)-c( ) Case.Number-1:n Missing.Values-ifelse(NAs0,Missing Value, ) missing.value-data.frame(Case.Number,Missing.Values) missing.values-missing.value[ which(Missing.Values=='Missing Value'),] list(CATEGORICAL DATA,DESCRIPTIVES=t(descriptives),CASE # OF MISSING VALUES=missing.values[,1]) } } dataset-data.frame(dataset) options(scipen=100) options(digits=2) lapply(dataset,find.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. David Winsemius, MD 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] kernel density plot
It is better to replace your later calls to plot with calls to lines instead, then you don't need to use par(new=T) which as you see tends to cause problems. -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of Muzna Alvi Sent: Sunday, April 03, 2011 4:56 AM To: r-help@r-project.org Subject: [R] kernel density plot I am using the following commands for plotting kernel density for three kinds of crops density(s22$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-t plot(t, xlim=c(-3,4), main=Net Income Distribution, axes=F, ylim=c(0,0.00035). xlab=Value in Rupees) par(new=T) density(s33$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-u plot(u, xlim=c(-3,4), axes=F, main=, col=red, ylim=c(0,0.00035)) par(new=T) density(s44$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-v plot(v, xlim=c(-3,4), col=blue, axes=F, main=, ylim=c(0,0.00035)) the problem is that in the graph that is drawn 1. the xlab gets hidden with the [N= and the bandwidth=] values 2. when i do par(new=T) this N and bandwidth value appears multiple times..overlapping each time and making the graph look untidy.. Is there any way of making these N and Bandwidth values not appear in the graph? Thanks -- -- [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] kernel density plot
i am sorry greg, can you explain that with an example? On Mon, Apr 4, 2011 at 12:11 AM, Greg Snow greg.s...@imail.org wrote: It is better to replace your later calls to plot with calls to lines instead, then you don't need to use par(new=T) which as you see tends to cause problems. -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of Muzna Alvi Sent: Sunday, April 03, 2011 4:56 AM To: r-help@r-project.org Subject: [R] kernel density plot I am using the following commands for plotting kernel density for three kinds of crops density(s22$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-t plot(t, xlim=c(-3,4), main=Net Income Distribution, axes=F, ylim=c(0,0.00035). xlab=Value in Rupees) par(new=T) density(s33$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-u plot(u, xlim=c(-3,4), axes=F, main=, col=red, ylim=c(0,0.00035)) par(new=T) density(s44$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-v plot(v, xlim=c(-3,4), col=blue, axes=F, main=, ylim=c(0,0.00035)) the problem is that in the graph that is drawn 1. the xlab gets hidden with the [N= and the bandwidth=] values 2. when i do par(new=T) this N and bandwidth value appears multiple times..overlapping each time and making the graph look untidy.. Is there any way of making these N and Bandwidth values not appear in the graph? Thanks -- -- [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. -- -- [[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] Help in splitting ists into sub-lists
Hi: here is one solution. Not so elegant but maybe it will give you some ideas: mylist1-rep(mylist,c(2,2)) a-matrix(c(index1,!index1,index2,!index2),ncol=4) mylist2-list() for (i in 1:4) { mylist2[[i]]-mylist1[[i]][,a[,i]] } mylist2 Andrija On Sun, Apr 3, 2011 at 6:56 PM, Axel Urbiz axel.ur...@gmail.com wrote: Dear List, Let's say I have a list whose components are 2 matrices (as exemplified in the mylist object below). I'd like to create a list with components being 4 matrices based on an logical index vector. is there a way to simplify what I'm doing to obtain the results in mylist2? I'd like something that would work on an arbitrary number of elements in mylist. mylist - list(matrix(1:9,3,3), matrix(10:18,3,3)) index1 - c(TRUE,FALSE,TRUE) index2- c(FALSE,TRUE,TRUE) mylist2 - list(mylist[[1]][,index1], mylist[[1]][,index1==FALSE], mylist[[2]][,index2], mylist[[2]][,index2==FALSE]) Thanks in advance, Axel. [[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. [[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] kernel density plot
It would be easier to test and example if we knew what s22, s33, and s44 were, but you can try: density(s22$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-t plot(t, xlim=c(-3,4), main=Net Income Distribution, axes=F, ylim=c(0,0.00035). xlab=Value in Rupees) density(s33$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-u lines(u, col=red) density(s44$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-v lines(v, col=blue) -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 From: Muzna Alvi [mailto:muzna.a...@gmail.com] Sent: Sunday, April 03, 2011 12:52 PM To: Greg Snow Cc: r-help@r-project.org Subject: Re: [R] kernel density plot i am sorry greg, can you explain that with an example? On Mon, Apr 4, 2011 at 12:11 AM, Greg Snow greg.s...@imail.orgmailto:greg.s...@imail.org wrote: It is better to replace your later calls to plot with calls to lines instead, then you don't need to use par(new=T) which as you see tends to cause problems. -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.orgmailto:greg.s...@imail.org 801.408.8111 -Original Message- From: r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.org [mailto:r-help-bounces@r-mailto:r-help-bounces@r- project.orghttp://project.org] On Behalf Of Muzna Alvi Sent: Sunday, April 03, 2011 4:56 AM To: r-help@r-project.orgmailto:r-help@r-project.org Subject: [R] kernel density plot I am using the following commands for plotting kernel density for three kinds of crops density(s22$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-t plot(t, xlim=c(-3,4), main=Net Income Distribution, axes=F, ylim=c(0,0.00035). xlab=Value in Rupees) par(new=T) density(s33$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-u plot(u, xlim=c(-3,4), axes=F, main=, col=red, ylim=c(0,0.00035)) par(new=T) density(s44$Net_income_Total.1, bw=nrd0,adjust=1, kernel=c(gaussian))-v plot(v, xlim=c(-3,4), col=blue, axes=F, main=, ylim=c(0,0.00035)) the problem is that in the graph that is drawn 1. the xlab gets hidden with the [N= and the bandwidth=] values 2. when i do par(new=T) this N and bandwidth value appears multiple times..overlapping each time and making the graph look untidy.. Is there any way of making these N and Bandwidth values not appear in the graph? Thanks -- -- [[alternative HTML version deleted]] __ R-help@r-project.orgmailto:R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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.
[R] setCoefTemplate
Hi everyone, I am trying to build a table putting standard errors horizontally. I haven't been able to do it. library(memisc) berkeley - aggregate(Table(Admit,Freq)~.,data=UCBAdmissions) berk0 - glm(cbind(Admitted,Rejected)~1,data=berkeley,family=binomial) berk1 - glm(cbind(Admitted,Rejected)~Gender,data=berkeley,family=binomial) berk2 - glm(cbind(Admitted,Rejected)~Gender+Dept,data=berkeley,family=binomial) setCoefTemplate(est.se=c(est = ($est:#)($se:#))) mtable(berk0,berk1,berk2, + coef.style=est.se, + summary.stats=c(Deviance,AIC,N)) Error in dim(ans) - newdims : dims [product 1] do not match the length of object [2] Thank you in advance. -- Sebastián Daza sebastian.d...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] zoo:rollapply by multiple grouping factors
# Hi there, # I am trying to apply a function over a moving-window for a large number of multivariate time-series that are grouped in a nested set of factors. I have spent a few days searching for solutions with no luck, so any suggestions are much appreciated. # The data I have are for the abundance dynamics of multiple species observed in multiple fixed plots at multiple sites. (I total I have 7 sites, ~3-5 plots/site, ~150 species/plot, for 60 time-steps each.) So my data look something like this: dat-data.frame(Site=rep(1), Plot=rep(c(rep(1,8),rep(2,8),rep(3,8)),1), Time=rep(c(1,1,2,2,3,3,4,4)), Sp=rep(1:2), Count=sample(24)) dat # Let the function I want to apply over a right-aligned window of w=2 time steps be: cv-function(x){sd(x)/mean(x)} w-2 # The final output I want would look something like this: Out-data.frame(dat,CV=round(c(NA,NA,runif(6,0,1),c(NA,NA,runif(6,0,1))),2)) # I could reshape and apply zoo:rollapply() to a given plot at a given site, and reshape again as follows: library(zoo) a-subset(dat,Site==1Plot==1) b-reshape(a[-c(1,2)],v.names='Count',idvar='Time',timevar='Sp',direction='wide') d-zoo(b[,-1],b[,1]) d out-rollapply(d, w, cv, na.pad=T, align='right') out # I would thereby have to loop through all my sites and plots which, although it deals with all species at once, still seems exceedingly inefficient. # So the question is, how do I use something like aggregate.zoo or tapply or even lapply to apply rollapply on each species' time series. # The closest I've come is the following two approaches: # First let: datx-list(Site=dat$Site,Plot=dat$Plot,Sp=dat$Sp) daty-dat$Count # Method 1. out1-tapply(seq(along=daty),datx,function(i,x=daty){ rollapply(zoo(x[i]), w, cv, na.pad=T, align='right') }) out1 out1[,,1] # Which works in that it gives me the right answers, but in a format from which I can't figure out how to get back into the format I want. # Method 2. fun-function(x){y-zoo(x);coredata(rollapply(y, w, cv,na.pad=T,align='right'))} out2-aggregate(daty,by=datx,fun) out2 # Which superficially works better, but again only in a format I can't figure out how to use because the output seems to be a mix of data.frame and lists. out2[1,4] out2[1,5] is.data.frame(out2) is.list(out2) # The situation is made more problematic by the fact that the time point of first survey can differ between plots (e.g., site1-plot3 may only start at time-point 3). As in... dat2-dat dat2-dat2[-which(dat2$Plot==3 dat2$Time3),] dat2 # I must therefore ensure that I'm keeping track of the true time associated with each value, not just the order of their occurences. This information is (seemingly) lost by both methods. datx-list(Site=dat2$Site,Plot=dat2$Plot,Sp=dat2$Sp) daty-dat2$Count # Method 1. out3-tapply(seq(along=daty),datx,function(i,x=daty){ rollapply(zoo(x[i]), w, cv, na.pad=T, align='right') }) out3 out3[1,3,1] time(out3[1,3,1]) # Method 2 out4-aggregate(daty,by=datx,fun) out4 time(out4[3,4]) # Am I going about this all wrong? Is there a different package to try? Any thoughts and suggestions are much appreciated! # R 2.12.2 GUI 1.36 Leopard build 32-bit (5691); zoo 1.6-4 # Thanks! # -mark -- ~~-- Ecology Evolutionary Biology University of California, Santa Cruz Long Marine Laboratory 100 Shaffer Road Santa Cruz, CA 95060-5730 Ph: 773-256-8645 Fax: 831-459-3383 http://people.ucsc.edu/~mnovak1/ ~~-- __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] another question on shapefiles and geom_point in ggplot2
Manuel: I changed your variable names from x to 'long' and y to 'lat' on the riqueza_out.csv file. The code below should do what you want. Also, since the legend title is kind of long, I broke it down into three lines so you can see more plot area. I am cc'ing the other groups so more people use it if needed. library(rgdal) library(ggplot2) library(sp) library(maptools) gpclibPermit() manuel - readOGR(dsn=., layer=AI_BIOTICA_010411_CRTM05) names(manuel);dim(manuel) slotNames(manuel) # look at the slot names # add the 'id' variable to the shapefile and use it to merge both files manuel@data$id = rownames(manuel@data) # convert shapefile to dataframe manuel.df - as.data.frame(manuel) # fortify to plot with ggplot2 manuel_fort - fortify(manuel,region=id) head(manuel_fort) # Merge shapefile and the as.dataframe shapefile manuel_merged - join(manuel_fort,manuel.df, by =id) head(manuel_merged) # Read in the csv file manuel_points - read.csv(riqueza_out.csv) head(manuel_points);dim(manuel_points) # fortify this one too for the points or else an error will ocurr manuel_points - fortify(manuel_points) manuel_points # Plot the shapefile and overlayed the points over it p - ggplot(manuel_merged, aes(long,lat,group=group)) + geom_polygon(aes(data=manuel_merged,fill=Area_Influ)) + geom_path(color=white) + theme_bw() # remove this if you don't want black and white background p + geom_point(data=manuel_points,aes(size=ACE,colour=ACE,group=NULL)) + scale_size(name = Número\nde\nespecies, breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20)) + scale_colour_gradientn(name = 'Número\nde\nespecies', colours = rainbow(6), breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20))+ xlab(Longitud) + ylab(Latitud) + opts(axis.text.x = theme_text(size = 8, vjust = 1)) + opts(axis.text.y = theme_text(size = 8, hjust = 1)) Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx From: Manuel Spínola mspinol...@gmail.com To: Felipe Carrillo mazatlanmex...@yahoo.com Sent: Sat, April 2, 2011 11:22:24 PM Subject: Re: another question on shapefiles and geom_point in ggplot2 No problem, thank you very much Felipe. Best, Manuel On 03/04/2011 12:19 a.m., Felipe Carrillo wrote: I meant to send you this one..Let me clean up the code a little bit and I will send it to you,,,do you mind if I send it to you in the morning? Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx From: Manuel Spínola mspinol...@gmail.com To: Felipe Carrillo mazatlanmex...@yahoo.com Sent: Sat, April 2, 2011 11:15:28 PM Subject: Re: another question on shapefiles and geom_point in ggplot2 Yes Felipe. That is the graph I was looking for. I got something closer but no like yours. How did you do it? Manuel On 03/04/2011 12:10 a.m., Felipe Carrillo wrote: I was able to open them,,I am attaching a picture of the graph I created..It's that what you had in mind? Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx From: Manuel Spínola mspinol...@gmail.com To: Felipe Carrillo mazatlanmex...@yahoo.com Sent: Sat, April 2, 2011 10:35:51 PM Subject: Re: another question on shapefiles and geom_point in ggplot2 It should be. I am sending them again. Manuel On 02/04/2011 10:23 p.m., Felipe Carrillo wrote: Manuel: I can't open the shapefile, is this the original one? Is the csv file the one that you are trying to overlay on top of the shapefile? Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx From: Manuel Spínola mspinol...@gmail.com To: Felipe Carrillo mazatlanmex...@yahoo.com Sent: Sat, April 2, 2011 6:14:09 PM Subject: Re: another question on shapefiles and geom_point in ggplot2 Files attached. On 02/04/2011 07:04 p.m., Felipe Carrillo wrote: If you want individual points overlayed on the shapefile, you need to add another variable to it before you fortify it. After you fortify merge both the fortified dataset and the original shapefile. Go ahead and post your shapefile to see if I can figure it out. Do you just want the points or want text also? Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx From: Manuel Spínola mspinol...@gmail.com To: Felipe Carrillo mazatlanmex...@yahoo.com Sent: Sat, April 2, 2011 5:24:02 PM Subject: Re: another question on shapefiles and geom_point in ggplot2 Hi Felipe, I did the same thing that I am trying know, attached is how it looks. Best, Manuel On 02/04/2011 06:09 p.m., Felipe Carrillo wrote: Manuel: I did something
Re: [R] R gui on windows how to force to always show the last line of output
i wonder why you want to change R but not your program itself, you just print the last sentence is ok? wuleiwong statistic@CSU On 03-Apr-2011 9:06 PM, lcn lcn...@gmail.com wrote: CTRL+L helps. 2011/4/2 stan zimine szim...@gmail.com Hi. Googled but did not found the answer for the following little issue. how to force R gui on windows (maybe a specific setting) to always show the last line of output in the window console. My program in R makes measurements every 5 mins in indefinite loop and prints results in the console. The problem: last messages are not visible, The scrolling bar of the gui console gets shorter. I.e. you have to scroll for the last messages. Thanks if anybody knows the sol to this prob. SZ __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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. [[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] Discretizing data rows into regular intervals
i suggest you'd better translate the.data format to another form like in SAS, the data can be a number count from a fixed time like 2005-01-01, then you can translate the datas to a sequence. On 03-Apr-2011 5:40 PM, Gabor Grothendieck ggrothendi...@gmail.com wrote: On Sat, Apr 2, 2011 at 9:31 PM, Linh Tran tra...@berkeley.edu wrote: Hi guys, I'd like to thank you ahead of time for any help that you can offer me. I'm kind of stuck trying to do this. I have a data frame with dates and values (note: only two columns shown): head(test) date value stop 1 01/02/05 100 12/01/07 2 07/16/05 200 12/01/07 3 12/20/05 150 12/01/07 4 04/01/06 250 12/01/07 5 10/01/06 10 12/01/07 What I need to do is create regularly spaced 3-month intervals (starting with the first observed date) with values that are closest to but recorded after the date created. I would stop at the stop date. So the result would look like: new_date value 1 01/02/05 100 2 04/02/05 100 3 07/02/05 100 4 10/02/05 200 5 01/02/06 150 6 04/02/06 250 7 07/02/06 250 8 10/02/06 10 9 01/02/07 10 etc etc etc 10/02/07 --- ## Final obs since next one would be 1/2/08 (after stop date) See question #13 in the zoo-faq vignette: http://cran.r-project.org/web/packages/zoo/index.html and note the existence of zoo's yearqtr class. -- 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. [[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] setCoefTemplate
Hi everyone, I am trying to build a table putting standard errors horizontally. I haven't been able to do it. library(memisc) berkeley - aggregate(Table(Admit,Freq)~.,data=UCBAdmissions) berk0 - glm(cbind(Admitted,Rejected)~1,data=berkeley,family=binomial) berk1 - glm(cbind(Admitted,Rejected)~Gender,data=berkeley,family=binomial) berk2 - glm(cbind(Admitted,Rejected)~Gender+Dept,data=berkeley,family=binomial) setCoefTemplate(est.se=c(est = ($est:#)($se:#))) mtable(berk0,berk1,berk2, + coef.style=est.se, + summary.stats=c(Deviance,AIC,N)) Error in dim(ans) - newdims : dims [product 1] do not match the length of object [2] Thank you in advance. -- Sebastián Daza sebastian.d...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] coefficients style
Hi everyone, I am trying to build a table putting standard errors horizontally. I haven't been able to do it. library(memisc) berkeley - aggregate(Table(Admit,Freq)~.,data=UCBAdmissions) berk0 - glm(cbind(Admitted,Rejected)~1,data=berkeley,family=binomial) berk1 - glm(cbind(Admitted,Rejected)~Gender,data=berkeley,family=binomial) berk2 - glm(cbind(Admitted,Rejected)~Gender+Dept,data=berkeley,family=binomial) setCoefTemplate(est.se=c(est = ($est:#)($se:#))) mtable(berk0,berk1,berk2, + coef.style=est.se, + summary.stats=c(Deviance,AIC,N)) Error in dim(ans) - newdims : dims [product 1] do not match the length of object [2] Thank you in advance. -- Sebastián Daza sebastian.d...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] style question
Hi everyone, I am trying to build a table putting standard errors horizontally. I haven't been able to do it. library(memisc) berkeley - aggregate(Table(Admit,Freq)~.,data=UCBAdmissions) berk0 - glm(cbind(Admitted,Rejected)~1,data=berkeley,family=binomial) berk1 - glm(cbind(Admitted,Rejected)~Gender,data=berkeley,family=binomial) berk2 - glm(cbind(Admitted,Rejected)~Gender+Dept,data=berkeley,family=binomial) setCoefTemplate(est.se=c(est = ($est:#)($se:#))) mtable(berk0,berk1,berk2, + coef.style=est.se, + summary.stats=c(Deviance,AIC,N)) Error in dim(ans) - newdims : dims [product 1] do not match the length of object [2] Thank you in advance. -- Sebastián Daza sebastian.d...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Homals package color function problem
Hello The Homals package and its plot options are excellent. However, I am unable to manipulate the colour in the plots. In a call such as: plot(mc_analysis, plot.dim = c(1,3), plot.type = jointplot, col = 1) this should be straightforward - but I can't seem to effect the plotted colours I have tried various combinations for col commands in other plot packages I know col = col = c('red', 'green') ccol = rcol = also the different ways of designating the colour col = 1, col = black, col = (#6698FF) I get no error messages, I have tried flushing R, restarting etc. the colours never change I'm under MacOS- 10.5.8, in R- 2.12.2 I'm sure it is something so obvious I'll crawl away and hide when I find it, but for now ideas anyone? Dylan -- View this message in context: http://r.789695.n4.nabble.com/Homals-package-color-function-problem-tp3423814p3423814.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] installing ncdf on Ubuntu 10.04
Hello I'm trying to install the ncdf package on a Ubuntu 10.04 laptop. Can you offer suggestions to install this package? install.packages(ncdf, repositories = TRUE) Installing package(s) into /home/steve/R/i486-pc-linux-gnu-library/2.12 (as lib is unspecified) --- Please select a CRAN mirror for use in this session --- Loading Tcl/Tk interface ... done Error in download.file(url, destfile, method, mode = wb, ...) : unused argument(s) (repositories = TRUE) Warning in download.packages(pkgs, destdir = tmpd, available = available, : download of package 'ncdf' failed install.packages(ncdf) Installing package(s) into /home/steve/R/i486-pc-linux-gnu-library/2.12 (as lib is unspecified) trying URL 'http://cran.mirrors.hoobly.com/src/contrib/ncdf_1.6.5.tar.gz' Content type 'application/x-gzip' length 78738 bytes (76 Kb) opened URL == downloaded 76 Kb * installing *source* package ncdf ... checking for nc-config... no checking for gcc... gcc checking whether the C compiler works... yes checking for C compiler default output file name... a.out checking for suffix of executables... checking whether we are cross compiling... no checking for suffix of object files... o checking whether we are using the GNU C compiler... yes checking whether gcc accepts -g... yes checking for gcc option to accept ISO C89... none needed checking how to run the C preprocessor... gcc -E checking for grep that handles long lines and -e... /bin/grep checking for egrep... /bin/grep -E checking for ANSI C header files... no checking for sys/types.h... no checking for sys/stat.h... no checking for stdlib.h... no checking for string.h... no checking for memory.h... no checking for strings.h... no checking for inttypes.h... no checking for stdint.h... no checking for unistd.h... no checking netcdf.h usability... no checking netcdf.h presence... no checking for netcdf.h... no configure: error: netcdf header netcdf.h not found ERROR: configuration failed for package ncdf * removing /home/steve/R/i486-pc-linux-gnu-library/2.12/ncdf The downloaded packages are in /tmp/Rtmpt709vH/downloaded_packages Warning message: In install.packages(ncdf) : installation of package 'ncdf' had non-zero exit status Thank you Steve [[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] Standard Error for Cointegration Results
Dear Sir/Madam, I have used ca.jo in urca package to identify the cointegration and cajorls to estimate the vecm. Althought both return the coefficients for long run relationship (or ect1 in cajorls), I am unable to find the standard error and t statistics. I spend some weeks to search around. I did find some similar enquiries before and answer provided Prof. Pfaff is to use vec2var. However, I still can't find corresponding info. after using vec2var function. Highly appreciate if someone could advise by using a step by step with codes. For example, in the paper of VAR, SVAR and SVEC Models: Implementation Within R Package vars, how to get the t-statistics / standard error for the beta in table 5, p.25? Thank you in advance for your effort in addressing my problem. Yours, Aries [[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] Structural equation modeling in R(lavaan,sem)
Daer all, I have a question concerning longitudinal data: When we have a longitudinal data and we have to do sem analysis there is in the package lavaan some functions,options in this package that help to do this or we can treat these data like non longitudinal data Thanks you a lot Antra EL MOUSSELLY Date: Tue, 29 Mar 2011 13:31:32 -0700 From: ml-node+3416220-556418632-225...@n4.nabble.com To: antr...@hotmail.com Subject: Re: Structural equation modeling in R(lavaan,sem) sem (the package) documentation is not intended to teach you how to do SEM (the technique) (there's very little R documentation that is intended to teach you how to do a particular statistical technique). There are several good books out there, but here's a free access journal article, which will help. http://www.biomedcentral.com/1756-0500/3/267 Might I also suggest you take a look at the semnet list, which is populated by practitioners of SEM. Jeremy On 29 March 2011 12:25, jouba [hidden email] wrote: Dear all, There is some where documentation to understand all indices in the output of the function sem(package lavaan ) ?? for example Chi-square test baseline model, Full model versus baseline model, Loglikelihood and Information Criteria, Root Mean Square Error of Approximation, Standardized Root Mean Square Residual⦠Th same question for the sem funtion (sem package) Thanks in advance for your help -- View this message in context: http://r.789695.n4.nabble.com/Structural-equation-modeling-in-R-lavaan-sem-tp3409642p3415954.html Sent from the R help mailing list archive at Nabble.com. __ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jeremy Miles Psychology Research Methods Wiki: www.researchmethodsinpsychology.com __ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. If you reply to this email, your message will be added to the discussion below:http://r.789695.n4.nabble.com/Structural-equation-modeling-in-R-lavaan-sem-tp3409642p3416220.html To unsubscribe from Structural equation modeling in R(lavaan,sem), click here. -- View this message in context: http://r.789695.n4.nabble.com/Structural-equation-modeling-in-R-lavaan-sem-tp3409642p3424029.html Sent from the R help mailing list archive at Nabble.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.
[R] R-project: plot 2 zoo objects (price series) that have some date mis-matches
I have 2 zoo objects - 1) Interest rate spread between 10-YR-US-Treasury and 2-YR-US-Treasury (object name = sprd) 2) SP 500 index (object name = spy) str(spy) ‘zoo’ series from 1976-06-01 to 2011-03-31 Data: num [1:8791] 99.8 100.2 100.1 99.2 98.6 ... Index: Class 'Date' num [1:8791] 2343 2344 2345 2346 2349 ... str(sprd) ‘zoo’ series from 1976-06-01 to 2011-03-31 Data: num [1:9088] 0.68 0.71 0.7 0.77 0.79 0.79 0.82 0.86 0.83 0.83 ... Index: Class 'Date' num [1:9088] 2343 2344 2345 2346 2349 ... Since there are NA data points in object 'sprd' I created another object that omits NA. The name of that object is 'sprdtmp'. str(sprdtmp) ‘zoo’ series from 1976-06-01 to 2011-03-31 Data: atomic [1:8704] 0.68 0.71 0.7 0.77 0.79 0.79 0.82 0.86 0.83 0.83 ... - attr(*, na.action)=Class 'omit' int [1:384] 25 70 95 111 118 128 149 190 224 260 ... Index: Class 'Date' num [1:8704] 2343 2344 2345 2346 2349 ... I want to plot both time series on the same plot with time/date on the x-axis and the axis label quarterly (or monthly). One problem is that the sprdtmp and spy objects do not have the same number of data points as there are times when the equity markets are closed while the interest rate markets are open. For the most part the dates overlap. Would this matter if I try to plot both on the same plot? And how would I go about plotting these objects in one plot. The second part of the plot requires that both are on different scales. I guess I could take a log of the sp and plot that along with the rate spread. But it would be nice for future reference how I could plot 2 series in one plot with 2 difference scales. I spent all night yesterday and this morning trying various options but I cant seem to get this to work. I have read through the ?plot, ?plot.zoo and ?axis documentation without any success. I would greatly appreciate if anyone can point me in the right direction. Thank you kindly. -- View this message in context: http://r.789695.n4.nabble.com/R-project-plot-2-zoo-objects-price-series-that-have-some-date-mis-matches-tp3423899p3423899.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] Function for finding NA's
aThanks David, After seeing the simplicity of your function versus the convoluted mess I worked up I now understand why it's not necessary to have a package to find NA's (and from what you said is a part of other packages such as Hmisc already). I am at the 2 1/2 month mark as an R user and have loads to learn. Simpler is better. Thanks David for your time and I will take the information you gave and put it to use in new situations. Tyler CC: r-help@r-project.org From: dwinsem...@comcast.net To: tyler_rin...@hotmail.com Subject: Re: [R] Function for finding NA's Date: Sun, 3 Apr 2011 14:19:40 -0400 On Apr 3, 2011, at 1:44 PM, Tyler Rinker wrote: Quick question, I tried to find a function in available packages to find NA's for an entire data set (or single variables) and report the row of missing values (NA's for each column). I searched the typical routes through the blogs and the help manuals for 15 minutes. Rather than spend any more time searching I created my own function to do this (probably in less time than it would have taken me to find the function). Now I still have the same question: Is this function (NAhunter I call it) already in existence? If so please direct me (because I'm sure they've written better code more efficiently). I highly doubt I'm this first person to want to find all the missing values in a data set so I assume there is a function for it but I just didn't spend enough time looking. If there is no existing function (big if here), is this something people feel is worthwhile for me to put into a package of some sort? I'm not sure that it would have occurred to people to include it in a package. Consider: getNa - function(dfrm) lapply(dfrm, function(x) which(is.na(x) ) ) cities long lat city pop 1 -58.38194 -34.59972 Buenos Aires NA 2 14.25000 40.8 NA NA getNa(cities) $long integer(0) $lat integer(0) $city [1] 2 $pop [1] 1 2 There are several packages with functions by the name `describe` that do most or all of rest of what you have proposed. I happen to use Harrell's Hmisc but the other versions should also be reviewed if you want to avoid re-inventing the wheel. -- David. Tyler Here's the code: NAhunter-function(dataset) { find.NA-function(variable) { if(is.numeric(variable)){ n-length(variable) mean-mean(variable, na.rm=T) median-median(variable, na.rm=T) sd-sd(variable, na.rm=T) NAs-is.na(variable) total.NA-sum(NAs) percent.missing-total.NA/n descriptives-data.frame(n,mean,median,sd,total.NA,percent.missing) rownames(descriptives)-c( ) Case.Number-1:n Missing.Values-ifelse(NAs0,Missing Value, ) missing.value-data.frame(Case.Number,Missing.Values) missing.values-missing.value[ which(Missing.Values=='Missing Value'),] list(NUMERIC DATA,DESCRIPTIVES=t(descriptives),CASE # OF MISSING VALUES=missing.values[,1]) } else{ n-length(variable) NAs-is.na(variable) total.NA-sum(NAs) percent.missing-total.NA/n descriptives-data.frame(n,total.NA,percent.missing) rownames(descriptives)-c( ) Case.Number-1:n Missing.Values-ifelse(NAs0,Missing Value, ) missing.value-data.frame(Case.Number,Missing.Values) missing.values-missing.value[ which(Missing.Values=='Missing Value'),] list(CATEGORICAL DATA,DESCRIPTIVES=t(descriptives),CASE # OF MISSING VALUES=missing.values[,1]) } } dataset-data.frame(dataset) options(scipen=100) options(digits=2) lapply(dataset,find.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. David Winsemius, MD West Hartford, CT [[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 do I modify uniroot function to return .0001 if error ?
I am calling the uniroot function from inside another function using these lines (last two lines of the function) : d - uniroot(k, c(.001, 250), tol=.05) return(d$root) The problem is that on occasion there's a problem with the values I'm passing to uniroot. In those instances uniroot stops and sends a message that it can't calculate the root because f.upper * f.lower is greater than zero. All I'd like to do in those cases is be able to set the return value of my calling function return(d$root) to .0001. But I'm not sure how to pull that off. I tried a few modifications to uniroot but so far no luck. For convenience, the uniroot function is shown below: uniroot - function (f, interval, ..., lower = min(interval), upper = max(interval), f.lower = f(lower, ...), f.upper = f(upper, ...), tol = .Machine$double.eps^0.25, maxiter = 1000) { if (!missing(interval) length(interval) != 2L) stop('interval' must be a vector of length 2) if (!is.numeric(lower) || !is.numeric(upper) || lower = upper) stop(lower upper is not fulfilled) if (is.na(f.lower)) stop(f.lower = f(lower) is NA) if (is.na(f.upper)) stop(f.upper = f(upper) is NA) if (f.lower * f.upper 0) stop(f.up * f.down 0) val - .Internal(zeroin2(function(arg) f(arg, ...), lower, upper, f.lower, f.upper, tol, as.integer(maxiter))) iter - as.integer(val[2L]) if (iter 0) { warning(_NOT_ converged in , maxiter, iterations) iter - maxiter } list(root = val[1L], f.root = f(val[1L], ...), iter = iter, estim.prec = val[3L]) } -- View this message in context: http://r.789695.n4.nabble.com/How-do-I-modify-uniroot-function-to-return-0001-if-error-tp3424092p3424092.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] Plotting MDS (multidimensional scaling)
The header metric mds was actually a leftover because I initially used cmdscale and did not bother changing it for the example. Thanks, Daniel -- View this message in context: http://r.789695.n4.nabble.com/Plotting-MDS-multidimensional-scaling-tp3422670p3424135.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] installing ncdf on Ubuntu 10.04
On Mon, Apr 4, 2011 at 6:35 AM, Steve Friedman friedman.st...@gmail.com wrote: Hello I'm trying to install the ncdf package on a Ubuntu 10.04 laptop. Can you offer suggestions to install this package? snip checking for netcdf.h... no configure: error: netcdf header netcdf.h not found ERROR: configuration failed for package ‘ncdf’ * removing ‘/home/steve/R/i486-pc-linux-gnu-library/2.12/ncdf’ Looks like you need the netCDF libraries. https://launchpad.net/ubuntu/+source/netcdf suggests that you need both the binary and a -dev package with headers. The ncdf binaries on Windows and Mac have been specially built to include the netCDF libraries, but the source package doesn't include them. -thomas -- Thomas Lumley Professor of Biostatistics University of Auckland __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Function for finding NA's
On Apr 3, 2011, at 3:46 PM, Tyler Rinker wrote: aThanks David, After seeing the simplicity of your function versus the convoluted mess I worked up I now understand why it's not necessary to have a package to find NA's (and from what you said is a part of other packages such as Hmisc already). I'm actually not aware that any of the `describe` variants will return the indices of NA's. In the case of real dataset such an object could be fairly large. It was the other descriptive functions that I said were probably already coded. I am at the 2 1/2 month mark as an R user and have loads to learn. Simpler is better. Thanks David for your time and I will take the information you gave and put it to use in new situations. You should also familiarize yourself with complete.cases() and the various functions that handle na.action parameters (linked from that help page). Note that complete.cases returns a logical vector (not the cases themselves) and is designed for indexing matrices or dataframes. Tyler CC: r-help@r-project.org From: dwinsem...@comcast.net To: tyler_rin...@hotmail.com Subject: Re: [R] Function for finding NA's Date: Sun, 3 Apr 2011 14:19:40 -0400 On Apr 3, 2011, at 1:44 PM, Tyler Rinker wrote: Quick question, I tried to find a function in available packages to find NA's for an entire data set (or single variables) and report the row of missing values (NA's for each column). I searched the typical routes through the blogs and the help manuals for 15 minutes. Rather than spend any more time searching I created my own function to do this (probably in less time than it would have taken me to find the function). Now I still have the same question: Is this function (NAhunter I call it) already in existence? If so please direct me (because I'm sure they've written better code more efficiently). I highly doubt I'm this first person to want to find all the missing values in a data set so I assume there is a function for it but I just didn't spend enough time looking. If there is no existing function (big if here), is this something people feel is worthwhile for me to put into a package of some sort? I'm not sure that it would have occurred to people to include it in a package. Consider: getNa - function(dfrm) lapply(dfrm, function(x) which(is.na(x) ) ) cities long lat city pop 1 -58.38194 -34.59972 Buenos Aires NA 2 14.25000 40.8 NA NA getNa(cities) $long integer(0) $lat integer(0) $city [1] 2 $pop [1] 1 2 There are several packages with functions by the name `describe` that do most or all of rest of what you have proposed. I happen to use Harrell's Hmisc but the other versions should also be reviewed if you want to avoid re-inventing the wheel. -- David. Tyler Here's the code: NAhunter-function(dataset) { find.NA-function(variable) { if(is.numeric(variable)){ n-length(variable) mean-mean(variable, na.rm=T) median-median(variable, na.rm=T) sd-sd(variable, na.rm=T) NAs-is.na(variable) total.NA-sum(NAs) percent.missing-total.NA/n descriptives- data.frame(n,mean,median,sd,total.NA,percent.missing) rownames(descriptives)-c( ) Case.Number-1:n Missing.Values-ifelse(NAs0,Missing Value, ) missing.value-data.frame(Case.Number,Missing.Values) missing.values-missing.value[ which(Missing.Values=='Missing Value'),] list(NUMERIC DATA,DESCRIPTIVES=t(descriptives),CASE # OF MISSING VALUES=missing.values[,1]) } else{ n-length(variable) NAs-is.na(variable) total.NA-sum(NAs) percent.missing-total.NA/n descriptives-data.frame(n,total.NA,percent.missing) rownames(descriptives)-c( ) Case.Number-1:n Missing.Values-ifelse(NAs0,Missing Value, ) missing.value-data.frame(Case.Number,Missing.Values) missing.values-missing.value[ which(Missing.Values=='Missing Value'),] list(CATEGORICAL DATA,DESCRIPTIVES=t(descriptives),CASE # OF MISSING VALUES=missing.values[,1]) } } dataset-data.frame(dataset) options(scipen=100) options(digits=2) lapply(dataset,find.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. David Winsemius, MD West Hartford, CT David Winsemius, MD 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] Quick question about duplicating vectors
Hello R users, I have a quick question, if I have a data set a, say a - c(1,2,3,4,5) and I want to get a new data set b, b - c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5) Could anyone tell me how I can obtain this? Thank you! Xiaoxi [[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] another question on shapefiles and geom_point in ggplot2
Thank you very much Felipe, Did you see the solution from ahmadou dicko? He doesn´t use gpclibPermit() I have another option but I cannot get the right fill for the id. See attached map. ai_biotica = readOGR(dsn=C:/ProyectosRespacial/ICE/SIG_Biotica_PHED, layer=AI_BIOTICA_010411_CRTM05) str(ai_biotica) # fortify to get the data fortify.ai_biotica - fortify.SpatialPolygonsDataFrame(ai_biotica,region='Area_Influ') names(fortify.ai_biotica) str(fortify.ai_biotica) levels(fortify.ai_biotica$group) # mapa ggplot(fortify.ai_biotica, aes(x = long, y=lat, group = group)) + geom_polygon(colour = black, fill = NA) geo = read.csv(riqueza_out.csv, sep = ,, header = T) names(geo) str(geo) summary(geo) # mapa con riqueza p = ggplot(geo, aes(x, y)) p + geom_point(aes(size = ACE, colour = ACE)) + theme_bw() + scale_size(name = Número de especies, breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20)) + scale_colour_gradientn(name = 'Número de especies', colours = heat.colors(10), breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20))+ xlab(Longitud) + ylab(Latitud) + opts(axis.text.x = theme_text(size = 8, vjust = 1)) + opts(axis.text.y = theme_text(size = 8, hjust = 1)) + geom_path(aes(x=long,y=lat,group=group, fill=id),data=fortify.ai_biotica) Best, Manuel On 03/04/2011 01:41 p.m., Felipe Carrillo wrote: Manuel: I changed your variable names from x to 'long' and y to 'lat' on the riqueza_out.csv file. The code below should do what you want. Also, since the legend title is kind of long, I broke it down into three lines so you can see more plot area. I am cc'ing the other groups so more people use it if needed. library(rgdal) library(ggplot2) library(sp) library(maptools) gpclibPermit() manuel - readOGR(dsn=., layer=AI_BIOTICA_010411_CRTM05) names(manuel);dim(manuel) slotNames(manuel) # look at the slot names # add the 'id' variable to the shapefile and use it to merge both files manuel@data$id mailto:manuel@data$id = rownames(manuel@data mailto:manuel@data) # convert shapefile to dataframe manuel.df - as.data.frame(manuel) # fortify to plot with ggplot2 manuel_fort - fortify(manuel,region=id) head(manuel_fort) # Merge shapefile and the as.dataframe shapefile manuel_merged - join(manuel_fort,manuel.df, by =id) head(manuel_merged) # Read in the csv file manuel_points - read.csv(riqueza_out.csv) head(manuel_points);dim(manuel_points) # fortify this one too for the points or else an error will ocurr manuel_points - fortify(manuel_points) manuel_points # Plot the shapefile and overlayed the points over it p - ggplot(manuel_merged, aes(long,lat,group=group)) + geom_polygon(aes(data=manuel_merged,fill=Area_Influ)) + geom_path(color=white) + theme_bw() # remove this if you don't want black and white background p + geom_point(data=manuel_points,aes(size=ACE,colour=ACE,group=NULL)) + scale_size(name = Número\nde\nespecies, breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20)) + scale_colour_gradientn(name = 'Número\nde\nespecies', colours = rainbow(6), breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20))+ xlab(Longitud) + ylab(Latitud) + opts(axis.text.x = theme_text(size = 8, vjust = 1)) + opts(axis.text.y = theme_text(size = 8, hjust = 1)) Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx *From:* Manuel Spínola mspinol...@gmail.com *To:* Felipe Carrillo mazatlanmex...@yahoo.com *Sent:* Sat, April 2, 2011 11:22:24 PM *Subject:* Re: another question on shapefiles and geom_point in ggplot2 No problem, thank you very much Felipe. Best, Manuel On 03/04/2011 12:19 a.m., Felipe Carrillo wrote: I meant to send you this one..Let me clean up the code a little bit and I will send it to you,,,do you mind if I send it to you in the morning? Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx *From:* Manuel Spínola mspinol...@gmail.com *To:* Felipe Carrillo mazatlanmex...@yahoo.com *Sent:* Sat, April 2, 2011 11:15:28 PM *Subject:* Re: another question on shapefiles and geom_point in ggplot2 Yes Felipe. That is the graph I was looking for. I got something closer but no like yours. How did you do it? Manuel On 03/04/2011 12:10 a.m., Felipe Carrillo wrote: I was able to open them,,I am attaching a picture of the graph I created..It's that what you had in mind? Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx *From:* Manuel Spínola mspinol...@gmail.com *To:* Felipe Carrillo mazatlanmex...@yahoo.com
Re: [R] style question
On May 1, 2008, at 1:46 PM, Sebastián Daza wrote: Hi everyone, I am trying to build a table putting standard errors horizontally. I haven't been able to do it. library(memisc) berkeley - aggregate(Table(Admit,Freq)~.,data=UCBAdmissions) berk0 - glm(cbind(Admitted,Rejected)~1,data=berkeley,family=binomial) berk1 - glm(cbind(Admitted,Rejected)~Gender,data=berkeley,family=binomial) berk2 - glm(cbind(Admitted,Rejected)~Gender +Dept,data=berkeley,family=binomial) setCoefTemplate(est.se=c(est = ($est:#)($se:#))) I'm not a skilled user of that package but just looking at the value of the last four leaves of the list you created makes me wonder if you meant to do something like the `ci.se.horizontal` variant? $ci.se.horizontal est se est ($est:#) (($se:#)) ci [($lwr:#) ($upr:#)] $ci.p est p lwr upr ($est:#) (($p:#)) [($lwr:#) ($upr:#)] $ci.p.horizontal est se est ($est:#) (($p:#)) ci [($lwr:#) ($upr:#)] $est.se # Your addition doesn't really look like the others in the list est ($est:#)($se:#) This runs without error: mtable(berk0,berk1,berk2, coef.style=ci.se.horizontal, summary.stats=c(Deviance,AIC,N)) On the other hand maybe you wanted this, (note a two item list): tt- setCoefTemplate(est.se=list(est = ($est:#), se=($se:#))) tt['est.se'] $est.se est se ($est:#) ($se:#) mtable(berk0,berk1,berk2, + coef.style=est.se, + summary.stats=c(Deviance,AIC,N)) Calls: berk0: glm(formula = cbind(Admitted, Rejected) ~ 1, family = binomial, data = berkeley) berk1: glm(formula = cbind(Admitted, Rejected) ~ Gender, family = binomial, data = berkeley) berk2: glm(formula = cbind(Admitted, Rejected) ~ Gender + Dept, family = binomial, data = berkeley) === berk0berk1berk2 --- (Intercept) -0.457 -0.2200.582 0.0310.0390.069 Gender: Female/Male-0.6100.100 0.0640.081 Dept: B/A -0.043 0.110 Dept: C/A -1.263 0.107 Dept: D/A -1.295 0.106 Dept: E/A -1.739 0.126 Dept: F/A -3.306 0.170 --- Deviance 877.056 783.607 20.204 AIC 947.996 856.547 103.144 N4526 4526 4526 === mtable(berk0,berk1,berk2, + coef.style=est.se, + summary.stats=c(Deviance,AIC,N)) Error in dim(ans) - newdims : dims [product 1] do not match the length of object [2] Thank you in advance. -- Sebastián Daza sebastian.d...@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. David Winsemius, MD 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] Avoiding loops in creating a coinvestment matrix
Hi, I am working on a dataset in which a number of venture capitalists invest in a number of firms. What I am creating is an asymmetric matrix M in which m(ij) is the volume (sum) of coinvestments of VC i with VC j (i.e., how much has VC i invested in companies that VC j also has investments in). The output should look like the coinvestments matrix produced with the code below. If possible I would like to avoid loops and optimize the code for speed because the real data is huge. If anybody has suggestions, I would be grateful. invest=c(20,50,40,30,10,20,20,30,40) vc=rep(c('A','B','C'),each=3) company=c('E','F','G','F','G','H','G','H','I') data=data.frame(vc,company,invest) data #data inv.mat=tapply(invest,list(vc,company),sum) inv.mat=replace(inv.mat,which(is.na(inv.mat)==T),0) inv.mat #investment matrix exist.mat=inv.mat0 coinvestments-matrix(0,nrow=length(unique(vc)),ncol=length(unique(vc))) for(i in unique(vc)){ for(j in unique(vc)){ i.is=which(unique(vc)==i) j.is=which(unique(vc)==j) i.invests=exist.mat[i,] j.invests=exist.mat[j,] which.i=which(i.invests==T) which.j=which(j.invests==T) i.invests.with.j=which.i[which.i%in%which.j] coinvestments[i.is,j.is]=sum(inv.mat[i.is,i.invests.with.j]) } } coinvestments system.time( for(i in unique(vc)){ for(j in unique(vc)){ i.is=which(unique(vc)==i) j.is=which(unique(vc)==j) i.invests=exist.mat[i,] j.invests=exist.mat[j,] which.i=which(i.invests==T) which.j=which(j.invests==T) i.invests.with.j=which.i[which.i%in%which.j] coinvestments[i.is,j.is]=sum(inv.mat[i.is,i.invests.with.j]) } } ) Thanks much, Daniel -- View this message in context: http://r.789695.n4.nabble.com/Avoiding-loops-in-creating-a-coinvestment-matrix-tp3424298p3424298.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] Quick question about duplicating vectors
On Apr 3, 2011, at 6:04 PM, Xiaoxi Gao wrote: Hello R users, I have a quick question, if I have a data set a, say a - c(1,2,3,4,5) You might as well learn to use more precise R terminology... that is a `vector`. and I want to get a new data set b, b - c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5) Could anyone tell me how I can obtain this? ?rep -- David Winsemius, MD 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] zoo:rollapply by multiple grouping factors
On Sun, Apr 3, 2011 at 11:58 AM, Mark Novak mnov...@ucsc.edu wrote: # Hi there, # I am trying to apply a function over a moving-window for a large number of multivariate time-series that are grouped in a nested set of factors. I have spent a few days searching for solutions with no luck, so any suggestions are much appreciated. # The data I have are for the abundance dynamics of multiple species observed in multiple fixed plots at multiple sites. (I total I have 7 sites, ~3-5 plots/site, ~150 species/plot, for 60 time-steps each.) So my data look something like this: dat-data.frame(Site=rep(1), Plot=rep(c(rep(1,8),rep(2,8),rep(3,8)),1), Time=rep(c(1,1,2,2,3,3,4,4)), Sp=rep(1:2), Count=sample(24)) dat # Let the function I want to apply over a right-aligned window of w=2 time steps be: cv-function(x){sd(x)/mean(x)} w-2 # The final output I want would look something like this: Out-data.frame(dat,CV=round(c(NA,NA,runif(6,0,1),c(NA,NA,runif(6,0,1))),2)) # I could reshape and apply zoo:rollapply() to a given plot at a given site, and reshape again as follows: library(zoo) a-subset(dat,Site==1Plot==1) b-reshape(a[-c(1,2)],v.names='Count',idvar='Time',timevar='Sp',direction='wide') d-zoo(b[,-1],b[,1]) d out-rollapply(d, w, cv, na.pad=T, align='right') out # I would thereby have to loop through all my sites and plots which, although it deals with all species at once, still seems exceedingly inefficient. # So the question is, how do I use something like aggregate.zoo or tapply or even lapply to apply rollapply on each species' time series. # The closest I've come is the following two approaches: # First let: datx-list(Site=dat$Site,Plot=dat$Plot,Sp=dat$Sp) daty-dat$Count # Method 1. out1-tapply(seq(along=daty),datx,function(i,x=daty){ rollapply(zoo(x[i]), w, cv, na.pad=T, align='right') }) out1 out1[,,1] # Which works in that it gives me the right answers, but in a format from which I can't figure out how to get back into the format I want. # Method 2. fun-function(x){y-zoo(x);coredata(rollapply(y, w, cv,na.pad=T,align='right'))} out2-aggregate(daty,by=datx,fun) out2 # Which superficially works better, but again only in a format I can't figure out how to use because the output seems to be a mix of data.frame and lists. out2[1,4] out2[1,5] is.data.frame(out2) is.list(out2) # The situation is made more problematic by the fact that the time point of first survey can differ between plots (e.g., site1-plot3 may only start at time-point 3). As in... dat2-dat dat2-dat2[-which(dat2$Plot==3 dat2$Time3),] dat2 # I must therefore ensure that I'm keeping track of the true time associated with each value, not just the order of their occurences. This information is (seemingly) lost by both methods. datx-list(Site=dat2$Site,Plot=dat2$Plot,Sp=dat2$Sp) daty-dat2$Count # Method 1. out3-tapply(seq(along=daty),datx,function(i,x=daty){ rollapply(zoo(x[i]), w, cv, na.pad=T, align='right') }) out3 out3[1,3,1] time(out3[1,3,1]) # Method 2 out4-aggregate(daty,by=datx,fun) out4 time(out4[3,4]) # Am I going about this all wrong? Is there a different package to try? Any thoughts and suggestions are much appreciated! # R 2.12.2 GUI 1.36 Leopard build 32-bit (5691); zoo 1.6-4 # Thanks! # -mark Try ave: dat$cv - ave(dat$Count, dat[c(Site, Plot, Sp)], FUN = function(x) rollapply(zoo(x), 2, cv, na.pad = TRUE, align = right)) -- 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] R-project: plot 2 zoo objects (price series) that have some date mis-matches
On Sun, Apr 3, 2011 at 1:57 PM, algotr8der algotr8...@gmail.com wrote: I have 2 zoo objects - 1) Interest rate spread between 10-YR-US-Treasury and 2-YR-US-Treasury (object name = sprd) 2) SP 500 index (object name = spy) str(spy) ‘zoo’ series from 1976-06-01 to 2011-03-31 Data: num [1:8791] 99.8 100.2 100.1 99.2 98.6 ... Index: Class 'Date' num [1:8791] 2343 2344 2345 2346 2349 ... str(sprd) ‘zoo’ series from 1976-06-01 to 2011-03-31 Data: num [1:9088] 0.68 0.71 0.7 0.77 0.79 0.79 0.82 0.86 0.83 0.83 ... Index: Class 'Date' num [1:9088] 2343 2344 2345 2346 2349 ... Since there are NA data points in object 'sprd' I created another object that omits NA. The name of that object is 'sprdtmp'. str(sprdtmp) ‘zoo’ series from 1976-06-01 to 2011-03-31 Data: atomic [1:8704] 0.68 0.71 0.7 0.77 0.79 0.79 0.82 0.86 0.83 0.83 ... - attr(*, na.action)=Class 'omit' int [1:384] 25 70 95 111 118 128 149 190 224 260 ... Index: Class 'Date' num [1:8704] 2343 2344 2345 2346 2349 ... I want to plot both time series on the same plot with time/date on the x-axis and the axis label quarterly (or monthly). One problem is that the sprdtmp and spy objects do not have the same number of data points as there are times when the equity markets are closed while the interest rate markets are open. For the most part the dates overlap. Would this matter if I try to plot both on the same plot? And how would I go about plotting these objects in one plot. The second part of the plot requires that both are on different scales. I guess I could take a log of the sp and plot that along with the rate spread. But it would be nice for future reference how I could plot 2 series in one plot with 2 difference scales. I spent all night yesterday and this morning trying various options but I cant seem to get this to work. I have read through the ?plot, ?plot.zoo and ?axis documentation without any success. I would greatly appreciate if anyone can point me in the right direction. Thank you kindly. Try this: plot(na.approx(cbind(z1, z2)), screen = 1) and in ?plot.zoo see the multivariate plotting example. where z1 and z2 are two zoo series. Omit screen = 1 if you want them in different panels. -- 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] another question on shapefiles and geom_point in ggplot2
Manuel: As far as I know one needs gpclibPermit() in order to fortify see this: Note: polygon geometry computations in maptools depend on the package gpclib, which has a restricted licence. It is disabled by default; to enable gpclib, type gpclibPermit() I am going to guess that ahmadou dicko doesn't show gpclibPermit() on his code because he loaded it with Rprofile or some other way. I tried to run his code without gpclibPermit() and it wouldn't let me fortify, so not sure how he did it. On the same note, your code didn't work when you were trying to join ai_biotica and ai_biotica@data because you are trying to join polygons and points (just a guess). Try to use fortify instead of fortify.SpatialPolygonsDataFrame. I tested with a different shapefile and it works for me. manuel@data$id = rownames(manuel@data) manuel_fort - fortify(manuel,region=id) manuel_merged - join(manuel_fort,manuel@data, by =id) Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx From: Manuel Spínola mspinol...@gmail.com To: Felipe Carrillo mazatlanmex...@yahoo.com Cc: r-h...@stat.math.ethz.ch; ggpl...@googlegroups.com Sent: Sun, April 3, 2011 1:51:02 PM Subject: Re: another question on shapefiles and geom_point in ggplot2 Thank you very much Felipe, Did you see the solution from ahmadou dicko? He doesn´t use gpclibPermit() I have another option but I cannot get the right fill for the id. See attached map. ai_biotica = readOGR(dsn=C:/ProyectosRespacial/ICE/SIG_Biotica_PHED, layer=AI_BIOTICA_010411_CRTM05) str(ai_biotica) # fortify to get the data fortify.ai_biotica - fortify.SpatialPolygonsDataFrame(ai_biotica,region='Area_Influ') names(fortify.ai_biotica) str(fortify.ai_biotica) levels(fortify.ai_biotica$group) # mapa ggplot(fortify.ai_biotica, aes(x = long, y=lat, group = group)) + geom_polygon(colour = black, fill = NA) geo = read.csv(riqueza_out.csv, sep = ,, header = T) names(geo) str(geo) summary(geo) # mapa con riqueza p = ggplot(geo, aes(x, y)) p + geom_point(aes(size = ACE, colour = ACE)) + theme_bw() + scale_size(name = Número de especies, breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20)) + scale_colour_gradientn(name = 'Número de especies', colours = heat.colors(10), breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20))+ xlab(Longitud) + ylab(Latitud) + opts(axis.text.x = theme_text(size = 8, vjust = 1)) + opts(axis.text.y = theme_text(size = 8, hjust = 1)) + geom_path(aes(x=long,y=lat,group=group, fill=id),data=fortify.ai_biotica) Best, Manuel On 03/04/2011 01:41 p.m., Felipe Carrillo wrote: Manuel: I changed your variable names from x to 'long' and y to 'lat' on the riqueza_out.csv file. The code below should do what you want. Also, since the legend title is kind of long, I broke it down into three lines so you can see more plot area. I am cc'ing the other groups so more people use it if needed. library(rgdal) library(ggplot2) library(sp) library(maptools) gpclibPermit() manuel - readOGR(dsn=., layer=AI_BIOTICA_010411_CRTM05) names(manuel);dim(manuel) slotNames(manuel) # look at the slot names # add the 'id' variable to the shapefile and use it to merge both files manuel@data$id = rownames(manuel@data) # convert shapefile to dataframe manuel.df - as.data.frame(manuel) # fortify to plot with ggplot2 manuel_fort - fortify(manuel,region=id) head(manuel_fort) # Merge shapefile and the as.dataframe shapefile manuel_merged - join(manuel_fort,manuel.df, by =id) head(manuel_merged) # Read in the csv file manuel_points - read.csv(riqueza_out.csv) head(manuel_points);dim(manuel_points) # fortify this one too for the points or else an error will ocurr manuel_points - fortify(manuel_points) manuel_points # Plot the shapefile and overlayed the points over it p - ggplot(manuel_merged, aes(long,lat,group=group)) + geom_polygon(aes(data=manuel_merged,fill=Area_Influ)) + geom_path(color=white) + theme_bw() # remove this if you don't want black and white background p + geom_point(data=manuel_points,aes(size=ACE,colour=ACE,group=NULL)) + scale_size(name = Número\nde\nespecies, breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20)) + scale_colour_gradientn(name = 'Número\nde\nespecies', colours = rainbow(6), breaks = c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20))+ xlab(Longitud) + ylab(Latitud) + opts(axis.text.x = theme_text(size = 8, vjust = 1)) + opts(axis.text.y = theme_text(size = 8, hjust = 1)) Felipe D. Carrillo Supervisory Fishery Biologist Department of the Interior US Fish Wildlife Service California, USA http://www.fws.gov/redbluff/rbdd_jsmp.aspx From: Manuel Spínola mspinol...@gmail.com To: Felipe Carrillo mazatlanmex...@yahoo.com Sent: Sat, April 2, 2011 11:22:24 PM Subject: Re: another question on shapefiles and geom_point in ggplot2 No problem, thank you very much Felipe. Best, Manuel On 03/04/2011 12:19 a.m.,
Re: [R] Quick question about duplicating vectors
Sorry for the imprecision. I haven't used r for a long time... CC: r-help@r-project.org From: dwinsem...@comcast.net To: rhel...@hotmail.com Subject: Re: [R] Quick question about duplicating vectors Date: Sun, 3 Apr 2011 18:17:55 -0400 On Apr 3, 2011, at 6:04 PM, Xiaoxi Gao wrote: Hello R users, I have a quick question, if I have a data set a, say a - c(1,2,3,4,5) You might as well learn to use more precise R terminology... that is a `vector`. and I want to get a new data set b, b - c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5) Could anyone tell me how I can obtain this? ?rep -- David Winsemius, MD West Hartford, CT [[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] power of 2 way ANOVA with interaction
I've been searching for an answer to this for a while but no joy. I have a simple 2-way ANOVA with an interaction. I'd like to determine the power of this test for each factor (factor A, factor B, and the A*B interaction). How can I do this in R? I used to do this with proc Glmpower in SAS, but I can find no analogue in R. [[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] replace last 3 characters of string
Hi, I would like to replace the last tree characters of the values of a certain column in a dataframe. This replacement should only take place if the last three characters correspond to the value /:/ and they should be replaced with (blank) I cannot perform a simple gsub because the characters /:/ might also be present somewhere else in the string values and then they should not be replaced. Could someone help me out on this one? Thx Bert [[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] replace last 3 characters of string
On Mon, Apr 04, 2011 at 01:39:34AM +0200, Bert Jacobs wrote: I would like to replace the last tree characters of the values of a certain column in a dataframe. This replacement should only take place if the last three characters correspond to the value /:/ and they should be replaced with (blank) I cannot perform a simple gsub because the characters /:/ might also be present somewhere else in the string values and then they should not be replaced. Keep reading up on regular expressions, this tends to pay off. Here we use the fact that you can achor a regexp to the end of a string: R aString - abc:def::: R gsub(:::$, , aString) [1] abc:def R Hth, Dirk -- Three out of two people have difficulties with fractions. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] replace last 3 characters of string
Thx I could imagine it was so simple.:) Bert -Original Message- From: Dirk Eddelbuettel [mailto:e...@master.debian.org] On Behalf Of Dirk Eddelbuettel Sent: 04 April 2011 01:39 To: Bert Jacobs Cc: r-help@r-project.org Subject: Re: [R] replace last 3 characters of string On Mon, Apr 04, 2011 at 01:39:34AM +0200, Bert Jacobs wrote: I would like to replace the last tree characters of the values of a certain column in a dataframe. This replacement should only take place if the last three characters correspond to the value /:/ and they should be replaced with (blank) I cannot perform a simple gsub because the characters /:/ might also be present somewhere else in the string values and then they should not be replaced. Keep reading up on regular expressions, this tends to pay off. Here we use the fact that you can achor a regexp to the end of a string: R aString - abc:def::: R gsub(:::$, , aString) [1] abc:def R Hth, Dirk -- Three out of two people have difficulties with fractions. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Linear Model with curve fitting parameter?
Steven: You are exactly right sorry I was confused. ### so log(y-intercept)+log(K) is a constant called b0 (is this right?) lm(log(Q)~log(A)+log(R)+log(S)-1) is fitting the model log(Q)=a*log(A)+r*log(R)+s*log(S) (no beta 0) and lm(log(Q)~log(A)+log(R)+log(S)) is fitting the model log(Q)=b0+a*log(A)+r*log(R)+s*log(S) ## These are the models I am trying to fit and if I have reasoned correctly above then I should be able to fit the below models similarly. manning log(Q)=log(b0)+log(K)+log(A)+r*log(R)+s*log(S) dingman log(Q)=log(b0)+log(K)+a*log(A)+r*log(R)+s*(log(S))^2 bjerklie log(Q)=log(b0)+log(K)+a*log(A)+r*log(R)+s*log(S) ### Thank you for all of your help! Stephen On Fri, Apr 1, 2011 at 2:44 PM, Steven McKinney smckin...@bccrc.ca wrote: -Original Message- From: stephen sefick [mailto:ssef...@gmail.com] Sent: April-01-11 5:44 AM To: Steven McKinney Cc: R help Subject: Re: [R] Linear Model with curve fitting parameter? Setting Z=Q-A would be the incorrect dimensions. I could Z=Q/A. I suspect this is confusion about what Q is. I was presuming that the Q in this following formula was log(Q) with Q from the original data. I have taken the log of the data that I have and this is the model formula without the K part lm(Q~offset(A)+R+S, data=x) If the model is Q=K*A*(R^r)*(S^s) then log(Q) = log(K) + log(A) + r*log(R) + s*log(S) Rearranging yields log(Q) - log(A) = log(K) + r*log(R) + s*log(S) so what I labeled 'Z' below is Z = log(Q) - log(A) = log(Q/A) so Z = log(K) + r*log(R) + s*log(S) and a linear model fit of Z ~ log(R) + log(S) will yield parameter estimates for the linear equation E(Z) = B0 + B1*log(R) + B2*log(S) (E(Z) = expected value of Z) so B0 estimate is an estimate of log(K) B1 estimate is an estimate of r B2 estimate is an estimate of s More details and careful notation will eventually lead to a reasonable description and analysis strategy. Best Steve McKinney Is fitting a nls model the same as fitting an ols? These data are hydraulic data from ~47 sites. To access predictive ability I am removing one site fitting a new model and then accessing the fit with a myriad of model assessment criteria. I should get the same answer with ols vs nls? Thank you for all of your help. Stephen On Thu, Mar 31, 2011 at 8:34 PM, Steven McKinney smckin...@bccrc.ca wrote: -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of stephen sefick Sent: March-31-11 3:38 PM To: R help Subject: [R] Linear Model with curve fitting parameter? I have a model Q=K*A*(R^r)*(S^s) A, R, and S are data I have and K is a curve fitting parameter. I have linearized as log(Q)=log(K)+log(A)+r*log(R)+s*log(S) I have taken the log of the data that I have and this is the model formula without the K part lm(Q~offset(A)+R+S, data=x) What is the formula that I should use? Let Z = Q - A for your logged data. Fitting lm(Z ~ R + S, data = x) should yield intercept parameter estimate = estimate for log(K) R coefficient parameter estimate = estimate for r S coefficient parameter estimate = estimate for s Steven McKinney Statistician Molecular Oncology and Breast Cancer Program British Columbia Cancer Research Centre Thanks for all of your help. I can provide a subset of data if necessary. -- Stephen Sefick | Auburn University | | Biological Sciences | | 331 Funchess Hall | | Auburn, Alabama | | 36849 | |___| | sas0...@auburn.edu | | http://www.auburn.edu/~sas0025 | |___| Let's not spend our time and resources thinking about things that are so little or so large that all they really do for us is puff us up and make us feel like gods. We are mammals, and have not exhausted the annoying little problems of being mammals. -K. Mullis A big computer, a complex algorithm and a long time does not equal science. -Robert Gentleman __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Stephen Sefick
[R] LIVE, ONLINE COURSE: Using R Software for Academic Research Analyses
The Information Institute (http://www.information-institute.org), a charitable, non-profit, educational and scientific organization, in conjunction with faculty from Virginia Commonwealth University (VCU), is offering live, interactive, synchronous online classes (one AM and one PM to reach all time zones) on the use of R for a variety of academic-research statistical analyses. The registration cost for the 14-hour, 5 week course is $225 (student); $295 (faculty); and $325 (practitioner). The informational (and registration) site for the AM version (AM by US Eastern Time) is here: https://www.regonline.com/R-asa-apr-sat-AM This course runs on five consecutive Saturdays from April 23 to May 21, from 11AM until 2PM ET. The informational (and registration) site for the PM version (PM by US Eastern Time) is here: https://www.regonline.com/R-inf-inst-apr This course runs on five consecutive Saturdays from April 23 to May 21, in the evening from 6PM until 9PM ET. The R Project for Statistical Computing provides a comprehensive environment for statistical analysis and graphics and is unrivaled in the availability of new, cutting-edge applications. R is a very powerful system for statistical computations and graphics, and runs on Windows, UNIX and Mac computers. It may be described as a combination of a statistics platform and a programming language. It is freely available for download from http://www.r-project.org/ . We designed this course specifically for researchers and statistical workers who want to learn a powerful new software platform for conducting a variety of essential academic research statistical analyses, including: (1) the statistical graphical displays available in R; (2) simple inference; (3) analysis of variance; (4) simple and multiple regression; (5) logistic regression; (6) analyzing longitudinal data; and (7) simultaneous inference and multiple comparisons. This course illustrates how to use R for common, research-oriented statistical analyses using generic data which is provided with or without the (optional) textbook. We will also use R Commander, the freely-available, menu-driven, statistics-oriented visual interface to R to illustrate how to utilize the various statistical functions. This live, interactive, synchronous online course is conducted in four separate, weekly 3-hour sessions. There is also a pre-class session to ensure that everyone has installed R and R Commander without problems. All classes are recorded and provided to participants so, if you miss a class, you will have the recording. Feel free to email ghub...@vcu.edu with any questions or for more information. Geoff Hubona Information Systems Department Virginia Commonwealth University [[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] Graph many points without hiding some
Thanks heaps for all your help, sorry for the late reply, I'm a bit overwhelmed by all the possibilities here! My variables have different lengths which makes automated randomisation difficult Peter, although not impossible. Thanks for the suggestion of 3D graphs Nick, I'm not sure how I can make it work with this particular project but I'll keep it in mind for the future. I'll try to get this working with transparancy in base graphics (thanks Claudia and Greg), but if that doesn't cut it I'll have to learn ggplot2 - which looks like a good idea anyway as I am very impressed with what it can do. Thankyou Dennis in particular for translating my code into ggplot, this will be a great help as I get started. Samuel On 1 April 2011 05:07, Greg Snow greg.s...@imail.org wrote: Just a note, Base graphics does support transparency as long as the device plotting to supports it. -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- project.org] On Behalf Of Dennis Murphy Sent: Thursday, March 31, 2011 1:36 AM To: Samuel Dennis Cc: R-help@r-project.org Subject: Re: [R] Graph many points without hiding some Hi: I can think of a couple: (1) size reduction of the points; (2) alpha transparency; (3) (1) + (2) From your original plot in base graphics, I reduced cex to 0.2 and it didn't look too bad: plot(rnorm(x,mean=19),rnorm(x),col=3,xlim=c(16,24), cex = 0.2) points(rnorm(x,mean=20),rnorm(x),col=1, cex = 0.2) points(rnorm(x,mean=21),rnorm(x),col=2, cex = 0.2) AFAIK, base graphics doesn't have alpha transparency available, but the ggplot2 package does. One approach is to adjust the alpha transparency on default size points; another is to combine reduced point size with alpha transparency. Here is your example rehashed for ggplot2. require(ggplot2) d - data.frame(x1 = rnorm(1, mean = 19), x2 = rnorm(1, mean = 20), x3 = rnorm(1, mean = 21), x = rnorm(1)) # Basically stacking x1 - x3, creating two new vars named variable and value dm - melt(d, id = 'x') # from reshape package, loads with ggplot2 # Alpha transparency is set to a low level with default point size, # but the colors in the legend are muted by the level of transparency ggplot(dm, aes(x = x, y = value, colour = variable)) + theme_bw() + geom_point(alpha = 0.05) + scale_colour_manual(values = c('x1' = 'black', 'x2' = 'red', 'x3' = 'green')) # A tradeoff is to reduce the point size and increase alpha a bit, but these changes will # also be reflected in the legend. ggplot(dm, aes(x = x, y = value, colour = variable)) + theme_bw() + geom_point(alpha = 0.15, size = 1) + scale_colour_manual(values = c('x1' = 'black', 'x2' = 'red', 'x3' = 'green')) You may well find the legend to be useless for this example, so to get rid of it, ggplot(dm, aes(x = x, y = value, colour = variable)) + theme_bw() + geom_point(alpha = 0.15, size = 1) + scale_colour_manual(values = c('x1' = 'black', 'x2' = 'red', 'x3' = 'green')) + opts(legend.position = 'none') The nice thing about the ggplot2 graph is that you can adjust the point size and alpha transparency to your tastes. The default point size is 2 and the default alpha = 1 (no transparency). HTH, Dennis On Wed, Mar 30, 2011 at 10:04 PM, Samuel Dennis sjdenn...@gmail.com wrote: I have a very large dataset with three variables that I need to graph using a scatterplot. However I find that the first variable gets masked by the other two, so the graph looks entirely different depending on the order of variables. Does anyone have any suggestions how to manage this? This code is an illustration of what I am dealing with: x - 1 plot(rnorm(x,mean=20),rnorm(x),col=1,xlim=c(16,24)) points(rnorm(x,mean=21),rnorm(x),col=2) points(rnorm(x,mean=19),rnorm(x),col=3) gives an entirely different looking graph to: x - 1 plot(rnorm(x,mean=19),rnorm(x),col=3,xlim=c(16,24)) points(rnorm(x,mean=20),rnorm(x),col=1) points(rnorm(x,mean=21),rnorm(x),col=2) despite being identical in all respects except for the order in which the variables are plotted. I have tried using pch=., however the colours are very difficult to discern. I have experimented with a number of other symbols with no real solution. The only way that appears to work is to iterate the plot with a for loop, and progressively add a few numbers from each variable, as below. However although I can do this simply with random numbers as I have done here, this is an extremely cumbersome method to use with real datasets.
Re: [R] How do I modify uniroot function to return .0001 if error ?
eric ericstrom at aol.com writes: I am calling the uniroot function from inside another function using these lines (last two lines of the function) : d - uniroot(k, c(.001, 250), tol=.05) return(d$root) The problem is that on occasion there's a problem with the values I'm passing to uniroot. In those instances uniroot stops and sends a message that it can't calculate the root because f.upper * f.lower is greater than zero. All I'd like to do in those cases is be able to set the return value of my calling function return(d$root) to .0001. But I'm not sure how to pull that off. I tried a few modifications to uniroot but so far no luck. Do not modify uniroot(). Use 'try' or 'tryCatch', for example e - try( d - uniroot(k, c(.001, 250), tol=.05), silent = TRUE ) if (class(e) == try-error) { return(0.0001) } else { return(d$root) } --Hans Werner __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] RGtk2: How to populate an GtkListStore data model?
hello all I am trying to learn how to use the RGtk2 package... so, my first problem is: I don't get the right way for populate my gtkListStore object! any help is welcome... because I am trying several day to mount the code... Thanks in advanced Cleber N. Borges --- # my testing code library(RGtk2) win - gtkWindowNew() datamodel - gtkListStoreNew('gchararray') treeview - gtkTreeViewNew() renderer - gtkCellRendererText() col_0 - gtkTreeViewColumnNewWithAttributes(title=TitleXXX, cell=renderer, text=Bar) nc_next - gtkTreeViewInsertColumn(object=treeview, column=col_0, position=0) gtkTreeViewSetModel( treeview, datamodel ) win$add( treeview ) # is there an alternative function for this? # iter - gtkTreeModelGetIterFirst( datamodel )[[2]] # this function don't give VALID iter # gtkListStoreIterIsValid( datamodel, iter ) result in FALSE iter - gtkListStoreInsert( datamodel, position=0 )[[2]] gtkListStoreIterIsValid( datamodel, iter ) # the help of this function say to terminated in -1 value # but -1 crash the R-pckage (or the gtk)... gtkListStoreSet(object=datamodel, iter=iter, 0, textFoo) # don't make any difference in the window... :-( R version 2.13.0 alpha (2011-03-27 r55091) Platform: i386-pc-mingw32/i386 (32-bit) locale: [1] LC_COLLATE=Portuguese_Brazil.1252 [2] LC_CTYPE=Portuguese_Brazil.1252 [3] LC_MONETARY=Portuguese_Brazil.1252 [4] LC_NUMERIC=C [5] LC_TIME=Portuguese_Brazil.1252 attached base packages: [1] stats graphics grDevices utils datasets methods [7] base other attached packages: [1] RGtk2_2.20.8 loaded via a namespace (and not attached): [1] tools_2.13.0 my gtk version == 2.16.2 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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 R object files on an RWeb server
Hello list: I have some R code/data sets that i'd like to make available to others through RWeb server use. The code/data were saved via save.image and should be made available by the following command sequence: con - url(https:// ) load(con) close(con) However, this is where the problem starts... For one thing, i am using google docs to host the R object file and google docs has secure https URL's, which apparently cannot be handled by R's url(). So my questions are these: 1) where can i host the R object file (for free) with a permanent hot-link and without https, that could be read in by the load() command? 2) is there a correct procedure for all this that can be used when working with Rweb? Packages are not the way to go as it seems one cannot install a package on an RWeb server! Thanks in advance, jose romero [[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] detect filetype (as in unix 'file')
Is there a way in R (in Linux) to detect the type of a file without invoking a shell? E.g to do this: system(file density.plot) density.plot: PDF document, version 1.4 but without using system()? I tried file() and file.info(), but both do display the information I am looking for. -- View this message in context: http://r.789695.n4.nabble.com/detect-filetype-as-in-unix-file-tp3424562p3424562.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] Ordering every row of a matrix while ignoring off diagonal elements
Sorry if this is a stupid question but I've been stuck on how to code so that I can order rows of a matrix without paying attention to the diagonal elements. Say, for example, I have a matrix d: d [,1] [,2] [,3] [,4] [1,] 0.00 2.384158 2.0065682 2.2998856 [2,] 2.384158 0.00 1.4599928 2.4333213 [3,] 2.006568 1.459993 0.000 0.9733285 [4,] 2.299886 0.00 0.9733285 0.000 Then I'd like ordered d to be like: [,1] [,2] [,3] [1,]342 [2,]314 [3,]421 [4,]231 So subject 1's smallest value is in column 3. Subject 2's second smallest value would be 3, etc. Note that subject 4 has two zeros (a tie) but if the diagonals are not in the equation, then the minimum value for this subject is from column 2. Right now I coded off diagonals as missing and then order it that way but I feel like it's cheating. Suggestions?? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] another question on shapefiles and geom_point in ggplot2
Hi all, 2011/4/4 Felipe Carrillo mazatlanmex...@yahoo.com: Manuel: As far as I know one needs gpclibPermit() in order to fortify see this: Note: polygon geometry computations in maptools depend on the package gpclib, which has a restricted licence. It is disabled by default; to enable gpclib, type gpclibPermit() I am going to guess that ahmadou dicko doesn't show gpclibPermit() on his code because he loaded it with Rprofile or some other way. I tried to run his code without gpclibPermit() and it wouldn't let me fortify, so not sure how he did it. On that specific point, Colin Arundel and Roger Bivant released the rgeos package on CRAN a few days [1]. This is a great achievement as it brings bindings to the GEOS C++ lib [2] - long story short, it makes the job the non-free [3] gpclib used to do. In its later release, maptools has an option to check if rgeos if present - if it is the case it is used instead of gpclib: library(maptools) Loading required package: foreign Loading required package: sp Loading required package: lattice Note: polygon geometry computations in maptools depend on the package gpclib, which has a restricted licence. It is disabled by default; to enable gpclib, type gpclibPermit() Checking rgeos availability as gpclib substitute: TRUE ?gpclibPermit Pierre [1] http://cran.r-project.org/web/packages/rgeos/ [2] http://trac.osgeo.org/geos/ [3] https://stat.ethz.ch/pipermail/r-sig-geo/2010-January/007400.html -- Scientist Landcare Research, New Zealand __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] Ordering every row of a matrix while ignoring off diagonal elements
I assume that this is what you did, and I would not call that cheating; it is just a reasonable way to solve the problem: x - as.matrix(read.table(textConnection( 0.00 2.384158 2.0065682 2.2998856 + 2.384158 0.00 1.4599928 2.4333213 + 2.006568 1.459993 0.000 0.9733285 + 2.299886 0.00 0.9733285 0.000))) closeAllConnections() # put NAs in diagonals diag(x) - NA # get the order x.ord - t(apply(x, 1, order)) # remove last column since this is where NAs sort x.ord[, -ncol(x.ord)] [,1] [,2] [,3] [1,]342 [2,]314 [3,]421 [4,]231 On Sun, Apr 3, 2011 at 11:51 PM, Vivian Shih v...@ucla.edu wrote: Sorry if this is a stupid question but I've been stuck on how to code so that I can order rows of a matrix without paying attention to the diagonal elements. Say, for example, I have a matrix d: d [,1] [,2] [,3] [,4] [1,] 0.00 2.384158 2.0065682 2.2998856 [2,] 2.384158 0.00 1.4599928 2.4333213 [3,] 2.006568 1.459993 0.000 0.9733285 [4,] 2.299886 0.00 0.9733285 0.000 Then I'd like ordered d to be like: [,1] [,2] [,3] [1,] 3 4 2 [2,] 3 1 4 [3,] 4 2 1 [4,] 2 3 1 So subject 1's smallest value is in column 3. Subject 2's second smallest value would be 3, etc. Note that subject 4 has two zeros (a tie) but if the diagonals are not in the equation, then the minimum value for this subject is from column 2. Right now I coded off diagonals as missing and then order it that way but I feel like it's cheating. Suggestions?? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Data Munger Guru What is the problem that you are trying to solve? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] replace last 3 characters of string
Will this do it for you: x - c('asdfasdf', 'sadfasdf/:/', 'sadf', 'asdf/:/') sub(/:/$, '', x) [1] asdfasdf sadfasdf sadf asdf On Sun, Apr 3, 2011 at 7:39 PM, Bert Jacobs bert.jac...@figurestofacts.be wrote: Hi, I would like to replace the last tree characters of the values of a certain column in a dataframe. This replacement should only take place if the last three characters correspond to the value /:/ and they should be replaced with (blank) I cannot perform a simple gsub because the characters /:/ might also be present somewhere else in the string values and then they should not be replaced. Could someone help me out on this one? Thx Bert [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Data Munger Guru What is the problem that you are trying to solve? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Structural equation modeling in R(lavaan,sem)
On 3 April 2011 12:38, jouba antr...@hotmail.com wrote: Daer all, I have a question concerning longitudinal data: When we have a longitudinal data and we have to do sem analysis there is in the package lavaan some functions,options in this package that help to do this or we can treat these data like non longitudinal data No, and (qualified) no. 1. There are (AFAIK) no options, functions that are specific to longitudinal data. 2. You don't treat these data as non-longitudinal data, you add parameters that are appropriate though. For example, look at the model shown on http://lavaan.ugent.be. dem60 and dem65 are two measures of the same construct at different timepoints, so there are correlations over time for each pair of measured variables that are measures of that construct - i.e. y1 ~~ y5 3. You would get much better answers on the SEM mailing list - semnet. You can join it here: http://www2.gsu.edu/~mkteer/semnet.html#Joining. Jeremy -- Jeremy Miles Psychology Research Methods Wiki: www.researchmethodsinpsychology.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.