Re: [R] Plotting MDS (multidimensional scaling)

2011-04-03 Thread Jari Oksanen
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

2011-04-03 Thread Alexander Engelhardt

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?

2011-04-03 Thread Krishna Kirti Das
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?

2011-04-03 Thread peter dalgaard

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

2011-04-03 Thread Gabor Grothendieck
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

2011-04-03 Thread Uwe Ligges



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

2011-04-03 Thread Muzna Alvi
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

2011-04-03 Thread mathijsdevaan
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

2011-04-03 Thread Prof Brian Ripley

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

2011-04-03 Thread Alexander Engelhardt

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

2011-04-03 Thread David Winsemius


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

2011-04-03 Thread David Winsemius


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

2011-04-03 Thread Mike Marchywka



 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

2011-04-03 Thread shahab
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

2011-04-03 Thread Duncan Murdoch

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

2011-04-03 Thread Alexander Engelhardt

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

2011-04-03 Thread John Kane
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?

2011-04-03 Thread lcn
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

2011-04-03 Thread Eduardo M. A. M.Mendes
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

2011-04-03 Thread lcn
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?

2011-04-03 Thread Duncan Murdoch

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

2011-04-03 Thread Ted Harding
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?

2011-04-03 Thread John Fox
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?

2011-04-03 Thread Krishna Kirti Das
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?

2011-04-03 Thread Krishna Kirti Das
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

2011-04-03 Thread mustafabinar
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?

2011-04-03 Thread Spencer Graves

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?

2011-04-03 Thread John Fox
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?

2011-04-03 Thread John Fox
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

2011-04-03 Thread John Kane
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

2011-04-03 Thread Samuel Le
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

2011-04-03 Thread David Winsemius


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

2011-04-03 Thread Axel Urbiz
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

2011-04-03 Thread Douglas Bates
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

2011-04-03 Thread Tyler Rinker

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

2011-04-03 Thread Dennis Murphy
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

2011-04-03 Thread David Winsemius


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

2011-04-03 Thread David Winsemius


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

2011-04-03 Thread Greg Snow
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

2011-04-03 Thread Muzna Alvi
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

2011-04-03 Thread andrija djurovic
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

2011-04-03 Thread Greg Snow
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

2011-04-03 Thread Sebastián Daza

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

2011-04-03 Thread Mark Novak

# 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

2011-04-03 Thread Felipe Carrillo
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

2011-04-03 Thread wulei wong
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

2011-04-03 Thread wulei wong
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

2011-04-03 Thread Sebastián Daza

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

2011-04-03 Thread Sebastián Daza

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

2011-04-03 Thread Sebastián Daza

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

2011-04-03 Thread dsg
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

2011-04-03 Thread Steve Friedman
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

2011-04-03 Thread Kin Ming Wong

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)

2011-04-03 Thread jouba

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

2011-04-03 Thread algotr8der
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

2011-04-03 Thread Tyler Rinker

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 ?

2011-04-03 Thread eric
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)

2011-04-03 Thread Daniel Malter
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

2011-04-03 Thread Thomas Lumley
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

2011-04-03 Thread David Winsemius


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

2011-04-03 Thread Xiaoxi Gao

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

2011-04-03 Thread Manuel Spínola

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

2011-04-03 Thread David Winsemius


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

2011-04-03 Thread Daniel Malter
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

2011-04-03 Thread David Winsemius


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

2011-04-03 Thread Gabor Grothendieck
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

2011-04-03 Thread Gabor Grothendieck
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

2011-04-03 Thread Felipe Carrillo
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

2011-04-03 Thread Xiaoxi Gao

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

2011-04-03 Thread Timothy Spier

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

2011-04-03 Thread Bert Jacobs
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

2011-04-03 Thread Dirk Eddelbuettel
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

2011-04-03 Thread Bert Jacobs
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?

2011-04-03 Thread stephen sefick
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

2011-04-03 Thread Geoffrey Hubona
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

2011-04-03 Thread Samuel Dennis
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 ?

2011-04-03 Thread Hans W Borchers
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?

2011-04-03 Thread Cleber N. Borges

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

2011-04-03 Thread jose romero
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')

2011-04-03 Thread Jeroen Ooms
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

2011-04-03 Thread Vivian Shih
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

2011-04-03 Thread Pierre Roudier
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

2011-04-03 Thread jim holtman
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

2011-04-03 Thread jim holtman
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)

2011-04-03 Thread Jeremy Miles
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.