Re: [R] Mixed Effects Model on Within-Subjects Design

2010-05-20 Thread Dave Deriso
Hi Dave,

Thank you for your helpful advice. I will take a look at the multicomp package.

I was wondering where the lme() function outputs the interaction
between condition*difficulty?

Below is the output to the code I had originally sent. Which one of
these is condition*difficulty?

Fixed effects: value ~ condition * diff
 Value Std.Error  DF   t-value p-value
(Intercept)   300109.95  9506.690 688 31.568289  0.
condition2 27717.65  9071.048 688  3.055617  0.0023
condition3-23718.72  9071.048 688 -2.614772  0.0091
diff50 56767.55  9071.048 688  6.258103  0.
diff75120031.80  9071.048 688 13.232408  0.
condition2:diff50 -45481.21 12828.399 688 -3.545354  0.0004
condition3:diff50   7333.37 12828.399 688  0.571651  0.5677
condition2:diff75 -38765.77 12828.399 688 -3.021871  0.0026
condition3:diff75  12919.59 12828.399 688  1.007109  0.3142

Also, why are diff25 and condition1 missing from the output??

Thanks again for your generous help!!!

Best,
Dave Deriso

On Wed, May 19, 2010 at 10:08 PM, David Atkins  wrote:
>
> Dave--
>
> Given that you want all comparisons among all means in your design, you
> won't get that directly in a call to lme (or lmer in lme4 package). Take a
> look at multcomp package and its vignettes, where I think you'll find what
> you're looking for.
>
> cheers, Dave
>
> --
> Dave Atkins, PhD
> Research Associate Professor
> Department of Psychiatry and Behavioral Science
> University of Washington
> datk...@u.washington.edu
>
> Center for the Study of Health and Risk Behaviors (CSHRB)
> 1100 NE 45th Street, Suite 300
> Seattle, WA  98105
> 206-616-3879
> http://depts.washington.edu/cshrb/
> (Mon-Wed)
>
> Center for Healthcare Improvement, for Addictions, Mental Illness,
>  Medically Vulnerable Populations (CHAMMP)
> 325 9th Avenue, 2HH-15
> Box 359911
> Seattle, WA 98104?
> 206-897-4210
> http://www.chammp.org
> (Thurs)
>
> Dear R Experts,
>
> I am attempting to run a mixed effects model on a within-subjects repeated
> measures design, but I am unsure if I am doing it properly. I was hoping
> that someone would be able to offer some guidance.
>
> There are 5 independent variables (subject, condition, difficulty,
> repetition) and 1 dependent measure (value). Condition and difficulty are
> fixed effects and have 3 levels each (1,2,3 and 25,50,75 respectively),
> while subject and repetition are random effects. Three repeated measurements
> (repetitions) were taken for each condition x difficulty pair for each
> subject, making this an entirely within-subject design.
>
>
>
> I would like an output that compares the significance of the 3 levels of
> difficulty for each condition, as well as the overall interaction of
> condition*difficulty. The ideal output would look like this:
>
> condition1:diff25 vs. condition1:diff50  p_value = 
> condition1:diff25 vs. condition1:diff75  p_value = 
> condition1:diff50 vs. condition1:diff75  p_value = 
>
> condition2:diff25 vs. condition1:diff50  p_value = 
> condition2:diff25 vs. condition1:diff75  p_value = 
> condition2:diff50 vs. condition1:diff75  p_value = 
>
> condition3:diff25 vs. condition1:diff50  p_value = 
> condition3:diff25 vs. condition1:diff75  p_value = 
> condition3:diff50 vs. condition1:diff75  p_value = 
>
> condition*diff  p_value = 
>
>
>
> Here is my code:
>
> #get the data
> study.data =read.csv("http://files.davidderiso.com/example_data.csv";,
> header=T)
> attach(study.data)
> subject = factor(subject)
> condition = factor(condition)
> diff = factor(diff)
> rep = factor(rep)
>
> #visualize whats happening
> interaction.plot(diff, condition, value, ylim=c(24,
> 45),ylab="value", xlab="difficulty", trace.label="condition")
>
> #compute the significance
> library(nlme)
> study.lme = lme(value~condition*diff,random=~1|subject/rep)
> summary(study.lme)
>
>
>
> Thank you so much for your generous help!!!
>
> Best,
> Dave Deriso
> UCSD Psychology
>
>        [[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] computer out of memory when using sigpathway

2010-05-20 Thread 山中程
Dear R users,
 
 
 
I am sorry to disturb you! But I really need your help for the usage of 
sigPathwy.
 
 
 
Actually, I want a sliding window analysis for possible chromosome expression 
pattern mining. My research microorganism is a plant pathogen, Gibberella zeae, 
and I first used SAS to divide locus number with 10, 20, 30, or 40 on the 
fungal chromosome according to their location. I really want to see whether 
among the continual 10, 20, 30, or 40 locus has some expression pattern that 
different from rest genes. Because I know sigPathway (R package, pathway 
analysis with microarray data) can do this kind of job. What I use SAS to do is 
to subset locus in arbitrary genes numbers, such as 10, 20, 30, 40, or so on, 
and I hope to use sigPathway to analysis whether these genes chromosome 
location have effect on its gene expression. 
 
When I used sigpathway to analyze my microarray data, it made my computer out 
of memory. I have tried the following R codes in several computer, but it 
always the same, even if it computed more than one day, it can not get any 
results. I also try to use memory.limit(size=NA), it dosen't work too.My 
computer is a cor duo 2G DDR2, I think it is large enough for my job. Would you 
please point out my problem and give me some suggestions? Thank you very much.
 
I attach my microarray data and R codes in the attachment, and I hope you can 
have a look. 
 
 
 
 #the following code is for annotation list initiation.
 
   
 
 setwd("C:/analysis data and codes")
 
   x <- read.table("chr1.txt",header=FALSE,sep="\t")
 
   attach(x)
 
   x$group <- paste(V2,V3,sep="_")
 
   group <- x$group
 
   y <- data.frame(group,V2,V3,V4)
 
   xx <- as.list(group)
 
   xx <- xx[!is.na(xx)]
 
   xx <- unlist(xx)
 
   xxUnique <- unique(xx)
 
   yy <- vector("list",length(xxUnique))
 
   for(i in 1:length(yy))
 
   {
 
MT <- "MT_lab"
 
yy[[i]] <- 
list(src=MT,title=xxUnique[i],probes=as.character(y[group==xxUnique[i],]$V4))
 
}
 
  
 
  
 
#the following code is for sigpathway analysis.
 
   
 
library(sigPathway)
 
YANG <- read.table("All microarray MT_LAB.txt",header=T,sep="\t")
 
attach(YANG)
 
Y <- data.frame(TF134_1_3DAK,TF134_2_3DAK,WT1_3DAK,WT2_3DAK,row.names=locus_no)
 
p <- c("1_trt","1_trt","0_norm","0_norm")
 
statList <- calcTStatFast(Y,p,ngroups=2)
 
hist(statList$pval,breaks=seq(0,1,0.025),xlab="p-value",ylab="Frequency",main="")
 
set.seed(1234)
 
YANG <- 
runSigPathway(yy,20,500,Y,p,nsim=100,weightType="constant",ngroup=2,npath=10,verbose=F,allpathway=F,alwaysUseRandomPerm=F)
 
#write.table(YANG$df.pathways[1:25, ])
 
write.table(YANG$df.pathways[1:25,],quote=F,sep="\t",file="chr1_sig.txt")
 
YANG$list.gPS[[1]] 
 
save.image("chr1_sig")__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Mixed Effects Model on Within-Subjects Design

2010-05-20 Thread Dave Deriso
Hi Dave,

Thank you for your helpful advice. I will take a look at the multicomp package.

I was wondering where the lme() function outputs the interaction
between condition*difficulty?

Below is the output to the code I had originally sent. Which one of
these is condition*difficulty?

Fixed effects: value ~ condition * diff
  Value Std.Error  DF   t-value p-value
(Intercept)   300109.95  9506.690 688 31.568289  0.
condition2 27717.65  9071.048 688  3.055617  0.0023
condition3-23718.72  9071.048 688 -2.614772  0.0091
diff50 56767.55  9071.048 688  6.258103  0.
diff75120031.80  9071.048 688 13.232408  0.
condition2:diff50 -45481.21 12828.399 688 -3.545354  0.0004
condition3:diff50   7333.37 12828.399 688  0.571651  0.5677
condition2:diff75 -38765.77 12828.399 688 -3.021871  0.0026
condition3:diff75  12919.59 12828.399 688  1.007109  0.3142

Also, why are diff25 and condition1 missing from the output??

Thanks again for your generous help!!!

Best,
Dave Deriso


On Wed, May 19, 2010 at 10:08 PM, David Atkins  wrote:
>
> Dave--
>
> Given that you want all comparisons among all means in your design, you won't 
> get that directly in a call to lme (or lmer in lme4 package). Take a look at 
> multcomp package and its vignettes, where I think you'll find what you're 
> looking for.
>
> cheers, Dave
>
> --
> Dave Atkins, PhD
> Research Associate Professor
> Department of Psychiatry and Behavioral Science
> University of Washington
> datk...@u.washington.edu
>
> Center for the Study of Health and Risk Behaviors (CSHRB)
> 1100 NE 45th Street, Suite 300
> Seattle, WA  98105
> 206-616-3879
> http://depts.washington.edu/cshrb/
> (Mon-Wed)
>
> Center for Healthcare Improvement, for Addictions, Mental Illness,
>  Medically Vulnerable Populations (CHAMMP)
> 325 9th Avenue, 2HH-15
> Box 359911
> Seattle, WA 98104?
> 206-897-4210
> http://www.chammp.org
> (Thurs)
>
> Dear R Experts,
>
> I am attempting to run a mixed effects model on a within-subjects repeated
> measures design, but I am unsure if I am doing it properly. I was hoping
> that someone would be able to offer some guidance.
>
> There are 5 independent variables (subject, condition, difficulty,
> repetition) and 1 dependent measure (value). Condition and difficulty are
> fixed effects and have 3 levels each (1,2,3 and 25,50,75 respectively),
> while subject and repetition are random effects. Three repeated measurements
> (repetitions) were taken for each condition x difficulty pair for each
> subject, making this an entirely within-subject design.
>
>
>
> I would like an output that compares the significance of the 3 levels of
> difficulty for each condition, as well as the overall interaction of
> condition*difficulty. The ideal output would look like this:
>
> condition1:diff25 vs. condition1:diff50  p_value = 
> condition1:diff25 vs. condition1:diff75  p_value = 
> condition1:diff50 vs. condition1:diff75  p_value = 
>
> condition2:diff25 vs. condition1:diff50  p_value = 
> condition2:diff25 vs. condition1:diff75  p_value = 
> condition2:diff50 vs. condition1:diff75  p_value = 
>
> condition3:diff25 vs. condition1:diff50  p_value = 
> condition3:diff25 vs. condition1:diff75  p_value = 
> condition3:diff50 vs. condition1:diff75  p_value = 
>
> condition*diff  p_value = 
>
>
>
> Here is my code:
>
> #get the data
> study.data =read.csv("http://files.davidderiso.com/example_data.csv";,
> header=T)
> attach(study.data)
> subject = factor(subject)
> condition = factor(condition)
> diff = factor(diff)
> rep = factor(rep)
>
> #visualize whats happening
> interaction.plot(diff, condition, value, ylim=c(24,
> 45),ylab="value", xlab="difficulty", trace.label="condition")
>
> #compute the significance
> library(nlme)
> study.lme = lme(value~condition*diff,random=~1|subject/rep)
> summary(study.lme)
>
>
>
> Thank you so much for your generous help!!!
>
> Best,
> Dave Deriso
> UCSD Psychology
>
>        [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] save in for loop

2010-05-20 Thread Ivan Calandra
Thanks to all of you for your answers!

Peter's is definitely the easiest :)
for (i in 1:4) {
   temp <- data.frame(a=(i+1):(i+10), b=LETTERS[(i+1):(i+10)])
   filename <- paste("file", i, sep="")
   assign(filename, temp)
   save(list=c(filename), file=paste(filename, ".rda", sep=""))
}

Tao, I don't understand why you have backslashes before "file" and after 
.rda. I guess it's something about regular expression, but I'm still 
very new to it.
eval(parse(text=paste("save(file", i, ", file=\"file", i, ".rda\")", 
sep="")))

Jorge, your solution does not work... I've just copy/pasted your code. 
My second great weakness is with the apply() family. So maybe I have to 
adjust some part of the code to my needs, but I'm unable to do it.
i <- 1:4
sapply(i, function(i) {
   x <- data.frame(a=(i+1):(i+10), b=LETTERS[(i+1):(i+10)])
save(x, file = paste("file", i, ".rda", sep=""))
} )


Anyway, everything's now fine!
Thanks again.
Ivan


Le 5/19/2010 20:29, Jorge Ivan Velez a écrit :
> Hi Ivan,
>
> How about this?
>
> i <- 1:4
> sapply(i,
> function(i){
> x <- data.frame(a=(i+1):(i+10), b=LETTERS[(i+1):(i+10)])
> save(x, file = paste("file", i, ".rda", sep=""))
> }
> )
>
> HTH,
> Jorge

Le 5/19/2010 22:20, Peter Ehlers a écrit :
> On 2010-05-19 12:05, Shi, Tao wrote:
>> Ivan,
>>
>> Try this:
>>
>> eval(parse(text=paste("save(file", i, ", file=\"file", i, 
>> ".RData\")", sep="")))
>>
>> ...Tao
>>
>
> Or just use 'list=' like this:
>
> for (i in 1:4) {
>   temp <- data.frame(a=(i+1):(i+10), b=LETTERS[(i+1):(i+10)])
>   filename <- paste("file", i, sep="")
>   assign(filename, temp)
>   save(list=c(filename), file=paste(filename, ".rda", sep=""))
> }
>
>  -Peter Ehlers
>
>>
>>
>> - Original Message 
>>> From: Ivan Calandra
>>> To: r-help@r-project.org
>>> Sent: Wed, May 19, 2010 7:56:44 AM
>>> Subject: [R] save in for loop
>>>
>>> Dear users,
>>
>> My problem concerns save() within a for loop.
>> Here is my
>>> code:
>>
>> for (i in 1:4) {
>> temp<- data.frame(a=(i+1):(i+10),
>>> b=LETTERS[(i+1):(i+10)])
>> filename<- paste("file", i, sep="")
>>
>>> assign(filename, temp)
>> save(filename, file=paste(filename, ".rda",
>>> sep=""))
>> }
>>
>> As you can see, save() doesn't work as I would like: (1)
>>> the object saved is called "filename" (instead of "file1", "file2", 
>>> etc), and
>>> (2) it of course contains only the name (as character) instead of the
>>> data.frame
>>
>> How can I fix it?
>>
> [snip]
>

-- 
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php


[[alternative HTML version deleted]]

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


[R] writing autocorrelation and partial auto correlation functions to a file

2010-05-20 Thread nuncio m
Dear All,
  I am very new to T.  I need to fit a ARIMA model to my time
series.  So I found the auto correlation functions and partial auto
correlation function in R. Now I want to save these valuse along with the
significance levels to a file.  How to do that?.  I tried some function in R
like write.table but returns an error  "cannot coerce class "acf" into a
data.frame".  The I tried "write"
but returned again an error "Error in cat(list(...), file, sep, fill,
labels, append) :
  argument 1 (type 'list') cannot be handled by 'cat'".
Where am I going wrong

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

[[alternative HTML version deleted]]

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


Re: [R] save in for loop

2010-05-20 Thread Ivan Calandra
Hi Dennis,

What is the problem with using eval(parse(text=...))? Is there a reason 
why not to use it?
Of course, in that case, "list" within save() is much easier and works 
perfectly.

In any case, thanks for your explanations :)
Ivan

Le 5/20/2010 10:26, Dennis Murphy a écrit :
> Hi:
>
> On Thu, May 20, 2010 at 12:43 AM, Ivan Calandra 
> mailto:ivan.calan...@uni-hamburg.de>> 
> wrote:
>
> Thanks to all of you for your answers!
>
> Peter's is definitely the easiest :)
> for (i in 1:4) {
>   temp <- data.frame(a=(i+1):(i+10), b=LETTERS[(i+1):(i+10)])
>   filename <- paste("file", i, sep="")
>   assign(filename, temp)
>   save(list=c(filename), file=paste(filename, ".rda", sep=""))
> }
>
>
>
> Tao, I don't understand why you have backslashes before "file" and after
>
> .rda. I guess it's something about regular expression, but I'm still
> very new to it.
> eval(parse(text=paste("save(file", i, ", file=\"file", i, ".rda\")",
> sep="")))
>
>
> eval(parse(text = ...)) is necessary sometimes, but usually not. This 
> is one of
> the 'not' cases. Re your question, the backslash is used to escape the 
> quote
> within the quote so that it is rendered properly when parsed/evaluated.
>
>
> Jorge, your solution does not work... I've just copy/pasted your code.
> My second great weakness is with the apply() family. So maybe I
> have to
> adjust some part of the code to my needs, but I'm unable to do it.
> i <- 1:4
> sapply(i, function(i) {
>   x <- data.frame(a=(i+1):(i+10), b=LETTERS[(i+1):(i+10)])
> save(x, file = paste("file", i, ".rda", sep=""))
> } )
>
>
> I tried a variation on this solution, and discovered that whenever you 
> load
> a .rda file, the name of the object is the same in all of them. It 
> 'works' in the
> sense that the right object is saved to each of the file*.rda; 
> however, when
> loaded, all the objects have the same name x, so in the end it doesn't 
> really
> work.
>
> Peter figured out that saving the object as a list was the key - that 
> way,
> the name of the object saved is the value of filename (file*) so that 
> when it
> is loaded back in, the object names are distinct.
>
> A more interesting question than it appeared at first...
>
> Dennis
>
>
> Anyway, everything's now fine!
> Thanks again.
> Ivan
>
>
> Le 5/19/2010 20:29, Jorge Ivan Velez a écrit :
> > Hi Ivan,
> >
> > How about this?
> >
> > i <- 1:4
> > sapply(i,
> > function(i){
> > x <- data.frame(a=(i+1):(i+10), b=LETTERS[(i+1):(i+10)])
> > save(x, file = paste("file", i, ".rda", sep=""))
> > }
> > )
> >
> > HTH,
> > Jorge
>
> Le 5/19/2010 22:20, Peter Ehlers a écrit :
> > On 2010-05-19 12:05, Shi, Tao wrote:
> >> Ivan,
> >>
> >> Try this:
> >>
> >> eval(parse(text=paste("save(file", i, ", file=\"file", i,
> >> ".RData\")", sep="")))
> >>
> >> ...Tao
> >>
> >
> > Or just use 'list=' like this:
> >
> > for (i in 1:4) {
> >   temp <- data.frame(a=(i+1):(i+10), b=LETTERS[(i+1):(i+10)])
> >   filename <- paste("file", i, sep="")
> >   assign(filename, temp)
> >   save(list=c(filename), file=paste(filename, ".rda", sep=""))
> > }
> >
> >  -Peter Ehlers
> >
> >>
> >>
> >> - Original Message 
> >>> From: Ivan Calandra >
> >>> To: r-help@r-project.org 
> >>> Sent: Wed, May 19, 2010 7:56:44 AM
> >>> Subject: [R] save in for loop
> >>>
> >>> Dear users,
> >>
> >> My problem concerns save() within a for loop.
> >> Here is my
> >>> code:
> >>
> >> for (i in 1:4) {
> >> temp<- data.frame(a=(i+1):(i+10),
> >>> b=LETTERS[(i+1):(i+10)])
> >> filename<- paste("file", i, sep="")
> >>
> >>> assign(filename, temp)
> >> save(filename, file=paste(filename, ".rda",
> >>> sep=""))
> >> }
> >>
> >> As you can see, save() doesn't work as I would like: (1)
> >>> the object saved is called "filename" (instead of "file1",
> "file2",
> >>> etc), and
> >>> (2) it of course contains only the name (as character) instead
> of the
> >>> data.frame
> >>
> >> How can I fix it?
> >>
> > [snip]
> >
>
> --
> Ivan CALANDRA
> PhD Student
> University of Hamburg
> Biozentrum Grindel und Zoologisches Museum
> Abt. Säugetiere
> Martin-Luther-King-Platz 3
> D-20146 Hamburg, GERMANY
> +49(0)40 42838 6231
> ivan.calan...@uni-hamburg.de 
>
> **
> http://www.for771.uni-bonn.de
> http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php
>
>
>[[alternative HTML version deleted]]
>
>
> __
> R-help@r-project.org  mailing l

[R] Deleting a file on a drive from within R

2010-05-20 Thread Sergey Goriatchev
Hello,

I have an Excel file on a drive and I extract data from it into R session.
Once I have extracted the data, I want to delete that Excel file from the drive.
Can I do that from within R, please?

Thank you for help!

Regards,
Sergey

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


Re: [R] Deleting a file on a drive from within R

2010-05-20 Thread Tal Galili
Check out:
?unlink





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




On Thu, May 20, 2010 at 11:46 AM, Sergey Goriatchev wrote:

> Hello,
>
> I have an Excel file on a drive and I extract data from it into R session.
> Once I have extracted the data, I want to delete that Excel file from the
> drive.
> Can I do that from within R, please?
>
> Thank you for help!
>
> Regards,
> Sergey
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Deleting a file on a drive from within R

2010-05-20 Thread Dave Deriso
file_to_delete = file.choose()
file.remove(file_to_delete)

Best,
Dave Deriso
UCSD Psychology

On Thu, May 20, 2010 at 1:46 AM, Sergey Goriatchev  wrote:
> Hello,
>
> I have an Excel file on a drive and I extract data from it into R session.
> Once I have extracted the data, I want to delete that Excel file from the 
> drive.
> Can I do that from within R, please?
>
> Thank you for help!
>
> Regards,
> Sergey
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] Deleting a file on a drive from within R

2010-05-20 Thread Sergey Goriatchev
Thank you, Dave!

Regards,
Sergey

On Thu, May 20, 2010 at 10:55, Dave Deriso  wrote:
> file_to_delete = file.choose()
> file.remove(file_to_delete)
>
> Best,
> Dave Deriso
> UCSD Psychology
>
> On Thu, May 20, 2010 at 1:46 AM, Sergey Goriatchev  wrote:
>> Hello,
>>
>> I have an Excel file on a drive and I extract data from it into R session.
>> Once I have extracted the data, I want to delete that Excel file from the 
>> drive.
>> Can I do that from within R, please?
>>
>> Thank you for help!
>>
>> Regards,
>> Sergey
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>



-- 
Simplicity is the last step of art./Bruce Lee
The more you know, the more you know you don't know. /Myself

I'm not young enough to know everything. /Oscar Wilde
Experience is one thing you can't get for nothing. /Oscar Wilde
When you are finished changing, you're finished. /Benjamin Franklin
Luck is where preparation meets opportunity. /George Patten

Kniven skärpes bara mot stenen.

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


