Re: [R] plotmath vector problem; full program enclosed

2010-07-07 Thread Paul Johnson
On Tue, Jul 6, 2010 at 12:41 PM, Duncan Murdoch
murdoch.dun...@gmail.com wrote:
 On 06/07/2010 10:54 AM, Paul Johnson wrote:

 Here's another example of my plotmath whipping boy, the Normal
 distribution.

 You want as.expression(b1), not expression(b1).  The latter means the
 expression consisting of the symbol b1.  The former means take the object
 stored in b1, and convert it to an expression..

 It's not perfect, because you'll end up with mu - -1.96sigma (i.e. two
 minus signs), but it's closer than what you had.

 Duncan Murdoch



Hi, Duncan and David

Thanks for looking.  I suspect from the comment you did not run the
code.  The expression examples I give do work fine already.  But I
have to explicitly put in values like 1.96 to make them work.  I'm
trying to avid that with substitute, which does work for b2, b3, b4,
b5, all but b1.  Why just one?

I'm uploading a picture of it so you can see for yourself:

http://pj.freefaculty.org/R/plotmathwrong.pdf

please look in the middle axis.

Why does only b1 not work, but the rest do?

-- 
Paul E. Johnson
Professor, Political Science
1541 Lilac Lane, Room 504
University of Kansas

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


Re: [R] Gray level mosaic plot with shading_Friendly

2010-07-07 Thread Achim Zeileis

On Tue, 6 Jul 2010, Michael Friendly wrote:


Michael Kubovy wrote:

Suppose we start with

data(Titanic)
mosaic(Titanic, shade = TRUE)

How do I combine the dashed box contours of shading_Friendly to indicate 
negative residuals, with three levels of gray: dark for abs(Pearson Resid) 
 4, lighter for 4  abs(Pearson Resid)  2, and lightest for bs(Pearson 
Resid)  2 ?




Do you mean [1] you want to plot positive residuals in color and negative in 
gray scale?

Or [2] to fold + and - residuals by shading all according to abs(resid), and
distinguishing + from - by the dashed box outlines?

In fact, I designed this coding scheme so that mosaic plots in color (with my 
blue - white - red scheme) would approximately do exactly what
you might want under [2], when rendered in B/W, since the fully saturated red 
and blue are close in  darkness in B/W.


And shading_hcl() has been written to do exactly what you want under [2]. 
While it is hard to come up with colors of different hues in HSV or HLS 
space that have the same brightness (aka lightness/luminance) and the same

colorfulness (aka chroma), this is easy in HCL.


Try
mosaic(Titanic, gp=shading_Friendly)
save as a jpg/png and try converting to B/W with an image program and see if 
this is good enough.


mosaic(Titanic, shade = TRUE)

is the same as

mosaic(Titanic, gp = shading_hcl)

which you can then modify to have different line types

mosaic(Titanic, gp = shading_hcl, gp_args = list(lty = 1:2))

If you print that on a grayscale printer you will see the same plot 
without any chroma, i.e.,


mosaic(Titanic, gp = shading_hcl, gp_args = list(lty = 1:2, c = 0))

The shading_hcl() function is introduced in Zeileis et al. (2007, JCGS), 
see ?shading_hcl, which provides more detailed references to HCL colors 
etc.


Best,
Z


Alternatively, write your own, shading_Kubovy, modeled on

shading_Friendly -
function (observed = NULL, residuals = NULL, expected = NULL,
   df = NULL, h = c(2/3, 0), lty = 1:2, interpolate = c(2, 4),
   eps = 0.01, line_col = black, ...)
{
   shading_hsv(observed = NULL, residuals = NULL, expected = NULL,
   df = NULL, h = h, v = 1, lty = lty, interpolate = interpolate,
   eps = eps, line_col = line_col, p.value = NA, ...)
}
environment: namespace:vcd
attr(,class)
[1] grapcon_generator

In the defaults, lty=1:2 is what distinguishes + and - for outline line type

hope this helps,
-Michael

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



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


[R] R C# (Mono)

2010-07-07 Thread bernde
Hello did anyone ever use C# in connection with R ?
i am looking into R extension but would like to use C# instead of C or C++
i wonder whether anyone has experience in particular with Mono for doing so 

many thanks in advance bernd

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


Re: [R] plotmath vector problem; full program enclosed

2010-07-07 Thread Allan Engelhardt


On 07/07/10 06:03, Paul Johnson wrote:
 [...]
 Hi, Duncan and David

 Thanks for looking.  I suspect from the comment you did not run the
 code.  The expression examples I give do work fine already.  But I
 have to explicitly put in values like 1.96 to make them work.  I'm
 trying to avid that with substitute, which does work for b2, b3, b4,
 b5, all but b1.  Why just one?

 I'm uploading a picture of it so you can see for yourself:

 http://pj.freefaculty.org/R/plotmathwrong.pdf

 please look in the middle axis.

 Why does only b1 not work, but the rest do?



Because you only had one (as.)expression in your original code, 
perhaps?  This version works for me:

### Filename: plotMathProblem.R
### Paul Johnson July 5, 2010
### email mepaulj...@ku.edu
### 2010-07-06 AE : Changes.

   sigma- 10.0
   mu- 4.0

   myx- seq( mu - 3.5*sigma,  mu+ 3.5*sigma, length.out=500)

   myDensity- dnorm(myx,mean=mu,sd=sigma)

   ### xpd needed to allow writing outside strict box of graph
   ### Need big bottom margin to add several x axes
   par(xpd=TRUE, ps=10, mar=c(18,2,2,2))

   plot(myx, myDensity, type=l, xlab=, ylab=Probability Density ,
main=myTitle1, axes=FALSE)
   axis(2, pos= mu - 3.6*sigma)
   axis(1, pos=0)

   lines(c(myx[1],myx[length(myx)]),c(0,0)) ### closes off axes


   addInteriorLine- function(x, m, sd){
 for (i in 1:(length(x))){
   lines( c(x[i],x[i]), c(0, dnorm(x[i],m=m,sd=sd)), lty= 14, lwd=.2)
 }
   }


   dividers- c(qnorm(0.025), -1, 0, 1, qnorm(0.975))
   addInteriorLine(mu+sigma*dividers, mu,sigma)

   # bquote creates an expression that text plotters can use
   t1-  bquote( mu== .(mu))
   mtext(bquote( mu == .(mu)), 1, at=mu, line=-1)


   addInteriorLabel- function(pos1, pos2, m, s){
 area- abs(100*( pnorm(m+pos1*s,m,s)-pnorm(m+pos2*s, m,s)))
 mid- m+0.5*(pos1+pos2)*s
 text(mid, 0.5*dnorm(mid,m,s),label=paste(round(area,2),%))
   }


   addInteriorLabel(dividers[1],dividers[2],  mu, sigma)
   addInteriorLabel(dividers[2],dividers[3],  mu, sigma)
   addInteriorLabel(dividers[3],dividers[4],  mu, sigma)
   addInteriorLabel(dividers[4],dividers[5],  mu, sigma)




### Following is problem point: axis will
### end up with correct labels, except for first point,
### where we end up with b1 instead of mu - 1.96*sigma.
b1- substitute( mu - d*sigma, list(d=*-round(dividers[1],2))*  )
b2- substitute( mu - sigma )
b3- substitute( mu )
b4- substitute( mu + sigma )
b5- substitute( mu + d*sigma, list(d=round(dividers[5],2)) )
##   plot(-20:50,-20:50,type=n,axes=F)
axis(1, line=4,at=mu+dividers*sigma,
labels=*as.expression(c(b1,b2,b3,b4,b5))*, padj=-1)




### This gets right result but have to hard code the dividers
   b1- expression( mu - 1.96*sigma )
   b2- expression( mu - sigma )
   b3- expression( mu )
   b4- expression( mu + sigma )
   b5- expression( mu + 1.96*sigma )
   axis(1, line=8,at=mu+dividers*sigma, labels=c(b1,b2,b3,b4,b5), padj=-1)




Allan

[[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] grayscale wireframe??

2010-07-07 Thread Allan Engelhardt
A standalone example is always appreciated (cf. the posting guide) but 
try and see if help(gray.colors, package=grDevices) is the sort of 
thing you are looking for.


Hope this helps

Allan

On 06/07/10 23:30, Marlin Keith Cox wrote:

I need grayscale formatting for a wireframe.
The only col.regions that I can find are color palettes are all colored:

rainbow(n, s = 1, v = 1, start = 0, end = max(1,n - 1)/n,
 gamma = 1, alpha = 1)
heat.colors(n, alpha = 1)
terrain.colors(n, alpha = 1)
topo.colors(n, alpha = 1)
cm.colors(n, alpha = 1)

The code follows:

X11()
library(lattice)
par(family=serif, cex=1.2)
wireframe(z ~ y*x, mat.df,
   drape = TRUE,col.regions=rainbow(100),
   zlab = list(Water mass error (%),rot=90), zlim=c(-50,180),
   xlab = list(Resistance error (%),rot=-9),
   ylab = list(Length error (%),rot=38),
   scales = list(arrows = FALSE),
screen = list(z = -35, x = -77, y = 10))





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


Re: [R] plotmath vector problem; full program enclosed

2010-07-07 Thread Allan Engelhardt

Ooops, I didn't convert this one to text right for the list.

b1- substitute( mu - d*sigma, list(d=*-round(dividers[1],2))*  )


should be

b1- substitute( mu - d*sigma, list(d=-round(dividers[1],2))  )


and similarly for

labels=*as.expression(c(b1,b2,b3,b4,b5))*, padj=-1)

read

labels=as.expression(c(b1,b2,b3,b4,b5)), padj=-1)

Apologies


Allan

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


Re: [R] ctree ordering nodes

2010-07-07 Thread Achim Zeileis

On Tue, 6 Jul 2010, Paras Sharma wrote:


Hello,

When using the ctree function, from library (party) what is the syntax to
order the Variables in the nodes in a specific way?

For example, how would I specify to make a binary come first, then a
continuous variable?


Not sure what you mean here. There is only one splitting variable in each 
node... The order in which the variables are tested is simply the order in 
which they appear in the model formula.


If you want to force ctree() to use a specific variable in the first split 
and a specific variable in the second split, this is currently not 
possible. ctree() always selects the locally optimal splitting variable.



Also is there a way to force ctree to show variables which are not
significant?


You can inspect the fitted ctree in various ways. See
  vignette(party, package = party)
for a few helpful examples, especially Section 5.3.

If you want to construct your own trees using different algorithms, the 
package partykit might be of interest which is currently under 
development on R-Forge.


hth,
Z

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


Re: [R] Visualization of coefficients

2010-07-07 Thread Tal Galili
Hello David,
Thanks to your posting I started looking at the function in the arm package.
 It appears this function is quite mature, and offers (for example) the
ability to easily overlap coefficients from several models.

I updated the post I published on the subject, so at the end of it I give an
example of comparing the coef of several models:
http://www.r-statistics.com/2010/07/visualization-of-regression-coefficients-in-r/

Thanks again for the pointer.

Best,
Tal




Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
--




On Wed, Jul 7, 2010 at 12:02 AM, David Atkins datk...@u.washington.eduwrote:



 FYI, there is already a function coefplot in the arm package; for example,
 compare:

  library(arm)
 Loading required package: MASS
 Loading required package: Matrix
 [snip]
 Attaching package: 'arm'

 The following object(s) are masked from 'package:coda':

