Re: [R] basic inquiry regarding write.csv

2009-02-17 Thread CJ Rubio

this is the code that i have so far:

> data.path <- file.path ("D:/documents/research/5 stations")
> setwd(data.path)
> getwd()

> data <- dir(".")
> num.files <- length(data)

>  for(i in 1:num.files) {
+  station.id <- substring(data[i], 1,8)
+  DF <- read.table(data[i], sep=",", blank.lines.skip = TRUE)
+  z <- zoo(DF[,4], as.Date(as.character(DF[,3]), "%m/%d/%Y"))
+  f <- function(x) time(x) [which.max(x)]
+  ix <- tapply(z, floor(as.yearmon(time(z))),f)
+  year <- (1988:2005)
+  date <- time(z[ix])
+  max.discharge <- coredata(z[ix])
+  data.info <- rbind(data.info, cbind(station.id, year, date,
max.discharge))
+  y <- split(data.info, data.info[station.id])
+  for (i in names(y)) {write.csv(y[[i]], file=paste(i, ".csv", sep=","))}
+ }

i have done the suggestion you gave me..  and there were no errors when i
run it.
hope this would further clarify my question.



So far, I cannot see any mistake, though the sep="" will be more
elegant that sep=",". Are you sure your working directory is
"E:/my_work_directory"? and is there any error msg?

BTW, a reproducible example will help to get better response from the list.

2009/2/18 CJ Rubio :
>
> thanks for your reply.. is there something wrong with the code i have?
> because it doesn't write the file in the directory that i am using...
>
> for (i in names(y))
>> {write.csv(y[[i]], file=paste(i, ".csv", sep=","))}
>
> thanks again..
>
>
>
> ronggui-2 wrote:
>>
>> If the "file" is a relative path, then it should be in the working
>> directory. Say, the working directory is E:/my_work_directory (of
>> course, you can get it by getwd()), and you export a data frame "a" to
>> csv by:
>> write.csv(a, file="a.csv"), then the file should be
>> "E:/my_work_directory/a.csv".
>>
>> Best
>>
>> 2009/2/18 CJ Rubio :
>>>
>>> i have a loop which looks likes this:
>>>
>>>
 data.info <- rbind(data.info, cbind(station.id, year, date,
 max.discharge))
>>> +  y <- split(data.info, data.info[station.id])
>>> +  for (i in names(y))
>>> {write.csv(y[[i]], file=paste(i, ".csv", sep=","))}
>>>
>>> i am wondering, where the file (which i am about to write in .csv
>>> format)
>>> will be saved? i looked at ?write.csv and it says there that :
>>>
>>> file   either a character string naming a file or a
>>> connection open for writing. "" indicates output to the console.
>>>
>>> correct me if i'm wrong but the way i undestand it is, i should have a
>>> file
>>> or a working directory where the .csv will be written.  if for example i
>>> have a working directory "E:/my_work_directory", how can i save this
>>> splitted files in the same directory?
>>>
>>> can anybody please enlighten me more with write.csv and the argument
>>> file?
>>>
>>> thanks..
>>>
>>> --
>>> View this message in context:
>>> http://www.nabble.com/basic-inquiry-regarding-write.csv-tp22073053p22073053.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.
>>>
>>
>>
>>
>> --
>> HUANG Ronggui, Wincent
>> Tel: (00852) 3442 3832
>> PhD Candidate
>> Dept of Public and Social Administration
>> City University of Hong Kong
>> Homepage: http://ronggui.huang.googlepages.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.
>>
>>
>
> --
> View this message in context:
> http://www.nabble.com/basic-inquiry-regarding-write.csv-tp22073053p22073393.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.
>



-- 
HUANG Ronggui, Wincent
Tel: (00852) 3442 3832
PhD Candidate
Dept of Public and Social Administration
City University of Hong Kong
Homepage: http://ronggui.huang.googlepages.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.



-- 
View this message in context: 
http://www.nabble.com/basic-inquiry-regarding-write.csv-tp22073053p22073615.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 

Re: [R] basic inquiry regarding write.csv

2009-02-17 Thread ronggui
So far, I cannot see any mistake, though the sep="" will be more
elegant that sep=",". Are you sure your working directory is
"E:/my_work_directory"? and is there any error msg?

BTW, a reproducible example will help to get better response from the list.

2009/2/18 CJ Rubio :
>
> thanks for your reply.. is there something wrong with the code i have?
> because it doesn't write the file in the directory that i am using...
>
> for (i in names(y))
>> {write.csv(y[[i]], file=paste(i, ".csv", sep=","))}
>
> thanks again..
>
>
>
> ronggui-2 wrote:
>>
>> If the "file" is a relative path, then it should be in the working
>> directory. Say, the working directory is E:/my_work_directory (of
>> course, you can get it by getwd()), and you export a data frame "a" to
>> csv by:
>> write.csv(a, file="a.csv"), then the file should be
>> "E:/my_work_directory/a.csv".
>>
>> Best
>>
>> 2009/2/18 CJ Rubio :
>>>
>>> i have a loop which looks likes this:
>>>
>>>
 data.info <- rbind(data.info, cbind(station.id, year, date,
 max.discharge))
>>> +  y <- split(data.info, data.info[station.id])
>>> +  for (i in names(y))
>>> {write.csv(y[[i]], file=paste(i, ".csv", sep=","))}
>>>
>>> i am wondering, where the file (which i am about to write in .csv format)
>>> will be saved? i looked at ?write.csv and it says there that :
>>>
>>> file   either a character string naming a file or a
>>> connection open for writing. "" indicates output to the console.
>>>
>>> correct me if i'm wrong but the way i undestand it is, i should have a
>>> file
>>> or a working directory where the .csv will be written.  if for example i
>>> have a working directory "E:/my_work_directory", how can i save this
>>> splitted files in the same directory?
>>>
>>> can anybody please enlighten me more with write.csv and the argument
>>> file?
>>>
>>> thanks..
>>>
>>> --
>>> View this message in context:
>>> http://www.nabble.com/basic-inquiry-regarding-write.csv-tp22073053p22073053.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.
>>>
>>
>>
>>
>> --
>> HUANG Ronggui, Wincent
>> Tel: (00852) 3442 3832
>> PhD Candidate
>> Dept of Public and Social Administration
>> City University of Hong Kong
>> Homepage: http://ronggui.huang.googlepages.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.
>>
>>
>
> --
> View this message in context: 
> http://www.nabble.com/basic-inquiry-regarding-write.csv-tp22073053p22073393.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.
>



-- 
HUANG Ronggui, Wincent
Tel: (00852) 3442 3832
PhD Candidate
Dept of Public and Social Administration
City University of Hong Kong
Homepage: http://ronggui.huang.googlepages.com/

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


[R] using stepAIC with negative binomial regression - error message help

2009-02-17 Thread t c
Dear List,
I am having problems running stepAIC with a negative binomial regression 
model.  I am working with data on manta ray abundance, using 20 predictor 
variables.  Predictors include variables for location (site), time (year, cos 
and sin of calendar day, length of day, percent lunar illumination), 
oceanography (sea surface temp mean and std, sea surface height mean and std), 
weather (cos wind direction, sin wind direction, wind speed, temperature, 
barometric pressure), and tides (dummy variables for high, falling, low, and 
rising).  So predictors are binary, categorical, and continuous variables.  I 
ran glm.nb to fit my full model containing all first order terms.  The model I 
ran was:
glm.nb.full stepAIC(glm.nb.full)
Start:  AIC=19240.46
mantas ~ site + year + cosday + sinday + daylength + lunarpercent + 
    sstmean + sststd + sshmean + sshstd + cosdir + sindir + spd + 
    temp + alt + tideht + high + falling + low + plankton
 
Error in dropterm.default(object, ...) : 
  number of rows in use has changed: remove missing values?
 
 
 
I know some of my predictors are missing values for certain dates.  I assume 
this is what “number of rows in use has changed” means.  Must I remove all rows 
that are missing values?  Or is there an option that does this for me?  Or am I 
completely off and there is something else going on?  This is my first time 
using R for statistics, so I am not sure if this is a mistake with my data or 
with my use of R.  I would appreciate any insight into what I am doing wrong.
 
Thanks,
 
Tim
 