Re: [R] Mixed Effects Model on Within-Subjects Design

2010-05-20 Thread ONKELINX, Thierry
Dear Dave,

I think you want this model.

lme(value~condition:diff - 1,random=~1|subject)

Note that I removed the replicate ID from the model. Include it in the model 
makes only sense if you can expect a similar replication effects the 
first/second/thirth time that a subject performs your test.

HTH,

Thierry


ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium

Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium

tel. + 32 54/436 185
thierry.onkel...@inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than 
asking him to perform a post-mortem examination: he may be able to say what the 
experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure 
that a reasonable answer can be extracted from a given body of data.
~ John Tukey
  

> -Oorspronkelijk bericht-
> Van: r-help-boun...@r-project.org 
> [mailto:r-help-boun...@r-project.org] Namens Dave Deriso
> Verzonden: donderdag 20 mei 2010 9:27
> Aan: David Atkins
> CC: r-help@r-project.org
> Onderwerp: Re: [R] Mixed Effects Model on Within-Subjects Design
> 
> Hi Dave,
> 
> Thank you for your helpful advice. I will take a look at the 
> multicomp package.
> 
> I was wondering where the lme() function outputs the
> interaction between condition*difficulty?
> 
> Below is the output to the code I had originally sent. Which 
> one of these is condition*difficulty?
> 
> Fixed effects: value ~ condition * diff
>   Value Std.Error  DF   t-value p-value
> (Intercept)   300109.95  9506.690 688 31.568289  0.
> condition2 27717.65  9071.048 688  3.055617  0.0023
> condition3-23718.72  9071.048 688 -2.614772  0.0091
> diff50 56767.55  9071.048 688  6.258103  0.
> diff75120031.80  9071.048 688 13.232408  0.
> condition2:diff50 -45481.21 12828.399 688 -3.545354  0.0004
> condition3:diff50   7333.37 12828.399 688  0.571651  0.5677
> condition2:diff75 -38765.77 12828.399 688 -3.021871  0.0026
> condition3:diff75  12919.59 12828.399 688  1.007109  0.3142
> 
> Also, why are diff25 and condition1 missing from the output??
> 
> Thanks again for your generous help!!!
> 
> Best,
> Dave Deriso
> 
> 
> On Wed, May 19, 2010 at 10:08 PM, David Atkins 
>  wrote:
> >
> > Dave--
> >
> > Given that you want all comparisons among all means in your 
> design, you won't get that directly in a call to lme (or lmer 
> in lme4 package). Take a look at multcomp package and its 
> vignettes, where I think you'll find what you're looking for.
> >
> > cheers, Dave
> >
> > --
> > Dave Atkins, PhD
> > Research Associate Professor
> > Department of Psychiatry and Behavioral Science University of 
> > Washington datk...@u.washington.edu
> >
> > Center for the Study of Health and Risk Behaviors (CSHRB) 
> 1100 NE 45th 
> > Street, Suite 300 Seattle, WA  98105
> > 206-616-3879
> > http://depts.washington.edu/cshrb/
> > (Mon-Wed)
> >
> > Center for Healthcare Improvement, for Addictions, Mental Illness,
> >  Medically Vulnerable Populations (CHAMMP)
> > 325 9th Avenue, 2HH-15
> > Box 359911
> > Seattle, WA 98104?
> > 206-897-4210
> > http://www.chammp.org
> > (Thurs)
> >
> > Dear R Experts,
> >
> > I am attempting to run a mixed effects model on a within-subjects 
> > repeated measures design, but I am unsure if I am doing it 
> properly. I 
> > was hoping that someone would be able to offer some guidance.
> >
> > There are 5 independent variables (subject, condition, difficulty,
> > repetition) and 1 dependent measure (value). Condition and 
> difficulty 
> > are fixed effects and have 3 levels each (1,2,3 and 25,50,75 
> > respectively), while subject and repetition are random 
> effects. Three 
> > repeated measurements
> > (repetitions) were taken for each condition x difficulty 
> pair for each 
> > subject, making this an entirely within-subject design.
> >
> >
> >
> > I would like an output that compares the significance of 
> the 3 levels 
> > of difficulty for each condition, as well as the overall 
> interaction 
> > of condition*difficulty. The ideal output would look like this:
> >
> > condition1:diff25 vs. condition1:diff50  p_value = 
> > condition1:diff25 vs. condition1:diff75  p_value = 
> > condition1:diff50 vs. condition1:diff75  p_value = 
> >
> > condition2:diff25 vs. condition1:diff50  p_value = 
> > condition2:diff25 vs. condition1:diff75  p_value = 
> > condition2:diff50 vs. condition1:diff75  p_value = 
> >
> > condition3:diff25 vs. condition1:diff50  p_value = 
> > condition3:diff25 vs. condition1:diff75  p_value = 
> > condition3:diff50 vs. condition1:diff75  p_value = 
> >
> >

Re: [R] Deleting a file on a drive from within R

2010-05-20 Thread Sergey Goriatchev
Thanks, Tal

Will check it out.

Regards,
Sergey

On Thu, May 20, 2010 at 10:51, Tal Galili  wrote:
> Check out:
> ?unlink
>
>
>
>
> Contact
> Details:---
> Contact me: tal.gal...@gmail.com |  972-52-7275845
> Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
> www.r-statistics.com (English)
> --
>
>
>
>
> On Thu, May 20, 2010 at 11:46 AM, Sergey Goriatchev 
> wrote:
>>
>> Hello,
>>
>> I have an Excel file on a drive and I extract data from it into R session.
>> Once I have extracted the data, I want to delete that Excel file from the
>> drive.
>> Can I do that from within R, please?
>>
>> Thank you for help!
>>
>> Regards,
>> Sergey
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>



-- 
Simplicity is the last step of art./Bruce Lee
The more you know, the more you know you don't know. /Myself

I'm not young enough to know everything. /Oscar Wilde
Experience is one thing you can't get for nothing. /Oscar Wilde
When you are finished changing, you're finished. /Benjamin Franklin
Luck is where preparation meets opportunity. /George Patten

Kniven skärpes bara mot stenen.

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


[R] Minimization problem

2010-05-20 Thread Fred
Dear R users,

I am trying to minimize two function simultaneously in R,

function(x)

minimize x[1],x[2],x[3]

mean(distribution(x1,x2,x3) ) - observed mean

std(distribution(x1,x2,x3)) - observed std

What I want to achieve is that simulated mean and standard deviation
of distribution related to x1 x2 x3  would be close to observed mean
and observed standard deviation.


is there any function in R can reach this?

Thank you for the help first .

F.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Compressed values on y-axis in effects plot

2010-05-20 Thread Simon Kiss
Dear colleagues, the code below generated the two effects plots that I  
have attached. I hope they are not stripped.


The original two models are as follows:
green_shift_mod=glm(green_shift ~ educ+party_id+educ:party_id,  
family=binomial, data=x)
carbon_tax_mod=glm(carbon_tax ~ educ+party_id+educ:party_id,  
family=binomial, data=x)


Then, I try to plot the effects of party_id by education for both models

It works well for carbon_tax_mod; but for green_shift_mod, effects  
plots the effects of party ID by education in a straight, horizontal  
line, with the values completely compressed.
I've looked through; all the variables included in the two models are  
identical save for the DV. And the DV's in both models are ordered  
factors.

Is any one familiar with this problem in effects plots?
Yours, Simon Kiss

