Re: [R] create a new data frame after comparing two columns of the previous data frame

2011-06-28 Thread Petr PIKAL
Hi

r-help-boun...@r-project.org napsal dne 28.06.2011 02:45:01:

 
 Right. You forgot the comma!

If you want to resolve this rather cryptic advice see

?[

Regards
Petr

 
 -- 
 David.
 
 On Jun 27, 2011, at 7:10 PM, Nanami wrote:
 
  it's a table; I read it from a file; I've tried to make it look 
  prettier:
 
  head(intra)
  chr   miRNA   start end   strand  ACC
  hsa_ID   region region_start region_end
  chr1 miRNA 1102484 1102578  +   ACC=MI342; ID=hsa- 
  mir-200b;
  exon  11024841102578
  chr1 miRNA 1103243 1103332  +   ACC=MI737; ID=hsa- 
  mir-200a;
  exon  11032431103332
  chr1 miRNA 1104385 1104467  +   ACC=MI0001641;  ID=hsa- 
  mir-429;
  exon  11043851104467
  chr1 miRNA 3044539 3044599  +   ACC=MI0015861; ID=hsa- 
  mir-4251;
  exon  30445393044599
  chr1 miRNA 3477260 3477354  -   ACC=MI0003556; ID=hsa- 
  mir-551a;
  exon  34772603477354
  chr1 miRNA 6489894 6489956  -   ACC=MI0015864; ID=hsa- 
  mir-4252;
  exon  64898946489956
 
 
  When I changed  with   I got this error:
 
  new - intra[(intra$start != intra$region_start )(intra$end !=
  intra$region_end)]
  Error in `[.data.frame`(intra, (intra$start != intra$region_start) 
  (intra$end !=  :
   undefined columns selected
 
  What I would like to do is to create a new table that only contains 
  the rows
  for which column start and region_start are not the same.
 
 
 
  --
  View this message in context: http://r.789695.n4.nabble.com/create-a-
 new-data-frame-after-comparing-two-columns-of-the-previous-data-frame-
 tp3628981p3629092.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.
 
 David Winsemius, MD
 West Hartford, CT
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

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


[R] data.frame: How to get the classes of all components and how to remove their factor structure?

2011-06-28 Thread Marius Hofert
Dear expeRts,

I have two questions concerning data frames:
(1) How can I apply the class function to each component in a data.frame? As 
you can see below, applying class to each column is not the right approach; 
applying it to each component seems bulky.
(2) After transforming the data frame a bit, the classes of certain components 
change to factor. How can I remove the factor structure?

Cheers,

Marius

x - c(2004:2010, 2002:2011, 2000:2011)
df - data.frame(x=x, group=c(rep(low,7), rep(middle,10), rep(high,12)), 
 y=x+100*runif(length(x))) 

## Question (1): why do the following lines do not give the same class?
apply(df, 2, class)
class(df$x)
class(df$group)
class(df$y)

df. - as.data.frame(xtabs(y ~ x + group, data=df))

class(df.$x)
class(df.$group)
class(df.$Freq)

## Question (2): how can I remove the factor structure from x?
df.$x - as.numeric(as.character(df.$x)) # seems bulky; note that 
as.numeric(df.$x) is not correct
class(df.$x)
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] means and error bars on xyplot for binary data

2011-06-28 Thread Louis Plough
Hi,
I have binary (0,1) data for a trait as my response variable, and
a dependent variable, genotype, with three classes (AA, AB, BB).

I would like to plot this data so that across the three genoytpes, even
though the points are all either 0 or 1, i want them to stack up or be seen
using 'jitter'.  So far I have been able to do this using xyplot {lattice}
(code below) but could not get the points to jitter or stack up on boxplot
or bwplot {lattice}.

 I would also like to add to the xyplot object, the mean of the values at
each of these classes  which will vary depending on how many data points are
at 0 and 1 for a given genotype.  I would also like to put error (i.e.
standard error) bars around these estimates.

I have tried using the points() function to put the mean at each of the
genotype classes, but the point ends up off the figure. Any ideas how to get
this going?

here is some example code.

 gtype-c(AA,AB,BB)
x-sample(gtype,20,replace=TRUE)
y-sample(c(0,1),20,replace=TRUE)
bin.data-data.frame(x,y)
xyplot(y~x, jitter.y=TRUE, jitter.x=TRUE,factor=.6, data=bin.data)

Then If I wanted to add the means to the plot, I would do this, which will
print the mean points on a box plot, but not an xyplot:

means1 - tapply(bin.data$y,bin.data$x,mean)
points(means1,col=red,pch=18)

Is there a way to get the means, and even error bars on the same xyplot?

Thanks for any help in advance,
Louis

[[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] Odp: data.frame: How to get the classes of all components and how to remove their factor structure?

2011-06-28 Thread Petr PIKAL
Hi

 Dear expeRts,
 
 I have two questions concerning data frames:
 (1) How can I apply the class function to each component in a 
data.frame? 
 As you can see below, applying class to each column is not the right 
 approach; applying it to each component seems bulky.
 (2) After transforming the data frame a bit, the classes of certain 
 components change to factor. How can I remove the factor structure?
 
 Cheers,
 
 Marius
 
 x - c(2004:2010, 2002:2011, 2000:2011)
 df - data.frame(x=x, group=c(rep(low,7), rep(middle,10), 
rep(high,12)), 
  y=x+100*runif(length(x))) 
 
 ## Question (1): why do the following lines do not give the same 
class?

from help page
?apply
Arguments
X
an array, including a matrix.

array is not a data frame

 apply(df, 2, class)
 class(df$x)
 class(df$group)
 class(df$y)

sapply(df, class)
x group y 
integer  factor numeric 


 
 df. - as.data.frame(xtabs(y ~ x + group, data=df))
 
 class(df.$x)
 class(df.$group)
 class(df.$Freq)
 
 ## Question (2): how can I remove the factor structure from x?
 df.$x - as.numeric(as.character(df.$x)) # seems bulky; note that 
 as.numeric(df.$x) is not correct

Actually it is correct in a sense it behaves as documented

?factor

Warning
The interpretation of a factor depends on both the codes and the levels 
attribute. Be careful only to compare factors with the same set of levels 
(in the same order). In particular, as.numeric applied to a factor is 
meaningless, and may happen by implicit coercion. To transform a factor f 
to approximately its original numeric values, as.numeric(levels(f))[f] is 
recommended and slightly more efficient than as.numeric(as.character(f)). 


Regards
Petr

 class(df.$x)





 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] Odp: data.frame: How to get the classes of all components and how to remove their factor structure?

2011-06-28 Thread Petr PIKAL
 
 Dear expeRts,
 
 I have two questions concerning data frames:
 (1) How can I apply the class function to each component in a 
data.frame? 
 As you can see below, applying class to each column is not the right 
 approach; applying it to each component seems bulky.
 (2) After transforming the data frame a bit, the classes of certain 
 components change to factor. How can I remove the factor structure?
 
 Cheers,
 
 Marius
 
 x - c(2004:2010, 2002:2011, 2000:2011)
 df - data.frame(x=x, group=c(rep(low,7), rep(middle,10), 
rep(high,12)), 
  y=x+100*runif(length(x))) 
 
 ## Question (1): why do the following lines do not give the same 
class?
 apply(df, 2, class)
 class(df$x)
 class(df$group)
 class(df$y)
 
 df. - as.data.frame(xtabs(y ~ x + group, data=df))
 
 class(df.$x)
 class(df.$group)
 class(df.$Freq)
 
 ## Question (2): how can I remove the factor structure from x?
 df.$x - as.numeric(as.character(df.$x)) # seems bulky; note that 

If you do it often you can

unfactor - function(x) as.numeric(as.character(x))
df.$x - unfactor(df.$x)

or you can use 
df. - as.data.frame(xtabs(y ~ x + group, data=df), 
stringsAsFactors=FALSE)
df.$x - as.numeric(df.$x)

But it seems to me that it is not much less bulkier.

Regards
Petr


 as.numeric(df.$x) is not correct
 class(df.$x)
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] testInstalledPackages

2011-06-28 Thread Uwe Ligges



On 27.06.2011 23:56, Cody Hamilton wrote:

Dear group,

When running the installation test:

testInstalledPackages(both,outDir='c:/Test')

I got the following message:

Running ‘testci.R’
comparing ‘testci.Rout’ to ‘testci.Rout.save’ ...
files differ in number of lines:


and then?



Please note the test does not result in 'OK' as do the other tests.  Is this a 
concern?


We do not know the package nor the output of the test / diffs. What 
answer do you expect now?


Best,
Uwe Ligges





Regards,
-Cody Hamilton

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

2011-06-28 Thread Uwe Ligges



On 27.06.2011 23:54, xin123620 wrote:

Dear R Users,

I was using R to import several years traffic data, but every time after
data successfully imported (no error or warning) when I tried to save this
workplace or tackle these data, R GUI would automatically shut down. Would
you have any ideas that how I should deal with this situation?  Thank you
very much.


- If this is not R-patched with updated package, update
- Try to make it reproducible and give us the relevant info so that we 
can reproduce.

- Tell what sessionInfo() says directly before you manage to crash R.

Uwe Ligges




Best,
Chengxin

--
View this message in context: 
http://r.789695.n4.nabble.com/R-GUI-shutdown-tp3628965p3628965.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] testInstalledPackages

2011-06-28 Thread Duncan Murdoch

On 27/06/2011 5:56 PM, Cody Hamilton wrote:

Dear group,

When running the installation test:

testInstalledPackages(both,outDir='c:/Test')

I got the following message:

Running ‘testci.R’
comparing ‘testci.Rout’ to ‘testci.Rout.save’ ...
files differ in number of lines:

Please note the test does not result in 'OK' as do the other tests.  Is this a 
concern?


Yes, you should look at the two files to determine why the length has 
changed.


Duncan Murdoch



Regards,
-Cody Hamilton

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Odp: data.frame: How to get the classes of all components and how to remove their factor structure?

2011-06-28 Thread Marius Hofert
Dear Petr,

thanks for your posts, they perfectly answered my questions.

Cheers,

Marius


On 2011-06-28, at 09:49 , Petr PIKAL wrote:

 
 Dear expeRts,
 
 I have two questions concerning data frames:
 (1) How can I apply the class function to each component in a 
 data.frame? 
 As you can see below, applying class to each column is not the right 
 approach; applying it to each component seems bulky.
 (2) After transforming the data frame a bit, the classes of certain 
 components change to factor. How can I remove the factor structure?
 
 Cheers,
 
 Marius
 
 x - c(2004:2010, 2002:2011, 2000:2011)
 df - data.frame(x=x, group=c(rep(low,7), rep(middle,10), 
 rep(high,12)), 
 y=x+100*runif(length(x))) 
 
 ## Question (1): why do the following lines do not give the same 
 class?
 apply(df, 2, class)
 class(df$x)
 class(df$group)
 class(df$y)
 
 df. - as.data.frame(xtabs(y ~ x + group, data=df))
 
 class(df.$x)
 class(df.$group)
 class(df.$Freq)
 
 ## Question (2): how can I remove the factor structure from x?
 df.$x - as.numeric(as.character(df.$x)) # seems bulky; note that 
 
 If you do it often you can
 
 unfactor - function(x) as.numeric(as.character(x))
 df.$x - unfactor(df.$x)
 
 or you can use 
 df. - as.data.frame(xtabs(y ~ x + group, data=df), 
 stringsAsFactors=FALSE)
 df.$x - as.numeric(df.$x)
 
 But it seems to me that it is not much less bulkier.
 
 Regards
 Petr
 
 
 as.numeric(df.$x) is not correct
 class(df.$x)
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] time series question

2011-06-28 Thread Albert-Jan Roskam
Hi,

I have a long data file with time data that change to wide format using 
reshape. 
The data contain Values and Factors. Some values are missing but can be 
obtained 
by multiplying value of year T-1 with Factor of year T. Sometimes, multiple 
succesive years have no values, so the calculated values need to be calculated 
sequentially. 


# sample data.
DF - data.frame(Var=rep(letters, 10), 
Fac=rep(runif(26), 10), 
Val=rep(runif(26), 10), 
Year=rep(2000:2009, 26))
DF[as.numeric(rownames(DF)) %% 3 == 0,Val] - NA # make some holes
DF2 - cast(melt(DF, id=c(Var, Year)), ... ~ variable + Year, fun=mean, 
na.rm=T)

# my attempt
library(reshape)
prev - grep(Val_, names(DF2)) - 1
this - grep(Fac_, names(DF2))
DF3 - DF2
DF3[, prev] - mapply(*, DF2[, this], DF2[, prev])

 
This doesn't work. Another option would be to use two loops for cols and rows, 
but I didn't get that to work either :-(

Suggestions for clean code, anyone?

Thank you in advance!

Cheers!!
Albert-Jan


~~
All right, but apart from the sanitation, the medicine, education, wine, public 
order, irrigation, roads, a fresh water system, and public health, what have 
the 
Romans ever done for us?
~~