traceplot


  data(Mroz, package = car)
  fm - glm(lfp ~ ., data = Mroz, family = binomial)
  coefplot(fm)

 with version below.

 cheeres, Dave

 
  detach(package:arm)

  coefplot - function(object, df = NULL, level = 0.95, parm = NULL,
 +labels = TRUE, xlab = Coefficient confidence intervals, ylab = ,
 +xlim = NULL, ylim = NULL,
 +las = 1, lwd = 1, lty = c(1, 2), pch = 19, col = 1,
 +length = 0, angle = 30, code = 3, ...)
 + {
 +cf - coef(object)
 +se - sqrt(diag(vcov(object)))
 +if(is.null(parm)) parm - seq_along(cf)
 +if(is.numeric(parm) | is.logical(parm)) parm - names(cf)[parm]
 +if(is.character(parm)) parm - which(names(cf) %in% parm)
 +cf - cf[parm]
 +se - se[parm]
 +k - length(cf)
 +
 +if(is.null(df)) {
 +  df - if(identical(class(object), lm)) df.residual(object) else 0
 +}
 +
 +critval - if(df  0  is.finite(df)) {
 +  qt((1 - level)/2, df = df)
 +} else {
 +  qnorm((1 - level)/2)
 +}
 +ci1 - cf + critval * se
 +ci2 - cf - critval * se
 +
 +lwd - rep(lwd, length.out = 2)
 +lty - rep(lty, length.out = 2)
 +pch - rep(pch, length.out = k)
 +col - rep(col, length.out = k)
 +
 +if(is.null(xlim)) xlim - range(c(0, min(ci1), max(ci2)))
 +if(is.null(ylim)) ylim - c(1 - 0.05 * k, 1.05 * k)
 +
 +if(isTRUE(labels)) labels - names(cf)
 +if(identical(labels, FALSE)) labels - 
 +labels - rep(labels, length.out = k)
 +
 +plot(0, 0, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab,
 +  axes = FALSE, type = n, las = las, ...)
 +arrows(ci1, 1:k, ci2, 1:k, lty = lty[1], lwd = lwd[1], col = col,
 +  length = length, angle = angle, code = code)
 +points(cf, 1:k, pch = pch, col = col)
 +abline(v = 0, lty = lty[2], lwd = lwd[2])
 +axis(1)
 +axis(2, at = 1:k, labels = labels, las = las)
 +box()
 + }
 
 
  coefplot(fm, parm = -1)




 Achim Zeileis wrote:

 I've thought about adding a plot() method for the coeftest() function in
 the lmtest package. Essentially, it relies on a coef() and a vcov()
 method being available - and that a central limit theorem holds. For
 releasing it as a general function in the package the code is still too
 raw, but maybe it's useful for someone on the list. Hence, I've included
 it below.

 An example would be to visualize all coefficients except the intercept for
 the Mroz data:

 data(Mroz, package = car)
 fm - glm(lfp ~ ., data = Mroz, family = binomial)
 coefplot(fm, parm = -1)

 hth,
 Z

 coefplot - function(object, df = NULL, level = 0.95, parm = NULL,
   labels = TRUE, xlab = Coefficient confidence intervals, ylab = ,
   xlim = NULL, ylim = NULL,
   las = 1, lwd = 1, lty = c(1, 2), pch = 19, col = 1,
   length = 0, angle = 30, code = 3, ...)
 {
   cf - coef(object)
   se - sqrt(diag(vcov(object)))
   if(is.null(parm)) parm - seq_along(cf)
   if(is.numeric(parm) | is.logical(parm)) parm - names(cf)[parm]
   if(is.character(parm)) parm - which(names(cf) %in% parm)
   cf - cf[parm]
   se - se[parm]
   k - length(cf)

   if(is.null(df)) {
 df - if(identical(class(object), lm)) df.residual(object) else 0
   }

   critval - if(df  0  is.finite(df)) {
 qt((1 - level)/2, df = df)
   } else {
 qnorm((1 - level)/2)
   }
   ci1 - cf + critval * se
   ci2 - cf - critval * se

   lwd - rep(lwd, length.out = 2)
   lty - rep(lty, length.out = 2)
   pch - rep(pch, length.out = k)
   col - rep(col, length.out = k)

   if(is.null(xlim)) xlim - range(c(0, min(ci1), max(ci2)))
   if(is.null(ylim)) ylim - c(1 - 0.05 * k, 1.05 * k)

   if(isTRUE(labels)) labels - names(cf)
   if(identical(labels, FALSE)) labels - 
   labels - rep(labels, length.out = k)

   plot(0, 0, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab,
 axes = FALSE, type = n, las = las, ...)
   arrows(ci1, 1:k, ci2, 

[R] Wavelet

2010-07-07 Thread nuncio m
Hi useRs,
   Is it possible to get MORLET wavelet in R
Thanks
nuncio

-- 
Nuncio.M
Research Scientist
National Center for Antarctic and Ocean research
Head land Sada
Vasco da Gamma
Goa-403804

[[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] A question about conducting crossed random effects in R

2010-07-07 Thread Paul Chatfield

There is a mixed effects e-mail list you might want to join for more in depth
discussion of these topics - you can subscribe here
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models.

In general the format for crossed effects would be
lmer(y~f1+f2+(1|r1)+(1|r2))   where f1, f2 are fixed effects and r1, r2 are
random.

I believe that the random lines correspond to those in SAS as below
random int/subject=r1; (1|r1) for random varying intercept with each
r1
random f1/subject=r1;  (f1|r1)for random varying slope with each r1

However, in your case, I suspect that the model might not run because
crossed random effects generally take time and with 3 way interactions
that's probably too much for R.  Either way, you can check that by building
your model up slowly and seeing if it runs for the simpler case.  Let me
know how far you get,

Paul

-- 
View this message in context: 
http://r.789695.n4.nabble.com/A-question-about-conducting-crossed-random-effects-in-R-tp2278443p2280578.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] How can I calculate Q-correlations?

2010-07-07 Thread Timo Stolz
I need to calculate Q-correlations, in order to quantify differences
between the profiles of my tested persons. Has anybody any experience
doing that? Which command/package targets Q-correlations?

Timo

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


Re: [R] Visualization of coefficients

2010-07-07 Thread Allan Engelhardt
Thanks Tal.  Nice summary on the web page.  I think the last example 
would be even better if it was a stand-alone piece of code (i.e.: with 
data) that we could run.  For example


library(arm)
data(Mroz, package = car)
M1-  glm(lfp ~ ., data = Mroz, family = binomial)
M2- bayesglm(lfp ~ ., data = Mroz, family = binomial)
M3-  glm(lfp ~ ., data = Mroz, family = binomial(probit))
coefplot(M2, xlim=c(-2, 6),intercept=TRUE)
coefplot(M1, add=TRUE, col.pts=red,  intercept=TRUE)
coefplot(M3, add=TRUE, col.pts=blue, intercept=TRUE, offset=0.2)


Allan

On 07/07/10 09:16, Tal Galili wrote:

Hello David,
Thanks to your posting I started looking at the function in the arm 
package.  It appears this function is quite mature, and offers (for 
example) the ability to easily overlap coefficients from several models.


I updated the post I published on the subject, so at the end of it I 
give an example of comparing the coef of several models:

http://www.r-statistics.com/2010/07/visualization-of-regression-coefficients-in-r/

Thanks again for the pointer.

Best,
Tal
[...]


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


Re: [R] Wavelet

2010-07-07 Thread Allan Engelhardt

Try

RSiteSearch(MORLET)

before you post.

Allan

On 07/07/10 09:38, nuncio m wrote:

Hi useRs,
Is it possible to get MORLET wavelet in R
Thanks
nuncio




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


Re: [R] Visualization of coefficients

2010-07-07 Thread Achim Zeileis

On Wed, 7 Jul 2010, Tal Galili wrote:


Hello David,
Thanks to your posting I started looking at the function in the arm package.
 It appears this function is quite mature, and offers (for example) the
ability to easily overlap coefficients from several models.


Re: more mature. arm's coefplot() is more flexible in certain respects, 
mine is more convenient in others. The overlay functionality is something 
arm's coefplot() is better in and it also as some further options 
(vertical vs. horizontal etc.). My coefplot() has the advantage that it 
does not need any modification as long as coef() and vcov() methods are 
available. Furthermore, level can specify the significance level 
(instead of always using one and two standard errors, respectively).

But it shouldn't be too hard to create a superset of all options.


I updated the post I published on the subject, so at the end of it I give an
example of comparing the coef of several models:
http://www.r-statistics.com/2010/07/visualization-of-regression-coefficient
s-in-r/


As Allan pointed out in his reply, something fully reproducible would be 
nice. Also, if you keep the example with quasi-complete separation, it 
would be worth pointing this out. (Because the maximum likelihood 
estimator is Infinity in this case.)


Finally, the Poisson model in comparison with the binomial models does not 
make much sense, I guess.


Best,
Z


Thanks again for the pointer.

Best,
Tal




Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
---
---




On Wed, Jul 7, 2010 at 12:02 AM, David Atkins datk...@u.washington.edu
wrote:


  FYI, there is already a function coefplot in the arm package;
  for example, compare:

   library(arm)
  Loading required package: MASS
  Loading required package: Matrix
  [snip]
  Attaching package: 'arm'

  The following object(s) are masked from 'package:coda':

     traceplot

   data(Mroz, package = car)
   fm - glm(lfp ~ ., data = Mroz, family = binomial)
 coefplot(fm)

with version below.

cheeres, Dave


 detach(package:arm)

 coefplot - function(object, df = NULL, level = 0.95, parm = NULL,
+    labels = TRUE, xlab = Coefficient confidence intervals, ylab =
,
+    xlim = NULL, ylim = NULL,
+    las = 1, lwd = 1, lty = c(1, 2), pch = 19, col = 1,
+    length = 0, angle = 30, code = 3, ...)
+ {
+    cf - coef(object)
+    se - sqrt(diag(vcov(object)))
+    if(is.null(parm)) parm - seq_along(cf)
+    if(is.numeric(parm) | is.logical(parm)) parm - names(cf)[parm]
+    if(is.character(parm)) parm - which(names(cf) %in% parm)
+    cf - cf[parm]
+    se - se[parm]
+    k - length(cf)
+
+    if(is.null(df)) {
+      df - if(identical(class(object), lm)) df.residual(object)
else 0
+    }
+
+    critval - if(df  0  is.finite(df)) {
+      qt((1 - level)/2, df = df)
+    } else {
+      qnorm((1 - level)/2)
+    }
+    ci1 - cf + critval * se
+    ci2 - cf - critval * se
+
+    lwd - rep(lwd, length.out = 2)
+    lty - rep(lty, length.out = 2)
+    pch - rep(pch, length.out = k)
+    col - rep(col, length.out = k)
+
+    if(is.null(xlim)) xlim - range(c(0, min(ci1), max(ci2)))
+    if(is.null(ylim)) ylim - c(1 - 0.05 * k, 1.05 * k)
+
+    if(isTRUE(labels)) labels - names(cf)
+    if(identical(labels, FALSE)) labels - 
+    labels - rep(labels, length.out = k)
+
+    plot(0, 0, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab,
+      axes = FALSE, type = n, las = las, ...)
+    arrows(ci1, 1:k, ci2, 1:k, lty = lty[1], lwd = lwd[1], col = col,
+      length = length, angle = angle, code = code)
+    points(cf, 1:k, pch = pch, col = col)
+    abline(v = 0, lty = lty[2], lwd = lwd[2])
+    axis(1)
+    axis(2, at = 1:k, labels = labels, las = las)
+    box()
+ }


 coefplot(fm, parm = -1)




Achim Zeileis wrote:

I've thought about adding a plot() method for the coeftest() function
in
the lmtest package. Essentially, it relies on a coef() and a vcov()
method being available - and that a central limit theorem holds. For
releasing it as a general function in the package the code is still
too
raw, but maybe it's useful for someone on the list. Hence, I've
included
it below.

An example would be to visualize all coefficients except the intercept
for
the Mroz data:

data(Mroz, package = car)
fm - glm(lfp ~ ., data = Mroz, family = binomial)
coefplot(fm, parm = -1)

hth,
Z

coefplot - function(object, df = NULL, level = 0.95, parm = NULL,
  labels = TRUE, xlab = Coefficient confidence intervals, ylab = ,
  xlim = NULL, ylim = NULL,
  las = 1, lwd = 1, lty = c(1, 2), pch = 19, col = 1,
  length = 0, angle = 30, code = 3, ...)
{
  cf - coef(object)
  se - sqrt(diag(vcov(object)))
  if(is.null(parm)) parm - seq_along(cf)
  

Re: [R] Wavelet

2010-07-07 Thread stephen sefick
I would even look at the packages page...  There are plenty of them.
Please, at least, do minimal research before clogging up this list.
Also, please read the posting guide that is appended to every mail to
this list.

On Wed, Jul 7, 2010 at 4:03 AM, Allan Engelhardt all...@cybaea.com wrote:
 Try

 RSiteSearch(MORLET)

 before you post.

 Allan

 On 07/07/10 09:38, nuncio m wrote:

 Hi useRs,
            Is it possible to get MORLET wavelet in R
 Thanks
 nuncio



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

| Auburn University   |
| Department of 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

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


Re: [R] Visualization of coefficients

2010-07-07 Thread Tal Galili
Hi Achim and Allan,
I updated the post with Allan's example (thanks Allan).

Achim, you wrote:
Finally, the Poisson model in comparison with the binomial models does not
make much sense, I guess.
I agree.  I wanted something to showcase the function on 3 models (with the
same predictors), and that's the easiest I could think of.  If you'd think
of a smarter example I'd be happy to incorporate it.

Best,
Tal



Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
--




On Wed, Jul 7, 2010 at 12:10 PM, Achim Zeileis achim.zeil...@uibk.ac.atwrote:

 On Wed, 7 Jul 2010, Tal Galili wrote:

  Hello David,
 Thanks to your posting I started looking at the function in the arm
 package.
  It appears this function is quite mature, and offers (for example) the
 ability to easily overlap coefficients from several models.


 Re: more mature. arm's coefplot() is more flexible in certain respects,
 mine is more convenient in others. The overlay functionality is something
 arm's coefplot() is better in and it also as some further options (vertical
 vs. horizontal etc.). My coefplot() has the advantage that it does not need
 any modification as long as coef() and vcov() methods are available.
 Furthermore, level can specify the significance level (instead of always
 using one and two standard errors, respectively).
 But it shouldn't be too hard to create a superset of all options.


  I updated the post I published on the subject, so at the end of it I give
 an
 example of comparing the coef of several models:

 http://www.r-statistics.com/2010/07/visualization-of-regression-coefficient
 s-in-r/


 As Allan pointed out in his reply, something fully reproducible would be
 nice. Also, if you keep the example with quasi-complete separation, it would
 be worth pointing this out. (Because the maximum likelihood estimator is
 Infinity in this case.)

 Finally, the Poisson model in comparison with the binomial models does not
 make much sense, I guess.

 Best,
 Z

  Thanks again for the pointer.

 Best,
 Tal




 Contact
 Details:---
 Contact me: tal.gal...@gmail.com |  972-52-7275845
 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
 www.r-statistics.com (English)

 ---
 ---




 On Wed, Jul 7, 2010 at 12:02 AM, David Atkins datk...@u.washington.edu
 wrote:


  FYI, there is already a function coefplot in the arm package;
  for example, compare:

   library(arm)
  Loading required package: MASS
  Loading required package: Matrix
  [snip]
  Attaching package: 'arm'

  The following object(s) are masked from 'package:coda':

 traceplot

   data(Mroz, package = car)
   fm - glm(lfp ~ ., data = Mroz, family = binomial)
  coefplot(fm)

 with version below.

 cheeres, Dave

 
  detach(package:arm)

  coefplot - function(object, df = NULL, level = 0.95, parm = NULL,
 +labels = TRUE, xlab = Coefficient confidence intervals, ylab =
 ,
 +xlim = NULL, ylim = NULL,
 +las = 1, lwd = 1, lty = c(1, 2), pch = 19, col = 1,
 +length = 0, angle = 30, code = 3, ...)
 + {
 +cf - coef(object)
 +se - sqrt(diag(vcov(object)))
 +if(is.null(parm)) parm - seq_along(cf)
 +if(is.numeric(parm) | is.logical(parm)) parm - names(cf)[parm]
 +if(is.character(parm)) parm - which(names(cf) %in% parm)
 +cf - cf[parm]
 +se - se[parm]
 +k - length(cf)
 +
 +if(is.null(df)) {
 +  df - if(identical(class(object), lm)) df.residual(object)
 else 0
 +}
 +
 +critval - if(df  0  is.finite(df)) {
 +  qt((1 - level)/2, df = df)
 +} else {
 +  qnorm((1 - level)/2)
 +}
 +ci1 - cf + critval * se
 +ci2 - cf - critval * se
 +
 +lwd - rep(lwd, length.out = 2)
 +lty - rep(lty, length.out = 2)
 +pch - rep(pch, length.out = k)
 +col - rep(col, length.out = k)
 +
 +if(is.null(xlim)) xlim - range(c(0, min(ci1), max(ci2)))
 +if(is.null(ylim)) ylim - c(1 - 0.05 * k, 1.05 * k)
 +
 +if(isTRUE(labels)) labels - names(cf)
 +if(identical(labels, FALSE)) labels - 
 +labels - rep(labels, length.out = k)
 +
 +plot(0, 0, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab,
 +  axes = FALSE, type = n, las = las, ...)
 +arrows(ci1, 1:k, ci2, 1:k, lty = lty[1], lwd = lwd[1], col = col,
 +  length = length, angle = angle, code = code)
 +points(cf, 1:k, pch = pch, col = col)
 +abline(v = 0, lty = lty[2], lwd = lwd[2])
 +axis(1)
 +axis(2, at = 1:k, labels = labels, las = las)
 +box()
 + }
 
 
  coefplot(fm, parm = -1)




 Achim Zeileis wrote:

 

Re: [R] Visualization of coefficients

2010-07-07 Thread Achim Zeileis

On Wed, 7 Jul 2010, Tal Galili wrote:


Hi Achim and Allan,I updated the post with Allan's example (thanks Allan).


Thanks!


Achim, you wrote:
Finally, the Poisson model in comparison with the binomial models does not
make much sense, I guess.
I agree.  I wanted something to showcase the function on 3 models (with the
same predictors), and that's the easiest I could think of.  If you'd think
of a smarter example I'd be happy to incorporate it.


You could generate Poisson data and then fit the binomial model to the 
threshold version of the response.


But I guess that would be a bit over the top. Also, one could argue in 
that case that a complementary log-log link should be employed.


Hence, I would simply say (verbally) that it works for baysglm, glm, lm, 
polr objects and that a default method is available which takes 
pre-computed coefficients and associated standard errors from any suitable 
model.


Best,
Z


Best,
Tal



Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
---
---




On Wed, Jul 7, 2010 at 12:10 PM, Achim Zeileis achim.zeil...@uibk.ac.at
wrote:
  On Wed, 7 Jul 2010, Tal Galili wrote:

Hello David,
Thanks to your posting I started looking at the
function in the arm package.
 It appears this function is quite mature, and
offers (for example) the
ability to easily overlap coefficients from several
models.


Re: more mature. arm's coefplot() is more flexible in certain
respects, mine is more convenient in others. The overlay functionality
is something arm's coefplot() is better in and it also as some further
options (vertical vs. horizontal etc.). My coefplot() has the
advantage that it does not need any modification as long as coef() and
vcov() methods are available. Furthermore, level can specify the
significance level (instead of always using one and two standard
errors, respectively).
But it shouldn't be too hard to create a superset of all options.

  I updated the post I published on the subject, so at the
  end of it I give an
  example of comparing the coef of several models:
http://www.r-statistics.com/2010/07/visualization-of-regression-coefficient

  s-in-r/


As Allan pointed out in his reply, something fully reproducible would
be nice. Also, if you keep the example with quasi-complete separation,
it would be worth pointing this out. (Because the maximum likelihood
estimator is Infinity in this case.)

Finally, the Poisson model in comparison with the binomial models does
not make much sense, I guess.

Best,
Z

Thanks again for the pointer.

Best,
Tal




Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il
(Hebrew) |
www.r-statistics.com (English)
---

---




On Wed, Jul 7, 2010 at 12:02 AM, David Atkins
datk...@u.washington.edu
wrote:


     FYI, there is already a function coefplot in the arm
package;
     for example, compare:

      library(arm)
     Loading required package: MASS
     Loading required package: Matrix
     [snip]
     Attaching package: 'arm'

     The following object(s) are masked from 'package:coda':

        traceplot

      data(Mroz, package = car)
      fm - glm(lfp ~ ., data = Mroz, family = binomial)
 coefplot(fm)

with version below.

cheeres, Dave


 detach(package:arm)

 coefplot - function(object, df = NULL, level = 0.95, parm =
NULL,
+    labels = TRUE, xlab = Coefficient confidence intervals,
ylab =
,
+    xlim = NULL, ylim = NULL,
+    las = 1, lwd = 1, lty = c(1, 2), pch = 19, col = 1,
+    length = 0, angle = 30, code = 3, ...)
+ {
+    cf - coef(object)
+    se - sqrt(diag(vcov(object)))
+    if(is.null(parm)) parm - seq_along(cf)
+    if(is.numeric(parm) | is.logical(parm)) parm -
names(cf)[parm]
+    if(is.character(parm)) parm - which(names(cf) %in% parm)
+    cf - cf[parm]
+    se - se[parm]
+    k - length(cf)
+
+    if(is.null(df)) {
+      df - if(identical(class(object), lm))
df.residual(object)
else 0
+    }
+
+    critval - if(df  0  is.finite(df)) {
+      qt((1 - level)/2, df = df)
+    } else {
+      qnorm((1 - level)/2)
+    }
+    ci1 - cf + critval * se
+    ci2 - cf - critval * se
+
+    lwd - rep(lwd, length.out = 2)
+    lty - rep(lty, length.out = 2)
+    pch - rep(pch, length.out = k)
+    col - rep(col, length.out = k)
+
+    if(is.null(xlim)) xlim - range(c(0, min(ci1), max(ci2)))
+    if(is.null(ylim)) ylim - c(1 - 0.05 * k, 1.05 * k)
+
+    if(isTRUE(labels)) labels - names(cf)
+    

Re: [R] Visualization of coefficients

2010-07-07 Thread Tal Galili
I Achim,
I retained the example (so to illustrate the use of the function) - but
pointed out to it's nonsensical nature.
Credit was mentioned to both you and Allan.

Thanks,
Tal

Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
--




On Wed, Jul 7, 2010 at 12:31 PM, Achim Zeileis achim.zeil...@uibk.ac.atwrote:

 On Wed, 7 Jul 2010, Tal Galili wrote:

  Hi Achim and Allan,I updated the post with Allan's example (thanks Allan).


 Thanks!


  Achim, you wrote:
 Finally, the Poisson model in comparison with the binomial models does
 not
 make much sense, I guess.
 I agree.  I wanted something to showcase the function on 3 models (with
 the
 same predictors), and that's the easiest I could think of.  If you'd think
 of a smarter example I'd be happy to incorporate it.


 You could generate Poisson data and then fit the binomial model to the
 threshold version of the response.

 But I guess that would be a bit over the top. Also, one could argue in that
 case that a complementary log-log link should be employed.

 Hence, I would simply say (verbally) that it works for baysglm, glm, lm,
 polr objects and that a default method is available which takes pre-computed
 coefficients and associated standard errors from any suitable model.

 Best,

 Z

  Best,
 Tal



 Contact
 Details:---
 Contact me: tal.gal...@gmail.com |  972-52-7275845
 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
 www.r-statistics.com (English)

 ---
 ---




 On Wed, Jul 7, 2010 at 12:10 PM, Achim Zeileis achim.zeil...@uibk.ac.at
 wrote:
  On Wed, 7 Jul 2010, Tal Galili wrote:

Hello David,
Thanks to your posting I started looking at the
function in the arm package.
 It appears this function is quite mature, and
offers (for example) the
ability to easily overlap coefficients from several
models.


 Re: more mature. arm's coefplot() is more flexible in certain
 respects, mine is more convenient in others. The overlay functionality
 is something arm's coefplot() is better in and it also as some further
 options (vertical vs. horizontal etc.). My coefplot() has the
 advantage that it does not need any modification as long as coef() and
 vcov() methods are available. Furthermore, level can specify the
 significance level (instead of always using one and two standard
 errors, respectively).
 But it shouldn't be too hard to create a superset of all options.

  I updated the post I published on the subject, so at the
  end of it I give an
  example of comparing the coef of several models:

 http://www.r-statistics.com/2010/07/visualization-of-regression-coefficient

  s-in-r/


 As Allan pointed out in his reply, something fully reproducible would
 be nice. Also, if you keep the example with quasi-complete separation,
 it would be worth pointing this out. (Because the maximum likelihood
 estimator is Infinity in this case.)

 Finally, the Poisson model in comparison with the binomial models does
 not make much sense, I guess.

 Best,
 Z

 Thanks again for the pointer.

 Best,
 Tal




 Contact
 Details:---
 Contact me: tal.gal...@gmail.com |  972-52-7275845
 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il
 (Hebrew) |
 www.r-statistics.com (English)

 ---

 ---




 On Wed, Jul 7, 2010 at 12:02 AM, David Atkins
 datk...@u.washington.edu
 wrote:


  FYI, there is already a function coefplot in the arm
 package;
  for example, compare:

   library(arm)
  Loading required package: MASS
  Loading required package: Matrix
  [snip]
  Attaching package: 'arm'

  The following object(s) are masked from 'package:coda':

 traceplot

   data(Mroz, package = car)
   fm - glm(lfp ~ ., data = Mroz, family = binomial)
  coefplot(fm)

 with version below.

 cheeres, Dave

 
  detach(package:arm)

  coefplot - function(object, df = NULL, level = 0.95, parm =
 NULL,
 +labels = TRUE, xlab = Coefficient confidence intervals,
 ylab =
 ,
 +xlim = NULL, ylim = NULL,
 +las = 1, lwd = 1, lty = c(1, 2), pch = 19, col = 1,
 +length = 0, angle = 30, code = 3, ...)
 + {
 +cf - coef(object)
 +se - sqrt(diag(vcov(object)))
 +if(is.null(parm)) parm - seq_along(cf)
 +if(is.numeric(parm) | is.logical(parm)) parm -
 names(cf)[parm]
 +if(is.character(parm)) parm - 

Re: [R] plotmath vector problem; full program enclosed

2010-07-07 Thread Duncan Murdoch

On 07/07/2010 1:03 AM, Paul Johnson wrote:

On Tue, Jul 6, 2010 at 12:41 PM, Duncan Murdoch
murdoch.dun...@gmail.com wrote:
  

On 06/07/2010 10:54 AM, Paul Johnson wrote:


Here's another example of my plotmath whipping boy, the Normal
distribution.

  

You want as.expression(b1), not expression(b1).  The latter means the
expression consisting of the symbol b1.  The former means take the object
stored in b1, and convert it to an expression..

It's not perfect, because you'll end up with mu - -1.96sigma (i.e. two
minus signs), but it's closer than what you had.

Duncan Murdoch





Hi, Duncan and David

Thanks for looking.  I suspect from the comment you did not run the
code.


And I'm certain from your comment that you didn't run my code, or read 
the explanation carefully.


Duncan Murdoch

  The expression examples I give do work fine already.  But I
have to explicitly put in values like 1.96 to make them work.  I'm
trying to avid that with substitute, which does work for b2, b3, b4,
b5, all but b1.  Why just one?
  
I'm uploading a picture of it so you can see for yourself:


http://pj.freefaculty.org/R/plotmathwrong.pdf

please look in the middle axis.

Why does only b1 not work, but the rest do?




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


[R] sqlQuery not functioning properly

2010-07-07 Thread Simon . Gingins
Hello,

I am using Rcommander to generate codes and when I go
through it to load my excel file, it works. But every time I
re-open my script and try to load it directly with the
script Rcommander generated, I get an error message.
Basically the code works one time and then no more.

Here's the script and the error message I get:

data - sqlQuery(channel = 2, select * from
[Questionnaire$])

[5] ERREUR:  
  ']' inattendu(e) dans pal - sqlQuery(channel = 2, select
* from [Questionnaire$]

I would be very grateful if someone could help me with that,
Thanks!

Simon

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


Re: [R] Break in the y-axis

2010-07-07 Thread Jim Lemon

On 07/07/2010 02:39 AM, beloitstudent wrote:


Thanks for the advice!  It has worked for the most part.  However, I am
still coming up with an error message when placing my break line in the axis
that I'm not sure what it means.  If you could help me out, that would be
fantastic...otherwise I might just have to see if I can add it on
powerpoint.  Here is the code you gave me and my script that doesn't work.
...
axis.break(axis=2, breakpos=50,style=slash,pos=-20,brw=0.02)*


This line should be:

axis.break(axis=2, breakpos=50,style=slash,pos=-20,brw=0.02)

That is, you have passed a character string instead of a number.

Jim

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


Re: [R] Help With ANOVA (corrected please ignore last email)

2010-07-07 Thread Dennis Murphy
Hi:

I'd suggest looking at the following plot (data in original post, copied
below):

library(lattice)
stripplot(Intensity ~ Group, data = zzzanova)

Some things stand out in this plot that merit attention.

As Josh Wiley pointed out in an earlier reply, the concentration of -4.60517
values
in this data set should be a concern. Of your 54 responses, 26 of them have
this value, spread across all six groups - in addition, this value is
concentrated
in the first few groups and becomes rarer in the later groups. Consider:

 with(zzzanova, table(Group))
Group
Group1 Group2 Group3 Group4 Group5 Group6
 9  8  9 10  9  9
 with(subset(zzzanova, Intensity  -1), table(Intensity, Group))
  Group
Intensity  Group1 Group2 Group3 Group4 Group5 Group6
  -4.60517  8  5  4  4  4  1

Also note that there is one other value with a negative intensity in group 6
(-0.024).
(Side note: Is negative intensity meaningful in the context of your
scientific problem?)

I'm curious about what this -4.60 value is supposed to represent. Is it a
missing
value code, as Josh inferred, a left-censoring value, or something else? The
reason
for asking is that its purpose could well have an impact on the type of
analysis
that is appropriate for these data:

*  if -4.60 is meant to be a missing value code, its inclusion greatly
inflates the
   degrees of freedom actually present in the data. Moreover, its
presence also
   inflates the variability both within and between groups, which
reduces the
   sensitivity of tests. Also observe from the strip plot above that if
-4.60 is a
   missing value code, then your non-missing data appear to increase in
variance
   with increasing group numbers; moreover, the imbalance in
observations between
   groups is more severe, which in turn makes the p-values of tests less
reliable
   for the non-missing data. (Even worse, the imbalance and increasing
variance
   don't appear to be independent if -4.60 represents a missing value.)
*  if -4.60 is meant to be something like a lower detection limit or
upper bound
   on a left-censored response, then one is artificially reducing
variability within
   groups and the p-values of tests turn out to be optimistic. There are
better ways
   to handle nondetects. One approach is outlined in Helsel's book
'Nondetects and
   Data Analysis', which is the reference of the R package NADA. Other
approaches
   to nondetect data are also available outside of R (e.g., the SCOUT
software at EPA).
   *   if -4.60 is meant to represent a lower bound on a (censored?)
response,
   then setting -4.60 as the response artificially increases the
variation in the
   data. In this case, one would be artificially left-truncating the
data, but for a
   different reason from that given in the immediately preceding point.

Josh also noted the difference in p-values of the between group test when
the groups
were factors as opposed to numeric. This is a common 'gotcha' in R - you
need to pay attention to
the classes of your inputs when fitting any statistical model. In the case
where the
groups comprise a factor, R performs a one-way ANOVA. [FWIW, I got the same
p-value as Josh for your 'complete' data (54 obs.)]. When the group variable
is numeric,
you are fitting a simple linear regression analysis, which implies that the
numeric
values of the groups are meaningful. There is a big difference in the
interpretation of
the two types of models.

HTH,
Dennis


On Tue, Jul 6, 2010 at 6:12 AM, Amit Patel amitrh...@yahoo.co.uk wrote:

 Sorry i had a misprint in the appendix code in the last email


 Hi I needed some help with ANOVA

 I have a problem with My ANOVA
 analysis. I have a dataset with a known ANOVA p-value, however I can
 not seem to re-create it in R.

 I have created a list (zzzanova) which contains
 1)Intensity Values
 2)Group Number (6 Different Groups)
 3)Sample Number (54 different samples)
 this is created by the script in Appendix 1

 I then conduct ANOVA with the command
  zzz.aov - aov(Intensity ~ Group, data = zzzanova)

 I get a p-value of
 Pr(F)1
 0.9483218

 The
 expected p-value is 0.00490 so I feel I maybe using ANOVA incorrectly
 or have put in a wrong formula. I am trying to do an ANOVA analysis
 across all 6 Groups. Is there something wrong with my formula. But I think
 I
 have made a mistake in the formula rather than anything else.




 APPENDIX 1

 datalist - c(-4.60517, -4.60517, -4.60517, -4.60517, -4.60517, -4.60517,
 -4.60517, 3.003749, -4.60517,
2.045314, 2.482557, -4.60517, -4.60517, -4.60517, -4.60517, 1.592743,
 -4.60517,
-4.60517, 0.91328, -4.60517, -4.60517, 1.827744, 2.457795, 0.355075,
 -4.60517, 2.39127,
2.016987, 2.319903, 1.146683, -4.60517, -4.60517, -4.60517, 1.846162,
 -4.60517, 2.121427, 1.973118,
-4.60517, 2.251568, -4.60517, 2.270724, 0.70338, 0.963816, -4.60517,
  0.023703, -4.60517,
2.043382, 1.070586, 

Re: [R] grouped logit regression

2010-07-07 Thread Corey Sparks

by grouped data are you saying that you have counts of outcomes and counts
of trials?
That is how I interpret the glogit in stata.  If that is the case you can
put your data into glm()
like this

fit-glm(nevents~xvars, weights=ntrials, family=binomial, data=yourdataset)
will fit the binomial regression model

summary(fit) 
will print the coefficients and model fit

In the future could you please read the posting guide and put in a data
example or some R code you have tried.
CS

-
Corey Sparks, PhD
Assistant Professor
Department of Demography and Organization Studies
University of Texas at San Antonio
501 West Durango Blvd
Monterey Building 2.270C
San Antonio, TX 78207
210-458-3166
corey.sparks 'at' utsa.edu
https://rowdyspace.utsa.edu/users/ozd504/www/index.htm
-- 
View this message in context: 
http://r.789695.n4.nabble.com/grouped-logit-regression-tp2280763p2280806.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] how to download files from ftp site

2010-07-07 Thread Benoit Wastine

Hi all,
I'm running R 2.10.1 on Windows XP and I'd like to read files from a ftp 
site.

Does anybody know how to do ?
Thanks
Benoit

--
Benoit Wastine
Laboratoire des Sciences du Climat et de l’Environnement (LSCE/IPSL)
CEA-CNRS-UVSQ
CE Saclay
Orme des merisiers
Bât 703 - Pte 13A
91191 Gif sur Yvette Cedex
France

Tel : 33 (0)1 69 08 21 97
Fax : 33 (0)1 69 08 77 16

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Odp: identifying odd or even number

2010-07-07 Thread Petr PIKAL
Hi

r-help-boun...@r-project.org napsal dne 01.07.2010 17:40:22:

 Hi, I run into problem when writing a syntax, I don't know syntax that 
will 
 return true or false if an integer is odd or even.

fff - function(x) as.logical(x%%2)

Regards
Petr



 Thanks
 
 OYEYEMI, Gafar Matanmi
 
 Department of Statistics
 
 University of Ilorin
 
 P.M.B 1515
 
 Ilorin, Kwara State
 
 Nigeria
 
 Tel: +2348052278655
 
 Tel: +2348068241885
 
 
 
[[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] numerical derivative R help

2010-07-07 Thread Setlhare Lekgatlhamang
Which package are you using?

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On Behalf Of Parminder Mankoo
Sent: Tuesday, July 06, 2010 10:50 PM
To: r-help@r-project.org
Subject: [R] numerical derivative R help
Importance: Low

I fit my CDF to sum of exponentials and now I want to take the numerical
derivative of this function to obtain probability density.I will really
appreciate your help reagrding the error messages I am getting which I
don't
understand.

*
*

 fitterma - function(xtime) {

a - -0.09144115

b - -0.01335756

c - -2.368057

d - -0.00600052

return(exp(a+b*xtime)+exp(c+d*xtime))

}


 numericDeriv(fitterma,xtime)


*Error in numericDeriv(fitterma, xtime) : *

*  cannot coerce type 'closure' to vector of type 'double'*

*
*

*Thanks,*

*parmee*

[[alternative HTML version deleted]]

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



DISCLAIMER:\ Sample Disclaimer added in a VBScript.\ ...{{dropped:3}}

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

2010-07-07 Thread Christofer Bogaso
Dear all, I have a date related question. Suppose I have a character string
March-2009, how I can convert it to a valid date object like
as.yearmon(2009-01-03) in the zoo package? Is there any possibility there?

Ans secondly is there any R function which will give the names of of all
months as LETTERS does?

Thanks for your time.

[[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] Need help in handling date

2010-07-07 Thread Whit Armstrong
?strptime

‘%B’ Full month name in the current locale.  (Also matches
  abbreviated name on input.)


On Wed, Jul 7, 2010 at 8:40 AM, Christofer Bogaso
bogaso.christo...@gmail.com wrote:
 Dear all, I have a date related question. Suppose I have a character string
 March-2009, how I can convert it to a valid date object like
 as.yearmon(2009-01-03) in the zoo package? Is there any possibility there?

 Ans secondly is there any R function which will give the names of of all
 months as LETTERS does?

 Thanks for your time.

        [[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] Need help in handling date

2010-07-07 Thread Marc Schwartz
On Jul 7, 2010, at 7:40 AM, Christofer Bogaso wrote:

 Dear all, I have a date related question. Suppose I have a character string
 March-2009, how I can convert it to a valid date object like
 as.yearmon(2009-01-03) in the zoo package? Is there any possibility there?
 
 Ans secondly is there any R function which will give the names of of all
 months as LETTERS does?
 
 Thanks for your time.


You may need to append a default day, so something like:

 as.Date(paste(March-2009, -15, sep = ), format = %B-%Y-%d)
[1] 2009-03-15

Otherwise on OSX, I get:

 as.Date(March-2009, format = %B-%Y)
[1] NA

The above behavior regarding missing components is system specific, as per the 
Note in ?as.Date. Also see ?strptime for the formatting options.

For the second part:

 month.name
 [1] January   February  March April May  
 [6] June  July  AugustSeptember October  
[11] November  December 

See ?LETTERS

HTH,

Marc Schwartz

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


Re: [R] how to download files from ftp site

2010-07-07 Thread Uwe Ligges

?download.file

Uwe Ligges


On 07.07.2010 14:05, Benoit Wastine wrote:

Hi all,
I'm running R 2.10.1 on Windows XP and I'd like to read files from a ftp
site.
Does anybody know how to do ?
Thanks
Benoit



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


Re: [R] Need help in handling date

2010-07-07 Thread Gabor Grothendieck
On Wed, Jul 7, 2010 at 8:40 AM, Christofer Bogaso
bogaso.christo...@gmail.com wrote:
 Dear all, I have a date related question. Suppose I have a character string
 March-2009, how I can convert it to a valid date object like
 as.yearmon(2009-01-03) in the zoo package? Is there any possibility there?

Try this.  The first returns a Date class object and the second
returns a character string:

 d - March-2009
 as.Date(as.yearmon(d, %B-%Y))
[1] 2009-03-01
 format(as.Date(as.yearmon(d, %B-%Y)))
[1] 2009-03-01

See R News 4/1 for more on times and dates.


 Ans secondly is there any R function which will give the names of of all
 months as LETTERS does?

 month.abb
 [1] Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
 month.name
 [1] January   February  March April May   June
 [7] July  AugustSeptember October   November  December

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


[R] can't open file

2010-07-07 Thread Sebastian Kruk
I have a text file log2.log encoded Ansi in Windows.

When I execute:

out - read.zoo(readLines(con - file(log2.log,
encoding=UCS-2LE)),FUN = as.chron)

have errors:

Error en file(file, rt) : no se puede abrir la conexión
Además: Mensajes de aviso perdidos
1: In file(file, rt) :
  sólo fue usado el primer elemento del argumento 'description'
2: In file(file, rt) :
  no fue posible abrir el archivo '#Software: Microsoft Internet
Information Services 5.0': No such file or directory

Why?

Thks,

Sebastián.

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


[R] Loading text file from data folder of package

2010-07-07 Thread Andrew Liu

Hello,

I am creating a package and in my vignette I would like to load a text 
file from the data folder of the package. Currently, I am doing the 
following:


filepath - paste(.libPaths(), pkgname, data, sample.txt, sep = /)
file(filepath)

Is there a better way of doing this?

Thanks,
Andrew

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

2010-07-07 Thread Andrew Leeser
I am looking for a way to trim leading and trailing spaces in a character 
string in R. For example:

   this is random text    

should become:

this is random text.

I have a short function to perform this task as follows:

trim - function(str){
   str - sub(^ +, , str)
   str - sub( +$, , str)
}

While this function does seem to work, I am curious if there's anything built 
into R that I can use instead, as that would be preferable.

Thanks in advance,
Andrew


  
[[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] Trimming in R

2010-07-07 Thread Henrique Dallazuanna
Try this:

gsub(^\\s+|\\s+$, ,this is random text)

On Wed, Jul 7, 2010 at 10:23 AM, Andrew Leeser aml05willi...@yahoo.comwrote:

 I am looking for a way to trim leading and trailing spaces in a character
 string in R. For example:

this is random text

 should become:

 this is random text.

 I have a short function to perform this task as follows:

 trim - function(str){
str - sub(^ +, , str)
str - sub( +$, , str)
 }

 While this function does seem to work, I am curious if there's anything
 built
 into R that I can use instead, as that would be preferable.

 Thanks in advance,
 Andrew



[[alternative HTML version deleted]]


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




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

[[alternative HTML version deleted]]

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


Re: [R] Loading text file from data folder of package

2010-07-07 Thread Henrique Dallazuanna
See

?system.file

On Wed, Jul 7, 2010 at 10:20 AM, Andrew Liu andrew.t@williams.eduwrote:

 Hello,

 I am creating a package and in my vignette I would like to load a text file
 from the data folder of the package. Currently, I am doing the following:

 filepath - paste(.libPaths(), pkgname, data, sample.txt, sep = /)
 file(filepath)

 Is there a better way of doing this?

 Thanks,
 Andrew

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




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

[[alternative HTML version deleted]]

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


Re: [R] can't open file

2010-07-07 Thread Henrique Dallazuanna
You've tried this?

out - read.zoo(log2.log, encoding=UCS-2LE, FUN = as.chron)

On Wed, Jul 7, 2010 at 10:16 AM, Sebastian Kruk residuo.so...@gmail.comwrote:

 I have a text file log2.log encoded Ansi in Windows.

 When I execute:

 out - read.zoo(readLines(con - file(log2.log,
 encoding=UCS-2LE)),FUN = as.chron)

 have errors:

 Error en file(file, rt) : no se puede abrir la conexión
 Además: Mensajes de aviso perdidos
 1: In file(file, rt) :
  sólo fue usado el primer elemento del argumento 'description'
 2: In file(file, rt) :
  no fue posible abrir el archivo '#Software: Microsoft Internet
 Information Services 5.0': No such file or directory

 Why?

 Thks,

 Sebastián.

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




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

[[alternative HTML version deleted]]

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


Re: [R] can't open file

2010-07-07 Thread Gabor Grothendieck
On Wed, Jul 7, 2010 at 9:16 AM, Sebastian Kruk residuo.so...@gmail.com wrote:
 I have a text file log2.log encoded Ansi in Windows.

 When I execute:

 out - read.zoo(readLines(con - file(log2.log,
 encoding=UCS-2LE)),FUN = as.chron)

 have errors:

 Error en file(file, rt) : no se puede abrir la conexión
 Además: Mensajes de aviso perdidos
 1: In file(file, rt) :
  sólo fue usado el primer elemento del argumento 'description'
 2: In file(file, rt) :
  no fue posible abrir el archivo '#Software: Microsoft Internet
 Information Services 5.0': No such file or directory

 Why?


The file argument of read.zoo is a character string giving the *name*
of the file, not the contents of the file as a vector of character
strings.  Alternately it can be a connection (undocumented but works)
or a data.frame so you likely want one of these:

read.zoo(file(log2.log, encoding=UCS-2LE), FUN = as.chron)
read.zoo(log2.log, FUN = as.chron)

See the examples section of ?read.zoo for more examples.

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


Re: [R] can't open file

2010-07-07 Thread Sebastian Kruk
If I edit(out):

structure(logical(0), .Dim = c(707L, 0L), .Dimnames = list(NULL,
NULL), index = structure(1:707, format = m/d/y, origin = structure(c(1,
1, 1970), .Names = c(month, day, year)), class = c(dates,
times)), class = zoo)

I tried:

z - read.zoo(textConnection(L), index = 1:2,FUN = function(x)
paste(x[,1], x[,2]))

I have de following error:

Error en x[, 1] : número incorreto de dimensiones

2010/7/7 Henrique Dallazuanna www...@gmail.com:

 You've tried this?

 out - read.zoo(log2.log, encoding=UCS-2LE, FUN = as.chron)

 On Wed, Jul 7, 2010 at 10:16 AM, Sebastian Kruk residuo.so...@gmail.com
 wrote:

 I have a text file log2.log encoded Ansi in Windows.

 When I execute:

 out - read.zoo(readLines(con - file(log2.log,
 encoding=UCS-2LE)),FUN = as.chron)

 have errors:

 Error en file(file, rt) : no se puede abrir la conexión
 Además: Mensajes de aviso perdidos
 1: In file(file, rt) :
  sólo fue usado el primer elemento del argumento 'description'
 2: In file(file, rt) :
  no fue posible abrir el archivo '#Software: Microsoft Internet
 Information Services 5.0': No such file or directory

 Why?

 Thks,

 Sebastián.

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



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


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


Re: [R] Trimming in R

2010-07-07 Thread David Winsemius


On Jul 7, 2010, at 9:23 AM, Andrew Leeser wrote:

I am looking for a way to trim leading and trailing spaces in a  
character

string in R. For example:

   this is random text

should become:

this is random text.

I have a short function to perform this task as follows:

trim - function(str){
   str - sub(^ +, , str)
   str - sub( +$, , str)
}

While this function does seem to work, I am curious if there's  
anything built

into R that I can use instead, as that would be preferable.


I have been using the trim function in package gdata. The code is  
quite similar to yours.


--

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] can't open file

2010-07-07 Thread Sebastian Kruk
I tried:

L - readLines(con - file(log2.log, encoding=UCS-2LE)
z - read.zoo(textConnection(L), index = 1:2,FUN = function(x)
paste(x[,1], x[,2]))

Error:
Error en x[, 1] : número incorreto de dimensiones

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


Re: [R] can't open file

2010-07-07 Thread Gabor Grothendieck
On Wed, Jul 7, 2010 at 9:43 AM, Sebastian Kruk residuo.so...@gmail.com wrote:
 I tried:

 L - readLines(con - file(log2.log, encoding=UCS-2LE)
 z - read.zoo(textConnection(L), index = 1:2,FUN = function(x)
 paste(x[,1], x[,2]))

 Error:
 Error en x[, 1] : número incorreto de dimensiones


Please provide a reproducible example like this:

library(zoo)
L - 4/4/200 10:10:10 3000
5/5/2000 12:12:12 4000
read.zoo(textConnection(L), index = 1:2, FUN = function(x) paste(x[,1], x[,2]))

which gives me
 4/4/200 10:10:10 5/5/2000 12:12:12
 3000  4000

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


Re: [R] can't open file

2010-07-07 Thread Duncan Murdoch

On 07/07/2010 9:16 AM, Sebastian Kruk wrote:

I have a text file log2.log encoded Ansi in Windows.
  


What Windows calls Ansi is called latin1 in R.  You said the 
encoding was UCS-2LE, which Windows calls Unicode.  Part of your 
problem might be this mismatched encoding.  Have you tried using 
encoding=latin1 when you read the file?


Duncan Murdoch

When I execute:

out - read.zoo(readLines(con - file(log2.log,
encoding=UCS-2LE)),FUN = as.chron)

have errors:

Error en file(file, rt) : no se puede abrir la conexión
Además: Mensajes de aviso perdidos
1: In file(file, rt) :
  sólo fue usado el primer elemento del argumento 'description'
2: In file(file, rt) :
  no fue posible abrir el archivo '#Software: Microsoft Internet
Information Services 5.0': No such file or directory

Why?

Thks,

Sebastián.

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



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

2010-07-07 Thread Kim Vanselow
Dear r-helpers,
I use function step.gam (package gam, T. Hastie) with several explanatory 
variables to build a model. Unfortunately, I obviously have too many variables. 
This message occurs on my 4 core 64bit machine with 8GB RAM in R2.11.1 for 
Windows (64bit build):

Error in array(FALSE, term.lengths) : 'dim' specifies too large an array

I read that this message occurs when running out of RAM. So, the question is:
Could it maybe help to run step.gam as a parallel process on more then 1 core 
and is parallel processing possible with this function?
I tried to manage it with snowfall, but with no success.

Any help would be greatly appreciated!

Thanks,
Kim

This is my step.gam code:

fit0 = gam( Vegetation ~ 1, family = binomial, data = d )
gam.scope = list(
~ 1 + ALTITUDE + s(ALTITUDE, 2),
~ 1 + SLOPE + s(SLOPE, 2),
~ 1 + SOUTH_EXPOSEDNESS + s(SOUTH_EXPOSEDNESS, 2),
~ 1 + WEST_EXPOSEDNESS + s(WEST_EXPOSEDNESS, 2),
~ 1 + VCURV + s(VCURV, 2),
~ 1 + HCURV + s(HCURV, 2),
~ 1 + PRECIPITATION_ANNUAL_MEAN + s(PRECIPITATION_ANNUAL_MEAN, 2),
~ 1 + PRECIPITATION_JULY + s(PRECIPITATION_JULY, 2),
~ 1 + PRECIPITATION_JANUARY + s(PRECIPITATION_JANUARY, 2),
~ 1 + DISTANCE_TO_ISOBATH + s(DISTANCE_TO_ISOBATH, 2),
~ 1 + VERTICAL_DISTANCE_TO_ISOBATH + 
s(VERTICAL_DISTANCE_TO_ISOBATH, 2),
~ 1 + HORIZONTAL_DISTANCE_TO_ISOBATH + 
s(HORIZONTAL_DISTANCE_TO_ISOBATH, 2),
~ 1 + UTM_EASTING + s(UTM_EASTING, 2),
~ 1 + UTM_NORTHING + s(UTM_NORTHING, 2),
~ 1 + DISTANCE_TO_SETTLEMENT + s(DISTANCE_TO_SETTLEMENT, 2),
~ 1 + Bd_1_Mean + s(Bd_1_Mean, 2),
~ 1 + Bd_2_Mean + s(Bd_2_Mean, 2),
~ 1 + Bd_3_Mean + s(Bd_3_Mean, 2),
~ 1 + Bd_4_Mean + s(Bd_4_Mean, 2),
~ 1 + Bd_5_Mean + s(Bd_5_Mean, 2),
~ 1 + Bd_1t_Mean + s(Bd_1t_Mean, 2),
~ 1 + Bd_2t_Mean + s(Bd_2t_Mean, 2),
~ 1 + Bd_3t_Mean + s(Bd_3t_Mean, 2),
~ 1 + Bd_4t_Mean + s(Bd_4t_Mean, 2),
~ 1 + Bd_5t_Mean + s(Bd_5t_Mean, 2),
~ 1 + NDVI_IR_R_Mean + s(NDVI_IR_R_Mean, 2),
~ 1 + NDVI_IR_Rededge_Mean + s(NDVI_IR_Rededge_Mean, 2),
~ 1 + NDVI_Rededge_R_Mean + s(NDVI_Rededge_R_Mean, 2),
~ 1 + ReTCI_Mean + s(ReTCI_Mean, 2),
~ 1 + NDVI_IR_R_t_Mean + s(NDVI_IR_R_t_Mean, 2),
~ 1 + NDVI_IR_Rededge_t_Mean + s(NDVI_IR_Rededge_t_Mean, 2),
~ 1 + NDVI_Rededge_R_t_Mean + s(NDVI_Rededge_R_t_Mean, 2),
~ 1 + ReTCI_t_Mean + s(ReTCI_t_Mean, 2))
gam.scope
fit = step.gam(fit0, scope = gam.scope, direction = both, trace = TRUE)

The error occurs with the function step.gam.

-- 
GMX DSL: Internet-, Telefon- und Handy-Flat ab 19,99 EUR/mtl.  
Bis zu 150 EUR Startguthaben inklusive! http://portal.gmx.net/de/go/dsl

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


Re: [R] Visualization of coefficients

2010-07-07 Thread Michael Friendly

Tal Galili wrote:

Hello David,
Thanks to your posting I started looking at the function in the arm package.
 It appears this function is quite mature, and offers (for example) the
ability to easily overlap coefficients from several models.

I updated the post I published on the subject, so at the end of it I give an
example of comparing the coef of several models:
http://www.r-statistics.com/2010/07/visualization-of-regression-coefficients-in-r/

Thanks again for the pointer.

Best,
Tal


Achim Zeileis wrote:

Re: more mature. arm's coefplot() is more flexible in certain respects, mine is more 
convenient in others. The overlay functionality is something arm's coefplot() is better 
in and it also as some further options (vertical vs. horizontal etc.). My coefplot() has 
the advantage that it does not need any modification as long as coef() and vcov() methods 
are available. Furthermore, level can specify the significance level (instead 
of always using one and two standard errors, respectively).
But it shouldn't be too hard to create a superset of all options.



@Tal:
For the example using library(arm) and the Mroz data, you posted the 
wrong image.  And loose the intercept in the example.


@Achim:
It would be worthwhile combining the generality of your version with the
overlay capability of the arm version, which is extremely useful for 
model comparison.  However, the arm version uses S4 methods.


-Michael


--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University  Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele StreetWeb:   http://www.datavis.ca
Toronto, ONT  M3J 1P3 CANADA

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

2010-07-07 Thread Petr PIKAL
Dear all

My question is more on statistics than on R, however it can be 
demonstrated by R. It is about pros and cons trying to find a relationship 
by aggregated data. I can have two variables which can be related and I 
measure them regularly during some time (let say a year) but I can not 
measure them in a same time - (e.g. I can not measure x and respective 
value of y, usually I have 3 or more values of x and only one value of y 
per day). 

I can make a aggregated values (let say quarterly). My questions are:

1.  Is such approach sound? Can I use it?
2.  What could be the problems
3.  Is there any other method to inspect variables which can be 
related but you can not directly measure them in a same time?

My opinion is, that it is not much sound to inspect aggregated values and 
there can be many traps especially if there are only few aggregated 
values. Below you can see my examples.

If you have some opinion on this issue, please let me know.

Best regards
Petr

Let us have a relation x/y

set.seed(555)
x - rnorm(120)
y - 5*x+3+rnorm(120)
plot(x, y)

As you can see there is clear relation which can be seen from plot. Now I 
make a factor for aggregation.

fac - rep(1:4,each=30)

xprum - tapply(x, fac, mean)
yprum - tapply(y, fac, mean)
plot(xprum, yprum)

Relationship is completely gone. Now let us make other fake data

xn - runif(120)*rep(1:4, each=30)
yn - runif(120)*rep(1:4, each=30)
plot(xn,yn)

There is no visible relation, xn and yn are independent but related to 
aggregation factor.

xprumn - tapply(xn, fac, mean)
yprumn - tapply(yn, fac, mean)
plot(xprumn, yprumn)

Here you can see perfect relation which is only due to aggregation factor.

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


Re: [R] Visualization of coefficients

2010-07-07 Thread Tal Galili
Thanks Michael - I now inserted the correct image for Achim example code.


Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
--




On Wed, Jul 7, 2010 at 5:19 PM, Michael Friendly frien...@yorku.ca wrote:

 Tal Galili wrote:

 Hello David,
 Thanks to your posting I started looking at the function in the arm
 package.
  It appears this function is quite mature, and offers (for example) the
 ability to easily overlap coefficients from several models.

 I updated the post I published on the subject, so at the end of it I give
 an
 example of comparing the coef of several models:

 http://www.r-statistics.com/2010/07/visualization-of-regression-coefficients-in-r/

 Thanks again for the pointer.

 Best,
 Tal


 Achim Zeileis wrote:

 Re: more mature. arm's coefplot() is more flexible in certain respects,
 mine is more convenient in others. The overlay functionality is something
 arm's coefplot() is better in and it also as some further options (vertical
 vs. horizontal etc.). My coefplot() has the advantage that it does not need
 any modification as long as coef() and vcov() methods are available.
 Furthermore, level can specify the significance level (instead of always
 using one and two standard errors, respectively).
 But it shouldn't be too hard to create a superset of all options.



 @Tal:
 For the example using library(arm) and the Mroz data, you posted the wrong
 image.  And loose the intercept in the example.

 @Achim:
 It would be worthwhile combining the generality of your version with the
 overlay capability of the arm version, which is extremely useful for model
 comparison.  However, the arm version uses S4 methods.

 -Michael


 --
 Michael Friendly Email: friendly AT yorku DOT ca
 Professor, Psychology Dept.
 York University  Voice: 416 736-5115 x66249 Fax: 416 736-5814
 4700 Keele StreetWeb:   http://www.datavis.ca
 Toronto, ONT  M3J 1P3 CANADA



[[alternative HTML version deleted]]

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


[R] help with Date class

2010-07-07 Thread Elisenda Vila
I am trying to work with the Date class which is written in S3 and I would
like to access to the elements of the class (for example the year). I've
tryed to do it for example like this:

as.Date(Sys.time)-w
w$year #Doesn't work
w[year] #is NA

I would like to know the correct way to acces to this value.

Thank you so much

-- 
Elisenda Vila

[[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] Different goodness of fit tests leads to contradictory conclusions

2010-07-07 Thread Kiyoshi Sasaki
I am trying to test goodness of fit for my legalistic regression using several 
options as shown below.  Hosmer-Lemeshow test (whose function I borrowed from 
a previous post), Hosmer–le Cessie omnibus lack of fit test (also borrowed 
from a previous post), Pearson chi-square test, and deviance test.  All the 
tests, except the deviance tests, produced p-values well above 0.05.  Would 
anyone please teach me why deviance test produce p-values such a small value 
(0.001687886) that suggest lack of fit, while other tests suggest good fit? Did 
I do something wrong?
 
Thank you for your time and help!
 
Kiyoshi
 
 
mod.fit - glm(formula = no.NA$repcnd ~  no.NA$svl, family = binomial(link = 
logit), data =  no.NA, na.action = na.exclude, control = list(epsilon = 
0.0001, maxit = 50, trace = F))
 
 # Option 1: Hosmer-Lemeshow test
 mod.fit - glm(formula = no.NA$repcnd ~  no.NA$svl, family = binomial(link = 
 logit), data =  no.NA, na.action = na.exclude, control = list(epsilon = 
 0.0001, maxit = 50, trace = F))
             
   hosmerlem - function (y, yhat, g = 10) 
{
cutyhat - cut(yhat, breaks = quantile(yhat, probs = seq(0, 1, 1/g)), 
include.lowest = T)
obs - xtabs(cbind(1 - y, y) ~ cutyhat)
expect - xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
chisq - sum((obs - expect)^2/expect)
P - 1 - pchisq(chisq, g - 2)
c(X^2 = chisq, Df = g - 2, P(Chi) = P)
}
 
 hosmerlem(no.NA$repcnd, fitted(mod.fit))
 X^2                       Df   
                        P(Chi) 
7.8320107            8.000            0.4500497
 
 
 # Option 2 - Hosmer–le Cessie omnibus lack of fit test:
 library(Design)
 lrm.GOF - lrm(formula = no.NA$repcnd ~  no.NA$svl, data =  no.NA, y = T, x 
 = T)
 resid(lrm.GOF,type = gof)
Sum of squared errors     Expected value|H0        
SD                     
Z                     P 
                 48.3487115           
               48.3017399            
              0.1060826     0.4427829     0.6579228
 
 # Option 3 - Pearson chi-square p-value:
 pp - sum(resid(lrm.GOF,typ=pearson)^2)
 1-pchisq(pp, mod.fit$df.resid)
[1] 0.4623282
 
 
 # Option 4 - Deviance (G^2) test:
 1-pchisq(mod.fit$deviance, mod.fit$df.resid)
[1] 0.001687886


  
[[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] can't open file

2010-07-07 Thread Sebastian Kruk
In both cases I obtain:

01/02/70
01/03/70
01/04/70
01/05/70
01/06/70
01/07/70
01/08/70
01/09/70
01/10/70
01/11/70
01/12/70
01/13/70
01/14/70
01/15/70
01/16/70
01/17/70
01/18/70
01/19/70
01/20/70
01/21/70
01/22/70
01/23/70
01/24/70
01/25/70
01/26/70
01/27/70
01/28/70
01/29/70
01/30/70
01/31/70
02/01/70
02/02/70
02/03/70
02/04/70
02/05/70
02/06/70
02/07/70
02/08/70
02/09/70
02/10/70
02/11/70
02/12/70
02/13/70
02/14/70
02/15/70
02/16/70
02/17/70
02/18/70
02/19/70
02/20/70
02/21/70
02/22/70
02/23/70
02/24/70
02/25/70
02/26/70
02/27/70
02/28/70
03/01/70
03/02/70
03/03/70
03/04/70
03/05/70
03/06/70
03/07/70
03/08/70
03/09/70
03/10/70
03/11/70
03/12/70
03/13/70
03/14/70
03/15/70
03/16/70
03/17/70
03/18/70
03/19/70
03/20/70
03/21/70
03/22/70
03/23/70
03/24/70
03/25/70
03/26/70
03/27/70
03/28/70
03/29/70
03/30/70
03/31/70
04/01/70
04/02/70
04/03/70
04/04/70
04/05/70
04/06/70
04/07/70
04/08/70
04/09/70
04/10/70
04/11/70
04/12/70
04/13/70
04/14/70
04/15/70
04/16/70
04/17/70
04/18/70
04/19/70
04/20/70
04/21/70
04/22/70
04/23/70
04/24/70
04/25/70
04/26/70
04/27/70
04/28/70
04/29/70
04/30/70
05/01/70
05/02/70
05/03/70
05/04/70
05/05/70
05/06/70
05/07/70
05/08/70
05/09/70
05/10/70
05/11/70
05/12/70
05/13/70
05/14/70
05/15/70
05/16/70
05/17/70
05/18/70
05/19/70
05/20/70
05/21/70
05/22/70
05/23/70
05/24/70
05/25/70
05/26/70
05/27/70
05/28/70
05/29/70
05/30/70
05/31/70
06/01/70
06/02/70
06/03/70
06/04/70
06/05/70
06/06/70
06/07/70
06/08/70


2010/7/7 Gabor Grothendieck ggrothendi...@gmail.com:
 On Wed, Jul 7, 2010 at 9:16 AM, Sebastian Kruk residuo.so...@gmail.com 
 wrote:
 I have a text file log2.log encoded Ansi in Windows.

 When I execute:

 out - read.zoo(readLines(con - file(log2.log,
 encoding=UCS-2LE)),FUN = as.chron)

 have errors:

 Error en file(file, rt) : no se puede abrir la conexión
 Además: Mensajes de aviso perdidos
 1: In file(file, rt) :
  sólo fue usado el primer elemento del argumento 'description'
 2: In file(file, rt) :
  no fue posible abrir el archivo '#Software: Microsoft Internet
 Information Services 5.0': No such file or directory

 Why?


 The file argument of read.zoo is a character string giving the *name*
 of the file, not the contents of the file as a vector of character
 strings.  Alternately it can be a connection (undocumented but works)
 or a data.frame so you likely want one of these:

 read.zoo(file(log2.log, encoding=UCS-2LE), FUN = as.chron)
 read.zoo(log2.log, FUN = as.chron)

 See the examples section of ?read.zoo for more examples.


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


Re: [R] help with Date class

2010-07-07 Thread David Winsemius


On Jul 7, 2010, at 7:25 AM, Elisenda Vila wrote:

I am trying to work with the Date class which is written in S3 and I  
would
like to access to the elements of the class (for example the year).  
I've

tryed to do it for example like this:

as.Date(Sys.time)-w


Throws an error  ... since Sys.time is a function


w$year #Doesn't work


Why should it? Dates are not stored as lists.

?Date


w[year] #is NA

I would like to know the correct way to acces to this value.


as.Date(Sys.time)-w
Error in as.Date.default(Sys.time) :
  do not know how to convert 'Sys.time' to class Date
 as.Date(Sys.time())-w
 w
[1] 2010-07-07
 format(w, %Y)
[1] 2010

--
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] xyplot: key inside the plot region / lme: confidence bands for predicted

2010-07-07 Thread Michael Friendly
No one replied to my second question: how to get standard errors or 
confidence intervals for the
estimated fixed effects from lme().AFAICS, intervals() only gives 
CIs for coefficients.


My working example is:

library(nlme)
library(lattice)

Ortho - Orthodont
Ortho$year - Ortho$age - 8  # make intercept = initial status

Ortho.mix1 - lme(distance ~ year * Sex, data=Ortho,
   random = ~ 1 + year | Subject, method=ML)
anova(Ortho.mix1)

# get predicted values
grid - expand.grid(year=0:6, Sex=c(Male, Female))
grid$age - grid$year+8  # plot vs. age
fm.mix1 -cbind(grid, distance = predict(Ortho.mix1, newdata = grid, 
level=0))


xyplot(distance ~ age, data=fm.mix1, groups=Sex, type=b,
   par.settings = list(superpose.symbol = list(cex = 1.2, pch=c(15,16))),
   auto.key=list(text=levels(fm.mix1$Sex), points = TRUE, x=0.05, 
y=0.9, corner=c(0,1)),

   main=Linear mixed model: predicted growth)

 intervals(Ortho.mix1)
Approximate 95% confidence intervals

Fixed effects:
  lower  est.  upper
(Intercept)21.607212 22.615625 23.6240383
year0.619660  0.784375  0.9490904
SexFemale  -3.041252 -1.406534  0.2281834
year:SexFemale -0.562889 -0.304830 -0.0467701
attr(,label)
[1] Fixed effects:

Random Effects:
 Level: Subject
  lower   est.upper
sd((Intercept))1.0710615  1.7045122 2.712600
sd(year)   0.0281228  0.1541351 0.844783
cor((Intercept),year) -0.9358505 -0.0311195 0.927652

Within-group standard error:
 lowerest.   upper
1.07081 1.31004 1.60272



Peter Ehlers wrote:

On 2010-07-02 9:37, Michael Friendly wrote:

I have two questions related to plotting predicted values for a linear
mixed model using xyplot:

1: With a groups= argument, I can't seem to get the key to appear inside
the xyplot. (I have the Lattice book,
but don't find an example that actually does this.)
2: With lme(), how can I generate confidence bands or prediction
intervals around the fitted values? Once
I get them, I'd like to add them to the plot.
Example:

library(nlme)
library(lattice)
Ortho - Orthodont
Ortho$year - Ortho$year - 8 # make intercept = initial status

Ortho.mix1 - lme(distance ~ year * Sex, data=Ortho,
random = ~ 1 + year | Subject, method=ML,
# correlation = corAR1 (form = ~ 1 | Subject)
)
Ortho.mix1

# predicted values
grid - expand.grid(year=0:6, Sex=c(Male, Female))
grid$age - grid$year+8 # plot vs. age
fm.mix1 -cbind(grid, distance = predict(Ortho.mix1, newdata = grid,
level=0))

xyplot(distance ~ age, data=fm.mix1, groups=Sex, type=b, pch=c(15,16),
cex=1.2,
auto.key=list(text=c(Male, Female), points = TRUE, x=8, y=26),
main=Linear mixed model: predicted growth)

Can someone help?

-Michael



Michael,

the x,y location should be specified in the unit square
(and you might want to set the 'corner' component).

 Ortho$year - Ortho$age - 8 ##fix typo
 xyplot(distance ~ age, data=fm.mix1, groups=Sex,
type=b, pch=c(15,16), cex=1.2,
auto.key=list(text=c(Male, Female),
  points = TRUE,
  x=0.1, y=0.8, corner=c(0,1)))

See description of 'key' in ?xyplot.

  -Peter Ehlers



--
Michael Friendly Email: friendly AT yorku DOT ca 
Professor, Psychology Dept.

York University  Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele StreetWeb:   http://www.datavis.ca
Toronto, ONT  M3J 1P3 CANADA

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


Re: [R] help with Date class

2010-07-07 Thread Dirk Eddelbuettel
On Wed, Jul 07, 2010 at 01:25:43PM +0200, Elisenda Vila wrote:
 I am trying to work with the Date class which is written in S3 and I would
 like to access to the elements of the class (for example the year). I've
 tryed to do it for example like this:
 
 as.Date(Sys.time)-w

w - Sys.Date()   # does the same, but in one step

 w$year #Doesn't work
 w[year] #is NA
 
 I would like to know the correct way to acces to this value.

wp - as.POSIXlt(w)# POSIXlt has the components
unclass(wp)# shows you all components
wp$year + 1900 # stored as year - 1900, see Unix manuals
wp$mon + 1 # stores as mon -1, see Unix manuals


Arguably, extractor functions would be of help here.

-- 
  Regards, Dirk

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


Re: [R] Different goodness of fit tests leads to contradictory conclusions

2010-07-07 Thread Joris Meys
The first two options are GOF-tests for ungrouped data, the latter two
can only be used for grouped data. According to my experience, the G^2
test is more influenced by this than the X^2 test (gives -wrongly-
significant deviations more easily when used for ungrouped data).

If you started from binary data, you can only use the Hosmer tests.

Cheers
Joris

On Wed, Jul 7, 2010 at 4:00 PM, Kiyoshi Sasaki skiyoshi2...@yahoo.com wrote:
 I am trying to test goodness of fit for my legalistic regression using 
 several options as shown below.  Hosmer-Lemeshow test (whose function I 
 borrowed from a previous post), Hosmer–le Cessie omnibus lack of fit test 
 (also borrowed from a previous post), Pearson chi-square test, and deviance 
 test.  All the tests, except the deviance tests, produced p-values well above 
 0.05.  Would anyone please teach me why deviance test produce p-values such a 
 small value (0.001687886) that suggest lack of fit, while other tests suggest 
 good fit? Did I do something wrong?

 Thank you for your time and help!

 Kiyoshi


 mod.fit - glm(formula = no.NA$repcnd ~  no.NA$svl, family = binomial(link = 
 logit), data =  no.NA, na.action = na.exclude, control = list(epsilon = 
 0.0001, maxit = 50, trace = F))

 # Option 1: Hosmer-Lemeshow test
 mod.fit - glm(formula = no.NA$repcnd ~  no.NA$svl, family = binomial(link = 
 logit), data =  no.NA, na.action = na.exclude, control = list(epsilon = 
 0.0001, maxit = 50, trace = F))

    hosmerlem - function (y, yhat, g = 10)
 {
 cutyhat - cut(yhat, breaks = quantile(yhat, probs = seq(0, 1, 1/g)), 
 include.lowest = T)
 obs - xtabs(cbind(1 - y, y) ~ cutyhat)
 expect - xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
 chisq - sum((obs - expect)^2/expect)
 P - 1 - pchisq(chisq, g - 2)
 c(X^2 = chisq, Df = g - 2, P(Chi) = P)
 }

 hosmerlem(no.NA$repcnd, fitted(mod.fit))
  X^2       Df       P(Chi)
 7.8320107    8.000    0.4500497


 # Option 2 - Hosmer–le Cessie omnibus lack of fit test:
 library(Design)
 lrm.GOF - lrm(formula = no.NA$repcnd ~  no.NA$svl, data =  no.NA, y = T, x 
 = T)
 resid(lrm.GOF,type = gof)
 Sum of squared errors Expected value|H0    SD 
 Z P
      48.3487115      48.3017399       
    0.1060826 0.4427829 0.6579228

 # Option 3 - Pearson chi-square p-value:
 pp - sum(resid(lrm.GOF,typ=pearson)^2)
 1-pchisq(pp, mod.fit$df.resid)
 [1] 0.4623282


 # Option 4 - Deviance (G^2) test:
 1-pchisq(mod.fit$deviance, mod.fit$df.resid)
 [1] 0.001687886



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





-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
joris.m...@ugent.be
---
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php

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


Re: [R] Gray level mosaic plot with shading_Friendly

2010-07-07 Thread Michael Kubovy
Dear Achim and Michael,

Thank you so much. Indeed, mosaic(Titanic, gp = shading_hcl, gp_args = list(lty 
= 1:2, c = 0)) does almost what I was looking for, except that for consistency 
and clarity, I would have expected the negative values on the legend to be be 
outlined with lty = 2.

Michael


On Jul 7, 2010, at 2:13 AM, Achim Zeileis wrote:

 On Tue, 6 Jul 2010, Michael Friendly wrote:
 
 Michael Kubovy wrote:
 Suppose we start with
 data(Titanic)
 mosaic(Titanic, shade = TRUE)
 How do I combine the dashed box contours of shading_Friendly to indicate 
 negative residuals, with three levels of gray: dark for abs(Pearson Resid) 
  4, lighter for 4  abs(Pearson Resid)  2, and lightest for bs(Pearson 
 Resid)  2 ?
 
 Do you mean [1] you want to plot positive residuals in color and negative in 
 gray scale?
 Or [2] to fold + and - residuals by shading all according to abs(resid), and
 distinguishing + from - by the dashed box outlines?
 
 In fact, I designed this coding scheme so that mosaic plots in color (with 
 my blue - white - red scheme) would approximately do exactly what
 you might want under [2], when rendered in B/W, since the fully saturated 
 red and blue are close in  darkness in B/W.
 
 And shading_hcl() has been written to do exactly what you want under [2]. 
 While it is hard to come up with colors of different hues in HSV or HLS space 
 that have the same brightness (aka lightness/luminance) and the same
 colorfulness (aka chroma), this is easy in HCL.
 
 Try
 mosaic(Titanic, gp=shading_Friendly)
 save as a jpg/png and try converting to B/W with an image program and see if 
 this is good enough.
 
 mosaic(Titanic, shade = TRUE)
 
 is the same as
 
 mosaic(Titanic, gp = shading_hcl)
 
 which you can then modify to have different line types
 
 mosaic(Titanic, gp = shading_hcl, gp_args = list(lty = 1:2))
 
 If you print that on a grayscale printer you will see the same plot without 
 any chroma, i.e.,
 
 mosaic(Titanic, gp = shading_hcl, gp_args = list(lty = 1:2, c = 0))
 
 The shading_hcl() function is introduced in Zeileis et al. (2007, JCGS), see 
 ?shading_hcl, which provides more detailed references to HCL colors etc.
 
 Best,
 Z
 
 Alternatively, write your own, shading_Kubovy, modeled on
 
 shading_Friendly -
 function (observed = NULL, residuals = NULL, expected = NULL,
   df = NULL, h = c(2/3, 0), lty = 1:2, interpolate = c(2, 4),
   eps = 0.01, line_col = black, ...)
 {
   shading_hsv(observed = NULL, residuals = NULL, expected = NULL,
   df = NULL, h = h, v = 1, lty = lty, interpolate = interpolate,
   eps = eps, line_col = line_col, p.value = NA, ...)
 }
 environment: namespace:vcd
 attr(,class)
 [1] grapcon_generator
 
 In the defaults, lty=1:2 is what distinguishes + and - for outline line type
 
 hope this helps,
 -Michael
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] interaction post hoc/ lme repeated measures

2010-07-07 Thread Doug Bourne
Hi Everyone,

I’m trying to figure out how to get R to analyze this experiment properly.  I 
have a series of subjects each with two legs.  Within each leg there are two 
bones that I am interested in.  There are also two treatments that I am 
interested in.  That results in four different combinations of treatments.  
Obviously, since the subjects only have two legs, they can’t receive each 
combination of treatment.  The groups are unbalanced so I can’t use aov.  As I 
understand it, lme should work but I am having a tough time figuring out what 
to use as a random term.  Is this correct?

x.lme - lme(effect~bone*treatment1*treatment2, random= ~1|subject, data=x)

anova(x.lme)

I’ve tried a number of different random terms and get very similar results (and 
even get similar results running aov) but I’d like to know what the proper way 
of doing it is.

My next question is how to do post hoc tests on this.  One website recommends 
something like this from multcomp: 
summary(glht(am2,linfct=mcp(myfactor=Tukey)))

But apparently this doesn’t work with interactions.  With my data, I get a 
significant three-way interaction.  If I run TukeyHSD after an aov (ignoring 
the repeated measures) I would get a series of p values for all the different 
combinations of treatments and bones.  Is there a way to do this with lme?

Thank you in advance for any help.

Doug Bourne
[[alternative HTML version deleted]]

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


Re: [R] how to display the clock time in the loop

2010-07-07 Thread Jack Luo
Thanks a bunch, it works.

On Thu, Jul 1, 2010 at 4:50 PM, Nikhil Kaza nikhil.l...@gmail.com wrote:

 explicit call to print usually works for me.

 library(audio)
 for (i in 1:5){
 wait(60)
 print(Sys.time())
 }


 On Jul 1, 2010, at 4:30 PM, Matt Shotwell wrote:

 Try to flush output after printing:

 cat(paste(Sys.time()),\n); flush(stdout())

 On Thu, 2010-07-01 at 16:17 -0400, Jack Luo wrote:

 Hi,


 I am doing some computation which is pretty time consuming, I want R to

 display CPU time after each iteration using the command Sys.time().
 However,

 I found that the code only began to display the CPU time after quite a
 while

 and several iterations have finished. Is there a way to ask R to display

 time right after each iteration is finished?


 Thanks,


 -Jun


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

 --
 Matthew S. Shotwell
 Graduate Student
 Division of Biostatistics and Epidemiology
 Medical University of South Carolina
 http://biostatmatt.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] F# vs. R

2010-07-07 Thread Sergey Goriatchev
Hello, everyone

F# is now public. Compiled code should run  faster than R.

Anyone has opinion on F# vs. R? Just curious

Best,
S

--
---
Kniven skärpes bara mot stenen.

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


Re: [R] relation in aggregated data

2010-07-07 Thread Joris Meys
You examples are pretty extreme... Combining 120 data points in 4
points is off course never going to give a result. Try :

fac - rep(1:8,each=15)
xprum - tapply(x, fac, mean)
yprum - tapply(y, fac, mean)
plot(xprum, yprum)

Relation is not obvious, but visible.

Yes, you lose information. Yes, your hypothesis changes. But in the
case you describe, averaging the x-values for every day (so you get an
average linked to 1 y value) seems like a possibility, given you take
that into account when formulating the hypothesis. Optimally, you
should take the standard error on the average into account for the
analysis, but this is complicated, often not done and in most cases
ignoring this issue is not influencing the results to that extent it
becomes important.

YMMV

Cheers

On Wed, Jul 7, 2010 at 4:24 PM, Petr PIKAL petr.pi...@precheza.cz wrote:
 Dear all

 My question is more on statistics than on R, however it can be
 demonstrated by R. It is about pros and cons trying to find a relationship
 by aggregated data. I can have two variables which can be related and I
 measure them regularly during some time (let say a year) but I can not
 measure them in a same time - (e.g. I can not measure x and respective
 value of y, usually I have 3 or more values of x and only one value of y
 per day).

 I can make a aggregated values (let say quarterly). My questions are:

 1.      Is such approach sound? Can I use it?
 2.      What could be the problems
 3.      Is there any other method to inspect variables which can be
 related but you can not directly measure them in a same time?

 My opinion is, that it is not much sound to inspect aggregated values and
 there can be many traps especially if there are only few aggregated
 values. Below you can see my examples.

 If you have some opinion on this issue, please let me know.

 Best regards
 Petr

 Let us have a relation x/y

 set.seed(555)
 x - rnorm(120)
 y - 5*x+3+rnorm(120)
 plot(x, y)

 As you can see there is clear relation which can be seen from plot. Now I
 make a factor for aggregation.

 fac - rep(1:4,each=30)

 xprum - tapply(x, fac, mean)
 yprum - tapply(y, fac, mean)
 plot(xprum, yprum)

 Relationship is completely gone. Now let us make other fake data

 xn - runif(120)*rep(1:4, each=30)
 yn - runif(120)*rep(1:4, each=30)
 plot(xn,yn)

 There is no visible relation, xn and yn are independent but related to
 aggregation factor.

 xprumn - tapply(xn, fac, mean)
 yprumn - tapply(yn, fac, mean)
 plot(xprumn, yprumn)

 Here you can see perfect relation which is only due to aggregation factor.

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




-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
joris.m...@ugent.be
---
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php

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

2010-07-07 Thread guox
We want to sum many vectors and numbers together as a vector if there is
at least one vector in the arguments.
For example, 1 + c(2,3) = c(3,4).
Since we are not sure arguments to sum, we are using sum function:
sum(v1,v2,...,n1,n2,..).
The problem is that sum returns the sum of all the values present in its
arguments:
sum(1,c(2,3))=6
sum(1,2,3)=6
We do not want to turn sum(v1,v2,...,n1,n2,..) to v1+v2+...+n1+n2+...
So do you know easy way to sum vectors (v1,v2,...) and numbers (n1,n2,...)
as a vector without using v1+v2+...+n1+n2+...? Thanks,

-james

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


Re: [R] Sum vectors and numbers

2010-07-07 Thread Erik Iverson



g...@ucalgary.ca wrote:

We want to sum many vectors and numbers together as a vector if there is
at least one vector in the arguments.
For example, 1 + c(2,3) = c(3,4).
Since we are not sure arguments to sum, we are using sum function:
sum(v1,v2,...,n1,n2,..).
The problem is that sum returns the sum of all the values present in its
arguments:
sum(1,c(2,3))=6
sum(1,2,3)=6
We do not want to turn sum(v1,v2,...,n1,n2,..) to v1+v2+...+n1+n2+...
So do you know easy way to sum vectors (v1,v2,...) and numbers (n1,n2,...)
as a vector without using v1+v2+...+n1+n2+...? Thanks,


 I must admit I've read this a couple times and am confused.  Can you 
give one example that shows exactly what you want, using R objects?


I guess my question is, what is wrong with 1 + c(2, 3)

Is this what you're looking for?

lapply(c(1,2,3), `+`, c(8,9,10))

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


Re: [R] Sum vectors and numbers

2010-07-07 Thread Gabor Grothendieck
On Wed, Jul 7, 2010 at 11:35 AM,  g...@ucalgary.ca wrote:
 We want to sum many vectors and numbers together as a vector if there is
 at least one vector in the arguments.
 For example, 1 + c(2,3) = c(3,4).
 Since we are not sure arguments to sum, we are using sum function:
 sum(v1,v2,...,n1,n2,..).
 The problem is that sum returns the sum of all the values present in its
 arguments:
 sum(1,c(2,3))=6
 sum(1,2,3)=6
 We do not want to turn sum(v1,v2,...,n1,n2,..) to v1+v2+...+n1+n2+...
 So do you know easy way to sum vectors (v1,v2,...) and numbers (n1,n2,...)
 as a vector without using v1+v2+...+n1+n2+...? Thanks,

Try this:


 L - list(1:3, 3, 4:6, 6)
 Reduce(+, L)
[1] 14 16 18

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


Re: [R] F# vs. R

2010-07-07 Thread Marc Schwartz
On Jul 7, 2010, at 10:31 AM, Sergey Goriatchev wrote:

 Hello, everyone
 
 F# is now public. Compiled code should run  faster than R.
 
 Anyone has opinion on F# vs. R? Just curious
 
 Best,
 S


The key time critical parts of R are written in compiled C and FORTRAN.

Of course, if you want to take the time to code and validate a Cox PH or mixed 
effects model in F# and then run them against R's coxph() or lme()/lmer() 
functions to test the timing, feel free...  :-)

So unless there is a pre-existing library of statistical and related 
functionality for F#, perhaps you need to reconsider your query.

Regards,

Marc Schwartz

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


[R] Appropriateness of survdiff {survival} for non-censored data

2010-07-07 Thread Frostygoat
I read through Harrington and Fleming (1982) but it is beyond my
statistical comprehension.  I have survival data for insects that have
a very finite expiration date.  I'm trying to test for differences in
survival distributions between different groups.  I understand that
the medical field is most often dealing with censored data and that
survival analysis, at least in the package survival, is largely built
around these conventions and differs from a classical biological
perspective.  For example, for lifetable analysis of insects there is
often no need to estimate survival using a Kaplan-Meir estimate
because it is relatively easy to follow a cohort of individuals
through the entire course of life.  Thus I question the
appropriateness of using survdiff in my analysis; I have exact data
yet I would be testing on the Kaplan-Meir estimate of these data in
survdiff.  Thanks for any help.

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


Re: [R] Trimming in R

2010-07-07 Thread Ken Takagi
David Winsemius dwinsemius at comcast.net writes:

 
 
 On Jul 7, 2010, at 9:23 AM, Andrew Leeser wrote:
 
  I am looking for a way to trim leading and trailing spaces in a  
  character
  string in R. For example:
 
 this is random text
 
  should become:
 
  this is random text.
 
  I have a short function to perform this task as follows:
 
  trim - function(str){
 str - sub(^ +, , str)
 str - sub( +$, , str)
  }
 
  While this function does seem to work, I am curious if there's  
  anything built
  into R that I can use instead, as that would be preferable.
 
 I have been using the trim function in package gdata. The code is  
 quite similar to yours.
 

You could also use the substr() function, assuming that you know a priori how
many spaces are before and after the text:

str =this is random text
nchar(str)
substr(str, 8, 26) # 7 spaces before and 8 after the text you want to keep.

HTH,
Ken

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


Re: [R] Different goodness of fit tests leads to contradictory conclusions

2010-07-07 Thread Kiyoshi Sasaki
Dear Joris,
 
Thank you for your prompt reply! I have a binary dependent variable (whether a 
snake is pregnant or not pregnant). Independent/predictor variable is the 
snake's body size. Each observation (row) of the data represents each snake. 
One column of the data contain '0' or '1' to indicate whether a snake is 
pregnant. Another column contain body size for each snake. So, if I understand 
correctly, I can use only Hosmer-Lemeshow test. Am I correct?
 
Thank you very much!
 
Kiyoshi

--- On Wed, 7/7/10, Joris Meys jorism...@gmail.com wrote:


From: Joris Meys jorism...@gmail.com
Subject: Re: [R] Different goodness of fit tests leads to contradictory 
conclusions
To: Kiyoshi Sasaki skiyoshi2...@yahoo.com
Cc: r-help@r-project.org
Date: Wednesday, July 7, 2010, 10:08 AM


The first two options are GOF-tests for ungrouped data, the latter two
can only be used for grouped data. According to my experience, the G^2
test is more influenced by this than the X^2 test (gives -wrongly-
significant deviations more easily when used for ungrouped data).

If you started from binary data, you can only use the Hosmer tests.

Cheers
Joris

On Wed, Jul 7, 2010 at 4:00 PM, Kiyoshi Sasaki skiyoshi2...@yahoo.com wrote:
 I am trying to test goodness of fit for my legalistic regression using 
 several options as shown below.  Hosmer-Lemeshow test (whose function I 
 borrowed from a previous post), Hosmer–le Cessie omnibus lack of fit test 
 (also borrowed from a previous post), Pearson chi-square test, and deviance 
 test.  All the tests, except the deviance tests, produced p-values well 
 above 0.05.  Would anyone please teach me why deviance test produce p-values 
 such a small value (0.001687886) that suggest lack of fit, while other tests 
 suggest good fit? Did I do something wrong?

 Thank you for your time and help!

 Kiyoshi


 mod.fit - glm(formula = no.NA$repcnd ~  no.NA$svl, family = binomial(link = 
 logit), data =  no.NA, na.action = na.exclude, control = list(epsilon = 
 0.0001, maxit = 50, trace = F))

 # Option 1: Hosmer-Lemeshow test
 mod.fit - glm(formula = no.NA$repcnd ~  no.NA$svl, family = binomial(link 
 = logit), data =  no.NA, na.action = na.exclude, control = list(epsilon = 
 0.0001, maxit = 50, trace = F))

    hosmerlem - function (y, yhat, g = 10)
 {
 cutyhat - cut(yhat, breaks = quantile(yhat, probs = seq(0, 1, 1/g)), 
 include.lowest = T)
 obs - xtabs(cbind(1 - y, y) ~ cutyhat)
 expect - xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
 chisq - sum((obs - expect)^2/expect)
 P - 1 - pchisq(chisq, g - 2)
 c(X^2 = chisq, Df = g - 2, P(Chi) = P)
 }

 hosmerlem(no.NA$repcnd, fitted(mod.fit))
  X^2                       Df   
                         P(Chi)
 7.8320107            8.000            0.4500497


 # Option 2 - Hosmer–le Cessie omnibus lack of fit test:
 library(Design)
 lrm.GOF - lrm(formula = no.NA$repcnd ~  no.NA$svl, data =  no.NA, y = T, 
 x = T)
 resid(lrm.GOF,type = gof)
 Sum of squared errors     Expected value|H0        
 SD                     
 Z                     P
                  48.3487115           
                48.3017399            
               0.1060826     0.4427829     0.6579228

 # Option 3 - Pearson chi-square p-value:
 pp - sum(resid(lrm.GOF,typ=pearson)^2)
 1-pchisq(pp, mod.fit$df.resid)
 [1] 0.4623282


 # Option 4 - Deviance (G^2) test:
 1-pchisq(mod.fit$deviance, mod.fit$df.resid)
 [1] 0.001687886



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





-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
joris.m...@ugent.be
---
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php



  
[[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] grayscale wireframe??

2010-07-07 Thread Sebastien Guyader

With grDevices package, I do the following to generate a greyscale:

  newcols - colorRampPalette(c(white, black)) #generates palette from
white to black
#OR
  newcols - colorRampPalette(c(grey90, grey10)) #generates palette
frome light to dark grey for better visibility

Then in the wireframe() command type the argument: col.regions=newcols(100)
-- 
View this message in context: 
http://r.789695.n4.nabble.com/grayscale-wireframe-tp2280281p2281127.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] how to define a function in R

2010-07-07 Thread Peter Ehlers

Chapter 10 of 'An Introduction to R' has the
suggestive title Writing your own functions.

Isn't that the *first* place you would look?

  -Peter Ehlers

On 2010-07-06 20:57, jd6688 wrote:


Thanks Joshua for the example which has been great help me to start.

Wish you the best
Jason Ding

On Tue, Jul 6, 2010 at 9:23 PM, Joshua Wiley-2 [via R]
ml-node+2280373-448579502-312...@n4.nabble.comml-node%2b2280373-448579502-312...@n4.nabble.com

wrote:



Hello,

As others have said, its hard to give specific advice without specific
needs, but that's okay; I made up some examples needs and some (rather
silly) code that might handle it.  Depending what you need to do, it
may help you get started.  I tried to explicitly name all the
arguments in any functions I used.

When I make gmail use basic text format instead of html, code is sent
poorly, so you can trundle off here to see the example, if you like.

http://gist.github.com/466164

Cheers,

Josh

On Tue, Jul 6, 2010 at 3:48 PM, jd6688[hidden 
email]http://user/SendEmail.jtp?type=nodenode=2280373i=0
wrote:



1. how to write a R script?
2.How to write a SAS like macro/generic process to process multiple files

by

using the same funstion in R?

Thanks in advance
--
View this message in context:

http://r.789695.n4.nabble.com/how-to-define-a-function-in-R-tp2280290p2280290.htmlhttp://r.789695.n4.nabble.com/how-to-define-a-function-in-R-tp2280290p2280290.html?by-user=t

Sent from the R help mailing list archive at Nabble.com.

__
[hidden email]http://user/SendEmail.jtp?type=nodenode=2280373i=1mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmlhttp://www.r-project.org/posting-guide.html

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





--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

__
[hidden email]http://user/SendEmail.jtp?type=nodenode=2280373i=2mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.htmlhttp://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] Adding two files into one and vlookup

2010-07-07 Thread raghu

I dont get any putput by merging.
merge(file1,file2)
[1] Date  Price
0 rows (or 0-length row.names)

file3-merge(file1,file2,by.file1='Date')
 file3
[1] Date  Price
0 rows (or 0-length row.names)

The above are the results. Nothing is coming in file3.

On 7/7/10, Peter Alspach-2 [via R]
ml-node+2280305-1773552212-309...@n4.nabble.com wrote:



 Tena koe

 Why is merge() not helpful?  From your description I would imagine

 merge(file1, file2, by='Date')

 would do what you require.

 HTH 

 Peter Alspach

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of raghu
 Sent: Wednesday, 7 July 2010 10:18 a.m.
 To: r-help@r-project.org
 Subject: [R] Adding two files into one and vlookup


 I have two files with dates and prices in each. The number of rows in
 each of
 them will differ. How do I create a new file which contains data from
 both
 these files? Cbind and merge are not helpful. For cbind because the
 rows are
 not the same replication occurs. Also if I have similar data how do I
 write
 a vlookup kind of function? I am giving an example below:
 Say Price1 file contains the following:
 Date Price
 2/3/2010   134.00
 3/3/2010   133.90
 4/3/2010   135.55

 And say price2 contains the following:
 Date  Price
 2/3/20102300
 3/3/20103200
 4/3/20101800
 5/3/20101900

 I want to take both these data together in a single file, and take the
 smaller vector (or matrix or dataframe??..i am new to R and still
 confused
 with the various objects) which is file1 (because it contains fewer
 rows )
 and vlookup prices in the second file basedon the dates on file1 and
 write
 three columns (date, price from 1 and price from2) in a new file. How
 do i
 do this please?

 Many thanks...
 R
 --
 View this message in context: http://r.789695.n4.nabble.com/Adding-two-
 files-into-one-and-vlookup-tp2280277p2280277.html
 Sent from the R help mailing list archive at Nabble.com.

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

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


 __
 View message @
 http://r.789695.n4.nabble.com/Adding-two-files-into-one-and-vlookup-tp2280277p2280305.html
 To start a new topic under R help, email
 ml-node+789696-608741344-309...@n4.nabble.com
 To unsubscribe from R help, click
  (link removed) 



-- 
'Raghu'

-- 
View this message in context: 
http://r.789695.n4.nabble.com/Adding-two-files-into-one-and-vlookup-tp2280277p2281026.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.


Re: [R] Adding two files into one and vlookup

2010-07-07 Thread Raghu
I was able to achieve the desired output. The issue I had was that
both the files contained the Date column in the data format. So, i
first converted this colum alone into numeric and then did a merge on
the files and it worked. The moral is that in any dataframe all
variable have to be of similar type and then merge works. Correct me
if i am wrong?

On 7/7/10, Raghu r.raghura...@gmail.com wrote:
 Say I have two files file and file2:

 file1 contains the following:
 Date  Price
 02/07/201053.96597903
 03/07/201056.92825807
 04/07/201039.27408645
 05/07/201042.59834151
 06/07/201070.68512383
 07/07/201010.92505265
 08/07/201052.12492249
 09/07/201049.88767957

 file2 contains the following:
 Date  Price
 03/07/20105.312006403
 04/07/2010673.0705924
 05/07/2010442.4679386
 06/07/2010851.9158985
 07/07/2010581.8592424

 I want to create a new file that should look like:
 Date  Price1Price2
 03/07/20105.31200640356.928
 04/07/2010673.070592439.274
 05/07/2010442.467938642.598
 06/07/2010851.915898570.685
 07/07/2010581.859242410.925

 Thx



 On 7/7/10, Erik Iverson er...@ccbr.umn.edu wrote:
 raghu wrote:
 I have two files with dates and prices in each. The number of rows in
 each
 of
 them will differ. How do I create a new file which contains data from
 both
 these files? Cbind and merge are not helpful. For cbind because the rows
 are
 not the same replication occurs. Also if I have similar data how do I
 write
 a vlookup kind of function? I am giving an example below:
 Say Price1 file contains the following:
 Date Price
 2/3/2010   134.00
 3/3/2010   133.90
 4/3/2010   135.55

 And say price2 contains the following:
 Date  Price
 2/3/20102300
 3/3/20103200
 4/3/20101800
 5/3/20101900

 I want to take both these data together in a single file, and take the
 smaller vector (or matrix or dataframe??..i am new to R and still
 confused
 with the various objects) which is file1 (because it contains fewer rows
 )
 and vlookup prices in the second file basedon the dates on file1 and
 write
 three columns (date, price from 1 and price from2) in a new file. How do
 i
 do this please?

 I think all this can be accomplished with merge.  Can you give
 reproducible
 examples as the posting guide suggests?

 Use read.table to read in your data into R objects, then use ?dput to
 give
 us
 the exact copies of the objects (probably data.frames by your example),
 and
 what
 output you want to have.  Being precise with the classes of objects
 you're
 working with is key, and ?dput is a great way to make sure we have the
 same
 objects as you.

 Another tip is common terminology. For instance, `vlookup` is not a term
 used in
 R, and many people will not know what it means.

 This way, everything is reproducible for us, and we can offer suggestions
 and
 show you what the exact output will be.  In short, making sure everyone
 is
 on
 the same page goes a long way when getting help from a mailing list.




 --
 'Raghu'



-- 
'Raghu'

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


Re: [R] Adding two files into one and vlookup

2010-07-07 Thread raghu

I tried and this error popped:
z1 - read.zoo(textConnection(Lines1), header = TRUE, format = fmt)
Error in textConnection(Lines1) : invalid 'text' argument

Also library(chron) returns error as it is an invalid library. Could you help?

On 7/7/10, Gabor Grothendieck [via R]
ml-node+2280348-527941373-309...@n4.nabble.com wrote:



 On Tue, Jul 6, 2010 at 8:26 PM, Gabor Grothendieck
 ggrothendi...@gmail.com wrote:
 On Tue, Jul 6, 2010 at 6:18 PM, raghu r.raghura...@gmail.com wrote:

 I have two files with dates and prices in each. The number of rows in
 each of
 them will differ. How do I create a new file which contains data from
 both
 these files? Cbind and merge are not helpful. For cbind because the rows
 are
 not the same replication occurs. Also if I have similar data how do I
 write
 a vlookup kind of function? I am giving an example below:
 Say Price1 file contains the following:
 Date Price
 2/3/2010   134.00
 3/3/2010   133.90
 4/3/2010   135.55

 And say price2 contains the following:
 Date  Price
 2/3/20102300
 3/3/20103200
 4/3/20101800
 5/3/20101900

 I want to take both these data together in a single file, and take the
 smaller vector (or matrix or dataframe??..i am new to R and still
 confused
 with the various objects) which is file1 (because it contains fewer rows
 )
 and vlookup prices in the second file basedon the dates on file1 and
 write
 three columns (date, price from 1 and price from2) in a new file. How do
 i
 do this please?


 Try this and for more read the three vignettes (pdf documents) in the
 zoo package and also read the R News 4/1 article on dates and times:

 Lines1 - Date Price
 2/3/2010   134.00
 3/3/2010   133.90
 4/3/2010   135.55

 Lines2 - Date  Price
 2/3/20102300
 3/3/20103200
 4/3/20101800
 5/3/20101900

 library(zoo)
 library(chron)
 z1 - read.zoo(textConnection(Lines1), header = TRUE, FUN = chron)
 z2 - read.zoo(textConnection(Lines2), header = TRUE, FUN = chron)

 I originally assumed that the dates were the usual month/day/year but
 looking at it again I suspect they are day/month/year so lets use Date
 class instead of chron replacing the last three statements with:

 fmt - %d/%m/%Y
 z1 - read.zoo(textConnection(Lines1), header = TRUE, format = fmt)
 z2 - read.zoo(textConnection(Lines2), header = TRUE, format = fmt)





 merge(z1, z2) # keep all rows in each
 merge(z1, z2, all = FALSE) # keep only rows in both
 merge(z1, z2, all = c(TRUE, FALSE)) # keep all rows in z1
 merge(z1, z2, all = c(FALSE, TRUE)) # keep all rows in z2


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


 __
 View message @
 http://r.789695.n4.nabble.com/Adding-two-files-into-one-and-vlookup-tp2280277p2280348.html
 To start a new topic under R help, email
 ml-node+789696-608741344-309...@n4.nabble.com
 To unsubscribe from R help, click
  (link removed) 



-- 
'Raghu'

-- 
View this message in context: 
http://r.789695.n4.nabble.com/Adding-two-files-into-one-and-vlookup-tp2280277p2281017.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.


Re: [R] Adding two files into one and vlookup

2010-07-07 Thread Raghu
Say I have two files file and file2:

file1 contains the following:
DatePrice
02/07/2010  53.96597903
03/07/2010  56.92825807
04/07/2010  39.27408645
05/07/2010  42.59834151
06/07/2010  70.68512383
07/07/2010  10.92505265
08/07/2010  52.12492249
09/07/2010  49.88767957

file2 contains the following:
DatePrice
03/07/2010  5.312006403
04/07/2010  673.0705924
05/07/2010  442.4679386
06/07/2010  851.9158985
07/07/2010  581.8592424

I want to create a new file that should look like:
Date  Price1Price2
03/07/2010  5.31200640356.928
04/07/2010  673.070592439.274
05/07/2010  442.467938642.598
06/07/2010  851.915898570.685
07/07/2010  581.859242410.925

Thx



On 7/7/10, Erik Iverson er...@ccbr.umn.edu wrote:
 raghu wrote:
 I have two files with dates and prices in each. The number of rows in each
 of
 them will differ. How do I create a new file which contains data from both
 these files? Cbind and merge are not helpful. For cbind because the rows
 are
 not the same replication occurs. Also if I have similar data how do I
 write
 a vlookup kind of function? I am giving an example below:
 Say Price1 file contains the following:
 Date Price
 2/3/2010   134.00
 3/3/2010   133.90
 4/3/2010   135.55

 And say price2 contains the following:
 Date  Price
 2/3/20102300
 3/3/20103200
 4/3/20101800
 5/3/20101900

 I want to take both these data together in a single file, and take the
 smaller vector (or matrix or dataframe??..i am new to R and still confused
 with the various objects) which is file1 (because it contains fewer rows )
 and vlookup prices in the second file basedon the dates on file1 and write
 three columns (date, price from 1 and price from2) in a new file. How do i
 do this please?

 I think all this can be accomplished with merge.  Can you give reproducible
 examples as the posting guide suggests?

 Use read.table to read in your data into R objects, then use ?dput to give
 us
 the exact copies of the objects (probably data.frames by your example), and
 what
 output you want to have.  Being precise with the classes of objects you're
 working with is key, and ?dput is a great way to make sure we have the same
 objects as you.

 Another tip is common terminology. For instance, `vlookup` is not a term
 used in
 R, and many people will not know what it means.

 This way, everything is reproducible for us, and we can offer suggestions
 and
 show you what the exact output will be.  In short, making sure everyone is
 on
 the same page goes a long way when getting help from a mailing list.




-- 
'Raghu'

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


Re: [R] Gray level mosaic plot with shading_Friendly

2010-07-07 Thread Achim Zeileis



On Wed, 7 Jul 2010, Michael Kubovy wrote:


Dear Achim and Michael,

Thank you so much. Indeed, mosaic(Titanic, gp = shading_hcl, gp_args = 
list(lty = 1:2, c = 0)) does almost what I was looking for, except that 
for consistency and clarity, I would have expected the negative values 
on the legend to be be outlined with lty = 2.


In the continuous legend, that is employed by default (legend_resbased), 
it is visually not very compelling to show line types as well. But you can 
set legend = legend_fixed which displays this information (but is less 
intuitive concerning the interval ranges).


Best,
Z


Michael


On Jul 7, 2010, at 2:13 AM, Achim Zeileis wrote:


On Tue, 6 Jul 2010, Michael Friendly wrote:


Michael Kubovy wrote:

Suppose we start with
data(Titanic)
mosaic(Titanic, shade = TRUE)
How do I combine the dashed box contours of shading_Friendly to indicate negative 
residuals, with three levels of gray: dark for abs(Pearson Resid)  4, lighter for 4 
 abs(Pearson Resid)  2, and lightest for bs(Pearson Resid)  2 ?


Do you mean [1] you want to plot positive residuals in color and negative in 
gray scale?
Or [2] to fold + and - residuals by shading all according to abs(resid), and
distinguishing + from - by the dashed box outlines?

In fact, I designed this coding scheme so that mosaic plots in color (with my 
blue - white - red scheme) would approximately do exactly what
you might want under [2], when rendered in B/W, since the fully saturated red 
and blue are close in  darkness in B/W.


And shading_hcl() has been written to do exactly what you want under [2]. While 
it is hard to come up with colors of different hues in HSV or HLS space that 
have the same brightness (aka lightness/luminance) and the same
colorfulness (aka chroma), this is easy in HCL.


Try
mosaic(Titanic, gp=shading_Friendly)
save as a jpg/png and try converting to B/W with an image program and see if 
this is good enough.


mosaic(Titanic, shade = TRUE)

is the same as

mosaic(Titanic, gp = shading_hcl)

which you can then modify to have different line types

mosaic(Titanic, gp = shading_hcl, gp_args = list(lty = 1:2))

If you print that on a grayscale printer you will see the same plot without any 
chroma, i.e.,

mosaic(Titanic, gp = shading_hcl, gp_args = list(lty = 1:2, c = 0))

The shading_hcl() function is introduced in Zeileis et al. (2007, JCGS), see 
?shading_hcl, which provides more detailed references to HCL colors etc.

Best,
Z


Alternatively, write your own, shading_Kubovy, modeled on

shading_Friendly -
function (observed = NULL, residuals = NULL, expected = NULL,
  df = NULL, h = c(2/3, 0), lty = 1:2, interpolate = c(2, 4),
  eps = 0.01, line_col = black, ...)
{
  shading_hsv(observed = NULL, residuals = NULL, expected = NULL,
  df = NULL, h = h, v = 1, lty = lty, interpolate = interpolate,
  eps = eps, line_col = line_col, p.value = NA, ...)
}
environment: namespace:vcd
attr(,class)
[1] grapcon_generator

In the defaults, lty=1:2 is what distinguishes + and - for outline line type

hope this helps,
-Michael

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] use sliding window to count substrings found in large string

2010-07-07 Thread Immanuel
Hello together,


I'm looking for advice on how to do some tests on strings.
What I want to do is the following:

(just an example, real strings/sequence are about 200-400 characters long)
given set of Strings:

String1 abcdefgh
String2 bcdefgop

use a sliding window of size x  to create an vector of all subsequences
of size x
found in the set (order matters! ).

Now create, for every string in the set, an vector containing the counts
on how often
each subsequence was found in this particular string.

 It would be great if someone could give me a vague outline on how to
start and which methods to work.
I did read through the man pages and goggled a lot, but still don't know
how to
approach this.

best regards,
Immanuel

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 - four.nines.cartesian.probability.grid

2010-07-07 Thread Charles Annis, P.E.
My routine (below) works OK but misbehaves if the on-screen plot is made
wider using the mouse.

The problem is caused by using 

par(usr)[1] - 0.07 * (par(usr)[2] - par(usr)[1]) 

to locate two items on the y-axis.  The rest of the labeling is controlled
by the line=0 parameter setting.  Of course resizing changes the absolute
plot width, and I'd rather have things with respect to the line= setting.
(Try it: Cut and paste the code, then resize the screen plot to be wider.)

Is there another way to do what I'm trying to do?

I'm using R version 2.11.1 (2010-05-31), running on an HP 64 bit Windows 7
machine.

Thanks



four.nines.cartesian.probability.grid -
function (X.data, x.title = X, y.title = Theoretical Normal CDF,
X.start, X.end)
{
probs - c(0.0001, 0.0003, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1,
0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.98, 0.99, 0.995, 0.998,
0.999, 0.9997, 0.)
z.vals - qnorm(probs)
y.min - z.vals[1]
y.max - z.vals[length(z.vals)]
npts - length(X.data)
X.mean - mean(X.data)
X.sd - sd(X.data)
X.001 - X.mean + X.sd * qnorm(0.001)
X.999 - X.mean + X.sd * qnorm(0.999)
plot(NA, axes = FALSE, xlim = c(X.start, X.end), ylim = c(qnorm(0.0001),
qnorm(0.)), xlab = , ylab = )
axis(side = 1, labels = TRUE, at = pretty(c(X.001, X.data, X.999), n=8),
line = 0, tick = TRUE, outer = FALSE)
mtext(text = x.title, side = 1, line = 2.3, adj = 0.5)

axis(side = 2, at = z.vals, labels = c(NA, NA,
as.character(probs)[-(1:2)]), 
line = 0, tck = -0.01, las = 1, cex.axis = 0.8)

x.loc - par(usr)[1] - 0.07 * (par(usr)[2] - par(usr)[1]) ## If plot
is resized wider this may be troublesome.
text(x = x.loc, y = z.vals[1], bquote(1* X *10^-4), cex = 0.75, xpd =
TRUE)
text(x = x.loc, y = z.vals[2], bquote(3* X *10^-4), cex = 0.75, xpd =
TRUE)

abline(h = qnorm(c(0.90, 0.95, 0.99, 0.999, 0.)), lty = 5, col =
gray)
abline(h = qnorm(c(0.10, 0.05, 0.01, 0.001, 0.0001)), lty = 5, col =
gray)  
mtext(text = y.title, side = 2, line = 3, at = (y.max + y.min)/2, adj =
0.5)
mtext(www.StatisticalEngineering.com, side = 1, line = 2.9, adj = 1,
cex = 0.7)
box()
}

windows(width = 6, height = 6, pointsize = 12)
par(mar = c(4, 4, 0.5, 1) + 0.1)
X.data - rnorm(n=100, mean = 8000, sd = 300)
four.nines.cartesian.probability.grid(X.data, X.start = 7000, X.end = 9000)
points(x = sort(X.data), y = qnorm1:length(X.data)) -
0.5)/length(X.data


Charles Annis, P.E.

charles.an...@statisticalengineering.com
561-352-9699
http://www.StatisticalEngineering.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] Calculate area under a curve

2010-07-07 Thread Greg Snow
The mean values theorem of integration (which I think typifies the differences 
in thinking between mathematicians and statisticians) says that the integral 
from a to b is equal to the average value of the curve between a and b times 
the distance from a to b (b - a).

I would be interested in how similar the areas computed using the other methods 
suggested compare to just averaging your data and multiplying by the distance.

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


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of suse
 Sent: Thursday, July 01, 2010 8:57 AM
 To: r-help@r-project.org
 Subject: [R] Calculate area under a curve
 
 
 Hi,
 
 I want to know the area under a curve, which is not given as a
 function, but
 as values in a time series. It is  not a smooth curve, but switches
 often
 between positive values and zero (the values describe the moisture in
 the
 soil over a year, one entry is one day). I already tried
 area.between.curves, but got only 0 as result. I guess, it doesn't work
 because of these multiple changes between 0 and positive values (most
 of the
 time, the values are 0 and in the certain case I tested, the positive
 values
 were only on single days; but for other values positive values last
 longer)
 I hope, someone can help me!
 
 THank you!
 --
 View this message in context: http://r.789695.n4.nabble.com/Calculate-
 area-under-a-curve-tp2275283p2275283.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-
 guide.html
 and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] Visualization of coefficients

2010-07-07 Thread David Winsemius


On Jul 7, 2010, at 10:19 AM, Michael Friendly wrote:


Tal Galili wrote:

Hello David,
Thanks to your posting I started looking at the function in the arm  
package.
It appears this function is quite mature, and offers (for example)  
the

ability to easily overlap coefficients from several models.
I updated the post I published on the subject, so at the end of it  
I give an

example of comparing the coef of several models:
http://www.r-statistics.com/2010/07/visualization-of-regression-coefficients-in-r/
Thanks again for the pointer.
Best,
Tal


Achim Zeileis wrote:
Re: more mature. arm's coefplot() is more flexible in certain  
respects, mine is more convenient in others. The overlay  
functionality is something arm's coefplot() is better in and it  
also as some further options (vertical vs. horizontal etc.). My  
coefplot() has the advantage that it does not need any modification  
as long as coef() and vcov() methods are available. Furthermore,  
level can specify the significance level (instead of always using  
one and two standard errors, respectively).

But it shouldn't be too hard to create a superset of all options.



@Tal:
For the example using library(arm) and the Mroz data, you posted the  
wrong image.  And loose the intercept in the example.


@Achim:
It would be worthwhile combining the generality of your version with  
the
overlay capability of the arm version, which is extremely useful for  
model comparison.  However, the arm version uses S4 methods.


To my reading coefplot somehow collects parameters from a list of  
model objects using S4 methods, but then passes these to  
coefplot.default which uses base graphics. There seems to be an  
implicit loop, ... perhaps some sort of S4 magic? ... that accumulates  
like-named coefficients in a .local(...) (function?) object until they  
all get reduced to class numeric, at which point they are passed to   
coefplot.default() which does not appear to be an S4 method.


(I am not particularly knowledgeable about S4 functions and  
dispatching, so any corrections or applifications to this account  
would be welcome.)


--
David.




-Michael


--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University  Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele StreetWeb:   http://www.datavis.ca
Toronto, ONT  M3J 1P3 CANADA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] use sliding window to count substrings found in large string

2010-07-07 Thread Gabor Grothendieck
On Wed, Jul 7, 2010 at 12:25 PM, Immanuel mane.d...@googlemail.com wrote:
 Hello together,


 I'm looking for advice on how to do some tests on strings.
 What I want to do is the following:

 (just an example, real strings/sequence are about 200-400 characters long)
 given set of Strings:

 String1 abcdefgh
 String2 bcdefgop

 use a sliding window of size x  to create an vector of all subsequences
 of size x
 found in the set (order matters! ).

 Now create, for every string in the set, an vector containing the counts
 on how often
 each subsequence was found in this particular string.

  It would be great if someone could give me a vague outline on how to
 start and which methods to work.
 I did read through the man pages and goggled a lot, but still don't know
 how to
 approach this.


Try this:

# generate an input string n long
set.seed(123)
n - 300
lets - paste(sample(letters[1:5], n, replace = TRUE), collapse = )

# get rolling k-length sequences and count
k - 3
table(substring(lets, 1:(n-k+1), k:n))

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


Re: [R] F# vs. R

2010-07-07 Thread Sergey Goriatchev
Hello, Marc

No, I do not want to validate Cox PH. :-)
I do use R daily, though right now I do not use the statistical part that much.

I just generally wonder if any R-user tried F# and his/her opinions.

Regards,
Sergey


On Wed, Jul 7, 2010 at 17:56, Marc Schwartz marc_schwa...@me.com wrote:
 On Jul 7, 2010, at 10:31 AM, Sergey Goriatchev wrote:

 Hello, everyone

 F# is now public. Compiled code should run  faster than R.

 Anyone has opinion on F# vs. R? Just curious

 Best,
 S


 The key time critical parts of R are written in compiled C and FORTRAN.

 Of course, if you want to take the time to code and validate a Cox PH or 
 mixed effects model in F# and then run them against R's coxph() or 
 lme()/lmer() functions to test the timing, feel free...  :-)

 So unless there is a pre-existing library of statistical and related 
 functionality for F#, perhaps you need to reconsider your query.

 Regards,

 Marc Schwartz





-- 
---
Kniven skärpes bara mot stenen.

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


Re: [R] nls + quasi-poisson distribution

2010-07-07 Thread Heather Turner
Dear Suresh,

The gnm package for generalized nonlinear models might be what you want
here. This allows you to specify nonlinear models with family=quasipoisson.

For an introduction to the package see the article in R News:
http://cran.r-project.org/doc/Rnews/Rnews_2007-2.pdf

If your model requires a custom nonlin function and you get stuck on
this, or you have any other queries about the package, I'm happy to be
contacted directly.

Best regards,

Heather


On 06/07/10 07:04, Suresh Krishna wrote:
 
 Hello R-helpers,
 
 I would like to fit a non-linear function to data (Discrete X axis,
 over-dispersed Poisson values on the Y axis).
 
 I found the functions gnlr in the gnlm package from Jim Lindsey: this
 can handle nonlinear regression equations for the parameters of Poisson
 and negative binomial distributions, among others. I also found the
 function nls2 in the software package accompanying the book Statistical
 tools for nonlinear regression by Huet et al: this can handle nonlinear
 regression with Poisson distributed Y-axis values.
 
 I was wondering if there was any other option: specifically, any option
 that handled nonlinear fitting with quasi-Poisson distributions (to
 handle the overdispersion).
 
 This is a very new area for me, and I am still trying to figure out the
 best way to do this, so I would appreciate any and all pointers.
 
 Thanks much, Suresh
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] possible to plot number line in R?

2010-07-07 Thread Greg Snow
Here is some code to get you started:


plot.new()
plot.window( c(0,10), c(-1, 1) )

axis(1, at=0:10, pos=0)

lines( c(2,2,5,5), c( -0.25, -0.5, -0.5, -0.25 ) )
text( 3, -0.6, 'Interval 1' )  # the plotrix package has a function for text in 
a box

lines( c(3,3,6,6), c( 0.1, 0.3, 0.3, 0.1) )
text( 4.5, 0.4, 'Interval 2' )

lines( c(4,4,7,7), c( -0.4, -0.6, -0.6, -0.4 ) )
text( 5.5, -0.7, 'Interval 3' )


pieces could be placed in functions to automate the parts that you want.


hope this helps,

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


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of Kroepfl, Julia (julia.kroe...@uni-graz.at)
 Sent: Friday, July 02, 2010 2:40 AM
 To: 'r-help@r-project.org'
 Subject: [R] possible to plot number line in R?
 
 Thank you very much for your answers, but I think I did not explain
 thoroughly enough what I needed. I attached a demo of the plot. I need
 the number line between 2 and 3, both values being shown on the line,
 interval values should be printed next to the dashes and lines should
 connect the intervals.
 
 Is this possible to do in R?
 
 
 
 Best Regards,
 
 Julia
 

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


Re: [R] use sliding window to count substrings found in large string

2010-07-07 Thread Immanuel
Hey,

big help, thanks!
One little question remains, if I create
more then one string and table ...
-

# generate an input string n long
set.seed(123)
n - 300
lets_1 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )
lets_2 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )


# get rolling k-length sequences and count
k - 3
table_1 -table(substring(lets_1, 1:(n-k+1), k:n))
table_2 -table(substring(lets_2, 1:(n-k+1), k:n))
---

is it possible to manipulate table_1 so that it contains zero entries
for all the substrings found in table_2 but not in table_1?

best regards
Immanuel






 Try this:

 # generate an input string n long
 set.seed(123)
 n - 300
 lets - paste(sample(letters[1:5], n, replace = TRUE), collapse = )

 # get rolling k-length sequences and count
 k - 3
 table(substring(lets, 1:(n-k+1), k:n))

   


[[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 on bar chart

2010-07-07 Thread Greg Snow
 fortune(197)

If anything, there should be a Law: Thou Shalt Not Even Think Of Producing A
Graph That Looks Like Anything From A Spreadsheet.
   -- Ted Harding (in a discussion about producing graphics)
  R-help (August 2007)

Filling graphics objects with lines dates back to the days when the only way to 
get decent graphics was with a pen plotter (device that actually had a 
mechanical arm that would move a pen over paper (or paper under the pen)) which 
could not easily fill areas (well you could draw the lines really densely, but 
the usual outcome was that you annoyed other users by taking a long time and 
using all the ink and often ended up with a hole in the paper instead of a nice 
fill.

Technology has progressed quite a bit since then, it is best if your graphs 
reflect that.  Tufte also discusses that such lines actually make graphs harder 
to read, can actually distort features, and cause other problems.

The fact that lattice graphics do not include a simple way to do this is a 
feature.

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


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of ppcrystal
 Sent: Friday, July 02, 2010 11:47 PM
 To: r-help@r-project.org
 Subject: [R] help on bar chart
 
 
 Hey guys,
 
 This is the bar chart that I am working on:
 
 library(lattice);
 data - data.frame(
   X1 = c(2300, 1300, 1300, 450),
   X2 = c(2110, 2220, 1100, 660),
   Y = factor(c(sample1, sample2, sample3, sample4))
   );
 barchart(
   Y ~ X1 + X2,
   data,
   stack = TRUE,
   horiz = TRUE,
   lwd = 1.5,
   xlab = expression(bold(Sample size)),
   col = colors()[c(24,1)],
   xlim = c(0,5000),
   xat = seq(0,5000,1000)
   );
 
 I wanted to make a bar chart that has hatching lines inside the bar:
 with
 sample 2 and 4 having vertical lines and sample 1 and 3 having
 horizontal
 lines, like the following (I kind of photoshopped the image to
 demonstrate
 what I wanted it to look like):
 
 http://r.789695.n4.nabble.com/file/n2277107/test.png
 
 Anyone knows how I can add hatching to the bar charts?
 
 Thanks very much for your time!!!
 --
 View this message in context: http://r.789695.n4.nabble.com/help-on-
 bar-chart-tp2277107p2277107.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-
 guide.html
 and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] use sliding window to count substrings found in large string

2010-07-07 Thread Gabor Grothendieck
On Wed, Jul 7, 2010 at 1:15 PM, Immanuel mane.d...@googlemail.com wrote:
 Hey,

 big help, thanks!
 One little question remains, if I create
 more then one string and table ...
 -

 # generate an input string n long
 set.seed(123)
 n - 300
 lets_1 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )
 lets_2 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )


 # get rolling k-length sequences and count
 k - 3
 table_1 -table(substring(lets_1, 1:(n-k+1), k:n))
 table_2 -table(substring(lets_2, 1:(n-k+1), k:n))
 ---

 is it possible to manipulate table_1 so that it contains zero entries
 for all the substrings found in table_2 but not in table_1?

 best regards
 Immanuel