quartz()
jpeg(filename="test.jpeg", type=c("quartz"))
plot(effect("educ:party_id", green_shift_mod, rug=TRUE),  
ylab="Probability of Disagreeing", xlab="Party ID", main="Probability  
of Disagreeing That The Green Shift Would Hurt The Economy")

dev.off()
quartz()
jpeg(filename="test2.jpeg", type=c("quartz"))
plot(effect("educ:party_id", carbon_tax_mod, rug=TRUE),  
ylab="Probability of Disagreeing", xlab="Party ID", main="Probability  
of Disagreeing That The Carbon Tax Would Hurt The Economy")

dev.off()



*
Simon J. Kiss, PhD
SSHRC and DAAD Post-Doctoral Fellow
John F. Kennedy Institute of North America Studies
Free University of Berlin
Lansstraße 7-9
14195 Berlin, Germany
Cell: +49 (0)1525-300-2812,
Web: http://www.jfki.fu-berlin.de/index.html








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


Re: [R] Mixed Effects Model on Within-Subjects Design

2010-05-20 Thread Dave Deriso
Hi Thierry,

Thank you so much for your response! I ran the model and I obtained
some strange results (see below). Is there a simple way to compute a
condition x difference interaction with the lme? Also, I read in the R
Book (Crawley, 2007) that repeated measures on the same day would be
temporal pseudoreplication. Being that I have 3 repeated measures for
each condition x difference pair, I assumed that each measurement
(variable "rep") would be a random effect variable along with subject.
Is this a correct assumption? If not, should I switch back to a
repeated measures ANOVA?

Thank you so much!!

Best,
Dave Deriso
UCSD Psychology

study.lme = lme(value~condition:diff - 1,random=~1|subject)
> summary(study.lme)
Linear mixed-effects model fit by REML
 Data: NULL
  AIC  BIClogLik
 19354.54 19405.71 -9666.272

Random effects:
 Formula: ~1 | subject
   (Intercept) Residual
StdDev:37786.52 59827.67

Fixed effects: value ~ condition:diff - 1
Value Std.Error  DF  t-value p-value
condition1:diff25 300110.0   9506.69 746 31.56829   0
condition2:diff25 327827.6   9506.69 746 34.48388   0
condition3:diff25 276391.2   9506.69 746 29.07334   0
condition1:diff50 356877.5   9506.69 746 37.53962   0
condition2:diff50 339113.9   9506.69 746 35.67108   0
condition3:diff50 340492.1   9506.69 746 35.81606   0
condition1:diff75 420141.8   9506.69 746 44.19432   0
condition2:diff75 409093.6   9506.69 746 43.03218   0
condition3:diff75 409342.6   9506.69 746 43.05837   0
 Correlation:
 cn1:25 cn2:25 cn3:25 cn1:50 cn2:50 cn3:50 cn1:75 cn2:75
condition2:diff25 0.545
condition3:diff25 0.545  0.545
condition1:diff50 0.545  0.545  0.545
condition2:diff50 0.545  0.545  0.545  0.545
condition3:diff50 0.545  0.545  0.545  0.545  0.545
condition1:diff75 0.545  0.545  0.545  0.545  0.545  0.545
condition2:diff75 0.545  0.545  0.545  0.545  0.545  0.545  0.545
condition3:diff75 0.545  0.545  0.545  0.545  0.545  0.545  0.545  0.545

Standardized Within-Group Residuals:
   Min  Q1 Med  Q3 Max
-6.15651930 -0.55808827 -0.02570532  0.52282057  5.06724039

Number of Observations: 783
Number of Groups: 29

On Thu, May 20, 2010 at 2:01 AM, ONKELINX, Thierry
 wrote:
> Dear Dave,
>
> I think you want this model.
>
> lme(value~condition:diff - 1,random=~1|subject)
>
> Note that I removed the replicate ID from the model. Include it in the model 
> makes only sense if you can expect a similar replication effects the 
> first/second/thirth time that a subject performs your test.
>
> HTH,
>
> Thierry
>
> 
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek
> team Biometrie & Kwaliteitszorg
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium
>
> Research Institute for Nature and Forest
> team Biometrics & Quality Assurance
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium
>
> tel. + 32 54/436 185
> thierry.onkel...@inbo.be
> www.inbo.be
>
> To call in the statistician after the experiment is done may be no more than 
> asking him to perform a post-mortem examination: he may be able to say what 
> the experiment died of.
> ~ Sir Ronald Aylmer Fisher
>
> The plural of anecdote is not data.
> ~ Roger Brinner
>
> The combination of some data and an aching desire for an answer does not 
> ensure that a reasonable answer can be extracted from a given body of data.
> ~ John Tukey
>
>
>> -Oorspronkelijk bericht-
>> Van: r-help-boun...@r-project.org
>> [mailto:r-help-boun...@r-project.org] Namens Dave Deriso
>> Verzonden: donderdag 20 mei 2010 9:27
>> Aan: David Atkins
>> CC: r-help@r-project.org
>> Onderwerp: Re: [R] Mixed Effects Model on Within-Subjects Design
>>
>> Hi Dave,
>>
>> Thank you for your helpful advice. I will take a look at the
>> multicomp package.
>>
>> I was wondering where the lme() function outputs the
>> interaction between condition*difficulty?
>>
>> Below is the output to the code I had originally sent. Which
>> one of these is condition*difficulty?
>>
>> Fixed effects: value ~ condition * diff
>>                       Value Std.Error  DF   t-value p-value
>> (Intercept)       300109.95  9506.690 688 31.568289  0.
>> condition2         27717.65  9071.048 688  3.055617  0.0023
>> condition3        -23718.72  9071.048 688 -2.614772  0.0091
>> diff50             56767.55  9071.048 688  6.258103  0.
>> diff75            120031.80  9071.048 688 13.232408  0.
>> condition2:diff50 -45481.21 12828.399 688 -3.545354  0.0004
>> condition3:diff50   7333.37 12828.399 688  0.571651  0.5677
>> condition2:diff75 -38765.77 12828.399 688 -3.021871  0.0026
>> condition3:diff75  12919.59 12828.399 688  1.007109  0.3142
>>
>> Also, why are diff25 and condition1 missing from the output??
>>
>> Thanks again for your generous help!!!
>>
>> Best,
>> Dave Deriso
>>
>>
>> On Wed, May 19, 2010 at 10:08 PM, David Atkins
>>  wrote:
>> >
>>

[R] Geneland error on unix: Error in MCMC(........ :, unused argument(s) (ploidy = 2, genotypes = geno)

2010-05-20 Thread Nevil Amos


I am receiving the above error ( full r session output below)  the 
script runs OK in windows.  and "genotypes" and "ploidy" are both 
correct arguments


any suggestions would be most welcome

Nevil Amos
MERG/ACB
Monash University School of Biological Sciences







> library(Geneland)
Loading required package: RandomFields
Loading required package: fields
Loading required package: spam
Package 'spam' is loaded. Spam version 0.21-0 (2010-03-13).
Type demo( spam) for some demos, help( Spam) for an overview
of this package.
Help for individual functions is optained by adding the
suffix '.spam' to the function name, e.g. 'help(chol.spam)'.

Attaching package: 'spam'

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

backsolve, forwardsolve, norm

 Try help(fields) for an overview of this library
fields web: http://www.image.ucar.edu/Software/Fields
Loading required package: mapproj
Loading required package: maps
Loading required package: snow
Loading required package: tcltk
Loading Tcl/Tk interface ... done
ooo
oGeneland is loaded   o
o o
o* Please *   o
o o
oRegister on  o
ohttp://www2.imm.dtu.dk/~gigu/Geneland/register.php   o
o o
oSee manual ono
ohttp://www2.imm.dtu.dk/~gigu/Geneland/#Manualo
o o
oType citation("Geneland") for a quick citation guide o
o o
oSee http://www2.imm.dtu.dk/~gigu/Geneland/#  o
ofor additional referenceso
o o
oThis is Geneland-3.2.1   o
o o
ooo
Warning message:
In fun(...) : no DISPLAY variable so Tk is not available
> DIRLIST<-c("Adult_ALL 
NO_ANW/")#,"Adult_Or_ANW/","Adult_Females/","Adult_Males/")

> for(d in DIRLIST){
+ theWd<- paste("/nfs/monash/home/namos/Rwork/",d,sep="")
+ setwd(theWd)
+ SPP.CODES <-"EYR"#c("BT","EYR","FH","SPP","STP")
+ for (sp in  SPP.CODES){
+ path.sp<- paste(theWd,sp,"/",sep="")
+ dir.create(path.sp)
+ GENO.TABLE<-paste(theWd,sp,"geno",sep="")
+ XY.TABLE<-paste(theWd,sp,"xy",sep="")
+ MINPOP=1
+ MAXPOP=25
+ INITPOP=1
+ NITS=50
+ THIN=NITS/1000
+ nrun <- 10
+ burnin <- 200
+ geno<-noquote(read.table(GENO.TABLE))
+ coord<-read.table(XY.TABLE)
+
+ ## Loop for multiple runs
+
+ for(irun in 1:nrun)
+ {
+ ## define path to MCMC directory
+
+ path.mcmc <- paste(path.sp,irun,"/",sep="")
+ dir.create(path.mcmc)
+ MCMC(coordinates=coord, ploidy=2, genotypes=geno, varnpop=T, 
npopmax=MAXPOP, npopinit=MINPOP, spatial=T, freq.model="Correlated", 
nit=NITS, thinning=500, rate.max=nrow(geno), delta.coord=100, 
path.mcmc=path.mcmc)

+ ## MCMC postprocessing
+ PostProcessChain(coordinates=coord,path.mcmc=path.mcmc,genotypes=geno, 
nxdom=50,nydom=50,burnin=burnin)

+ }
+ ## Computing average posterior probability
+
+ lpd <- rep(NA,nrun)
+ for(irun in 1:nrun)
+ {
+ path.lpd <- paste(path.mcmc,"log.posterior.density.txt",sep="")
+ lpd[irun] <- mean(scan(path.lpd)[-(1:burnin)])
+ }
+ ## Runs sorted by decreasing average posterior probability:
+ order(lpd,decreasing=TRUE)
+
+ #pdf("Number of pops.pdf",((210-10)/25.4),((297-10)/25.4))
+ #theWd<-getwd()
+ ##setwd(theWd)
+ #par(mfrow = c(5, 2),cex=0.25)
+ #for(irun in 1:nrun)
+ #
+ #{
+ #
+ ###Below is code form the Geneland Plotnpop function
+ #fileparam <- paste(path.mcmc, "parameters.txt", sep = "")
+ #param <- as.matrix(read.table(fileparam))
+ #thinning <- as.numeric(param[param[, 1] == "thinning", 3])
+ #filenpop <- paste(path.mcmc, "populations.numbers.txt", sep = "")
+ #npop <- scan(filenpop)
+ # sub <- -(1:burnin)
+ #
+ #plot((1:length(npop)) * thinning, npop, type = "l",ylab = "Number of 
classes", xlab = "Index of MCMC iteration Whole 
chain",main=paste(path.mcmc,"\n","Run:", irun, sep = " "), ylim = c(1, 
max(npop)+ 0.5))
+ #hist(npop[sub], plot = TRUE, prob = TRUE, breaks = seq(0.5,max(npop) 
+ 0.5, 1), xlab = paste("Nb. of pop. along the chain (after a burnin of 
", burnin, "x", thinning, "i t.)", sep = ""), main = "Number of 
populations long the chain after burnin")

+ #
+ #
+ #}
+ #dev.off()
+ }
+ }
Error in MCMC(coordinates = coord, ploidy = 2, genotypes = geno, varnpop 
= T,  :

  unused argument(s) (ploidy = 2, genotypes = geno)
Execution halted

_

[R] Reading results of commands in Microsoft Word typed in the terminal window, A question from a Blind R user.

2010-05-20 Thread Faiz Rasool
Hi all, 

I would like to read the results of the commands type in the terminal window in 
Microsoft Word. As a blind user my options are somewhat limited and are time 
consuming if I want to see the results of the commands that I have type 
earlier. for example if  my first two commands were
 x<-c(1,2,3,4,5)
mean(x)
and I have typed ten more commands after the first two commands it is not easy 
for me to see that what was the result of mean(x)
but if I can somehow divert the results of the commands to Microsoft Word it is 
comparatively easy for me to see what was the result of mean(x) and what were 
the results of other commands. One another advantage of diverting R's output to 
Microsoft Word for me is that from there they can be easily copied into 
assignments as well.

Any ideas and suggestions are appreciated.

faiz.







[[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] Mixed Effects Model on Within-Subjects Design

2010-05-20 Thread ONKELINX, Thierry
Dear Dave,

The model is correct. But it is not testing the hypotheses you want. The null 
hypothesis is that each level of the interaction is zero, which is clearly not 
the fact.

I misread you question, therefore my advise was not accurate. Since you are 
explicitly interested in comparing all levels of the interaction, I would 
suggest that you create a new variable that defines the interaction and then 
use a multiple comparison test.

library(nlme)
library(multcomp)
study.data <- read.csv("http://files.davidderiso.com/example_data.csv";, 
header=T)
study.data$combination <- with(study.data, factor(condition):factor(diff))
study.lme <- lme(value ~ combination - 1,random=~1|subject, data = study.data)
summary(glht(study.lme, linfct = mcp(combination = "Tukey")))

HTH,

Thierry


ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium

Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium

tel. + 32 54/436 185
thierry.onkel...@inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than 
asking him to perform a post-mortem examination: he may be able to say what the 
experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure 
that a reasonable answer can be extracted from a given body of data.
~ John Tukey
  

> -Oorspronkelijk bericht-
> Van: der...@gmail.com [mailto:der...@gmail.com] Namens Dave Deriso
> Verzonden: donderdag 20 mei 2010 11:24
> Aan: ONKELINX, Thierry
> CC: r-help@r-project.org
> Onderwerp: Re: [R] Mixed Effects Model on Within-Subjects Design
> 
> Hi Thierry,
> 
> Thank you so much for your response! I ran the model and I 
> obtained some strange results (see below). Is there a simple 
> way to compute a condition x difference interaction with the 
> lme? Also, I read in the R Book (Crawley, 2007) that repeated 
> measures on the same day would be temporal pseudoreplication. 
> Being that I have 3 repeated measures for each condition x 
> difference pair, I assumed that each measurement (variable 
> "rep") would be a random effect variable along with subject.
> Is this a correct assumption? If not, should I switch back to 
> a repeated measures ANOVA?
> 
> Thank you so much!!
> 
> Best,
> Dave Deriso
> UCSD Psychology
> 
> study.lme = lme(value~condition:diff - 1,random=~1|subject)
> > summary(study.lme)
> Linear mixed-effects model fit by REML
>  Data: NULL
>   AIC  BIClogLik
>  19354.54 19405.71 -9666.272
> 
> Random effects:
>  Formula: ~1 | subject
>(Intercept) Residual
> StdDev:37786.52 59827.67
> 
> Fixed effects: value ~ condition:diff - 1
> Value Std.Error  DF  t-value p-value
> condition1:diff25 300110.0   9506.69 746 31.56829   0
> condition2:diff25 327827.6   9506.69 746 34.48388   0
> condition3:diff25 276391.2   9506.69 746 29.07334   0
> condition1:diff50 356877.5   9506.69 746 37.53962   0
> condition2:diff50 339113.9   9506.69 746 35.67108   0
> condition3:diff50 340492.1   9506.69 746 35.81606   0
> condition1:diff75 420141.8   9506.69 746 44.19432   0
> condition2:diff75 409093.6   9506.69 746 43.03218   0
> condition3:diff75 409342.6   9506.69 746 43.05837   0
>  Correlation:
>  cn1:25 cn2:25 cn3:25 cn1:50 cn2:50 cn3:50 
> cn1:75 cn2:75
> condition2:diff25 0.545
> condition3:diff25 0.545  0.545
> condition1:diff50 0.545  0.545  0.545
> condition2:diff50 0.545  0.545  0.545  0.545 
> condition3:diff50 0.545  0.545  0.545  0.545  0.545
> condition1:diff75 0.545  0.545  0.545  0.545  0.545  0.545
> condition2:diff75 0.545  0.545  0.545  0.545  0.545  0.545  0.545
> condition3:diff75 0.545  0.545  0.545  0.545  0.545  0.545  
> 0.545  0.545
> 
> Standardized Within-Group Residuals:
>Min  Q1 Med  Q3 Max
> -6.15651930 -0.55808827 -0.02570532  0.52282057  5.06724039
> 
> Number of Observations: 783
> Number of Groups: 29
> 
> On Thu, May 20, 2010 at 2:01 AM, ONKELINX, Thierry
>  wrote:
> > Dear Dave,
> >
> > I think you want this model.
> >
> > lme(value~condition:diff - 1,random=~1|subject)
> >
> > Note that I removed the replicate ID from the model. 
> Include it in the model makes only sense if you can expect a 
> similar replication effects the first/second/thirth time that 
> a subject performs your test.
> >
> > HTH,
> >
> > Thierry
> >
> > 
> --
> > --
> > ir. Thierry Onkelinx
> > Instituut voor natuur- en bosonderzoek team Biometrie & 
> Kwaliteitszorg 
> > Gaverstraat 4 9500 Geraardsbergen Belgium
> >
> > Research Institute for Nature and Forest team Biometrics & Quality 
> > Assur

[R] Sweave and uttf-8 under Windows XP

2010-05-20 Thread Sven Garbade
Hi list,

I need to process a Rnw file and and a csv file (both are encoded in
UTF-8) under Windows XP  (R Version 2.11.0, i386, mingw32).

I can source and run the Rnw file:

> Stangle("Bericht.Rnw")
Writing to file Bericht.R
> source(file("Bericht.R", encoding="UTF-8"))

which runs fine, but running Sweave() failed:

> Sweave("Bericht.Rnw")
Writing to file Bericht.tex
Processing code chunks ...
 1 : term verbatim
 2 : term verbatim

Fehler:  chunk 2
Error in parse(text = chunk) : Unerwartetes Eingabe in:
"
data$X6.1.Staats.angehö"
>

My first idea was to run Sweave like Sweave(file("Bericht.Rnw",
encoding="UTF-8")), but Sweave seems to expect a file name as a
character rather a connection. The Rnw and csv file were created under
Linux in a UTF-8 locale and were further processed by xelatex.

Are there any ideas to get Sweave working? I would be nice when the
file encoding can be left in UTF-8.

Thanks, Sven

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


Re: [R] Geneland error on unix: Error in MCMC(........ :, unused argument(s) (ploidy = 2, genotypes = geno)

2010-05-20 Thread Jombart, Thibaut
Hello,

you may want to contact directly the maintainer of Geneland for that kind of 
issue. In any case, this post would be best suited for R-sig-genetics:
https://stat.ethz.ch/mailman/listinfo/r-sig-genetics

Best regards,

Thibaut

--
##
Dr Thibaut JOMBART
MRC Centre for Outbreak Analysis and Modelling
Department of Infectious Disease Epidemiology
Imperial College - Faculty of Medicine
St Mary’s Campus
Norfolk Place
London W2 1PG
United Kingdom
Tel. : 0044 (0)20 7594 3658
t.jomb...@imperial.ac.uk
http://sites.google.com/site/thibautjombart/
http://adegenet.r-forge.r-project.org/


From: Nevil Amos [nevil.a...@gmail.com]
Sent: 20 May 2010 10:34
To: r-help@r-project.org; Jombart, Thibaut
Subject: Geneland error on unix:  Error in MCMC( :,  unused argument(s) 
(ploidy = 2, genotypes = geno)

I am receiving the above error ( full r session output below)  the
script runs OK in windows.  and "genotypes" and "ploidy" are both
correct arguments

any suggestions would be most welcome

Nevil Amos
MERG/ACB
Monash University School of Biological Sciences







 > library(Geneland)
Loading required package: RandomFields
Loading required package: fields
Loading required package: spam
Package 'spam' is loaded. Spam version 0.21-0 (2010-03-13).
Type demo( spam) for some demos, help( Spam) for an overview
of this package.
Help for individual functions is optained by adding the
suffix '.spam' to the function name, e.g. 'help(chol.spam)'.

Attaching package: 'spam'

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

 backsolve, forwardsolve, norm

  Try help(fields) for an overview of this library
fields web: http://www.image.ucar.edu/Software/Fields
Loading required package: mapproj
Loading required package: maps
Loading required package: snow
Loading required package: tcltk
Loading Tcl/Tk interface ... done
ooo
oGeneland is loaded   o
o o
o* Please *   o
o o
oRegister on  o
ohttp://www2.imm.dtu.dk/~gigu/Geneland/register.php   o
o o
oSee manual ono
ohttp://www2.imm.dtu.dk/~gigu/Geneland/#Manualo
o o
oType citation("Geneland") for a quick citation guide o
o o
oSee http://www2.imm.dtu.dk/~gigu/Geneland/#  o
ofor additional referenceso
o o
oThis is Geneland-3.2.1   o
o o
ooo
Warning message:
In fun(...) : no DISPLAY variable so Tk is not available
 > DIRLIST<-c("Adult_ALL
NO_ANW/")#,"Adult_Or_ANW/","Adult_Females/","Adult_Males/")
 > for(d in DIRLIST){
+ theWd<- paste("/nfs/monash/home/namos/Rwork/",d,sep="")
+ setwd(theWd)
+ SPP.CODES <-"EYR"#c("BT","EYR","FH","SPP","STP")
+ for (sp in  SPP.CODES){
+ path.sp<- paste(theWd,sp,"/",sep="")
+ dir.create(path.sp)
+ GENO.TABLE<-paste(theWd,sp,"geno",sep="")
+ XY.TABLE<-paste(theWd,sp,"xy",sep="")
+ MINPOP=1
+ MAXPOP=25
+ INITPOP=1
+ NITS=50
+ THIN=NITS/1000
+ nrun <- 10
+ burnin <- 200
+ geno<-noquote(read.table(GENO.TABLE))
+ coord<-read.table(XY.TABLE)
+
+ ## Loop for multiple runs
+
+ for(irun in 1:nrun)
+ {
+ ## define path to MCMC directory
+
+ path.mcmc <- paste(path.sp,irun,"/",sep="")
+ dir.create(path.mcmc)
+ MCMC(coordinates=coord, ploidy=2, genotypes=geno, varnpop=T,
npopmax=MAXPOP, npopinit=MINPOP, spatial=T, freq.model="Correlated",
nit=NITS, thinning=500, rate.max=nrow(geno), delta.coord=100,
path.mcmc=path.mcmc)
+ ## MCMC postprocessing
+ PostProcessChain(coordinates=coord,path.mcmc=path.mcmc,genotypes=geno,
nxdom=50,nydom=50,burnin=burnin)
+ }
+ ## Computing average posterior probability
+
+ lpd <- rep(NA,nrun)
+ for(irun in 1:nrun)
+ {
+ path.lpd <- paste(path.mcmc,"log.posterior.density.txt",sep="")
+ lpd[irun] <- mean(scan(path.lpd)[-(1:burnin)])
+ }
+ ## Runs sorted by decreasing average posterior probability:
+ order(lpd,decreasing=TRUE)
+
+ #pdf("Number of pops.pdf",((210-10)/25.4),((297-10)/25.4))
+ #theWd<-getwd()
+ ##setwd(theWd)
+ #par(mfrow = c(5, 2),cex=0.25)
+ #for(irun in 1:nrun)
+ #
+ #{
+ #
+ ###Below is code form the Geneland Plotnpop function
+ #fileparam <- paste(path.mcmc, "parameters.txt", sep = "")
+ #param <- as.matrix(read.table(fi

Re: [R] Minimization problem

2010-05-20 Thread Bernardo Rangel Tura
On Thu, 2010-05-20 at 01:35 -0700, Fred wrote:
> Dear R users,
> 
> I am trying to minimize two function simultaneously in R,
> 
> function(x)
> 
> minimize x[1],x[2],x[3]
> 
> mean(distribution(x1,x2,x3) ) - observed mean
> 
> std(distribution(x1,x2,x3)) - observed std
> 
> What I want to achieve is that simulated mean and standard deviation
> of distribution related to x1 x2 x3  would be close to observed mean
> and observed standard deviation.
> 
> 
> is there any function in R can reach this?
> 
> Thank you for the help first .
> 
> F.

Hi!
Do you need use optim, something like this

test <- function(parameters){
m.error <- mean(distribution(x1,x2,x3) ) - observed mean
m.sd <- std(distribution(x1,x2,x3)) - observed std
res <- cbind(m.error,sd.error)
return(res)
}
-- 
Bernardo Rangel Tura, M.D,MPH,Ph.D
National Institute of Cardiology
Brazil

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


Re: [R] Sweave and uttf-8 under Windows XP

2010-05-20 Thread Sven Garbade
Great, thanks, Sveb

2010/5/20 Duncan Murdoch :
> Sven Garbade wrote:
>>
>> Hi list,
>>
>> I need to process a Rnw file and and a csv file (both are encoded in
>> UTF-8) under Windows XP  (R Version 2.11.0, i386, mingw32).
>>
>
> Take a look at this post:
>
> http://tolstoy.newcastle.edu.au/R/e10/help/10/05/4889.html
>
> which discussed this issue recently.
>
> Duncan Murdoch
>>
>> I can source and run the Rnw file:
>>
>>
>>>
>>> Stangle("Bericht.Rnw")
>>>
>>
>> Writing to file Bericht.R
>>
>>>
>>> source(file("Bericht.R", encoding="UTF-8"))
>>>
>>
>> which runs fine, but running Sweave() failed:
>>
>>
>>>
>>> Sweave("Bericht.Rnw")
>>>
>>
>> Writing to file Bericht.tex
>> Processing code chunks ...
>>  1 : term verbatim
>>  2 : term verbatim
>>
>> Fehler:  chunk 2
>> Error in parse(text = chunk) : Unerwartetes Eingabe in:
>> "
>> data$X6.1.Staats.angehö"
>>
>> My first idea was to run Sweave like Sweave(file("Bericht.Rnw",
>> encoding="UTF-8")), but Sweave seems to expect a file name as a
>> character rather a connection. The Rnw and csv file were created under
>> Linux in a UTF-8 locale and were further processed by xelatex.
>>
>> Are there any ideas to get Sweave working? I would be nice when the
>> file encoding can be left in UTF-8.
>>
>> Thanks, Sven
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>

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


[R] R 2.11.1 Scheduled for May 31

2010-05-20 Thread Peter Dalgaard
This is to announce that we plan to release R version 2.11.1 on Monday,
May 31, 2010.

Those directly involved should review the generic schedule at
http://developer.r-project.org/release-checklist.html

The source tarballs will be made available daily (barring build
troubles) via

http://cran.r-project.org/src/base-prerelease/

  For the R Core Team
  Peter Dalgaard

-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.com

___
r-annou...@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-announce

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

2010-05-20 Thread Maulik Shah
May I ask for the help on how to plot Item characteristic curve (ICC) one
item each. I have fit 3 parameter model for a test with 40 items.Currently I
get ICC for all the 40 items in a single graph making it difficult to
interpret.

Thanks and Regards,
-- 
Maulik Shah

[[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] Sweave and uttf-8 under Windows XP

2010-05-20 Thread Duncan Murdoch

Sven Garbade wrote:

Hi list,

I need to process a Rnw file and and a csv file (both are encoded in
UTF-8) under Windows XP  (R Version 2.11.0, i386, mingw32).
  


Take a look at this post:

http://tolstoy.newcastle.edu.au/R/e10/help/10/05/4889.html

which discussed this issue recently.

Duncan Murdoch

I can source and run the Rnw file:

  

Stangle("Bericht.Rnw")


Writing to file Bericht.R
  

source(file("Bericht.R", encoding="UTF-8"))



which runs fine, but running Sweave() failed:

  

Sweave("Bericht.Rnw")


Writing to file Bericht.tex
Processing code chunks ...
 1 : term verbatim
 2 : term verbatim

Fehler:  chunk 2
Error in parse(text = chunk) : Unerwartetes Eingabe in:
"
data$X6.1.Staats.angehö"
  


My first idea was to run Sweave like Sweave(file("Bericht.Rnw",
encoding="UTF-8")), but Sweave seems to expect a file name as a
character rather a connection. The Rnw and csv file were created under
Linux in a UTF-8 locale and were further processed by xelatex.

Are there any ideas to get Sweave working? I would be nice when the
file encoding can be left in UTF-8.

Thanks, Sven

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

2010-05-20 Thread Mohan L
Dear All,

I have data some thisng like this :

> data <- read.csv(file='ipsample.csv',sep=',' , header=TRUE)

> data
  State  Jan  Feb  Mar  Apr  May Jun
1   AAA11022   0
2   BBB 1298 1195 1212 1244 1158 845
3  CCC 00012   1
4   DDD5   11   17   15   10   9
5   EEE   18   28   27   23   23  16
6   FFF   68  152  184  135  111  86

from this data frame, I took "Jan" as base and calculating weightage like
this :

> basemonth.sum <- sum(data[[2]])

> basemonth.sum
[1] 1390

> basemonth.data <- data[[2]]

>  basemonth.data
[1]1 129805   18   68

> weightage <- basemonth.data / basemonth.sum

> weightage
[1] 0.0007194245 0.9338129496 0.00 0.0035971223 0.0129496403
[6] 0.0489208633


The above is the weightage for base month "Jan".  Now I need to calculate
weighted states data. What I need to do is :
(((Feb[i]-Jan[1])*weightage)+Jan[1]) for all column. The "Jan" column is
fixed. I need to do the calculation in all the column Feb, Mar etc...


 StateJan   Feb
Mar
1   AAA1(((Feb[1]-Jan[1])*weightage[1])+Jan[1])
(((Mar[1]-Jan[1])*weightage[1])+Jan[1])
2   BBB 1298 (((Feb[2]-Jan[2])*weightage[2])+Jan[2])
(((Mar[2]-Jan[1])*weightage[2])+Ja[1])
3   CCC 0
4   DDD5
5   EEE   18
6   FFF   68

I am struggling with this . I have framed a logic using for loop. But it
seems me very bad logic.  Any help will be greatly appreciated.

Thanks & Rg
Mohan L

[[alternative HTML version deleted]]

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


[R] [R-pkgs] New package: `lavaan' for latent variable analysis (including structural equation modeling)

2010-05-20 Thread Yves Rosseel

Dear R-users,

A new package called `lavaan' (for latent variable analysis) has been 
uploaded to CRAN. The current version of lavaan (0.3-1) can be used for 
path analysis, confirmatory factor analysis, structural equation 
modeling, and growth curve modeling.


More information can be found on the website: http://lavaan.org

Some notable features of lavaan:

- the 'lavaan model syntax' allows users to express their models in a 
compact, elegant and useR-friendly way


- lavaan is robust and reliable: there are no convergence problems and 
numerical results are very close (if not identical) to the commercial 
package Mplus


- many different estimators are available: ML, GLS, WLS, robust ML using 
Satorra-Bentler corrections, and FIML for data with missing values


- full support for meanstructures and multiple groups

- user friendly output including standardized solutions, fit measures, 
modification indices and more


To get a first impression of how the 'lavaan model syntax' looks like, 
below is the full R code for fitting a SEM model:


## begin R Code ##

library(lavaan)

# The industrialization and Political Democracy Example
# Bollen (1989), page 332

model <- '
  # latent variable definitions
 ind60 =~ x1 + x2 + x3
 dem60 =~ y1 + y2 + y3 + y4
 dem65 =~ y5 + y6 + y7 + y8

  # regressions
dem60 ~ ind60
dem65 ~ ind60 + dem60

  # residual correlations
y1 ~~ y5
y2 ~~ y4 + y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8
'

fit <- sem(model, data=PoliticalDemocracy)
summary(fit, fit.measures=TRUE)

## end R code ##


Please feel free to contact me directly with questions and comments.

Best,

Yves Rosseel.



--
Yves Rosseel -- http://www.da.ugent.be
Department of Data Analysis, Ghent University
Henri Dunantlaan 1, B-9000 Gent, Belgium
-

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

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


Re: [R] writing autocorrelation and partial auto correlation functions to a file

2010-05-20 Thread Glen Barnett
The problem is that the acf function (like many R functions) returns a list
containing many different things. For example, I have a short series in
the vector z:

> acz <- acf(z)
> str(acz)
List of 6
 $ acf   : num [1:11, 1, 1] 1 -0.0668 -0.7401 0.0627 0.5954 ...
 $ type  : chr "correlation"
 $ n.used: int 11
 $ lag   : num [1:11, 1, 1] 0 1 2 3 4 5 6 7 8 9 ...
 $ series: chr "z"
 $ snames: NULL
 - attr(*, "class")= chr "acf"

Let's say you are interested in extracting the acf and lag values from that

> acd <- data.frame(cbind(lag=acz$lag,acf=acz$acf))
> acd
   lag acf
10  1.
21 -0.06676557
32 -0.74006273
43  0.06266564
54  0.59542482
65 -0.06786556
76 -0.41370293
87  0.03464621
98  0.21887326
10   9 -0.07485743
11  10 -0.04835570

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 could I restrict and reordered data.frames?

2010-05-20 Thread Csima Gabriella
Dear Everyone,

 

I 've just begun to use the library ncdf and I would like to compare 
meteorological observational data with forecast data, so to make verification. 
The netcdf files I'm using contain data of many different parameters in many 
different stations.  I could read easily that I needed, but naturally I do not 
need the data of all the stations. On the other hand, the order of the stations 
is not the same in the observation files and in the forecast files.

Let's take that I have a list of those stations (with station numbers) where I 
would like to make the verification. I read the observations in all the 
possible stations and I receive a data.frame (first column with the station 
numbers, second column with - let' say - the temperature data...and naturally 
we can have more columns with different parameters). I make the same with the 
forecast data, as I wrote the orders of the station numbers in the two dataset 
are different, and naturally there are some stations that you can find in one 
data.frame but not in the other.



How could I make (or rewrite) my two data.frames (observation and forecast), 
where the first coulumn is totally the same as in the station list (even the 
order of the stations)??

 

For example, I have this data.frame as obsesrvation:



12866   14.4

12844   14.1

12843   16.5
12860   14.9
128519.8

12846   15.3





...and  have this data.frame as forecast:



12830   12.808611
12836   12.725081
12843   15.241580
12844   15.185887
12846   13.723515
128518.498717
12860   15.715260
12866   14.262023
12870   12.968392


...and my list of the stations I needed 



12836   
12843
12846
12860
12870   




Thank you very much for your help or suggestions in advance!

Gabriella Csima

csim...@met.hu 

 

[[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] pie3D

2010-05-20 Thread n.via...@libero.it
Dear list,
I would like to know if is it possible to add a legend to a pie3D??
Thanks for your attention!

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


Re: [R] pie3D

2010-05-20 Thread Jim Lemon

On 05/20/2010 09:06 PM, n.via...@libero.it wrote:

Dear list,
I would like to know if is it possible to add a legend to a pie3D??
Thanks for your attention!


Hi n.vialma,
You sure can. Try adding this to the example:

legend(-1,1,c("Haters","Opposers","Apathizers","Lovers"),fill=rainbow(4))

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.


[R] Overlap of leaf labels

2010-05-20 Thread Ayesha Jadoon
Hi,

I have tried looking at the archives but havent found any answer that works
till now (Sorry if i have missed anything)


I am a newbie to R and i am trying to carry out hierarchical clustering
using hclust -> as.dendrogram and then plotting the results as a dendrogram
using the plot function plot(object).

My question is :

In the function "plot", can one decrease the leaf label size to make them
readable and clear? I am including over 380 proteins in my dendrogram
and each leaf has a label which are currently overlapping and not
decipherable?


Thanks


Ayesha

[[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] Overlap of leaf labels

2010-05-20 Thread Ivan Calandra

Hi,
I think that one of the cex arguments in par() can be what you're 
looking for. But since I've never plotted any dendrogram, I don't know 
which one, if any.

HTH,
Ivan

Le 5/20/2010 14:08, Ayesha Jadoon a écrit :

Hi,

I have tried looking at the archives but havent found any answer that works
till now (Sorry if i have missed anything)


I am a newbie to R and i am trying to carry out hierarchical clustering
using hclust ->  as.dendrogram and then plotting the results as a dendrogram
using the plot function plot(object).

My question is :

In the function "plot", can one decrease the leaf label size to make them
readable and clear? I am including over 380 proteins in my dendrogram
and each leaf has a label which are currently overlapping and not
decipherable?


Thanks


Ayesha

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

   


--
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php

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


[R] use object within rda file in for loop

2010-05-20 Thread Ivan Calandra

Dear users,

I would like to process all the lists from all *.rda files that I have 
in one folder.

Up to now, I can load all the *.rda files without any problem.

The problem is when I want to access the list saved within each *.rda 
file (only one list per rda file).


Here is my code:
fpath <- "D:/R"
listnames <- list.files(path=fpath, pattern=glob2rx("*.rda"), 
full.names=FALSE)

for (i in 1:length(listnames)) {
  load(paste(fpath, listnames[i], sep="/"))
  z <- list_in_listnames[i]  ## here is my problem
  **do something on z**
}

It might be really simple, but listnames is a character vector and I 
cannot find how to store the values within each element into z.
I think there would be a function to get and use an object with a given 
pattern in its name. The pattern itself is no problem, the problem is 
the function.
I've tried to look on RSiteSearch but, probably because I couldn't 
figure out the keywords, haven't found anything helpful.

Maybe I'm just on the wrong path to do it or missed something obvious...

I hope my question is clear. If not, please let me know what would help 
you to understand.


Thanks in advance
Ivan

--
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php

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


[R] logarithmic integrals in R?

2010-05-20 Thread Oliver Kullmann
Hello,

I couldn't find information on whether the logarithmic integrals

Li_m(x) = integral_0^x log(t)^(-m) dt

for x >= 0 are available in R?

Best wishes

Oliver

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 plot the sample size on top of the bars of a barchart plot?

2010-05-20 Thread Jabba

> 
> Dear Marco,
>  
> I was trying using lattice barchart() but your suggestion will do the
> trick.
>  
> Many thanks
>  
> Javier

Lattice is great but a bit hard for me to learn and customize. I think
you should use ltext() in this case or edit panel.barchart() by hand,
but I don't know how.

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


Re: [R] Re : Manipulating Data Frames

2010-05-20 Thread Henrique Dallazuanna
Try this:

(data[,-(1:2)] - data[,2]) * prop.table(data[,2]) + data[,2]


On Thu, May 20, 2010 at 7:38 AM, Mohan L  wrote:

> Dear All,
>
> I have data some thisng like this :
>
> > data <- read.csv(file='ipsample.csv',sep=',' , header=TRUE)
>
> > data
>  State  Jan  Feb  Mar  Apr  May Jun
> 1   AAA11022   0
> 2   BBB 1298 1195 1212 1244 1158 845
> 3  CCC 00012   1
> 4   DDD5   11   17   15   10   9
> 5   EEE   18   28   27   23   23  16
> 6   FFF   68  152  184  135  111  86
>
> from this data frame, I took "Jan" as base and calculating weightage like
> this :
>
> > basemonth.sum <- sum(data[[2]])
>
> > basemonth.sum
> [1] 1390
>
> > basemonth.data <- data[[2]]
>
> >  basemonth.data
> [1]1 129805   18   68
>
> > weightage <- basemonth.data / basemonth.sum
>
> > weightage
> [1] 0.0007194245 0.9338129496 0.00 0.0035971223 0.0129496403
> [6] 0.0489208633
>
>
> The above is the weightage for base month "Jan".  Now I need to calculate
> weighted states data. What I need to do is :
> (((Feb[i]-Jan[1])*weightage)+Jan[1]) for all column. The "Jan" column is
> fixed. I need to do the calculation in all the column Feb, Mar etc...
>
>
>  StateJan   Feb
>Mar
> 1   AAA1(((Feb[1]-Jan[1])*weightage[1])+Jan[1])
> (((Mar[1]-Jan[1])*weightage[1])+Jan[1])
> 2   BBB 1298 (((Feb[2]-Jan[2])*weightage[2])+Jan[2])
> (((Mar[2]-Jan[1])*weightage[2])+Ja[1])
> 3   CCC 0
> 4   DDD5
> 5   EEE   18
> 6   FFF   68
>
> I am struggling with this . I have framed a logic using for loop. But it
> seems me very bad logic.  Any help will be greatly appreciated.
>
> Thanks & Rg
> Mohan L
>
>[[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



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

[[alternative HTML version deleted]]

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


Re: [R] computer out of memory when using sigpathway

2010-05-20 Thread Martin Morgan
On 05/19/2010 09:15 PM, ?? wrote:
> Dear R users,
> 
> 
> 
> I am sorry to disturb you! But I really need your help for the usage
> of sigPathwy.

Please ask questions about Bioconductor packages on the Bioconductor
mailing list, as other Bioconductor users are the ones mostly likely to
have experience with this package.

  http://bioconductor.org/docs/mailList.html

Include the package maintainer in your email as they are in the best
position to provide informed advice.

  > packageDescription("sigPathways")$Maintainer

Include the output of sessionInfo() in your email, so that the specific
versions of R and packages are unambiguous.

  > sessionInfo()

Provide a minimal and self-contained example, so that those wishing to
help can easily reproduce your problem. This requires additional work on
your part, but will often reward by pointing out where the problem is.
Sending large attachments is NOT recommended; instead, use example data
from the package, from an appropriate experiment data package, or simulated.

When reporting an error, copy and paste the output of the relevant
section of your script. This provides additional clues to where exactly
a problem occurs.

After the error occurs, call

  > traceback()

to discover where during the computation the problem was.

Run

  > example(runSigPathway)

to see an example, and compare your input data to the input in the
example. For instance, is your 'yy' structured the same as 'G' in the
example?

Hope this leads to a solution,

Martin


> 
> 
> 
> Actually, I want a sliding window analysis for possible chromosome
> expression pattern mining. My research microorganism is a plant
> pathogen, Gibberella zeae, and I first used SAS to divide locus
> number with 10, 20, 30, or 40 on the fungal chromosome according to
> their location. I really want to see whether among the continual 10,
> 20, 30, or 40 locus has some expression pattern that different from
> rest genes. Because I know sigPathway (R package, pathway analysis
> with microarray data) can do this kind of job. What I use SAS to do
> is to subset locus in arbitrary genes numbers, such as 10, 20, 30,
> 40, or so on, and I hope to use sigPathway to analysis whether these
> genes chromosome location have effect on its gene expression.
> 
> When I used sigpathway to analyze my microarray data, it made my
> computer out of memory. I have tried the following R codes in several
> computer, but it always the same, even if it computed more than one
> day, it can not get any results. I also try to use
> memory.limit(size=NA), it dosen't work too.My computer is a cor duo
> 2G DDR2, I think it is large enough for my job. Would you please
> point out my problem and give me some suggestions? Thank you very
> much.
> 
> I attach my microarray data and R codes in the attachment, and I hope
> you can have a look.
> 
> 
> 
> #the following code is for annotation list initiation.
> 
> 
> 
> setwd("C:/analysis data and codes")
> 
> x <- read.table("chr1.txt",header=FALSE,sep="\t")
> 
> attach(x)
> 
> x$group <- paste(V2,V3,sep="_")
> 
> group <- x$group
> 
> y <- data.frame(group,V2,V3,V4)
> 
> xx <- as.list(group)
> 
> xx <- xx[!is.na(xx)]
> 
> xx <- unlist(xx)
> 
> xxUnique <- unique(xx)
> 
> yy <- vector("list",length(xxUnique))
> 
> for(i in 1:length(yy))
> 
> {
> 
> MT <- "MT_lab"
> 
> yy[[i]] <-
> list(src=MT,title=xxUnique[i],probes=as.character(y[group==xxUnique[i],]$V4))
>
>  }
> 
> 
> 
> 
> 
> #the following code is for sigpathway analysis.
> 
> 
> 
> library(sigPathway)
> 
> YANG <- read.table("All microarray MT_LAB.txt",header=T,sep="\t")
> 
> attach(YANG)
> 
> Y <-
> data.frame(TF134_1_3DAK,TF134_2_3DAK,WT1_3DAK,WT2_3DAK,row.names=locus_no)
>
>  p <- c("1_trt","1_trt","0_norm","0_norm")
> 
> statList <- calcTStatFast(Y,p,ngroups=2)
> 
> hist(statList$pval,breaks=seq(0,1,0.025),xlab="p-value",ylab="Frequency",main="")
>
>  set.seed(1234)
> 
> YANG <-
> runSigPathway(yy,20,500,Y,p,nsim=100,weightType="constant",ngroup=2,npath=10,verbose=F,allpathway=F,alwaysUseRandomPerm=F)
>
>  #write.table(YANG$df.pathways[1:25, ])
> 
> write.table(YANG$df.pathways[1:25,],quote=F,sep="\t",file="chr1_sig.txt")
>
>  YANG$list.gPS[[1]]
> 
> save.image("chr1_sig")
> 
> 
> 
> __ R-help@r-project.org
> mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do
> read the posting guide http://www.R-project.org/posting-guide.html 
> and provide commented, minimal, self-contained, reproducible code.
Rr

-- 
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793

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


Re: [R] Re : Adding column sum to new row in data frame

2010-05-20 Thread David Winsemius


On May 20, 2010, at 1:01 AM, Joshua Wiley wrote:


Dear Mohan,

First, I would like to modify my code slightly to:

data <- rbind(data,data.frame(State="Total",t(apply(data[,-1], 2, sum,
na.rm=TRUE

This actually will add a 7th level to your factor automatically.  The
reason I wanted to change from using c() to data.frame() is that if
one uses c(), all the columns are converted to character (this has to
do with different methods for rbind, see ?rbind particularly the
Details and Value section which describe the different methods for
rbind and what its behavior will be if it is using the default
method).  This may not be an issue, but it would hamper any subsequent
calculations you may wish to perform on your data.



Actually the coercion to a single element type occurs as soon as you  
use c()


> c("a", 1:10)
 [1] "a"  "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10"

Whether rbind might also coerce would depend on the class of arguments  
it is supplied. If you have created an argument with c() that is all  
character because that is the "least common denominator", and then try  
to rbind it to a dataframe, the numeric columns will character-ized.

--

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 could I restrict and reordered data.frames?

2010-05-20 Thread jim holtman
? merge

> obs <- read.table(textConnection("12866   14.4
+
+ 12844   14.1
+
+ 12843   16.5
+ 12860   14.9
+ 128519.8
+
+ 12846   15.3"), col.names=c('station', 'obs'))
> fore <- read.table(textConnection("12830   12.808611
+ 12836   12.725081
+ 12843   15.241580
+ 12844   15.185887
+ 12846   13.723515
+ 128518.498717
+ 12860   15.715260
+ 12866   14.262023
+ 12870   12.968392"), col.names=c('station', 'fore'))
> closeAllConnections()
>
> # use 'merge' to group them
>
> x <- merge(obs, fore, by="station", all=TRUE)
>
> x
  station  obs  fore
1   12830   NA 12.808611
2   12836   NA 12.725081
3   12843 16.5 15.241580
4   12844 14.1 15.185887
5   12846 15.3 13.723515
6   12851  9.8  8.498717
7   12860 14.9 15.715260
8   12866 14.4 14.262023
9   12870   NA 12.968392
>

You can then order by the first column.

On Thu, May 20, 2010 at 6:22 AM, Csima Gabriella  wrote:
> Dear Everyone,
>
>
>
> I 've just begun to use the library ncdf and I would like to compare 
> meteorological observational data with forecast data, so to make 
> verification. The netcdf files I'm using contain data of many different 
> parameters in many different stations.  I could read easily that I needed, 
> but naturally I do not need the data of all the stations. On the other hand, 
> the order of the stations is not the same in the observation files and in the 
> forecast files.
>
> Let's take that I have a list of those stations (with station numbers) where 
> I would like to make the verification. I read the observations in all the 
> possible stations and I receive a data.frame (first column with the station 
> numbers, second column with - let' say - the temperature data...and naturally 
> we can have more columns with different parameters). I make the same with the 
> forecast data, as I wrote the orders of the station numbers in the two 
> dataset are different, and naturally there are some stations that you can 
> find in one data.frame but not in the other.
>
>
>
> How could I make (or rewrite) my two data.frames (observation and forecast), 
> where the first coulumn is totally the same as in the station list (even the 
> order of the stations)??
>
>
>
> For example, I have this data.frame as obsesrvation:
>
>
>
> 12866   14.4
>
> 12844   14.1
>
> 12843   16.5
> 12860   14.9
> 12851    9.8
>
> 12846   15.3
>
>
>
>
>
> ...and  have this data.frame as forecast:
>
>
>
> 12830   12.808611
> 12836   12.725081
> 12843   15.241580
> 12844   15.185887
> 12846   13.723515
> 12851    8.498717
> 12860   15.715260
> 12866   14.262023
> 12870   12.968392
>
>
> ...and my list of the stations I needed
>
>
>
> 12836
> 12843
> 12846
> 12860
> 12870
>
>
>
>
> Thank you very much for your help or suggestions in advance!
>
> Gabriella Csima
>
> csim...@met.hu
>
>
>
>        [[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
Cincinnati, OH
+1 513 646 9390

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.


[R] sort a data.frame

2010-05-20 Thread Yuan Jian
Hello,
 
I have a dataframe:
dd <- data.frame(b = c("chr2", "chr1", "chr15", "chr13"),  
  x = c("A", "D", "A", "C"), y = c(8, 3, 9, 9), 
   z = c(1, 1, 1, 2)) 
 
>dd
  b x y z
1  chr2 A 8 1
2  chr1 D 3 1
3 chr15 A 9 1
4 chr13 C 9 2

Now I want to sort them according column "b", but only its number is considered:
  b x y z
1  chr1 D 3 1
2 chr13 C 9 2
3 chr15 A 9 1
4  chr2 A 8 1

thanks
jian
 
 
 


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

2010-05-20 Thread GenTOstats

Hey all

I am relatively new to analysis in R 
(and am a geneticist so apologies for rubbish questions ),

I have a question regarding repeated measures analysis

this was what i have done so far 
1)ph<-make.rm(constant=c("sid","pheno1","snp","sex","age"),repeated=c("measure1","measure2"),data=file)

2) i want to analyze 
 aov( repdat ~ pheno1 + snp + pheno1:snp + repdat:snp + repdat:pheno1 +
repdat:pheno1:snp + age + sex )

when i use SPSS (GLM repeated measures) it automatically generates the above
model for me)

But i read that I need to add in the sample ID term together with the Error
term to get the correct sum of squares 

3) Does this mean i need 

repdat ~ pheno1 + snp + pheno1:snp + repdat:snp + repdat:pheno1 +
repdat:pheno1:snp + age + sex 
 + sid + Error(sid/repdat + pheno1 + snp)

to get the correct anova summary?

Thanks

GS



-- 
View this message in context: 
http://r.789695.n4.nabble.com/Repeated-measures-in-R-tp2224326p2224326.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] Re : Manipulating Data Frames

2010-05-20 Thread David Winsemius


On May 20, 2010, at 6:38 AM, Mohan L wrote:


Dear All,

I have data some thisng like this :


data <- read.csv(file='ipsample.csv',sep=',' , header=TRUE)



data

 State  Jan  Feb  Mar  Apr  May Jun
1   AAA11022   0
2   BBB 1298 1195 1212 1244 1158 845
3  CCC 00012   1
4   DDD5   11   17   15   10   9
5   EEE   18   28   27   23   23  16
6   FFF   68  152  184  135  111  86

from this data frame, I took "Jan" as base and calculating weightage  
like

this :


basemonth.sum <- sum(data[[2]])



basemonth.sum

[1] 1390


basemonth.data <- data[[2]]



basemonth.data

[1]1 129805   18   68


weightage <- basemonth.data / basemonth.sum



weightage

[1] 0.0007194245 0.9338129496 0.00 0.0035971223 0.0129496403
[6] 0.0489208633


The above is the weightage for base month "Jan".  Now I need to  
calculate

weighted states data. What I need to do is :
(((Feb[i]-Jan[1])*weightage)+Jan[1]) for all column. The "Jan"  
column is

fixed. I need to do the calculation in all the column Feb, Mar etc...


data[, 3:7]*(
data[ , 2]/sum(data[ , 2]) )

Gives the reweighted estimates. You could easily cbind them to data[ ,  
1:2]





StateJan   Feb
   Mar
1   AAA1(((Feb[1]-Jan[1])*weightage[1])+Jan[1])
(((Mar[1]-Jan[1])*weightage[1])+Jan[1])
2   BBB 1298 (((Feb[2]-Jan[2])*weightage[2])+Jan[2])
(((Mar[2]-Jan[1])*weightage[2])+Ja[1])
3   CCC 0
4   DDD5
5   EEE   18
6   FFF   68

I am struggling with this . I have framed a logic using for loop.  
But it

seems me very bad logic.  Any help will be greatly appreciated.

Thanks & Rg
Mohan L

[[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] use object within rda file in for loop

2010-05-20 Thread Ivan Calandra

I've found the answer:
get() is exactly what I need, I completely forgot about this function.

There might be a better way, but that works for me.
Ivan

Le 5/20/2010 14:13, Ivan Calandra a écrit :

Dear users,

I would like to process all the lists from all *.rda files that I have 
in one folder.

Up to now, I can load all the *.rda files without any problem.

The problem is when I want to access the list saved within each *.rda 
file (only one list per rda file).


Here is my code:
fpath <- "D:/R"
listnames <- list.files(path=fpath, pattern=glob2rx("*.rda"), 
full.names=FALSE)

for (i in 1:length(listnames)) {
  load(paste(fpath, listnames[i], sep="/"))
  z <- list_in_listnames[i]  ## here is my problem
  **do something on z**
}

It might be really simple, but listnames is a character vector and I 
cannot find how to store the values within each element into z.
I think there would be a function to get and use an object with a 
given pattern in its name. The pattern itself is no problem, the 
problem is the function.
I've tried to look on RSiteSearch but, probably because I couldn't 
figure out the keywords, haven't found anything helpful.

Maybe I'm just on the wrong path to do it or missed something obvious...

I hope my question is clear. If not, please let me know what would 
help you to understand.


Thanks in advance
Ivan



--
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php

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


Re: [R] col allocation is not right

2010-05-20 Thread jim holtman
Look at the RColorBrewer package to understand how you can create the
colors you want.  Also do "?colorRampPalette"

On Wed, May 19, 2010 at 4:24 PM, Changbin Du  wrote:
> Thanks, Phil!
>
> Does that mean only eight colors can be used in R plot?
> THe following codes works for me, but, if I use the number, it does not
> work.
>
> plot(svm.auc, col=2, main="ROC curves comparing classification performance\n
> of six machine learning models")
> legend(0.5, 0.6, c(ns, nb, nr, nt, nl,ne), c(2,3,4,5,6,"black")) # Draw a
> legend.
>
> plot(bo.auc, col=3, add=T) # add=TRUE draws on the existing chart
> plot(rf.auc, col=4, add=T)
> plot(tree.auc, col=5, add=T)
> plot(nn.auc, col=6, add=T)
> plot(en.auc, col="black",lty=4,lwd=3, add=T)
>
>
>
> On Wed, May 19, 2010 at 11:30 AM, Phil Spector 
> wrote:
>
>> Changbin -
>>   Please take a look at the help file for the function
>> "palette", which is how R maps col= numbers to colors.
>> Also look at the default output of that function:
>>
>>  palette()
>>>
>> [1] "black"   "red"     "green3"  "blue"    "cyan"    "magenta" "yellow"
>> [8] "gray"
>>
>> You might get better results by specifying the colors that you
>> want directly.
>>
>>                                        - Phil Spector
>>                                         Statistical Computing Facility
>>                                         Department of Statistics
>>                                         UC Berkeley
>>                                         spec...@stat.berkeley.edu
>>
>>
>>
>> On Wed, 19 May 2010, Changbin Du wrote:
>>
>>  plot(svm.auc, col=2, main="ROC curves comparing classification
>>> performance\n
>>> of six machine learning models")
>>> legend(0.5, 0.6, c(ns, nb, nr, nt, nl,ne), 2:6, 9) # Draw a legend.
>>>
>>> plot(bo.auc, col=3, add=T) # add=TRUE draws on the existing chart
>>> plot(rf.auc, col=4, add=T)
>>> plot(tree.auc, col=5, add=T)
>>> plot(nn.auc, col=6, add=T)
>>> plot(en.auc, col=9,lty="dotted",lwd=3, add=T)
>>>
>>> Hi, Dear community,
>>>
>>> I am use the above codes to draw plot, but find that col =9 is not used by
>>> the R.  Instead, it use the col=2 when plot en.auc. WHy this happens and
>>> how
>>> to check the col allocation in R?
>>>
>>> Thanks!
>>>
>>>
>>> --
>>> Sincerely,
>>> Changbin
>>> --
>>>
>>>        [[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.
>>>
>>>
>
>
> --
> Sincerely,
> Changbin
> --
>
>        [[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
Cincinnati, OH
+1 513 646 9390

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] sort a data.frame

2010-05-20 Thread jim holtman
> dd
  b x y z
1  chr2 A 8 1
2  chr1 D 3 1
3 chr15 A 9 1
4 chr13 C 9 2
> # add column with just numbers
> dd$sort <- as.integer(gsub("\\D+", "", dd$b))
> dd[order(dd$sort),]  # notice it is a numeric, not character order
  b x y z sort
2  chr1 D 3 11
1  chr2 A 8 12
4 chr13 C 9 2   13
3 chr15 A 9 1   15
>


On Thu, May 20, 2010 at 8:28 AM, Yuan Jian  wrote:
> Hello,
>
> I have a dataframe:
> dd <- data.frame(b = c("chr2", "chr1", "chr15", "chr13"),
>   x = c("A", "D", "A", "C"), y = c(8, 3, 9, 9),
>    z = c(1, 1, 1, 2))
>
>>dd
>   b x y z
> 1  chr2 A 8 1
> 2  chr1 D 3 1
> 3 chr15 A 9 1
> 4 chr13 C 9 2
>
> Now I want to sort them according column "b", but only its number is 
> considered:
>   b x y z
> 1  chr1 D 3 1
> 2 chr13 C 9 2
> 3 chr15 A 9 1
> 4  chr2 A 8 1
>
> thanks
> jian
>
>
>
>
>
>
>        [[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
Cincinnati, OH
+1 513 646 9390

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] writing autocorrelation and partial auto correlation functions to a file

2010-05-20 Thread David Winsemius


On May 20, 2010, at 3:38 AM, nuncio m wrote:


Dear All,
 I am very new to T.  I need to fit a ARIMA model to my time
series.  So I found the auto correlation functions and partial auto
correlation function in R. Now I want to save these valuse along  
with the
significance levels to a file.  How to do that?.  I tried some  
function in R
like write.table but returns an error  "cannot coerce class "acf"  
into a

data.frame".  The I tried "write"
but returned again an error "Error in cat(list(...), file, sep, fill,
labels, append) :
 argument 1 (type 'list') cannot be handled by 'cat'".
Where am I going wrong


If you are just interested in getting the same information that you  
see on the console into a text file then:


?capture.output

--

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] Overlap of leaf labels

2010-05-20 Thread David Winsemius


On May 20, 2010, at 8:17 AM, Ivan Calandra wrote:


Hi,
I think that one of the cex arguments in par() can be what you're  
looking for. But since I've never plotted any dendrogram, I don't  
know which one, if any.


?par

Appears that the first effort should be to use cex.lab = 0.3 or some  
such.




HTH,
Ivan

Le 5/20/2010 14:08, Ayesha Jadoon a écrit :

Hi,

I have tried looking at the archives but havent found any answer  
that works

till now (Sorry if i have missed anything)


I am a newbie to R and i am trying to carry out hierarchical  
clustering
using hclust ->  as.dendrogram and then plotting the results as a  
dendrogram

using the plot function plot(object).

My question is :

In the function "plot", can one decrease the leaf label size to  
make them

readable and clear? I am including over 380 proteins in my dendrogram
and each leaf has a label which are currently overlapping and not
decipherable?


Thanks


Ayesha

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




--
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php

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


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] sort a data.frame

2010-05-20 Thread Nikhil Kaza

Try this.

dd[order(gsub("chr","",dd$b)),]

You need regular expressions if chr is not the only characterstring  
that is prepended to the numbers.

look for
?strsplit


Nikhil Kaza
University of North Carolina
nikhil.l...@gmail.com



On May 20, 2010, at 8:28 AM, Yuan Jian wrote:


Hello,

I have a dataframe:
dd <- data.frame(b = c("chr2", "chr1", "chr15", "chr13"),
  x = c("A", "D", "A", "C"), y = c(8, 3, 9, 9),
   z = c(1, 1, 1, 2))


dd

  b x y z
1  chr2 A 8 1
2  chr1 D 3 1
3 chr15 A 9 1
4 chr13 C 9 2

Now I want to sort them according column "b", but only its number is  
considered:

  b x y z
1  chr1 D 3 1
2 chr13 C 9 2
3 chr15 A 9 1
4  chr2 A 8 1

thanks
jian






[[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] sort a data.frame

2010-05-20 Thread Mohamed Lajnef

The solution given by Jim is more correct

dd[order(as.numeric(substr(dd$b,4,5))),]

regards
M
jim holtman a écrit :

dd


  b x y z
1  chr2 A 8 1
2  chr1 D 3 1
3 chr15 A 9 1
4 chr13 C 9 2
  

# add column with just numbers
dd$sort <- as.integer(gsub("\\D+", "", dd$b))
dd[order(dd$sort),]  # notice it is a numeric, not character order


  b x y z sort
2  chr1 D 3 11
1  chr2 A 8 12
4 chr13 C 9 2   13
3 chr15 A 9 1   15
  



On Thu, May 20, 2010 at 8:28 AM, Yuan Jian  wrote:
  

Hello,

I have a dataframe:
dd <- data.frame(b = c("chr2", "chr1", "chr15", "chr13"),
  x = c("A", "D", "A", "C"), y = c(8, 3, 9, 9),
   z = c(1, 1, 1, 2))



dd
  

  b x y z
1  chr2 A 8 1
2  chr1 D 3 1
3 chr15 A 9 1
4 chr13 C 9 2

Now I want to sort them according column "b", but only its number is considered:
  b x y z
1  chr1 D 3 1
2 chr13 C 9 2
3 chr15 A 9 1
4  chr2 A 8 1

thanks
jian






   [[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] Overlap of leaf labels

2010-05-20 Thread Ivan Calandra

Hi again,

Have you also tried cex.lab? And values <0.8, maybe 0.5 or even 0.3?

Output as Pdf/Ps shouldn't be a problem, or is it?

I cannot help you more without a reproducible example (read the posting 
guide, it's linked at the bottom of every email from the list), though 
even with it, I'm not sure I would be able to help you (but others would!)


HTH,
Ivan

PS: for more answers, don't write offlist!

Le 5/20/2010 15:14, ayesha.2.jad...@googlemail.com a écrit :

I have tried it but it didn't work! It doesn't decrease the label size 
significantly when I put

Cex=0.8 as an argument for plot()

I am hoping to export the output to pdf/ps format for visualisation as they 
allow zooming in and out!

I am really stuck! :(

Thanks!

Ayesha Jadoon

Kings College London
FWB, SE19NH


Sent using BlackBerry® from Orange

-Original Message-
From: Ivan Calandra
Date: Thu, 20 May 2010 14:17:56
To:
Subject: Re: [R] Overlap of leaf labels

Hi,
I think that one of the cex arguments in par() can be what you're
looking for. But since I've never plotted any dendrogram, I don't know
which one, if any.
HTH,
Ivan

Le 5/20/2010 14:08, Ayesha Jadoon a écrit :
   

Hi,

I have tried looking at the archives but havent found any answer that works
till now (Sorry if i have missed anything)


I am a newbie to R and i am trying to carry out hierarchical clustering
using hclust ->   as.dendrogram and then plotting the results as a dendrogram
using the plot function plot(object).

My question is :

In the function "plot", can one decrease the leaf label size to make them
readable and clear? I am including over 380 proteins in my dendrogram
and each leaf has a label which are currently overlapping and not
decipherable?


Thanks


Ayesha

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


 
   


--
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php

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


Re: [R] Re : Re : Nomogram with multiple interactions (package rms)

2010-05-20 Thread Frank E Harrell Jr

On 05/20/2010 01:42 AM, Marc Carpentier wrote:

Thank you for your responses, but I don't think you're right about the doc...
I carefully looked at it before posting and ran the examples, looked in 
Vanderbilt Biostat doc, and just looked again example(nomogram) :
1st example : categorical*continous : two axes for each sex
f<- lrm(y ~ lsp(age,50)+sex*rcs(cholesterol,4)+blood.pressure)


Hi Marc,

My apologies; I misread my own example.  This will take some digging 
into the code.  If you have time to do this before I do, code change 
suggestions welcomed.


Frank




2nd : continous*continous : one "age" axe for each specified value of 
cholesterol
g<- lrm(y ~ sex + rcs(age,3)*rcs(cholesterol,3))

3rd and 4th : categorical*continous : two axes for each sex (4th with fun)
f<- psm(Surv(d.time,death) ~ sex*age, dist='lognormal')

5th : categorical*continous : two axes for each sex (with fun)
g<- lrm(Y ~ age+rcs(cholesterol,4)*sex)

I'm desperately trying to represent a case of categorical*(continous+continous) 
:
f2<- cph(Surv(d.time,death) ~ sex*(rcs(cholesterol,4)+blood.pressure)
The best solution I can think of is to draw one nomogram for each sex :
Assuming 'male' is the ref level of sex :
1st nomogram : one axe for rcs(cholesterol,4), one axe for blood.pressure
2nd nomogram : one axe for sex:rcs(cholesterol,4), one axe for 
sex:blood.pressure, both shifted because of the sex own effect.
(I badly draw it in my previous mail)
I didn't see any example of this "adjustement" of nomogram to 'male' or 
'female'...

I hope I gave a clearer explanation and I'm not wrong about this unmentioned 
case.

Marc




- Message d'origine 
De : Frank E Harrell Jr
À : Marc Carpentier
Cc : r-help-request Mailing List
Envoyé le : Jeu 20 mai 2010, 0h 55min 32s
Objet : Re: Re : [R] Nomogram with multiple interactions (package rms)

On 05/19/2010 04:36 PM, Marc Carpentier wrote:

I'm sorry. I don't understand the "omit" solution, and maybe I mislead you with 
my explanation.

With the data from the "f" exemple of nomogram() :
Let's declare :
f2<- cph(Surv(d.time,death) ~ sex*(age+blood.pressure))
I guess the best (and maybe the only) way to represent it with a nomogram is to 
plot two nomograms (I couldn't find better).
Is there a way to have :

Nomogram1 : "male" :
- points 1-100 ---
- age ("men") ---
- blood.pressure ("men") ---
- linear predictor ---

And nomogram2 : "female" :
- points 1-100 ---
- age ("female") ---
- blood.pressure ("female") ---
- linear predictor ---

As I said I tried and failed (nomogram() still wants me to define
interact=list(...)) with :
plot(nomorgam(f2, adj.to=list(sex="male")) #and "female" for the other one

Marc


I think the documentation tells you how to do this.  But you failed to
look at the output from example(nomogram).  In one of the examples two
continuous predictors have two axes each, with male and female in close
proximity.  Or maybe I'm just missing your point.

Frank





- Message d'origine 
De : Frank E Harrell Jr
À : Marc Carpentier; r-help-request Mailing 
List
Envoyé le : Mer 19 mai 2010, 22h 28min 51s
Objet : Re: [R] Nomogram with multiple interactions (package rms)

On 05/19/2010 03:17 PM, Marc Carpentier wrote:

Dear list, I'm facing the following problem : A cox model with my sex
variable interacting with several continuous variables :
cph(S~sex*(x1+x2+x3)) And I'd like to make a nomogram. I know it's a
bit tricky and one mights argue that nomogram is not a good a
choice... I could use the parameter
interact=list(sex=("male","female"),x1=c(a,b,c))... but with rcs or
pol transformations of x1, x2 and x3, the choice of the
categorization (a,b,c,...) is arbitrary and the nomogram not so
useful... Considering that sex is the problem, I thought I could draw
two nomograms, one for each sex... based on one model. These would be
great. Do you think it's possible ?


Yes, you can specify constant predictors not to draw with the omit=
argument.  But try first to draw everything.  Shouldn't you just get 2
axes each for x1 x2 x3?



Taking the exemple of the help of nomogram() (package "rms") : f<-
psm(Surv(d.time,death) ~ sex*age, dist=if(.R.)'lognormal' else
'gaussian')


Drop the if(.R.) which was just corrected in the documentation.  Use
dist='lognormal'

Frank



Let's add the previously defined blood.pressure effect with an
interaction with sex too (with cph) : f2<- cph(Surv(d.time,death) ~
sex*(age+blood.pressure))

I thought of the parameter adt.to : plot(nomorgam(f2,
adj.to=list(sex="male")) #and "female" for the other one

But nomogram() still wants me to define interact=list(...) Thanks for
any advice you might have (with adj.to or any alternative...)

Marc Carpentier










--
Frank E Harrell Jr   Professor and ChairmanSchool of Medicine
 Department of Biostatistics   Vanderbilt University

__

Re: [R] sort a data.frame

2010-05-20 Thread Jorge Ivan Velez
Hi Yuan,

Try

dd[order(as.numeric(gsub("[^0-9]", "", dd$b))), ]

HTH,
Jorge

On Thu, May 20, 2010 at 8:28 AM, Yuan Jian <> wrote:

> Hello,
>
> I have a dataframe:
> dd <- data.frame(b = c("chr2", "chr1", "chr15", "chr13"),
>   x = c("A", "D", "A", "C"), y = c(8, 3, 9, 9),
>z = c(1, 1, 1, 2))
>
> >dd
>   b x y z
> 1  chr2 A 8 1
> 2  chr1 D 3 1
> 3 chr15 A 9 1
> 4 chr13 C 9 2
>
> Now I want to sort them according column "b", but only its number is
> considered:
>   b x y z
> 1  chr1 D 3 1
> 2 chr13 C 9 2
> 3 chr15 A 9 1
> 4  chr2 A 8 1
>
> thanks
> jian
>
>
>
>
>
>
>[[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] Minimization problem

2010-05-20 Thread Ravi Varadhan
That would not work in "optim".  The objective function must take in a
vector of parameters and return a scalar value.  Your objective function
does not return a scalar.  

This would work:

test <- function(x, ...){
m.error <- mean(distribution(x)) - observed.mean
sd.error <- std(distribution(x)) - observed.std
res <- c(m.error, sd.error)
return(sum(res^2))  # returns a scalar
}

Ravi.

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Bernardo Rangel Tura
Sent: Thursday, May 20, 2010 6:02 AM
To: r-help@r-project.org
Subject: Re: [R] Minimization problem

On Thu, 2010-05-20 at 01:35 -0700, Fred wrote:
> Dear R users,
> 
> I am trying to minimize two function simultaneously in R,
> 
> function(x)
> 
> minimize x[1],x[2],x[3]
> 
> mean(distribution(x1,x2,x3) ) - observed mean
> 
> std(distribution(x1,x2,x3)) - observed std
> 
> What I want to achieve is that simulated mean and standard deviation
> of distribution related to x1 x2 x3  would be close to observed mean
> and observed standard deviation.
> 
> 
> is there any function in R can reach this?
> 
> Thank you for the help first .
> 
> F.

Hi!
Do you need use optim, something like this

test <- function(parameters){
m.error <- mean(distribution(x1,x2,x3) ) - observed mean
m.sd <- std(distribution(x1,x2,x3)) - observed std
res <- cbind(m.error,sd.error)
return(res)
}
-- 
Bernardo Rangel Tura, M.D,MPH,Ph.D
National Institute of Cardiology
Brazil

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

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


Re: [R] Overlap of leaf labels

2010-05-20 Thread David Winsemius
I got an offlist response saying my advice was not correct, not  
surprising in the absence of a reproducible example on which to test.  
Looking at the help page, which seems the sensible place for anyone to  
start, we see that the second plotting example for dendrograms uses  
lab.cex and it need to be in a list offered to nodePar.


plot(dend1, nodePar = list(lab.cex = 0.3))

So I think my advice was correct up to the wider interpretation of  
"... or some such."


--
David.


On May 20, 2010, at 9:15 AM, David Winsemius wrote:



On May 20, 2010, at 8:17 AM, Ivan Calandra wrote:


Hi,
I think that one of the cex arguments in par() can be what you're  
looking for. But since I've never plotted any dendrogram, I don't  
know which one, if any.


?par

Appears that the first effort should be to use cex.lab = 0.3 or some  
such.




HTH,
Ivan

Le 5/20/2010 14:08, Ayesha Jadoon a écrit :

Hi,

I have tried looking at the archives but havent found any answer  
that works

till now (Sorry if i have missed anything)


I am a newbie to R and i am trying to carry out hierarchical  
clustering
using hclust ->  as.dendrogram and then plotting the results as a  
dendrogram

using the plot function plot(object).

My question is :

In the function "plot", can one decrease the leaf label size to  
make them
readable and clear? I am including over 380 proteins in my  
dendrogram

and each leaf has a label which are currently overlapping and not
decipherable?


Thanks


Ayesha

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




--
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php

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


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.


David Winsemius, MD
West Hartford, CT

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


Re: [R] use object within rda file in for loop

2010-05-20 Thread Henrik Bengtsson
On Thu, May 20, 2010 at 3:06 PM, Ivan Calandra
 wrote:
> I've found the answer:
> get() is exactly what I need, I completely forgot about this function.
>
> There might be a better way, but that works for me.

I R.utils, there are saveObject() and loadObject(), which allows you
to save and load objects without having to worry about their names.
For example:

library("R.utils");
x <- list(a=2, b=3:40, fcn=rnorm);
saveObject(x, "objs.Rbin");

Then later you can load the list object that 'x' held as whatever you
want, say, 'y', e.g.

library("R.utils");
y <- loadObject("objs.Rbin");

This way you don't "contaminate" the working environment with
variables named according to the file, which sometimes can be a
surprise (since you don't always know what the file contain).  Except
from that it works just like save()/load().

/Henrik

> Ivan
>
> Le 5/20/2010 14:13, Ivan Calandra a écrit :
>>
>> Dear users,
>>
>> I would like to process all the lists from all *.rda files that I have in
>> one folder.
>> Up to now, I can load all the *.rda files without any problem.
>>
>> The problem is when I want to access the list saved within each *.rda file
>> (only one list per rda file).
>>
>> Here is my code:
>> fpath <- "D:/R"
>> listnames <- list.files(path=fpath, pattern=glob2rx("*.rda"),
>> full.names=FALSE)
>> for (i in 1:length(listnames)) {
>>  load(paste(fpath, listnames[i], sep="/"))
>>  z <- list_in_listnames[i]  ## here is my problem
>>  **do something on z**
>> }
>>
>> It might be really simple, but listnames is a character vector and I
>> cannot find how to store the values within each element into z.
>> I think there would be a function to get and use an object with a given
>> pattern in its name. The pattern itself is no problem, the problem is the
>> function.
>> I've tried to look on RSiteSearch but, probably because I couldn't figure
>> out the keywords, haven't found anything helpful.
>> Maybe I'm just on the wrong path to do it or missed something obvious...
>>
>> I hope my question is clear. If not, please let me know what would help
>> you to understand.
>>
>> Thanks in advance
>> Ivan
>>
>
> --
> Ivan CALANDRA
> PhD Student
> University of Hamburg
> Biozentrum Grindel und Zoologisches Museum
> Abt. Säugetiere
> Martin-Luther-King-Platz 3
> D-20146 Hamburg, GERMANY
> +49(0)40 42838 6231
> ivan.calan...@uni-hamburg.de
>
> **
> http://www.for771.uni-bonn.de
> http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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


Re: [R] Overlap of leaf labels

2010-05-20 Thread Ayesha Jadoon
Hi

Sorry to be honest, i hit reply before i realised!


Thanks a bunch!
It worked! I can visualise the labels!

Ayesha Jadoon

PhD student
PSD, FWB
Kings college London
SE1 9NH
UK
On Thu, May 20, 2010 at 2:38 PM, David Winsemius wrote:

> I got an offlist response saying my advice was not correct, not surprising
> in the absence of a reproducible example on which to test. Looking at the
> help page, which seems the sensible place for anyone to start, we see that
> the second plotting example for dendrograms uses lab.cex and it need to be
> in a list offered to nodePar.
>
> plot(dend1, nodePar = list(lab.cex = 0.3))
>
> So I think my advice was correct up to the wider interpretation of "... or
> some such."
>
> --
> David.
>
>
>
> On May 20, 2010, at 9:15 AM, David Winsemius wrote:
>
>
>> On May 20, 2010, at 8:17 AM, Ivan Calandra wrote:
>>
>> Hi,
>>> I think that one of the cex arguments in par() can be what you're looking
>>> for. But since I've never plotted any dendrogram, I don't know which one, if
>>> any.
>>>
>>
>> ?par
>>
>> Appears that the first effort should be to use cex.lab = 0.3 or some such.
>>
>>
>> HTH,
>>> Ivan
>>>
>>> Le 5/20/2010 14:08, Ayesha Jadoon a écrit :
>>>
 Hi,

 I have tried looking at the archives but havent found any answer that
 works
 till now (Sorry if i have missed anything)


 I am a newbie to R and i am trying to carry out hierarchical clustering
 using hclust ->  as.dendrogram and then plotting the results as a
 dendrogram
 using the plot function plot(object).

 My question is :

 In the function "plot", can one decrease the leaf label size to make
 them
 readable and clear? I am including over 380 proteins in my dendrogram
 and each leaf has a label which are currently overlapping and not
 decipherable?


 Thanks


 Ayesha

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



>>> --
>>> Ivan CALANDRA
>>> PhD Student
>>> University of Hamburg
>>> Biozentrum Grindel und Zoologisches Museum
>>> Abt. Säugetiere
>>> Martin-Luther-King-Platz 3
>>> D-20146 Hamburg, GERMANY
>>> +49(0)40 42838 6231
>>> ivan.calan...@uni-hamburg.de
>>>
>>> **
>>> http://www.for771.uni-bonn.de
>>> http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php
>>>
>>> __
>>> R-help@r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>
>> 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.
>>
>
> 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.
>

[[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] use object within rda file in for loop

2010-05-20 Thread Ivan Calandra

Hi Henrik,

Thank you for this suggestion, it does sound easier indeed!

Ivan

Le 5/20/2010 15:50, Henrik Bengtsson a écrit :

On Thu, May 20, 2010 at 3:06 PM, Ivan Calandra
  wrote:
   

I've found the answer:
get() is exactly what I need, I completely forgot about this function.

There might be a better way, but that works for me.
 

I R.utils, there are saveObject() and loadObject(), which allows you
to save and load objects without having to worry about their names.
For example:

library("R.utils");
x<- list(a=2, b=3:40, fcn=rnorm);
saveObject(x, "objs.Rbin");

Then later you can load the list object that 'x' held as whatever you
want, say, 'y', e.g.

library("R.utils");
y<- loadObject("objs.Rbin");

This way you don't "contaminate" the working environment with
variables named according to the file, which sometimes can be a
surprise (since you don't always know what the file contain).  Except
from that it works just like save()/load().

/Henrik

   

Ivan

Le 5/20/2010 14:13, Ivan Calandra a écrit :
 

Dear users,

I would like to process all the lists from all *.rda files that I have in
one folder.
Up to now, I can load all the *.rda files without any problem.

The problem is when I want to access the list saved within each *.rda file
(only one list per rda file).

Here is my code:
fpath<- "D:/R"
listnames<- list.files(path=fpath, pattern=glob2rx("*.rda"),
full.names=FALSE)
for (i in 1:length(listnames)) {
  load(paste(fpath, listnames[i], sep="/"))
  z<- list_in_listnames[i]  ## here is my problem
  **do something on z**
}

It might be really simple, but listnames is a character vector and I
cannot find how to store the values within each element into z.
I think there would be a function to get and use an object with a given
pattern in its name. The pattern itself is no problem, the problem is the
function.
I've tried to look on RSiteSearch but, probably because I couldn't figure
out the keywords, haven't found anything helpful.
Maybe I'm just on the wrong path to do it or missed something obvious...

I hope my question is clear. If not, please let me know what would help
you to understand.

Thanks in advance
Ivan

   

--
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php

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

 
   


--
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php

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


[R] finding euclidean proximate points in two datasets

2010-05-20 Thread Alexander Shenkin
Hello all,

I've been pouring through the various spatial packages, but haven't come
across the right thing yet.

Given a set of points in 2-d space X, i'm trying to find the subset of
points in Y proximate to each point in X.  Furthermore, the proximity
threshold of each point in X differs (X$threshold).  I've constructed
this myself already, but it's horrificly slow with a dataset of 40k+
points in one set, and a 700 in the other.

A very inefficient example of what I'm looking for:

for (pt in X$idx) {
proximity[i] = euclidian_dist(X[pt]$x, X[pt]$y, Y$x, Y$y) <
X$threshold
i = i+1
}

Perhaps crossdist() in spatstat is what I should use, and then code a
comparison with X$threshold after the cross-distances are computed.
However, I was wondering if there was another tool I should be
considering.  Any and all thoughts are very welcome.  Thanks in advance.

Thanks,
Allie
-- 
Alexander Shenkin
PhD Candidate
School of Natural Resources and Environment
University of Florida

http://snre.ufl.edu/people/students.asp

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


Re: [R] finding euclidean proximate points in two datasets

2010-05-20 Thread David Winsemius


On May 20, 2010, at 10:02 AM, Alexander Shenkin wrote:


Hello all,

I've been pouring through the various spatial packages, but haven't  
come

across the right thing yet.


There is a SIG for such questions.



Given a set of points in 2-d space X, i'm trying to find the subset of
points in Y proximate to each point in X.  Furthermore, the proximity
threshold of each point in X differs (X$threshold).  I've constructed
this myself already, but it's horrificly slow with a dataset of 40k+
points in one set, and a 700 in the other.

A very inefficient example of what I'm looking for:


Not really a reproducible example. If euclidean_dist is a function ,  
then it is not one in any of the packages I have installed.




   for (pt in X$idx) {
   proximity[i] = euclidian_dist(X[pt]$x, X[pt]$y, Y$x, Y$y) <
X$threshold
i = i+1
   }



Have you considered first creating a subset of candidate points that  
are within "threshold" of each reference point on both coordinates.  
That might sidestep a lot of calculations on points that are easily  
eliminated on a single comparison. Then you could calculate distances  
within that surviving subset of points. On average that should give  
you an over 50% "hit rate":


> (4/3)*pi*0.5^3
[1] 0.5235988




Perhaps crossdist() in spatstat is what I should use, and then code a
comparison with X$threshold after the cross-distances are computed.
However, I was wondering if there was another tool I should be
considering.  Any and all thoughts are very welcome.  Thanks in  
advance.


Thanks,
Allie
--
Alexander Shenkin
PhD Candidate
School of Natural Resources and Environment
University of Florida

--
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] [Off topic?] Time dependent Cox model fitting and validation

2010-05-20 Thread Jabba

DeaR users.


These days i'm working on fitting an extended Cox model with
time-dependent covariables and possibly time-varying effects. My
data are in counting process format as described in Therneau&Grambsh's
`Modeling Survival Data', page 68. I'm trying to follow Harrell's
`Regression Modeling Strategies' advices for the choice of my  model.
This study aims to the development of a prognostic model, so it'is
primary predictive.

I have to do stepwise model selection and provide a measure of
predictive accuracy. I'm using rms's cph and validate function with
bw=TRUE option.



1. Is validate good at resampling from a counting process format
database? Or should i use a somewhat modified version?

2. Why fastbw(fit,"aic") and step(fit) don't select the
same model? step() appears to stop first. I can't manage to get the
stopping rule in the help files.

3. cph seems to be a bit less "permissive" than coxph in parsing the
model formula. Particularly i have some difficulty in modeling
interactions between covariables and time. Am I totally misguided? Is
there any reference on this topic?

Now a theoretical one:

4. Is it somewhat sensible to use cox.zph() and schoenfeld
residuals to investigate which time dependent variables could need a
time interaction parameter for estimating a time-varying effect?


Thanks in advance for any advice.

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

2010-05-20 Thread Albert-Jan Roskam
Hi,
 
Is there a built in function that returns a character vector of all subclasses 
of a given superclass? 
showClass(Class = "SomeClass") contains the info that I want, but I don't know 
how to access it.
getSubClasses <- function(superClass) return(setdiff(getClasses(.GlobalEnv), 
superClass)) wíll only work if the global enviroment isn;t filled with other 
class definitions.

I've seen functions that somewhat resemble it, but are not quite the same. 
Thanks 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] Time problems (POSIXct)

2010-05-20 Thread Research

Hello,

I have a zoo time series object (say a)  with the following time 
stamp/format




[1] "1950-01-03 GMT" "1950-01-04 GMT" "1950-01-05 GMT" "1950-01-06 GMT"
[5] "1950-01-09 GMT" "1950-01-10 GMT"


and another (say b) with


[1] "1950-01-05" "1950-01-06" "1950-01-09" "1950-01-10" "1950-01-11"
[6] "1950-01-12"


I want to merge these series   but when I try:


> head(merge(a, b ))
a   b
1950-01-03 02:00:00   16.66   NA
1950-01-04 02:00:00   16.85   NA
1950-01-05 02:00:00   16.93   NA
1950-01-06 02:00:00   16.98   NA
1950-01-09 02:00:00   17.08   NA
1950-01-10 02:00:00   17.03   NA
Warning message:
In merge.zoo(a, b) :
  Index vectors are of different classes: POSIXt Date
>

Anybody can help?  I tried as.POSIXct(a,tz="GMT") but it doesn't seem to 
work


Many thanks in advance!

Costas
  _
platform   i386-pc-mingw32
arch   i386
os mingw32
system i386, mingw32
status
major  2
minor  10.1
year   2009
month  12
day14
svn rev50720
language   R
version.string R version 2.10.1 (2009-12-14)
>

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


Re: [R] getSubClasses()?

2010-05-20 Thread Martin Morgan
On 05/20/2010 07:47 AM, Albert-Jan Roskam wrote:
> Hi,
> 
> Is there a built in function that returns a character vector of all
> subclasses of a given superclass? showClass(Class = "SomeClass")
> contains the info that I want, but I don't know how to access it. 

Hi Albert-Jan --

Maybe

> slotNames(class(getClass("vector")))
 [1] "slots"  "contains"   "virtual""prototype"  "validity"
 [6] "access" "className"  "package""subclasses" "versionKey"
[11] "sealed"
> names(getClass("vector")@subclasses)
 [1] "logical""numeric""character"
 [4] "complex""integer""raw"
 [7] "expression" "list"   "structure"
[10] "array"  "matrix" "signature"
[13] "ObjectsWithPackage" "mts""ordered"
[16] "namedList"  "listOfMethods"

Martin

> getSubClasses <- function(superClass)
> return(setdiff(getClasses(.GlobalEnv), superClass)) wíll only work if
> the global enviroment isn;t filled with other class definitions.
> 
> I've seen functions that somewhat resemble it, but are not quite the
> same. Thanks 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.


-- 
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793

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


Re: [R] Time problems (POSIXct)

2010-05-20 Thread Cedrick W. Johnson

try:

as.POSIXct(datestr, format="%y-%m-%d")

for the first timeseries...

-c


On 5/20/2010 10:48 AM, Research wrote:

Hello,

I have a zoo time series object (say a) with the following time
stamp/format



[1] "1950-01-03 GMT" "1950-01-04 GMT" "1950-01-05 GMT" "1950-01-06 GMT"
[5] "1950-01-09 GMT" "1950-01-10 GMT"


and another (say b) with


[1] "1950-01-05" "1950-01-06" "1950-01-09" "1950-01-10" "1950-01-11"
[6] "1950-01-12"


I want to merge these series but when I try:


 > head(merge(a, b ))
a b
1950-01-03 02:00:00 16.66 NA
1950-01-04 02:00:00 16.85 NA
1950-01-05 02:00:00 16.93 NA
1950-01-06 02:00:00 16.98 NA
1950-01-09 02:00:00 17.08 NA
1950-01-10 02:00:00 17.03 NA
Warning message:
In merge.zoo(a, b) :
Index vectors are of different classes: POSIXt Date
 >

Anybody can help? I tried as.POSIXct(a,tz="GMT") but it doesn't seem to
work

Many thanks in advance!

Costas
_
platform i386-pc-mingw32
arch i386
os mingw32
system i386, mingw32
status
major 2
minor 10.1
year 2009
month 12
day 14
svn rev 50720
language R
version.string R version 2.10.1 (2009-12-14)
 >

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Antwort: Re: Multiple language output - Correct in RGui, wrong in .txt after sink()

2010-05-20 Thread mark . redshaw
Dear Brian,
many thanks for reply it helped a great deal. But firstly an apology for 
not providing the "at a minimum",my (newbe) error, sorry for making you 
guess.
The info is/was:
> sessionInfo()
R version 2.10.0 (2009-10-26) 
i386-pc-mingw32 

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252 
LC_MONETARY=German_Germany.1252 LC_NUMERIC=C 
[5] LC_TIME=German_Germany.1252 

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

I am "stuck" using the system I have , no way for me to change that and 
still get access to the data I need. The email is from Lotus Notes and I 
also have no idea how this works with encoding but think it is responsible 
for the problems with the Korean text.
For reference the fonts I have used are MS Mincho and Arial Unicode MS 
which both seen to have all the characters I need in the languages I am 
using.

I have taken your suggestions and tested and had positive results.

With the connection that you suggested, I found that on my system that the 
solution was:
con <- file("output.txt",open = "a",encoding="UTF-8")

I also started looking the locale following you hint and found that by 
changing this together with the connection I could get the output in the 
form i wanted in all languages including Korean
I am, as you say, just thankful for the "miracle" of this working at all.
It may not be optimal, and ideas from you or others would be welcome, but 
it works!
many thanks
Mark

""
RM_EN <- c("Alfalfa hay","Alfalfa meal","Alfalfa silage")
RM_DE <- c("Luzerneheu","Lurzernegrünmehl","Luzernesilage")
RM_RU <- c("Люцерновое сено","Люцерновая 
травяная мука","Люцерновый 
сенаж")
RM_CN <- c("苜蓿干草","苜蓿草粉","苜蓿青贮")
RM_JP <- c("アルファルファ乾草","アルファルファ 
ミール","アルファルファ サ
イレージ")
RM_KR <- c("알팔파 건초","알팔파 박","알팔파 사일리지")

RMLANG <- data.frame(RM_EN,RM_DE,RM_RU,RM_CN,RM_JP,RM_KR)
nrm <- NROW(RMLANG)

con <- file("output.txt",open = "a",encoding="UTF-8")
for(i in 1:nrm)
{
cat("English", as.character(RMLANG$RM_EN[i]), file=con,"\n",sep="")
cat("German", as.character(RMLANG$RM_DE[i]), file=con,"\n",sep="")
Sys.setlocale("LC_ALL","Chinese_CHN")
cat("Chinese", as.character(RMLANG$RM_CN[i]), file=con,"\n",sep="")
Sys.setlocale("LC_ALL","Japanese")
cat("Japanese", as.character(RMLANG$RM_JP[i]), file=con,"\n",sep="")
Sys.setlocale("LC_ALL","Korean")
cat("Korean", as.character(RMLANG$RM_KR[i]), file=con,"\n",sep="")
Sys.setlocale("LC_ALL","German")
}
close(con)
""

Dr. Mark Redshaw 
Animal Nutrition Services 
Evonik Degussa GmbH, HN-M-AN, Rodenbacher Chaussee 4, 63457 Hanau, Germany 

Tel: +49 61 81 59 6788 
www.aminoacidsandmore.com 




Prof Brian Ripley  
19.05.2010 23:35

An
mark.reds...@evonik.com
Kopie
r-help@r-project.org
Thema
Re: [R] Multiple language output  - Correct in RGui, wrong in .txt after 
sink()






You haven't given us the 'at a minimum' information asked for in the 
the posting guide (but we can guess you are using Windows), nor do we 
know the intended encoding of this email (I see no encoding in the 
header as it reached me, but it seems sensible viewed as UTF-8). And 
the absence of basic information does make it *really* hard to help 
here -- this reply is my third guess at what might be happening.

We also do not know the font you are using in RGui, but I am 
not aware of any Windows font which covers correctly Russian and CJK.
However, it is not just a question of knowing the font name: different 
versions of Windows, including different language-specific versions, 
have different fonts with the same name.

RGui (since about R 2.7.0) works in UCS-2 encoding.  Sink files work 
in the locale's encoding (another of the pieces of information you did 
not tell us, but on Windows it is 8-bit or specific to one of 
Simplified Chinese, Traditional Chinese, Japanese or Korean -- I'd 
guess from your address it was CP1252, but it *is* part of the 'at a 
minimum').  So whereas R can store non-native strings in UTF-8 
(provided you get them in as such), it can only output them if told 
how to: the designer of RGui did so but you in using 
sink('output.txt') did not.

cat+sink is an inefficient way to write to a file: try using the file= 
argument on an opened connection.  And you can set the encoding on 
that connection.  I really don't know what you meant by 'the 
characters as I expect': in a file they have to be in *some* encoding 
and you are not looking at bits but as a representation in some 
unspecified file viewer.  One possibility is that you meant UCS-2 
(what Windows tends incorrectly to call 'Unicode' files), in which 
case you can use something like

con <- file("foo", encoding="UCS-2LE")
cat(..., file=con)
...
close(con)

You can use a connection with sink() too.

Think of it more as a miracle (and much unappreciated hard work and 
inspired des

Re: [R] Time problems (POSIXct)

2010-05-20 Thread Gabor Grothendieck
The warning message does tell you exactly what was wrong.   You are
trying to merge zoo objects that have two different index classes.
The first is "POSIXct" and the second is "Date" class.  Next time
please provide your data in reproducible form
as shown here:

> library(zoo)
> z1 <- zoo(1:6, as.Date(c("1950-01-05", "1950-01-06", "1950-01-09",
+   "1950-01-10", "1950-01-11", "1950-01-12")))
> z2 <- zoo(c(16.66, 16.85, 16.93, 16.98, 17.08, 17.03),
+ as.POSIXct(c("1950-01-03, GMT", "1950-01-04, GMT", "1950-01-05, GMT",
+   "1950-01-06, GMT", "1950-01-09 GMT", "1950-01-10 GMT")))
>
> merge(z1, z2 = aggregate(z2, as.Date, identity))
   z1z2
1950-01-03 NA 16.66
1950-01-04 NA 16.85
1950-01-05  1 16.93
1950-01-06  2 16.98
1950-01-09  3 17.08
1950-01-10  4 17.03
1950-01-11  5NA
1950-01-12  6NA
>
>
> # or
> z2.Date <- z2
> time(z2.Date) <- as.Date(time(z2.Date))
>
> merge(z1, z2.Date)
   z1 z2.Date
1950-01-03 NA   16.66
1950-01-04 NA   16.85
1950-01-05  1   16.93
1950-01-06  2   16.98
1950-01-09  3   17.08
1950-01-10  4   17.03
1950-01-11  5  NA
1950-01-12  6  NA


On Thu, May 20, 2010 at 10:48 AM, Research  wrote:
> Hello,
>
> I have a zoo time series object (say a)  with the following time
> stamp/format
>
>
>
> [1] "1950-01-03 GMT" "1950-01-04 GMT" "1950-01-05 GMT" "1950-01-06 GMT"
> [5] "1950-01-09 GMT" "1950-01-10 GMT"
>
>
> and another (say b) with
>
>
> [1] "1950-01-05" "1950-01-06" "1950-01-09" "1950-01-10" "1950-01-11"
> [6] "1950-01-12"
>
>
> I want to merge these series   but when I try:
>
>
>> head(merge(a, b ))
>                    a                   b
> 1950-01-03 02:00:00       16.66   NA
> 1950-01-04 02:00:00       16.85   NA
> 1950-01-05 02:00:00       16.93   NA
> 1950-01-06 02:00:00       16.98   NA
> 1950-01-09 02:00:00       17.08   NA
> 1950-01-10 02:00:00       17.03   NA
> Warning message:
> In merge.zoo(a, b) :
>  Index vectors are of different classes: POSIXt Date
>>
>
> Anybody can help?  I tried as.POSIXct(a,tz="GMT") but it doesn't seem to
> work
>
> Many thanks in advance!
>
> Costas
>              _
> platform       i386-pc-mingw32
> arch           i386
> os             mingw32
> system         i386, mingw32
> status
> major          2
> minor          10.1
> year           2009
> month          12
> day            14
> svn rev        50720
> language       R
> version.string R version 2.10.1 (2009-12-14)
>>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] getSubClasses()?

2010-05-20 Thread Albert-Jan Roskam
Hi Martin,
 
Great, the second solution is precisely what I've been looking for. Thanks a 
lot!

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?
~~

--- On Thu, 5/20/10, Martin Morgan  wrote:


From: Martin Morgan 
Subject: Re: [R] getSubClasses()?
To: "Albert-Jan Roskam" 
Cc: "R Mailing List" 
Date: Thursday, May 20, 2010, 4:58 PM


On 05/20/2010 07:47 AM, Albert-Jan Roskam wrote:
> Hi,
> 
> Is there a built in function that returns a character vector of all
> subclasses of a given superclass? showClass(Class = "SomeClass")
> contains the info that I want, but I don't know how to access it. 

Hi Albert-Jan --

Maybe

> slotNames(class(getClass("vector")))
[1] "slots"      "contains"   "virtual"    "prototype"  "validity"
[6] "access"     "className"  "package"    "subclasses" "versionKey"
[11] "sealed"
> names(getClass("vector")@subclasses)
[1] "logical"            "numeric"            "character"
[4] "complex"            "integer"            "raw"
[7] "expression"         "list"               "structure"
[10] "array"              "matrix"             "signature"
[13] "ObjectsWithPackage" "mts"                "ordered"
[16] "namedList"          "listOfMethods"

Martin

> getSubClasses <- function(superClass)
> return(setdiff(getClasses(.GlobalEnv), superClass)) wíll only work if
> the global enviroment isn;t filled with other class definitions.
> 
> I've seen functions that somewhat resemble it, but are not quite the
> same. Thanks 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.


-- 
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793



  
[[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] finding euclidean proximate points in two datasets

2010-05-20 Thread Alexander Shenkin
On 5/20/2010 9:18 AM, David Winsemius wrote:
> 
> On May 20, 2010, at 10:02 AM, Alexander Shenkin wrote:
> 
>> Hello all,
>>
>> I've been pouring through the various spatial packages, but haven't come
>> across the right thing yet.
> 
> There is a SIG for such questions.

thanks - just joined it.

>> Given a set of points in 2-d space X, i'm trying to find the subset of
>> points in Y proximate to each point in X.  Furthermore, the proximity
>> threshold of each point in X differs (X$threshold).  I've constructed
>> this myself already, but it's horrificly slow with a dataset of 40k+
>> points in one set, and a 700 in the other.
>>
>> A very inefficient example of what I'm looking for:
> 
> Not really a reproducible example. If euclidean_dist is a function ,
> then it is not one in any of the packages I have installed.

it's not reproducible - i'll make a better effort to include
reproducible code in the future.  and that function is just one i would
have written, but didn't want to clog the email with code.  Anyway, here
is a reproducible example:

X = data.frame(x=c(1,2,3), y=c(2,3,1), threshold=c(1,2,4))
Y = data.frame(x=c(5,2,3,4,2,5,2,3), y=c(5,2,2,4,1,2,3,1))
proximate=list()
i=1
for (pt in 1:length(X$x)) {
proximate[[i]] <- sqrt((X[pt,]$x - Y$x)^2 + (X[pt,]$y - Y$y)^2) >
X[pt,]$threshold
i=i+1
}
proximate

>>
>>for (pt in X$idx) {
>>proximity[i] = euclidian_dist(X[pt]$x, X[pt]$y, Y$x, Y$y) <
>> X$threshold
>> i = i+1
>>}
>>
> 
> Have you considered first creating a subset of candidate points that are
> within "threshold" of each reference point on both coordinates. That
> might sidestep a lot of calculations on points that are easily
> eliminated on a single comparison. Then you could calculate distances
> within that surviving subset of points. On average that should give you
> an over 50% "hit rate":
> 
>> (4/3)*pi*0.5^3
> [1] 0.5235988

That's a nice idea.  I'll still be waiting quite a while while my
machine cranks, but not as long.  Still - I suspect there would be much
bigger gains if there were tailored packages.  I'll re-post over on
sig-geo about that.  thanks.

>> Perhaps crossdist() in spatstat is what I should use, and then code a
>> comparison with X$threshold after the cross-distances are computed.
>> However, I was wondering if there was another tool I should be
>> considering.  Any and all thoughts are very welcome.  Thanks in advance.
>>
>> Thanks,
>> Allie
>> -- 
>> Alexander Shenkin
>> PhD Candidate
>> School of Natural Resources and Environment
>> University of Florida

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] esthetics --- extending the lm command to fixed effects?

2010-05-20 Thread ivo welch
dear R wizards:

not important.  more a curiosity or esthetics question.

is there a way to extend the standard lm command, so that it takes a new
argument that handles fixed effects?   right now, I have (provided to me
from an expert---I would have never figured this one out):

   diffid <- function(h,id) {
   id <- as.factor(id)[, drop=TRUE]
   apply(as.matrix(h), 2, function(x) x - tapply(x,id,mean)[id]
   }

which is used as

 r= lm( diffid(y, firmid) ~ diffid(x, firmid ) )

it works, but it would be much nicer if I could just write

r= lm( y ~ x + z, fixed.effects=firmid )

does this already exists as a package?  or has someone figured out how to
program this?

as I wrote---this is a curiosity question, not a substance question.

regards,

/iaw

Ivo Welch (ivo.we...@brown.edu, ivo.we...@gmail.com)

[[alternative HTML version deleted]]

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


[R] RODBC: owerwrite into a named range in Excel

2010-05-20 Thread Jay
Hello,

Let's say that I have a data frame of n numbers I want to transfer
into a Excel spreadsheet. I have opened the conection to the file
using ODBC, and I can query the content of these n cells without
problem. However, how do I transfer my new values to these cells?
I.e., overwite them.

Should I use sqlSave() or sqlUpdate()?

Using the update I get the error: "cannot update ‘data_001’ without
unique column"

sqlUpdate(connection_name, my_new_data_frame,
"name_of_the_range_in_excel")

Let's say that the range in Excel is in E10:E20. (if it matters)


BR,
Jay

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


Re: [R] offset in gam and spatial scale of variables

2010-05-20 Thread Lucia Rueda

Hi,

Thanks for the inputs. I talked to my coworker, who has been the one doing
the analysis. Perhaps I wasn't making myself clear about the “differences in
spatial scales”.  Here is what he says:

"The truth is that measuring scales (i.e all area related variable are
measured in m2) and spatial definition of initial cartography are
homogeneous among extracted variables. But all variables (ie. sum of the
total rocky bottom in the surrounding area) are computed for each different
integration areas (buffer) (i.e in an area of 40squaremeters around the
sample, in an area of 80m2, …).
The question is then if we can build a model that includes variables
measured at different buffers (for example a model that includes 3
variables:  1.-  the amount of rocky bottom in an area of 80m2 ; 2- the
amount of sandy bottom in an area of 200m2; and the mean depth calculated in
a surrounding area of 50m2) considering that each variable may be expressing
different ecological processes. I believe that if there is not an ecological
constrain in the interpretation of the variables (and their ecological
effect over the specie), including them in a model is correct, unless there
is not a mathematical constrain."

Also, about the spatial correlation I thought from what I've read so far
that I had to build the model and then check if there was spatial
correlation in the residuals since they are supposed to be i.i.d. And if it
turns out that they are then I have to do something about it like gamm, gee,
sar, car, etc.

Cheers,

Lucia


-- 
View this message in context: 
http://r.789695.n4.nabble.com/offset-in-gam-and-spatial-scale-of-variables-tp483p2224528.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] factor

2010-05-20 Thread RockO

Try, ofr a factor:

> x <-c("A","B","C")
> x <- factor(x)
> levels(x)
[1] "A" "B" "C"
> x <- factor(x,levels=levels(x)[c(3,2,1)])
> levels(x)
[1] "C" "B" "A"

Rock, DRF
-- 
View this message in context: 
http://r.789695.n4.nabble.com/factor-tp879983p2224639.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] Using svychisq inside user-defined function

2010-05-20 Thread Thomas Lumley

On Wed, 19 May 2010, Sabatier, Jennifer F. (CDC/OID/NCHHSTP) wrote:


Hi R-help,

Yes, this is my second request for assistance in a single day

I am attempting to use svychisq() inside a function  I made.  The goal
of this function is to produce a table of summary statistics that I can
later output to EXCEL (simple frequencies and sample sizes from regular
crosstabulation on dataset "data" but the chi-square using survey
methods on "audit").

Here's my code (I can't supply data for you as I am not that
sophisticated and the real data is not cleared for public consumption -
I really apologize):



You could have checked to see if the same problem occurs using one of the 
built-in data sets, but it turns out that the diagnosis is fairly 
straightforward. Also, since you use the function 'crosstab' without saying 
which package you got it from, it might be hard to get your code to run anyway.

Your problem is that you are passing a bare name and trying to substitute it 
into a formula, and that isn't how R function arguments or formulas work.

Simplifying the problem to the essentials, your function would work only if

myFormulaMakingFunction<-function(svyX){
   ~svyX+SEX
}

myFormulaMakingFunction(con)

returned ~con+SEX.  It actually returns ~svyX+SEX

There are various ways around this.  My preferred one would be

myPreferredWay<-function(formula){
 update(formula, ~.+SEX)
}

myPreferredWay(~con)

which does return ~con+SEX.

Then, as an example:

data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)

myIllustrativeExample<-function(X,formula){
   tab<-table(X,apiclus1$stype)
   chisq<-svychisq(update(formula, ~.+stype), statistic="adjWald",design=dclus1)
   list(tab,chisq)
}

myIllustrativeExample(apiclus1$comp.imp, ~comp.imp)


Now, this can be improved: it appears that the vector argument is supposed to 
be the same variable as the formula argument, so we shouldn't make the user 
supply both of them.

mySupererogatoryEffort <- function(formula, design){
X<-model.frame(formula, model.frame(design))[[1]]
   tab<-table(X, model.frame(design)$stype)
   chisq<-svychisq(update(formula, ~.+stype), statistic="adjWald",design=design)
   list(table=tab, F=chisq$statistic, df=chisq$parameter, p=chisq$p.value)
}

mySupererogatoryEffort(~comp.imp, dclus1)


As a final note, there is no round= argument to svychisq().  There is a round= 
argument to svytable(), which is documented on the same help page, but 4 is not 
one of the allowed values.

 -thomas


# create my svydesign object

audit <- svydesign(id~id, strata=~field, weights=~wt, data=data,
fpc=~AllocProportion)

# my function to create my table

mkMyCrossTable <- function(X, svyX, T) {

tbl <- crosstab(X, data$SEX, prop.c=TRUE)
tbl <- data.frame(cbind(tbl$t, tbl$prop.col))
tbl$var <- rownames(tbl)

chisq <- svychisq(~svyX + SEX, design=audit, statistic="adjWald",
round=4)
chisq <- data.frame(do.call("cbind", chisq)
chisq <- data.frame(chisq[,3])

Table <- data.frame(tbl$var,
   paste(formatC(tbl$X0.1*100, format="f", digits=1),
"%", sep=""),
  tbl$X0,
  paste(formatC(tbl$X1.1*100, format="f",
digits=1), "%", sep=""),
  tbl$X1,
  chisq[1])
Table[2: length(Table[,1]), 6] <- NA
Table <- NAToUnkown(Table, unknown = " ")
Colnames(Table) <- c(T, "Male (%)", "Male (n)", "Female (%)", "Female
(n)", "p-value")
Table

}

con3 <- mkMyCrossTable(data$con, con, "Constituency")




Thomas Lumley   Assoc. Professor, Biostatistics
tlum...@u.washington.eduUniversity of Washington, Seattle

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Svy function doesn't work nested in user-defined function

2010-05-20 Thread Sabatier, Jennifer F. (CDC/OID/NCHHSTP)
 
Hi R-help,

I posted about this late yesterday but got no response. I may have put
TMI in the original request.  Not to mention I couldn't cut and paste
yesterday because I was working R off a non-network computer while
asking for help on a network computer.


Essentially, I have this user-defined function:

test <- function(X){
  

chisq <- svychisq(~X + SEX, design=audit,
statisic="adjWald", round=4)

}

test(con)



"con" is a data variable in my design object audit.

When I just run: 

chisq <- svychisq(~X + SEX, design=audit, statisic="adjWald", round=4)

It works just fine.  

It's only when it's nested in the function test() that it falls apart.  

I get this error:

Error in `[.data.frame`(design$variables, , as.character(rows)) : 
  undefined columns selected


How can I solve this problem?

Thanks,

Jen (obviously a new R-user, as of late 2009)

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


Re: [R] esthetics --- extending the lm command to fixed effects?

2010-05-20 Thread Thomas Lumley

On Thu, 20 May 2010, ivo welch wrote:


dear R wizards:

not important.  more a curiosity or esthetics question.

is there a way to extend the standard lm command, so that it takes a new
argument that handles fixed effects?   right now, I have (provided to me
from an expert---I would have never figured this one out):

  diffid <- function(h,id) {
  id <- as.factor(id)[, drop=TRUE]
  apply(as.matrix(h), 2, function(x) x - tapply(x,id,mean)[id]
  }


Simpler would be

   diffid<-function(h,id){ h-ave(h,id)}


which is used as

r= lm( diffid(y, firmid) ~ diffid(x, firmid ) )

it works, but it would be much nicer if I could just write

   r= lm( y ~ x + z, fixed.effects=firmid )

does this already exists as a package?  or has someone figured out how to
program this?


I would just have used lm(y~x+z+factor(firmid)).  Admittedly, you get a whole 
bunch of uninteresting coefficients in the output, but it's not that hard to 
subset them out.

There are two implementation of this in Bill Venables' course notes on advanced 
programming. I think they are also in 'S Programming', but I can't find my copy 
right now.  These were motivated by computational problems: the full design 
matrix for the linear model was too large for memory at the time (last century).


As a final note, I would strongly discourage
   r= lm( y ~ x + z, fixed.effects=firmid )
as a specification, and would argue for
   r= lm( y ~ x + z, fixed.effects=~firmid )

I think the ability to have some subset of the arguments in a modelling call 
silently treated as formulas was a bad decision, although it must have looked 
user-friendly at the time.

 -thomas

Thomas Lumley   Assoc. Professor, Biostatistics
tlum...@u.washington.eduUniversity of Washington, Seattle

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

2010-05-20 Thread kaveh vakili
Dear List,

I'd like to call pyhton function from within R. I tried installing the latest 
version of RSPython:

wget http://www.omegahat.org/RSPython/RSPython_0.7-1.tar.gz
R CMD INSTALL --clean RSPython_0.7-1.tar.gz

I get a compile error (posted below). 

Did anyone else run against this ? Is there a solution ?

checking for python... /usr/bin/python
Python version 2.6
Using threads
checking for gcc... gcc
checking for C compiler default output file name... a.out
checking whether the C compiler works... yes
checking whether we are cross compiling... no
checking for suffix of executables... 
checking for suffix of object files... o
checking whether we are using the GNU C compiler... yes
checking whether gcc accepts -g... yes
checking for gcc option to accept ISO C89... none needed
R version 2
Looking for libR.so in lib/
configure: creating ./config.status
config.status: creating Makefile
config.status: creating src/Makevars
config.status: creating cleanup
config.status: creating inst/scripts/RPython.csh
config.status: creating inst/scripts/RPython.bsh
** libs
gcc -std=gnu99 -I/usr/share/R/include -I../inst/include 
-I/usr/include/python2.6 
-D_R_=1 -DUSE_R=1 -fpic  -g -O2 -c GeneralConverters.c -o 
GeneralConverters.o
In file included from ../inst/include/UserConverters.h:4,
 from GeneralConverters.c:1:
../inst/include/RPythonModule.h:4:20: error: Python.h: No such file or directory
In file included from ../inst/include/UserConverters.h:4,
 from GeneralConverters.c:1:
../inst/include/RPythonModule.h:14: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or 
‘__attribute__’ before ‘*’ token
../inst/include/RPythonModule.h:15: error: expected ‘)’ before ‘*’ token
../inst/include/RPythonModule.h:17: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or 
‘__attribute__’ before ‘*’ token
In file included from GeneralConverters.c:1:
../inst/include/UserConverters.h:13: error: expected ‘)’ before ‘*’ token
../inst/include/UserConverters.h:15: error: expected ‘)’ before ‘*’ token
../inst/include/UserConverters.h:17: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or 
‘__attribute__’ before ‘FromTargetConverterMatch’
../inst/include/UserConverters.h:18: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or 
‘__attribute__’ before ‘FromTargetConverter’
../inst/include/UserConverters.h:24: error: expected specifier-qualifier-list 
before ‘FromPythonConverterMatch’
../inst/include/UserConverters.h:35: error: expected declaration specifiers or 
‘...’ before ‘PyClassObject’
In file included from GeneralConverters.c:1:
../inst/include/UserConverters.h:37: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or 
‘__attribute__’ before ‘*’ token
../inst/include/UserConverters.h:40: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or 
‘__attribute__’ before ‘ToTargetConverter’
../inst/include/UserConverters.h:46: error: expected specifier-qualifier-list 
before ‘ToPythonConverter’
../inst/include/UserConverters.h:63: error: expected ‘)’ before ‘*’ token
../inst/include/UserConverters.h:64: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or 
‘__attribute__’ before ‘*’ token
../inst/include/UserConverters.h:74: error: expected ‘)’ before ‘Rf_match’
../inst/include/UserConverters.h:78: error: expected declaration specifiers or 
‘...’ before ‘ToPythonConverter’
GeneralConverters.c: In function ‘removeFromTargetConverterByIndex’:
GeneralConverters.c:23: error: ‘RSFromPythonConverter’ has no member named 
‘next’
GeneralConverters.c:29: error: ‘RSFromTargetConverter’ has no member named 
‘next’
GeneralConverters.c:36: error: ‘RSFromTargetConverter’ has no member named 
‘next’
GeneralConverters.c:36: error: ‘RSFromTargetConverter’ has no member named 
‘next’
GeneralConverters.c: In function ‘removeToTargetConverterByIndex’:
GeneralConverters.c:53: error: ‘RSToPythonConverter’ has no member named ‘next’
GeneralConverters.c:59: error: ‘RSToTargetConverter’ has no member named ‘next’
GeneralConverters.c:66: error: ‘RSToTargetConverter’ has no member named ‘next’
GeneralConverters.c:66: error: ‘RSToTargetConverter’ has no member named ‘next’
GeneralConverters.c: In function ‘removeFromTargetConverterByDescription’:
GeneralConverters.c:85: error: ‘RSFromTargetConverter’ has no member named 
‘next’
GeneralConverters.c:85: error: ‘RSFromTargetConverter’ has no member named 
‘next’
GeneralConverters.c:87: error: ‘RSFromTargetConverter’ has no member named 
‘next’
GeneralConverters.c:94: error: ‘RSFromTargetConverter’ has no member named 
‘next’
GeneralConverters.c: In function ‘RPython_removeConverter’:
GeneralConverters.c:128: warning: implicit declaration of function ‘free’
GeneralConverters.c:128: warning: incompatible implicit declaration of built-in 
function ‘free’
GeneralConverters.c:135: error: ‘RSToTargetConverter’ has no member named 
‘description’
GeneralConverters.c:138: warning: incompatible implicit declaration of built-in 
function ‘free’
GeneralConverters.c: In function ‘getNumConverters’:
GeneralConverters.c:160: error: ‘RSFromTargetConverter’ has no member named 
‘next’
Gener

[R] Comparing three groups, data: present, absent

2010-05-20 Thread Mächler Marc Jacques
Dear R-Experts, Dear friends of dung.

I have a statistical Problem, to which nobody I asked could give me an answer. 
Maybe you can.

I was in the African-Savanna and made a Dung-Monitoring. This means I walked 
randomly over the field and for every Dung-Event I found I noted following 
parameters:

Species (Waterbuck, Giraffe, Reedbuck); 
Age-Category (less than a week, more than a week & less than a month, more than 
a month);
Termites present (0 or 1);

Now I want to compare for every species, if the percentage of Termites attack 
differ between the different Age-Categories

and if the percentage of attacks of Species differ between the same  
Age-Categories.

I only found a lot of very complicated models for "Presence / Absence". But I 
do not need to find a treshold, I just want to compare groups. 

My Question now:

- What test do I have to use to compare groups.
- Do I have to transform the data (Termite.Attack which is 0 and 1)?

I really appreciate every hint.

Thanks a lot, 

Dung-Master Marc-Jacques from the swiss federal institute of technology 

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


Re: [R] esthetics --- extending the lm command to fixed effects?

2010-05-20 Thread ivo welch
hi thomas---

thanks for the answer.  the problem with the "+factor(fmid)" is not just
that it provides uninteresting coefficients and that it eats more memory,
but that it is also MUCH slower when there are (hundred of) thousands of
fixed effects.

Does Bill Venables describe how to do extend the lm() function?  I googled
"course notes on advanced programming Venables", but did not find it.  Do
you have a better link?   (hopefully, this is a short explanation---I know
the algorithm.  I want to learn how to coax it into an lm statement.)

regards,

/iaw

Ivo Welch (ivo.we...@brown.edu, ivo.we...@gmail.com)


I would just have used lm(y~x+z+factor(firmid)).  Admittedly, you get a
> whole bunch of uninteresting coefficients in the output, but it's not that
> hard to subset them out.
>
> There are two implementation of this in Bill Venables' course notes on
> advanced programming. I think they are also in 'S Programming', but I can't
> find my copy right now.  These were motivated by computational problems: the
> full design matrix for the linear model was too large for memory at the time
> (last century).
>
>
> As a final note, I would strongly discourage
>
>   r= lm( y ~ x + z, fixed.effects=firmid )
> as a specification, and would argue for
>
>   r= lm( y ~ x + z, fixed.effects=~firmid )
>
> I think the ability to have some subset of the arguments in a modelling
> call silently treated as formulas was a bad decision, although it must have
> looked user-friendly at the time.
>
> -thomas
>
> Thomas Lumley   Assoc. Professor, Biostatistics
> tlum...@u.washington.eduUniversity of Washington, Seattle
>
>

[[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] RSpython Ubuntu

2010-05-20 Thread kaveh vakili

kaveh vakili  ulb.ac.be> writes:

> 
> Dear List,
> 
> I'd like to call pyhton function from within R. I tried installing the latest 
> version of RSPython:
> 
> wget http://www.omegahat.org/RSPython/RSPython_0.7-1.tar.gz
> R CMD INSTALL --clean RSPython_0.7-1.tar.gz
> 
> I get a compile error (posted below). 
> 
> Did anyone else run against this ? Is there a solution ?


Sorry, solved,
i needed pyhton headers.

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

2010-05-20 Thread Luc Villandre

Dear R-users,

I've recently switched to Ubuntu and I've decided to use Kate to edit my 
R code. I really like how Kate allows one to simply pipe their code to 
the terminal. However, I would find it even better if I could actually 
get the Kate console to display error messages (or any console outputs 
in fact) in red and my input in a different color. The instance of xterm 
that Ubuntu seems to run by default supports color coding. Can the Kate 
console be set up to do just the same?


Although this is clearly an issue related to Kate, I was thinking that 
some people on this list might also use Kate in conjunction with R and 
might have already inquired about this.


I thank you very much for your help.

Luc Villandré

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


Re: [R] offset in gam and spatial scale of variables

2010-05-20 Thread Joris Meys
On Thu, May 20, 2010 at 3:20 PM, Lucia Rueda  wrote:

>
> Hi,
>
> Thanks for the inputs. I talked to my coworker, who has been the one doing
> the analysis. Perhaps I wasn't making myself clear about the “differences
> in
> spatial scales”.  Here is what he says:
>
> "The truth is that measuring scales (i.e all area related variable are
> measured in m2) and spatial definition of initial cartography are
> homogeneous among extracted variables. But all variables (ie. sum of the
> total rocky bottom in the surrounding area) are computed for each different
> integration areas (buffer) (i.e in an area of 40squaremeters around the
> sample, in an area of 80m2, …).
> The question is then if we can build a model that includes variables
> measured at different buffers (for example a model that includes 3
> variables:  1.-  the amount of rocky bottom in an area of 80m2 ; 2- the
> amount of sandy bottom in an area of 200m2; and the mean depth calculated
> in
> a surrounding area of 50m2) considering that each variable may be
> expressing
> different ecological processes. I believe that if there is not an
> ecological
> constrain in the interpretation of the variables (and their ecological
> effect over the specie), including them in a model is correct, unless there
> is not a mathematical constrain."
>

If you look upon it that way, you might indeed consider using them in
different buffers, but as you said, you should be able to interprete them in
an ecological way. I'd be surprised if depth and bottom have a different
effect-scale, as they both are related to the territorium of the animal.
Plus, you cannot conclude anything from the difference in deviance
explained. You can't say anything about the homerange or so based on the
observation that more deviance is explained when looking on a scale of 200m2
for example. So if you have good ecological reasons to include them, you
can, but if it's merely because on one scale they explain more of the
deviance, I still believe it is a very dangerous approach...


> Also, about the spatial correlation I thought from what I've read so far
> that I had to build the model and then check if there was spatial
> correlation in the residuals since they are supposed to be i.i.d. And if it
> turns out that they are then I have to do something about it like gamm,
> gee,
> sar, car, etc.
>
That's an approach that is often used. In essence, that's true. Correlation
between the raw data can be due to cocorrelation with some other factor in
space or time. But a pre-analysis of correlations and autocorrelations can
tell you already quite some. In any case, you always have to check the
residuals after the model building. My main point was that using the
correlation will definitely influence the significance of the parameters.

Anyway, good luck with it. I learnt pretty fast that as long as you can
explain what you're doing and why you're doing it, there's a big grey zone
between right and wrong. Otherwise it wouldn't be statistics, would it? ;-)

Cheers
Joris

>
> Cheers,
>
> Lucia
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/offset-in-gam-and-spatial-scale-of-variables-tp483p2224528.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.
>



-- 
Joris Meys
Statistical Consultant

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

Coupure Links 653
B-9000 Gent

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

[[alternative HTML version deleted]]

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


[R] writing function

2010-05-20 Thread arnaud Gaboury
Dear group,

I am trying to write functions, but as a beginner, everything is not so
obvious.
Let's say I want the results in a list of elemts like this :
tot1, tot2, etc

Here is a function:

toto <-
function(x,y)

{

for(i in x:y){

paste(c("tot",i),collapse="")<-(i*2)

}
}

If I type this : 
>toto(1,5)
I get this message error:
Error in paste(c("tot", i), collapse = "") <- (i * 2) : 
  target of assignment expands to non-language object

How can I write it to get the result I want (i.e tot1, tot2... with tot1=2,
tot2=4...) in my environment?

TY for any help

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


Re: [R] writing function

2010-05-20 Thread Tal Galili
Try this:

paste("tot", 4:16, sep = "")


Or:

func <- function(x,y)
{
 paste("tot", x:y, sep = "")
}
func(4,16)




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




On Thu, May 20, 2010 at 6:54 PM, arnaud Gaboury wrote:

> Dear group,
>
> I am trying to write functions, but as a beginner, everything is not so
> obvious.
> Let's say I want the results in a list of elemts like this :
> tot1, tot2, etc
>
> Here is a function:
>
> toto <-
> function(x,y)
>
> {
>
> for(i in x:y){
>
> paste(c("tot",i),collapse="")<-(i*2)
>
> }
> }
>
> If I type this :
> >toto(1,5)
> I get this message error:
> Error in paste(c("tot", i), collapse = "") <- (i * 2) :
>  target of assignment expands to non-language object
>
> How can I write it to get the result I want (i.e tot1, tot2... with tot1=2,
> tot2=4...) in my environment?
>
> TY for any help
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[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] esthetics --- extending the lm command to fixed effects?

2010-05-20 Thread Achim Zeileis

On Thu, 20 May 2010, ivo welch wrote:


hi thomas---

thanks for the answer.  the problem with the "+factor(fmid)" is not just
that it provides uninteresting coefficients and that it eats more memory,
but that it is also MUCH slower when there are (hundred of) thousands of
fixed effects.


There is also the plm() function in the "plm" package which provides fixed 
effects models (among many other models for panel data).


However, if I recall correctly, it internally employs the full regressor 
matrix, not the demeaned one. More details are explained in the vignette, 
though:

  vignette("plm", package = "plm")

hth,
Z


Does Bill Venables describe how to do extend the lm() function?  I googled
"course notes on advanced programming Venables", but did not find it.  Do
you have a better link?   (hopefully, this is a short explanation---I know
the algorithm.  I want to learn how to coax it into an lm statement.)

regards,

/iaw

Ivo Welch (ivo.we...@brown.edu, ivo.we...@gmail.com)


I would just have used lm(y~x+z+factor(firmid)).  Admittedly, you get a

whole bunch of uninteresting coefficients in the output, but it's not that
hard to subset them out.

There are two implementation of this in Bill Venables' course notes on
advanced programming. I think they are also in 'S Programming', but I can't
find my copy right now.  These were motivated by computational problems: the
full design matrix for the linear model was too large for memory at the time
(last century).


As a final note, I would strongly discourage

  r= lm( y ~ x + z, fixed.effects=firmid )
as a specification, and would argue for

  r= lm( y ~ x + z, fixed.effects=~firmid )

I think the ability to have some subset of the arguments in a modelling
call silently treated as formulas was a bad decision, although it must have
looked user-friendly at the time.

-thomas

Thomas Lumley   Assoc. Professor, Biostatistics
tlum...@u.washington.eduUniversity of Washington, Seattle




[[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] sqldf: issues with natural joins

2010-05-20 Thread Nick Switanek
Hello,

I'm having trouble discovering what's going wrong with my use of natural
joins via sqldf.

Following the instructions under 4i at http://code.google.com/p/sqldf/,
which discusses creating indices to speed joins, I have been only unreliably
able to get natural joins to work.

For example,

> Tid <- c('AES 01-01-02 10:58:00', 'AES 01-01-02 11:53:00', 'AES 01-01-05
10:58:00', 'AES 01-01-11 12:30:00')
> A <- data.frame(Tid, dfName = 'a')
> B <- data.frame(Tid = Tid[2:4], dfName = 'b')
> C <- data.frame(Tid = Tid[1:3], dfName = 'c')

# then use the sqldf library
> library(sqldf)
> sqldf()

# to create indices on the Tid variable shared across data.frames
> sqldf('create index indA on A(Tid)')
> sqldf('create index indB on B(Tid)')
> sqldf('create index indC on C(Tid)')

# check to make sure everything is there
> sqldf('select * from sqlite_master')

# doing a natural join (implicitly on Tid)
# does not give the expected joins
> sqldf('select * from main.A natural join main.B')
[1] TiddfName
<0 rows> (or 0-length row.names)
> sqldf('select * from main.A natural join main.C')
[1] TiddfName
<0 rows> (or 0-length row.names)
> sqldf('select * from main.B natural join main.C')
[1] TiddfName
<0 rows> (or 0-length row.names)

# even using a where clause (which doesn't have the efficiency qualities I
need the indexed natural joins for) is problematic, setting values of the
dfName variable incorrectly for the data from C
> sqldf('select * from main.B b, main.C c where b.Tid = c.Tid')
Tid dfName   Tid dfName
1 AES 01-01-02 11:53:00  b AES 01-01-02 11:53:00  b
2 AES 01-01-05 10:58:00  b AES 01-01-05 10:58:00  b

I'm grateful for your guidance on what I'm doing wrong with the natural join
in sqldf.

many thanks,
Nick

[[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] [Off topic?] Time dependent Cox model fitting and validation

2010-05-20 Thread Marco Barbàra

this is self-reply only to insert my real name

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


[R] Question about difftime()

2010-05-20 Thread Stella Pachidi
Dear R experts,

I have a question about the result of difftime() function: Does it
take into account the different number of days in each month. In my
example, I have the following:

> firstDay
[1] "2010-02-20"
> lastDay
[1] "2010-05-20 16:00:00"
> difftime(lastDay,firstDay,units='days')
Time difference of 89.625 days
>

When I count the days I get 88 days from 20/02/2010 to 20/05/2010
consequently the difference in days should be 87.

On the contrary, difftime gives a higher number, so I doubt whether it
takes into account the fact that february has 28 days (or 29). Could
you please help?

Thank you very much in advance.

Kind regards,
Stella

--
Stella Pachidi
Master in Business Informatics student
Utrecht University

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


Re: [R] p-values < 2.2e-16 not reported

2010-05-20 Thread Will Eagle

On 2010-05-20 08:52, Shi, Tao wrote:

Will,

  I'm wondering if you have any
insights after looking at the cor.test source code.  It seems to be fine to me, as the p 
value is either calculated by "your first method" or a
.C code.

...Tao


Dear Tao,

I think the described problem of p-values < 2.2e-16 not being reported 
(or set to zero, respectively) applies at least to the cor.test(, 
method="pearson", alternative=c("greater","two.sided")) function, but 
also other statistical test functions. You find the critical piece of 
code with methods(cor.test) and getAnywhere(cor.test.default) in the 
following lines:

...
p <- pt(STATISTIC, df) # line 27
...
PVAL <- switch(alternative,  # lines 144-145, reformatted
less = p,
greater = 1 - p,
two.sided = 2 * min(p, 1 - p))

This means that if you calculate a pearson correlation (which is the 
default) with option alternativ="two.sided" (which is the default) or 
alternative="greater", an additive operation of the 2nd kind is 
performed, which causes that a very small p-value is set to zero. This 
problem should not apply to alternative="less" since the p-value has 
been generated without an additive operation with the pt() function.


What I still find confusing is that:

> pt(10,100,lower=FALSE)
[1] 4.950844e-17
> pt(10,100,lower=TRUE)   # should give 1 - 4.950844e-17
[1] 1
> 1-pt(10,100,lower=TRUE) # should give 4.950844e-17
[1] 0

So it seems that this distribution function(s) only give the exact value 
if the option lower.tail=FALSE is used. Therefore, also cor.test(..., 
method="pearson", alternative="less") might not report p-values < 
2.2e-16. I find this behavior of R very unintuitive.


The .C code routines are by the way only used for methods=c("kendall", 
"spearman") as far as I have seen.


Best, Will

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


Re: [R] finding euclidean proximate points in two datasets

2010-05-20 Thread David Winsemius


On May 20, 2010, at 11:12 AM, Alexander Shenkin wrote:


On 5/20/2010 9:18 AM, David Winsemius wrote:


On May 20, 2010, at 10:02 AM, Alexander Shenkin wrote:


Hello all,

I've been pouring through the various spatial packages, but  
haven't come

across the right thing yet.


There is a SIG for such questions.


thanks - just joined it.

Given a set of points in 2-d space X, i'm trying to find the  
subset of
points in Y proximate to each point in X.  Furthermore, the  
proximity
threshold of each point in X differs (X$threshold).  I've  
constructed

this myself already, but it's horrificly slow with a dataset of 40k+
points in one set, and a 700 in the other.

A very inefficient example of what I'm looking for:


Not really a reproducible example. If euclidean_dist is a function ,
then it is not one in any of the packages I have installed.


it's not reproducible - i'll make a better effort to include
reproducible code in the future.  and that function is just one i  
would
have written, but didn't want to clog the email with code.  Anyway,  
here

is a reproducible example:

X = data.frame(x=c(1,2,3), y=c(2,3,1), threshold=c(1,2,4))
Y = data.frame(x=c(5,2,3,4,2,5,2,3), y=c(5,2,2,4,1,2,3,1))
proximate=list()
i=1
for (pt in 1:length(X$x)) {
   proximate[[i]] <- sqrt((X[pt,]$x - Y$x)^2 + (X[pt,]$y - Y$y)^2) >
X[pt,]$threshold
   i=i+1
}
proximate


This might be a more R-ish way of approaching the problem. It is not  
always the case that apply functions will be faster than looping, but  
they may be more compact and they get you ready for thinking about  
indexed solutions which would be generally faster.


> apply(X, 1, function(.x) {  (.x[1] - Y$x )^2+(.x[2] - Y$y)^2  
<.x[3]^2 } )


  [,1]  [,2]  [,3]
[1,] FALSE FALSE FALSE
[2,] FALSE  TRUE  TRUE
[3,] FALSE  TRUE  TRUE
[4,] FALSE FALSE  TRUE
[5,] FALSE FALSE  TRUE
[6,] FALSE FALSE  TRUE
[7,] FALSE  TRUE  TRUE
[8,] FALSE FALSE  TRUE

Transposing lets you check that the results are the same as your  
solution modulo yours and mine being slightly different classes:


> t( apply(X, 1, function(.x) {  (.x[1] - Y$x )^2+(.x[2] - Y$y)^2  
<.x[3]^2 } ) )

  [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
[1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[2,] FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE
[3,] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

> proximate
[[1]]
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

[[2]]
[1] FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE

[[3]]
[1] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE





  for (pt in X$idx) {
  proximity[i] = euclidian_dist(X[pt]$x, X[pt]$y, Y$x, Y$y) <
X$threshold
   i = i+1
  }



Have you considered first creating a subset of candidate points  
that are

within "threshold" of each reference point on both coordinates. That
might sidestep a lot of calculations on points that are easily
eliminated on a single comparison. Then you could calculate distances
within that surviving subset of points. On average that should give  
you

an over 50% "hit rate":


(4/3)*pi*0.5^3

[1] 0.5235988


That's a nice idea.  I'll still be waiting quite a while while my
machine cranks, but not as long.  Still - I suspect there would be  
much

bigger gains if there were tailored packages.  I'll re-post over on
sig-geo about that.  thanks.

Perhaps crossdist() in spatstat is what I should use, and then  
code a

comparison with X$threshold after the cross-distances are computed.
However, I was wondering if there was another tool I should be
considering.  Any and all thoughts are very welcome.  Thanks in  
advance.


Thanks,
Allie
--
Alexander Shenkin
PhD Candidate
School of Natural Resources and Environment
University of Florida


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] writing function

2010-05-20 Thread Ivan Calandra

Hi,

the problem is not about the environment, but about the
paste(c("tot", i), collapse = "")
which is not recognized as an object.
Maybe assign() could do the trick

You could also do it this way (though it's not exactly what you want, 
but it might be better):

toto <- function(x,y){
  tot <- vector(mode="list", length=(y-x+1))   ##create an empty list 
of correct length

  for (i in x:y) {
tot[[i]] <- (i*2)   ##store each value into each element of the list
  }
  return(tot)  ##return the list
}

HTH,
Ivan




Le 5/20/2010 17:54, arnaud Gaboury a écrit :

Dear group,

I am trying to write functions, but as a beginner, everything is not so
obvious.
Let's say I want the results in a list of elemts like this :
tot1, tot2, etc

Here is a function:

toto<-
function(x,y)

{

for(i in x:y){

paste(c("tot",i),collapse="")<-(i*2)

}
}

If I type this :
   

toto(1,5)
 

I get this message error:
Error in paste(c("tot", i), collapse = "")<- (i * 2) :
   target of assignment expands to non-language object

How can I write it to get the result I want (i.e tot1, tot2... with tot1=2,
tot2=4...) in my environment?

TY for any help

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

   


--
Ivan CALANDRA
PhD Student
University of Hamburg
Biozentrum Grindel und Zoologisches Museum
Abt. Säugetiere
Martin-Luther-King-Platz 3
D-20146 Hamburg, GERMANY
+49(0)40 42838 6231
ivan.calan...@uni-hamburg.de

**
http://www.for771.uni-bonn.de
http://webapp5.rrz.uni-hamburg.de/mammals/eng/mitarbeiter.php

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


Re: [R] About the breakpoint when making heatmap with lots of variables

2010-05-20 Thread Shi, Tao
The "breakpoint" you mentioned is irrelevant here.  Clustering 15672 
genes (I assume this is a microarray data) requires lots of memory.  I 
suggest you either filter your gene list down to thousands or just plot 
the column dendrogram without showing the heatmap.

plot(hclust(dist(x)))

...Tao




- Original Message 
> From: Meng Wu 
> To: r-help@r-project.org
> Sent: Wed, May 19, 2010 7:31:29 PM
> Subject: [R] About the breakpoint when making heatmap with lots of variables
> 
> HI,All

I am trying to create a heatmap with 24 samples with 15672 
> varibles, I read
in the table in R, and then made it as a matrix, then try to 
> create the
heatmap using heatmap(x,...)

However, I received the error 
> message as:
> heatmap(t(x))
Error: cannot allocate vector of size 936.8 
> Mb
R(2925,0xa0b16500) malloc: *** mmap(size=982261760) failed (error 
> code=12)
*** error: can't allocate region
*** set a breakpoint in 
> malloc_error_break to debug
R(2925,0xa0b16500) malloc: *** 
> mmap(size=982261760) failed (error code=12)
*** error: can't allocate 
> region
*** set a breakpoint in malloc_error_break to 
> debug
R(2925,0xa0b16500) malloc: *** mmap(size=982261760) failed (error 
> code=12)
*** error: can't allocate region
*** set a breakpoint in 
> malloc_error_break to debug
R(2925,0xa0b16500) malloc: *** 
> mmap(size=982261760) failed (error code=12)
*** error: can't allocate 
> region
*** set a breakpoint in malloc_error_break to 
> debug
R(2925,0xa0b16500) malloc: *** mmap(size=982261760) failed (error 
> code=12)
*** error: can't allocate region
*** set a breakpoint in 
> malloc_error_break to debug
R(2925,0xa0b16500) malloc: *** 
> mmap(size=982261760) failed (error code=12)
*** error: can't allocate 
> region
*** set a breakpoint in malloc_error_break to 
> debug
R(2925,0xa0b16500) malloc: *** mmap(size=982261760) failed (error 
> code=12)
*** error: can't allocate region
*** set a breakpoint in 
> malloc_error_break to debug


It seems like I do not have enough 
> memory, how can I set a breakpoint 
> as
suggested?

Thanks,

Meng


> [[alternative HTML version 
> deleted]]

__

> ymailto="mailto:R-help@r-project.org"; 
> href="mailto:R-help@r-project.org";>R-help@r-project.org mailing list

> href="https://stat.ethz.ch/mailman/listinfo/r-help"; target=_blank 
> >https://stat.ethz.ch/mailman/listinfo/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] sqldf: issues with natural joins

2010-05-20 Thread Gabor Grothendieck
There are two problems:

1. A natural join will join all columns with the same names in the two
tables and that includes not only Tid but also dfName and since there
are no rows that have the same Tid and dfName the result has zero
rows.

2. the heuristic it uses fails when you retrieve the same column name
from multiple tables so use method = "raw" to turn off the heuristic.
The heuristic will be improved to cover this case in the future.
Read FAQ #1 on the home page:
http://code.google.com/p/sqldf/#1._How_does_sqldf_handle_classes_and_factors?

This should work:

> sqldf('select * from main.A join main.B using(Tid)', method = "raw")
 Tid dfName dfName
1  AES 01-01-02 11:53:00  a  b
2 AES 01-01-05\n10:58:00  a  b
3  AES 01-01-11 12:30:00  a  b

This works too as the double dfName no longer exists to confuse the heuristic:

names(B)[2] <- "dfNameB"
sqldf('select * from main.A join main.B using(Tid)')



On Thu, May 20, 2010 at 12:04 PM, Nick Switanek  wrote:
> Hello,
>
> I'm having trouble discovering what's going wrong with my use of natural
> joins via sqldf.
>
> Following the instructions under 4i at http://code.google.com/p/sqldf/,
> which discusses creating indices to speed joins, I have been only unreliably
> able to get natural joins to work.
>
> For example,
>
>> Tid <- c('AES 01-01-02 10:58:00', 'AES 01-01-02 11:53:00', 'AES 01-01-05
> 10:58:00', 'AES 01-01-11 12:30:00')
>> A <- data.frame(Tid, dfName = 'a')
>> B <- data.frame(Tid = Tid[2:4], dfName = 'b')
>> C <- data.frame(Tid = Tid[1:3], dfName = 'c')
>
> # then use the sqldf library
>> library(sqldf)
>> sqldf()
>
> # to create indices on the Tid variable shared across data.frames
>> sqldf('create index indA on A(Tid)')
>> sqldf('create index indB on B(Tid)')
>> sqldf('create index indC on C(Tid)')
>
> # check to make sure everything is there
>> sqldf('select * from sqlite_master')
>
> # doing a natural join (implicitly on Tid)
> # does not give the expected joins
>> sqldf('select * from main.A natural join main.B')
> [1] Tid    dfName
> <0 rows> (or 0-length row.names)
>> sqldf('select * from main.A natural join main.C')
> [1] Tid    dfName
> <0 rows> (or 0-length row.names)
>> sqldf('select * from main.B natural join main.C')
> [1] Tid    dfName
> <0 rows> (or 0-length row.names)
>
> # even using a where clause (which doesn't have the efficiency qualities I
> need the indexed natural joins for) is problematic, setting values of the
> dfName variable incorrectly for the data from C
>> sqldf('select * from main.B b, main.C c where b.Tid = c.Tid')
>                    Tid dfName                   Tid dfName
> 1 AES 01-01-02 11:53:00      b AES 01-01-02 11:53:00      b
> 2 AES 01-01-05 10:58:00      b AES 01-01-05 10:58:00      b
>
> I'm grateful for your guidance on what I'm doing wrong with the natural join
> in sqldf.
>
> many thanks,
> Nick

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


Re: [R] Question about difftime()

2010-05-20 Thread Marc Schwartz
On May 20, 2010, at 11:08 AM, Stella Pachidi wrote:

> Dear R experts,
> 
> I have a question about the result of difftime() function: Does it
> take into account the different number of days in each month. In my
> example, I have the following:
> 
>> firstDay
> [1] "2010-02-20"
>> lastDay
> [1] "2010-05-20 16:00:00"
>> difftime(lastDay,firstDay,units='days')
> Time difference of 89.625 days
>> 
> 
> When I count the days I get 88 days from 20/02/2010 to 20/05/2010
> consequently the difference in days should be 87.
> 
> On the contrary, difftime gives a higher number, so I doubt whether it
> takes into account the fact that february has 28 days (or 29). Could
> you please help?
> 
> Thank you very much in advance.
> 
> Kind regards,
> Stella


Try recounting:

8 days remaining in February this year (not a leap year), after the 20th.
31 days in March
30 days in April
20 days in May

Total?  

HTH,

Marc Schwartz

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


Re: [R] Kate terminal window

2010-05-20 Thread Vojtěch Zeisek
Dne Čt 20. května 2010 17:50:42 Luc Villandre napsal(a):
> Dear R-users,
> 
> I've recently switched to Ubuntu and I've decided to use Kate to edit my
> R code. I really like how Kate allows one to simply pipe their code to
> the terminal. However, I would find it even better if I could actually
> get the Kate console to display error messages (or any console outputs
> in fact) in red and my input in a different color. The instance of xterm
> that Ubuntu seems to run by default supports color coding. Can the Kate
> console be set up to do just the same?

Hello,
I do not know how it is in Kate, but You can use Rkward (the package is in 
Ubuntu's repositories). It is GUI for R. It uses Kate part as integrated text 
editor (so it behaves as Kate and it shares its configuration) and it has log, 
error console and all needed stuff...
Well, might be, it is not what You wished to hear, but it works and I hope it 
helps little bit. :-)
Best regards,
Vojtěch Zeisek

> Although this is clearly an issue related to Kate, I was thinking that
> some people on this list might also use Kate in conjunction with R and
> might have already inquired about this.
> 
> I thank you very much for your help.
> 
> Luc Villandré
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
>  http://www.R-project.org/posting-guide.html and provide commented,
>  minimal, self-contained, reproducible code.
-- 
Vojtěch Zeisek

Komunita openSUSE GNU/Linuxu /
Community of the openSUSE GNU/Linux

http://www.opensuse.org/
http://web.natur.cuni.cz/~zeisek/


signature.asc
Description: This is a digitally signed message part.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Strange behaviour when using diff with POSIXt and POSIXlt objects

2010-05-20 Thread Julian Burgos
Dear list,

I´m calculating time differences between series of time stamps and I noticed
something odd:

If I do this...

> time1=strptime("2009 05 31 22 57 00",format="%Y %m %d %H %M")
> time2=strptime("2009 05 31 23 07 00",format="%Y %m %d %H %M")
>
> diff(c(time1,time2),units="mins")
Time difference of 10 mins

.. I get the correct response in minutes.  But if I try the same thing with
different values, say..

> time3=strptime("2009 06 01 00 47 00",format="%Y %m %d %H %M")
> time4=strptime("2009 06 01 00 57 00",format="%Y %m %d %H %M")
>
> diff(c(time3,time4))
Time difference of NA secs

...which is not what I´m looking for. The difference should also be 10
minutes.

I burned a few neurons (and searched the documentation) and I cannot figure
why this happens.  Any ideas?

All the best,

Julian

Julian Mariano Burgos
Hafrannsóknastofnunin/Marine Research Institute
Skúlagata 4, 121 Reykjavík, Iceland
Sími/Telephone : +354-5752037
Bréfsími/Telefax:  +354-5752001
Netfang/Email: jul...@hafro.is, jmbur...@uw.edu

[[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] Courses***More R courses scheduled for June 2010 by XLSolutions Corp

2010-05-20 Thread Sue Turner
We've scheduled more R courses for June 2010 in Washington DC,
San Francisco, Seattle, New York City and other USA cities

http://www.xlsolutions-corp.com/Rcourses

May - June 2010

*** R/S: Programming Essentials
*** R Fundamentals and Programming Techniques
*** R/S-PLUS Functions by Example
*** R/S Systems: Advanced Programming


More on website

http://www.xlsolutions-corp.com/Rcourses

Ask for group discount and reserve your seat Now - Earlybird Rates.
Payment due after the class! Email Sue Turner: sue at
xlsolutions-corp.com

Phone: 206-686-1578

Please let us know if you and your colleagues are interested in this
class to take advantage of group discount. Register now to secure your
seat.

Cheers,
Elvis Miller, PhD
Manager Training.
XLSolutions Corporation
206 686 1578
www.xlsolutions-corp.com
elvis at xlsolutions-corp.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] sqldf: issues with natural joins

2010-05-20 Thread Gabor Grothendieck
Although that works I had meant to write:

> names(B)[2] <- "dfNameB"
> # ... other commands
> sqldf('select * from main.A natural join main.B')

so that now only Tid is in common so the natural join just picks it up
and also the heuristic works again since we no longer retrieve
duplicate column names.

On Thu, May 20, 2010 at 12:32 PM, Gabor Grothendieck
 wrote:
> There are two problems:
>
> 1. A natural join will join all columns with the same names in the two
> tables and that includes not only Tid but also dfName and since there
> are no rows that have the same Tid and dfName the result has zero
> rows.
>
> 2. the heuristic it uses fails when you retrieve the same column name
> from multiple tables so use method = "raw" to turn off the heuristic.
> The heuristic will be improved to cover this case in the future.
> Read FAQ #1 on the home page:
> http://code.google.com/p/sqldf/#1._How_does_sqldf_handle_classes_and_factors?
>
> This should work:
>
>> sqldf('select * from main.A join main.B using(Tid)', method = "raw")
>                     Tid dfName dfName
> 1  AES 01-01-02 11:53:00      a      b
> 2 AES 01-01-05\n10:58:00      a      b
> 3  AES 01-01-11 12:30:00      a      b
>
> This works too as the double dfName no longer exists to confuse the heuristic:
>
> names(B)[2] <- "dfNameB"
> sqldf('select * from main.A join main.B using(Tid)')
>
>
>
> On Thu, May 20, 2010 at 12:04 PM, Nick Switanek  wrote:
>> Hello,
>>
>> I'm having trouble discovering what's going wrong with my use of natural
>> joins via sqldf.
>>
>> Following the instructions under 4i at http://code.google.com/p/sqldf/,
>> which discusses creating indices to speed joins, I have been only unreliably
>> able to get natural joins to work.
>>
>> For example,
>>
>>> Tid <- c('AES 01-01-02 10:58:00', 'AES 01-01-02 11:53:00', 'AES 01-01-05
>> 10:58:00', 'AES 01-01-11 12:30:00')
>>> A <- data.frame(Tid, dfName = 'a')
>>> B <- data.frame(Tid = Tid[2:4], dfName = 'b')
>>> C <- data.frame(Tid = Tid[1:3], dfName = 'c')
>>
>> # then use the sqldf library
>>> library(sqldf)
>>> sqldf()
>>
>> # to create indices on the Tid variable shared across data.frames
>>> sqldf('create index indA on A(Tid)')
>>> sqldf('create index indB on B(Tid)')
>>> sqldf('create index indC on C(Tid)')
>>
>> # check to make sure everything is there
>>> sqldf('select * from sqlite_master')
>>
>> # doing a natural join (implicitly on Tid)
>> # does not give the expected joins
>>> sqldf('select * from main.A natural join main.B')
>> [1] Tid    dfName
>> <0 rows> (or 0-length row.names)
>>> sqldf('select * from main.A natural join main.C')
>> [1] Tid    dfName
>> <0 rows> (or 0-length row.names)
>>> sqldf('select * from main.B natural join main.C')
>> [1] Tid    dfName
>> <0 rows> (or 0-length row.names)
>>
>> # even using a where clause (which doesn't have the efficiency qualities I
>> need the indexed natural joins for) is problematic, setting values of the
>> dfName variable incorrectly for the data from C
>>> sqldf('select * from main.B b, main.C c where b.Tid = c.Tid')
>>                    Tid dfName                   Tid dfName
>> 1 AES 01-01-02 11:53:00      b AES 01-01-02 11:53:00      b
>> 2 AES 01-01-05 10:58:00      b AES 01-01-05 10:58:00      b
>>
>> I'm grateful for your guidance on what I'm doing wrong with the natural join
>> in sqldf.
>>
>> many thanks,
>> Nick
>

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

2010-05-20 Thread Thomas Stewart
I know this was asked and answered in 2005, I just want to make sure the
answer still holds in 2010.

Question:

If I have the mixed-model with

lmer( y ~ x1 + x2 + (1 | id) + (z1 | w ) , ... )

there is no way for me to specify a correlation structure for the random
effects? Right?

If I want to specify correlation structure then I need to use lme in the
nlme package?  Right?  Crossed Random effects are coded with the pdBlocked
syntax?  Right?

Thanks,
-tgs

[[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] Strange behaviour when using diff with POSIXt and POSIXlt objects

2010-05-20 Thread jim holtman
Please provide information as to what version you are using;  works fine for me:

> time3=strptime("2009 06 01 00 47 00",format="%Y %m %d %H %M")
> time4=strptime("2009 06 01 00 57 00",format="%Y %m %d %H %M")
>
> diff(c(time3,time4))
Time difference of 10 mins
>

I have version 2.10.1

On Thu, May 20, 2010 at 12:36 PM, Julian Burgos  wrote:
> Dear list,
>
> I´m calculating time differences between series of time stamps and I noticed
> something odd:
>
> If I do this...
>
>> time1=strptime("2009 05 31 22 57 00",format="%Y %m %d %H %M")
>> time2=strptime("2009 05 31 23 07 00",format="%Y %m %d %H %M")
>>
>> diff(c(time1,time2),units="mins")
> Time difference of 10 mins
>
> .. I get the correct response in minutes.  But if I try the same thing with
> different values, say..
>
>> time3=strptime("2009 06 01 00 47 00",format="%Y %m %d %H %M")
>> time4=strptime("2009 06 01 00 57 00",format="%Y %m %d %H %M")
>>
>> diff(c(time3,time4))
> Time difference of NA secs
>
> ...which is not what I´m looking for. The difference should also be 10
> minutes.
>
> I burned a few neurons (and searched the documentation) and I cannot figure
> why this happens.  Any ideas?
>
> All the best,
>
> Julian
>
> Julian Mariano Burgos
> Hafrannsóknastofnunin/Marine Research Institute
> Skúlagata 4, 121 Reykjavík, Iceland
> Sími/Telephone : +354-5752037
> Bréfsími/Telefax:  +354-5752001
> Netfang/Email: jul...@hafro.is, jmbur...@uw.edu
>
>        [[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
Cincinnati, OH
+1 513 646 9390

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.


  1   2   >