[[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] lattice multiple y-scale possible?

2011-06-28 Thread Steve Corsi

Hi
I am attempting to use the lattice bwplot function to generate boxplots 
of numerous parameters (1-panel/parameter) by site (x-axis). The 
parameters have quite different ranges of values, so it would be best to 
have a separate y-axis range for each panel. Below is a basic example of 
what I am trying to do. As is seen, the y-axes need to be scaled 
individually to make this useful. Any information on how to do this 
would be much appreciated:


rm(mydat)
rm(tempdf)

for (i in 1:5){
for (ii in 1:5){
dat - sample(1:20)*10^ii
tempdf - data.frame(dat)
tempdf$parameter - paste(parameter ,ii,sep=)
tempdf$site - paste(site,i,sep=)
if(!exists(mydat)) {mydat - tempdf
  }else {mydat - rbind(mydat,tempdf)}
  }
}

bwplot(dat~site|parameter,data=mydat,
  layout=c(1,5),
  cex=2,
  xlab=Site Name,
  ylab=,
  labels=levels(mydat$site), scales=list(tick.number=list(4)))

thanks
Steve

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


[R] Running R from windows command prompt

2011-06-28 Thread siddharth arun
1. I have a R program in a file say functions.R.
I load the functions.R file the R using source(function.R) and then call
functionsf1(), f2() etc. which are declared and defined within function.R
 file.
I also need to load a couple of R libraries using library() before I can
use f1(), f2() etc.

My question is can I acheive all this (i.e. calling function f1() and f2())
from the windows prompt without opening the R environment ? If yes then how?


2. Also, Is there any way to scan strings directly. Like scan() function
only scans numerical values. Is there any way to scan strings?

-- 
Siddharth Arun,
4th Year Undergraduate student
Industrial Engineering and Management,
IIT Kharagpur

[[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] Running R from windows command prompt

2011-06-28 Thread Uwe Ligges



On 28.06.2011 11:54, siddharth arun wrote:

1. I have a R program in a file say functions.R.
I load the functions.R file the R using source(function.R) and then call
functionsf1(), f2() etc. which are declared and defined within function.R
  file.
I also need to load a couple of R libraries using library() before I can
use f1(), f2() etc.

My question is can I acheive all this (i.e. calling function f1() and f2())
from the windows prompt without opening the R environment ? If yes then how?



Put all the code into a file, e.g. foo.R, and run

R CMD BATCH foo.R from the windows command shell.



2. Also, Is there any way to scan strings directly. Like scan() function
only scans numerical values. Is there any way to scan strings?


Yes, it is called scan() which is for arbitrary data, including 
character(). See ?scan. Otherwise, you may want to look into ?readLines


Uwe Ligges

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

2011-06-28 Thread suman pal
Hello,
I basically want to use R-help, and post some problems which I am facing. The 
Ref is a well known Genome Biology paper Bioconductor: open software 
development for computational biology and bioinformatics by Robert C Gentleman 
et al., 2004. Generating Heatmaps till Fig2 is working so I think esetSel is 
not the problem..
However, for generating the Figure 3, for GO annotations the following command 
is not working:
 ll - mget(geneNames(esetSel),hgu95av2LOCUSID)   
#it is displaying error messages  Error in mget(geneNames(esetSel), 
hgu95av2LOCUSID) : object 'hgu95av2LOCUSID' not foundand also geneNames not 
found try featureNames instead
Hence I cant proceed to the next set of commands provided in the paper which 
are as follows:
 ll - unique(unlist(ll)) mf - as.data.frame(GOHyperG(ll))[, 1:3] mf - 
 mf[mf$pvalue  0.01, ]
sincerelySumanCCMB Hyderabad.

[[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] gam confidence interval (package mgcv)

2011-06-28 Thread Simon Wood
not sure if I'm missing something here, but since you are using a log 
link, isn't the ratio you are looking for given by the `treatmentB' 
parameter in the summary (independent of X)


 summary(gfit)
[snip]

Parametric coefficients:
Estimate Std. Error t value Pr(|t|)
(Intercept)  1.831830.03885  47.152  2e-16 ***
treatmentB   0.445130.05567   7.996  1.4e-09 ***
---
[snip]

Let mu = E(y), and b be a binary indicator for treatment B. You want
mu|b=1/mu|b=0

log(mu|b=1) = intercept + treatmentB + s(X)
log(mu|b=0) = intercept + s(X)

= log(mu|b=1) - log(mu|b=0) = treatmentB

so mu|b=1/mu|b = exp(treatmentB)

So you can get the required interval by finding and interval for 
treatment B and exponentiating...


tB - coef(gfit)[2]
se.tB - sqrt(vcov(gfit)[2,2])
exp(c(tB - 2*se.tB,tB+2*se.tB))

On 06/28/2011 03:45 AM, Remko Duursma wrote:

Dear R-helpers,

I am trying to construct a confidence interval on a prediction of a
gam fit. I have the Wood (2006) book, and section 5.2.7 seems relevant
but I am not able to apply that to this, different, problem.

Any help is appreciated!

Basically I have a function Y = f(X) for two different treatments A
and B.  I am interested in the treatment ratios : Y(treatment = B) /
Y(treatment = A) as a function of X, including a confidence interval
for this treatment ratio (because we are testing this ratio against
some value, across the range of X).

The X values that Y is measured at differs between the treatments, but
the ranges are similar.


# Reproducible problem:
X1- runif(20, 0.5, 4)
X2- runif(20, 0.5, 4)

Y1- 20*exp(-0.5*X1) + rnorm(20)
Y2- 30*exp(-0.5*X2) + rnorm(20)

# Look at data:
plot(X1, Y1, pch=19, col=blue, ylim=c(0,max(Y1,Y2)), xlim=c(0,5))
points(X2, Y2, pch=19, col=red)

# Full dataset
dfr- data.frame(X=c(X1,X2), Y=c(Y1,Y2), treatment=c(rep(A,20),rep(B,20)))

# Fit gam
# I use a gamma family here although it is not necessary: in the real
problem it is, though.
gfit- gam(Y ~ treatment + s(X), data=dfr, family=Gamma(link=log))

# I am interested in the relationship:
# Y(treatment ==B) / Y(treatment==A) as a function of X, with a
confidence interval!

Do I just do a bootstrap here? Or is there a more appropriate method?

Thanks a lot for any help.

greetings,
Remko





-
Remko Duursma
Research Lecturer

Hawkesbury Institute for the Environment
University of Western Sydney
Hawkesbury Campus, Richmond

Mobile: +61 (0)422 096908
www.remkoduursma.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] Running R from windows command prompt

2011-06-28 Thread siddharth arun
Thanks for your help.

I tried the way you mentioned for my first question. But I am not getting
any results.
Can you please explain in detail the process through which I can run a R
code from windows command prompt.

2011/6/28 Uwe Ligges lig...@statistik.tu-dortmund.de



 On 28.06.2011 11:54, siddharth arun wrote:

 1. I have a R program in a file say functions.R.
 I load the functions.R file the R using source(function.R) and then
 call
 functionsf1(), f2() etc. which are declared and defined within
 function.R
  file.
 I also need to load a couple of R libraries using library() before I can
 use f1(), f2() etc.

 My question is can I acheive all this (i.e. calling function f1() and
 f2())
 from the windows prompt without opening the R environment ? If yes then
 how?



 Put all the code into a file, e.g. foo.R, and run

 R CMD BATCH foo.R from the windows command shell.



  2. Also, Is there any way to scan strings directly. Like scan() function
 only scans numerical values. Is there any way to scan strings?


 Yes, it is called scan() which is for arbitrary data, including
 character(). See ?scan. Otherwise, you may want to look into ?readLines

 Uwe Ligges




-- 
Siddharth Arun,
4th Year Undergraduate student
Industrial Engineering and Management,
IIT Kharagpur

[[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] problem with corrgram function

2011-06-28 Thread Niels Janssen

Dear list,

I have a problem with the corrgram function. It does not seem to 
color large negative correlations, while the same correlation, if 
positive, provides no problems. Is this a bug?


require(corrgram)   
a = seq(1,100)
b = -jitter(seq(1,100), 80)
cor(a,b) # r about -.96
c=as.data.frame(cbind(a,b))
	corrgram(c, order=NULL, lower.panel=panel.pie,upper.panel=NULL, 
text.panel=panel.txt) # no color


c$b = -1*c$b # flip direction of correlation
cor(c$a, c$b) # r now about +.96
	corrgram(c, order=NULL, lower.panel=panel.pie,upper.panel=NULL, 
text.panel=panel.txt) #no problem with color.


Thanks!

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


Re: [R] lattice multiple y-scale possible?

2011-06-28 Thread Dennis Murphy
Hi:

bwplot(dat~site|parameter,data=mydat,
 layout=c(1,5),
 cex=2,
 xlab=Site Name,
 ylab=,
 labels=levels(mydat$site),
 scales=list(tick.number=list(4), rot = c(0, 90),
 y = list(relation = 'free')))

Does that work?

Dennis

On Mon, Jun 27, 2011 at 10:46 PM, Steve Corsi srco...@usgs.gov wrote:
 Hi
 I am attempting to use the lattice bwplot function to generate boxplots of
 numerous parameters (1-panel/parameter) by site (x-axis). The parameters
 have quite different ranges of values, so it would be best to have a
 separate y-axis range for each panel. Below is a basic example of what I am
 trying to do. As is seen, the y-axes need to be scaled individually to make
 this useful. Any information on how to do this would be much appreciated:

 rm(mydat)
 rm(tempdf)

 for (i in 1:5){
 for (ii in 1:5){
 dat - sample(1:20)*10^ii
 tempdf - data.frame(dat)
 tempdf$parameter - paste(parameter ,ii,sep=)
 tempdf$site - paste(site,i,sep=)
 if(!exists(mydat)) {mydat - tempdf
  }else {mydat - rbind(mydat,tempdf)}
  }
 }

 bwplot(dat~site|parameter,data=mydat,
      layout=c(1,5),
      cex=2,
      xlab=Site Name,
      ylab=,
      labels=levels(mydat$site), scales=list(tick.number=list(4)))

 thanks
 Steve

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


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

2011-06-28 Thread siddharth arun
I want to find the error in the forecasting values.
I used the function accuracy() under forecast library. But, I didn't
understand how it is calculating the errors?

Can anybody help?

-- 
Siddharth Arun,
4th Year Undergraduate student
Industrial Engineering and Management,
IIT Kharagpur

[[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] Creating a Polar Plot with expanding points as radius increases

2011-06-28 Thread Jim Lemon

On 06/28/2011 08:35 AM, Patrick Jemison wrote:

I'd like to create a polar plot similar to those created by the polarFreq
function in the openair package.  However, this package seems to be specific
to wind speed and direction, and requires a ws (wind speed) and a wd
(wind direction) column.  My data is unrelated to wind speed, but I'd like
to be able to get a plot that does what polarFreq's plots do; I'd like to
have points that expand as their distance from the center increases.

http://www.oga-lab.net/RGM2/func.php?rd_id=openair:polarFreq

Without this property, my plot looks more like spokes on a wheel, or similar
to the polar.plot function in the plotrix package:
http://www.oga-lab.net/RGM2/func.php?rd_id=plotrix:polar.plot

Is there a package or function available that allows the creation of a plot
as described above for a generic set of polar data?  I spent a few hours
today trying to find something I could use, but no luck yet.

I'm using data of the form: x = r, y = theta, and z = a value to plot at the
corresponding r,theta (as a comparison, for the polarFreq plot, z would be
the frequency that is indicated by various colors in the plot).


Hi Patrick,
polar.plot (or other plots in the radial.plot family) currently don't 
plot sectors or fill the (what the heck is the intersection of a sector 
and an annulus called?) to illustrate a third value (say, proportion) in 
addition to direction and intensity. However, it can be done, and I 
might just have a shot at it if I can scratch up an hour or two sometime 
in the next couple of years. Well, maybe sooner than that. Look out for 
radial.pie.


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] Kernel Density Estimation at manually specified points

2011-06-28 Thread Stephen Ellison
May be a kludge, but it might be simpler to write your own density function for 
a few specified points.

For example

my.density - function(x, bw = 'nrd0', at) {
x-na.omit(x)

#
#Borrowed from density.default for compatibility

if (is.character(bw)) {
if (length(x)  2) 
stop(need at least 2 points to select a bandwidth automatically)
bw - switch(tolower(bw), nrd0 = bw.nrd0(x), nrd = bw.nrd(x), 
ucv = bw.ucv(x), bcv = bw.bcv(x), sj = , `sj-ste` = bw.SJ(x, 
method = ste), `sj-dpi` = bw.SJ(x, method = dpi), 
stop(unknown bandwidth rule))
}
##
at - matrix(at, ncol=1)
y - apply(at, 1, FUN=function(a, x, bw) sum(dnorm(a, x, bw)/length(x)), 
x=x, bw=bw )
return(list(x=at, y=y))

}

x-rnorm(500)
plot(density(x))
at-seq(-2, 2, 0.25)
points(my.density(x, at=at))

 

 -Original Message-
 From: r-help-boun...@r-project.org 
 [mailto:r-help-boun...@r-project.org] On Behalf Of Carsten Harlaß
 Sent: 27 June 2011 20:19
 To: dcarl...@tamu.edu
 Cc: r-help@r-project.org
 Subject: Re: [R] Kernel Density Estimation at manually 
 specified points
 
 Hello,
 
 thanks for your response.
 
 Of course I already thought about this simple solution of 
 the problem.
 But I think this is not a nice workaround.
 
 If I understand the manual correctly, density() already makes 
 an approximation. But this approximation is just evaluated at 
 the wrong points.
 
 See:
 
 Details
 
 The algorithm used in density.default disperses the mass of 
 the empirical distribution function over a regular grid of at 
 least 512 points and then uses the fast Fourier transform to 
 convolve this approximation with a discretized version of the 
 kernel and then uses linear approximation to evaluate the 
 density at the specified points.
 
 and
 
 nthe number of equally spaced points at which the 
 density is to be
 estimated. When n  512, it is rounded up to a power of 2 
 during the calculations (as fft is used) and the final result 
 is interpolated by approx. So it almost always makes sense to 
 specify n as a power of two. 
 
 So this workaround means putting a second approximation on 
 top of an other approximation. That's not nice. Or is it?
 
 With kind regards
 
 Carsten
 
 Am 27.06.2011 19:16, schrieb David L Carlson:
  Look at ?approx. For your example (of course your random 
 numbers give 
  different results):
  
  approx(f$x, f$y, c(-2, -1, 0, 1, 2))
  $x
  [1] -2 -1  0  1  2
  
  $y
  [1] 0.03757113 0.19007982 0.31941779 0.37066592 0.10227509
  
  approx gives NA's if you try to interpolate outside the 
 bounds of the data.
  --
  David L Carlson
  Associate Professor of Anthropology
  Texas AM University
  College Station, TX 77843-4352
  
  
  -Original Message-
  From: r-help-boun...@r-project.org 
  [mailto:r-help-boun...@r-project.org] On Behalf Of Carsten Harlaß
  Sent: Sunday, June 26, 2011 7:02 PM
  To: r-help@r-project.org
  Subject: [R] Kernel Density Estimation at manually specified points
  
  Hello,
  
  my name is Carsten. This ist my first post to R-help mailing list.
  
  I estimate densities with the function density out of the package 
  stats.
  
  A simplified example:
  
  
  #generation of test data
  n=10
  z = rnorm(n)
  
  #density estimation
  f=density(z,kernel=epanechnikov,n=n)
  
  #evaluation
  print(f$y[5])
  
  Here I can only evaluate the estimation at given points. 
 These points 
  are determined by the parameter n. By default they are equidistant 
  distributed on the interesting interval.
  
  But I need to evaluate the estimation (the estimated densitiy 
  function) at manually specified points. For example I want 
 to compute f(z[i]).
  This means I am interested in the estimated density at a the 
  observation z[i].
  
  Does anyone know how I can compute this? I think this is an 
 ordinary 
  task so I would be surprised if R can not manage this. But 
 even after 
  a long search I have found nothing.
  
  Thanks in advance
  
  Carsten Harlaß
  
  --
  Carsten Harlaß
  Aachen University of Applied Sciences
  Campus Jülich
  E-Mail: carsten_harl...@web.de
  
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide 
  http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
  
 
  --
  Carsten Harlass
  Aachen University of Applied Sciences
  Campus Juelich
  E-Mail: carsten_harl...@web.de
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

Re: [R] Running R from windows command prompt

2011-06-28 Thread Mike Marchywka




 To: lig...@statistik.tu-dortmund.de
 CC: r-help@r-project.org
 Subject: Re: [R] Running R from windows command prompt
 
 Thanks for your help.
 
 I tried the way you mentioned for my first question. But I am not getting
 any results.
 Can you please explain in detail the process through which I can run a R
 code from windows command prompt.

While your problem likely has nothing to do with windoh's, I would suggest
you go get cygwin ( see google) and use that. I and probably others
have lots of scripts that work there and on linux. Bash scripts are
a better way to proceed if you want to do more than a test case and
they integrate with lots of other existing things. 


 
 2011/6/28 Uwe Ligges lig...@statistik.tu-dortmund.de
 
 
 
  On 28.06.2011 11:54, siddharth arun wrote:
 
  1. I have a R program in a file say functions.R.
  I load the functions.R file the R using source(function.R) and then
  call
  functionsf1(), f2() etc. which are declared and defined within
  function.R
   file.
  I also need to load a couple of R libraries using library() before I can
  use f1(), f2() etc.
 
  My question is can I acheive all this (i.e. calling function f1() and
  f2())
  from the windows prompt without opening the R environment ? If yes then
  how?
 
 
 
  Put all the code into a file, e.g. foo.R, and run
 
  R CMD BATCH foo.R from the windows command shell.
 
 
 
   2. Also, Is there any way to scan strings directly. Like scan() function
  only scans numerical values. Is there any way to scan strings?
 
 
  Yes, it is called scan() which is for arbitrary data, including
  character(). See ?scan. Otherwise, you may want to look into ?readLines
 
  Uwe Ligges
 
 
 
 
 -- 
 Siddharth Arun,
 4th Year Undergraduate student
 Industrial Engineering and Management,
 IIT Kharagpur
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
  
[[alternative HTML version deleted]]

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


Re: [R] Running R from windows command prompt

2011-06-28 Thread Uwe Ligges



On 28.06.2011 13:56, Mike Marchywka wrote:






To: lig...@statistik.tu-dortmund.de
CC: r-help@r-project.org
Subject: Re: [R] Running R from windows command prompt

Thanks for your help.

I tried the way you mentioned for my first question. But I am not getting
any results.
Can you please explain in detail the process through which I can run a R
code from windows command prompt.



Either save your results using a function call in the foo.R file or just 
look into the foo.Rout file R generated for you.




While your problem likely has nothing to do with windoh's, I would suggest
you go get cygwin ( see google) and use that.


Sure you can do, but note that cygwin is not officially supported.



I and probably others
have lots of scripts that work there and on linux. Bash scripts are
a better way to proceed if you want to do more than a test case and
they integrate with lots of other existing things.


R CMD BATCH foo.R behaves almost the same way in both Windows command 
shell and the bash.


Best,
Uwe Ligges








2011/6/28 Uwe Liggeslig...@statistik.tu-dortmund.de




On 28.06.2011 11:54, siddharth arun wrote:


1. I have a R program in a file say functions.R.
I load the functions.R file the R using source(function.R) and then
call
functionsf1(), f2() etc. which are declared and defined within
function.R
  file.
I also need to load a couple of R libraries using library() before I can
use f1(), f2() etc.

My question is can I acheive all this (i.e. calling function f1() and
f2())
from the windows prompt without opening the R environment ? If yes then
how?




Put all the code into a file, e.g. foo.R, and run

R CMD BATCH foo.R from the windows command shell.



  2. Also, Is there any way to scan strings directly. Like scan() function

only scans numerical values. Is there any way to scan strings?



Yes, it is called scan() which is for arbitrary data, including
character(). See ?scan. Otherwise, you may want to look into ?readLines

Uwe Ligges





--
Siddharth Arun,
4th Year Undergraduate student
Industrial Engineering and Management,
IIT Kharagpur

[[alternative HTML version deleted]]

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




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


[R] Date format issue.

2011-06-28 Thread Ashim Kapoor
Dear R helpers,

I have 2 questions : -

1. My excel sheet has a column with dates like 01/03/1980 which is formatted
as 03/80 when I read this into R it reads as Mar-80. How can I read it in
the source format ?

2.

 v-c(Mar-80)
 as.Date(v,format=%b-%y)
[1] NA


Could someone please tell me where I am wrong in this date format conversion
?

Thank you all very much.

[[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] Date format issue.

2011-06-28 Thread jim holtman
You get the NA since it is indeterminate as to the date; paste on a 1
for the day

 v-c(Mar-80)
 as.Date(paste(v, '1'),format=%b-%y %d)
[1] 1980-03-01


On Tue, Jun 28, 2011 at 8:04 AM, Ashim Kapoor ashimkap...@gmail.com wrote:
 Dear R helpers,

 I have 2 questions : -

 1. My excel sheet has a column with dates like 01/03/1980 which is formatted
 as 03/80 when I read this into R it reads as Mar-80. How can I read it in
 the source format ?

 2.

 v-c(Mar-80)
 as.Date(v,format=%b-%y)
 [1] NA


 Could someone please tell me where I am wrong in this date format conversion
 ?

 Thank you all very much.

        [[alternative HTML version deleted]]

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




-- 
Jim Holtman
Data Munger Guru

What is the problem that you are trying to solve?

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


Re: [R] Date format issue.

2011-06-28 Thread Ashim Kapoor
  You get the NA since it is indeterminate as to the date; paste on a 1

 for the day


 Alright Jim,
Many thanks.

[[alternative HTML version deleted]]

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


[R] Using Match in a lookup table

2011-06-28 Thread Michael Pearmain
Hi All,

I'm having a few problems using match and a lookup table, previous Googling
show numerous solutions to matching a lookup table to a dataset,
My situation is slightly different as i have multiple lookup tables, (that i
cannot merge - for integrity reasons) that i wish to match against my data,
and each of these files is large, so lots of for / if conditions are not
ideal. (withstanding that i have multiple files of course)


For example,
I have data:

 v - c(foo, foo, bar, bar, baz)
 id - c(1,2)
 id2 - c(3)
 name - c(foo, bar)
 name2 - c(baz)
 df1 - data.frame(id, name)
 df2 - data.frame(id2, name2)

 v - df1$id[match(v,df1$name)]
 v
[1]  1  1  2  2 NA

So here i actually want to return
 v
[1]  1  1  2  2 baz

so next time i can run
v - df2$id[match(v,df2$name)]

and return:
 v
[1]  1  1  2  2 3

Any help very much appreciated

Mike

[[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] coxph() - unexpected result using Crawley's seedlings data (The R Book)

2011-06-28 Thread Jacob Brogren
Hi,

I ran the example on pp. 799-800 from Machael Crawley's The R Book using 
package survival v. 2.36-5, R 2.13.0 and RStudio 0.94.83. The model is a Cox's 
Proportional Hazards model. The result was quite different compared to the R 
Book. I have compared my code to the code in the book but can not find any 
differences in the function call. My results are attached as well as a link to 
the results presented in the book (link to Google Books).

When running the examples on pp. 797-799 I can't detect any differences in 
results so I don't think there are errors in the data set or in the creation of 
the status variable.

-
Original from the R Book:
http://books.google.com/books?id=8D4HVx0apZQClpg=PA799ots=rQgd_8ofeSdq=r%20coxph%20crawleypg=PA799#v=onepageqf=false

-
My result:
 summary(model1)
Call:
coxph(formula = Surv(death, status) ~ strata(cohort) * gapsize, 
data = seedlings)

  n= 60, number of events= 60 

coef exp(coef)  se(coef)  z 
Pr(|z|)
gapsize-0.001893  0.998109  0.593372 -0.003
0.997
gapsize:strata(cohort)cohort=September  0.717407  2.049112  0.860807  0.833
0.405

   exp(coef) exp(-coef) lower .95 upper .95
gapsize   0.9981  1.0020.3120 3.193
gapsize:strata(cohort)cohort=September2.0491  0.4880.379211.074

Rsquare= 0.022   (max possible= 0.993 )
Likelihood ratio test= 1.35  on 2 df,   p=0.5097
Wald test= 1.32  on 2 df,   p=0.5178
Score (logrank) test = 1.33  on 2 df,   p=0.514

Anyone have an idea why this is occurring?

Kind Regards

Jacob
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] chull increase number of points

2011-06-28 Thread pete
  

Dear R-help, 

I am using the chull function to create a convex
hull of a series of about 20,000 data points. 

A 
pain is temporary, glory is forever! 
Powered by Linux. www.linux.org
Scanned for viruses using ClamAV. www.clamav.net.



[[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] Error in library (nls)

2011-06-28 Thread matotope
Hi everybody,
I'm not very experienced with R software. I have used it several times for
some of the population genetics analyses. I have problem with executing one
of the script. The script is created by another software called Gimlet and
it is aimed to calculate rarefaction curve in R software.
However, when I try to execute the script, it says: *Error in library(nls) :
there is no package called 'nls'.*

I used the script previously on another PC about 2 or 3 years ago without
problems. I tried to reinstall R but no help. I'm using the version 2.13.0. 

Any help would be kindly appreciated (I'm in hurry :o...)

matotope

--
View this message in context: 
http://r.789695.n4.nabble.com/Error-in-library-nls-tp3630143p3630143.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] connecting R and PostgreSQL

2011-06-28 Thread Aya74656
Dear R-helpers,

I'm an absolute beginner using both R and PostgreSQL, but now I have to work 
with both programs. I need to connect R and my Postgres-database, but every 
attempt so far has failed (I tried using the RpgSQL package as well as 
RdbiPgSQL, the first, following this manual 
(http://code.google.com/p/rpostgresql/) didn't find any drivers for the 
database (step no. 1) whereas the second doesn't work with R version 2.13). 

Could someone please be so kind to either provide a step-by-step instruction on 
how to make this connection work or direct me to a manual?

Thanks in advance.

Yours sincerly,

Marie
--

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


Re: [R] Error in library (nls)

2011-06-28 Thread Ben Bolker
matotope matotope at gmail.com writes:

 
 Hi everybody,
 I'm not very experienced with R software. I have used it several times for
 some of the population genetics analyses. I have problem with executing one
 of the script. The script is created by another software called Gimlet and
 it is aimed to calculate rarefaction curve in R software.
 However, when I try to execute the script, it says: *Error in library(nls) :
 there is no package called 'nls'.*
 
 I used the script previously on another PC about 2 or 3 years ago without
 problems. I tried to reinstall R but no help. I'm using the version 2.13.0. 
 
 Any help would be kindly appreciated (I'm in hurry :o...)

  nls has been moved to the stats package.
  You don't need to load the (no longer existent) 'nls' package to use it.
  Try commenting out that line in the script.

  Ben Bolker

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


Re: [R] chull increase number of points

2011-06-28 Thread Ben Bolker
pete pete at nevill.uk.net writes:

 I am using the chull function to create a convex
 hull of a series of about 20,000 data points. 

  And your question is ... ?

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


Re: [R] Date format issue.

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 8:04 AM, Ashim Kapoor wrote:


Dear R helpers,

I have 2 questions : -

1. My excel sheet has a column with dates like 01/03/1980 which is  
formatted

as 03/80 when I read this into R it reads as Mar-80.


You should be able to change that from the format menu in Excel. If  
you use a custom format of -mm-dd and then save as csv you get  
dates in %Y-%m-%d acceptable format for importing, namely your date  
which on the data entry line appears as 3/1/80 appears in the cell as  
1980-03-01, and likewise in the text output the same way, and could be  
input using a colClasses argument of Date.



How can I read it in
the source format ?

2.


v-c(Mar-80)
as.Date(v,format=%b-%y)

[1] NA




Could someone please tell me where I am wrong in this date format  
conversion

?


--

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] coxph() - unexpected result using Crawley's seedlings data (The R Book)

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 6:51 AM, Jacob Brogren wrote:


Hi,

I ran the example on pp. 799-800 from Machael Crawley's The R Book  
using package survival v. 2.36-5, R 2.13.0 and RStudio 0.94.83. The  
model is a Cox's Proportional Hazards model. The result was quite  
different compared to the R Book. I have compared my code to the  
code in the book but can not find any differences in the function  
call. My results are attached as well as a link to the results  
presented in the book (link to Google Books).


Shouldn't this instead go to Mr Crawley? Despite his unfortunately  
successful efforts to give the impression that his book is  
authoritative and somehow official, it is neither. I have never seen a  
contribution to the R effort from Mr. Crawley.


--
David.




When running the examples on pp. 797-799 I can't detect any  
differences in results so I don't think there are errors in the data  
set or in the creation of the status variable.


-
Original from the R Book:
http://books.google.com/books?id=8D4HVx0apZQClpg=PA799ots=rQgd_8ofeSdq=r%20coxph%20crawleypg=PA799#v 
=onepageqf=false


-
My result:

summary(model1)

Call:
coxph(formula = Surv(death, status) ~ strata(cohort) * gapsize,
   data = seedlings)

 n= 60, number of events= 60

   coef exp(coef)   
se(coef)  z Pr(|z|)
gapsize-0.001893  0.998109  0.593372  
-0.0030.997
gapsize:strata(cohort)cohort=September  0.717407  2.049112   
0.860807  0.8330.405


  exp(coef) exp(-coef) lower .95  
upper .95
gapsize   0.9981  1.002 
0.3120 3.193
gapsize:strata(cohort)cohort=September2.0491  0.488 
0.379211.074


Rsquare= 0.022   (max possible= 0.993 )
Likelihood ratio test= 1.35  on 2 df,   p=0.5097
Wald test= 1.32  on 2 df,   p=0.5178
Score (logrank) test = 1.33  on 2 df,   p=0.514

Anyone have an idea why this is occurring?

Kind Regards

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


David Winsemius, MD
West Hartford, CT

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


[R] A masked function is a masked function by any other name

2011-06-28 Thread Juan Carlos Borrás
Dear all,
It looks like I do not grasp the concept of masked functions enough as
to solve this trivial problem.
The code that replicates the problem (a source code tree that realizes
a R package actually) is under github so one can call it clone it
easily from the command line (though more experienced users will spot
the problem by browsing through the package code):
git clone http://jcbor...@github.com/jcborras/rseedpkg.git

rseedpkg builds and installs with the usual sequence:
R CMD build rseedpkg
R CMD INSTALL rseedpkg_0.01-1.tar.gz

Last but not least one can test it from the command line:
Rscript --verbose --default-packages=testthat,log4r -e
test_package('rseedpkg')

The thing is that if one changes the call log4r:::debug() to plain
debug() in R/f1.r and R/f2.r then one ends up calling base:::debug()
and not log4r:::debug() even though the former should be masked by the
later as log4r is a package dependency of my dummy rseedpkg. And
that's something that I cannot understand...

Thanks in advance
jcb!

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

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 6:18 AM, Michael Pearmain wrote:


Hi All,

I'm having a few problems using match and a lookup table, previous  
Googling

show numerous solutions to matching a lookup table to a dataset,
My situation is slightly different as i have multiple lookup tables,  
(that i
cannot merge - for integrity reasons) that i wish to match against  
my data,
and each of these files is large, so lots of for / if conditions are  
not

ideal. (withstanding that i have multiple files of course)


For example,
I have data:


v - c(foo, foo, bar, bar, baz)
id - c(1,2)
id2 - c(3)
name - c(foo, bar)
name2 - c(baz)
df1 - data.frame(id, name)
df2 - data.frame(id2, name2)



v - df1$id[match(v,df1$name)]
v

[1]  1  1  2  2 NA


A numeric vector.


So here i actually want to return

v

[1]  1  1  2  2 baz


Not possibly a numeric vector.



so next time i can run
v - df2$id[match(v,df2$name)]

and return:

v

[1]  1  1  2  2 3


You need to keep track of the successful matches in df1 and then ypu  
probably want to interleave them with matches in df2. Perhaps instead  
use ifelse to do the work for you:


 ifelse(!is.na(match(v,df1$name)), match(v,df1$name),  
match(v,df2$name2) )

[1] 1 1 2 2 1




Any help very much appreciated

Mike

[[alternative HTML version deleted]]

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


David Winsemius, MD
West Hartford, CT

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


Re: [R] A masked function is a masked function by any other name

2011-06-28 Thread Duncan Murdoch

On 28/06/2011 9:38 AM, Juan Carlos Borrás wrote:

Dear all,
It looks like I do not grasp the concept of masked functions enough as
to solve this trivial problem.
The code that replicates the problem (a source code tree that realizes
a R package actually) is under github so one can call it clone it
easily from the command line (though more experienced users will spot
the problem by browsing through the package code):
git clone http://jcbor...@github.com/jcborras/rseedpkg.git

rseedpkg builds and installs with the usual sequence:
R CMD build rseedpkg
R CMD INSTALL rseedpkg_0.01-1.tar.gz

Last but not least one can test it from the command line:
Rscript --verbose --default-packages=testthat,log4r -e
test_package('rseedpkg')

The thing is that if one changes the call log4r:::debug() to plain
debug() in R/f1.r and R/f2.r then one ends up calling base:::debug()
and not log4r:::debug() even though the former should be masked by the
later as log4r is a package dependency of my dummy rseedpkg. And
that's something that I cannot understand...


If you are using ::: (three colons), then you may be looking into the 
unexported functions in log4r.  The only normal way to see unexported 
functions is to use three colons.


Duncan Murdoch

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


Re: [R] Error in library (nls)

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 8:09 AM, matotope wrote:


Hi everybody,
I'm not very experienced with R software. I have used it several  
times for
some of the population genetics analyses. I have problem with  
executing one
of the script. The script is created by another software called  
Gimlet and

it is aimed to calculate rarefaction curve in R software.
However, when I try to execute the script, it says: *Error in  
library(nls) :

there is no package called 'nls'.*


Did you check to see whether nls was installed under 2.13.0? (actually  
you did check with the library call and it isn't.)




I used the script previously on another PC about 2 or 3 years ago  
without
problems. I tried to reinstall R but no help. I'm using the version  
2.13.0.


Any help would be kindly appreciated (I'm in hurry :o...)

matotope

--
View this message in context: 
http://r.789695.n4.nabble.com/Error-in-library-nls-tp3630143p3630143.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.


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] coxph() - unexpected result using Crawley's seedlings data (The R Book)

2011-06-28 Thread Robert A LaBudde

Did you create the 'status' variable the way indicated on p. 797?

Frequently with Surv() it pays to use syntax such as Surv(death, 
status==1) to make a clear logical statement of what is an event 
(status==1) vs. censored.


PS. Next time include head(seedlings) and str(seedlings) to make 
clear what you are using as data.



At 06:51 AM 6/28/2011, Jacob Brogren wrote:

Hi,

I ran the example on pp. 799-800 from Machael Crawley's The R Book 
using package survival v. 2.36-5, R 2.13.0 and RStudio 0.94.83. The 
model is a Cox's Proportional Hazards model. The result was quite 
different compared to the R Book. I have compared my code to the 
code in the book but can not find any differences in the function 
call. My results are attached as well as a link to the results 
presented in the book (link to Google Books).


When running the examples on pp. 797-799 I can't detect any 
differences in results so I don't think there are errors in the data 
set or in the creation of the status variable.


-
Original from the R Book:
http://books.google.com/books?id=8D4HVx0apZQClpg=PA799ots=rQgd_8ofeSdq=r%20coxph%20crawleypg=PA799#v=onepageqf=false

-
My result:
 summary(model1)
Call:
coxph(formula = Surv(death, status) ~ strata(cohort) * gapsize,
data = seedlings)

  n= 60, number of events= 60

coef 
exp(coef)  se(coef)  z Pr(|z|)
gapsize-0.001893  0.998109  0.593372 
-0.0030.997
gapsize:strata(cohort)cohort=September  0.717407  2.049112  0.860807 
 0.8330.405


   exp(coef) exp(-coef) lower 
.95 upper .95
gapsize   0.9981  1.002 
0.3120 3.193
gapsize:strata(cohort)cohort=September2.0491  0.488 
0.379211.074


Rsquare= 0.022   (max possible= 0.993 )
Likelihood ratio test= 1.35  on 2 df,   p=0.5097
Wald test= 1.32  on 2 df,   p=0.5178
Score (logrank) test = 1.33  on 2 df,   p=0.514

Anyone have an idea why this is occurring?

Kind Regards

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



Robert A. LaBudde, PhD, PAS, Dpl. ACAFS  e-mail: r...@lcfltd.com
Least Cost Formulations, Ltd.URL: http://lcfltd.com/
824 Timberlake Drive Tel: 757-467-0954
Virginia Beach, VA 23464-3239Fax: 757-467-2947

Vere scire est per causas scire

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


Re: [R] connecting R and PostgreSQL

2011-06-28 Thread Gabor Grothendieck
On Tue, Jun 28, 2011 at 8:15 AM,  aya74...@gmx.de wrote:
 Dear R-helpers,

 I'm an absolute beginner using both R and PostgreSQL, but now I have to work 
 with both programs. I need to connect R and my Postgres-database, but every 
 attempt so far has failed (I tried using the RpgSQL package as well as 
 RdbiPgSQL, the first, following this manual 
 (http://code.google.com/p/rpostgresql/) didn't find any drivers for the 
 database (step no. 1) whereas the second doesn't work with R version 2.13).

 Could someone please be so kind to either provide a step-by-step instruction 
 on how to make this connection work or direct me to a manual?


For RpgSQL there is installation info if you click on the Installation
link here:

http://cran.r-project.org/web/packages/RpgSQL/index.html

and the same info is in the INSTALL file in the package.

-- 
Statistics  Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

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


Re: [R] About the covariant

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 12:57 AM, Lao Meng wrote:


Thanks David for your reply.
You said a single slope and intercept are estimated for each  
variable.Actually I can only get one intercept no matter how many


Sorry. You are right. You get individual slopes (and differences for  
factors) reference to a single intercept (unless you use different  
formula specification)




variables exist,but a slope for each variable.

Since the regression is done via:lm(CD4 ~ time + gender + income)
It seems that the explanatory variable(time) and the two  
covariants(gender,income) are treated in the same way,but I think  
explanatory variable and covariant should be treated differently


I do not understand what you are saying when you use the word  
'differently' and increasing the number of times you say it is not  
improving communication.



although I don't know how to do it.
Also,they are not both numeric,if gender are F(Female) and  
M(Male),and income are L(Low),M(median),H(High).


Yes. discrete, unordered factors can have associated estimated  
effects, which will be differences from the intercept level. The  
intercept in you case would probably be Female/High, since the default  
ordering of factor levels is alphabetic. How are these multiple  
question arising? Are you in the middle of an introductory regression  
class?


--
David




2011/6/28 David Winsemius dwinsem...@comcast.net

On Jun 27, 2011, at 10:02 PM, Lao Meng wrote:

Hi all,I have some questions about the covariants of regression.

My target: To explore the trend of CD4 level through a period of time.

Response variable: CD4 count
Explanatory variable:time

Also, the demology information is available,such as  
gender,occupation,income

level...

Q1,Are these variables of demology information called covariant?
Q2,How can I correct the impact of covariant so that I can get the
corrected result of CD4's change through the time period?
Q3,How to treat the covariants in regression?I've looked up to many  
papers

of R on regression,which treat the covariant in the same

way as the Explanatory variable,like following:
lm(CD4 ~ time + gender + income)

Yes that seems pretty standard practice. It does, of course, force  
the relationships to a) be linear and b) means that a single slope  
and intercept are estimated for each variable, neither of a} or b}  
assumptions may be true.



From above expression of regression,it's obvious that the response  
variables

and covariants are treated the same way,

In what sense are you making that claim? True they are both numeric,  
but what else are you saying?


--
David

but acturally

they are totally different.



Thanks for your help.

My best.

   [[alternative HTML version deleted]]

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

David Winsemius, MD
West Hartford, CT




David Winsemius, MD
West Hartford, CT

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


Re: [R] Creating a Polar Plot with expanding points as radius increases

2011-06-28 Thread Patrick Jemison
Hi Jim,

That sounds pretty great!  I am happy to have contributed a stimulus
for action to be taken in further developing the tools.  I'll keep an
eye out for your update.

Thanks,
Patrick

On Tue, Jun 28, 2011 at 6:29 AM, Jim Lemon j...@bitwrit.com.au wrote:
 On 06/28/2011 08:35 AM, Patrick Jemison wrote:

 I'd like to create a polar plot similar to those created by the polarFreq
 function in the openair package.  However, this package seems to be
 specific
 to wind speed and direction, and requires a ws (wind speed) and a wd
 (wind direction) column.  My data is unrelated to wind speed, but I'd like
 to be able to get a plot that does what polarFreq's plots do; I'd like to
 have points that expand as their distance from the center increases.

 http://www.oga-lab.net/RGM2/func.php?rd_id=openair:polarFreq

 Without this property, my plot looks more like spokes on a wheel, or
 similar
 to the polar.plot function in the plotrix package:
 http://www.oga-lab.net/RGM2/func.php?rd_id=plotrix:polar.plot

 Is there a package or function available that allows the creation of a
 plot
 as described above for a generic set of polar data?  I spent a few hours
 today trying to find something I could use, but no luck yet.

 I'm using data of the form: x = r, y = theta, and z = a value to plot at
 the
 corresponding r,theta (as a comparison, for the polarFreq plot, z would be
 the frequency that is indicated by various colors in the plot).

 Hi Patrick,
 polar.plot (or other plots in the radial.plot family) currently don't plot
 sectors or fill the (what the heck is the intersection of a sector and an
 annulus called?) to illustrate a third value (say, proportion) in addition
 to direction and intensity. However, it can be done, and I might just have a
 shot at it if I can scratch up an hour or two sometime in the next couple of
 years. Well, maybe sooner than that. Look out for radial.pie.

 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 required for GO Annotation problem

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 5:34 AM, suman pal wrote:


Hello,
I basically want to use R-help, and post some problems which I am  
facing. The Ref is a well known Genome Biology paper Bioconductor:  
open software development for computational biology and  
bioinformatics by Robert C Gentleman et al., 2004. Generating  
Heatmaps till Fig2 is working so I think esetSel is not the problem..
However, for generating the Figure 3, for GO annotations the  
following command is not working:

ll - mget(geneNames(esetSel),hgu95av2LOCUSID)
#it is displaying error messages  Error in mget(geneNames(esetSel),  
hgu95av2LOCUSID) : object 'hgu95av2LOCUSID' not foundand also  
geneNames not found try featureNames instead


You seem to be missing a data object. A google search produces this as  
the third hit that suggest you are attempting to run outdate ... well  
it is 7 years old ... code:


http://www.rforge.net/bin/results/compton/2.1/data/annotation/hgu95av2-00check.html

Hence I cant proceed to the next set of commands provided in the  
paper which are as follows:
ll - unique(unlist(ll)) mf - as.data.frame(GOHyperG(ll))[, 1:3]  
mf - mf[mf$pvalue  0.01, ]

sincerelySumanCCMB Hyderabad.

This would seem a natural question for the BoConductor mailing list...  
after you do a bit more basic searching.




[[alternative HTML version deleted]]

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


David Winsemius, MD
West Hartford, CT

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


Re: [R] Using Match in a lookup table

2011-06-28 Thread David Winsemius

On Jun 28, 2011, at 10:29 AM, Michael Pearmain wrote:

 Thanks for the idea David,
 My problem comes from having (say) upto 10 different match files, so  
 nested ifelse whilst would work doesn't seem and elegant solution,

Well, you cannot  mix text and numeric modes in a vector as you  
appeared to expect.

There may be a limit on the nesting if ifelse's as well. ISTR 7 but  
have not been able to document that vague memory.

It seems that separate matching runs followed by some logic that  
cbind()-ed and then collapsed the results perhaps by just dropping  
all the NA's row-wise.

-- 
David.



 However if needs must..

 Mike

 On 28 June 2011 14:39, David Winsemius dwinsem...@comcast.net wrote:

 On Jun 28, 2011, at 6:18 AM, Michael Pearmain wrote:

 Hi All,

 I'm having a few problems using match and a lookup table, previous  
 Googling
 show numerous solutions to matching a lookup table to a dataset,
 My situation is slightly different as i have multiple lookup tables,  
 (that i
 cannot merge - for integrity reasons) that i wish to match against  
 my data,
 and each of these files is large, so lots of for / if conditions are  
 not
 ideal. (withstanding that i have multiple files of course)


 For example,
 I have data:

 v - c(foo, foo, bar, bar, baz)
 id - c(1,2)
 id2 - c(3)
 name - c(foo, bar)
 name2 - c(baz)
 df1 - data.frame(id, name)
 df2 - data.frame(id2, name2)

 v - df1$id[match(v,df1$name)]
 v
 [1]  1  1  2  2 NA

 A numeric vector.


 So here i actually want to return
 v
 [1]  1  1  2  2 baz

 Not possibly a numeric vector.



 so next time i can run
 v - df2$id[match(v,df2$name)]

 and return:
 v
 [1]  1  1  2  2 3

 You need to keep track of the successful matches in df1 and then ypu  
 probably want to interleave them with matches in df2. Perhaps  
 instead use ifelse to do the work for you:

  ifelse(!is.na(match(v,df1$name)), match(v,df1$name),  
 match(v,df2$name2) )

 [1] 1 1 2 2 1



 Any help very much appreciated

 Mike

[[alternative HTML version deleted]]

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

 David Winsemius, MD
 West Hartford, CT



David Winsemius, MD
West Hartford, CT


[[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] chull increase number of points

2011-06-28 Thread pete
 

Dear R-help, 

I am using the chull function to create a convex hull of a series of about
20,000 data points. 

A 
[[alternative HTML version deleted]]

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


Re: [R] lattice multiple y-scale possible?

2011-06-28 Thread Steven R Corsi

Dennis

Adding in y = list(relation = 'free') to the scales argument worked  
very well.


Thanks!
Steve

===
Steven R. CorsiPhone: (608) 821-3835
Research Hydrologist   email: srco...@usgs.gov
U.S. Geological Survey
Wisconsin Water Science Center
8505 Research Way
Middleton, WI 53562
===


On 6/28/2011 6:23 AM, Dennis Murphy wrote:

Hi:

bwplot(dat~site|parameter,data=mydat,
  layout=c(1,5),
  cex=2,
  xlab=Site Name,
  ylab=,
  labels=levels(mydat$site),
  scales=list(tick.number=list(4), rot = c(0, 90),
  y = list(relation = 'free')))

Does that work?

Dennis

On Mon, Jun 27, 2011 at 10:46 PM, Steve Corsisrco...@usgs.gov  wrote:

Hi
I am attempting to use the lattice bwplot function to generate boxplots of
numerous parameters (1-panel/parameter) by site (x-axis). The parameters
have quite different ranges of values, so it would be best to have a
separate y-axis range for each panel. Below is a basic example of what I am
trying to do. As is seen, the y-axes need to be scaled individually to make
this useful. Any information on how to do this would be much appreciated:

rm(mydat)
rm(tempdf)

for (i in 1:5){
for (ii in 1:5){
dat- sample(1:20)*10^ii
tempdf- data.frame(dat)
tempdf$parameter- paste(parameter ,ii,sep=)
tempdf$site- paste(site,i,sep=)
if(!exists(mydat)) {mydat- tempdf
  }else {mydat- rbind(mydat,tempdf)}
  }
}

bwplot(dat~site|parameter,data=mydat,
  layout=c(1,5),
  cex=2,
  xlab=Site Name,
  ylab=,
  labels=levels(mydat$site), scales=list(tick.number=list(4)))

thanks
Steve

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



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


Re: [R] Error in library (nls)

2011-06-28 Thread matotope
Hi Ben,
thanks for your advise. I put the comment in front of the line with nls. Now
the message doesn't appear but the script won't produce the result files
which I need (the result files are empty). 
When executing the script, R asks for number of iterations and after that it
states Read 20 items but nothing happens. Probably also some other part of
the script needs to be amended.

Here is the script:

#library (nls)
options (digits = 6, show.error.messages=FALSE)

cat (Enter the number of iterations: )
niter-as.integer (readLines(n = 1))
cat (Parameter a (asymptote) for Kohn's equation, file=Param_a_Kohn.txt,
append=TRUE, fill=TRUE)
cat (Number of iterations: ,niter, file=Param_a_Kohn.txt, append=TRUE,
fill=TRUE)
cat ( Niter  Min Mean Max Sd,
file=Param_a_Kohn.txt,append=TRUE, fill=TRUE)
cat (Parameter b for Kohn's equation, file=Param_b_Kohn.txt,
append=TRUE, fill=TRUE)
cat (Number of iterations: ,niter, file=Param_b_Kohn.txt, append=TRUE,
fill=TRUE)
cat ( Niter  Min Mean Max Sd,
file=Param_b_Kohn.txt,append=TRUE, fill=TRUE)
cat (Parameter a (asymptote) for Eggert's equation,
file=Param_a_Eggert.txt, append=TRUE, fill=TRUE)
cat (Number of iterations: ,niter, file=Param_a_Eggert.txt, append=TRUE,
fill=TRUE)
cat ( Niter  Min Mean Max Sd,
file=Param_a_Eggert.txt,append=TRUE, fill=TRUE)
cat (Parameter b for Eggert's equation, file=Param_b_Eggert.txt,
append=TRUE, fill=TRUE)
cat (Number of iterations: ,niter, file=Param_b_Eggert.txt, append=TRUE,
fill=TRUE)
cat ( Niter  Min Mean Max Sd,
file=Param_b_Eggert.txt,append=TRUE, fill=TRUE)

#construction of the dataset to estimate the parameter from output file from
GEMINI
createdata-function(n, niter=500) {
concent-function(echa) apply(matrix(1
:length(echa),ncol=1),1,flocal1,pop=echa)
flocal1-function(pop,j) length(unique(pop[1:j]))
indsel-length (n)
if (indsel=1) stop (n1 expected)
faesel-sum(n)
echa-rep(1:indsel,n)
return(echa,faesel)
}


#function to estimate and print the results
calculate-function(data){
dataset-createdata(data, niter)
coeffaKohn-rep(0, niter)
coeffbKohn-rep(0, niter)
niterKohn-0
coeffaEggert-rep(0, niter)
coeffbEggert-rep(0, niter)
niterEggert-0

draw-array(0,c(niter+1, dataset$faesel))
draw1-array(0,c(niter+1, dataset$faesel))
draw2-array(0,c(niter+1, dataset$faesel))
draw3-array(0,c(niter+1, dataset$faesel))

for (k in 1:niter) {
w-rep(0,dataset$faesel)
w-w+concent(sample(dataset$echa,dataset$faesel,replace=F))
faesel-length (w)
dxy-cbind.data.frame(1:faesel,w)
names(dxy)-c(x,y)

draw[k,1:faesel]-w[1:faesel]

ok-0
while(ok5){
tmp1-NA
tmp2-NA
tmp-try(coef(nls1-nls(y~a*x/(b+x), data=dxy,
start=list(a=max(dxy$y),b=1
tmp1-tmp [1]
tmp2-tmp [2]
if (is.character(tmp1)) {coeffaKohn[k]-NA
ok-ok+1}
else {coeffaKohn[k]-tmp1
niterKohn-niterKohn + 1
draw1[k,1:faesel]-predict(nls1)
ok-5}
if (is.character(tmp2)) {coeffbKohn[k]-NA}
else {coeffbKohn[k]-tmp2}
}

ok-0
while(ok5){
tmp1-NA
tmp2-NA
tmp-try(coef(nls3-nls(y~a*(1-exp(b*x)), data=dxy,
start=list(a=max(dxy$y),b=-1
tmp1-tmp [1]
tmp2-tmp [2]
if (is.character(tmp1)) {coeffaEggert[k]-NA
ok-ok+1}
else {coeffaEggert[k]-tmp1
niterEggert-niterEggert + 1
draw3[k,1:faesel]-predict(nls3)
ok-5}
if (is.character(tmp2)) {coeffbEggert[k]-NA}
else {coeffbEggert[k]-tmp2}
}
}

x11()
plot(draw[1,], main=Number of unique genotypes against number of feces
analysed, sub=observed= circles; mean of observed= black line; Kohn's eq=
red line; Eggert's eq= blue line, xlab=Number of feces analysed,
ylab=Number of unique genotypes)
for (k in 2:niter) {
lines(draw[k,], type=o)}
for (k in 1:faesel) {
draw[niter+1,k]-mean(draw[1:niter,k])
draw1[niter+1,k]-mean(draw1[1:niter,k])
draw2[niter+1,k]-mean(draw2[1:niter,k])
draw3[niter+1,k]-mean(draw3[1:niter,k])
}
lines(draw[niter+1,], type=l, col=black, lwd=2)
lines(draw1[niter+1,], type=l, col=red, lwd=2)
lines(draw3[niter+1,], type=l, col=blue, lwd=2)

minaKohn-min(coeffaKohn, na.rm=TRUE)
maxaKohn-max(coeffaKohn, na.rm=TRUE)
meanaKohn-mean(coeffaKohn, na.rm=TRUE)
sdaKohn-sd(coeffaKohn, na.rm=TRUE)
minbKohn-min(coeffbKohn, na.rm=TRUE)
maxbKohn-max(coeffbKohn, na.rm=TRUE)
meanbKohn-mean(coeffbKohn, na.rm=TRUE)
sdbKohn-sd(coeffaKohn, na.rm=TRUE)
minaEggert-min(coeffaEggert, na.rm=TRUE)
maxaEggert-max(coeffaEggert, na.rm=TRUE)
meanaEggert-mean(coeffaEggert, na.rm=TRUE)
sdaEggert-sd(coeffaEggert, na.rm=TRUE)
minbEggert-min(coeffbEggert, na.rm=TRUE)
maxbEggert-max(coeffbEggert, na.rm=TRUE)
meanbEggert-mean(coeffbEggert, na.rm=TRUE)
sdbEggert-sd(coeffbEggert, na.rm=TRUE)

pop-

cat(pop,niterKohn, minaKohn, meanaKohn, maxaKohn, sdaKohn, sep=  , file =
Param_a_Kohn.txt,append=TRUE, fill = TRUE)
cat(pop,niterKohn, minbKohn, meanbKohn, maxbKohn, sdbKohn, sep=  , file =
Param_b_Kohn.txt,append=TRUE, fill = TRUE)

cat(pop,niterEggert, minaEggert, meanaEggert, maxaEggert, sdaEggert, sep= 
, file = Param_a_Eggert.txt,append=TRUE, fill = TRUE)
cat(pop,niterEggert, minbEggert, meanbEggert, maxbEggert, 

Re: [R] Using Match in a lookup table

2011-06-28 Thread Michael Pearmain
Thanks for the idea David,
My problem comes from having (say) upto 10 different match files, so nested
ifelse whilst would work doesn't seem and elegant solution,

However if needs must..

Mike

On 28 June 2011 14:39, David Winsemius dwinsem...@comcast.net wrote:


 On Jun 28, 2011, at 6:18 AM, Michael Pearmain wrote:

  Hi All,

 I'm having a few problems using match and a lookup table, previous
 Googling
 show numerous solutions to matching a lookup table to a dataset,
 My situation is slightly different as i have multiple lookup tables, (that
 i
 cannot merge - for integrity reasons) that i wish to match against my
 data,
 and each of these files is large, so lots of for / if conditions are not
 ideal. (withstanding that i have multiple files of course)


 For example,
 I have data:

  v - c(foo, foo, bar, bar, baz)
 id - c(1,2)
 id2 - c(3)
 name - c(foo, bar)
 name2 - c(baz)
 df1 - data.frame(id, name)
 df2 - data.frame(id2, name2)


  v - df1$id[match(v,df1$name)]
 v

 [1]  1  1  2  2 NA


 A numeric vector.


 So here i actually want to return

 v

 [1]  1  1  2  2 baz


 Not possibly a numeric vector.



 so next time i can run
 v - df2$id[match(v,df2$name)]

 and return:

 v

 [1]  1  1  2  2 3


 You need to keep track of the successful matches in df1 and then ypu
 probably want to interleave them with matches in df2. Perhaps instead use
 ifelse to do the work for you:

  ifelse(!is.na(match(v,df1$**name)), match(v,df1$name),
 match(v,df2$name2) )

 [1] 1 1 2 2 1



 Any help very much appreciated

 Mike

[[alternative HTML version deleted]]

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


 David Winsemius, MD
 West Hartford, CT



[[alternative HTML version deleted]]

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


Re: [R] Error in library (nls)

2011-06-28 Thread matotope
Hi,
I installed the 2.8.0 version and it seems to be working fine now.
Thanks for the help guys!
matotope

--
View this message in context: 
http://r.789695.n4.nabble.com/Error-in-library-nls-tp3630143p3630552.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] problem with rJava : same as message from wwreith on Mon, 27 Jun 2011

2011-06-28 Thread Olivier ETERRADOSSI

Hi R-listers,

I run R version 2.13.0, on a i386-pc-mingw32/i386 platform under Windows XP.
I perform daily update of my installed packages.
I've got the most recent Java for this platform installed, as I am told 
by the JavaUpdater when checking.


(and I've read the posting guide, so that I'm not going to get fired by 
Uwe as poor wwreith [joke])...


But :

since one of the last daily updates (cannot tell exactly when because I 
spent some days without calling library(rJava), I get this error message 
(translated from french for the first part, in english for the second) :

Error : .onLoad failed in loadNamespace() for 'rJava', details :
 Call : fun(...)


error : No CurrentVersion entry in 'Software\JavaSoft\Java Runtime 
Environment'! Try re-installing Java and make sure R and Java have 
matching architectures.


So it seems to be quite the same as wwreith.

Can somebody help ?

Thanks, all the best to everybody. Olivier

--
Olivier ETERRADOSSI
Maître-Assistant
animateur du groupe Sensomines (Institut Carnot M.I.N.E.S)
-
CMGD
Pôle Matériaux Polymères Avancés
axe Propriétés Psycho-Sensorielles des Matériaux
-
Ecole des Mines d'Alès
Hélioparc, 2 av. P. Angot, F-64053 PAU CEDEX 9
tel std: +33 (0)5.59.30.54.25
tel direct: +33 (0)5.59.30.90.35
fax: +33 (0)5.59.30.63.68
http://www.mines-ales.fr
e-mail : olivier.eterrado...@mines-ales.fr

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


Re: [R] coxph() - unexpected result using Crawley's seedlings data (The R Book)

2011-06-28 Thread Jacob Brogren
Hi,

sorry about that; here is the full output - data set, structure, model and 
result.

Cheers

Jacob

 seedlings
  cohort death gapsize status
1  September 7  0.5889  1
2  September 3  0.6869  1
3  September12  0.1397  1
4  September 1  0.1921  1
5  September 4  0.2798  1
6  September 2  0.2607  1
7  September 6  0.9467  1
8  September 6  0.6375  1
9  September 8  0.1527  1
10 September 3  0.8237  1
11 September 1  0.5979  1
12 September 1  0.2914  1
13 September 3  0.5053  1
14 September 5  0.4714  1
15 September 2  0.6041  1
16 September 8  0.8812  1
17 September 4  0.8416  1
18 September 1  0.0577  1
19 September 2  0.9034  1
20 September 2  0.4348  1
21 September 7  0.9878  1
22 September11  0.1486  1
23 September14  0.5003  1
24 September 1  0.8507  1
25 September10  0.8187  1
26 September14  0.0291  1
27 September 1  0.3785  1
28 September 4  0.8384  1
29 September 2  0.8351  1
30 September 2  0.9674  1
31   October 1  0.6943  1
32   October 1  0.2591  1
33   October 2  0.7397  1
34   October 2  0.4663  1
35   October14  0.9115  1
36   October 5  0.1750  1
37   October 1  0.5628  1
38   October 8  0.2681  1
39   October 5  0.6967  1
40   October 2  0.7020  1
41   October 4  0.7971  1
42   October 3  0.4047  1
43   October 5  0.0498  1
44   October10  0.0364  1
45   October 9  0.4080  1
46   October 1  0.6226  1
47   October11  0.3002  1
48   October 3  0.8111  1
49   October21  0.4894  1
50   October 1  0.0375  1
51   October 4  0.2560  1
52   October 9  0.2168  1
53   October 8  0.7437  1
54   October 1  0.9082  1
55   October 3  0.9496  1
56   October 9  0.1040  1
57   October 9  0.8691  1
58   October16  0.9502  1
59   October 6  0.0790  1
60   October 1  0.5658  1
 str(seedlings)
'data.frame':   60 obs. of  4 variables:
 $ cohort : Factor w/ 2 levels October,September: 2 2 2 2 2 2 2 2 2 2 ...
 $ death  : int  7 3 12 1 4 2 6 6 8 3 ...
 $ gapsize: num  0.589 0.687 0.14 0.192 0.28 ...
 $ status : num  1 1 1 1 1 1 1 1 1 1 ...
 model1 - coxph(Surv(death,status)~strata(cohort)*gapsize,data=seedlings)
 summary(model1)
Call:
coxph(formula = Surv(death, status) ~ strata(cohort) * gapsize, 
data = seedlings)

  n= 60, number of events= 60 

coef exp(coef)  se(coef)  z 
Pr(|z|)
gapsize-0.001893  0.998109  0.593372 -0.003
0.997
gapsize:strata(cohort)cohort=September  0.717407  2.049112  0.860807  0.833
0.405

   exp(coef) exp(-coef) lower .95 upper .95
gapsize   0.9981  1.0020.3120 3.193
gapsize:strata(cohort)cohort=September2.0491  0.4880.379211.074

Rsquare= 0.022   (max possible= 0.993 )
Likelihood ratio test= 1.35  on 2 df,   p=0.5097
Wald test= 1.32  on 2 df,   p=0.5178
Score (logrank) test = 1.33  on 2 df,   p=0.514

28 jun 2011 kl. 15.48 skrev Robert A LaBudde:

 Did you create the 'status' variable the way indicated on p. 797?
 
 Frequently with Surv() it pays to use syntax such as Surv(death, status==1) 
 to make a clear logical statement of what is an event (status==1) vs. 
 censored.
 
 PS. Next time include head(seedlings) and str(seedlings) to make clear what 
 you are using as data.
 
 
 At 06:51 AM 6/28/2011, Jacob Brogren wrote:
 Hi,
 
 I ran the example on pp. 799-800 from Machael Crawley's The R Book using 
 package survival v. 2.36-5, R 2.13.0 and RStudio 0.94.83. The model is a 
 Cox's Proportional Hazards model. The result was quite different compared to 
 the R Book. I have compared my code to the code in the book but can not find 
 any differences in the function call. My results are attached as well as a 
 link to the results presented in the book (link to Google Books).
 
 When running the examples on pp. 797-799 I can't detect any differences in 
 results so I don't think there are errors in the data set or in the creation 
 of the status variable.
 
 -
 Original from the R Book:
 http://books.google.com/books?id=8D4HVx0apZQClpg=PA799ots=rQgd_8ofeSdq=r%20coxph%20crawleypg=PA799#v=onepageqf=false
 
 -
 My result:
  summary(model1)
 Call:
 coxph(formula = Surv(death, status) ~ strata(cohort) * gapsize,
data = seedlings)
 
  n= 60, number of events= 60
 
coef exp(coef)  se(coef)  z 
 Pr(|z|)
 gapsize-0.001893  0.998109  0.593372 -0.003  
   0.997
 

Re: [R] help required for GO Annotation problem

2011-06-28 Thread James W. MacDonald

Hi Suman,

On 6/28/2011 10:02 AM, David Winsemius wrote:


On Jun 28, 2011, at 5:34 AM, suman pal wrote:


Hello,
I basically want to use R-help, and post some problems which I am
facing. The Ref is a well known Genome Biology paper Bioconductor:
open software development for computational biology and
bioinformatics by Robert C Gentleman et al., 2004. Generating
Heatmaps till Fig2 is working so I think esetSel is not the problem..
However, for generating the Figure 3, for GO annotations the following
command is not working:

ll - mget(geneNames(esetSel),hgu95av2LOCUSID)

#it is displaying error messages Error in mget(geneNames(esetSel),
hgu95av2LOCUSID) : object 'hgu95av2LOCUSID' not foundand also
geneNames not found try featureNames instead


As David notes, this question is better asked on the BioC listserv. And 
it is really old code. The LOCUSID object has been deprecated for a long 
time now, replaced by ENTREZID, and as the warning notes, geneNames() 
has been deprecated in favor of featureNames.


So try

ll - mget(featureNames(esetSel), hgu95av2ENTREZID)

also GOHyperG() is gone as well, replaced by S4 methods that might be a 
bit less easy to figure out. I don't have the book you are working 
through, but you would be better served if you were to use the vignettes 
that come with the packages you are using instead (or at least as a way 
to get updated examples of how to use the package).


Anyway, please direct all future questions about BioC packages to the 
BioC listserv.


Best,

Jim




You seem to be missing a data object. A google search produces this as
the third hit that suggest you are attempting to run outdate ... well it
is 7 years old ... code:

http://www.rforge.net/bin/results/compton/2.1/data/annotation/hgu95av2-00check.html



Hence I cant proceed to the next set of commands provided in the paper
which are as follows:

ll - unique(unlist(ll)) mf - as.data.frame(GOHyperG(ll))[, 1:3]
mf - mf[mf$pvalue  0.01, ]

sincerelySumanCCMB Hyderabad.


This would seem a natural question for the BoConductor mailing list...
after you do a bit more basic searching.



[[alternative HTML version deleted]]

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


David Winsemius, MD
West Hartford, CT

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


--
James W. MacDonald, M.S.
Biostatistician
Douglas Lab
University of Michigan
Department of Human Genetics
5912 Buhl
1241 E. Catherine St.
Ann Arbor MI 48109-5618
734-615-7826
**
Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues 


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


Re: [R] What does class call mean? How do I make class formula into a call?

2011-06-28 Thread Rita Carreira


Thank you Bert and Prof. Ripley for your feedback. I did read the language 
documentation and it was not entirely clear to me, but I'm one of those people 
that has to read and digest something before it clicks. However, I did realize 
that the issue with calland formula was not the real reason why my program 
did not work. The real reason was much more trivial: I put the arguments inside 
the systemfit function out of order. Eventually, I figured it out. The good 
thing about this is that I learned about the existence of the R language 
documentation. 
Thank you again both!
Rita
=
If you think education is expensive, try ignorance.--Derek Bok



 Date: Sat, 25 Jun 2011 05:59:40 +0100
 From: rip...@stats.ox.ac.uk
 To: gunter.ber...@gene.com
 CC: ritacarre...@hotmail.com
 Subject: Re: [R] What does class call mean? How do I make class formula 
 into a call?

 This is really a misleading subject: it is already a call! From
 ?class

 Many R objects have a ‘class’ attribute, a character vector giving
 the names of the classes from which the object _inherits_. If the
 object does not have a class attribute, it has an implicit class,
 ‘matrix’, ‘array’ or the result of ‘mode(x)’ (except that
 integer vectors have implicit class ‘integer’).

 So, simply remove the class if you want the mode: but anything which
 needs to know this is call will be looking at the mode and not the
 class.

  zz - ~x
  class(zz)
 [1] formula
  mode(zz)
 [1] call

 And see ?mode and ?call. Formulae and calls which are not formulae
 are completely different: you cannot coerce one to the other.


 On Fri, 24 Jun 2011, Bert Gunter wrote:

  Well, this is kind of complicated. The first place you should go for
  help is not this list, but the R docs. Specfically ?call. This
  assumes familiarity with R's (S3) class system and language structure,
  however.. For this, I suggest ?UseMethod and consulting the R Language
  Definition Manual.
 
  Perhaps some brave soul on this list will attempt a short explanation
  in reply. But I am not (s)he.
 
  Cheers,
  Bert
 
  Oh -- as for specific suggestions, I think you need to do what the
  posting guide asks and provide a minimal reproducible example to give
  people a clearer idea of what's going on.
 
  On Fri, Jun 24, 2011 at 2:58 PM, Rita Carreira ritacarre...@hotmail.com 
  wrote:
 
  I have a list called tabs that I would like to have the same
  structure as my list eqSystem. The two look like they have the
  same format but they are different because when I look at their
  attributes, class(eqSystem[[1]]) is call but class(tabs[[1]]) is
  formula. I want to have class(tabs[[1]]) as a call too. So what
  does call mean? And how do I make an object of type formula be
  of type call?
  Thank you so much!!!--Rita
  class(tabs)
  [1] list
  class(tabs[1])
  [1] list
  class(tabs[[1]])
  [1] formula class(eqSystem)
  [1] list
  class(eqSystem[1])
  [1] list
  class(eqSystem[[1]])
  [1] call
 
 
  Rita
  =
  If you think education is expensive, try ignorance.--Derek Bok
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide 
  http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 
 
 
  --
  Men by nature long to get on to the ultimate truths, and will often
  be impatient with elementary studies or fight shy of them. If it were
  possible to reach the ultimate truths without the elementary studies
  usually prefixed to them, these would not be preparatory studies but
  superfluous diversions.
 
  -- Maimonides (1135-1204)
 
  Bert Gunter
  Genentech Nonclinical Biostatistics
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 

 --
 Brian D. Ripley, rip...@stats.ox.ac.uk
 Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
 University of Oxford, Tel: +44 1865 272861 (self)
 1 South Parks Road, +44 1865 272866 (PA)
 Oxford OX1 3TG, UK Fax: +44 1865 272595
  
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] problem with corrgram function

2011-06-28 Thread Eik Vettorazzi
yes it is.
and a correlation of 0 isn't exactly white (#FF) either.
have a look at the panel.pie function.

the crucial part is

ncol - 14
pal - col.corrgram(ncol)
col.ind - round(ncol * (corr + 1)/2)

so an correlation near -1 maps to an index 0, which isn't a proper index
in R.
Alter these lines to
 ncol - 15 #so 0 becomes #FF
 pal - col.corrgram(ncol)
 col.ind - round((ncol-1) * (corr + 1)/2)+1

hth.


Am 28.06.2011 13:11, schrieb Niels Janssen:
 Dear list,
 
 I have a problem with the corrgram function. It does not seem to
 color large negative correlations, while the same correlation, if
 positive, provides no problems. Is this a bug?
 
 require(corrgram)   
 a = seq(1,100)
 b = -jitter(seq(1,100), 80)
 cor(a,b) # r about -.96
 c=as.data.frame(cbind(a,b))
 corrgram(c, order=NULL, lower.panel=panel.pie,upper.panel=NULL,
 text.panel=panel.txt) # no color
 
 c$b = -1*c$b # flip direction of correlation
 cor(c$a, c$b) # r now about +.96
 corrgram(c, order=NULL, lower.panel=panel.pie,upper.panel=NULL,
 text.panel=panel.txt) #no problem with color.
 
 Thanks!
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

-- 
Eik Vettorazzi

Department of Medical Biometry and Epidemiology
University Medical Center Hamburg-Eppendorf

Martinistr. 52
20246 Hamburg

T ++49/40/7410-58243
F ++49/40/7410-57790

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


Re: [R] coxph() - unexpected result using Crawley's seedlings data (The R Book)

2011-06-28 Thread Jacob Brogren
All,

I rerun once again and managed to reproduce the results from the text book. 
Made no changes to the code. Could it be some problem with convergence?

Anyhow, now it works!

Cheers

Jacob

ps. I find The R Book very useful ds.

28 jun 2011 kl. 15.48 skrev Robert A LaBudde:

 Did you create the 'status' variable the way indicated on p. 797?
 
 Frequently with Surv() it pays to use syntax such as Surv(death, status==1) 
 to make a clear logical statement of what is an event (status==1) vs. 
 censored.
 
 PS. Next time include head(seedlings) and str(seedlings) to make clear what 
 you are using as data.
 
 
 At 06:51 AM 6/28/2011, Jacob Brogren wrote:
 Hi,
 
 I ran the example on pp. 799-800 from Machael Crawley's The R Book using 
 package survival v. 2.36-5, R 2.13.0 and RStudio 0.94.83. The model is a 
 Cox's Proportional Hazards model. The result was quite different compared to 
 the R Book. I have compared my code to the code in the book but can not find 
 any differences in the function call. My results are attached as well as a 
 link to the results presented in the book (link to Google Books).
 
 When running the examples on pp. 797-799 I can't detect any differences in 
 results so I don't think there are errors in the data set or in the creation 
 of the status variable.
 
 -
 Original from the R Book:
 http://books.google.com/books?id=8D4HVx0apZQClpg=PA799ots=rQgd_8ofeSdq=r%20coxph%20crawleypg=PA799#v=onepageqf=false
 
 -
 My result:
  summary(model1)
 Call:
 coxph(formula = Surv(death, status) ~ strata(cohort) * gapsize,
data = seedlings)
 
  n= 60, number of events= 60
 
coef exp(coef)  se(coef)  z 
 Pr(|z|)
 gapsize-0.001893  0.998109  0.593372 -0.003  
   0.997
 gapsize:strata(cohort)cohort=September  0.717407  2.049112  0.860807  0.833  
   0.405
 
   exp(coef) exp(-coef) lower .95 upper 
 .95
 gapsize   0.9981  1.002 0.3120 3.193
 gapsize:strata(cohort)cohort=September2.0491  0.488 0.379211.074
 
 Rsquare= 0.022   (max possible= 0.993 )
 Likelihood ratio test= 1.35  on 2 df,   p=0.5097
 Wald test= 1.32  on 2 df,   p=0.5178
 Score (logrank) test = 1.33  on 2 df,   p=0.514
 
 Anyone have an idea why this is occurring?
 
 Kind Regards
 
 Jacob
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 
 Robert A. LaBudde, PhD, PAS, Dpl. ACAFS  e-mail: r...@lcfltd.com
 Least Cost Formulations, Ltd.URL: http://lcfltd.com/
 824 Timberlake Drive Tel: 757-467-0954
 Virginia Beach, VA 23464-3239Fax: 757-467-2947
 
 Vere scire est per causas scire
 
 

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


Re: [R] Running R from windows command prompt

2011-06-28 Thread Arun Kumar Saha
Hi Siddharth, many experts already answered your query, however I
would like to share how I run R in command prompt:

1. open command prompt
2. change working directory: cd C:\\R-2.13.0\bin\i386
(put the entire path here, however many people might find this step
weird, you can have better management setting window's path variable
appropriately)
3. type R.exe

You can use R within command prompt with same efficiency. However most
awkward thing I find in this process is you can never copy-paste any
code. So everything you need to type there manually!

HTH
_

Arun Kumar Saha, FRM
QUANTITATIVE RISK AND HEDGE CONSULTING SPECIALIST
Visit me at: http://in.linkedin.com/in/ArunFRM

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 masked function is a masked function by any other name

2011-06-28 Thread Juan Carlos Borrás
Ops!
Thank-you Duncan for clarifying the 2 vs. 3 colon difference and a
couple of other things.
Working like a charm now.
Cheers,
jcb!

 If you are using ::: (three colons), then you may be looking into the
 unexported functions in log4r.  The only normal way to see unexported
 functions is to use three colons.

 Duncan Murdoch


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


Re: [R] Running R from windows command prompt

2011-06-28 Thread Robert Baer

Subject: Re: [R] Running R from windows command prompt

Hi Siddharth, many experts already answered your query, however I
would like to share how I run R in command prompt:

1. open command prompt
2. change working directory: cd C:\\R-2.13.0\bin\i386
(put the entire path here, however many people might find this step
weird, you can have better management setting window's path variable
appropriately)
3. type R.exe

You can use R within command prompt with same efficiency. However most
awkward thing I find in this process is you can never copy-paste any
code. So everything you need to type there manually!
---
Actually, you can 'cut and paste' from at Windows Cmd prompt.  It is done by 
clicking the
C: icon in the upper left of the command window, choosing edit, and 'copy' 
or 'paste'

as desired.

Rob

HTH
_

Arun Kumar Saha, FRM
QUANTITATIVE RISK AND HEDGE CONSULTING SPECIALIST
Visit me at: http://in.linkedin.com/in/ArunFRM

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] How do I output all the R-squares of an SUR? summary(fitSUR$eq[[1:4]])$r.squared does not work

2011-06-28 Thread StellathePug
Greetings R Users,
I have a system of equations for which I would like to output all the
R-squares. Assume there are four equations in my system, the only way I
found to output all the R-squares is by calling them out one by one as this:

summary(fitSUR$eq[[1]])$r.squared
summary(fitSUR$eq[[2]])$r.squared
summary(fitSUR$eq[[3]])$r.squared
summary(fitSUR$eq[[4]])$r.squared

But isn't there a way of making this automatic, for example, something more
along the way of 

summary(fitSUR$eq[[1:4]])$r.squared ?

The above does not work, by the way.

I have attached below my sample program. Thanks!
Rita


 SAMPLE PROGRAM 

YX-as.data.frame(matrix(rnorm(280),ncol=14, nrow=20))   ##
generate variables
names(YX) -c(paste(Y, 1:4, sep=), paste(X, 1:10, sep=)) ##
assign variables' names

library(systemfit)

## EQUATIONS:
EQ1 - Y1 ~ X1 + X2 + X4 + X7 + X10  ## equation 1 formula
EQ2 - Y2 ~ X2 + X3 + X5 + X8 + X10  ## equation 2 formula
EQ3 - Y3 ~ X5 + X6 + X7 + X9## equation 3 formula
EQ4 - Y4 ~ X1 + X3 + X4 + X6 + X9   ## equation 4 formula

eqSystem -list(form1 = EQ1, form2 = EQ2, form3 = EQ3, form4 = EQ4)

fitSUR - systemfit(eqSystem, method =SUR, data=YX)


## How do I out put all the R-squares of the system without having to type a
single line for each?

summary(fitSUR$eq[[1]])$r.squared
summary(fitSUR$eq[[2]])$r.squared
summary(fitSUR$eq[[3]])$r.squared
summary(fitSUR$eq[[4]])$r.squared

summary(fitSUR$eq[[1]])$adj.r.squared
summary(fitSUR$eq[[2]])$adj.r.squared
summary(fitSUR$eq[[3]])$adj.r.squared
summary(fitSUR$eq[[4]])$adj.r.squared

--
View this message in context: 
http://r.789695.n4.nabble.com/How-do-I-output-all-the-R-squares-of-an-SUR-summary-fitSUR-eq-1-4-r-squared-does-not-work-tp3630601p3630601.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 do I output all the R-squares of an SUR? summary(fitSUR$eq[[1:4]])$r.squared does not work

2011-06-28 Thread Achim Zeileis

On Tue, 28 Jun 2011, StellathePug wrote:


Greetings R Users,
I have a system of equations for which I would like to output all the
R-squares. Assume there are four equations in my system, the only way I
found to output all the R-squares is by calling them out one by one as this:

summary(fitSUR$eq[[1]])$r.squared
summary(fitSUR$eq[[2]])$r.squared
summary(fitSUR$eq[[3]])$r.squared
summary(fitSUR$eq[[4]])$r.squared


You can abbreviate that to:

  sapply(fitSUR$eq, function(x) summary(x)$r.squared)

Instead of calling summary(fitSUR$eq[[1]]), you can also look at 
summary(fitSUR)$eq[[1]] which leads to identical output. Hence you could 
also do


  sapply(summary(fitSUR)$eq, [[, r.squared)

depending on which you find more intuitive.

hth,
Z


But isn't there a way of making this automatic, for example, something more
along the way of

summary(fitSUR$eq[[1:4]])$r.squared ?

The above does not work, by the way.

I have attached below my sample program. Thanks!
Rita


 SAMPLE PROGRAM 

YX-as.data.frame(matrix(rnorm(280),ncol=14, nrow=20))   ##
generate variables
names(YX) -c(paste(Y, 1:4, sep=), paste(X, 1:10, sep=)) ##
assign variables' names

library(systemfit)

## EQUATIONS:
EQ1 - Y1 ~ X1 + X2 + X4 + X7 + X10  ## equation 1 formula
EQ2 - Y2 ~ X2 + X3 + X5 + X8 + X10  ## equation 2 formula
EQ3 - Y3 ~ X5 + X6 + X7 + X9## equation 3 formula
EQ4 - Y4 ~ X1 + X3 + X4 + X6 + X9   ## equation 4 formula

eqSystem -list(form1 = EQ1, form2 = EQ2, form3 = EQ3, form4 = EQ4)

fitSUR - systemfit(eqSystem, method =SUR, data=YX)


## How do I out put all the R-squares of the system without having to type a
single line for each?

summary(fitSUR$eq[[1]])$r.squared
summary(fitSUR$eq[[2]])$r.squared
summary(fitSUR$eq[[3]])$r.squared
summary(fitSUR$eq[[4]])$r.squared

summary(fitSUR$eq[[1]])$adj.r.squared
summary(fitSUR$eq[[2]])$adj.r.squared
summary(fitSUR$eq[[3]])$adj.r.squared
summary(fitSUR$eq[[4]])$adj.r.squared

--
View this message in context: 
http://r.789695.n4.nabble.com/How-do-I-output-all-the-R-squares-of-an-SUR-summary-fitSUR-eq-1-4-r-squared-does-not-work-tp3630601p3630601.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.


[R] extracting data

2011-06-28 Thread Ana Kolar
Let's say I have an original data set which is called A and data extracted from 
this original data set, called B. Based on these A and B data set I would like 
to get data set C which includes all the remaining data from the data set A 
after we exclude data of the data set B.

Any idea how to do this?
[[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] extracting data

2011-06-28 Thread Sarah Goslee
Hi Ana,

On Tue, Jun 28, 2011 at 12:28 PM, Ana Kolar annako...@yahoo.com wrote:
 Let's say I have an original data set which is called A and data extracted 
 from this original data set, called B. Based on these A and B data set I 
 would like to get data set C which includes all the remaining data from the 
 data set A after we exclude data of the data set B.

 Any idea how to do this?

Yes. Several.

But to know which one to suggest, I need to know more about your data.

How about a toy example, so the list members can see your index
variables, etc? Or how you created the subset B, and why you can't
just use the opposite of that procedure?

Sarah
-- 
Sarah Goslee
http://www.functionaldiversity.org

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


Re: [R] extracting data

2011-06-28 Thread Ana Kolar
Hi Sarah,

Thank you for your response. Here is a toy example:


library(MatchIt)
data(lalonde)

A-lalonde
f-treat ~ age + I(age^2) + educ + I(educ^2) + black + hispan +
    married + nodegree + re74 + I(re74^2) + re75 + I(re75^2)
m-nearest
m.out.base - matchit(formula=f, data=A, method=m)

B - match.data(m.out.base)

An - nrow(A)
Bn - nrow(B)

Cn - An - Bn
C - ??





From: Sarah Goslee sarah.gos...@gmail.com
To: Ana Kolar annako...@yahoo.com
Cc: R r-help@r-project.org
Sent: Tuesday, 28 June 2011, 18:44
Subject: Re: [R] extracting data

Hi Ana,

On Tue, Jun 28, 2011 at 12:28 PM, Ana Kolar annako...@yahoo.com wrote:
 Let's say I have an original data set which is called A and data extracted 
 from this original data set, called B. Based on these A and B data set I 
 would like to get data set C which includes all the remaining data from the 
 data set A after we exclude data of the data set B.

 Any idea how to do this?

Yes. Several.

But to know which one to suggest, I need to know more about your data.

How about a toy example, so the list members can see your index
variables, etc? Or how you created the subset B, and why you can't
just use the opposite of that procedure?

Sarah
-- 
Sarah Goslee
http://www.functionaldiversity.org



[[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] cumulative incidence plot vs survival plot

2011-06-28 Thread array chip
Thank you Alan! Now I sort of understand what it means by competing risk! So in 
cuminc() function, the argument fstatus should be coded like: 0=censored, 
1=event of interest, 2=event of competing risk. Then the function will 
calculate 
CI for each of the 2 types of events (event of interest and event of competing 
risk), am I correct?

What about running regular Cox regression for recurrence? any problem there? 
for 
example, need to take into competing risk as well or regular Cox regression is 
still fine?

Thanks!

John



- Original Message 
From: alanm (Alan Mitchell) al...@crab.org
To: array chip arrayprof...@yahoo.com; David Winsemius 
dwinsem...@comcast.net
Cc: r-help@r-project.org
Sent: Tue, June 28, 2011 9:20:22 AM
Subject: RE: [R] cumulative incidence plot vs survival plot

John,

Since death precludes recurrence, censoring deaths would violate the KM
estimator assumption that additional follow-up would eventually lead to
an event.  If your goal is to estimate the probability of recurrence,
then you want CI with deaths as a competing risk.  The cuminc function
in the cmprsk package is a great place to start.  

Gooley has a great paper on the difference between CI and 1-KM (See
Statistics in Medicine, 18, 695-706 (1999)).  

HTH,

Alan Mitchell, MSc
Biostatistician
al...@crab.org

-Original Message-
From: array chip [mailto:arrayprof...@yahoo.com] 
Sent: Monday, June 27, 2011 2:04 PM
To: David Winsemius
Cc: r-help@r-project.org
Subject: Re: [R] cumulative incidence plot vs survival plot

Hi David,

Thanks for responding, and plain text ...(didn't realized I was in rich
text).

The endpoint is disease recurrence, I was producing a regular KM plot of
recurrence-free probability. Then someone recommend using cumulative
incidence is preferred because death was censored in the dataset. I did
a little googling, I found CI was used often in the context of competing
risk. I am totally new to competing risk and trying to understand what
competing risk means and why CI is preferred than KM survival in this
context. If you could share your thoughts helping me to understand,
greatly appreciated.

Searched archive, found people talking about cmprsk package for
estimating and plotting CI. would that be the same as the code you
suggested: plot(time,
cumsum(dead))

Thanks very much!

John





From: David Winsemius dwinsem...@comcast.net

Cc: r-help@r-project.org
Sent: Mon, June 27, 2011 1:45:35 PM
Subject: Re: [R] cumulative incidence plot vs survival plot


On Jun 27, 2011, at 4:31 PM, array chip wrote:

 Hi, I am wondering if anyone can explain to me if cumulative incidence

 (CI) is just 1 minus kaplan-Meier survival?

First tell us what you think CI is defined as. I suspect it is not the
same. The 

KM estimator is cumulative product of (alive-n(dead))/alive so is the
product of 

interval survival probabilities. I doubt that your definition of CI has
a similar denominator.


 Under what circumstance, you should use
 cumulative incidence vs KM survival? If the relationship is just CI =
 1-survival, then what difference it makes to use one vs. the other?
 
 And in R how I can draw a cumulative incidence plot.

plot(time, cumsum(dead)) ...?

 I know I can make a
 Kaplan-Meier survival plot using plot(survfit()), for example:
 
 fit-survfit(Surv(time,status)~group,data=data)
 plot(fit, col=1:2)
 
 How to draw CI plot then?

As above. Specify what you are seeking.

There is a well-defined relationship between S(t) and the cumulative
hazard. 
Maybe you should do a little study of those terms in texts regarding
survival 
analysis.

 Thanks very much!
 
 John
 [[alternative HTML version deleted]]

Isn't it time you learned to post in plain text?

--
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] Saved EPS does not match screen when using bquote(.(i))

2011-06-28 Thread John Kruschke
Many thanks to Uwe Ligges, Peter Ehlers, and Dennis Murphy for suggesting
work-arounds for this bug. Because the suggestions are work-arounds, rather
than actually correcting the bug, I have opted simply to copy and paste the
plotting commands a few times with the subscripts specified as constants
instead of as evaluated variables (and to use the easy and familiar savePlot
command for EPS output). For future cases that have more iterations through
a loop, I may resort to the suggested workarounds. Thanks again.

John K. Kruschke, Professor


On Thu, Jun 23, 2011 at 12:27 PM, Peter Ehlers ehl...@ucalgary.ca wrote:


 I think that there may be a problem with the way
 bquote(), for() and savePlot() play together in
 the OP's example (multiple plots on a windows device;
 bquote using the loop index).

 Here's a version using replayPlot():

 ## show mu with subscripts 4 and 9:
  x11()
  par(mfrow = c(2,1))
  for (i in c(4, 9)) {
plot(0, 0, main = bquote(mu[.(i)]))
  }

 ## now record and replay:
  z - recordPlot()
  replayPlot(z)
 ## both subscripts are now 9


 Simple workarounds are:

 1. As Uwe and Dennis show, for saving to file,
 use the appropriate device.

 2. For recalling plots (e.g. with recording turned on),
 a) use substitute() instead of bquote() or
 b) insert something like
   i - i
 before the plot() call.

 sessionInfo()
 R version 2.13.0 Patched (2011-05-24 r55981)
 Platform: i386-pc-mingw32/i386 (32-bit)

 locale:
 [1] LC_COLLATE=English_Canada.1252  LC_CTYPE=English_Canada.1252
 [3] LC_MONETARY=English_Canada.**1252 LC_NUMERIC=C
 [5] LC_TIME=English_Canada.1252

 attached base packages:
 [1] stats graphics  grDevices utils datasets  methods   base

 Peter Ehlers



 On 2011-06-22 22:14, Dennis Murphy wrote:

 Hi:

 As Uwe suggested...

 pdf('testgraph.pdf')
 layout( matrix( 1:2 , nrow=2 ) )
 for ( i in 1:2 ) {
   plot( 0 , 0 , xlab=bquote(mu[.(i)]) )
  }
 dev.off()

 postscript('testgraph.ps')
 layout( matrix( 1:2 , nrow=2 ) )
 for ( i in 1:2 ) {
   plot( 0 , 0 , xlab=bquote(mu[.(i)]) )
  }
 dev.off()

 png('testgraph.png')
 layout( matrix( 1:2 , nrow=2 ) )
 for ( i in 1:2 ) {
   plot( 0 , 0 , xlab=bquote(mu[.(i)]) )
  }
 dev.off()

 The three graphs look the same (although the PS graph is rotated to
 landscape while the other two are portrait). The main point is that
 mu_1 and mu_2 show up correctly in the two panels in all three graphs
 (at least on my viewers).

 The following thread from last January describes some of the problems
 that certain viewers have with Greek letters, which appear to be
 viewer and platform dependent:

 http://r-project.markmail.org/**search/?q=pdf%20incorrect#**
 query:pdf%20incorrect+page:2+**mid:egmb6utulrxgcznw+state:**resultshttp://r-project.markmail.org/search/?q=pdf%20incorrect#query:pdf%20incorrect+page:2+mid:egmb6utulrxgcznw+state:results

 I'm guessing that I've seen about a half dozen or so similar posts in
 this forum over the past year and a half, so you can check the list
 archives for related problems.

 HTH,
 Dennis

 On Wed, Jun 22, 2011 at 8:07 PM, John 
 Kruschkejohnkruschke@gmail.**comjohnkrusc...@gmail.com
  wrote:

 Here's a fairly minimal-case example in which the saved EPS does not
 match
 the screen. The error comes when using bquote(.(i)) instead of bquote(1),
 as
 demonstrated by the two minimally different cases below. Very strange.
 Any
 clues as to why?

 # begin ---

 # Version A. X axis labels have subscripts as constants. EPS is correct.
 windows()
 layout( matrix( 1:2 , nrow=2 ) )
 plot( 0 , 0 , xlab=bquote(mu[1]) )
 plot( 0 , 0 , xlab=bquote(mu[2]) )
 savePlot( file=SavePlotTestA.eps , type=eps ) # Axis labels are
 correct
 in EPS.

 # Version B. X axis labels have subscripts as variable index. EPS is
 wrong!
 windows()
 layout( matrix( 1:2 , nrow=2 ) )
 for ( i in 1:2 ) {
  plot( 0 , 0 , xlab=bquote(mu[.(i)]) )
 }
 savePlot( file=SavePlotTestB.eps , type=eps ) # X-AXIS OF PLOT 1 IS
 WRONG IN EPS.

 #-- end -

 Thanks!

 John K. Kruschke, Professor
 http://www.indiana.edu/%**7Ekruschke/**DoingBayesianDataAnalysis/http://www.indiana.edu/%7Ekruschke/DoingBayesianDataAnalysis/
 


 2011/6/22 Uwe 
 Liggeslig...@statistik.tu-**dortmund.delig...@statistik.tu-dortmund.de
 



 On 22.06.2011 13:50, John Kruschke wrote:

  The error happens when using the savePlot() command, like this:
 savePlot( file=TestSavePlot.eps , type=eps )
 savePlot( file=TestSavePlot.jpg , type=jpg )


 Well, plot directly into a device, for postscript:

 postscript(estSavePlot.eps, additionalArguments .)
 plot(1:10)
 dev.off()

 Uwe Ligges


  The images in the two saved files are not the same, with the JPG being

 correct but the EPS being wrong.

 When you suggest starting separate devices explicitly, what do you
 mean?
 (I've skimmed through the results of ??device, but can't make sense of
 it.)
 Thank you!

 John K. Kruschke, Professor


 2011/6/22 Uwe 
 

[R] how to print = in plot title

2011-06-28 Thread array chip
Hi, how can I print = (I mean the symbol of just one character) in the main 
title of a plot?

for example:

plot(1:10, main=paste(x =, x))

where variable x is some number generated on the fly.

Thanks

John

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

2011-06-28 Thread Bert Gunter
This is highly system dependent: what character do you intend to use
for this 2 character representation? Hence, you need to follow the
posting guide and give the at a minimum system info. ?sessionInfo

-- Bert

On Tue, Jun 28, 2011 at 10:25 AM, array chip arrayprof...@yahoo.com wrote:
 Hi, how can I print = (I mean the symbol of just one character) in the main
 title of a plot?

 for example:

 plot(1:10, main=paste(x =, x))

 where variable x is some number generated on the fly.

 Thanks

 John

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




-- 
Men by nature long to get on to the ultimate truths, and will often
be impatient with elementary studies or fight shy of them. If it were
possible to reach the ultimate truths without the elementary studies
usually prefixed to them, these would not be preparatory studies but
superfluous diversions.

-- Maimonides (1135-1204)

Bert Gunter
Genentech Nonclinical Biostatistics

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


Re: [R] extracting data

2011-06-28 Thread Peter Ehlers

On 2011-06-28 09:54, Ana Kolar wrote:

Hi Sarah,

Thank you for your response. Here is a toy example:


library(MatchIt)
data(lalonde)

A-lalonde
f-treat ~ age + I(age^2) + educ + I(educ^2) + black + hispan +
 married + nodegree + re74 + I(re74^2) + re75 + I(re75^2)
m-nearest
m.out.base- matchit(formula=f, data=A, method=m)

B- match.data(m.out.base)

An- nrow(A)
Bn- nrow(B)

Cn- An - Bn
C- ??


Can't you just use

 idx - setdiff(rownames(A), rownames(B))
 C - A[idx, ]

Peter Ehlers








From: Sarah Gosleesarah.gos...@gmail.com
To: Ana Kolarannako...@yahoo.com
Cc: Rr-help@r-project.org
Sent: Tuesday, 28 June 2011, 18:44
Subject: Re: [R] extracting data

Hi Ana,

On Tue, Jun 28, 2011 at 12:28 PM, Ana Kolarannako...@yahoo.com  wrote:

Let's say I have an original data set which is called A and data extracted from 
this original data set, called B. Based on these A and B data set I would like 
to get data set C which includes all the remaining data from the data set A 
after we exclude data of the data set B.

Any idea how to do this?


Yes. Several.

But to know which one to suggest, I need to know more about your data.

How about a toy example, so the list members can see your index
variables, etc? Or how you created the subset B, and why you can't
just use the opposite of that procedure?

Sarah
--
Sarah Goslee
http://www.functionaldiversity.org




[[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] cumulative incidence plot vs survival plot

2011-06-28 Thread array chip
Alan,

Let's say that I code censoring as 0, recurrence as 1 for fstat and 
death/competing risk as 2. If a patient did not have recurrence and lost 
follow-up at 2 years in terms of recurrence monitoring, but he also died at 5 
years. How should I code this patient? I think I still code this patient as 0 
(censoring) because lost-of-followup occurred before death, am I correct?

Thanks very much!

John





- Original Message 
From: alanm (Alan Mitchell) al...@crab.org
To: array chip arrayprof...@yahoo.com; David Winsemius 
dwinsem...@comcast.net
Cc: r-help@r-project.org
Sent: Tue, June 28, 2011 9:20:22 AM
Subject: RE: [R] cumulative incidence plot vs survival plot

John,

Since death precludes recurrence, censoring deaths would violate the KM
estimator assumption that additional follow-up would eventually lead to
an event.  If your goal is to estimate the probability of recurrence,
then you want CI with deaths as a competing risk.  The cuminc function
in the cmprsk package is a great place to start.  

Gooley has a great paper on the difference between CI and 1-KM (See
Statistics in Medicine, 18, 695-706 (1999)).  

HTH,

Alan Mitchell, MSc
Biostatistician
al...@crab.org

-Original Message-
From: array chip [mailto:arrayprof...@yahoo.com] 
Sent: Monday, June 27, 2011 2:04 PM
To: David Winsemius
Cc: r-help@r-project.org
Subject: Re: [R] cumulative incidence plot vs survival plot

Hi David,

Thanks for responding, and plain text ...(didn't realized I was in rich
text).

The endpoint is disease recurrence, I was producing a regular KM plot of
recurrence-free probability. Then someone recommend using cumulative
incidence is preferred because death was censored in the dataset. I did
a little googling, I found CI was used often in the context of competing
risk. I am totally new to competing risk and trying to understand what
competing risk means and why CI is preferred than KM survival in this
context. If you could share your thoughts helping me to understand,
greatly appreciated.

Searched archive, found people talking about cmprsk package for
estimating and plotting CI. would that be the same as the code you
suggested: plot(time,
cumsum(dead))

Thanks very much!

John





From: David Winsemius dwinsem...@comcast.net

Cc: r-help@r-project.org
Sent: Mon, June 27, 2011 1:45:35 PM
Subject: Re: [R] cumulative incidence plot vs survival plot


On Jun 27, 2011, at 4:31 PM, array chip wrote:

 Hi, I am wondering if anyone can explain to me if cumulative incidence

 (CI) is just 1 minus kaplan-Meier survival?

First tell us what you think CI is defined as. I suspect it is not the
same. The 

KM estimator is cumulative product of (alive-n(dead))/alive so is the
product of 

interval survival probabilities. I doubt that your definition of CI has
a similar denominator.


 Under what circumstance, you should use
 cumulative incidence vs KM survival? If the relationship is just CI =
 1-survival, then what difference it makes to use one vs. the other?
 
 And in R how I can draw a cumulative incidence plot.

plot(time, cumsum(dead)) ...?

 I know I can make a
 Kaplan-Meier survival plot using plot(survfit()), for example:
 
 fit-survfit(Surv(time,status)~group,data=data)
 plot(fit, col=1:2)
 
 How to draw CI plot then?

As above. Specify what you are seeking.

There is a well-defined relationship between S(t) and the cumulative
hazard. 
Maybe you should do a little study of those terms in texts regarding
survival 
analysis.

 Thanks very much!
 
 John
 [[alternative HTML version deleted]]

Isn't it time you learned to post in plain text?

--
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] how to print = in plot title

2011-06-28 Thread array chip
Thanks Bert. here is info:

R version 2.12.2 (2011-02-25)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United 
States.1252LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C   LC_TIME=English_United States.1252   
 


attached base packages:
[1] splines   stats graphics  grDevices datasets  utils methods   
base 


other attached packages:
[1] rms_3.3-0   Hmisc_3.8-3 cmprsk_2.2-2survival_2.36-5 
rcom_2.2-3.1rscproxy_1.3-1 


loaded via a namespace (and not attached):
[1] cluster_1.13.3  grid_2.12.2 lattice_0.19-17 tools_2.12.2 





- Original Message 
From: Bert Gunter gunter.ber...@gene.com
To: array chip arrayprof...@yahoo.com
Cc: R r-help@r-project.org
Sent: Tue, June 28, 2011 10:32:29 AM
Subject: Re: [R] how to print = in plot title

This is highly system dependent: what character do you intend to use
for this 2 character representation? Hence, you need to follow the
posting guide and give the at a minimum system info. ?sessionInfo

-- Bert

On Tue, Jun 28, 2011 at 10:25 AM, array chip arrayprof...@yahoo.com wrote:
 Hi, how can I print = (I mean the symbol of just one character) in the main
 title of a plot?

 for example:

 plot(1:10, main=paste(x =, x))

 where variable x is some number generated on the fly.

 Thanks

 John

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




-- 
Men by nature long to get on to the ultimate truths, and will often
be impatient with elementary studies or fight shy of them. If it were
possible to reach the ultimate truths without the elementary studies
usually prefixed to them, these would not be preparatory studies but
superfluous diversions.

-- Maimonides (1135-1204)

Bert Gunter
Genentech Nonclinical Biostatistics

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


Re: [R] how to print = in plot title

2011-06-28 Thread Peter Ehlers

On 2011-06-28 10:25, array chip wrote:

Hi, how can I print = (I mean the symbol of just one character) in the main
title of a plot?

for example:

plot(1:10, main=paste(x=, x))

where variable x is some number generated on the fly.


  x - 2.718
  plot(0, 0)
  title(bquote( x %=% .(x) ))

?plotmath

Peter Ehlers



Thanks

John

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


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


Re: [R] how to print = in plot title

2011-06-28 Thread array chip
Thank Peter! How do I make the title in bold font?

John



- Original Message 
From: Peter Ehlers ehl...@ucalgary.ca
To: array chip arrayprof...@yahoo.com
Cc: R r-help@r-project.org
Sent: Tue, June 28, 2011 10:52:07 AM
Subject: Re: [R] how to print = in plot title

On 2011-06-28 10:25, array chip wrote:
 Hi, how can I print = (I mean the symbol of just one character) in the main
 title of a plot?

 for example:

 plot(1:10, main=paste(x=, x))

 where variable x is some number generated on the fly.

   x - 2.718
   plot(0, 0)
   title(bquote( x %=% .(x) ))

?plotmath

Peter Ehlers


 Thanks

 John

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

2011-06-28 Thread Cody Hamilton
Hello Duncan,

testci.R is a test function for the survival package.  I compared the .Rout and 
.Rout.save files by eyeball (I'm on a Windows 7 machine, so I can't use the 
diff function).  The only differences I found were in the file headers.

This is the header from the .Rout file:


R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-pc-mingw32/i386 (32-bit)

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

  Natural language support but running in an English locale

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

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


Here is the header from the .Rout.save file:


R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

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

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

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

As you can see, the header for the .Rout.save file is shorter resulting in a 
difference in length of three lines.  This appears to be nothing serious to 
worry about, but the the testInstalledPackages function will be run by our IT 
team (who has no experience with R) as part of the IQ testing, and I am worried 
that they may be concerned by the 'files differ in numbers of lines:' message.  
Is there anything I can do to avoid the message?

Regards,
   -Cody


--- On Tue, 6/28/11, Duncan Murdoch murdoch.dun...@gmail.com wrote:

 From: Duncan Murdoch murdoch.dun...@gmail.com
 Subject: Re: [R] testInstalledPackages
 To: Cody Hamilton cody.sh...@yahoo.com
 Cc: r-help@r-project.org
 Date: Tuesday, June 28, 2011, 1:01 AM
 On 27/06/2011 5:56 PM, Cody Hamilton
 wrote:
  Dear group,
 
  When running the installation test:
 
  testInstalledPackages(both,outDir='c:/Test')
 
  I got the following message:
 
  Running ‘testci.R’
  comparing ‘testci.Rout’ to ‘testci.Rout.save’
 ...
  files differ in number of lines:
 
  Please note the test does not result in 'OK' as do the
 other tests.  Is this a concern?
 
 Yes, you should look at the two files to determine why the
 length has 
 changed.
 
 Duncan Murdoch
 
 
  Regards,
      -Cody Hamilton
 
  __
  R-help@r-project.org
 mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained,
 reproducible code.
 


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


Re: [R] how to print = in plot title

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 1:52 PM, Peter Ehlers wrote:


On 2011-06-28 10:25, array chip wrote:
Hi, how can I print = (I mean the symbol of just one character)  
in the main

title of a plot?

for example:

plot(1:10, main=paste(x=, x))

where variable x is some number generated on the fly.


 x - 2.718
 plot(0, 0)
 title(bquote( x %=% .(x) ))


I think John wants the mathematical symbol. As was pointed out in a  
question last week, the `=` plotmath symbol needs to be flanked by  
operands. Non-printing operands can be created with the phantom  
function:


title(main=expression(phantom()=phantom()) )

Contrary to Gunters's comment, this is probably going to work on all  
the three major OS platforms. It depends only on whether there is a  
Symbol font mapped to the output device.


?plotmath


Yes. The details are there.



Peter Ehlers




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] how to print = in plot title

2011-06-28 Thread Bert Gunter
For the record, I was wrong -- using plotmath's less than or equal
does NOT require platform info. However, I was unsure if you meant
that it was some kind of an arrow you wanted to render, a clear
misinterpretation on my part.

-- Bert

On Tue, Jun 28, 2011 at 10:57 AM, array chip arrayprof...@yahoo.com wrote:
 Thank Peter! How do I make the title in bold font?

 John



 - Original Message 
 From: Peter Ehlers ehl...@ucalgary.ca
 To: array chip arrayprof...@yahoo.com
 Cc: R r-help@r-project.org
 Sent: Tue, June 28, 2011 10:52:07 AM
 Subject: Re: [R] how to print = in plot title

 On 2011-06-28 10:25, array chip wrote:
 Hi, how can I print = (I mean the symbol of just one character) in the 
 main
 title of a plot?

 for example:

 plot(1:10, main=paste(x=, x))

 where variable x is some number generated on the fly.

   x - 2.718
   plot(0, 0)
   title(bquote( x %=% .(x) ))

 ?plotmath

 Peter Ehlers


 Thanks

 John

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




-- 
Men by nature long to get on to the ultimate truths, and will often
be impatient with elementary studies or fight shy of them. If it were
possible to reach the ultimate truths without the elementary studies
usually prefixed to them, these would not be preparatory studies but
superfluous diversions.

-- Maimonides (1135-1204)

Bert Gunter
Genentech Nonclinical Biostatistics

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


Re: [R] cumulative incidence plot vs survival plot

2011-06-28 Thread alanm (Alan Mitchell)
John,

Since death precludes recurrence, censoring deaths would violate the KM
estimator assumption that additional follow-up would eventually lead to
an event.  If your goal is to estimate the probability of recurrence,
then you want CI with deaths as a competing risk.  The cuminc function
in the cmprsk package is a great place to start.  

Gooley has a great paper on the difference between CI and 1-KM (See
Statistics in Medicine, 18, 695-706 (1999)).  

HTH,

Alan Mitchell, MSc
Biostatistician
al...@crab.org

-Original Message-
From: array chip [mailto:arrayprof...@yahoo.com] 
Sent: Monday, June 27, 2011 2:04 PM
To: David Winsemius
Cc: r-help@r-project.org
Subject: Re: [R] cumulative incidence plot vs survival plot

Hi David,

Thanks for responding, and plain text ...(didn't realized I was in rich
text).

The endpoint is disease recurrence, I was producing a regular KM plot of
recurrence-free probability. Then someone recommend using cumulative
incidence is preferred because death was censored in the dataset. I did
a little googling, I found CI was used often in the context of competing
risk. I am totally new to competing risk and trying to understand what
competing risk means and why CI is preferred than KM survival in this
context. If you could share your thoughts helping me to understand,
greatly appreciated.

Searched archive, found people talking about cmprsk package for
estimating and plotting CI. would that be the same as the code you
suggested: plot(time,
cumsum(dead))

Thanks very much!

John





From: David Winsemius dwinsem...@comcast.net

Cc: r-help@r-project.org
Sent: Mon, June 27, 2011 1:45:35 PM
Subject: Re: [R] cumulative incidence plot vs survival plot


On Jun 27, 2011, at 4:31 PM, array chip wrote:

 Hi, I am wondering if anyone can explain to me if cumulative incidence

 (CI) is just 1 minus kaplan-Meier survival?

First tell us what you think CI is defined as. I suspect it is not the
same. The 

KM estimator is cumulative product of (alive-n(dead))/alive so is the
product of 

interval survival probabilities. I doubt that your definition of CI has
a similar denominator.


 Under what circumstance, you should use
 cumulative incidence vs KM survival? If the relationship is just CI =
 1-survival, then what difference it makes to use one vs. the other?
 
 And in R how I can draw a cumulative incidence plot.

plot(time, cumsum(dead)) ...?

 I know I can make a
 Kaplan-Meier survival plot using plot(survfit()), for example:
 
 fit-survfit(Surv(time,status)~group,data=data)
 plot(fit, col=1:2)
 
 How to draw CI plot then?

As above. Specify what you are seeking.

There is a well-defined relationship between S(t) and the cumulative
hazard. 
Maybe you should do a little study of those terms in texts regarding
survival 
analysis.

 Thanks very much!
 
 John
 [[alternative HTML version deleted]]

Isn't it time you learned to post in plain text?

--
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] Specify ID variable in daisy{cluster}

2011-06-28 Thread adlynch
Hi David -
I wanted to thank you for your response!  This was exactly what I needed.
Sincerely,
Alicia

On Fri, Jun 17, 2011 at 3:33 PM, David L Carlson [via R]
ml-node+3606287-1326677984-245...@n4.nabble.com wrote:
 You need to use hhid as the rownames for housing.cluster rather than
 including it as a variable in the data.frame:

 housing.cluster -data.frame(htypec1, afforcr1, resyrc1, crowdcc1, chprbos1)
 rownames(housing.cluster) - hhid

 Then it will not be included in the cluster analysis but will be used to
 label the dendrogram.

 -Original Message-
 From: [hidden email] [mailto:[hidden email]] On
 Behalf Of adlynch
 Sent: Thursday, June 16, 2011 12:35 PM
 To: [hidden email]
 Subject: [R] Specify ID variable in daisy{cluster}

 Hi All - I am using the daisy function from the cluster library to create a
 dissimilarity matrix.  I'm going to use that matrix to run a cluster
 analysis.  My participants are identified with the variable, hhid.  However,
 when I try to keep hhid in the dataset that I use to create the
 dissimilarity matrix, daisy uses it to create the matrix rather than
 ignoring it as an ID variable.  I need to have the ID variable so I can
 later on identify which cluster each participant was classified as.  Any
 thoughts would be much appreciated!

 housing.cluster -data.frame(hhid, htypec1, afforcr1, resyrc1, crowdcc1,
 chprbos1)
 housingdiss - daisy(housing.cluster, metric=gower)


 --
 View this message in context:
 http://r.789695.n4.nabble.com/Specify-ID-variable-in-daisy-cluster-tp3603136
 p3603136.html
 Sent from the R help mailing list archive at Nabble.com.

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

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


 
 If you reply to this email, your message will be added to the discussion
 below:
 http://r.789695.n4.nabble.com/Specify-ID-variable-in-daisy-cluster-tp3603136p3606287.html
 To unsubscribe from Specify ID variable in daisy{cluster}, click here.



-- 
Alicia Doyle Lynch, Ph.D.
Boston College, Lynch School of Education
Department of Counseling, Developmental, and Educational Psychology
Chestnut Hill, MA 02467
Phone: (617) 552-4534
e-mail: doyl...@bc.edu


--
View this message in context: 
http://r.789695.n4.nabble.com/Specify-ID-variable-in-daisy-cluster-tp3603136p3630875.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] Need help on nnet

2011-06-28 Thread arun
Hi Georg,

I am new to R and I am curious if there is a simple way to do the feature
selection you described:

feature selection is essentially an exhaustive approach which tries
every possible subset of your predictors, trains a network and sees what
the prediction error is. The subset which is best (lowest error) is then
chosen in the end. It normally (as a side-effect) also gives you something
like an importance ranking of the variables when using backward or forward
feature selection. But be careful of interactions between variables. 

Is it an option with nnet or should I use leaps in conjunction with nnet ?

Thanks,
Arun



--
View this message in context: 
http://r.789695.n4.nabble.com/Need-help-on-nnet-tp3081744p3630984.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] renaming multiple columns + interpolating temperature series

2011-06-28 Thread tomtomme
Greetings R Users, 
I´m new to R but at least managed to read in multiple files:

filenames - list.files(path=getwd())  
numfiles - length(filenames) 
for (all_temp in c(1:numfiles)) {
  filenames[all_temp] - paste(filenames[all_temp],sep=)
 
assign(gsub([.]ASC$,temp,filenames[all_temp]),read.delim2(filenames[all_temp],
fileEncoding=ISO-8859-15, skip = 4))
  }

Now I want to change the column names on the fly within the above loop. How?
I only found out for one file:

  colnames(w01_10temp) - c(date, time, temp, na)

I want then to lineary interpolate date, time and temp from the
original 5 to 1 second interval for all the files, like:

old:
date timetemp   na
1   22.05.1116:00:0023.653  NA
2   22.05.1116:00:0523.541  NA
...

new:
date timetemp   na
1   22.05.1116:00:0023.653  NA
2   22.05.1116:00:0123.631  NA
3   22.05.1116:00:0223.609  NA
 ...

I already found out about the zoo package, but could not really figure out
how to use it correctly...
Thanks for any help.

--
View this message in context: 
http://r.789695.n4.nabble.com/renaming-multiple-columns-interpolating-temperature-series-tp3630927p3630927.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] Function unfold package RcmdrPlugin.survival

2011-06-28 Thread Rapha
Dear all,

I am using the function ‘unfold’ from the ‘RcmdrPlugin.survival’ to convert
my time-varying covariates dataset from wide to long. I managed to have it
working for my data.

However, the problem I have is that the observations after an event, won’t
be dropped from the dataset. For example, see the dataframe below: the event
occurs at 1.2 (event.time=1), but the 1.3 to 1.6 will remain in the dataset.
I did not find in the 'unfold' function an option to drop them, but I was
probably not looking well. 

From the Rossi dataset example, I saw that observations following event were
dropped, and I understood this might be because the values of the
time-varying covariates are NA, after the event.
In my case, all time-varying covariates still have values, even after the
event, because they were extracted from a very large number of environmental
raster files; so I guess the function ‘unfold’ sees them as observations
being at risk again.

Do you know/is there any way, to get rid of these remaining observations
(1.3 to 1.6)?


start stop event.time id_cell month rvf_bin lake month_rain month_veg
1.1 010   1 2   1  2.5 2321
1.2 121   1 2   1  2.5 5666
1.3 230   1 2   1  2.5  311
1.4 340   1 2   1  2.5 1221
1.5 450   1 2   1  2.5 1120
1.6 560   1 2   1  2.5 2135
2.1 010   2NA   0  2.5 1222

Many thanks for your advice,

Best regards,

Rapha

--
View this message in context: 
http://r.789695.n4.nabble.com/Function-unfold-package-RcmdrPlugin-survival-tp3630809p3630809.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 do I output all the R-squares of an SUR? summary(fitSUR$eq[[1:4]])$r.squared does not work

2011-06-28 Thread StellathePug
 sapply(fitSUR$eq, function(x) summary(x)$r.squared)
 You can abbreviate that to:
 sapply(summary(fitSUR)$eq, [[, r.squared)

This is fantastic! Thanks so much. 

I had a hunch that it would be something related to the apply family but I
am still not very good at using it. Thank you immensely for the two
examples, they made my day!!!

Rita

--
View this message in context: 
http://r.789695.n4.nabble.com/How-do-I-output-all-the-R-squares-of-an-SUR-summary-fitSUR-eq-1-4-r-squared-does-not-work-tp3630601p3630900.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-Installation on Unix -- Make: Don't know how to make #. Stop.

2011-06-28 Thread Zhou, Hong
Hi, all,

./configure was run successfully on my HP-UX ia64 server with exit=0, but when 
type make at prompt, get this error Make: Don't know how to make #.  Stop.  
Does anyone has any clues about this message? Thank you very much!
#make
Rmath.h is unchanged
`libRblas.sl' is up to date.
/app/R/R-2.13.0/lib/libRblas.sl is unchanged
`libbz2.a' is up to date.
`libpcre.a' is up to date.
`libz.a' is up to date.
../../../src/include/libintl.h is unchanged
../../../include/libintl.h is unchanged
`charsetalias.h' is up to date.
`libintl.a' is up to date.
`libtre.a' is up to date.
`liblzma.a' is up to date.
Make: Don't know how to make #.  Stop.
*** Error exit code 1
Stop.
Hong

[[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] Loop through each subject

2011-06-28 Thread Edward Patzelt
R help -

I am attempting to write a script that has multiple subjects in 1 data file.
 Each subject has multiple rows with columns as variables.  Here is my code,
I am having problem executing it on each unique subject id (dat$Subject).

   getwd()

setwd(/Users/edwardpatzelt/Desktop/Neuroimaging/MERIT/SRRT/merge)


dat - read.table(test2.txt, header = TRUE, na.strings = NA,
stringsAsFactors = FALSE, sep = \t)


for(i in 1:length(dat))

 {

 for (i in 1:)dat[(unique(dat$Subject)),)]

 {

colg - dat[grep(Green, dat$CueProbe),]

   colg - data.frame(colg$SRRTCue.OnsetTime/1000, (colg$SRRTFix2.OnsetTime-
colg$SRRTCue.OnsetTime)/1000, (ifelse((colg$SRRTProbe.ACC == 1 | colg$Probe==
+), 1, 0)))

   colr - dat[grep(Red, dat$CueProbe),]

   colr - data.frame(colr$SRRTCue.OnsetTime/1000, (colr$SRRTFix2.OnsetTime-
colr$SRRTCue.OnsetTime)/1000, (ifelse((colr$SRRTProbe.ACC == 1 | colr$Probe==
+), 1, 0)))

 write.table(colg, file  = paste(dat$Subject[[1]], sep = \t, append =
green.txt), col.names = FALSE, row.names = FALSE)

  write.table(colr, file  = paste(dat$Subject[[1]], sep = \t, append =
red.txt), col.names = FALSE, row.names = FALSE)

 }

  }






-- 
Edward H. Patzelt
Research Assistant – TRiCAM Lab
University of Minnesota – Psychology/Psychiatry
VA Medical Center
Office: S355 Elliot Hall - Twin Cities Campus
Phone: 612-626-0072  Email: patze...@umn.edu

Please consider the environment before printing this email
www.psych.umn.edu/research/tricam

[[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 print = in plot title

2011-06-28 Thread array chip
Thank you David and Bert. 

x-3plot(1:10)
title(bquote( x = .(x) ))

would do what I want. But I also want the title printed in bold font. so I tried

x-3
plot(1:10)
title(bquote(bold(x = .(x)) ))

But this did not print the less than equal to symbol and the number 3 (from 
variable x) in bold. Anyway to solve that?

Thanks again!

John




- Original Message 
From: David Winsemius dwinsem...@comcast.net
To: Peter Ehlers ehl...@ucalgary.ca
Cc: array chip arrayprof...@yahoo.com; R r-help@r-project.org
Sent: Tue, June 28, 2011 11:04:07 AM
Subject: Re: [R] how to print = in plot title


On Jun 28, 2011, at 1:52 PM, Peter Ehlers wrote:

 On 2011-06-28 10:25, array chip wrote:
 Hi, how can I print = (I mean the symbol of just one character) in the 
main
 title of a plot?
 
 for example:
 
 plot(1:10, main=paste(x=, x))
 
 where variable x is some number generated on the fly.
 
  x - 2.718
  plot(0, 0)
  title(bquote( x %=% .(x) ))

I think John wants the mathematical symbol. As was pointed out in a question 
last week, the `=` plotmath symbol needs to be flanked by operands. 
Non-printing operands can be created with the phantom function:

title(main=expression(phantom()=phantom()) )

Contrary to Gunters's comment, this is probably going to work on all the three 
major OS platforms. It depends only on whether there is a Symbol font mapped to 
the output device.
 
 ?plotmath

Yes. The details are there.

 
 Peter Ehlers
 


David Winsemius, MD
West Hartford, CT

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


[R] Axes labels, greek letters and spaces

2011-06-28 Thread Sam Albers
Hello all,

I can't seem to figure how to use a greek character in expression() in
plot() labels without adding a space. So for example below when plotting
this out

x-1:10
plot(x,x^2, xlab=expression(Chlorophyll~italic(a)~mu~g~cm^-2))

the axis label read as  μ g cm^-2 because I have space there with a tilda.

But if I remove the tilda then my units are mug cm^-2.

Can anyone recommend a way that I can modify the axis label to look for like
this: μg cm^-2

Thanks in advance!

Sam

[[alternative HTML version deleted]]

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


Re: [R] Loop through each subject

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 1:44 PM, Edward Patzelt wrote:


R help -

I am attempting to write a script that has multiple subjects in 1  
data file.
Each subject has multiple rows with columns as variables.  Here is  
my code,
I am having problem executing it on each unique subject id (dat 
$Subject).


One problem that I see is that you are calling all of your files the  
same thing (i.e. overwriting earlier results. Why aren't you using the  
loop index in the naming process?


(And aeppnd is a logical argument in write.table.)

?write.table

--
David.


  getwd()

setwd(/Users/edwardpatzelt/Desktop/Neuroimaging/MERIT/SRRT/merge)


dat - read.table(test2.txt, header = TRUE, na.strings = NA,
stringsAsFactors = FALSE, sep = \t)


for(i in 1:length(dat))

{

for (i in 1:)dat[(unique(dat$Subject)),)]

{

   colg - dat[grep(Green, dat$CueProbe),]

  colg - data.frame(colg$SRRTCue.OnsetTime/1000, (colg 
$SRRTFix2.OnsetTime-
colg$SRRTCue.OnsetTime)/1000, (ifelse((colg$SRRTProbe.ACC == 1 | colg 
$Probe==

+), 1, 0)))

  colr - dat[grep(Red, dat$CueProbe),]

  colr - data.frame(colr$SRRTCue.OnsetTime/1000, (colr 
$SRRTFix2.OnsetTime-
colr$SRRTCue.OnsetTime)/1000, (ifelse((colr$SRRTProbe.ACC == 1 | colr 
$Probe==

+), 1, 0)))

write.table(colg, file  = paste(dat$Subject[[1]], sep = \t,  
append =

green.txt), col.names = FALSE, row.names = FALSE)

 write.table(colr, file  = paste(dat$Subject[[1]], sep = \t,  
append =

red.txt), col.names = FALSE, row.names = FALSE)

}

 }






--
Edward H. Patzelt
Research Assistant – TRiCAM Lab
University of Minnesota – Psychology/Psychiatry
VA Medical Center
Office: S355 Elliot Hall - Twin Cities Campus
Phone: 612-626-0072  Email: patze...@umn.edu

Please consider the environment before printing this email
www.psych.umn.edu/research/tricam

[[alternative HTML version deleted]]

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


David Winsemius, MD
West Hartford, CT

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


Re: [R] chull increase number of points

2011-06-28 Thread Ben Bolker
 pete at nevill.uk.net writes:

 I am using the chull function to create a convex hull of a series of about
 20,000 data points. 

  You already posted this statement (not a question).  One more try?  (You
might as well read the posting guide while you're at it -- please refrain
from sending your e-mail in HTML format if you can avoid it.)


  Ben Bolker

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

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 2:19 PM, array chip wrote:


Thank you David and Bert.

x-3plot(1:10)
title(bquote( x = .(x) ))

would do what I want. But I also want the title printed in bold  
font. so I tried


x-3
plot(1:10)
title(bquote(bold(x = .(x)) ))

But this did not print the less than equal to symbol and the  
number 3 (from

variable x) in bold. Anyway to solve that?


Nope. The Symbol font has no bold type face. Bolding of numbers _can_  
be done inside bquote with bold(as.character(.(x)) )


--
David


Thanks again!

John




- Original Message 
From: David Winsemius dwinsem...@comcast.net
To: Peter Ehlers ehl...@ucalgary.ca
Cc: array chip arrayprof...@yahoo.com; R r-help@r-project.org
Sent: Tue, June 28, 2011 11:04:07 AM
Subject: Re: [R] how to print = in plot title


On Jun 28, 2011, at 1:52 PM, Peter Ehlers wrote:


On 2011-06-28 10:25, array chip wrote:
Hi, how can I print = (I mean the symbol of just one character)  
in the

main

title of a plot?

for example:

plot(1:10, main=paste(x=, x))

where variable x is some number generated on the fly.


x - 2.718
plot(0, 0)
title(bquote( x %=% .(x) ))


I think John wants the mathematical symbol. As was pointed out in a  
question

last week, the `=` plotmath symbol needs to be flanked by operands.
Non-printing operands can be created with the phantom function:

title(main=expression(phantom()=phantom()) )

Contrary to Gunters's comment, this is probably going to work on all  
the three
major OS platforms. It depends only on whether there is a Symbol  
font mapped to

the output device.


?plotmath


Yes. The details are there.



Peter Ehlers




David Winsemius, MD
West Hartford, CT


David Winsemius, MD
West Hartford, CT

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


Re: [R] Axes labels, greek letters and spaces

2011-06-28 Thread David Winsemius

The * operator can be used for a non-space separation.

expression(Chlorophyll*italic(a)~mu*g~cm^-2)

--  
David.


On Jun 28, 2011, at 2:21 PM, Sam Albers wrote:


Hello all,

I can't seem to figure how to use a greek character in expression() in
plot() labels without adding a space. So for example below when  
plotting

this out

x-1:10
plot(x,x^2, xlab=expression(Chlorophyll~italic(a)~mu~g~cm^-2))

the axis label read as  μ g cm^-2 because I have space there with a  
tilda.


But if I remove the tilda then my units are mug cm^-2.

Can anyone recommend a way that I can modify the axis label to look  
for like

this: μg cm^-2

Thanks in advance!

Sam.


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] Error in library (nls)

2011-06-28 Thread Ben Bolker
matotope matotope at gmail.com writes:

 I installed the 2.8.0 version and it seems to be working fine now.

  Hard to see from a brief glance at the code why it would not work
in 2.13.0 as well, but glad you got the problem solved.  If you will
need to maintain and extend this code in the future it may be worth
figuring out what breaks between 2.8.0 and 2.13.0 and fixing it,
as it is often difficult to get advice on the list about old versions
of R (as releases are bi-annual, 2.8.0 should be about 2.5 years old,
which is antiquated by R community standards ...)

  Ben Bolker

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


Re: [R] cumulative incidence plot vs survival plot

2011-06-28 Thread alanm (Alan Mitchell)
So in cuminc() function, the argument fstatus should be coded like:
0=censored, 1=event of interest, 2=event of competing risk. Then the
function will calculate CI for each of the 2 types of events (event of
interest and event of competing risk), am I correct?
Correct.

What about running regular Cox regression for recurrence? any problem
there? for example, need to take into competing risk as well or regular
Cox regression is still fine?
This is a little trickier.  I would strongly suggest reading up on this
before doing any analyses.  There are a couple of different ways to do
this, but the crr function in cmprsk will perform a competing risks
regression.


Alan Mitchell, MSc
Ph: (206) 839-1708
al...@crab.org

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


Re: [R] coxph() - unexpected result using Crawley's seedlings data (The R Book)

2011-06-28 Thread Ben Bolker
Jacob Brogren jacob at brogren.nu writes:

 
 All,
 
 I rerun once again and managed to reproduce the results from the text book. 
 Made no changes to the code. Could
 it be some problem with convergence?

  It is possible, but *extremely* unlikely, to get non-deterministic
results from R (i.e. running the same code twice from an identical
state and getting different answers [run from a clean R session to be
absolutely sure]) -- this only happens if there is a deep, C-level bug
in the internals of the code, which is unlikely in a piece of core
functionality like coxph().  It is slightly more likely, but still
unlikely, that running the code changes the state of the R session in
a subtle way that makes it run differently the second time in a row.
By far the most likely situation is that you have made some minor
change in the state (i.e. you redefined some variable) that allows
the code to reproduce the results in the book.  (I know I've done this
many times, even when I was initially fairly certain that I hadn't
changed anything.)

  In defense of Crawley: he's written a book that is very useful to
a lot of R users, even if advanced R users sometimes find it a bit
sloppy in places. I agree that he hasn't contributed much back
to the community (except for helping a large population
of beginning users to learn how to use R, which is non-trivial),
but I don't think he has gone out of his way to claim any kind of
official status (other than calling his book The R Book).

  Ben Bolker

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

2011-06-28 Thread array chip
David,

I tried your suggestion, still not working:

x-3
plot(1:10)
title(bquote(bold(x) = bold(as.character(.(x)) )))


It prints x=as.character(3) as title

Thanks

John



- Original Message 
From: David Winsemius dwinsem...@comcast.net
To: array chip arrayprof...@yahoo.com
Cc: Peter Ehlers ehl...@ucalgary.ca; R r-help@r-project.org
Sent: Tue, June 28, 2011 11:30:23 AM
Subject: Re: [R] how to print = in plot title


On Jun 28, 2011, at 2:19 PM, array chip wrote:

 Thank you David and Bert.
 
 x-3plot(1:10)
 title(bquote( x = .(x) ))
 
 would do what I want. But I also want the title printed in bold font. so I 
tried
 
 x-3
 plot(1:10)
 title(bquote(bold(x = .(x)) ))
 
 But this did not print the less than equal to symbol and the number 3 (from
 variable x) in bold. Anyway to solve that?

Nope. The Symbol font has no bold type face. Bolding of numbers _can_ be done 
inside bquote with bold(as.character(.(x)) )

--David
 
 Thanks again!
 
 John
 
 
 
 
 - Original Message 
 From: David Winsemius dwinsem...@comcast.net
 To: Peter Ehlers ehl...@ucalgary.ca
 Cc: array chip arrayprof...@yahoo.com; R r-help@r-project.org
 Sent: Tue, June 28, 2011 11:04:07 AM
 Subject: Re: [R] how to print = in plot title
 
 
 On Jun 28, 2011, at 1:52 PM, Peter Ehlers wrote:
 
 On 2011-06-28 10:25, array chip wrote:
 Hi, how can I print = (I mean the symbol of just one character) in the
 main
 title of a plot?
 
 for example:
 
 plot(1:10, main=paste(x=, x))
 
 where variable x is some number generated on the fly.
 
 x - 2.718
 plot(0, 0)
 title(bquote( x %=% .(x) ))
 
 I think John wants the mathematical symbol. As was pointed out in a question
 last week, the `=` plotmath symbol needs to be flanked by operands.
 Non-printing operands can be created with the phantom function:
 
 title(main=expression(phantom()=phantom()) )
 
 Contrary to Gunters's comment, this is probably going to work on all the three
 major OS platforms. It depends only on whether there is a Symbol font mapped 
to
 the output device.
 
 ?plotmath
 
 Yes. The details are there.
 
 
 Peter Ehlers
 
 
 
 David Winsemius, MD
 West Hartford, CT

David Winsemius, MD
West Hartford, CT

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


Re: [R] cumulative incidence plot vs survival plot

2011-06-28 Thread array chip
Thank you Alan again!  Hope you could also share your thought on my another 
email about the coding of  censoring before death..

Thanks again!

John



- Original Message 
From: alanm (Alan Mitchell) al...@crab.org
To: array chip arrayprof...@yahoo.com; David Winsemius 
dwinsem...@comcast.net
Cc: r-help@r-project.org
Sent: Tue, June 28, 2011 11:40:52 AM
Subject: RE: [R] cumulative incidence plot vs survival plot

So in cuminc() function, the argument fstatus should be coded like:
0=censored, 1=event of interest, 2=event of competing risk. Then the
function will calculate CI for each of the 2 types of events (event of
interest and event of competing risk), am I correct?
Correct.

What about running regular Cox regression for recurrence? any problem
there? for example, need to take into competing risk as well or regular
Cox regression is still fine?
This is a little trickier.  I would strongly suggest reading up on this
before doing any analyses.  There are a couple of different ways to do
this, but the crr function in cmprsk will perform a competing risks
regression.


Alan Mitchell, MSc
Ph: (206) 839-1708
al...@crab.org

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


Re: [R] how to print = in plot title

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 2:42 PM, array chip wrote:


David,

I tried your suggestion, still not working:

x-3
plot(1:10)
title(bquote(bold(x) = bold(as.character(.(x)) )))


It prints x=as.character(3) as title


Kids these days! No innovative spirit. Why in my day... they would be  
jumping with the obvious fix to an elder's error.


x-3
plot(1:10)
title( bquote(bold(x) = bold(.(as.character(x)) )))

--
David.


Thanks

John



- Original Message 
From: David Winsemius dwinsem...@comcast.net
To: array chip arrayprof...@yahoo.com
Cc: Peter Ehlers ehl...@ucalgary.ca; R r-help@r-project.org
Sent: Tue, June 28, 2011 11:30:23 AM
Subject: Re: [R] how to print = in plot title


On Jun 28, 2011, at 2:19 PM, array chip wrote:


Thank you David and Bert.

x-3plot(1:10)
title(bquote( x = .(x) ))

would do what I want. But I also want the title printed in bold  
font. so I

tried

x-3
plot(1:10)
title(bquote(bold(x = .(x)) ))

But this did not print the less than equal to symbol and the  
number 3 (from

variable x) in bold. Anyway to solve that?


Nope. The Symbol font has no bold type face. Bolding of numbers  
_can_ be done

inside bquote with bold(as.character(.(x)) )

--David


Thanks again!

John




- Original Message 
From: David Winsemius dwinsem...@comcast.net
To: Peter Ehlers ehl...@ucalgary.ca
Cc: array chip arrayprof...@yahoo.com; R r-help@r-project.org
Sent: Tue, June 28, 2011 11:04:07 AM
Subject: Re: [R] how to print = in plot title


On Jun 28, 2011, at 1:52 PM, Peter Ehlers wrote:


On 2011-06-28 10:25, array chip wrote:
Hi, how can I print = (I mean the symbol of just one  
character) in the

main

title of a plot?

for example:

plot(1:10, main=paste(x=, x))

where variable x is some number generated on the fly.


x - 2.718
plot(0, 0)
title(bquote( x %=% .(x) ))


I think John wants the mathematical symbol. As was pointed out in a  
question

last week, the `=` plotmath symbol needs to be flanked by operands.
Non-printing operands can be created with the phantom function:

title(main=expression(phantom()=phantom()) )

Contrary to Gunters's comment, this is probably going to work on  
all the three
major OS platforms. It depends only on whether there is a Symbol  
font mapped

to

the output device.


?plotmath


Yes. The details are there.



Peter Ehlers




David Winsemius, MD
West Hartford, CT


David Winsemius, MD
West Hartford, CT


David Winsemius, MD
West Hartford, CT

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


Re: [R] Loop through each subject

2011-06-28 Thread Edward Patzelt
Huh, not sure I understand your feedback.  There is only 1 file that is 1
large dataframe.  I want to execute the commands on each subject in the
dataframe and ouput their respective files with the subject # and file type
appended (10_green.txt).

How would I do this loop index?
- Show quoted text -

On Tue, Jun 28, 2011 at 1:27 PM, David Winsemius dwinsem...@comcast.netwrote:


 On Jun 28, 2011, at 1:44 PM, Edward Patzelt wrote:

  R help -

 I am attempting to write a script that has multiple subjects in 1 data
 file.
 Each subject has multiple rows with columns as variables.  Here is my
 code,
 I am having problem executing it on each unique subject id (dat$Subject).


 One problem that I see is that you are calling all of your files the same
 thing (i.e. overwriting earlier results. Why aren't you using the loop index
 in the naming process?

 (And aeppnd is a logical argument in write.table.)

 ?write.table

 --
 David.


  getwd()

 setwd(/Users/edwardpatzelt/**Desktop/Neuroimaging/MERIT/**SRRT/merge)


 dat - read.table(test2.txt, header = TRUE, na.strings = NA,
 stringsAsFactors = FALSE, sep = \t)


 for(i in 1:length(dat))

 {

 for (i in 1:)dat[(unique(dat$Subject)),)**]

 {

   colg - dat[grep(Green, dat$CueProbe),]

  colg - data.frame(colg$SRRTCue.**OnsetTime/1000,
 (colg$SRRTFix2.OnsetTime-
 colg$SRRTCue.OnsetTime)/1000, (ifelse((colg$SRRTProbe.ACC == 1 |
 colg$Probe==
 +), 1, 0)))

  colr - dat[grep(Red, dat$CueProbe),]

  colr - data.frame(colr$SRRTCue.**OnsetTime/1000,
 (colr$SRRTFix2.OnsetTime-
 colr$SRRTCue.OnsetTime)/1000, (ifelse((colr$SRRTProbe.ACC == 1 |
 colr$Probe==
 +), 1, 0)))

write.table(colg, file  = paste(dat$Subject[[1]], sep = \t, append =
 green.txt), col.names = FALSE, row.names = FALSE)

  write.table(colr, file  = paste(dat$Subject[[1]], sep = \t, append =
 red.txt), col.names = FALSE, row.names = FALSE)

}

  }






 --
 Edward H. Patzelt
 Research Assistant – TRiCAM Lab
 University of Minnesota – Psychology/Psychiatry
 VA Medical Center
 Office: S355 Elliot Hall - Twin Cities Campus
 Phone: 612-626-0072  Email: patze...@umn.edu

 Please consider the environment before printing this email
 www.psych.umn.edu/research/**tricamhttp://www.psych.umn.edu/research/tricam

[[alternative HTML version deleted]]

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


 David Winsemius, MD
 West Hartford, CT




-- 
Edward H. Patzelt
Research Assistant – TRiCAM Lab
University of Minnesota – Psychology/Psychiatry
VA Medical Center
Office: S355 Elliot Hall - Twin Cities Campus
Phone: 612-626-0072  Email: patze...@umn.edu

Please consider the environment before printing this email
www.psych.umn.edu/research/tricam

[[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 print = in plot title

2011-06-28 Thread array chip
Thank you David for the obvious fix and sorry for my lacking innovative spirit, 
:-)

Really, I know I haven't grasped the essence of plotting these math symbols in 
R.

John



- Original Message 
From: David Winsemius dwinsem...@comcast.net
To: array chip arrayprof...@yahoo.com
Cc: Peter Ehlers ehl...@ucalgary.ca; R r-help@r-project.org
Sent: Tue, June 28, 2011 11:50:02 AM
Subject: Re: [R] how to print = in plot title


On Jun 28, 2011, at 2:42 PM, array chip wrote:

 David,

 I tried your suggestion, still not working:

 x-3
 plot(1:10)
 title(bquote(bold(x) = bold(as.character(.(x)) )))


 It prints x=as.character(3) as title

Kids these days! No innovative spirit. Why in my day... they would be  
jumping with the obvious fix to an elder's error.

x-3
plot(1:10)
title( bquote(bold(x) = bold(.(as.character(x)) )))

-- 
David.

 Thanks

 John



 - Original Message 
 From: David Winsemius dwinsem...@comcast.net
 To: array chip arrayprof...@yahoo.com
 Cc: Peter Ehlers ehl...@ucalgary.ca; R r-help@r-project.org
 Sent: Tue, June 28, 2011 11:30:23 AM
 Subject: Re: [R] how to print = in plot title


 On Jun 28, 2011, at 2:19 PM, array chip wrote:

 Thank you David and Bert.

 x-3plot(1:10)
 title(bquote( x = .(x) ))

 would do what I want. But I also want the title printed in bold  
 font. so I
 tried

 x-3
 plot(1:10)
 title(bquote(bold(x = .(x)) ))

 But this did not print the less than equal to symbol and the  
 number 3 (from
 variable x) in bold. Anyway to solve that?

 Nope. The Symbol font has no bold type face. Bolding of numbers  
 _can_ be done
 inside bquote with bold(as.character(.(x)) )

 --David

 Thanks again!

 John




 - Original Message 
 From: David Winsemius dwinsem...@comcast.net
 To: Peter Ehlers ehl...@ucalgary.ca
 Cc: array chip arrayprof...@yahoo.com; R r-help@r-project.org
 Sent: Tue, June 28, 2011 11:04:07 AM
 Subject: Re: [R] how to print = in plot title


 On Jun 28, 2011, at 1:52 PM, Peter Ehlers wrote:

 On 2011-06-28 10:25, array chip wrote:
 Hi, how can I print = (I mean the symbol of just one  
 character) in the
 main
 title of a plot?

 for example:

 plot(1:10, main=paste(x=, x))

 where variable x is some number generated on the fly.

 x - 2.718
 plot(0, 0)
 title(bquote( x %=% .(x) ))

 I think John wants the mathematical symbol. As was pointed out in a  
 question
 last week, the `=` plotmath symbol needs to be flanked by operands.
 Non-printing operands can be created with the phantom function:

 title(main=expression(phantom()=phantom()) )

 Contrary to Gunters's comment, this is probably going to work on  
 all the three
 major OS platforms. It depends only on whether there is a Symbol  
 font mapped
 to
 the output device.

 ?plotmath

 Yes. The details are there.


 Peter Ehlers



 David Winsemius, MD
 West Hartford, CT

 David Winsemius, MD
 West Hartford, CT

David Winsemius, MD
West Hartford, CT

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


Re: [R] means and error bars on xyplot for binary data

2011-06-28 Thread Peter Ehlers

On 2011-06-28 00:12, Louis Plough wrote:

Hi,
I have binary (0,1) data for a trait as my response variable, and
a dependent variable, genotype, with three classes (AA, AB, BB).

I would like to plot this data so that across the three genoytpes, even
though the points are all either 0 or 1, i want them to stack up or be seen
using 'jitter'.  So far I have been able to do this using xyplot {lattice}
(code below) but could not get the points to jitter or stack up on boxplot
or bwplot {lattice}.

  I would also like to add to the xyplot object, the mean of the values at
each of these classes  which will vary depending on how many data points are
at 0 and 1 for a given genotype.  I would also like to put error (i.e.
standard error) bars around these estimates.

I have tried using the points() function to put the mean at each of the
genotype classes, but the point ends up off the figure. Any ideas how to get
this going?

here is some example code.


gtype-c(AA,AB,BB)

x-sample(gtype,20,replace=TRUE)
y-sample(c(0,1),20,replace=TRUE)
bin.data-data.frame(x,y)
xyplot(y~x, jitter.y=TRUE, jitter.x=TRUE,factor=.6, data=bin.data)

Then If I wanted to add the means to the plot, I would do this, which will
print the mean points on a box plot, but not an xyplot:

means1- tapply(bin.data$y,bin.data$x,mean)
points(means1,col=red,pch=18)

Is there a way to get the means, and even error bars on the same xyplot?


I don't know why you want to do what you're doing, but it's
pretty easy with lattice (or with ggplot2). But do use
stripplot instead of xyplot. Define an appropriate panel
function to add the mean points and error bars:

  library(lattice)
  library(Hmisc) ## for error bars, using panel.xYplot

  mn - with(bin.data, tapply(y, x, mean))
  M - Cbind(mn, lower = mn - .1, upper = mn + .1) ## NB: capital 'C'

  stripplot(jitter(y, factor =.6) ~ x,
data = bin.data,
ylab = ,
panel = function(...) {
  panel.stripplot(..., jitter=TRUE, factor=1.2, pch=1, cex=2)
  panel.xYplot(seq_along(mn), M, pch=16, cex=2, col=red)
}
  )

Peter Ehlers



Thanks for any help in advance,
Louis

[[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] how to print = in plot title

2011-06-28 Thread David Winsemius


On Jun 28, 2011, at 2:52 PM, array chip wrote:

Thank you David for the obvious fix and sorry for my lacking  
innovative spirit,

:-)

Really, I know I haven't grasped the essence of plotting these math  
symbols in

R.


A lot of people (including me) have had trouble understanding to to  
write proper plotmath code. Here's my 10 cent version of an  
intro...things that were not clear to me on first, second, or third  
reading of the help(plotmath) page. (I am not saying these things are  
not in the page, just that they were not clear to me after multiple  
efforts at digestion. Searching r-help archives for worked examples  
was the most productive learning strategy.)


The stuff on left side of the table in the plotmath help page  
represent the special words and symbols and they need to be  
separated with the proper operators, since spaces are ignored (as  
actually happens in R parsing as well but we are generally shielded  
form that fact). Words other than these specials can be entered  
without quotes since expressions are not evaluated.  The list of greek  
letters is not complete. Not all of the script greeks work.  The *  
and ~ operators are the usual way to separate plotmath sub- 
expressions with either no-space or space displayed respectively.  
Commas are also special in separating individual expression values, so  
don't use them unless you want a multiple expression value (or need  
them inside a plotmath 'paste' or 'list'). You generally do not want  
multiple values of expression vector when constructing 'main' or 'sub'  
arguments, but may need them when constructing labels for plot axes.   
Only use quotes when they are really needed. You need to construct the  
arguments so that infix operations have flanking operands and notice  
that the plotmath paste function has no sep argument.





John



- Original Message 
From: David Winsemius dwinsem...@comcast.net
To: array chip arrayprof...@yahoo.com
Cc: Peter Ehlers ehl...@ucalgary.ca; R r-help@r-project.org
Sent: Tue, June 28, 2011 11:50:02 AM
Subject: Re: [R] how to print = in plot title


On Jun 28, 2011, at 2:42 PM, array chip wrote:


David,

I tried your suggestion, still not working:

x-3
plot(1:10)
title(bquote(bold(x) = bold(as.character(.(x)) )))


It prints x=as.character(3) as title


Kids these days! No innovative spirit. Why in my day... they would be
jumping with the obvious fix to an elder's error.

x-3
plot(1:10)
title( bquote(bold(x) = bold(.(as.character(x)) )))

--
David.


Thanks

John



- Original Message 
From: David Winsemius dwinsem...@comcast.net
To: array chip arrayprof...@yahoo.com
Cc: Peter Ehlers ehl...@ucalgary.ca; R r-help@r-project.org
Sent: Tue, June 28, 2011 11:30:23 AM
Subject: Re: [R] how to print = in plot title


On Jun 28, 2011, at 2:19 PM, array chip wrote:


Thank you David and Bert.

x-3plot(1:10)
title(bquote( x = .(x) ))

would do what I want. But I also want the title printed in bold
font. so I
tried

x-3
plot(1:10)
title(bquote(bold(x = .(x)) ))

But this did not print the less than equal to symbol and the
number 3 (from
variable x) in bold. Anyway to solve that?


Nope. The Symbol font has no bold type face. Bolding of numbers
_can_ be done
inside bquote with bold(as.character(.(x)) )

--David


Thanks again!

John




- Original Message 
From: David Winsemius dwinsem...@comcast.net
To: Peter Ehlers ehl...@ucalgary.ca
Cc: array chip arrayprof...@yahoo.com; R r-help@r-project.org
Sent: Tue, June 28, 2011 11:04:07 AM
Subject: Re: [R] how to print = in plot title


On Jun 28, 2011, at 1:52 PM, Peter Ehlers wrote:


On 2011-06-28 10:25, array chip wrote:

Hi, how can I print = (I mean the symbol of just one
character) in the

main

title of a plot?

for example:

plot(1:10, main=paste(x=, x))

where variable x is some number generated on the fly.


x - 2.718
plot(0, 0)
title(bquote( x %=% .(x) ))


I think John wants the mathematical symbol. As was pointed out in a
question
last week, the `=` plotmath symbol needs to be flanked by operands.
Non-printing operands can be created with the phantom function:

title(main=expression(phantom()=phantom()) )

Contrary to Gunters's comment, this is probably going to work on
all the three
major OS platforms. It depends only on whether there is a Symbol
font mapped

to

the output device.


?plotmath


Yes. The details are there.



Peter Ehlers




David Winsemius, MD
West Hartford, CT


David Winsemius, MD
West Hartford, CT


David Winsemius, MD
West Hartford, CT


David Winsemius, MD
West Hartford, CT

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


  1   2   >