Tim Clark
PhD Candidate
Department of Zoology
University of Hawaii
Honolulu, HI  96816


  
[[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] basic inquiry regarding write.csv

2009-02-17 Thread CJ Rubio

thanks for your reply.. is there something wrong with the code i have?
because it doesn't write the file in the directory that i am using...

for (i in names(y))
> {write.csv(y[[i]], file=paste(i, ".csv", sep=","))}

thanks again..



ronggui-2 wrote:
> 
> If the "file" is a relative path, then it should be in the working
> directory. Say, the working directory is E:/my_work_directory (of
> course, you can get it by getwd()), and you export a data frame "a" to
> csv by:
> write.csv(a, file="a.csv"), then the file should be
> "E:/my_work_directory/a.csv".
> 
> Best
> 
> 2009/2/18 CJ Rubio :
>>
>> i have a loop which looks likes this:
>>
>>
>>> data.info <- rbind(data.info, cbind(station.id, year, date,
>>> max.discharge))
>> +  y <- split(data.info, data.info[station.id])
>> +  for (i in names(y))
>> {write.csv(y[[i]], file=paste(i, ".csv", sep=","))}
>>
>> i am wondering, where the file (which i am about to write in .csv format)
>> will be saved? i looked at ?write.csv and it says there that :
>>
>> file   either a character string naming a file or a
>> connection open for writing. "" indicates output to the console.
>>
>> correct me if i'm wrong but the way i undestand it is, i should have a
>> file
>> or a working directory where the .csv will be written.  if for example i
>> have a working directory "E:/my_work_directory", how can i save this
>> splitted files in the same directory?
>>
>> can anybody please enlighten me more with write.csv and the argument
>> file?
>>
>> thanks..
>>
>> --
>> View this message in context:
>> http://www.nabble.com/basic-inquiry-regarding-write.csv-tp22073053p22073053.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.
>>
> 
> 
> 
> -- 
> HUANG Ronggui, Wincent
> Tel: (00852) 3442 3832
> PhD Candidate
> Dept of Public and Social Administration
> City University of Hong Kong
> Homepage: http://ronggui.huang.googlepages.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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/basic-inquiry-regarding-write.csv-tp22073053p22073393.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] basic inquiry regarding write.csv

2009-02-17 Thread ronggui
If the "file" is a relative path, then it should be in the working
directory. Say, the working directory is E:/my_work_directory (of
course, you can get it by getwd()), and you export a data frame "a" to
csv by:
write.csv(a, file="a.csv"), then the file should be
"E:/my_work_directory/a.csv".

Best

2009/2/18 CJ Rubio :
>
> i have a loop which looks likes this:
>
>
>> data.info <- rbind(data.info, cbind(station.id, year, date,
>> max.discharge))
> +  y <- split(data.info, data.info[station.id])
> +  for (i in names(y))
> {write.csv(y[[i]], file=paste(i, ".csv", sep=","))}
>
> i am wondering, where the file (which i am about to write in .csv format)
> will be saved? i looked at ?write.csv and it says there that :
>
> file   either a character string naming a file or a
> connection open for writing. "" indicates output to the console.
>
> correct me if i'm wrong but the way i undestand it is, i should have a file
> or a working directory where the .csv will be written.  if for example i
> have a working directory "E:/my_work_directory", how can i save this
> splitted files in the same directory?
>
> can anybody please enlighten me more with write.csv and the argument file?
>
> thanks..
>
> --
> View this message in context: 
> http://www.nabble.com/basic-inquiry-regarding-write.csv-tp22073053p22073053.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.
>



-- 
HUANG Ronggui, Wincent
Tel: (00852) 3442 3832
PhD Candidate
Dept of Public and Social Administration
City University of Hong Kong
Homepage: http://ronggui.huang.googlepages.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] Fwd: Percentiles/Quantiles with Weighting

2009-02-17 Thread Dieter Menne
Stavros Macrakis  alum.mit.edu> writes:

> 
> Some minor improvements and corrections below
> 
> # Simple weighted quantile
> #
> # v  A vector of sortable observations
> # w A numeric vector of positive weights
> # p  The quantile 0<=p<=1
> #
> # Nothing fancy: no interpolation etc.; NA cases not thought through
> 
>  wquantile <- function(v,w=rep(1,length(v)),p=.5)

You could compare it with wtd.quantile in Hmisc.

Dieter

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

2009-02-17 Thread CJ Rubio

i have a loop which looks likes this:


> data.info <- rbind(data.info, cbind(station.id, year, date,
> max.discharge))
+  y <- split(data.info, data.info[station.id])
+  for (i in names(y))
 {write.csv(y[[i]], file=paste(i, ".csv", sep=","))}

i am wondering, where the file (which i am about to write in .csv format)
will be saved? i looked at ?write.csv and it says there that :

file   either a character string naming a file or a
connection open for writing. "" indicates output to the console. 

correct me if i'm wrong but the way i undestand it is, i should have a file
or a working directory where the .csv will be written.  if for example i
have a working directory "E:/my_work_directory", how can i save this
splitted files in the same directory?

can anybody please enlighten me more with write.csv and the argument file?

thanks..

-- 
View this message in context: 
http://www.nabble.com/basic-inquiry-regarding-write.csv-tp22073053p22073053.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] Possible Cause of Segmentation Fault

2009-02-17 Thread Moumita Das
Hi All,

If you have already finished reading my previous emails regarding
segmentation fault , please have a look at this .I think this may help you
to diagnose the reason for the segmentation fault and help me,because i
don't understand much.

Rather than running the script using  the command "
source("new_regression.R") ", what I did was ,simply typed in the  commands
in R-prompt  and the results were:



> drv<-MySQL()

> drv



> dbConnect(drv, user="xyz", password="xyz",dbname =xyz_database, host =
xyz.com)

Error in mysqlNewConnection(drv, ...) : *object "xyz.com" not found*

* *

* *

Thanks

Moumita

[[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] Added system Info:--Segmentation Fault occured while connecting to the database

2009-02-17 Thread Moumita Das
Oops !! I had also included just one library in my script  i.e RMySQL.So
sorry for the inconvenience:(
-- 
Thanks
Moumita

[[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] Added system Info:--Segmentation Fault occured while connecting to the database

2009-02-17 Thread Moumita Das
Hi All,
Wanted to add some more information ,regarding my problem.
configuration of teh OS and R:---
Linux 2.6.18-6-686
> R.Version()
$platform
[1] "i486-pc-linux-gnu"

$arch
[1] "i486"

$os
[1] "linux-gnu"

$system
[1] "i486, linux-gnu"

$status
[1] "Patched"

$major
[1] "2"

$minor
[1] "4.0"

$year
[1] "2006"

$month
[1] "11"

$day
[1] "25"

$`svn rev`
[1] "39997"

$language
[1] "R"

$version.string
[1] "R version 2.4.0 Patched (2006-11-25 r39997)"



-- 
Thanks
Moumita

[[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] Segmentation Fault occured while connecting to the database

2009-02-17 Thread Moumita Das
Hi All,
Can anyone help me please?I don't know much about segmentation faults.I
understand what it is,but why my script's throwing the error i don't know.

This is my main function:

*main<-function()*

*{*

*dbName<-"xyz_database"*

*hostName<-"xyz.com"*

*con<-myDbconnect(dbName,hostName) *

**

*#Fetching exact sub group names*

*sub_grp_exact_num<-getSbGrpExactNum(con)*

*sub_grp_exact_num_data<-sub_grp_exact_num$matrix*

*sub_grp_exact_num_data_size<-sub_grp_exact_num$dim*

there's is some other code following the abovei thought
this much code is enough to explain my problem..error occurs when the
myDbconnect function is called at teh third line
**

*}*

myDbconnect  function is here:--

*myDbconnect<-function(dbName, hostName)*

*{*

*print("myDbconnect ")*

*drv<-MySQL()*

*#print(drv)*


*con <- dbConnect(drv, user="xyz", password="xyz",dbname =
dbName, host = hostName)*

*#return(con)*

*}*

When  print("myDbconnect print") is the first line of the myDbconnect
function  "myDbconnect print" gets printed.Though followed by  error
messages,shown below:-

*> source("new_regression.R")*

*Loading required package: DBI*

*[1] "myDbconnect print"*

* *** caught segfault 

*address 0x55, cause 'memory not mapped'*

* *

*Traceback:*

* 1: .Call("RS_MySQL_newConnection", drvId, con.params, groups,
default.file, PACKAGE = .MySQLPkgName)*

* 2: mysqlNewConnection(drv, ...)*

* 3: .class1(object)*

* 4: .class1(object)*

* 5: is(object, Cl)*

* 6: .valueClassTest(standardGeneric("dbConnect"), "DBIConnection",
"dbConnect")*

* 7: dbConnect(drv, user = "xyz", password = "xyz", dbname = dbName, host
= hostName)*

* 8: myDbconnect(dbName, hostName)*

* 9: main()*

*10: eval.with.vis(expr, envir, enclos)*

*11: eval.with.vis(ei, envir)*

*12: source("new_regression.R")*

* *

*Possible actions:*

*1: abort (with core dump)*

*2: normal R exit*

*3: exit R without saving workspace*

*4: exit R saving workspace*

*Selection:*

Now  if I change the position of the print statement  from first line to
fourth line after  dbConnect function:--

myDbconnect<-function(dbName, hostName)

{

drv<-MySQL()

#print(drv)

con <- dbConnect(drv, user="xyz", password="xyz",dbname =
dbName, host = hostName)

print("myDbconnect  print")

#return(con)

}

In this case the text  ("myDbconnect  print") is not printed ,otherwise what
I see is same as above i.e.:--

> source("new_regression.R")

Loading required package: DBI

*** caught segfault ***

address 0x55, cause 'memory not mapped'

Traceback:

 1: .Call("RS_MySQL_newConnection", drvId, con.params, groups,
default.file,
PACKAGE = .MySQLPkgName)

 2: mysqlNewConnection(drv, ...)

 3: .class1(object)

 4: .class1(object)

 5: is(object, Cl)

 6: .valueClassTest(standardGeneric("dbConnect"), "DBIConnection",
"dbConnect")

 7: dbConnect(drv, user = "xyz", password = "xyz", dbname = dbName, host
= hostName)

 8: myDbconnect(dbName, hostName)

 9: main()

10: eval.with.vis(expr, envir, enclos)

11: eval.with.vis(ei, envir)

12: source("new_regression.R")



Possible actions:

1: abort (with core dump)

2: normal R exit

3: exit R without saving workspace

4: exit R saving workspace

Selection:

*CONCLUSION:-- So the problem arises in the* *dbconnect function* *inside
myDbconnect  function.*

*i.e this line :--*

con <- dbConnect(drv, user="xyz", password="xyz",dbname = dbName, host =
hostName)

It's a segfault. A segmentation fault occurs when a program *attempts to
access a memory location that it is not allowed to access*, or *attempts to
access a memory location in a way that is not allowed* (for example,
attempting to write to a read-only location, or to overwrite part of the
operating system).Don't understand how I am trying to do any of these
incorrect operations.

Moreover  *traceback()  *is a function in R that is used for debugging ;
prints the call stack of the last uncaught error, i.e., the sequence of
calls that lead to the error.I have not used it anywhere in  my script
but  still
it shows the errors.



P.S:--- I  found these links,buti don't understand,where am I going wrong
:--

http://tolstoy.newcastle.edu.au/R/e2/devel/06/12/1436.html

http://www.ens.gu.edu.au/ROBERTK/R/HELP/00a/0111.HTML

http://tolstoy.newcastle.edu.au/R/e2/devel/06/12/1438.html


Any help is welcome..Even if anybody doesn't know to solve the problem,is
there anyother way ,i can get some more information,out of my script
regarding the cause of the error.do i need to make any changes,somewhere
else other than the script,to get rid of this problem.

-- 
Thanks in advance
Moumita


-- 
Thanks
Moumita

[[alternative HTML version deleted]]

___

[R] Index-G1 error

2009-02-17 Thread mauede
I am using some functions from package clusterSim to evaluate the best clusters 
layout.
Here is the features vector I am using to cluater 12 signals:

> alpha.vec
 [1] 0.8540039 0.8558350 0.8006592 0.8066406 0.8322754 0.8991699 0.8212891
 [8] 0.8815918 0.9050293 0.9174194 0.8613281 0.8425293

In the following I pasted an excerpt of my program:


library(clusterSim)

 setwd("C:/Documents and 
Settings/Monville/SpAn-Tests/Signals-Classification-Dir/Classification-by-Signal-Dir")

## LOAD DONOHO ALPHA FILE 
 features.tab <- read.table("Signals-Alpha-and-Zero-Moments.txt",header=T)

## EXTRACT SIGNALS ID
 insig <- features.tab[,"SignalID"]

## GENERATE DISTANCE MATRIX LABELS
 rwn <- NULL
 for(i in 1:length(insig)) {
   rwn <- c(rwn,as.character(insig[i]))
 }

## EXTRACT FEATURES VECTOR
 alpha.vec <- features.tab[,"Signal_Donoho_Alpha"]

## CALCULATE DISTANCE MATRIX - EUCLIDEAN DISTANCE
## ATTACH SIGNAL LABELS
 dist.mat <- dist(as.matrix(alpha.vec),method = "euclidean")
 yy <- as.matrix(dist.mat)
 rownames(yy) <- rwn

 nc.min <- 2#MINIMUM CLUSTERS NUM.
 nc.max <- length(alpha.vec) -1 #MAXIMUM CLUSTERS NUM.

## - PAM - INDEX.G1 
--
 res <- array(0,c(nc.max - nc.min +1,2))
 res[,1] <- nc.min:nc.max
 clusters <- NULL
 for (nc in nc.min:nc.max) {  
   cl <- pam(dist.mat,nc,diss=TRUE)
   res[nc-nc.min+1,2] <- G1 <- 
index.G1(as.matrix(alpha.vec),cl$cluster,d=dist.mat,centrotypes="medoids")
   clusters <- rbind(clusters, cl$cluster)
 }
###

I get the following error whenever I use  index.G1:

> nc <- 2
>  cl <- pam(dist.mat,nc,diss=TRUE)
> cl
Medoids:
 ID  
[1,]  5 5
[2,]  9 9
Clustering vector:
 [1] 1 1 1 1 1 2 1 2 2 2 1 1
Objective function:
 build   swap 
0.01873272 0.01620992 

Available components:
[1] "medoids""id.med" "clustering" "objective"  "isolation" 
[6] "clusinfo"   "silinfo""diss"   "call"  
> ndex.G1(as.matrix(alpha.vec),cl$cluster,d=dist.mat,centrotypes="medoids")
Error: could not find function "ndex.G1"
> index.G1(as.matrix(alpha.vec),cl$cluster,d=dist.mat,centrotypes="medoids")
Error in centers[i, ] <- .medoid(x[cl == i, ], d[cl == i, cl == i]) : 
  number of items to replace is not a multiple of replacement length


Please, notice that the same code work if I call index.G2, index.G3, index.S
Whay  it doesn't work with index.G1 ?

Thank you very much.
Maura 








tutti i telefonini TIM!


[[alternative HTML version deleted]]

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


[R] How to: C# / R interface

2009-02-17 Thread Baggett, Jonathan W
This is my situation:
I have a significant amount of data, and need to send it in pieces to R. I need 
R to return certain parameters for further use. 
I am sending files from C# (that are being queried from a database) into R.  
Currently I am trying to use the R(D)-Com package to figure out how to do this. 
 Along with sending in the file(s) to R, I will be sending the expressions for 
calculating distribution parameters in R as well.  I do not really know how to 
go about doing this, and there doesn't appear to be much literature online.  
Please let me know if you have any suggestions, sample code, references, etc.   

Thank you,

Jonathan Baggett
Georgia Institute of Technology
Industrial and Systems Engineering
---

"The price of freedom is constant vigilance." -Patrick Henry

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

2009-02-17 Thread Stavros Macrakis
Some minor improvements and corrections below

# Simple weighted quantile
#
# v  A vector of sortable observations
# w A numeric vector of positive weights
# p  The quantile 0<=p<=1
#
# Nothing fancy: no interpolation etc.; NA cases not thought through

 wquantile <- function(v,w=rep(1,length(v)),p=.5)
  {
if ( !is.numeric(w) || length(v) != length(w) )
  stop("Values and weights must be equal-length vectors")
if ( !is.numeric(p) || any( p<0 | p>1) )
  stop("Quantiles must be 0<=p<=1")
if ( min(w) < 0 ) stop("Weights must be non-negative numbers")
ranking <- order(v)
sumw <- cumsum(w[ranking])
plist <- sumw / sumw[ length(sumw) ]
sapply(p, function(p) v [ ranking [ which.max( plist >= p ) ] ] )
  }

I would appreciate any comments people have on this -- whether
correctness, efficiency, style, 

 -s

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


Re: [R] Help with rgl

2009-02-17 Thread Yihui Xie
"Chinese extend a helping hand to Russians who happen to be in Brazil
about a package written in Germany," which gladdened an American.
Trotsky would be even more proud  -- and amazed!! :-)

Regards,
Yihui
--
Yihui Xie 
Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086
Mobile: +86-15810805877
Homepage: http://www.yihui.name
School of Statistics, Room 1037, Mingde Main Building,
Renmin University of China, Beijing, 100872, China



On Wed, Feb 18, 2009 at 11:07 AM, roger koenker
 wrote:
> Why I love R  [Number  6]:
>
> Chinese  extend a helping hand to Russians who happen to be in Brazil
> about a package written in Germany.   Trotsky would be proud  -- and amazed!
>
> url:www.econ.uiuc.edu/~rogerRoger Koenker
> email   rkoen...@uiuc.edu   Department of Economics
> vox:217-333-4558University of Illinois
> fax:217-244-6678Champaign, IL 61820
>
>
> On Feb 17, 2009, at 8:37 PM, Yihui Xie wrote:
>
>> (1) you'll need ImageMagick installed to use the command "convert" to
>> convert image sequences into GIF animations; see ?movie3d
>> (2) "viewport" is read only!! see ?open3d carefully
>>
>> Regards,
>> Yihui
>> --
>> Yihui Xie 
>> Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086
>> Mobile: +86-15810805877
>> Homepage: http://www.yihui.name
>> School of Statistics, Room 1037, Mingde Main Building,
>> Renmin University of China, Beijing, 100872, China
>>
>>
>>
>> On Tue, Feb 17, 2009 at 2:45 AM, Iuri Gavronski  wrote:
>>>
>>> Hi,
>>>
>>> I don't know much about the RGL package, and I have read the
>>> documentation and tried some parameters, with no luck... I would like
>>> to generate a movie from a 3D object (code below), where the vortex A
>>> is closer to the observer, and then the object rotates and the B
>>> vortex gets closer. I would like to capture this movie to a file.
>>>
>>> By the way, I am not being able to insert unicode text with text3d.
>>>
>>> rgl 0.82, R 2.8.1, Windows Vista.
>>>
>>> Any help would be appreciated.
>>>
>>> Code follows:
>>>
>>> library(rgl)
>>> open3d()
>>>
>>> coord.1=c(0,100,0)
>>> coord.2=c(100,100,0)
>>> coord.3=c(100,0,0)
>>> coord.4=c(0,0,0)
>>> coord.5=c(50,50,70)
>>>
>>> pyrcolor="red"
>>> triangles3d(rbind(coord.1,coord.4,coord.5),color=pyrcolor)
>>> triangles3d(rbind(coord.1,coord.2,coord.5),color=pyrcolor)
>>> triangles3d(rbind(coord.2,coord.3,coord.5),color=pyrcolor)
>>> triangles3d(rbind(coord.3,coord.4,coord.5),color=pyrcolor)
>>> quads3d(rbind(coord.1,coord.2,coord.3,coord.4),color=pyrcolor)
>>>
>>> vertices = LETTERS[1:5]
>>> text3d(coord.1,text=vertices[1],adj=1,color="blue")
>>> text3d(coord.2,text=vertices[2],adj=0,color="blue")
>>> text3d(coord.3,text=vertices[3],adj=0,color="blue")
>>> text3d(coord.4,text=vertices[4],adj=1,color="blue")
>>> text3d(coord.5,text=vertices[5],adj=0,color="blue")
>>>
>>> # couldn't make this work...
>>> #open3d(viewport=c(0,0,686,489))
>>> #par3d(zoom = 1.157625)
>>>
>>> filename = "piramide.png"
>>> rgl.snapshot(filename)
>>>
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/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] Plotting Binned Data

2009-02-17 Thread Gundala Viswanath
Dear all,

I have a binned data that looks like this:

> dat
  (-1,9]   (9,19]  (19,29]  (29,39]  (39,49]  (49,59]  (59,69]  (69,79]
10063374   79   1643443
 (79,89]  (89,99]
   62

I tried to plot a histogram overlayed with curve.
With the following snippet:

library(lattice)
pdf("myfile.pdf")

hist(dat)
lines(dat,col="red")
dev.off()
__ END__

However the number of bins it plot is not the same
with "dat" (i.e. 10 bins). It only gives 6 bins.
The actual output can be seen in this URL

http://docs.google.com/Doc?id=dcvdrfrh_5cm5qkchw

How can I fix the code so that it gives a exact plot
to "dat" above, with same number of bins and its
respective frequency?

- Gundala Viswanath
Jakarta - Indonesia

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


Re: [R] Help with rgl

2009-02-17 Thread roger koenker

Why I love R  [Number  6]:

Chinese  extend a helping hand to Russians who happen to be in Brazil
about a package written in Germany.   Trotsky would be proud  -- and  
amazed!


url:www.econ.uiuc.edu/~rogerRoger Koenker
email   rkoen...@uiuc.edu   Department of Economics
vox:217-333-4558University of Illinois
fax:217-244-6678Champaign, IL 61820


On Feb 17, 2009, at 8:37 PM, Yihui Xie wrote:


(1) you'll need ImageMagick installed to use the command "convert" to
convert image sequences into GIF animations; see ?movie3d
(2) "viewport" is read only!! see ?open3d carefully

Regards,
Yihui
--
Yihui Xie 
Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086
Mobile: +86-15810805877
Homepage: http://www.yihui.name
School of Statistics, Room 1037, Mingde Main Building,
Renmin University of China, Beijing, 100872, China



On Tue, Feb 17, 2009 at 2:45 AM, Iuri Gavronski  wrote:

Hi,

I don't know much about the RGL package, and I have read the
documentation and tried some parameters, with no luck... I would like
to generate a movie from a 3D object (code below), where the vortex A
is closer to the observer, and then the object rotates and the B
vortex gets closer. I would like to capture this movie to a file.

By the way, I am not being able to insert unicode text with text3d.

rgl 0.82, R 2.8.1, Windows Vista.

Any help would be appreciated.

Code follows:

library(rgl)
open3d()

coord.1=c(0,100,0)
coord.2=c(100,100,0)
coord.3=c(100,0,0)
coord.4=c(0,0,0)
coord.5=c(50,50,70)

pyrcolor="red"
triangles3d(rbind(coord.1,coord.4,coord.5),color=pyrcolor)
triangles3d(rbind(coord.1,coord.2,coord.5),color=pyrcolor)
triangles3d(rbind(coord.2,coord.3,coord.5),color=pyrcolor)
triangles3d(rbind(coord.3,coord.4,coord.5),color=pyrcolor)
quads3d(rbind(coord.1,coord.2,coord.3,coord.4),color=pyrcolor)

vertices = LETTERS[1:5]
text3d(coord.1,text=vertices[1],adj=1,color="blue")
text3d(coord.2,text=vertices[2],adj=0,color="blue")
text3d(coord.3,text=vertices[3],adj=0,color="blue")
text3d(coord.4,text=vertices[4],adj=1,color="blue")
text3d(coord.5,text=vertices[5],adj=0,color="blue")

# couldn't make this work...
#open3d(viewport=c(0,0,686,489))
#par3d(zoom = 1.157625)

filename = "piramide.png"
rgl.snapshot(filename)



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


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


Re: [R] How to connect R and WinBUGS/OpenBUGS/LinBUGS in Linux in Feb. 2009

2009-02-17 Thread Paul Heinrich Dietrich

Hi Uwe,
Thank you for your guidance.  I have installed R2WinBUGS and WinBUGS14 under
wine.  Using ?bugs for help, it tells me:

useWINE: logical; attempt to use the Wine emulator to run 'WinBUGS',
defaults to 'FALSE' on Windows, and 'TRUE' otherwise. Not available in
S-PLUS.

WINE: character, path to 'wine' binary file, it is tried hard (by a guess
and the utilities 'which' and 'locate')  to get the information
automatically if not given.

newWINE: Use new versions of Wine that have 'winepath' utility

WINEPATH: character, path to 'winepath' binary file, it is tried hard (by a
guess and the utilities 'which' and 'locate')  to get the information
automatically if not given.

..and the following code is a simple Bayesian version of a t-test...

  Directory Paths  
MyModelPath <- "/home/me/Compound/R/WinBUGS/"
MyBUGSPath <- "/home/me/.wine/drive_c/Program Files/WinBUGS14/"
MyModelFile <- paste(MyModelPath, "model.bug", sep="")
WINEPATH <- "/usr/bin/wine"

  Create Data Set  
# Here is some fake data
n_draws <- 50
x <- round(runif(n_draws, 1, 2))
y <- ifelse(x == 1, rnorm(n_draws, 1, 1), rnorm(n_draws, 1.2, 0.8))
MyData <- as.data.frame(cbind(y, x))
y.n <- NROW(MyData$y)
x.j <- length(unique(x))
summary(MyData)

##  Format Data for WinBUGS  ##
MyBUGSData <- list(y=MyData$y, x=MyData$x, n=y.n, x.j=x.j)
MyBUGSData

##  WinBUGS Model File  ###
library(R2WinBUGS)
cat("model
{
for (i in 1:n)
{
y[i] ~ dnorm(mu[i], tau)
mu[i] <- alpha + beta[x[i]]
}
### STZ (Sum-To-Zero) Constraints
beta[1] <- -sum(beta[2:x.j])
### Priors
alpha ~ dnorm(0.0, 1.0E-4)
for (i in 2:x.j)
{
beta[i] ~ dnorm(0.0, 1.0E-4)
}
tau ~ dgamma(0.01, 0.01)
precision <- sqrt(1/tau)
}",
file=MyModelFile)
file.show(MyModelFile)

#  WinBUGS Model  #
MyModel <- bugs(MyBUGSData, inits=NULL,
model.file=MyModelFile, 
parameters.to.save=c("alpha", "beta", "precision"), 
n.chains=3, n.iter=2000, n.burnin=1000, n.thin=1, codaPkg=TRUE,
bugs.directory = MyBUGSPath, working.directory=MyModelPath, 
useWINE=TRUE, WINEPATH=WINEPATH, debug=TRUE)

The output says:

ERROR:  
  cannot open the connection

I'm wondering if I've misinterpreted how to set my paths with wine, because
I can go to the following path, double-click on WinBUGS14.exe, and open it
just fine: /home/me/.wine/drive_c/Program Files/WinBUGS14/

I can also go to Applications > Wine > Browse C:\ Drive and navigate to
WinBUGS.

Please help if I've done something wrong.  Thanks.
-- 
View this message in context: 
http://www.nabble.com/How-to-connect-R-and-WinBUGS-OpenBUGS-LinBUGS-in-Linux-in-Feb.-2009-tp22058716p22069602.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] FW: Can't access CRAN

2009-02-17 Thread Gabor Grothendieck
Maybe the suggestions here help:

http://cran.r-project.org/bin/windows/base/rw-FAQ.html#The-Internet-download-functions-fail_002e

On Tue, Feb 17, 2009 at 8:19 PM, Charles Annis, P.E.
 wrote:
> Dear R-Helpers:
>
> I'm running R version 2.8.1 on a 1 year old HP Pavilion with a AMD Athelon
> 64 Dual Core Processor 4800+ 2.50GHz chip, Vista Home Premium SP1, and 2 GB
> RAM.  I can't add packages (even running as Administrator) because I can't
> connect to a CRAN mirror.
>
>> chooseCRANmirror()
>
> (After nearly a minute I see the menu, and after any selection I immediately
> see)
>
> Warning message:
> In open.connection(con, "r") :
>  unable to connect to 'cran.r-project.org' on port 80.
>>
>
> BUT - I also have R 2.6.2 installed on this machine.  It works perfectly.
>
> Consulting the FAQ I tried changing permissions but no-go, plus the
> permissions on R 2.6.2 and R 2.8.1 appear to be identical.
>
> I'm baffled.  Can anyone illuminate this situation?
>
> Thanks
>
> Charles Annis, P.E.
>
> charles.an...@statisticalengineering.com
> phone: 561-352-9699
> eFax:  614-455-3265
> http://www.StatisticalEngineering.com
>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 randomly eliminate half the entries in a vector?

2009-02-17 Thread Esmail Bonakdarian

Gene Leynes wrote:

This is my first help post, hope it works!

Just check out the "sample" function
At the command line type:
?sample

I think it will be pretty clear from the documentation.


Yes, most excellent suggestion and quite helpful!

Thanks,
Esmail

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 randomly eliminate half the entries in a vector?

2009-02-17 Thread Gene Leynes
This is my first help post, hope it works!

Just check out the "sample" function
At the command line type:
?sample

I think it will be pretty clear from the documentation.

On Tue, Feb 17, 2009 at 9:13 PM, Esmail Bonakdarian wrote:

> (sorry if this is a duplicate-problems with posting at my end)
> 
> Hello all,
>
> I need some help with a nice R-idiomatic and efficient solution to a
> small problem.
>
> Essentially, I am trying to eliminate randomly half of the entries in
> a vector that contains index values into some other vectors.
>
> More details:
>
> I am working with two strings/vectors of 0s and 1s. These will contain
> about 200 elements (always the same number for both)
>
> I want to:
>
> 1. determines the locations of where the two strings differ
>
>--> easy using xor(s1, s2)
>
> 2. *randomly* selects *half* of those positions
>
>--> not sure how to do this. I suppose the result would be
>a list of index positions of size sum(xor(s1, s2))/2
>
> 3. exchange (flip) the bits in those random positions for both strings
>
>--> I have something that seems to do that, but it doesn't look
>slick and I wonder how efficient it is.
>
> Mostly I need help for #2, but will happily accept suggestions for #3,
> or for that matter anything that looks odd.
>
> Below my partial solution .. the HUX function is what I am trying
> to finish if someone can point me in the right direction.
>
> Thanks
> Esmail
> --
>
> rm(list=ls())
>
> 
> # create a binary vector of size "len"
> #
> create_bin_Chromosome <- function(len)
> {
>   sample(0:1, len, replace=T)
> }
>
> 
> # HUX - half uniform crossover
> #
> # 1. determines the locations of where the two strings
> #differ (easy xor)
> #
> # 2. randomly selects half of those positions
> #
> # 3. exchanges (flips) the bits in those positions for
> #both
> #
> HUX <- function(b1, b2)
> {
>   # 1. find differing bits
>   r=xor(b1, b2)
>
>   # positions where bits differ
>   different = which(r==TRUE)
>
>   cat("\nhrp: ", different, "\n")
>   # 2. ??? how to do this best so that each time
>   #a different half subset is selected? I.e.,
>   #sum(r)/2 positions.
>
>   # 3. this flips *all* positions, should really only flip
>   #half of them (randomly selected half)
>   new_b1 = b1
>   new_b2 = b2
>
>   for(i in different)  # should contain half the entries (randomly)
>   {
> new_b1[i] = b2[i]
> new_b2[i] = b1[i]
>   }
>
>   result <- matrix(c(new_b1, new_b2), 2, LEN, byrow=T)
>   result
> }
>
> LEN = 5
> b1=create_bin_Chromosome(LEN)
> b2=create_bin_Chromosome(LEN)
>
> cat(b1, "\n")
> cat(b2, "\n")
>
> idx=HUX(b1, b2)
> cat("\n\n")
> cat(idx[1,], "\n")
> cat(idx[2,], "\n")
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[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] rbind: number of columns of result is not a multiple of vector length (arg 1)

2009-02-17 Thread jim holtman
You could 'split' the dataframe and then write out each element: something like

x <- split(data.info, data.info$station.id)
for (i in names(x)) write.csv(x[[i]], file=paste(i, ".csv", sep=""))

On Tue, Feb 17, 2009 at 10:11 PM, CJ Rubio  wrote:
>
> it works perfectly thank you for your help
>
> what if i want to seperate each stations data and save in a .csv file??
> while i'm asking you these i'm also finding some ways to do so..
>
> thank you again.
>
>
>
>
>
> jholtman wrote:
>>
>> try using:
>>
>> data.info <- rbind(data.info, cbind(station.id, year, date,
>> max.discharge))
>>
>> On Tue, Feb 17, 2009 at 9:26 PM, CJ Rubio  wrote:
>>>
>>> i have the following constructed and running very well,, thanks to Gabor
>>> Grothendieck for his help.
>>>
data.info <- c("station.id", "year", "date", "max.discharge")

 for(i in 1:num.files) {
>>> + station.id <- substring(data[i], 1,8)
>>> + DF <- read.table(data[i], sep=",", blank.lines.skip = TRUE)
>>> + z <- zoo(DF[,4], as.Date(as.character(DF[,3]), "%m/%d/%Y"))
>>> + f <- function(x) time(x) [which.max(x)]
>>> + ix <- tapply(z, floor(as.yearmon(time(z))),f)
>>> + year <- (1988:2005)
>>> + date <- time(z[ix])
>>> + max.discharge <- coredata(z[ix])
>>> + data.info <- rbind(data.info, c(station.id, year, date, max.discharge))
>>> + }
>>>
>>> my problem with my code occurs in the part where I arrange my results..
>>> after running this code, i get this warning:
>>>
>>> Warning message:
>>> In rbind(data.info, c(station.id, year, date, max.discharge)) :
>>>  number of columns of result is not a multiple of vector length (arg 1)
>>>
>>>
>>> i can't figure out what to do to produce the result i wanted:
>>> (for each station, it should look like this:)
>>>
>>> data.info "station.id" "year"   "date"  "max.discharge"
>>>   "01014000" 1988   "1988-11-07"   4360
>>>   "01014000" 1989   "1989-05-13" 2
>>>   "01014000" 1990   "1990-10-25"   9170
>>>   "01014000" 1991   "1991-04-22" 12200
>>>   "01014000" 1992   "1992-03-29" 11800
>>>
>>>   "01014000" 2005   "2005-04-04" 22100
>>>
>>> thanks in advence for your help..
>>> --
>>> View this message in context:
>>> http://www.nabble.com/rbind%3A-number-of-columns-of-result-is-not-a-multiple-of-vector-length-%28arg-1%29-tp22070942p22070942.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.
>>>
>>
>>
>>
>> --
>> 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.
>>
>>
>
> --
> View this message in context: 
> http://www.nabble.com/rbind%3A-number-of-columns-of-result-is-not-a-multiple-of-vector-length-%28arg-1%29-tp22070942p22071324.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.
>



-- 
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] how to randomly eliminate half the entries in a vector?

2009-02-17 Thread Esmail Bonakdarian
(sorry if this is a duplicate-problems with posting at my end)

Hello all,

I need some help with a nice R-idiomatic and efficient solution to a
small problem.

Essentially, I am trying to eliminate randomly half of the entries in
a vector that contains index values into some other vectors.

More details:

I am working with two strings/vectors of 0s and 1s. These will contain
about 200 elements (always the same number for both)

I want to:

1. determines the locations of where the two strings differ

--> easy using xor(s1, s2)

2. *randomly* selects *half* of those positions

--> not sure how to do this. I suppose the result would be
a list of index positions of size sum(xor(s1, s2))/2

3. exchange (flip) the bits in those random positions for both strings

--> I have something that seems to do that, but it doesn't look
slick and I wonder how efficient it is.

Mostly I need help for #2, but will happily accept suggestions for #3,
or for that matter anything that looks odd.

Below my partial solution .. the HUX function is what I am trying
to finish if someone can point me in the right direction.

Thanks
Esmail
--

rm(list=ls())


# create a binary vector of size "len"
#
create_bin_Chromosome <- function(len)
{
   sample(0:1, len, replace=T)
}


# HUX - half uniform crossover
#
# 1. determines the locations of where the two strings
#differ (easy xor)
#
# 2. randomly selects half of those positions
#
# 3. exchanges (flips) the bits in those positions for
#both
#
HUX <- function(b1, b2)
{
   # 1. find differing bits
   r=xor(b1, b2)

   # positions where bits differ
   different = which(r==TRUE)

   cat("\nhrp: ", different, "\n")
   # 2. ??? how to do this best so that each time
   #a different half subset is selected? I.e.,
   #sum(r)/2 positions.

   # 3. this flips *all* positions, should really only flip
   #half of them (randomly selected half)
   new_b1 = b1
   new_b2 = b2

   for(i in different)  # should contain half the entries (randomly)
   {
 new_b1[i] = b2[i]
 new_b2[i] = b1[i]
   }

   result <- matrix(c(new_b1, new_b2), 2, LEN, byrow=T)
   result
}

LEN = 5
b1=create_bin_Chromosome(LEN)
b2=create_bin_Chromosome(LEN)

cat(b1, "\n")
cat(b2, "\n")

idx=HUX(b1, b2)
cat("\n\n")
cat(idx[1,], "\n")
cat(idx[2,], "\n")

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


Re: [R] rbind: number of columns of result is not a multiple of vector length (arg 1)

2009-02-17 Thread CJ Rubio

it works perfectly thank you for your help

what if i want to seperate each stations data and save in a .csv file??
while i'm asking you these i'm also finding some ways to do so..

thank you again.





jholtman wrote:
> 
> try using:
> 
> data.info <- rbind(data.info, cbind(station.id, year, date,
> max.discharge))
> 
> On Tue, Feb 17, 2009 at 9:26 PM, CJ Rubio  wrote:
>>
>> i have the following constructed and running very well,, thanks to Gabor
>> Grothendieck for his help.
>>
>>>data.info <- c("station.id", "year", "date", "max.discharge")
>>>
>>> for(i in 1:num.files) {
>> + station.id <- substring(data[i], 1,8)
>> + DF <- read.table(data[i], sep=",", blank.lines.skip = TRUE)
>> + z <- zoo(DF[,4], as.Date(as.character(DF[,3]), "%m/%d/%Y"))
>> + f <- function(x) time(x) [which.max(x)]
>> + ix <- tapply(z, floor(as.yearmon(time(z))),f)
>> + year <- (1988:2005)
>> + date <- time(z[ix])
>> + max.discharge <- coredata(z[ix])
>> + data.info <- rbind(data.info, c(station.id, year, date, max.discharge))
>> + }
>>
>> my problem with my code occurs in the part where I arrange my results..
>> after running this code, i get this warning:
>>
>> Warning message:
>> In rbind(data.info, c(station.id, year, date, max.discharge)) :
>>  number of columns of result is not a multiple of vector length (arg 1)
>>
>>
>> i can't figure out what to do to produce the result i wanted:
>> (for each station, it should look like this:)
>>
>> data.info "station.id" "year"   "date"  "max.discharge"
>>   "01014000" 1988   "1988-11-07"   4360
>>   "01014000" 1989   "1989-05-13" 2
>>   "01014000" 1990   "1990-10-25"   9170
>>   "01014000" 1991   "1991-04-22" 12200
>>   "01014000" 1992   "1992-03-29" 11800
>>
>>   "01014000" 2005   "2005-04-04" 22100
>>
>> thanks in advence for your help..
>> --
>> View this message in context:
>> http://www.nabble.com/rbind%3A-number-of-columns-of-result-is-not-a-multiple-of-vector-length-%28arg-1%29-tp22070942p22070942.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.
>>
> 
> 
> 
> -- 
> 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.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/rbind%3A-number-of-columns-of-result-is-not-a-multiple-of-vector-length-%28arg-1%29-tp22070942p22071324.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] rbind: number of columns of result is not a multiple of vector length (arg 1)

2009-02-17 Thread jim holtman
try using:

data.info <- rbind(data.info, cbind(station.id, year, date, max.discharge))

On Tue, Feb 17, 2009 at 9:26 PM, CJ Rubio  wrote:
>
> i have the following constructed and running very well,, thanks to Gabor
> Grothendieck for his help.
>
>>data.info <- c("station.id", "year", "date", "max.discharge")
>>
>> for(i in 1:num.files) {
> + station.id <- substring(data[i], 1,8)
> + DF <- read.table(data[i], sep=",", blank.lines.skip = TRUE)
> + z <- zoo(DF[,4], as.Date(as.character(DF[,3]), "%m/%d/%Y"))
> + f <- function(x) time(x) [which.max(x)]
> + ix <- tapply(z, floor(as.yearmon(time(z))),f)
> + year <- (1988:2005)
> + date <- time(z[ix])
> + max.discharge <- coredata(z[ix])
> + data.info <- rbind(data.info, c(station.id, year, date, max.discharge))
> + }
>
> my problem with my code occurs in the part where I arrange my results..
> after running this code, i get this warning:
>
> Warning message:
> In rbind(data.info, c(station.id, year, date, max.discharge)) :
>  number of columns of result is not a multiple of vector length (arg 1)
>
>
> i can't figure out what to do to produce the result i wanted:
> (for each station, it should look like this:)
>
> data.info "station.id" "year"   "date"  "max.discharge"
>   "01014000" 1988   "1988-11-07"   4360
>   "01014000" 1989   "1989-05-13" 2
>   "01014000" 1990   "1990-10-25"   9170
>   "01014000" 1991   "1991-04-22" 12200
>   "01014000" 1992   "1992-03-29" 11800
>
>   "01014000" 2005   "2005-04-04" 22100
>
> thanks in advence for your help..
> --
> View this message in context: 
> http://www.nabble.com/rbind%3A-number-of-columns-of-result-is-not-a-multiple-of-vector-length-%28arg-1%29-tp22070942p22070942.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.
>



-- 
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] Help with rgl

2009-02-17 Thread Yihui Xie
(1) you'll need ImageMagick installed to use the command "convert" to
convert image sequences into GIF animations; see ?movie3d
(2) "viewport" is read only!! see ?open3d carefully

Regards,
Yihui
--
Yihui Xie 
Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086
Mobile: +86-15810805877
Homepage: http://www.yihui.name
School of Statistics, Room 1037, Mingde Main Building,
Renmin University of China, Beijing, 100872, China



On Tue, Feb 17, 2009 at 2:45 AM, Iuri Gavronski  wrote:
> Hi,
>
> I don't know much about the RGL package, and I have read the
> documentation and tried some parameters, with no luck... I would like
> to generate a movie from a 3D object (code below), where the vortex A
> is closer to the observer, and then the object rotates and the B
> vortex gets closer. I would like to capture this movie to a file.
>
> By the way, I am not being able to insert unicode text with text3d.
>
> rgl 0.82, R 2.8.1, Windows Vista.
>
> Any help would be appreciated.
>
> Code follows:
>
> library(rgl)
> open3d()
>
> coord.1=c(0,100,0)
> coord.2=c(100,100,0)
> coord.3=c(100,0,0)
> coord.4=c(0,0,0)
> coord.5=c(50,50,70)
>
> pyrcolor="red"
> triangles3d(rbind(coord.1,coord.4,coord.5),color=pyrcolor)
> triangles3d(rbind(coord.1,coord.2,coord.5),color=pyrcolor)
> triangles3d(rbind(coord.2,coord.3,coord.5),color=pyrcolor)
> triangles3d(rbind(coord.3,coord.4,coord.5),color=pyrcolor)
> quads3d(rbind(coord.1,coord.2,coord.3,coord.4),color=pyrcolor)
>
> vertices = LETTERS[1:5]
> text3d(coord.1,text=vertices[1],adj=1,color="blue")
> text3d(coord.2,text=vertices[2],adj=0,color="blue")
> text3d(coord.3,text=vertices[3],adj=0,color="blue")
> text3d(coord.4,text=vertices[4],adj=1,color="blue")
> text3d(coord.5,text=vertices[5],adj=0,color="blue")
>
> # couldn't make this work...
> #open3d(viewport=c(0,0,686,489))
> #par3d(zoom = 1.157625)
>
> filename = "piramide.png"
> rgl.snapshot(filename)
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] rbind: number of columns of result is not a multiple of vector length (arg 1)