Turn them into factors with the appropriate levels before counting
them with table:

# generate an input string n long
set.seed(123)
n - 300
lets_1 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )
lets_2 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )

# get rolling k-length sequences and count
k - 3
s1 - substring(lets_1, 1:(n-k+1), k:n)
s2 - substring(lets_2, 1:(n-k+1), k:n)
levs - sort(unique(union(s1, s2)))
table(factors(s1, levs))
table(factors(s2, levs))

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


Re: [R] use sliding window to count substrings found in large string

2010-07-07 Thread Gabor Grothendieck
On Wed, Jul 7, 2010 at 1:25 PM, Gabor Grothendieck
ggrothendi...@gmail.com wrote:
 On Wed, Jul 7, 2010 at 1:15 PM, Immanuel mane.d...@googlemail.com wrote:
 Hey,

 big help, thanks!
 One little question remains, if I create
 more then one string and table ...
 -

 # generate an input string n long
 set.seed(123)
 n - 300
 lets_1 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )
 lets_2 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )


 # get rolling k-length sequences and count
 k - 3
 table_1 -table(substring(lets_1, 1:(n-k+1), k:n))
 table_2 -table(substring(lets_2, 1:(n-k+1), k:n))
 ---

 is it possible to manipulate table_1 so that it contains zero entries
 for all the substrings found in table_2 but not in table_1?

 best regards
 Immanuel


 Turn them into factors with the appropriate levels before counting
 them with table:

 # generate an input string n long
 set.seed(123)
 n - 300
 lets_1 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )
 lets_2 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )

 # get rolling k-length sequences and count
 k - 3
 s1 - substring(lets_1, 1:(n-k+1), k:n)
 s2 - substring(lets_2, 1:(n-k+1), k:n)
 levs - sort(unique(union(s1, s2)))
 table(factors(s1, levs))
 table(factors(s2, levs))