2009-02-17 Thread CJ Rubio

i have the following constructed and running very well,, thanks to Gabor
Grothendieck for his help.

>data.info <- c("station.id", "year", "date", "max.discharge")
>
> for(i in 1:num.files) {
+ station.id <- substring(data[i], 1,8)
+ DF <- read.table(data[i], sep=",", blank.lines.skip = TRUE)
+ z <- zoo(DF[,4], as.Date(as.character(DF[,3]), "%m/%d/%Y"))
+ f <- function(x) time(x) [which.max(x)]
+ ix <- tapply(z, floor(as.yearmon(time(z))),f)
+ year <- (1988:2005)
+ date <- time(z[ix])
+ max.discharge <- coredata(z[ix])
+ data.info <- rbind(data.info, c(station.id, year, date, max.discharge))
+ }

my problem with my code occurs in the part where I arrange my results..
after running this code, i get this warning:

Warning message:
In rbind(data.info, c(station.id, year, date, max.discharge)) :
  number of columns of result is not a multiple of vector length (arg 1)


i can't figure out what to do to produce the result i wanted:
(for each station, it should look like this:)

data.info "station.id" "year"   "date"  "max.discharge" 
   "01014000" 1988   "1988-11-07"   4360 
   "01014000" 1989   "1989-05-13" 2  
   "01014000" 1990   "1990-10-25"   9170
   "01014000" 1991   "1991-04-22" 12200
   "01014000" 1992   "1992-03-29" 11800

   "01014000" 2005   "2005-04-04" 22100

thanks in advence for your help..
-- 
View this message in context: 
http://www.nabble.com/rbind%3A-number-of-columns-of-result-is-not-a-multiple-of-vector-length-%28arg-1%29-tp22070942p22070942.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] how to randomly eliminate half the entries in a vector?

2009-02-17 Thread Jeremiah Rounds

Here is what I got for script through your third question:

 

set.seed(1)

x1 = rbinom(200,1,.5)

x2 = rbinom(200,1,.5)

differ = x1 != x2

differ.indexes = (1:length(x1))[differ == TRUE]

#you were unclear if you want to round up or round down on odd index of 
differ.indexes

n = floor( length(differ.indexes)/2)

#sampling without replacement

random.indexes = sample(differ.indexes,n )

swapping = x1[random.indexes] #with 1s and 0s you can do this without this 
variable.  

x1[random.indexes] = x2[random.indexes] 

x2[random.indexes] = swapping

 

 

Good luck,

Jeremiah
 
> Date: Tue, 17 Feb 2009 20:17:51 -0500
> From: esmail...@gmail.com
> To: r-help@r-project.org
> Subject: [R] how to randomly eliminate half the entries in a vector?
> 
> Hello all,
> 
> I need some help with a nice R-idiomatic and efficient solution to a
> small problem.
> 
> Essentially, I am trying to eliminate randomly half of the entries in
> a vector that contains index values into some other vectors.
> 
> 
> More details:
> 
> I am working with two strings/vectors of 0s and 1s. These will contain
> about 200 elements (always the same number for both)
> 
> I want to:
> 
> 1. determines the locations of where the two strings differ
> 
> --> easy using xor(s1, s2)
> 
> 2. *randomly* selects *half* of those positions
> 
> --> not sure how to do this. I suppose the result would be
> a list of index positions of size sum(xor(s1, s2))/2
> 
> 
> 3. exchange (flip) the bits in those random positions for both strings
> 
> --> I have something that seems to do that, but it doesn't look
> slick and I wonder how efficient it is.
> 
> Mostly I need help for #2, but will happily accept suggestions for #3,
> or for that matter anything that looks odd.
> 
> Below my partial solution .. the HUX function is what I am trying
> to finish if someone can point me in the right direction.
> 
> Thanks
> Esmail
> --
> 
> 
> 
> rm(list=ls())
> 
> 
> # create a binary vector of size "len"
> #
> create_bin_Chromosome <- function(len)
> {
> sample(0:1, len, replace=T)
> }
> 
> 
> 
> 
> # HUX - half uniform crossover
> #
> # 1. determines the locations of where the two strings
> # differ (easy xor)
> #
> # 2. randomly selects half of those positions
> #
> # 3. exchanges (flips) the bits in those positions for
> # both
> #
> HUX <- function(b1, b2)
> {
> # 1. find differing bits
> r=xor(b1, b2)
> 
> # positions where bits differ
> different = which(r==TRUE)
> 
> cat("\nhrp: ", different, "\n")
> # 2. ??? how to do this best so that each time
> # a different half subset is selected? I.e.,
> # sum(r)/2 positions.
> 
> 
> # 3. this flips *all* positions, should really only flip
> # half of them (randomly selected half)
> new_b1 = b1
> new_b2 = b2
> 
> for(i in different) # should contain half the entries (randomly)
> {
> new_b1[i] = b2[i]
> new_b2[i] = b1[i]
> }
> 
> result <- matrix(c(new_b1, new_b2), 2, LEN, byrow=T)
> result
> }
> 
> 
> 
> LEN = 5
> b1=create_bin_Chromosome(LEN)
> b2=create_bin_Chromosome(LEN)
> 
> cat(b1, "\n")
> cat(b2, "\n")
> 
> idx=HUX(b1, b2)
> cat("\n\n")
> cat(idx[1,], "\n")
> cat(idx[2,], "\n")
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

_

 go.

[[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] cumsum vs. sum

2009-02-17 Thread Gabor Grothendieck
Check out sum.exact and cumsum.exact in the caTools package.

> library(caTools)
Loading required package: bitops
> x <- 1/(12:14)
> sum(x) - cumsum(x)[3]
[1] 2.775558e-17
> sum.exact(x) - cumsum.exact(x)[3]
[1] 0


On Tue, Feb 17, 2009 at 5:12 PM, Stavros Macrakis  wrote:
> I recently traced a bug of mine to the fact that cumsum(s)[length(s)]
> is not always exactly equal to sum(s).
>
> For example,
>
> x<-1/(12:14)
> sum(x) - cumsum(x)[3]  => 2.8e-17
>
> Floating-point addition is of course not exact, and in particular is
> not associative, so there are various possible reasons for this.
> Perhaps sum uses clever summing tricks to get more accurate results?
> In some quick experiments, it does seem to get more accurate results
> than cumsum.
>
> It might be worth documenting.
>
> -s
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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] FW: Can't access CRAN

2009-02-17 Thread Charles Annis, P.E.
Dear R-Helpers:

I'm running R version 2.8.1 on a 1 year old HP Pavilion with a AMD Athelon
64 Dual Core Processor 4800+ 2.50GHz chip, Vista Home Premium SP1, and 2 GB
RAM.  I can't add packages (even running as Administrator) because I can't
connect to a CRAN mirror.

> chooseCRANmirror()

(After nearly a minute I see the menu, and after any selection I immediately
see)

Warning message:
In open.connection(con, "r") : 
  unable to connect to 'cran.r-project.org' on port 80.
> 

BUT - I also have R 2.6.2 installed on this machine.  It works perfectly.

Consulting the FAQ I tried changing permissions but no-go, plus the
permissions on R 2.6.2 and R 2.8.1 appear to be identical.

I'm baffled.  Can anyone illuminate this situation?

Thanks

Charles Annis, P.E.

charles.an...@statisticalengineering.com
phone: 561-352-9699
eFax:  614-455-3265
http://www.StatisticalEngineering.com
 

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


[R] how to randomly eliminate half the entries in a vector?

2009-02-17 Thread Esmail Bonakdarian

Hello all,

I need some help with a nice R-idiomatic and efficient solution to a
small problem.

Essentially, I am trying to eliminate randomly half of the entries in
a vector that contains index values into some other vectors.


More details:

I am working with two strings/vectors of 0s and 1s. These will contain
about 200 elements (always the same number for both)

I want to:

1. determines the locations of where the two strings differ

   --> easy using xor(s1, s2)

2. *randomly* selects *half* of those positions

   --> not sure how to do this. I suppose the result would be
   a list of index positions of size sum(xor(s1, s2))/2


3. exchange (flip) the bits in those random positions for both strings

   --> I have something that seems to do that, but it doesn't look
   slick and I wonder how efficient it is.

Mostly I need help for #2, but will happily accept suggestions for #3,
or for that matter anything that looks odd.

Below my partial solution .. the HUX function is what I am trying
to finish if someone can point me in the right direction.

Thanks
Esmail
--



rm(list=ls())


# create a binary vector of size "len"
#
create_bin_Chromosome <- function(len)
{
  sample(0:1, len, replace=T)
}




# HUX - half uniform crossover
#
# 1. determines the locations of where the two strings
#differ (easy xor)
#
# 2. randomly selects half of those positions
#
# 3. exchanges (flips) the bits in those positions for
#both
#
HUX <- function(b1, b2)
{
  # 1. find differing bits
  r=xor(b1, b2)

  # positions where bits differ
  different = which(r==TRUE)

  cat("\nhrp: ", different, "\n")
  # 2. ??? how to do this best so that each time
  #a different half subset is selected? I.e.,
  #sum(r)/2 positions.


  # 3. this flips *all* positions, should really only flip
  #half of them (randomly selected half)
  new_b1 = b1
  new_b2 = b2

  for(i in different)  # should contain half the entries (randomly)
  {
new_b1[i] = b2[i]
new_b2[i] = b1[i]
  }

  result <- matrix(c(new_b1, new_b2), 2, LEN, byrow=T)
  result
}



LEN = 5
b1=create_bin_Chromosome(LEN)
b2=create_bin_Chromosome(LEN)

cat(b1, "\n")
cat(b2, "\n")

idx=HUX(b1, b2)
cat("\n\n")
cat(idx[1,], "\n")
cat(idx[2,], "\n")

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


Re: [R] ylim in plot.ts?

2009-02-17 Thread Gabor Grothendieck
Check out:

https://www.stat.math.ethz.ch/pipermail/r-help/2009-February/187562.html

On Tue, Feb 17, 2009 at 6:45 PM, Oliver Bandel
 wrote:
> Hello,
>
>
> I tried to use ylim=c(x,y) in a plot.ts().
>
> This was ignored.
>
> How can I achieve it to create such graphics?
>
>
> Ciao,
>   Oliver Bandel
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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


Re: [R] :How to insert a column in a data frame

2009-02-17 Thread Rolf Turner


On 18/02/2009, at 12:51 PM, Laura Rodriguez Murillo wrote:


Hi dear list,

I wonder if somebody can help me with this. I have a text file with
300 rows and around 30 columns and I need to insert a column that
has the number 1 in every row. This new column should be placed
between columns 6 and 7.

As an example: I would want to insert a column (consiting just of the
number 1 in every row, as in column 6) between columns 6 and 7.

847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A  
T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A  
T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A  
T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A  
T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A  
T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A  
T T


So the new table would look like this:

847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C  
A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C  
A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C  
A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C  
A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C  
A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C  
A T T


You can also look at this problem as duplicating column 6.


(1) y <- cbind(x[,1:6],clyde=rep(1,300),x[,7:ncol(x)])

(2) y <- as.data.frame(c(x[1:6],list(clyde=rep(1,300)),x[7:length(x)]))

where x is your original data stored as a matrix or data frame in  
method (1),
and as a data frame in method (2), and y is the new data set with the  
extra

column of 1's.

Or, if you really are just duplicating column 6:

y <- x[,c(1:6,6,7:ncol(x))]
names(y)[7] <- "clyde"# For consistency with other methods.

These work with your toy example.  I haven't tried them out with data  
sets of

the size of your real data set.

cheers,

Rolf Turner

##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

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


[R] Normal cdf modified function

2009-02-17 Thread Fernando Saldanha
I wonder if an R package would have a function that calculates the following.

Let Y be a normal multivariate function. For example, let Y have 4
dimensions. I want to calculate

P(Y1 < Z1, Y2 < Z2, Y3 > Z3, Y4 > Z4).