That should be factor, not factors:

table(factor(s1, levs))
table(factor(s2, levs))

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


Re: [R] using objects from different workspaces

2010-07-07 Thread Greg Snow
Check out the high performance computing task view on CRAN.

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


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of Martin Dózsa
 Sent: Sunday, July 04, 2010 5:12 PM
 To: r-help@r-project.org
 Subject: [R] using objects from different workspaces
 
 Hi all,
 
 I have the following problem: I need to run a large number of
 simulations,
 therefore I use many computers. After the computations I need to make
 some
 operations with the obtained results (such as weighted average, sum,
 etc.).
 
 My question is, how is it possible to combine the output of several R
 sessions?
 
 My objects are quite complex (multi-dimensional arrays), therefore
 export to
 some simple csv file seems not to be a solution. I also need to do this
 operation several times, so it would be good if it could be
 automatized.
 
 Thank you in advance,
 Martin
 
   [[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] F# vs. R

2010-07-07 Thread Bernardo Rangel Tura
On Wed, 2010-07-07 at 17:31 +0200, Sergey Goriatchev wrote:
 Hello, everyone
 
 F# is now public. Compiled code should run faster than R.
 
 Anyone has opinion on F# vs. R? Just curious
 
 Best,
 S


Sergey,

F# is public, but is not open source. 

F# run in windows but run in AIX, linux, MAC, UNIX etc?

Compiled code should run faster than R, but is precise?

Compiled code should run faster than R, but is reliable?

Compiled code should run faster than R, but have 2.440 packages for
extend your capacities?

Compiled code should run faster than R, but critical code is in C o
FORTRAN

So I think the F# is not a good alternative, if your concern is velocity
dou you look Littler

http://code.google.com/p/littler/

-- 
Bernardo Rangel Tura, M.D,MPH,Ph.D
National Institute of Cardiology
Brazil

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


Re: [R] Data Labels in a barchart (Lattice or otherwise)

2010-07-07 Thread Greg Snow
 fortune(197)

If anything, there should be a Law: Thou Shalt Not Even Think Of Producing A
Graph That Looks Like Anything From A Spreadsheet.
   -- Ted Harding (in a discussion about producing graphics)
  R-help (August 2007)

Also read the discussion started with:
http://tolstoy.newcastle.edu.au/R/e2/help/07/08/22858.html



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


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of RaoulD
 Sent: Sunday, July 04, 2010 9:44 PM
 To: r-help@r-project.org
 Subject: [R] Data Labels in a barchart (Lattice or otherwise)
 
 
 Hi,
 
 Can anyone please help me with how I could add labels with the value
 for
 each bar in a barchart? (similar to how data labels can be added in
 Excel) I
 have done a lot of searching but havent been lucky.
 
 Thanks,
 Raoul
 --
 View this message in context: http://r.789695.n4.nabble.com/Data-
 Labels-in-a-barchart-Lattice-or-otherwise-tp2278027p2278027.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-
 guide.html
 and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] use sliding window to count substrings found in large string