There are R functions to do the calculation if all the inequalities
are of the type "<" (the cdf). But is there an R function where the
two types of inequalities ("<" and ">") can be mixed? (The user would
have to specify the set of indexes with inequalities of the type ">")

Thanks for any suggestions.

FS

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


Re: [R] Python and R

2009-02-17 Thread Esmail Bonakdarian
On Tue, Feb 17, 2009 at 6:05 PM, Barry Rowlingson
 wrote:
> 2009/2/17 Esmail Bonakdarian :

>  When I need to use the two together, it's easiest with 'rpy'. This
> lets you call R functions from python, so you can do:
>
>  from rpy import r
>  r.hist(z)

wow .. that is pretty straight forward, I'll have to check out rpy for sure.

> to get a histogram of the values in a python list 'z'. There are some
> complications converting structured data types between the two but
> they can be overcome, and apparently are handled better with the next
> generation Rpy2 (which I've not got into yet). Google for rpy for
> info.

Will do!

>> Is there much of a performance hit either way? (as both are interpreted
>> languages)
>
>  Not sure what you mean here. Do you mean is:
>
>  R> sum(x)
>
> faster than
>
> Python> sum(x)
>
> and how much worse is:
>
> Python> from rpy import r
> Python> r.sum(x)
>

Well, I have a program written in R which already takes quite a while
to run. I was
just wondering if I were to rewrite most of the logic in Python - the
main thing I use
in R are its regression facilities - if it would speed things up. I
suspect not since
both of them are interpreted, and the bulk of the time is taken up by
R's regression
calls.

Esmail

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


[R] :How to insert a column in a data frame

2009-02-17 Thread Laura Rodriguez Murillo
Hi dear list,

I wonder if somebody can help me with this. I have a text file with
300 rows and around 30 columns and I need to insert a column that
has the number 1 in every row. This new column should be placed
between columns 6 and 7.

As an example: I would want to insert a column (consiting just of the
number 1 in every row, as in column 6) between columns 6 and 7.

847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T

So the new table would look like this:

847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T
847 847 0 0 2 1 1 G G T C A G T T C G A C C C G C A A A G G G G A C A T T

You can also look at this problem as duplicating column 6.

Thank you for your help!


Laura R Murillo
Columbia University
New York

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

2009-02-17 Thread Oliver Bandel

Hello,


I tried to use ylim=c(x,y) in a plot.ts().

This was ignored.

How can I achieve it to create such graphics?


Ciao,
   Oliver Bandel

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Saving data to a MS Access Database

2009-02-17 Thread Bill Cunliffe
I have the following dataframe:

 

ad <- data.frame(dates, av, sn$SectorName)

colnames(ad) <- c("Date", "Value", "Tag")

 

which has data (rows 10 to 20, for example) as follows:

 

 DateValue   Tag

10 2008-01-16-0.20875ConDisc

11 2008-01-17-0.21125ConDisc

12 2008-01-18-0.08875ConDisc

13 2008-01-21-0.10875ConDisc

14 2008-01-220.02125 ConDisc

15 2008-01-230.09500 ConDisc

16 2008-01-240.02500 ConDisc

17 2008-01-250.03500 ConDisc

18 2008-01-280.06625 ConDisc

19 2008-01-290.20500 ConDisc

20 2008-01-300.11625 ConDisc

 

I want to import this data in into an existing MS Access database which has
the following properties:

 

Name: tblAD

ID: Autonumber, primary key

Date: Date

Value: Double

Tag: String

 

I am trying this by using the following:

 

result <- sqlUpdate(channel, ad, "tblAD")

 

But I get the following error:

> result <- sqlUpdate(channel, ad, "tblAD")

Error in sqlUpdate(channel, ad, "tblAD") :

cannot update 'tblAD' without unique column

 

 

Can anyone suggest where I am going wrong?

 


[[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] Python and R

2009-02-17 Thread Esmail Bonakdarian
Hello!

On Tue, Feb 17, 2009 at 5:58 PM, Warren Young  wrote:
>
> Esmail Bonakdarian wrote:
>>
>> I am just wondering if any of you are doing most of your scripting
>> with Python instead of R's programming language and then calling
>> the relevant R functions as needed?
>
> No, but if I wanted to do such a thing, I'd look at Sage: http://sagemath.org/

ah .. thanks for the pointer, I had not heard of Sage, I was just
starting to look at
SciPy.

> It'll give you access to a lot more than just R.
>
>> Is there much of a performance hit either way? (as both are interpreted
>> languages)
>
> Are you just asking, or do you have a particular execution time goal, which 
> if exceeded
> would prevent doing this?  I ask because I suspect it's the former, and fast 
> enough is fast
> enough.

I put together a large'ish R program last year, but I think I would be
happier if I could code
it in say Python - but I would rather not do that at the expense of
execution time.

Thanks again for telling me about Sage.

Esmail

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


Re: [R] Python and R

2009-02-17 Thread Barry Rowlingson
2009/2/17 Esmail Bonakdarian :
> Hello all,
>
> I am just wondering if any of you are doing most of your scripting
> with Python instead of R's programming language and then calling
> the relevant R functions as needed?

 I tend to use R in its native form for data analysis and modelling,
and python for all my other programming needs (gui stuff with PyQt4,
web stuff, text processing etc etc).

> And if so, what is your experience with this and what sort of
> software/library  do you use in combination with Python to be able
> to access R's functionality.

 When I need to use the two together, it's easiest with 'rpy'. This
lets you call R functions from python, so you can do:

 from rpy import r
 r.hist(z)

to get a histogram of the values in a python list 'z'. There are some
complications converting structured data types between the two but
they can be overcome, and apparently are handled better with the next
generation Rpy2 (which I've not got into yet). Google for rpy for
info.

> Is there much of a performance hit either way? (as both are interpreted
> languages)

 Not sure what you mean here. Do you mean is:

 R> sum(x)

faster than

Python> sum(x)

and how much worse is:

Python> from rpy import r
Python> r.sum(x)

?

 Knuth's remark on premature optimization applies, as ever

Barry

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


Re: [R] Python and R

2009-02-17 Thread Warren Young

Esmail Bonakdarian wrote:


I am just wondering if any of you are doing most of your scripting
with Python instead of R's programming language and then calling
the relevant R functions as needed?


No, but if I wanted to do such a thing, I'd look at Sage: 
http://sagemath.org/


It'll give you access to a lot more than just R.


Is there much of a performance hit either way? (as both are interpreted
languages)


Are you just asking, or do you have a particular execution time goal, 
which if exceeded would prevent doing this?  I ask because I suspect 
it's the former, and fast enough is fast enough.


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


Re: [R] frequency table for multiple variables

2009-02-17 Thread Hans Ekbrand
On Tue, Feb 17, 2009 at 10:00:40AM -0600, Marc Schwartz wrote:
> on 02/17/2009 09:06 AM Hans Ekbrand wrote:
> > Hi r-help!
> > 
> > Consider the following data-frame:
> > 
> >var1 var2 var3 
> > 1 314 
> > 2 223 
> > 3 223 
> > 4 44   NA 
> > 5 435 
> > 6 223 
> > 7 343 
> > 
> > How can I get R to convert this into the following?
> > 
> > Value 1  2  3  4  5 
> > var1  0  3  2  2  0
> > var2  1  3  1  2  0 
> > var3  0  0  4  1  1
> 
> 
> > t(sapply(DF, function(x) table(factor(x, levels = 1:5
>  1 2 3 4 5
> var1 0 3 2 2 0
> var2 1 3 1 2 0
> var3 0 0 4 1 1
> 
> 
> The key is to turn each column into a factor with explicitly defined
> common levels for tabulation. This enables the table result to have a
> consistent format across each column, allowing for a matrix to be
> created, rather than a list.

Thanks alot, Marc. Neat and efficient, just what I wanted.

BTW, before I saw that you actually included code, I tried on my own,
and wrote this:

my.count <- function(data.frame, levels) {
  result.df <- data.frame(matrix(nrow=length(data.frame),ncol=levels))
  for (i in 1:length(data.frame)) {
result.df[i,] <- table(factor(data.frame[[i]], levels = c(1:levels)))
  }
  result.df
}

which produces the same result. I take this to be a an instructive
example of unnecessary use of for-loops in R.

-- 
Hans Ekbrand (http://sociologi.cjb.net) 
Q. What is that strange attachment in this mail?
A. My digital signature, see www.gnupg.org for info on how you could
   use it to ensure that this mail is from me and has not been
   altered on the way to you.


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


Re: [R] Percentiles/Quantiles with Weighting

2009-02-17 Thread Stavros Macrakis
Here is one kind of weighted quantile function.

The basic idea is very simple:

wquantile <- function( v, w, p )
  {
v <- v[order(v)]
w <- w[order(v)]
v [ which.max( cumsum(w) / sum(w) >= p ) ]
  }

With some more error-checking and general clean-up, it looks like this:

# Simple weighted quantile
#
# v  A numeric vector of observations
# w  A numeric vector of positive weights
# p  The probability 0<=p<=1
#
# Nothing fancy: no interpolation etc.

# Basic idea

wquantile <- function( v, w, p )
  {
v <- v[order(v)]
w <- w[order(v)]
v [ which.max( cumsum(w) / sum(w) >= p ) ]
  }

# Simple weighted quantile
#
# v  A numeric vector of observations
# w  A numeric vector of positive weights
# p  The probability 0<=p<=1
#
# Nothing fancy: no interpolation etc.

wquantile <- function(v,w=rep(1,length(v)),p=.5)
   {
 if (!is.numeric(v) || !is.numeric(w) || length(v) != length(w))
   stop("Values and weights must be equal-length numeric vectors")
 if ( !is.numeric(p) || any( p<0 | p>1 ) )
   stop("Quantiles must be 0<=p<=1")
 ranking <- order(v)
 sumw <- cumsum(w[ranking])
 if ( is.na(w[1]) || w[1]<0 ) stop("Weights must be non-negative numbers")
 plist <- sumw/sumw[length(sumw)]
 sapply(p, function(p) v [ ranking [ which.max( plist >= p ) ] ])
   }

I would appreciate any comments people have on this -- whether
correctness, efficiency, style, 

  -s


On Tue, Feb 17, 2009 at 11:57 AM, Brigid Mooney  wrote:
> Hi All,
>
> I am looking at applications of percentiles to time sequenced data.  I had
> just been using the quantile function to get percentiles over various
> periods, but am more interested in if there is an accepted (and/or
> R-implemented) method to apply weighting to the data so as to weigh recent
> data more heavily.
>
> I wrote the following function, but it seems quite inefficient, and not
> really very flexible in its applications - so if anyone has any suggestions
> on how to look at quantiles/percentiles within R while also using a
> weighting schema, I would be very interested.
>
> Note - this function supposes the data in X is time-sequenced, with the most
> recent (and thus heaviest weighted) data at the end of the vector
>
> WtPercentile <- function(X=rnorm(100), pctile=seq(.1,1,.1))
> {
>  Xprime <- NA
>
>  for(i in 1:length(X))
>  {
>Xprime <- c(Xprime, rep(X[i], times=i))
>  }
>
>  print("Percentiles:")
>  print(quantile(X, pctile))
>  print("Weighted:")
>  print(Xprime)
>  print("Weighted Percentiles:")
>  print(quantile(Xprime, pctile, na.rm=TRUE))
> }
>
> WtPercentile(1:10)
> WtPercentile(rnorm(10))
>
>[[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] Python and R

2009-02-17 Thread Esmail Bonakdarian

Hello all,

I am just wondering if any of you are doing most of your scripting
with Python instead of R's programming language and then calling
the relevant R functions as needed?

And if so, what is your experience with this and what sort of
software/library  do you use in combination with Python to be able
to access R's functionality.

Is there much of a performance hit either way? (as both are interpreted
languages)

Thanks,
Esmail

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


Re: [R] matrix output

2009-02-17 Thread Duncan Murdoch

On 17/02/2009 5:31 PM, phoebe kong wrote:

Hi friends,

I have questions about printing a pretty big size matrix.

As you could see from below, the matrix wasn't showed in R at full size
(11X11), but it was cut partly into three smaller matrices (11X4,11X4,11X3).
I'm wondering if there is a way to show the whole matrix with dimension
11X11, do you know how to make it?


If you increase the output width, e.g.

options(width=1)

won't wrap (but it'll probably be too wide to be useful).  You can also 
reduced the number of decimal places, e.g.


SYrounded <- round(SY, 3)




If R really couldn't fit the full big matrix at once, what about output the
FULL matrix to a .pdf document? I have been wondering for a long time if we
could output something other than graphic, like data frame or text, to a
.pdf file.


The main way to produce nice text in a PDF document through R is to use 
Sweave and LaTeX.  Too much explanation needed for a simple example here.


Duncan Murdoch



Thanks in advance for your patience and answer.

SY

Sldur Slonset Waso  Sboutn
Sldur1  -0.6744252  -0.08427312  -0.2871798
Slonset327   1.000   0.14257353   0.1981339
Waso   327 327.000   1.   0.5104723
Sboutn 327 327.000 327.   1.000
MidSleep   327 327.000 327. 327.000
SleepEff   327 327.000 327. 327.000
NumTrials   22  22.000  22.  22.000
MeanTrials  22  22.000  22.  22.000
MedianTrials22  22.000  22.  22.000
NS3all5tage296 296.000 296. 296.000
HGsoc118   325 325.000 325. 325.000

 MidSleepSleepEff   NumTrials
MeanTrials
Sldur -0.14512244   0.6721889  0.26482213
0.26256352
Slonset0.73991223  -0.5613362 -0.09429701
-0.02540937
Waso   0.08729977  -0.6836098  0.18577075
0.26369283
Sboutn 0.06169839  -0.4895902  0.10897798
0.07058159
MidSleep   1.  -0.1498673 -0.17786561
0.13043478
SleepEff 327.   1.000 -0.01637493
-0.20835686
NumTrials 22.  22.000  1.
-0.16768424
MeanTrials22.  22.000 46.
1.
MedianTrials  22.  22.000 46.
46.
NS3all5tage  296. 296.000 44.
44.
HGsoc118 325. 325.000 44.
44.

 MedianTrials  NS3all5tage HGsoc118
Sldur   0.1857708   0.07849234  -0.03751973
Slonset 0.2919255  -0.14858206   0.05323562
Waso0.3856578  -0.08148054   0.09341454
Sboutn  0.2557877  -0.05049218   0.06847118
MidSleep0.4172784  -0.09344423   0.05287818
SleepEff   -0.3088650   0.12099185  -0.08453191
NumTrials  -0.1108233  -0.06779422  -0.13164200
MeanTrials  0.9203207  -0.07625088   0.10584919
MedianTrials1.000  -0.10894996   0.13615222
NS3all5tage44.000   1.  -0.10554137
HGsoc118   44.000 632.   1.

[[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] matrix output

2009-02-17 Thread phoebe kong
Hi friends,

I have questions about printing a pretty big size matrix.

As you could see from below, the matrix wasn't showed in R at full size
(11X11), but it was cut partly into three smaller matrices (11X4,11X4,11X3).
I'm wondering if there is a way to show the whole matrix with dimension
11X11, do you know how to make it?

If R really couldn't fit the full big matrix at once, what about output the
FULL matrix to a .pdf document? I have been wondering for a long time if we
could output something other than graphic, like data frame or text, to a
.pdf file.

Thanks in advance for your patience and answer.

SY

Sldur Slonset Waso  Sboutn
Sldur1  -0.6744252  -0.08427312  -0.2871798
Slonset327   1.000   0.14257353   0.1981339
Waso   327 327.000   1.   0.5104723
Sboutn 327 327.000 327.   1.000
MidSleep   327 327.000 327. 327.000
SleepEff   327 327.000 327. 327.000
NumTrials   22  22.000  22.  22.000
MeanTrials  22  22.000  22.  22.000
MedianTrials22  22.000  22.  22.000
NS3all5tage296 296.000 296. 296.000
HGsoc118   325 325.000 325. 325.000

 MidSleepSleepEff   NumTrials
MeanTrials
Sldur -0.14512244   0.6721889  0.26482213
0.26256352
Slonset0.73991223  -0.5613362 -0.09429701
-0.02540937
Waso   0.08729977  -0.6836098  0.18577075
0.26369283
Sboutn 0.06169839  -0.4895902  0.10897798
0.07058159
MidSleep   1.  -0.1498673 -0.17786561
0.13043478
SleepEff 327.   1.000 -0.01637493
-0.20835686
NumTrials 22.  22.000  1.
-0.16768424
MeanTrials22.  22.000 46.
1.
MedianTrials  22.  22.000 46.
46.
NS3all5tage  296. 296.000 44.
44.
HGsoc118 325. 325.000 44.
44.

 MedianTrials  NS3all5tage HGsoc118
Sldur   0.1857708   0.07849234  -0.03751973
Slonset 0.2919255  -0.14858206   0.05323562
Waso0.3856578  -0.08148054   0.09341454
Sboutn  0.2557877  -0.05049218   0.06847118
MidSleep0.4172784  -0.09344423   0.05287818
SleepEff   -0.3088650   0.12099185  -0.08453191
NumTrials  -0.1108233  -0.06779422  -0.13164200
MeanTrials  0.9203207  -0.07625088   0.10584919
MedianTrials1.000  -0.10894996   0.13615222
NS3all5tage44.000   1.  -0.10554137
HGsoc118   44.000 632.   1.

[[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] using sapply to apply function to some columns of a dataframe

2009-02-17 Thread Duncan Murdoch

On 17/02/2009 4:42 PM, mwestp...@worldbank.org wrote:

Hello:

I would like to sum every x columns of a dataframe for each row.  For instance,
if x is 10, then for dataframe df, this function will sum the first ten elements
together and then the next ten:

sapply(list(colnames(df)[1:10], colnames(df)[11:20]),function(x)apply( df[,x],
1, sum))

If the number of columns is quite large (1000's), then manually entering the
list above is not practical.  Any suggestions?

I would also like to do a variant of the above, where I sum every nth element.


I think the easiest way to do this is to convert the dataframe into an 
array with 3 indices, and sum over one of them.  For example:


rows <- 20
cols <- 120

df <- matrix(1:(rows*cols), rows, cols)
# in your case, df <- as.matrix( df )

arr <- array( df, c(rows, 10, cols/10))
sums <- apply( arr, c(1,3), sum)

Duncan Murdoch

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


[R] cumsum vs. sum

2009-02-17 Thread Stavros Macrakis
I recently traced a bug of mine to the fact that cumsum(s)[length(s)]
is not always exactly equal to sum(s).

For example,

 x<-1/(12:14)
 sum(x) - cumsum(x)[3]  => 2.8e-17

Floating-point addition is of course not exact, and in particular is
not associative, so there are various possible reasons for this.
Perhaps sum uses clever summing tricks to get more accurate results?
In some quick experiments, it does seem to get more accurate results
than cumsum.

It might be worth documenting.

 -s

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


Re: [R] Subset Regression Package

2009-02-17 Thread Hans W. Borchers
Alex Roy  gmail.com> writes:

> 
> Dear all ,
>   Is there any subset regression (subset selection
> regression) package in R other than "leaps"?


Lars and Lasso are other 'subset selection' methods, see the corresponding
packages 'lars' and 'lasso2' and its description in The Elements of Statistical
Learning.
Also, 'dr', "Methods for dimension reduction for regression", or  'relaimpo',
"Relative importance of regressors in linear models", can be considered.


> Thanks and regards
> 
> Alex
>

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


[R] using sapply to apply function to some columns of a dataframe

2009-02-17 Thread mwestphal

Hello:

I would like to sum every x columns of a dataframe for each row.  For instance,
if x is 10, then for dataframe df, this function will sum the first ten elements
together and then the next ten:

sapply(list(colnames(df)[1:10], colnames(df)[11:20]),function(x)apply( df[,x],
1, sum))

If the number of columns is quite large (1000's), then manually entering the
list above is not practical.  Any suggestions?

I would also like to do a variant of the above, where I sum every nth element.


Thanks,

Michael

--
Michael I. Westphal, PhD
Africa Region Water Resources (AFTWR)
South Asia Region Sustainable Development (SASSD)
World Development Report 2010: "Development in a Changing Climate"
www.worldbank.org/wdr2010
Room J6-007 (mail stop: J6-603)
Tel: 202.473.1217
The World Bank
1818  H St NW, Washington DC 20433, USA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Cross classified or Multiple membership or Hierarchical (3 level ) logistic models using Umacs

2009-02-17 Thread Luwis Tapiwa Diya
Dear R users,

I would like to fit cross classified or multiple membership logistic models
or a 3 level hierarchical logistic model using the Umacs package. Can anyone
advise me on how to proceed or better point me to examples of  how its done.

Regards,

-- 
Luwis Diya,
Leuven Biostatistics and Statistical Bioinformatics Centre (L-BioStat),
Kapucijnenvoer 35 blok d - bus 7001,
3000 Leuven,
Belgium

Tel: +32 16 336886 or +32 16 336892
Fax: +32 16 337015

[[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] Whitening Time Series

2009-02-17 Thread Pele

Hi Bob - your suggesting worked out great... Many thanks!  

Also, thanks everyone for the other suggestions!


Bob McCall wrote:
> 
> Look in the package "forecast" for the function "Arima". It will do what
> you want. It's different than arima function in the stats package.
> Bob
> 
> Pele wrote:
>> 
>> Hi R users,
>> 
>> I am doing cross correlation analysis on  2 time series (call them
>> y-series and x-series) where I need the use the model developed on the
>> x-series to prewhiten the yseries..  Can someone point me to a
>> function/filter in R that would allow me to do that? 
>> 
>> Thanks in advance for any help!
>> 
> 
> 

-- 
View this message in context: 
http://www.nabble.com/Whitening-Time-Series-tp22041765p22066246.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] Subset Regression Package

2009-02-17 Thread Ben Bolker



Alex Roy wrote:
> 
> Dear all ,
>   Is there any subset regression (subset selection
> regression) package in R other than "leaps"?
> 
> 

RSiteSearch("{subset regression}") doesn't turn up much other than
special-purpose tools for ARMA models etc..  What does leaps not do that you
need it to do?

  Ben Bolker

-- 
View this message in context: 
http://www.nabble.com/Subset-Regression-Package-tp22054068p22066217.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] printing out the summary for lm into a txt file

2009-02-17 Thread kayj


this is the error message that I was getting 

 "In sink() ... : no sink to remove"


i got it to work . thanks for the help


David Winsemius wrote:
> 
> " did not work."
> 
> Might that mean errors? Care to share?
> 
> Running that through my wetware R interpreter, I think I am seeing you  
> ask for creation of models with variable outcomes from a the first 100  
> columns of "data". And then you are specifying a formula that includes  
> a weird mixture of 3 vectors, data[,"x1"],x2 x3,... plus all of the  
> data.frame, "data"? Have you tested this loop without the sink's to  
> see if it produces anything meaningful?
> 
> Once you test the process, then issue the sink(filename, append=TRUE)  
> once, before the loop, do your analyses in the loop, and then close  
> with sink() after the loop.
> 
> -- 
> David Winsemius
> On Feb 17, 2009, at 12:52 PM, kayj wrote:
> 
>>
>> Hi All,
>>
>>
>> I am trying to run several linear regressions and print out the  
>> summay and
>> the anova reslts on the top of
>> each other for each model. Below is a sample progarm that did not  
>> work. is
>> it possible to print the
>> anova below the summary of lm in one file?
>>
>> thanks for your help
>>
>>
>>
>> ##
>>
>> data<-read.table("data.txt", header=T, sep='\t')
>>
>> for (i in 1:100){
>>
>> y<-data[,i]
>>
>> lm.model<-lm(y~data$x1+data$x2+data$x3+data)
>>
>> sink("results.txt", append=T)
>>
>> s<-summary(lm.Model)
>> print(s)
>> sink()
>>
>> an<-anova(lm.Model)
>> print(an)
>> sink()
>>
>> }
>> -- 
>> View this message in context:
>> http://www.nabble.com/printing-out-the-summary-for-lm-into-a-txt-file-tp22062643p22062643.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/printing-out-the-summary-for-lm-into-a-txt-file-tp22062643p22065365.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] printing out the summary for lm into a txt file

2009-02-17 Thread kayj

Thanks a lot for your help, it worked.




Greg Snow-2 wrote:
> 
> The sink() command stops the sinking, so you send the lm output to the
> file, then stop the sinking before printing out the anova result.  So the
> simplest thing to try is to put the first sink (with the filename and
> append=T) before you start the loop, remove all calls to sink within the
> loop, then do a single sink() after the loop ends.
> 
> Hope this helps,
> 
> -- 
> Gregory (Greg) L. Snow Ph.D.
> Statistical Data Center
> Intermountain Healthcare
> greg.s...@imail.org
> 801.408.8111
> 
> 
>> -Original Message-
>> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
>> project.org] On Behalf Of kayj
>> Sent: Tuesday, February 17, 2009 10:52 AM
>> To: r-help@r-project.org
>> Subject: [R] printing out the summary for lm into a txt file
>> 
>> 
>> Hi All,
>> 
>> 
>> I am trying to run several linear regressions and print out the summay
>> and
>> the anova reslts on the top of
>> each other for each model. Below is a sample progarm that did not work.
>> is
>> it possible to print the
>> anova below the summary of lm in one file?
>> 
>> thanks for your help
>> 
>> 
>> 
>> ##
>> 
>> data<-read.table("data.txt", header=T, sep='\t')
>> 
>> for (i in 1:100){
>> 
>> y<-data[,i]
>> 
>> lm.model<-lm(y~data$x1+data$x2+data$x3+data)
>> 
>> sink("results.txt", append=T)
>> 
>> s<-summary(lm.Model)
>> print(s)
>> sink()
>> 
>> an<-anova(lm.Model)
>> print(an)
>> sink()
>> 
>> }
>> --
>> View this message in context: http://www.nabble.com/printing-out-the-
>> summary-for-lm-into-a-txt-file-tp22062643p22062643.html
>> Sent from the R help mailing list archive at Nabble.com.
>> 
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-
>> guide.html
>> and provide commented, minimal, self-contained, reproducible code.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: 
http://www.nabble.com/printing-out-the-summary-for-lm-into-a-txt-file-tp22062643p22065292.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] Percentiles/Quantiles with Weighting

2009-02-17 Thread roger koenker


url:www.econ.uiuc.edu/~rogerRoger Koenker
emailrkoen...@uiuc.eduDepartment of Economics
vox: 217-333-4558University of Illinois
fax:   217-244-6678Champaign, IL 61820



On Feb 17, 2009, at 1:58 PM, Brigid Mooney wrote:

Thanks for pointing me to the quantreg package as a resource.  I was  
hoping to ask be able to address one quick follow-up question...


I get slightly different variants between using the rq funciton with  
formula = mydata ~ 1 as I would if I ran the same data using the  
quantile function.

Example:

mydata <- (1:10)^2/2
pctile <- seq(.59, .99, .1)

quantile(mydata, pctile)
59%69%79%89%99%
20.015 26.075 32.935 40.595 49.145

rq(mydata~1, tau=pctile)
Call:
rq(formula = mydata ~ 1, tau = pctile)
Coefficients:
tau= 0.59 tau= 0.69 tau= 0.79 tau= 0.89 tau= 0.99
(Intercept)18  24.532  40.550
Degrees of freedom: 10 total; 9 residual

Is it correct to assume this is due to the different accepted  
methods of calculating quantiles?  If so, do you know where I would  
be able to see the algorithms used in these functions?  I'm not  
finding it in the documentation for function rq, and am new enough  
to R that I don't know where those references would generally be.




Yes,  quantile() in base R documents 9 varieties of quantiles, 2 more  
than

William Empson's famous 7 Types of Ambiguity.  In quantreg the function
rq() finds a solution to an underlying optimization problem and  
doesn't ask

any further into the nature of the ambiguity -- it does often produce a
warning indicating that there may be more than one solution.  The  
default
base R quantile is interpolated, while the default rq() with method =  
"br"
using the simplex algorithm finds an order statistic, typically.  If  
you prefer

something more like interpolation, you can try rq() with method = "fn"
which is using an interior point algorithm and when there are multiple
solutions it tends to produce something more like the centroid of the
solution set.  I hope that this helps.




On Tue, Feb 17, 2009 at 12:29 PM, roger koenker   
wrote:

http://www.nabble.com/weighted-quantiles-to19864562.html#a19865869

gives one possibility...

url:www.econ.uiuc.edu/~rogerRoger Koenker
emailrkoen...@uiuc.eduDepartment of Economics
vox: 217-333-4558University of Illinois
fax:   217-244-6678Champaign, IL 61820




On Feb 17, 2009, at 10:57 AM, Brigid Mooney wrote:

Hi All,

I am looking at applications of percentiles to time sequenced data.   
I had

just been using the quantile function to get percentiles over various
periods, but am more interested in if there is an accepted (and/or
R-implemented) method to apply weighting to the data so as to weigh  
recent

data more heavily.

I wrote the following function, but it seems quite inefficient, and  
not
really very flexible in its applications - so if anyone has any  
suggestions

on how to look at quantiles/percentiles within R while also using a
weighting schema, I would be very interested.

Note - this function supposes the data in X is time-sequenced, with  
the most

recent (and thus heaviest weighted) data at the end of the vector

WtPercentile <- function(X=rnorm(100), pctile=seq(.1,1,.1))
{
 Xprime <- NA

 for(i in 1:length(X))
 {
  Xprime <- c(Xprime, rep(X[i], times=i))
 }

 print("Percentiles:")
 print(quantile(X, pctile))
 print("Weighted:")
 print(Xprime)
 print("Weighted Percentiles:")
 print(quantile(Xprime, pctile, na.rm=TRUE))
}

WtPercentile(1:10)
WtPercentile(rnorm(10))

   [[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] Percentiles/Quantiles with Weighting

2009-02-17 Thread Brigid Mooney
Thanks for pointing me to the quantreg package as a resource.  I was hoping
to ask be able to address one quick follow-up question...

I get slightly different variants between using the rq funciton with formula
= mydata ~ 1 as I would if I ran the same data using the quantile function.

Example:

mydata <- (1:10)^2/2
pctile <- seq(.59, .99, .1)

quantile(mydata, pctile)
59%69%79%89%99%
20.015 26.075 32.935 40.595 49.145

rq(mydata~1, tau=pctile)
Call:
rq(formula = mydata ~ 1, tau = pctile)
Coefficients:
tau= 0.59 tau= 0.69 tau= 0.79 tau= 0.89 tau= 0.99
(Intercept)18  24.532  40.550
Degrees of freedom: 10 total; 9 residual

Is it correct to assume this is due to the different accepted methods of
calculating quantiles?  If so, do you know where I would be able to see the
algorithms used in these functions?  I'm not finding it in the documentation
for function rq, and am new enough to R that I don't know where those
references would generally be.




On Tue, Feb 17, 2009 at 12:29 PM, roger koenker  wrote:

> http://www.nabble.com/weighted-quantiles-to19864562.html#a19865869
>
> gives one possibility...
>
> url:www.econ.uiuc.edu/~rogerRoger Koenker
> emailrkoen...@uiuc.eduDepartment of Economics
> vox: 217-333-4558University of Illinois
> fax:   217-244-6678Champaign, IL 61820
>
>
>
>
> On Feb 17, 2009, at 10:57 AM, Brigid Mooney wrote:
>
>   Hi All,
>>
>> I am looking at applications of percentiles to time sequenced data.  I had
>> just been using the quantile function to get percentiles over various
>> periods, but am more interested in if there is an accepted (and/or
>> R-implemented) method to apply weighting to the data so as to weigh recent
>> data more heavily.
>>
>> I wrote the following function, but it seems quite inefficient, and not
>> really very flexible in its applications - so if anyone has any
>> suggestions
>> on how to look at quantiles/percentiles within R while also using a
>> weighting schema, I would be very interested.
>>
>> Note - this function supposes the data in X is time-sequenced, with the
>> most
>> recent (and thus heaviest weighted) data at the end of the vector
>>
>> WtPercentile <- function(X=rnorm(100), pctile=seq(.1,1,.1))
>> {
>>  Xprime <- NA
>>
>>  for(i in 1:length(X))
>>  {
>>   Xprime <- c(Xprime, rep(X[i], times=i))
>>  }
>>
>>  print("Percentiles:")
>>  print(quantile(X, pctile))
>>  print("Weighted:")
>>  print(Xprime)
>>  print("Weighted Percentiles:")
>>  print(quantile(Xprime, pctile, na.rm=TRUE))
>> }
>>
>> WtPercentile(1:10)
>> WtPercentile(rnorm(10))
>>
>>[[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] printing out the summary for lm into a txt file

2009-02-17 Thread Greg Snow
The sink() command stops the sinking, so you send the lm output to the file, 
then stop the sinking before printing out the anova result.  So the simplest 
thing to try is to put the first sink (with the filename and append=T) before 
you start the loop, remove all calls to sink within the loop, then do a single 
sink() after the loop ends.

Hope this helps,

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


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of kayj
> Sent: Tuesday, February 17, 2009 10:52 AM
> To: r-help@r-project.org
> Subject: [R] printing out the summary for lm into a txt file
> 
> 
> Hi All,
> 
> 
> I am trying to run several linear regressions and print out the summay
> and
> the anova reslts on the top of
> each other for each model. Below is a sample progarm that did not work.
> is
> it possible to print the
> anova below the summary of lm in one file?
> 
> thanks for your help
> 
> 
> 
> ##
> 
> data<-read.table("data.txt", header=T, sep='\t')
> 
> for (i in 1:100){
> 
> y<-data[,i]
> 
> lm.model<-lm(y~data$x1+data$x2+data$x3+data)
> 
> sink("results.txt", append=T)
> 
> s<-summary(lm.Model)
> print(s)
> sink()
> 
> an<-anova(lm.Model)
> print(an)
> sink()
> 
> }
> --
> View this message in context: http://www.nabble.com/printing-out-the-
> summary-for-lm-into-a-txt-file-tp22062643p22062643.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] creating a map

2009-02-17 Thread Greg Snow
It should be in the datasets package that is automatically loaded with R (at 
least my copy), try ?state and you should see the help for it and others.

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

From: Alina Sheyman [mailto:alina...@gmail.com]
Sent: Tuesday, February 17, 2009 11:02 AM
To: Greg Snow
Cc: r-help@r-project.org
Subject: Re: [R] creating a map

Thanks Greg,
  do you know where i can find the sate.center dataset that you mention?
On Tue, Feb 17, 2009 at 12:28 PM, Greg Snow 
mailto:greg.s...@imail.org>> wrote:
You need to give the symbols function the locations where you want the centers 
of the circles to be.  Some datesets with map information also have centers of 
the states that you can use, for the USA, there is the state.center dataset 
that may work for you, or the maptools package function get.Pcent will compute 
a center for polygons (there are probably other similar functions in other 
packages).

For adding circles to a map of the USA, you may want to look at the state.vbm 
data in the TeachingDemos package (works with maptools package), but you will 
need to computer the centers of the polygons, they don't match state.center or 
others.

Hope this helps,

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


> -Original Message-
> From: r-help-boun...@r-project.org 
> [mailto:r-help-boun...@r-
> project.org] On Behalf Of Alina Sheyman
> Sent: Tuesday, February 17, 2009 9:53 AM
> To: r-help@r-project.org
> Subject: [R] creating a map
>
> I'm trying to create a fairly basic map using R. What i want to get is
> the
> map of the country with circles representing a count of students in
> each
> state.
> What I've done so far is as following -
>   map("state")
>
> symbols(data1$count,circles=log(data1$count)*3,fg=col,bg=col,add=T,inch
> es=F)
>
> this gives me the map of the country, but one that's not populated by
> my
> counts.
> Does anyone know what I'm doing wrong?
>
> Also, if anyone can recommend a good reference for creating maps in R,
> I'd
> really appreciate that.
>
> thank you
>
>   [[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] Overdispersion with binomial distribution

2009-02-17 Thread Ben Bolker
 Thanks for the clarification.
  I actually had MASS open to that page while
I was composing my reply but forgot to mention
it (trying to do too many things at once) ...

  Ben Bolker

Prof Brian Ripley wrote:
> On Tue, 17 Feb 2009, Ben Bolker wrote:
> 
>> Jessica L Hite/hitejl/O/VCU  vcu.edu> writes:
>>
>>> I am attempting to run a glm with a binomial model to analyze proportion
>>> data.
>>> I have been following Crawley's book closely and am wondering if there is
>>> an accepted standard for how much is too much overdispersion? (e.g. change
>>> in AIC has an accepted standard of 2).
> 
>>  In principle, in the null case (i.e. data are really binomial)
>> the deviance is  chi-squared distributed with the df equal
>> to the residual df.
> 
> *Approximately*, provided the expected counts are not near or below 
> one.  See MASS §7.5 for an analysis of the size of the approximation 
> errors (which can be large and in both directions).
> 
> Given that I once had a consulting job where the over-dispersion was 
> causing something close ot panic and was entirely illusory, the lack 
> of the 'approximately' can have quite serious consequences.
> 
>>  For example:
>>
>> example(glm)
>> deviance(glm.D93) ## 5.13
>> summary(glm.D93)$dispersion ## 1 (by definition)
>> dfr <- df.residual(glm.D93)
>> deviance(glm.D93)/dfr ## 1.28
>> d2 <- sum(residuals(glm.D93,"pearson")^2) ## 5.17
>> (disp2 <- d2/dfr)  ## 1.293
>>
>> gg2 <- update(glm.D93,family=quasipoisson)
>> summary(gg2)$dispersion  ## 1.293, same as above
>>
>> pchisq(d2,df=dfr,lower.tail=FALSE)
>>
>> all.equal(coef(glm.D93),coef(gg2)) ## TRUE
>>
>> se1 <- coef(summary(glm.D93))[,"Std. Error"]
>> se2 <- coef(summary(gg2))[,"Std. Error"]
>> se2/se1
>>
>> # (Intercept)outcome2outcome3  treatment2  treatment3
>> #   1.1372341.1372341.1372341.1372341.137234
>>
>> sqrt(disp2)
>> # [1] 1.137234
>>
>>> My code and output are below, given the example in the book, these data are
>>> WAY overdispersed .do I mention this and go on or does this signal the
>>> need to try a different model? If so, any suggestions on the type of
>>> distribution (gamma or negative binomial ?)?
>>  Way overdispersed may indicate model lack of fit.  Have
>> you examined residuals/data for outliers etc.?
>>
>>  quasibinomial should be fine, or you can try beta-binomial
>> (see the aod package) ...
>>
>>
>>> attach(Clutch2)
>>>  y<-cbind(Total,Size-Total)
>>> glm1<-glm(y~Pred,"binomial")
>>> summary(glm1)
>>>
>>> Call:
>>> glm(formula = y ~ Pred, family = "binomial")
>>>
>>> Deviance Residuals:
>>> Min   1Q   Median   3Q  Max
>>> -9.1022  -2.7899  -0.4781   2.6058   8.4852
>>>
>>> Coefficients:
>>> Estimate Std. Error z value Pr(>|z|)
>>> (Intercept)  1.350950.06612  20.433  < 2e-16 ***
>>> PredF   -0.348110.11719  -2.970  0.00297 **
>>> PredSN  -3.291560.10691 -30.788  < 2e-16 ***
>>> PredW   -1.464510.12787 -11.453  < 2e-16 ***
>>> PredWF  -0.564120.13178  -4.281 1.86e-05 ***
>>> ---
>>>  the output for residual deviance and df does not change even when I
>>> use quasibinomial, is this ok?  #
>>  That's as expected.
>>
>>>  library(MASS)
>>  you don't really need MASS for quasibinomial.
>>
 glm2<-glm(y~Pred,"quasibinomial")
 summary(glm2)
>>> Call:
>>> glm(formula = y ~ Pred, family = "quasibinomial")
>>>
>>> Deviance Residuals:
>>> Min   1Q   Median   3Q  Max
>>> -9.1022  -2.7899  -0.4781   2.6058   8.4852
>>>
>>> Coefficients:
>>> Estimate Std. Error t value Pr(>|t|)
>>> (Intercept)   1.3510 0.2398   5.633 1.52e-07 ***
>>> PredF-0.3481 0.4251  -0.819  0.41471
>>> PredSN   -3.2916 0.3878  -8.488 1.56e-13 ***
>>> PredW-1.4645 0.4638  -3.157  0.00208 **
>>> PredWF   -0.5641 0.4780  -1.180  0.24063
>>> ---
>>> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
>>>
>>> (Dispersion parameter for quasibinomial family taken to be 13.15786)
>>>
>>> Null deviance: 2815.5  on 108  degrees of freedom
>>> Residual deviance: 1323.5  on 104  degrees of freedom
>>>   (3 observations deleted due to missingness)
>>> AIC: NA
>>>
>>> Number of Fisher Scoring iterations: 5
>>>
 anova(glm2,test="F")
>>> Analysis of Deviance Table
>>>
>>> Model: quasibinomial, link: logit
>>>
>>> Response: y
>>>
>>> Terms added sequentially (first to last)
>>>
>>>   Df Deviance Resid. Df Resid. Dev  F   Pr(>F)
>>> NULL108 2815.5
>>> Pred   4   1492.0   104 1323.5 28.349 6.28e-16 ***
>>> ---
>>> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
 model1<-update(glm2,~.-Pred)
 anova(glm2,model1,test="F")
>>> Analysis of Deviance Table
>>>
>>> Model 1: y ~ Pred
>>> Model 2: y ~ 1
>>>   Resid. Df Resid. Dev  Df Deviance  F   Pr(>F)
>>> 1   104 1323.5
>>> 2   108 2815.5  -4  -1492.0 28.349 6.28e-16 ***
>>> ---
>>> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.0

Re: [R] printing out the summary for lm into a txt file

2009-02-17 Thread David Winsemius

" did not work."

Might that mean errors? Care to share?

Running that through my wetware R interpreter, I think I am seeing you  
ask for creation of models with variable outcomes from a the first 100  
columns of "data". And then you are specifying a formula that includes  
a weird mixture of 3 vectors, data[,"x1"],x2 x3,... plus all of the  
data.frame, "data"? Have you tested this loop without the sink's to  
see if it produces anything meaningful?


Once you test the process, then issue the sink(filename, append=TRUE)  
once, before the loop, do your analyses in the loop, and then close  
with sink() after the loop.


--
David Winsemius
On Feb 17, 2009, at 12:52 PM, kayj wrote:



Hi All,


I am trying to run several linear regressions and print out the  
summay and

the anova reslts on the top of
each other for each model. Below is a sample progarm that did not  
work. is

it possible to print the
anova below the summary of lm in one file?

thanks for your help



##

data<-read.table("data.txt", header=T, sep='\t')

for (i in 1:100){

y<-data[,i]

lm.model<-lm(y~data$x1+data$x2+data$x3+data)

sink("results.txt", append=T)

s<-summary(lm.Model)
print(s)
sink()

an<-anova(lm.Model)
print(an)
sink()

}
--
View this message in context: 
http://www.nabble.com/printing-out-the-summary-for-lm-into-a-txt-file-tp22062643p22062643.html
Sent from the R help mailing list archive at Nabble.com.

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


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


Re: [R] Overdispersion with binomial distribution

2009-02-17 Thread Prof Brian Ripley

On Tue, 17 Feb 2009, Ben Bolker wrote:


Jessica L Hite/hitejl/O/VCU  vcu.edu> writes:


I am attempting to run a glm with a binomial model to analyze proportion
data.
I have been following Crawley's book closely and am wondering if there is
an accepted standard for how much is too much overdispersion? (e.g. change
in AIC has an accepted standard of 2).



 In principle, in the null case (i.e. data are really binomial)
the deviance is  chi-squared distributed with the df equal
to the residual df.


*Approximately*, provided the expected counts are not near or below 
one.  See MASS §7.5 for an analysis of the size of the approximation 
errors (which can be large and in both directions).


Given that I once had a consulting job where the over-dispersion was 
causing something close ot panic and was entirely illusory, the lack 
of the 'approximately' can have quite serious consequences.



 For example:

example(glm)
deviance(glm.D93) ## 5.13
summary(glm.D93)$dispersion ## 1 (by definition)
dfr <- df.residual(glm.D93)
deviance(glm.D93)/dfr ## 1.28
d2 <- sum(residuals(glm.D93,"pearson")^2) ## 5.17
(disp2 <- d2/dfr)  ## 1.293

gg2 <- update(glm.D93,family=quasipoisson)
summary(gg2)$dispersion  ## 1.293, same as above

pchisq(d2,df=dfr,lower.tail=FALSE)

all.equal(coef(glm.D93),coef(gg2)) ## TRUE

se1 <- coef(summary(glm.D93))[,"Std. Error"]
se2 <- coef(summary(gg2))[,"Std. Error"]
se2/se1

# (Intercept)outcome2outcome3  treatment2  treatment3
#   1.1372341.1372341.1372341.1372341.137234

sqrt(disp2)
# [1] 1.137234


My code and output are below, given the example in the book, these data are
WAY overdispersed .do I mention this and go on or does this signal the
need to try a different model? If so, any suggestions on the type of
distribution (gamma or negative binomial ?)?


 Way overdispersed may indicate model lack of fit.  Have
you examined residuals/data for outliers etc.?

 quasibinomial should be fine, or you can try beta-binomial
(see the aod package) ...



attach(Clutch2)
 y<-cbind(Total,Size-Total)
glm1<-glm(y~Pred,"binomial")
summary(glm1)

Call:
glm(formula = y ~ Pred, family = "binomial")

Deviance Residuals:
Min   1Q   Median   3Q  Max
-9.1022  -2.7899  -0.4781   2.6058   8.4852

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.350950.06612  20.433  < 2e-16 ***
PredF   -0.348110.11719  -2.970  0.00297 **
PredSN  -3.291560.10691 -30.788  < 2e-16 ***
PredW   -1.464510.12787 -11.453  < 2e-16 ***
PredWF  -0.564120.13178  -4.281 1.86e-05 ***
---
 the output for residual deviance and df does not change even when I
use quasibinomial, is this ok?  #


 That's as expected.


 library(MASS)


 you don't really need MASS for quasibinomial.


glm2<-glm(y~Pred,"quasibinomial")
summary(glm2)


Call:
glm(formula = y ~ Pred, family = "quasibinomial")

Deviance Residuals:
Min   1Q   Median   3Q  Max
-9.1022  -2.7899  -0.4781   2.6058   8.4852

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.3510 0.2398   5.633 1.52e-07 ***
PredF-0.3481 0.4251  -0.819  0.41471
PredSN   -3.2916 0.3878  -8.488 1.56e-13 ***
PredW-1.4645 0.4638  -3.157  0.00208 **
PredWF   -0.5641 0.4780  -1.180  0.24063
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

(Dispersion parameter for quasibinomial family taken to be 13.15786)

Null deviance: 2815.5  on 108  degrees of freedom
Residual deviance: 1323.5  on 104  degrees of freedom
  (3 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 5


anova(glm2,test="F")

Analysis of Deviance Table

Model: quasibinomial, link: logit

Response: y

Terms added sequentially (first to last)

  Df Deviance Resid. Df Resid. Dev  F   Pr(>F)
NULL108 2815.5
Pred   4   1492.0   104 1323.5 28.349 6.28e-16 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

model1<-update(glm2,~.-Pred)
anova(glm2,model1,test="F")

Analysis of Deviance Table

Model 1: y ~ Pred
Model 2: y ~ 1
  Resid. Df Resid. Dev  Df Deviance  F   Pr(>F)
1   104 1323.5
2   108 2815.5  -4  -1492.0 28.349 6.28e-16 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

coef(glm2)

(Intercept)   PredF  PredSN   PredW  PredWF
  1.3509550  -0.3481096  -3.2915601  -1.4645097  -0.5641223

Thanks
Jessica
hitejl  vcu.edu



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

[R] multiple levels of nesting with trellis plots

2009-02-17 Thread R User R User
Hi guys,
I have a tricky problem that I'd appreciate your help with.

I have two categorical variables, say varA and varB and an associated
frequency Freq for combinations of the levels of varA and varB. This was
created with a table() call.

I'd now like to make panel plots of the frequency. I can do a panel plot by
varA:
barchart(data$Freq ~ data$varB | data$varA)

How do I achieve two levels of nesting for say a third categorical variable
varC?
That is, how do I specify 2 'given's in the the panel plot, and preferably
arranged the panels vertically within the 'inner' panel/

thanks very much,
Richie

[[alternative HTML version deleted]]

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


Re: [R] how to control the overall shape of a figure?

2009-02-17 Thread Oliver
Thanks very much, exactly what I need.

Oliver

On Feb 16, 10:36 pm, Marc Schwartz  wrote:
> on 02/16/2009 07:51 PM Oliver wrote:
>
>
>
> > hi,
>
> > I am a R beginner. One thing I notice is that when do graphing is,
>
> > if I want to draw two figures in a row such as this:
>
> > par(mfrow(1, 2))
> > plot(...)
> > plot(...)
>
> > Each figure inside will be rectangle instead of the familiar square
> > shape.
> > Though you can drag the edge the window to resize it. I would have
> > prefer this can be done automatically ... also when I do the pdf
> > export for example. Is this possible?
>
> > The thing is if you do par (mfrow(2,2)), then you got all 4 figures in
> > perfect square shape ... but one row two figures are pretty common to
> > me.
>
> > thanks for help
>
> > Oliver
>
> The overall shape of the plot region is controlled by par("pty"), which
> is set to "m" by default, for 'maximal' plotting region.
>
> Set this to "s" to create a square plot region. For example:
>
> par(mfrow = c(1, 2), pty = "s")
> plot(1)
> plot(1)
>
> In order to do this with a PDF file, ensure that the overall plot
> dimensions are in the proper relative aspect ratio. For example:
>
> pdf(file = "SquarePlots.pdf", height = 4, width = 8)
> par(mfrow = c(1, 2), pty= "s")
> plot(1)
> plot(1)
> dev.off()
>
> This will create 2 square plots, side-by-side, within an overall plot
> size of 4x8. You can further adjust the plot margins and other
> characteristics as may be required.
>
> See ?par and ?Devices for more information regarding these
> customizations and options to set size arguments for specific devices.
>
> HTH,
>
> Marc Schwartz
>
> __
> r-h...@r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guidehttp://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] printing out the summary for lm into a txt file

2009-02-17 Thread kayj

Hi All,


I am trying to run several linear regressions and print out the summay and
the anova reslts on the top of 
each other for each model. Below is a sample progarm that did not work. is
it possible to print the
anova below the summary of lm in one file?

thanks for your help



##

data<-read.table("data.txt", header=T, sep='\t')

for (i in 1:100){

y<-data[,i]

lm.model<-lm(y~data$x1+data$x2+data$x3+data)

sink("results.txt", append=T)

s<-summary(lm.Model)
print(s)
sink()

an<-anova(lm.Model)
print(an)
sink()

} 
-- 
View this message in context: 
http://www.nabble.com/printing-out-the-summary-for-lm-into-a-txt-file-tp22062643p22062643.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] Error with "make" with R-devel

2009-02-17 Thread Lana Schaffer
Hi,
I am getting an error compiling the R-devel on a suse architecture
64-bit architecture.  The cp attribute is sending 'trusted.lov'
and an error.  This is a sample of the output:

>  make[3]: Entering directory
>`/lustre/people/schaffer/R-devel/src/library/base'
>building package 'base'
>make[4]: Entering directory
>`/lustre/people/schaffer/R-devel/src/library/base'
>all.R is unchanged
>cmp: EOF on ../../../library/base/R/base
>cp: getting attribute `trusted.lov' of `all.R': Operation not permitted
>make[4]: *** [mkR] Error 1
>make[4]: Leaving directory
>`/lustre/people/schaffer/R-devel/src/library/base'
>make[3]: *** [all] Error 2

Is there someone who can modify the Makefile so that I am able
to compile R-devel?

Lana Schaffer
Biostatistics/Informatics
The Scripps Research Institute
DNA Array Core Facility
La Jolla, CA 92037
(858) 784-2263
(858) 784-2994
schaf...@scripps.edu 

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


Re: [R] creating a map

2009-02-17 Thread Alina Sheyman
Thanks Greg,
  do you know where i can find the sate.center dataset that you mention?

On Tue, Feb 17, 2009 at 12:28 PM, Greg Snow  wrote:

> You need to give the symbols function the locations where you want the
> centers of the circles to be.  Some datesets with map information also have
> centers of the states that you can use, for the USA, there is the
> state.center dataset that may work for you, or the maptools package function
> get.Pcent will compute a center for polygons (there are probably other
> similar functions in other packages).
>
> For adding circles to a map of the USA, you may want to look at the
> state.vbm data in the TeachingDemos package (works with maptools package),
> but you will need to computer the centers of the polygons, they don't match
> state.center or others.
>
> Hope this helps,
>
> --
> Gregory (Greg) L. Snow Ph.D.
> Statistical Data Center
> Intermountain Healthcare
> greg.s...@imail.org
> 801.408.8111
>
>
> > -Original Message-
> > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> > project.org] On Behalf Of Alina Sheyman
> > Sent: Tuesday, February 17, 2009 9:53 AM
> > To: r-help@r-project.org
> > Subject: [R] creating a map
> >
> > I'm trying to create a fairly basic map using R. What i want to get is
> > the
> > map of the country with circles representing a count of students in
> > each
> > state.
> > What I've done so far is as following -
> >   map("state")
> >
> > symbols(data1$count,circles=log(data1$count)*3,fg=col,bg=col,add=T,inch
> > es=F)
> >
> > this gives me the map of the country, but one that's not populated by
> > my
> > counts.
> > Does anyone know what I'm doing wrong?
> >
> > Also, if anyone can recommend a good reference for creating maps in R,
> > I'd
> > really appreciate that.
> >
> > thank you
> >
> >   [[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] Overdispersion with binomial distribution

2009-02-17 Thread Ben Bolker
Jessica L Hite/hitejl/O/VCU  vcu.edu> writes:

> I am attempting to run a glm with a binomial model to analyze proportion
> data.
> I have been following Crawley's book closely and am wondering if there is
> an accepted standard for how much is too much overdispersion? (e.g. change
> in AIC has an accepted standard of 2).

  In principle, in the null case (i.e. data are really binomial)
the deviance is  chi-squared distributed with the df equal
to the residual df.

  For example:

example(glm)
deviance(glm.D93) ## 5.13
summary(glm.D93)$dispersion ## 1 (by definition)
dfr <- df.residual(glm.D93)
deviance(glm.D93)/dfr ## 1.28
d2 <- sum(residuals(glm.D93,"pearson")^2) ## 5.17
(disp2 <- d2/dfr)  ## 1.293

gg2 <- update(glm.D93,family=quasipoisson)
summary(gg2)$dispersion  ## 1.293, same as above

pchisq(d2,df=dfr,lower.tail=FALSE)

all.equal(coef(glm.D93),coef(gg2)) ## TRUE

se1 <- coef(summary(glm.D93))[,"Std. Error"]
se2 <- coef(summary(gg2))[,"Std. Error"]
se2/se1

# (Intercept)outcome2outcome3  treatment2  treatment3 
#   1.1372341.1372341.1372341.1372341.137234 

sqrt(disp2)
# [1] 1.137234

> My code and output are below, given the example in the book, these data are
> WAY overdispersed .do I mention this and go on or does this signal the
> need to try a different model? If so, any suggestions on the type of
> distribution (gamma or negative binomial ?)?

  Way overdispersed may indicate model lack of fit.  Have
you examined residuals/data for outliers etc.?  

  quasibinomial should be fine, or you can try beta-binomial
(see the aod package) ...


> attach(Clutch2)
>  y<-cbind(Total,Size-Total)
> glm1<-glm(y~Pred,"binomial")
> summary(glm1)
> 
> Call:
> glm(formula = y ~ Pred, family = "binomial")
> 
> Deviance Residuals:
> Min   1Q   Median   3Q  Max
> -9.1022  -2.7899  -0.4781   2.6058   8.4852
> 
> Coefficients:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept)  1.350950.06612  20.433  < 2e-16 ***
> PredF   -0.348110.11719  -2.970  0.00297 **
> PredSN  -3.291560.10691 -30.788  < 2e-16 ***
> PredW   -1.464510.12787 -11.453  < 2e-16 ***
> PredWF  -0.564120.13178  -4.281 1.86e-05 ***
> ---
>  the output for residual deviance and df does not change even when I
> use quasibinomial, is this ok?  #

  That's as expected.
 
>  library(MASS)

  you don't really need MASS for quasibinomial.

> > glm2<-glm(y~Pred,"quasibinomial")
> > summary(glm2)
> 
> Call:
> glm(formula = y ~ Pred, family = "quasibinomial")
> 
> Deviance Residuals:
> Min   1Q   Median   3Q  Max
> -9.1022  -2.7899  -0.4781   2.6058   8.4852
> 
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept)   1.3510 0.2398   5.633 1.52e-07 ***
> PredF-0.3481 0.4251  -0.819  0.41471
> PredSN   -3.2916 0.3878  -8.488 1.56e-13 ***
> PredW-1.4645 0.4638  -3.157  0.00208 **
> PredWF   -0.5641 0.4780  -1.180  0.24063
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> 
> (Dispersion parameter for quasibinomial family taken to be 13.15786)
> 
> Null deviance: 2815.5  on 108  degrees of freedom
> Residual deviance: 1323.5  on 104  degrees of freedom
>   (3 observations deleted due to missingness)
> AIC: NA
> 
> Number of Fisher Scoring iterations: 5
> 
> > anova(glm2,test="F")
> Analysis of Deviance Table
> 
> Model: quasibinomial, link: logit
> 
> Response: y
> 
> Terms added sequentially (first to last)
> 
>   Df Deviance Resid. Df Resid. Dev  F   Pr(>F)
> NULL108 2815.5
> Pred   4   1492.0   104 1323.5 28.349 6.28e-16 ***
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> > model1<-update(glm2,~.-Pred)
> > anova(glm2,model1,test="F")
> Analysis of Deviance Table
> 
> Model 1: y ~ Pred
> Model 2: y ~ 1
>   Resid. Df Resid. Dev  Df Deviance  F   Pr(>F)
> 1   104 1323.5
> 2   108 2815.5  -4  -1492.0 28.349 6.28e-16 ***
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> > coef(glm2)
> (Intercept)   PredF  PredSN   PredW  PredWF
>   1.3509550  -0.3481096  -3.2915601  -1.4645097  -0.5641223
> 
> Thanks
> Jessica
> hitejl  vcu.edu

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


Re: [R] creating a map

2009-02-17 Thread Greg Snow
You need to give the symbols function the locations where you want the centers 
of the circles to be.  Some datesets with map information also have centers of 
the states that you can use, for the USA, there is the state.center dataset 
that may work for you, or the maptools package function get.Pcent will compute 
a center for polygons (there are probably other similar functions in other 
packages).

For adding circles to a map of the USA, you may want to look at the state.vbm 
data in the TeachingDemos package (works with maptools package), but you will 
need to computer the centers of the polygons, they don't match state.center or 
others.

Hope this helps,

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


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
> project.org] On Behalf Of Alina Sheyman
> Sent: Tuesday, February 17, 2009 9:53 AM
> To: r-help@r-project.org
> Subject: [R] creating a map
> 
> I'm trying to create a fairly basic map using R. What i want to get is
> the
> map of the country with circles representing a count of students in
> each
> state.
> What I've done so far is as following -
>   map("state")
> 
> symbols(data1$count,circles=log(data1$count)*3,fg=col,bg=col,add=T,inch
> es=F)
> 
> this gives me the map of the country, but one that's not populated by
> my
> counts.
> Does anyone know what I'm doing wrong?
> 
> Also, if anyone can recommend a good reference for creating maps in R,
> I'd
> really appreciate that.
> 
> thank you
> 
>   [[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] Percentiles/Quantiles with Weighting

2009-02-17 Thread David Winsemius
I do know that Harrell's Quantile function in the Hmisc package will  
allow quantile estimates from models. Whether it is general enough to  
extend to time series, I have no experience and cannot say.


--
David Winsemius


On Feb 17, 2009, at 11:57 AM, Brigid Mooney wrote:


Hi All,

I am looking at applications of percentiles to time sequenced data.   
I had

just been using the quantile function to get percentiles over various
periods, but am more interested in if there is an accepted (and/or
R-implemented) method to apply weighting to the data so as to weigh  
recent

data more heavily.

I wrote the following function, but it seems quite inefficient, and  
not
really very flexible in its applications - so if anyone has any  
suggestions

on how to look at quantiles/percentiles within R while also using a
weighting schema, I would be very interested.

Note - this function supposes the data in X is time-sequenced, with  
the most

recent (and thus heaviest weighted) data at the end of the vector

WtPercentile <- function(X=rnorm(100), pctile=seq(.1,1,.1))
{
 Xprime <- NA

 for(i in 1:length(X))
 {
   Xprime <- c(Xprime, rep(X[i], times=i))
 }

 print("Percentiles:")
 print(quantile(X, pctile))
 print("Weighted:")
 print(Xprime)
 print("Weighted Percentiles:")
 print(quantile(Xprime, pctile, na.rm=TRUE))
}

WtPercentile(1:10)
WtPercentile(rnorm(10))

[[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] merging files with different structures

2009-02-17 Thread Grant Gillis
> Hello list,
>
> I am sorry for the previous half post.  I accidentily hit send.  Thanks
> again in advance for any help.
>
> I have many (approx 20) files that I have merged.  Each data set contains
> rows for individuals and data in 2 - 5 columns (depending upon which data
> set).  The individuals in each data set are not necessarily the same and are
> all duplicated (different data in columns) across sheets I am trying to
> merge.  I have used the merge function
>
For example
>
> d1<-read.csv("AlleleReport.csv")
> d2<-read.csv("AlleleReport.csv")
>
> m1 <- merge(d1, d2,  by = c("IND", intersect(colnames(d1), colnames(d2))),
> all = TRUE)
> m2 <- merge(m1, d3,  by = c("IND", intersect(colnames(m1), colnames(d3))),
> all = TRUE)



My problem is that when the data is merged it looks something like


Ind   L1 L1.1L2 L2.1L3   L3.1
a 12  13  NA NA  NA  NA
a  NA NA 22 43   34   45
b  14  1545   64   NA  NA
b   NANANA  NA 99  84

Is there a way that I can merge the rows for each individual?

Cheers


Grant







>
>
>
>
>
>
>
>

[[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] spss-file problem with foreign 0.8-32

2009-02-17 Thread Prof Brian Ripley
It is a issue specific to 0.8-32 and some files (most likely those 
with some (not all) Windows codepages declared).


We are trying to collect together some examples, and will update 
foreign accordingly later in the week.


On Tue, 17 Feb 2009, Harry Haupt wrote:



Hi,
after updating to foreign version 0.8-32, I experienced the following error
when I tried to load a SPSS file:

Fehler in inherits(x, "factor") : objekt "cp" nicht gefunden
Zusätzlich: Warning message:
In read.spss("***l.sav", use.value.labels = TRUE, to.data.frame = TRUE) :
 ***.sav: File-indicated character representation code (1252) looks like a
Windows codepage

It worked without problems with earlier versions.

Thanks for any clues,
best,
Harry
--
View this message in context: 
http://www.nabble.com/spss-file-problem-with-foreign-0.8-32-tp22059259p22059259.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.



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


[R] merging files with different structures

2009-02-17 Thread Grant Gillis
Hello list,

Thanks in advance for any help.

I have many (approx 20) files that I have merged.  For example

d1<-read.csv("AlleleReport.csv")
d2<-read.csv("AlleleReport.csv")

m1 <- merge(d1, d2,  by = c("IND", intersect(colnames(d1), colnames(d2))),
all = TRUE)
m2 <- merge(m1, d3,  by = c("IND", intersect(colnames(m1), colnames(d3))),
all = TRUE)

[[alternative HTML version deleted]]

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


Re: [R] creating a map

2009-02-17 Thread David Winsemius

Two places that have worked examples leap to mind:

--- Sarkar's online accompaniment to his book:
http://lmdvr.r-forge.r-project.org/figures/figures.html
Thumbing through the hard copy I see Figure 6.5 might of interest.

--- Addicted to R's graphics gallery:
http://addictedtor.free.fr/graphiques/search.php?q=map&engine=RGG

--  
David Winsemius


On Feb 17, 2009, at 11:53 AM, Alina Sheyman wrote:

I'm trying to create a fairly basic map using R. What i want to get  
is the
map of the country with circles representing a count of students in  
each

state.
What I've done so far is as following -
 map("state")
symbols 
(data1$count,circles=log(data1$count)*3,fg=col,bg=col,add=T,inches=F)


this gives me the map of the country, but one that's not populated  
by my

counts.
Does anyone know what I'm doing wrong?

Also, if anyone can recommend a good reference for creating maps in  
R, I'd

really appreciate that.

thank you

[[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] Percentiles/Quantiles with Weighting

2009-02-17 Thread Brigid Mooney
Hi All,

I am looking at applications of percentiles to time sequenced data.  I had
just been using the quantile function to get percentiles over various
periods, but am more interested in if there is an accepted (and/or
R-implemented) method to apply weighting to the data so as to weigh recent
data more heavily.

I wrote the following function, but it seems quite inefficient, and not
really very flexible in its applications - so if anyone has any suggestions
on how to look at quantiles/percentiles within R while also using a
weighting schema, I would be very interested.

Note - this function supposes the data in X is time-sequenced, with the most
recent (and thus heaviest weighted) data at the end of the vector

WtPercentile <- function(X=rnorm(100), pctile=seq(.1,1,.1))
{
  Xprime <- NA

  for(i in 1:length(X))
  {
Xprime <- c(Xprime, rep(X[i], times=i))
  }

  print("Percentiles:")
  print(quantile(X, pctile))
  print("Weighted:")
  print(Xprime)
  print("Weighted Percentiles:")
  print(quantile(Xprime, pctile, na.rm=TRUE))
}

WtPercentile(1:10)
WtPercentile(rnorm(10))

[[alternative HTML version deleted]]

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


Re: [R] How to connect R and WinBUGS/OpenBUGS/LinBUGS in Linux in Feb. 2009

2009-02-17 Thread Uwe Ligges



Paul Heinrich Dietrich wrote:

Hi all,
I've managed to get JAGS working on my Ubuntu Hardy Linux with a 32-bit
computer and AMD processors using R 2.8.1.  JAGS is great.  I've read that
JAGS is the fastest, but that hasn't been my experience.  At any rate, I
have more experience with WinBUGS under Windows and would like a version of
that working as well.

It seems like I've read a lot on the subject and tried a lot, but haven't
managed to get BUGS to work yet.  The most success I've had is to install
WinBUGS or OpenBUGS using this method:
http://www.math.aau.dk/~slb/kurser/bayes-08/install.html

What you also need to know is that you need to open Wine and add a drive. 
Although Z is recommended, I haven't been able to specify it, but have

gotten a D drive to work, using:

wine D:/opt/OpenBUGS/winbugs.exe

Using this method, OpenBUGS opens.  Now, to be able to open it with R.  I've
read all sorts of discussions about BRugs (which is no longer on CRAN, but
old versions can still be found), rbugs, and R2WinBUGS (which I'm used to
using on Windows with WinBUGS).  Some people say R2WinBUGS cannot run
OpenBUGS on Linux, some claim they've done it (I think).  It seems the same
thing with everything else.  I've tried making the linbugs and cbugs file
recommended elsewhere online.  It's all very confusing.


For short: It is quite unlikely that BRugs / OpenBUGS (which is called 
LinBUGS under Linux) works natively under your Linux (although it might 
work under very specific settings). BRugs is available for Windows users 
from the "CRAN extras" repsository maintained by Brian Ripley. We moved 
it in order to meet GPL compliance issues.


Hence a standard recommendation is to use R2WinBUGS under native R under 
Linux with WinBUGS running under wine. R2WinBUGS can use wine to do so. 
See the help page ?bugs once you have loaded R2WinBUGS.


Best wishes,
Uwe Ligges


Can someone show a method that works currently, along with some sample code? 
I'm also new to Linux, and confused by path conventions.  For example, in

rbugs, it shows an example of a path such as
"/var/scratch/jyan/wine-20040408/wine", and I don't see how to modify this. 
I have no /var/scratch to begin with, and think Wine is installed in

/home/me/.wine...(I don't have Linux in front of me right now).

Please help.  Thanks.


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


Re: [R] ARIMA models

2009-02-17 Thread Gabor Grothendieck
see auto.arima in the forecast package.

On Tue, Feb 17, 2009 at 10:20 AM, emj83  wrote:
>
> is there some sort of R function which can advise me of the best ARIMA(p,q,r)
> model to use based on the Schwarz criterion e.g for e.g p=0-5, q =0, r=0-5
> or for example p+r< 5???
>
> or is this something I will have to write my own code for?
>
> Thanks Emma
> --
> View this message in context: 
> http://www.nabble.com/ARIMA-models-tp22059382p22059382.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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


[R] creating a map

2009-02-17 Thread Alina Sheyman
I'm trying to create a fairly basic map using R. What i want to get is the
map of the country with circles representing a count of students in each
state.
What I've done so far is as following -
  map("state")
 symbols(data1$count,circles=log(data1$count)*3,fg=col,bg=col,add=T,inches=F)

this gives me the map of the country, but one that's not populated by my
counts.
Does anyone know what I'm doing wrong?

Also, if anyone can recommend a good reference for creating maps in R, I'd
really appreciate that.

thank you

[[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] Chromatogram deconvolution and peak matching

2009-02-17 Thread Katharine Mullen
Hoi Bart,

I think you're right that ALS should be applicable to this problem.
Unfortunately in writing I see that there is a bug when the spectra are
NOT constrained to nonnegative values (the package has been used to my
knowledge only in fitting multiway mass spectra thus far, where this
constraint is always applied); I'll have a patched version soon.

I'd be interested in hearing about where the manual could be improved,
too.

#2D chromatogram generation

time <- seq(0,20,by=0.05)
f <- function(x,rt) dnorm((x-rt),mean=0,sd=rt/35)
C1 <- matrix(c( f(time,6.1), f(time,5.6), f(time,15)), nrow=401,ncol=3)
C2 <- matrix(c( f(time,2.1), f(time,4), f(time,8)), nrow=401,ncol=3)

#spectrum generation
spectra <- function(x,a,b,c,d,e) a + b*(x-e) + c*((x-e)^2) + d*((x-e)^3)
x <- 220:300
s1 <- spectra(x,(-194.2),2.386,(-0.009617),(1.275e-05),0)
s2 <- spectra(x,(-1.054e02),1.3,(-5.239e-03),(6.927e-06),-20)
s3 <- spectra(x,(-194.2),2.386,(-0.009617),(1.275e-05),20)
spec <- matrix(c(s1,s2,s3), nrow=81,ncol=3)

## data matrix 1
chrom1 <- C1 %*% t(spec)

## data matrix 2
chrom2 <- C2 %*% t(spec)

require(ALS)

## you want to recover 2 (401 x 3) matrices C1 and C2 and a
## (82 x 3) matrix E such that
## chrom1 = C1 E^T, chrom2 = C2 E^T

## some starting guess for elution profiles
## these need to be pretty good
Cstart1 <- matrix(c( f(time,4), f(time,6), f(time,12)), nrow=401,ncol=3)
Cstart2 <- matrix(c( f(time,3), f(time,5), f(time,10)), nrow=401,ncol=3)

xx <- als( CList = list(Cstart1, Cstart2), PsiList = list(chrom1, chrom2),
  S=matrix(1, nrow=81, ncol=3), x=time, x2=x, uniC=TRUE,
  normS=TRUE, nonnegS=TRUE)

par(mfrow=c(3,1))

matplot(time, xx$CList[[1]], col=2, type="l",
main="elution profiles chrom1", lty=1,
ylab="amplitude")

matplot(time, C1, col=1, add=TRUE, type="l", lty=1)

legend(10,2,legend=c("real", "estimated"), col=1:2, lty=1)

matplot(time, xx$CList[[2]], col=2, type="l",
main="elution profiles chrom2", lty=1,
ylab="amplitude")

matplot(time, C2, col=1, add=TRUE, type="l", lty=1)

legend(10,6,legend=c("real", "estimated"), col=1:2, lty=1)

matplot(x, xx$S, col=2, type="l", main="spectra", lty=1,
ylab="amplitude")

matplot(x, spec, col=1, add=TRUE, type="l", lty=1)

legend(270,.95,legend=c("real", "estimated"), col=1:2, lty=1)


On Tue, 17 Feb 2009, bartjoosen wrote:

>
> Hi,
>
> I'm trying to match peaks between chromatographic runs.
> I'm able to match peaks when they are chromatographed with the same method,
> but not when there are different methods are used and spectra comes in to
> play.
>
> While searching I found the ALS package which should be usefull for my
> application, but I couldn't figure it out.
>
> I made some dummy chroms with R, which mimic my actual datasets, to play
> with, but after looking at the manuals of ALS, I'm affraid I can't get the
> job done. Can someone put me on the right way?
>
> Here is my code to generate the dummy chroms, which also plots the 2 chroms
> and the spectra of the 3 peaks:
>
> #2D chromatogram generation
> par(mfrow=c(3,1))
> time <- seq(0,20,by=0.05)
> f <- function(x,rt) dnorm((x-rt),mean=0,sd=rt/35)
> c1 <- f(time,6.1)
> c2 <- f(time,5.6)
> c3 <- f(time,15)
> plot(c1+c2+c3~time,type="l",main="chrom1")
>
> #spectrum generation
> spectra <- function(x,a,b,c,d,e) a + b*(x-e) + c*((x-e)^2) + d*((x-e)^3)
> x <- 220:300
> s1 <- spectra(x,(-194.2),2.386,(-0.009617),(1.275e-05),0)
> s2 <- spectra(x,(-1.054e02),1.3,(-5.239e-03),(6.927e-06),-20)
> s3 <- spectra(x,(-194.2),2.386,(-0.009617),(1.275e-05),20)
>
> chrom1.tot <-
> data.frame(time,outer(c1,s1,"*")+outer(c2,s2,"*")+outer(c2,s2,"*"))
> names(chrom.tot)[-1] <- x
>
> #generation of chromatogram 2
> c1 <- f(time,2.1)
> c2 <- f(time,4)
> c3 <- f(time,8)
> plot(c1+c2+c3~time,type="l",main="chrom2")
>
> chrom2.tot <-
> data.frame(time,outer(c1,s1,"*")+outer(c2,s2,"*")+outer(c2,s2,"*"))
> names(chrom.tot)[-1] <- x
>
> plot(s1~x,type="l",main="spectra")
> lines(s2~x,col=2)
> lines(s3~x,col=3)
>
> Thanks for your time
>
> Kind Regards
>
> Bart
> --
> View this message in context: 
> http://www.nabble.com/Chromatogram-deconvolution-and-peak-matching-tp22057592p22057592.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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


Re: [R] Comparison of age categories using contrasts

2009-02-17 Thread Mark Difford

Hi Dylan, Chuck,

Mark Difford wrote:
>> Coming to your question [?] about how to generate the kind of contrasts
>> that Patrick wanted 
>> using contrast.Design. Well, it is not that straightforward, though I may
>> have missed 
>> something in the documentation to the function. In the past I have
>> cobbled them together 
>> using the kind of hack given below:

Here is a much simpler route, and a correction to my earlier posting (helped
by a little off-list prompting from Prof. Harrell). All that is required to
get successive difference (i.e. sdif) contrasts using contrast.Design is the
following:


contrast(l, a=list(f=c("a","b","c")), b=list(f=c("b","c","d")))

## Run a full example
set.seed(10101010) 
x <- rnorm(100) 
y1 <- x[1:25] * 2 + rnorm(25, mean=1) 
y2 <- x[26:50] * 2.6 + rnorm(25, mean=1.5) 
y3 <- x[51:75] * 2.9 + rnorm(25, mean=5) 
y4 <- x[76:100] * 3.5 + rnorm(25, mean=5.5) 

d <- data.frame(x=x, y=c(y1,y2,y3,y4), f=factor(rep(letters[1:4], each=25)))

dd <- datadist(d) 
options(datadist='dd') 
l <- ols(y ~ x * f, data=d)

## compare with multcomp by setting x=0
summary(glht(l, linfct=mcp(f="Seq")))
contrast(l, b=list(x=0, f=c("a","b","c")), a=list(x=0, f=c("b","c","d")))

Regards, and apologies for any confusion, Mark.


Mark Difford wrote:
> 
> Hi Dylan,
> 
>>> Am I trying to use contrast.Design() for something that it was not
>>> intended for? ...
> 
> I think Prof. Harrell's main point had to do with how interactions are
> handled. You can also get the kind of contrasts that Patrick was
> interested in via multcomp. If we do this using your artificial data set
> we see that the contrasts are the same as those got by fitting the model
> using contr.sdif, but a warning is generated about interactions in the
> model &c. [see example code]
> 
> Part of Prof. Harrell's "system" is that in generating contrasts via
> contr.Design an appropriate value is automatically chosen for the
> interacting variable (in this case the median value of x) so that a
> reasonable default set of contrasts is calculated. This can of course be
> changed.
> 
> Coming to your question [?] about how to generate the kind of contrasts
> that Patrick wanted using contrast.Design. Well, it is not that
> straightforward, though I may have missed something in the documentation
> to the function. In the past I have cobbled them together using the kind
> of hack given below:
> 
> ## Exampe code
> x <- rnorm(100) 
> y1 <- x[1:25] * 2 + rnorm(25, mean=1)
> y2 <- x[26:50] * 2.6 + rnorm(25, mean=1.5)
> y3 <- x[51:75] * 2.9 + rnorm(25, mean=5)
> y4 <- x[76:100] * 3.5 + rnorm(25, mean=5.5)
> 
> 
> d <- data.frame(x=x, y=c(y1,y2,y3,y4), f=factor(rep(letters[1:4],
> each=25))) 
> 
> 
> # now try with contrast.Design(): 
> library(multcomp)
> dd <- datadist(d) 
> options(datadist='dd') 
> l <- ols(y ~ x * f, data=d)
> 
> ## model fitted using successive difference contrasts
> library(MASS)
> T.lm <- lm(y ~ x * f, contrasts=list(f="contr.sdif"), data=d)
> summary(T.lm)
> 
> ## show the warning: model fitted using ols() and treatment contrasts
> summary(glht(l, linfct=mcp(f="Seq")))
> 
> ## "custom" successive difference contrasts using contrast.Design
> TFin <- {}
> for (i in 1:(length(levels(d$f))-1)) {
> tcont <- contrast(l, a=list(f=levels(d$f)[i]),
> b=list(f=levels(d$f)[i+1]))
> TFin <- rbind(TFin, tcont)
> rownames(TFin)[i] <-  paste(levels(d$f)[i], levels(d$f)[i+1], sep="-")
> }
> TFin[,1:9]
> 
> Regards, Mark.
> 
> 
> 
> Dylan Beaudette-2 wrote:
>> 
>> On Mon, Feb 16, 2009 at 5:28 PM, Patrick Giraudoux
>>  wrote:
>>> Greg Snow a écrit :
 One approach is to create your own contrasts matrix:


> mycmat <- diag(8)
> mycmat[ row(mycmat) == col(mycmat) + 1 ] <- -1
> mycmati <- solve(mycmat)
> contrasts(agefactor) <- mycmati[,-1]
>

 Now when you use agefactor, the intercept will be the first age group
 and the slopes will be the differences between the pairs of groups
 (make sure that the order of the levels of agefactor is correct).

 The difference between this method and the contr.sdif function in MASS
 is how the intercept will end up being interpreted (and the dimnames).

 Hope this helps,


>>>
>>> Actually, exactly what I needed including the reference to contr.sdif in
>>> MASS I did not spot before (although I am a faithful reader of the
>>> yellow book... but so many things still escape to me). Again thanks a
>>> lot.
>>>
>>> Patrick
>>>
>> 
>> Glad you were able to solve your problem. Frank replied earlier with
>> the suggestion to use contrast.Design() to perform the tests after the
>> initial model has been fit. Perhaps I am a little denser than normal,
>> but I cannot figure out how to apply contrast.Design() to a similar
>> model with several levels of a factor. Example data:
>> 
>> # need these
>> library(lattice)
>> library(Design)
>> 
>> # replicate an important experimental dataset
>> set.seed(10101010)
>> x <- rnorm(100)
>> y

Re: [R] Survival-Analysis: How to get numerical values from survfit (and not just a plot)?

2009-02-17 Thread Eik Vettorazzi

Hi Bernhard,
I'm wondering what you will expect to get in "dividing" two proportional 
survival curves from a fitted cox model.
Anyway, you can provide a newdata object to the survfit function 
containing any combination of cofactors you are interested in and then 
use summary, eg:


fit <- coxph( Surv(futime,fustat)~resid.ds+rx+ecog.ps,data=ovarian)
summary(survfit( fit,newdata=data.frame(rx=1, ecog.ps=2, resid.ds=0)))

hth.

Bernhard Reinhardt schrieb:

Hi!

I came across R just a few days ago since I was looking for a toolbox 
for cox-regression.


I´ve read
"Cox Proportional-Hazards Regression for Survival Data
Appendix to An R and S-PLUS Companion to Applied Regression" from John 
Fox.


As described therein plotting survival-functions works well 
(plot(survfit(model))). But I´d like to do some manipulation with the 
survival-functions before plotting them e.g. dividing one 
survival-function by another.


survfit() only returns an object of the following structure

 n events median 0.9LCL 0.9UCL
55.000 55.000  1.033  0.696  1.637

Can you tell me how I can calculate a survival- or baseline-function 
out of these values and how I extract the values from the object? I´m 
sure the calculation is done by the corresponding plot-routine, but I 
couldn´t find that one either.


Regards

Bernhard

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html

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


--
Eik Vettorazzi
Institut für Medizinische Biometrie und Epidemiologie
Universitätsklinikum Hamburg-Eppendorf

Martinistr. 52
20246 Hamburg

T ++49/40/42803-8243
F ++49/40/42803-7790

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


[R] How to connect R and WinBUGS/OpenBUGS/LinBUGS in Linux in Feb. 2009

2009-02-17 Thread Paul Heinrich Dietrich

Hi all,
I've managed to get JAGS working on my Ubuntu Hardy Linux with a 32-bit
computer and AMD processors using R 2.8.1.  JAGS is great.  I've read that
JAGS is the fastest, but that hasn't been my experience.  At any rate, I
have more experience with WinBUGS under Windows and would like a version of
that working as well.

It seems like I've read a lot on the subject and tried a lot, but haven't
managed to get BUGS to work yet.  The most success I've had is to install
WinBUGS or OpenBUGS using this method:
http://www.math.aau.dk/~slb/kurser/bayes-08/install.html

What you also need to know is that you need to open Wine and add a drive. 
Although Z is recommended, I haven't been able to specify it, but have
gotten a D drive to work, using:

wine D:/opt/OpenBUGS/winbugs.exe

Using this method, OpenBUGS opens.  Now, to be able to open it with R.  I've
read all sorts of discussions about BRugs (which is no longer on CRAN, but
old versions can still be found), rbugs, and R2WinBUGS (which I'm used to
using on Windows with WinBUGS).  Some people say R2WinBUGS cannot run
OpenBUGS on Linux, some claim they've done it (I think).  It seems the same
thing with everything else.  I've tried making the linbugs and cbugs file
recommended elsewhere online.  It's all very confusing.

Can someone show a method that works currently, along with some sample code? 
I'm also new to Linux, and confused by path conventions.  For example, in
rbugs, it shows an example of a path such as
"/var/scratch/jyan/wine-20040408/wine", and I don't see how to modify this. 
I have no /var/scratch to begin with, and think Wine is installed in
/home/me/.wine...(I don't have Linux in front of me right now).

Please help.  Thanks.
-- 
View this message in context: 
http://www.nabble.com/How-to-connect-R-and-WinBUGS-OpenBUGS-LinBUGS-in-Linux-in-Feb.-2009-tp22058716p22058716.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] ARIMA models

2009-02-17 Thread emj83

is there some sort of R function which can advise me of the best ARIMA(p,q,r)
model to use based on the Schwarz criterion e.g for e.g p=0-5, q =0, r=0-5
or for example p+r< 5???

or is this something I will have to write my own code for?

Thanks Emma
-- 
View this message in context: 
http://www.nabble.com/ARIMA-models-tp22059382p22059382.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] R scripts and parameters

2009-02-17 Thread Duncan Murdoch

On 2/17/2009 10:55 AM, mau...@alice.it wrote:

A couple of weeks ago I asked how it is possible to run an R script (not a 
function) passing some parameters.
Someone suggested the function "commandArgs()".
I read the on-line help and found no clarifying example. Therefore I do not 
know how to use it appropriately.
I noticed this function returns the pathname of the R executable which is not 
what I need.

I meant to ask if it is possible to open an R session and launch a script 
passing parameters that the script can retrieve and use itself.
Just like in C I can run a program and call it with  some arguments 


Example_Prog A B C
  
The program "Example_Prog" can acess its own arguments through the data structures "argc" an "argv".


How can I launch an R script simulating the above mechanism ? 
Shall I use source ("script-name") ?

Where are the arguments to be passed, as part of the source call ?
Is the function "commandArgs" to be places as one of the first code lines of 
the script in order to access its own arguments ?
Is there any "commandArgs" usage example ?



Gabor gave you a solution from within R.  If you want to run a script 
from the command line, then use commandArgs(TRUE).  For example, put 
this into the file test.R:


commandArgs(TRUE)

(The TRUE says you only want to see the trailing arguments, not 
everything else on the command line.)


Then from the command line, do

Rscript test.R A B C

and you'll see the output

[1] "A" "B" "C"

Duncan Murdoch

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


Re: [R] spss-file problem with foreign 0.8-32

2009-02-17 Thread Harry Haupt


Peter Dalgaard wrote:
> 
> Yes, something in the logic appears to have gotten garbled.
> 
> It's in this part of read,spss:
> 
> if (is.character(reencode)) {
> cp <- reencode
> reencode <- TRUE
> }
> else if (codepage <= 500 || codepage >= 2000) {
> attr(rval, "codepage") <- NULL
> reencode <- FALSE
> }
> else if (m <- match(cp, knownCP, 0L))
> cp <- names(knownCP)[m]
> 
> if you get to the 2nd "else if" then cp is unset. Possible the attempted
> match should be of codepage? But it still looks wrong: Why restrict to
> codepages between 500 and 2000 when knownCP contains several values
> above 1???
> 
> A workaround is to set reencode="ascii" (or whatever is relevant).
> 
> -- 
>O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
>   c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
>  (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
> ~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

Dear Peter,
thanks, reencode="ascii" fixed it (and leaves just the warning message,
which seems to have no effect).
Best,
Harry

-
---
Centre for Statistics
Bielefeld University, Germany
-- 
View this message in context: 
http://www.nabble.com/spss-file-problem-with-foreign-0.8-32-tp22059259p22060774.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] R scripts and parameters

2009-02-17 Thread Gabor Grothendieck
Try this:

A <- 1
B <- 2
C <- 3
source("myfile.R")

Now the code in myfile can access A, B and C.

On Tue, Feb 17, 2009 at 10:55 AM,   wrote:
> A couple of weeks ago I asked how it is possible to run an R script (not a 
> function) passing some parameters.
> Someone suggested the function "commandArgs()".
> I read the on-line help and found no clarifying example. Therefore I do not 
> know how to use it appropriately.
> I noticed this function returns the pathname of the R executable which is not 
> what I need.
>
> I meant to ask if it is possible to open an R session and launch a script 
> passing parameters that the script can retrieve and use itself.
> Just like in C I can run a program and call it with  some arguments
>
>> Example_Prog A B C
>
> The program "Example_Prog" can acess its own arguments through the data 
> structures "argc" an "argv".
>
> How can I launch an R script simulating the above mechanism ?
> Shall I use source ("script-name") ?
> Where are the arguments to be passed, as part of the source call ?
> Is the function "commandArgs" to be places as one of the first code lines of 
> the script in order to access its own arguments ?
> Is there any "commandArgs" usage example ?
>
> Thank you very much in advance.
> Maura
>
>
>
>
>
>
> tutti i telefonini TIM!
>
>
>[[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] Survival-Analysis: How to get numerical values from survfit (and not just a plot)?

2009-02-17 Thread David Winsemius
I seriously doubt that a survfit object could only contain that  
information. I suspect that you are erroneously thinking that what  
print.survfit offers is the entire story.


What does str(survfit(, data=) ) show you?

> data(aml)
> aml.mdl <- survfit(Surv(time, status) ~ x, data=aml)
# this is with survfit.Design since I load Hmisc/Design by default now.
> str(aml.mdl)
List of 18
 $ n: int 23
 $ time : num [1:20] 9 13 18 23 28 31 34 45 48 161 ...
 $ n.risk   : num [1:20] 11 10 8 7 6 5 4 3 2 1 ...
 $ n.event  : num [1:20] 1 1 1 1 0 1 1 0 1 0 ...
 $ surv : num [1:20] 0.909 0.818 0.716 0.614 0.614 ...
 $ type : chr "right"
 $ ntimes.strata: Named int [1:2] 10 10
  ..- attr(*, "names")= chr [1:2] "1" "2"
 $ strata   : Named num [1:2] 10 10
  ..- attr(*, "names")= chr [1:2] "Maintained" "Nonmaintained"
 $ strata.all   : Named int [1:2] 11 12
  ..- attr(*, "names")= chr [1:2] "Maintained" "Nonmaintained"
 $ std.err  : num [1:20] 0.0953 0.1421 0.1951 0.2487 0.2487 ...
 $ upper: num [1:20] 0.987 0.951 0.899 0.835 0.835 ...
 $ lower: num [1:20] 0.508 0.447 0.35 0.266 0.266 ...
 $ conf.type: chr "log-log"
 $ conf.int : num 0.95
 $ maxtime  : num 161
 $ units: chr "Day"
 $ time.label   : chr "time"
 $ call : language survfit(formula = Surv(time, status) ~ x,  
data = aml)

 - attr(*, "class")= chr "survfit"
>

I also don't think survfit returns a Cox model.


On Feb 17, 2009, at 10:37 AM, Bernhard Reinhardt wrote:


Hi!

I came across R just a few days ago since I was looking for a  
toolbox for cox-regression.


I´ve read
"Cox Proportional-Hazards Regression for Survival Data
Appendix to An R and S-PLUS Companion to Applied Regression" from  
John Fox.


As described therein plotting survival-functions works well  
(plot(survfit(model))). But I´d like to do some manipulation with  
the survival-functions before plotting them e.g. dividing one  
survival-function by another.


survfit() only returns an object of the following structure

n events median 0.9LCL 0.9UCL
55.000 55.000  1.033  0.696  1.637

Can you tell me how I can calculate a survival- or baseline-function  
out of these values and how I extract the values from the object? I 
´m sure the calculation is done by the corresponding plot-routine,  
but I couldn´t find that one either.


Regards

Bernhard

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

2009-02-17 Thread SLCMSR

Hi,

is it possible to wrap a text using Hershey fonts? "\n" does not work!

Thanks in advance,
Martina

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


Re: [R] R crash after fGarch update

2009-02-17 Thread John Kerpel
Prof Ripley:

Many thanks - it did indeed say it cannot find fGarch after I tried your
advice - but a completely clean re-install did the trick.

John

On Tue, Feb 17, 2009 at 1:09 AM, Prof Brian Ripley wrote:

> Start R with --vanilla, or rename youe saved workspace (.RData).
> Then
>
> library(fGarch)
> load(".RData")  # or whatever you renamed it to.
>
> This will either work or (more ikely) tell you it cannot find fGarch or a
> package it depends on).
>
> On Mon, 16 Feb 2009, John Kerpel wrote:
>
> Hi folks!
>> After updating my packages my R seems to have completely crashed as will
>> not
>> start up - even after I installed 2.8.1 from 2.8.0.
>>
>
> You haven't told us your OS: I am guesing Windows.
>
>  I get the following:
>>
>> Fatal error: unable to restore saved data in .Rdata
>>
>> Error in loadNamespeace(name): there is no package called fGarch
>>
>> But I do have a package called fGarch.
>>
>> After I hit ok, it crashes and exits.  I cannot use any functionality at
>> all.  What do I do?
>>
>> John
>>
>>[[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.
>>
>>
> --
> Brian D. Ripley,  rip...@stats.ox.ac.uk
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel:  +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax:  +44 1865 272595
>

[[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] frequency table for multiple variables

2009-02-17 Thread Marc Schwartz
on 02/17/2009 09:06 AM Hans Ekbrand wrote:
> Hi r-help!
> 
> Consider the following data-frame:
> 
>var1 var2 var3 
> 1 314 
> 2 223 
> 3 223 
> 4 44   NA 
> 5 435 
> 6 223 
> 7 343 
> 
> How can I get R to convert this into the following?
> 
> Value 1  2  3  4  5 
> var1  0  3  2  2  0
> var2  1  3  1  2  0 
> var3  0  0  4  1  1


> t(sapply(DF, function(x) table(factor(x, levels = 1:5
 1 2 3 4 5
var1 0 3 2 2 0
var2 1 3 1 2 0
var3 0 0 4 1 1


The key is to turn each column into a factor with explicitly defined
common levels for tabulation. This enables the table result to have a
consistent format across each column, allowing for a matrix to be
created, rather than a list.

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.


[R] R scripts and parameters

2009-02-17 Thread mauede
A couple of weeks ago I asked how it is possible to run an R script (not a 
function) passing some parameters.
Someone suggested the function "commandArgs()".
I read the on-line help and found no clarifying example. Therefore I do not 
know how to use it appropriately.
I noticed this function returns the pathname of the R executable which is not 
what I need.

I meant to ask if it is possible to open an R session and launch a script 
passing parameters that the script can retrieve and use itself.
Just like in C I can run a program and call it with  some arguments 

> Example_Prog A B C
  
The program "Example_Prog" can acess its own arguments through the data 
structures "argc" an "argv".

How can I launch an R script simulating the above mechanism ? 
Shall I use source ("script-name") ?
Where are the arguments to be passed, as part of the source call ?
Is the function "commandArgs" to be places as one of the first code lines of 
the script in order to access its own arguments ?
Is there any "commandArgs" usage example ?

Thank you very much in advance.
Maura






tutti i telefonini TIM!


[[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] spss-file problem with foreign 0.8-32

2009-02-17 Thread Peter Dalgaard
Harry Haupt wrote:
> Hi,
> after updating to foreign version 0.8-32, I experienced the following error
> when I tried to load a SPSS file:
> 
> Fehler in inherits(x, "factor") : objekt "cp" nicht gefunden
> Zusätzlich: Warning message:
> In read.spss("***l.sav", use.value.labels = TRUE, to.data.frame = TRUE) :
>   ***.sav: File-indicated character representation code (1252) looks like a
> Windows codepage
> 
> It worked without problems with earlier versions.
> 
> Thanks for any clues,
> best,
> Harry

Yes, something in the logic appears to have gotten garbled.

It's in this part of read,spss:

if (is.character(reencode)) {
cp <- reencode
reencode <- TRUE
}
else if (codepage <= 500 || codepage >= 2000) {
attr(rval, "codepage") <- NULL
reencode <- FALSE
}
else if (m <- match(cp, knownCP, 0L))
cp <- names(knownCP)[m]

if you get to the 2nd "else if" then cp is unset. Possible the attempted
match should be of codepage? But it still looks wrong: Why restrict to
codepages between 500 and 2000 when knownCP contains several values
above 1???

A workaround is to set reencode="ascii" (or whatever is relevant).

-- 
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - (p.dalga...@biostat.ku.dk)  FAX: (+45) 35327907

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


Re: [R] Survival-Analysis: How to get numerical values from survfit (and not just a plot)?

2009-02-17 Thread Arthur Allignol

Hi,

See ?survfit.object

if fit is the object you get using survfit,
fit$surv will give you the survival probability.

Best,
arthur

Bernhard Reinhardt wrote:

Hi!

I came across R just a few days ago since I was looking for a toolbox 
for cox-regression.


I´ve read
"Cox Proportional-Hazards Regression for Survival Data
Appendix to An R and S-PLUS Companion to Applied Regression" from John Fox.

As described therein plotting survival-functions works well 
(plot(survfit(model))). But I´d like to do some manipulation with the 
survival-functions before plotting them e.g. dividing one 
survival-function by another.


survfit() only returns an object of the following structure

 n events median 0.9LCL 0.9UCL
55.000 55.000  1.033  0.696  1.637

Can you tell me how I can calculate a survival- or baseline-function out 
of these values and how I extract the values from the object? I´m sure 
the calculation is done by the corresponding plot-routine, but I 
couldn´t find that one either.


Regards

Bernhard

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Survival-Analysis: How to get numerical values from survfit (and not just a plot)?

2009-02-17 Thread Bernhard Reinhardt

Hi!

I came across R just a few days ago since I was looking for a toolbox 
for cox-regression.


I´ve read
"Cox Proportional-Hazards Regression for Survival Data
Appendix to An R and S-PLUS Companion to Applied Regression" from John Fox.

As described therein plotting survival-functions works well 
(plot(survfit(model))). But I´d like to do some manipulation with the 
survival-functions before plotting them e.g. dividing one 
survival-function by another.


survfit() only returns an object of the following structure

 n events median 0.9LCL 0.9UCL
55.000 55.000  1.033  0.696  1.637

Can you tell me how I can calculate a survival- or baseline-function out 
of these values and how I extract the values from the object? I´m sure 
the calculation is done by the corresponding plot-routine, but I 
couldn´t find that one either.


Regards

Bernhard

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


Re: [R] joining "one-to-many"

2009-02-17 Thread Monica Pisica


Ok,
 
I feel properly ashamed. I suppose my "real" data is a little bit different 
than my toy data (although i don't know how) because i did try the merge 
function as simple as merge(t1, t2) and did not work. Maybe a reset of my 
session will solve my problems and more coffee my confusion.
 
Again, thanks for your help,
 
Monica
 

> Date: Tue, 17 Feb 2009 10:09:17 -0500
> Subject: Re: [R] joining "one-to-many"
> From: ggrothendi...@gmail.com
> To: pisican...@hotmail.com
> CC: r-help@r-project.org
>
> Try merge(t1, t2)
>
>
> On Tue, Feb 17, 2009 at 9:33 AM, Monica Pisica wrote:
>>
>> Hello list,
>>
>> I am wondering if a joining "one-to-many" can be done a little bit easier. I 
>> tried merge function but I was not able to do it, so I end up using for and 
>> if.
>>
>> Suppose you have a table with locations, each location repeated several 
>> times, and some attributes at that location. The second table has the same 
>> locations, but only once with a different set of attributes. I would like to 
>> add the second set of attributes to the first table.
>>
>> Example:
>>
>> set.seed <- 123
>> loc <- c(rep("L1", 3), rep("L2", 5), rep("L3", 2))
>> val1 <- round(rnorm(10),2)
>> val2 <- c("a", "b", "c", "a", "b", "d", "f", "e", "b", "e")
>> t1 <- data.frame(loc, val1, val2)
>> t2 <- data.frame(loc=c("L1","L2","L3"), val3 = c("m", "n", "p"), val4 = 
>> c(25, 67, 48))
>>
>> # join one-to-many
>>
>> n <- nrow(t1)
>> m <- nrow(t2)
>> t1$val3 <- rep(1, n)
>> t1$val4 <- rep(1, n)
>>
>> for (i in 1:n) {
>> for (j in 1:m){
>> if (t1$loc[i]==t2$loc[j]) {
>> t1$val3[i] <- as.character(t2$val3[j])
>> t1$val4[i] <- t2$val4[j]
>> }
>> }
>> }
>>
>> Desired result:
>>
>> t1
>> loc val1 val2 val3 val4
>> 1 L1 -0.41 a m 25
>> 2 L1 -0.69 b m 25
>> 3 L1 0.36 c m 25
>> 4 L2 1.11 a n 67
>> 5 L2 0.15 b n 67
>> 6 L2 -0.80 d n 67
>> 7 L2 -0.08 f n 67
>> 8 L2 -1.01 e n 67
>> 9 L3 -1.01 b p 48
>> 10 L3 -2.50 e p 48
>>
>>
>> This code works OK but it is slow if the data frames are actually bigger 
>> than my little example. I hope somebody knows of a better way of doing these 
>> type of things.
>>
>> Thanks,
>>
>> Monica
>> _
>>
>>
>> 22009
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/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] Uninstall question

2009-02-17 Thread Sundar Dorai-Raj
This is on the Mac FAQ:

http://cran.cnr.berkeley.edu/bin/macosx/RMacOSX-FAQ.html#How-can-R-for-Mac-OS-X-be-uninstalled_003f

HTH,

--sundar

On Tue, Feb 17, 2009 at 7:17 AM, ANJAN PURKAYASTHA
 wrote:
> I need to uninstall R 2.7.1 from my Mac. What is the best way to uninstall
> it? Simply delete the R icon in the Applications folder?
> Or is it more involved?
> TIA,
> Anjan
>
> --
> =
> anjan purkayastha, phd
> bioinformatics analyst
> whitehead institute for biomedical research
> nine cambridge center
> cambridge, ma 02142
>
> purkayas [at] wi [dot] mit [dot] edu
> 703.740.6939
>
>[[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] joining "one-to-many"

2009-02-17 Thread ONKELINX, Thierry
Hi Monica,

merge(t1, t2) works on your example. So why don't you use merge?

HTH,

Thierry




ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and 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 Monica Pisica
Verzonden: dinsdag 17 februari 2009 15:33
Aan: R help project
Onderwerp: [R] joining "one-to-many"


Hello list,

I am wondering if a joining "one-to-many" can be done a little bit
easier. I tried merge function but I was not able to do it, so I end up
using for and if.

Suppose you have a table with locations, each location repeated several
times, and some attributes at that location. The second table has the
same locations, but only once with a different set of attributes. I
would like to add the second set of attributes to the first table.

Example:

set.seed <- 123
loc <- c(rep("L1", 3), rep("L2", 5), rep("L3", 2))
val1 <- round(rnorm(10),2)
val2 <- c("a", "b", "c", "a", "b", "d", "f", "e", "b", "e")
t1 <- data.frame(loc, val1, val2)
t2 <- data.frame(loc=c("L1","L2","L3"), val3 = c("m", "n", "p"), val4 =
c(25, 67, 48))

# join one-to-many

n <- nrow(t1)
m <- nrow(t2)
t1$val3 <- rep(1, n)
t1$val4 <- rep(1, n)

for (i in 1:n) {
for (j in 1:m){
if (t1$loc[i]==t2$loc[j]) {
t1$val3[i] <- as.character(t2$val3[j])
t1$val4[i] <- t2$val4[j]
}
}
}

Desired result:

t1
   loc  val1 val2 val3 val4
1   L1 -0.41am   25
2   L1 -0.69bm   25
3   L1  0.36cm   25
4   L2  1.11an   67
5   L2  0.15bn   67
6   L2 -0.80dn   67
7   L2 -0.08fn   67
8   L2 -1.01en   67
9   L3 -1.01bp   48
10  L3 -2.50ep   48


This code works OK but it is slow if the data frames are actually bigger
than my little example. I hope somebody knows of a better way of doing
these type of things.

Thanks,

Monica
_


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

Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message 
and any annex are purely those of the writer and may not be regarded as stating 
an official position of INBO, as long as the message is not confirmed by a duly 
signed document.

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

2009-02-17 Thread ANJAN PURKAYASTHA
I need to uninstall R 2.7.1 from my Mac. What is the best way to uninstall
it? Simply delete the R icon in the Applications folder?
Or is it more involved?
TIA,
Anjan

-- 
=
anjan purkayastha, phd
bioinformatics analyst
whitehead institute for biomedical research
nine cambridge center
cambridge, ma 02142

purkayas [at] wi [dot] mit [dot] edu
703.740.6939

[[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] spss-file problem with foreign 0.8-32

2009-02-17 Thread Harry Haupt

Hi,
after updating to foreign version 0.8-32, I experienced the following error
when I tried to load a SPSS file:

Fehler in inherits(x, "factor") : objekt "cp" nicht gefunden
Zusätzlich: Warning message:
In read.spss("***l.sav", use.value.labels = TRUE, to.data.frame = TRUE) :
  ***.sav: File-indicated character representation code (1252) looks like a
Windows codepage

It worked without problems with earlier versions.

Thanks for any clues,
best,
Harry
-- 
View this message in context: 
http://www.nabble.com/spss-file-problem-with-foreign-0.8-32-tp22059259p22059259.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] joining "one-to-many"

2009-02-17 Thread Gabor Grothendieck
Try merge(t1, t2)


On Tue, Feb 17, 2009 at 9:33 AM, Monica Pisica  wrote:
>
> Hello list,
>
> I am wondering if a joining "one-to-many" can be done a little bit easier. I 
> tried merge function but I was not able to do it, so I end up using for and 
> if.
>
> Suppose you have a table with locations, each location repeated several 
> times, and some attributes at that location. The second table has the same 
> locations, but only once with a different set of attributes. I would like to 
> add the second set of attributes to the first table.
>
> Example:
>
> set.seed <- 123
> loc <- c(rep("L1", 3), rep("L2", 5), rep("L3", 2))
> val1 <- round(rnorm(10),2)
> val2 <- c("a", "b", "c", "a", "b", "d", "f", "e", "b", "e")
> t1 <- data.frame(loc, val1, val2)
> t2 <- data.frame(loc=c("L1","L2","L3"), val3 = c("m", "n", "p"), val4 = c(25, 
> 67, 48))
>
> # join one-to-many
>
> n <- nrow(t1)
> m <- nrow(t2)
> t1$val3 <- rep(1, n)
> t1$val4 <- rep(1, n)
>
> for (i in 1:n) {
>for (j in 1:m){
>if (t1$loc[i]==t2$loc[j]) {
>t1$val3[i] <- as.character(t2$val3[j])
>t1$val4[i] <- t2$val4[j]
>}
>}
> }
>
> Desired result:
>
> t1
>   loc  val1 val2 val3 val4
> 1   L1 -0.41am   25
> 2   L1 -0.69bm   25
> 3   L1  0.36cm   25
> 4   L2  1.11an   67
> 5   L2  0.15bn   67
> 6   L2 -0.80dn   67
> 7   L2 -0.08fn   67
> 8   L2 -1.01en   67
> 9   L3 -1.01bp   48
> 10  L3 -2.50ep   48
>
>
> This code works OK but it is slow if the data frames are actually bigger than 
> my little example. I hope somebody knows of a better way of doing these type 
> of things.
>
> Thanks,
>
> Monica
> _
>
>
> 22009
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/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.


  1   2   >