2010-07-07 Thread Immanuel
Hey,
saved my day.
Now can watch the football semi-final
thanks
 Turn them into factors with the appropriate levels before counting
 them with table:

 # generate an input string n long
 set.seed(123)
 n - 300
 lets_1 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )
 lets_2 - paste(sample(letters[1:5], n, replace = TRUE), collapse = )

 # get rolling k-length sequences and count
 k - 3
 s1 - substring(lets_1, 1:(n-k+1), k:n)
 s2 - substring(lets_2, 1:(n-k+1), k:n)
 levs - sort(unique(union(s1, s2)))
 table(factor(s1, levs))
 table(factor(s2, levs))



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] subscripts for panel.superpose in lattice

2010-07-07 Thread László Sándor
Hi,

I am trying to superimpose (overlay) regression lines to scatter plots
by groups with xyplot (dysfunctional code below). However, my call of
panel.superpose breaks down because of the subscripts requirement. I
tried to research the documentation and examples, but I cannot figure
out how to make xyplot plug subscripts to a panel... call. Could you
have a look? It would be greatly appreciated.

Thank you,

Laszlo

scatter_contrast - function(depvar,bins,cutvar,cutvarname = NULL,
yvarlab = NULL,xvarlab =
NULL,nbins=20,maxbins=100,yrange=c(0,9),plottitle=NULL,legendtitle=NULL)
 {
library('lattice')
library('grid')
trellis.par.set(
plot.symbol = list(cex = 1.5,col=rgb(26,71,111,max=255)),
superpose.symbol = list(cex = rep(1,
times=7),pch=c(15:21),col=c(rgb(26,71,111,max=255),
rgb(144,53,59,max=255),rgb(85,117,47,max=255),#ff,orange,#00ff00,brown),fill=c(rgb(26,71,111,max=255),
  
rgb(144,53,59,max=255),rgb(85,117,47,max=255),#ff,orange,#00ff00,brown)),
plot.line = list(cex = 1,lwd=2,col=rgb(26,71,111,max=255)),
superpose.line = list(cex = rep(1,
times=7),lwd=rep(2,times=7),col=c(rgb(26,71,111,max=255),
rgb(144,53,59,max=255),rgb(85,117,47,max=255),#ff,orange,#00ff00,brown)),
reference.line = list(col=rgb(234,242,243,max=255)),
add.line = list(rgb(85,117,47,max=255),lwd=2),
grid.pars = list(col=black),#rgb(234,242,243,max=255)),
superpose.polygon = 
list(col=black),#rgb(234,242,243,max=255)),
fontsize = list(text=16),   
par.xlab.text = list(cex = 0.8),
par.ylab.text = list(cex = 0.8),
)
if (length(unique(bins))maxbins) bins - binning(bins,nbins)
temp - summary(cutvar)
cut - 1*(cutvar  temp[3])
if (length(unique(!is.na(cutvar)))6) cut - cutvar
legval1 - names(data.frame(cutvar))
xl - names(data.frame(bins))
leg - paste(unique(sort(cut)))
legval2 - leg[1:(length(unique(leg))-1)]
if (length(na.omit(cutvar)) == length(cutvar)) legval2 - leg
ht - depvar[[1]]
if (is.na(ht)) ht - 0
if (ht == hist) {
bins - 
replace(bins,binsquantile(bins,0.99),quantile(bins,0.99))
bins - 
replace(bins,binsquantile(bins,0.05),quantile(bins,0.05))
mes - 
aggregate(matrix(1,length(bins),1),list(bins,cut),sum,na.rm = TRUE)
cnt - aggregate(matrix(1,length(bins),1),list(cut),sum,na.rm = 
TRUE)
mes$x - mes$V1
mes$V1 - NULL
mes - merge(mes,cnt,by.x = Group.2,by.y = Group.1)
mes$x - mes$x/mes$V1 }
else {mes - aggregate(depvar,list(bins,cut),mean,na.rm = TRUE)}
if (yrange[2] == 9) {
ran - (max(mes$x)- min(mes$x))*0.1
yrange = c(min(mes$x)-ran,max(mes$x)+ran)}
xyplot(x ~ Group.1,groups = Group.2,data = mes, subscripts =TRUE,type = 
p,
auto.key = list(cex=0.7,cex.title=0.7,title=legendtitle,space =
bottom,points = FALSE,lines=TRUE,columns=2),
#   key = simpleKey(paste(cutvarname,leg),points = FALSE,lines=TRUE 
),
xlab = xvarlab, ylab = yvarlab,ylim=yrange,
aspect=fill,
panel = function(...) {
panel.grid(h = -1, v = 0)
panel.xyplot(...)
#   panel.abline(lm(depvar ~ bins))
panel.superpose(bins,depvar,subscripts,groups = 
cutvar,type=nr)
panel.axis(side = left, outside = TRUE,tck = 
-1, line.col = 1)
panel.axis(side = bottom, outside = TRUE,tck 
= -1, line.col = 1)},
main = plottitle,
par.settings = list(axis.line = list(col = 
transparent)),
axis = function(side, ...) {
if (side == left) {
grid.lines(x = c(0, 0), y = c(0, 
1),default.units =
npc,gp=gpar(col='black')) }

else if (side == bottom) {
grid.lines(x = c(0, 1), y = c(0, 
0),default.units =
npc,gp=gpar(col='black')) }
axis.default(side = side, ...)
}
)
}

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


Re: [R] export VTK from R : impossible to write data as float

2010-07-07 Thread yepri

So, the problem was that R exports only double sized floats (double), and
Paraview requires single sized floats. The solution was just to write : 

writeBin(data,bfile_celldata,endian=swap,size=4) 

Special Thanks to Sebastian Gibb :
http://www.mail-archive.com/r-help@r-project.org/msg100994.html

alex
-- 
View this message in context: 
http://r.789695.n4.nabble.com/export-VTK-from-R-impossible-to-write-data-as-float-tp2278756p2281217.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] [R-pkgs] ggplot2 version 0.8.8

2010-07-07 Thread Hadley Wickham
ggplot2 

ggplot2 is a plotting system for R, based on the grammar of graphics,
which tries to take the good parts of base and lattice graphics and
avoid bad parts. It takes care of many of the fiddly details
that make plotting a hassle (like drawing legends) as well as
providing a powerful model of graphics that makes it easy to produce
complex multi-layered graphics.

To install or update, run:
install.packages(c(ggplot2, plyr))

Find out more at http://had.co.nz/ggplot2, and check out the nearly 500
examples of ggplot in use.  If you're interested, you can also sign up to
the ggplot2 mailing list at http://groups.google.com/group/ggplot2, or track
development at  http://github.com/hadley/ggplot2

ggplot2 0.8.8 (2010-07-02) 

This version fixes the following 23 bugs:

* coord_equal finally works as expected (thanks to continued prompting
from Jean-Olivier Irisson)
* coord_equal renamed to coord_fixed to better represent capabilities
* coord_polar and coord_polar: new munching system that uses distances
(as defined by the coordinate system) to figure out how many pieces
each segment should be broken in to (thanks to prompting from
Jean-Olivier Irisson)
* fix ordering bug in facet_wrap (thanks to bug report by Frank Davenport)
* geom_errorh correctly responds to height parameter outside of aes
* geom_hline and geom_vline will not impact legend when used for fixed
intercepts
* geom_hline/geom_vline: intercept values not set quite correctly
which caused a problem in conjunction with transformed scales
(reported by Seth Finnegan)
* geom_line: can now stack lines again with position = stack (fixes #74)
* geom_segment: arrows now preserved in non-Cartesian coordinate
system (fixes #117)
* geom_smooth now deals with missing values in the same way as
geom_line (thanks to patch from Karsten Loesing)
* guides: check all axis labels for expressions (reported by Benji Oswald)
* guides: extra 0.5 line margin around legend (fixes #71)
* guides: non-left legend positions now work once more (thanks to
patch from Karsten Loesing)
* label_bquote works with more expressions (factors now cast to
characters, thanks to Baptiste Auguie for bug report)
* scale_color: add missing US spellings
* stat: panels with no non-missing values trigged errors with some
statistics. (reported by Giovanni Dall'Olio)
* stat: statistics now also respect layer parameter inherit.aes
(thanks to bug report by Lorenzo Isella and investigation by Brian
Diggs)
* stat_bin no longer drops 0-count bins by default
* stat_bin: fix small bug when dealing with single bin with NA
position (reported by John Rauser)
* stat_binhex: uses range of data from scales when computing binwidth
so hexes are the same size in all facets (thanks to Nicholas Lewin-Koh
for the bug report)
* stat_qq has new dparam parameter for specifying distribution
parameters (thanks to Yunfeng Zhang for the bug report)
* stat_smooth now uses built-in confidence interval (with small sample
correction) for linear models (thanks to suggestion by Ian Fellows)
* stat_spoke: correctly calculate stat_spoke (cos and sin were
flipped, thanks to Jean-Olivier Irisson for bug report and fix)


-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
http://had.co.nz/

___
R-packages mailing list
r-packa...@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-packages

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


Re: [R] how to draw the legend about color from 3d picture

2010-07-07 Thread Sebastien Guyader

It may not help the original poster, but here's a solution based on what Greg
said above:

# Load plotrix
library(plotrix)

# Create a new layout to divide the graphics in 2, the first one (displaying
the persp() graph) being 4 times larger than the second one (displying the
legend)

layout(matrix(c(1,2),1,2,byrow=T),widths=c(4,1))

# Setup the color palette, here for exemple a 20-step heatmap from white
through yellow to red colors
grey.colors - colorRampPalette( c(white, yellow, black) ) 
color - grey.colors(20)

# Run the persp() command with right definition of facets and colors
zfacet - z[-1,-1] + z[-1,-ncol(z)] + z[-nrow(z),-1] + z[-nrow(z),-ncol(z)]
persp(x, y, z, col=color[cut(zfacet, nbcol)], theta=...)

# Stetup and display the legend, 
col.labels - c(0.0,0.2,0.4,0.6,0.8,1.0) # For z values between
0 and 1
plot(0:10,type=n,axes=FALSE,xlab=,ylab=) # Blank plot required for the
legend to be added to
color.legend(1,1.5,9,8,col.labels,color,gradient=y,align=rb)

-- 
View this message in context: 
http://r.789695.n4.nabble.com/how-to-draw-the-legend-about-color-from-3d-picture-tp866475p2281321.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] plotmath vector problem; full program enclosed

2010-07-07 Thread Paul Johnson
On Wed, Jul 7, 2010 at 5:41 AM, Duncan Murdoch murdoch.dun...@gmail.com wrote:



 You want as.expression(b1), not expression(b1).  The latter means
 the
 expression consisting of the symbol b1.  The former means take the
 object
 stored in b1, and convert it to an expression..


Thanks to Duncan and Allen, who pointed out that I was not even
reading my own code carefully.  I apologize for trying your patience.

Before I stop swinging at this one, I still want to bother everybody
about one thing, which really was the original question, but I did not
know the words to ask it.

The full code below is a working example that works, but I don't
understand why. Focus on these two commands that produce 2 axes.  Both
produce the desired output, but, as far as I can see, they should not!

1:
   axis(1, line=6, at=mu+dividers*sigma,
labels=as.expression(c(b1,b2,b3,b4,b5), padj=-1))


2:
   axis(1, line=9, at=mu+dividers*sigma,
labels=c(as.expression(b1),b2,b3,b4,b5), padj=-1)

This second one shouldn't work, I think.
It has as.expression on only the first element, and yet they all come
out right. Is there a spill over effect?

My original question should not have asked why b1 does not print
correctly when I do this:

   axis(1, line=9, at=mu+dividers*sigma,
labels=c(expression(b1),b2,b3,b4,b5), padj=-1)

but the correct question should have been why do b2, b3, b4 , and b5
get processed properly into plot math even though they are not
expressions??

???
pj

### Filename: plotMathProblem.R
### Paul Johnson July 7, 2010
### email me paulj...@ku.edu

  sigma - 10.0
  mu - 4.0

  myx - seq( mu - 3.5*sigma,  mu+ 3.5*sigma, length.out=500)

  myDensity - dnorm(myx,mean=mu,sd=sigma)

  ### xpd needed to allow writing outside strict box of graph
  ### Need big bottom margin to add several x axes
  par(xpd=TRUE, ps=10, mar=c(18,2,2,2))

  plot(myx, myDensity, type=l, xlab=, ylab=Probability Density ,
main=myTitle1, axes=FALSE)
  axis(2, pos= mu - 3.6*sigma)
  axis(1, pos=0)

  lines(c(myx[1],myx[length(myx)]),c(0,0)) ### closes off axes


  addInteriorLine - function(x, m, sd){
for (i in 1:(length(x))){
  lines( c(x[i],x[i]), c(0, dnorm(x[i],m=m,sd=sd)), lty= 14, lwd=.2)
}
  }


  dividers - c(qnorm(0.025), -1, 0, 1, qnorm(0.975))
  addInteriorLine(mu+sigma*dividers, mu,sigma)

  # bquote creates an expression that text plotters can use
  t1 -  bquote( mu== .(mu))
  mtext(bquote( mu == .(mu)), 1, at=mu, line=-1)


  addInteriorLabel - function(pos1, pos2, m, s){
area - abs(100*( pnorm(m+pos1*s,m,s)-pnorm(m+pos2*s, m,s)))
mid - m+0.5*(pos1+pos2)*s
text(mid, 0.5*dnorm(mid,m,s),label=paste(round(area,2),%))
  }


  addInteriorLabel(dividers[1],dividers[2],  mu, sigma)
  addInteriorLabel(dividers[2],dividers[3],  mu, sigma)
  addInteriorLabel(dividers[3],dividers[4],  mu, sigma)
  addInteriorLabel(dividers[4],dividers[5],  mu, sigma)

   b1 - substitute( mu - d*sigma, list(d=round(dividers[1],2)) )
   b2 - substitute( mu - sigma )
   b3 - substitute( mu )
   b4 - substitute( mu + sigma )
   b5 - substitute( mu + d*sigma, list(d=round(dividers[5],2)) )
   axis(1, line=4, at=mu+dividers*sigma, labels=c(b1,b2,b3,b4,b5), padj=-1)


   axis(1, line=6, at=mu+dividers*sigma,
labels=as.expression(c(b1,b2,b3,b4,b5), padj=-1))



   axis(1, line=9, at=mu+dividers*sigma,
labels=c(as.expression(b1),b2,b3,b4,b5), padj=-1)


-- 
Paul E. Johnson
Professor, Political Science
1541 Lilac Lane, Room 504
University of Kansas

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

2010-07-07 Thread Sarah Lewis
I am pleased to announce to agenda for next weeks LondonR meeting:

LondonR meeting - 13th July 2010

Date:    Tuesday 13th July 2010
Time:    6pm - 9pm
Venue:  The Shooting Star
    125 - 129 Middlesex Street
    E1 7JF
    (Nearest Tubes- Liverpool Street, Moorgate or Bank) 

Agenda

6.00pm Mango Solutions   Introduction

6.20pm Chris Campbell Image analysis using R

7.00pm Matthew Dowle News from data.table 1.4.1 and 1.5

7.40pm     Andy Nicholls   How I'm selling R at GSK

8.30pm Drinks and Networking 

Please register at lond...@mango-solutions.com 

LondonR is a free event.

For more information and past presentations, please visit www.londonr.org
Sarah Lewis

Hadley Wickham, Creator of ggplot2 - first time teaching in the UK. 1st - 2nd  
November 2010. 
To book your seat please go to http://mango-solutions.com/news.html 

T: +44 (0)1249 767700 Ext: 200
F: +44 (0)1249 767707
M: +44 (0)7746 224226
www.mango-solutions.com
Unit 2 Greenways Business Park 
Bellinger Close
Chippenham
Wilts
SN15 1BN
UK 


LEGAL NOTICE
This message is intended for the use o...{{dropped:9}}

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


[R] Weired problem when passing arguments using ...?

2010-07-07 Thread thmsfuller...@gmail.com
Hello All,

I'm trying to pass the argument col.names to write.csv using '...'.
But I got the following warnings. Maybe it is very simple. But I'm not
sure what I am wrong. Could you please help point to me what the
problem is?


#
fun=function(x, ...) {
  fr=parent.frame()
  tmp=get(x, envir=fr)
  write.csv(
  tmp
  , file=paste(x, '.csv', sep='')
  , ...
  )
}

f=data.frame(x=1:10,y=letters[1:10])

fun('f', col.names=F)


 fun('f', col.names=F)
Warning message:
In write.csv(tmp, file = paste(x, .csv, sep = ), ...) :
  attempt to set 'col.names' ignored

-- 
Tom

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


[R] how to process this in R

2010-07-07 Thread jd6688

Here are what i am going to accomplish: 

I have 400 files named as xxx.txt. the content of the file looks like the
following: 

namecount 

1. aaa 100 
2. bbb2000 
3. ccc300 
4. ddd   3000 

 
more that 1000 rows in each files. 

these are the areas i need help: 
1. how can i only read in the files with the string patterns ggg or fff as
part of the file names? 
  for instance, I only need the file names with the ggg or fff in it 
 x_ggg_y_1.txt 
 _fff__xxx.txt 

i don't need to read in the files, such as _aaa_.txt 

2.how cam rename the files: 

  for instance: x_ggg_y_1.txt==changed to ggg1a.txt 


3.after the files read in, how can i only keep the rows with the aaa and
bbb, everything elses show be removed from the files, but the files still
remain the same file name? 

   for instance, in the x_ggg_y_1.txt file, it shouls looks like: 
 namecount 

1. aaa100 
2. bbb2000 
3. aaa300 
4. bbb400 


Thanks so lot, I am very new to R, I am looking forward to any helps from
you. 

-- 
View this message in context: 
http://r.789695.n4.nabble.com/how-to-process-this-in-R-tp2281283p2281283.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] Batch files process and String parsing

2010-07-07 Thread jd6688

Here are what i am going to accomplish:

I have 400 files named as xxx.txt. the content of the file looks like the
following:

namecount

1. aaa 100
2. bbb2000
3. ccc300
4. ddd   3000


more that 1000 rows in each files.

these are the areas i need help:
1. how can i only read in the files with the string patterns ggg or fff as
part of the file names?
  for instance, I only need the file names with the ggg or fff in it
 x_ggg_y_1.txt
 _fff__xxx.txt

i don't need to read in the files, such as _aaa_.txt

2.how cam rename the files:

  for instance: x_ggg_y_1.txt==changed to ggg1a.txt


3.after the files read in, how can i only keep the rows with the aaa and
bbb, everything elses show be removed from the files, but the files still
remain the same file name?

   for instance, in the x_ggg_y_1.txt file, it shouls looks like:
 namecount

1. aaa100
2. bbb2000
3. aaa300
4. bbb400


Thanks so lot, I am very new to R, I am looking forward to any helps from
you.
   
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Batch-files-process-and-String-parsing-tp2281208p2281208.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] forcing a zero level in contr.sum

2010-07-07 Thread Bond, Stephen
I need to use contr.sum and observe that some levels are not statistically 
different from the overall mean of zero.
What is the proper way of forcing the zero estimate? It seems the column 
corresponding to that level should become a column of zeros.
Is there a way to achieve that without me constructing the design matrix?
Thank you.

Stephen Bond


[[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] xyplot of function only shows diagonal plots

2010-07-07 Thread DrHatch
I expected this script to show nine panels, each with a plot of the 
function.

But when I run it, only the diagonal panels have content.
Executing  print(trel_1) shows what I expected in the upper left corner.
I am using R ver. 2.11.0 and lattice ver. 0.18-8 on Windows XP under 
Eclipse and StatET.

I have searched the documentation, but found no answer. What am I missing?

--- simple.R ---

# Test xyplot layout for a function
#
require(lattice)

   f - function(x,a,b){return( a * sin(x)+ b)}
   trellis.par.set(theme = canonical.theme(Windows))
   x - seq(-4,4, by=0.1)
# case 1x1
   a - c( 1 )
   b - c( 3 )
   trel_1 - xyplot(f(x,a,b) ~ x | a * b, ,  type = l,  
   main = 1x1, ylab = f(x), xlab = x,

   xlim = c(-6,6), ylim = c(-6,6),   )
# case 3x3
   a - c(1,2,3)
   b - c(1,2,3)
   trel_3 - xyplot(f(x,a,b) ~ x | a * b, ,  type = l,  
   main = 3x3, ylab = f(x), xlab = x,

   xlim = c(-6,6), ylim = c(-6,6),   )

print(trel_3)
#print(trel_1)

--- end of simple.R ---

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


  1   2   >